J. Phys. Chem. 1996, 100, 13561-13568
13561
Simulation of Geminate Recombination in the Presence of an Absorbing Plane N. J. B. Green* and R. D. Spencer-Smith Department of Chemistry, King’s College London, Strand, London WC2R 2LS, U.K. ReceiVed: February 19, 1996; In Final Form: May 3, 1996X
A theoretical investigation is presented of the competition between recombination of a geminate pair and reaction with an absorbing plane surface. Full random flights simulations are compared with the predictions of the independent pairs approximation for a variety of spatial configurations and, where possible, analytic solutions. Differences resulting from the competition are manifested as an underestimate of the amount of recombination and average hit-hit separation in the independent pairs relative to the random flights. These discrepancies are rationalized in terms of a correlation between the pair separation and the closest distance of the pair to the plane. The exact solution of a related one-dimensional problem is presented in order to give further insight into the correlation effect.
1. Introduction A distinct series of processes follows the energy loss from a high kinetic energy electron as it passes through a medium. The energy of the particle is transferred to molecules resulting either in ionization or in excitation. The ionized electrons may go on to ionize and excite further molecules of the media until insufficient energy for ionization or excitation is left. The ionized and excited molecules may then fragment and/or relax back to their ground state. These processes result in small groups of species lying in the wake of the primary electron, collectively called a radiation track. The ensuing kinetics are highly non-homogeneous, since they are strongly influenced by the initial spatial distribution of reactive particles in the track, and have been intensively studied, both by computer simulation1-6 and by experiment.7-12 Experiments can be either direct, where a sample is irradiated and one of the species is observed directly, or indirect, such as measurement of yields in scavenging experiments. Indirect measurements are then related back to the underlying kinetics using a theory or an empirical parametrization.10,13-15 The scope of direct measurements is currently limited by several features: only a few species can be observed in sufficiently low concentrations for experiments, and many processes of interest take place on a time scale that is too short to probe directly using pulse radiolysis. Additionally, practical experiments are often unable to afford precise control over the physical properties of the system. This is particularly true in respect to forming specific species in given orientations and positions relative to each other and, as in this paper, to macroscopic features of the system. Computer simulations avoid these problems, but at the expense of simulating a model system that may or may not correspond accurately to reality. Initial particle distributions are laid down, and the ensuing kinetics are simulated. The challenge is to reproduce the totality of experimental measurements and to gain insights inaccessible by practical means, while retaining a physically and chemically realistic model, thus enabling reliable conclusions to be drawn from the model about the very fast initial kinetics. Work by this group3,4,16-20 has principally centered on development of stochastic techniques, such as the analytic master equation (ME), and Monte Carlo methods, such as random flights (RF) and independent reaction times (IRT) simulations. X
Abstract published in AdVance ACS Abstracts, June 15, 1996.
S0022-3654(96)00487-X CCC: $12.00
These studies have provided insight into the relationship between observables (e.g. kinetics and product yields) and the underlying track structure, and the IRT method is widely used in conjunction with Monte Carlo track structure simulations for modeling radiation effects. However, many problems of interest involve competition between recombination in the track and reaction of track particles with other species. For most situations the required modifications to the techniques are reasonably simple;14 however there are problems where scavenging by a large, almost static moiety is important, such as a surface (radiation enhancement of corrosion), a colloidal particle, or a macromolecule. For these problems it is of interest to be able to model not only the competition kinetics but also the spatial distribution of reactions on the macroscopic moiety; for example, several groups are currently applying stochastic simulation methods to model radical attack on DNA. As yet, no suitable practical experiments have been designed and implemented that can give this type of information. The close proximity of a macroscopic entity will affect the chemistry and kinetics of the species in solution, both by restricting the space available for free diffusion and by acting as a giant scavenger, competing with recombination. It is very likely that theories and simulation methods based on the independent pairs approximation, such as the IRT method, will be applied to such problems in the near future. The purpose of this paper is 2-fold, therefore: firstly to describe how the independent pairs method can be generalized to include a large reactive entity and secondly to test whether the independent pairs approximation (which is also the basis of Smoluchowski theory, for example) is likely to be accurate for such systems or whether there are other qualitative effects that need to be considered more deeply. Although the most important reactive entity of this type is DNA, it is too complex for a test of this sort, because there is no unique measure of distance between events, for example. We have therefore chosen a simple geometry, the plane. Future investigations will extend the theory to cylindrical and spherical moieties. Reaction with a flat plane might be thought of as a crude model of radical attack on a metal surface and can be applied anywhere that the modeled system is so large in relation to the radicals that it is locally flat and two-dimensional. Section 2 describes the modifications necessary to include the reactive plane in the IRT and RF simulation methods, paying special attention to the modeling of the hit positions, which has not been attempted before. The simulators are tested against analytic solutions where possible and compared against each © 1996 American Chemical Society
13562 J. Phys. Chem., Vol. 100, No. 32, 1996
Green and Spencer-Smith
other, searching for possible correlation effects. The main result of this paper, presented in section 3, is that the competition between recombination and reaction with the plane cannot be described fully in the IRT method, where all reactive distances are assumed to diffuse independently of one another. Some interesting three-body correlation effects are observed, and attempts are made to explain their origin and to analyze them, including the exact solution of a one-dimensional system which shows similar effects. Finally a summary of the findings and suggestions for future work are given. 2. Calculation Details Since this investigation is a first attempt to look for correlations and competition effects, only a small number of reactive particles have been used in the calculations (1 or 2), and they are placed in specific configurations relative to the plane at the start of each realization. The simulation is then run forward in time, and the results are averaged over 105-107 repetitions. Particles are removed with probability 1 if they hit the plane z ) 0. In other words, reaction with the plane is assumed to be diffusion-controlled. The details of the basic IRT and RF codes have been discussed fully elsewhere,4 and apart from a general description only the alterations resulting from the introduction of the surface will be covered here. (a) Random Flights. In the random flights (RF) trajectory simulations time is discretized. At each time step δt the particles, which are assumed uncharged, take a random flight sampled from a spherical normal distribution whose variance in each dimension is 2Dδt. This procedure is equivalent to sampling from the Green function for the diffusion equation. If two particles encounter during δt a reaction is counted. Encounter is determined in two ways: either two particle centers are found to lie within a given distance at the end of the time step, or a survival probability is calculated for each pair conditional on their relative positions at the start and end of the time step.21 The method has been exhaustively tested for systems where exact or numerical solutions are known.4 To this scheme were added two similar checks for reaction with the plane. The first, proximity, worked simply by inspecting the coordinates of each particle at the end of the time step to see if they had passed through the plane. The second check incorporates the possibility of the particles crossing and then recrossing the plane in a single time step. A survival probability is calculated conditional on the distance from the particle to the plane at the start and end of δt:21
Ω ) 1 - exp(-z1z2/Dδt)
(1)
where Ω is the probability that a particle does not hit the plane z ) 0 given that it diffuses from z1 to z2 in a time step δt with a diffusion constant D. If a reaction with the plane occurs during δt, then it is necessary to generate x and y coordinates for the position of the hit. It is possible to use the Brownian bridge process18 to obtain an explicit representation of the probability density of the hit position conditional on the values of z1 and z2 and the fact that a hit has taken place. However, it is hardly necessary to go to this level of sophistication in a random flights simulation, since the time step depends on both the interparticle distance and the distance of the closest particle to the plane and is therefore very small under conditions where a hit is likely; we therefore present an approximate method for positioning the hit. Consider the probability density for the x coordinate at any given instant t + δt (0 e e 1), in the time period from t to
Figure 1. Method for estimating plane hit positions for reaction (a) by proximity and (b) on the Brownian bridge.
t + δt, during which the x coordinate diffuses from x1 to x2. The required density is a normal distribution with mean x1 + (x2 - x1) and a variance of 2(1 - )Dδt.21 The y coordinate is independent of x and has an analogous density. Note that the mean values of the x and y coordinates can be obtained by simple linear interpolation between (x1,y1) and (x2,y2). The density of the hit position is independent of the sign of z2, because if the plane is hit at some time in the time step, then there is an equal probability of subsequent diffusion to z2 and to -z2 (this is known as the reflection principle in the theory of diffusion22). In the RF simulation we therefore generate the hit position as follows: construct the linear interpolation from (x1,y1,z1) to (x2,y2,-|z2|) and place the hit at the point where this line crosses the plane z ) 0. This procedure will underestimate the width of the distribution in a given time step, but the time steps are so small that the results will not be significantly affected. The approximation is summarized in Figure 1. (b) Independent Reaction Times. The independent reaction times method is also a stochastic simulation, but it does not follow the particle trajectories explicitly. This makes it much faster, at the expense of introducing the independent pairs approximation, in which the relative motions of all pairs of particles (and of all particles relative to the plane) are assumed to be independent. Models based on this approximation have been studied extensively and generally shown to agree well with the more accurate random flights treatment.17-20 After the initial configuration has been laid down all the interparticle separations are calculated. For each reactive pair the separation distance is used to generate a random reaction time from the marginal probability distribution for that pair conditional on its separation. The method for generating reaction times for radical recombination has been described elsewhere.4 To include reaction with the plane an additional reaction time tp is generated for each particle, conditional on its initial distance from the plane z.
Simulation of Geminate Recombination
tp )
1 z2 4D (erfc-1(U))2
J. Phys. Chem., Vol. 100, No. 32, 1996 13563
(2)
where D is the diffusion constant and U {0 e U < 1} is a uniform random number. These times are added to the list of radical-radical reaction times generated from the pairwise separations. The first reaction occurs at the minimum of these times, and all the reaction times generated for the reacted particle(s) are removed from further consideration. The second reaction occurs at the minimum remaining time and so on. If a particle survives reaction until tp, it is deemed to react with the plane. Since the diffusion of each particle relative to the plane is assumed independent of the relative motion of the radical pairs, the hit position is generated from a two-dimensional Gaussian distribution of circular symmetry with the appropriate variance, 2Dtp. This is the exact result for the joint density of the x and y coordinates of a free particle at tp.22 3. Results The results presented here are divided into two principal sets. Tests with one and two particles are used to validate the accuracy of the methods; these can be thought of as control experiments. Then simulations with two particles are used to provide information on the competition between recombination and reaction with the plane. Test Simulations. A series of simulations with a single particle were carried out in order to test the RF and IRT programs for a system that is analytically soluble. The particle was started at 5, 10, 20, 30, 40, or 50 Å from the plane and was allowed to diffuse normally with a diffusion constant of 0.7 Å2 ps-1. The particle has an infinite time asymptotic reaction probability of 1. However, the realizations were halted at 105 ps if no reaction had occurred. The results are shown in Figure 2a, where it can be seen that the kinetics predicted by the two simulation methods are in excellent agreement. They also agree almost exactly with the analytic solution, which is not shown.
( )
Wp ) erfc
z x4Dt
(3)
In addition to the kinetics the probability distribution of hits on the plane was also simulated. Figure 2b contains the simulated marginal distribution of the x coordinate of the hit for the case when the particle starts at z ) 20 Å from the plane. These data are compared with the analytical solution for this density, which is obtained by integrating the joint probability density of x, y, and t over the domain 0 < t < ∞, -∞ < y < ∞.
p(x) ) ∫0
∫-∞∞
∞
z
x4πDt3
e-z /4Dt 2
1 -(x2+y2)/4Dt dy dt (4) e 4πDt
which gives the well-known result, a Cauchy distribution (ref 22, p 175).
p(x) )
z π(x + z2) 2
(5)
Again the agreement is excellent, even though the simulation was not actually run to infinite time. The survival probability of the particle at the cutoff time of 100 ns is about 0.05. If the particle is started farther from the plane, the distribution broadens, with a consequent fall in the signal-to-noise ratio,
Figure 2. Results for one particle starting at 5, 10, 20, 30, 40, and 50 Å from the plane. (a) The fraction of particles that have hit the plane up to a given time. Lines represent IRT simulations; symbols, RF simulations. The analytic solution (eq 3) is indistinguishable from the simulated curves. (b) Marginal density of the x coordinate of the hit for one particle starting 20 Å from the plane. The diamonds represent the RF simulation, and the dashed lines the IRT simulation. The analytic solution at infinite time (eq 5) is indistinguishable from the simulated densities.
but the RF and IRT curves still lie very close to one another. These results verify that the simulators are working correctly, particularly the method of generating the hit positions on the plane in the RF simulator. The two-particle simulations yield two additional sets of data. The first is the recombination kinetics of the species. The second is the relative distribution of pairs of hits on the plane, when both particles encounter the plane in the same realization. To verify the correct functioning of these parts of the programs, two tests were made. The first was to check the recombination kinetics by starting the particles far enough from the plane (about 1000 Å) to ensure that there was no competition from the plane. The agreement is excellent between both the current RF and IRT codes that include the plane and the analytic solution for the time-dependent reaction probability:
Wr )
( )
r-a a erfc r x4D′t
(6)
In each case there is no recorded reaction with the plane. To test the simulated hit pair separation density, one of the particles is placed touching the plane. This immediately reacts with the plane before any diffusion can occur, with no chance of recombination, thus fixing the position of one of the hits. The second particle is initially placed along the axis perpendicular to the plane passing through the first particle’s initial
13564 J. Phys. Chem., Vol. 100, No. 32, 1996
Green and Spencer-Smith
Figure 3. Distribution of hit separations for pairs of particles hitting the plane with one particle fixed at the origin and one particle starting at 10 or 50 Å from the plane. Again the simulated density is indistinguishable from the analytic solution (eq 7).
position. If the second particle strikes the plane during the time of the simulation, then the distance between the two hits is recorded. The predicted distribution of the separation between the hits is easily calculated by a method similar to that of eq 4, integrating the joint probability density in cylindrical polar coordinates over the domain 0 < t < ∞, 0 < θ < 2π, yielding
p(r) )
rz (r + z2)3/2 2
(7)
and the simulation results, once again devoid of any possible correlation effects, agree to within statistical significance (Figure 3). Two Particle Simulations. Results are now presented from simulations with two particles in which there is real competition between recombination and reaction with the plane. In the first set the particles start 10 Å apart and the pair is oriented parallel to the plane; typical results are illustrated in Figure 4a-d. It can be seen that the IRT method consistently underestimates the amount of recombination and overestimates the number of plane hits in all cases. The difference decreases as the initial distance from the plane increases, corresponding to a reduction of the ability of the plane to interfere with the recombination. This, in conjunction with the test results above, suggests that the presence of the plane leads to a three-body correlation effect. We believe that the origin of this correlation can be explained in the following way. It is shown in Appendix A that the infinitesimal covariance of the pair distance and the z coordinate of the A particle is given by
(zB - zA) r
σ2(r,zA) ) -2D
(8)
(This is the coefficient in the cross term of a diffusion equation for the joint density of r and zA.) A similar formula applies to the covariance of r and zB, with zA and zB interchanged. If zm denotes the smaller of zA and zB, i.e. the z coordinate of the particle closer to the plane, then the infinitesimal covariance of zm and r is given by
|zB - zA| r
σ2(r,zm) ) -2D
(9)
It is clear that zm and r are negatively correlated. An increase in one tends (statistically) to be accompanied by a decrease in the other. In trajectories that lead to recombination r tends to decrease as time proceeds, and so on these trajectories zm is
more likely to increase than in the whole ensemble of possible trajectories. Reaction with the plane results from first-passage of zm to 0, and therefore trajectories where the pair moves toward recombination tend (statistically) to stay away from the plane and so are less likely to hit the plane than would be expected if zm and r were independent, as assumed in the IRT simulation. The negative correlation is a statistical effect (there is no force) keeping the particles away from the plane as they approach one and other or Vice Versa. As a result of this negative correlation between plane hitting and recombination, it is clear that the plane is less effective at competing with recombination than expected on the basis of independent pairs, in accord with the simulation results. A related effect can be seen in the hit pair separations in Figure 4d. In the IRT simulation there is a higher probability of a small hit-hit distance than in the RF simulation. The reason for this can be found in the same negative correlation: trajectories where the plane is hit are more likely to be those where the particles diffuse away from one another than in the measure of all possible trajectories. Trajectories where the particles diffuse toward one another are more likely to recombine, and as argued above, the plane is less efficient at scavenging them than expected on the basis of the independent pairs approximation. Since the plane preferentially picks up pairs which diffuse apart, and this correlation is absent in the IRT model, the distance between the hits is typically larger than found in the IRT simulation. We have also investigated spatial configurations in which the pair is initially arranged perpendicular to the plane. The particles are again 10 Å apart but are placed on an axis orthogonal to the plane. The results presented in Figure 5 are indexed by the distance from the plane to the center of mass of the particles. For example, in the 5 Å run the A particle (which always starts nearer to the plane in the simulations) actually touches the plane and the B particle is 10 Å further away. It can be seen from Figure 5a-d that there is now a larger difference between the RF and IRT methods than for the simulations where the particles started parallel to the plane. It can be seen from the 5 Å run, which is equivalent to a singleparticle run because the A particle reacts instantly, that the simulators are still working properly. For the true two-particle systems the IRT method underestimates the recombination and overestimates the plane hits as before. The discrepancy becomes smaller as the initial distance from the plane is increased. Perhaps the most surprising result can be seen in Figure 5b, where the A plane hits for the IRT and RF cross over. Qualitatively, these results are not surprising in the light of the discussion above. The negative correlation between zm and r is maximized by placing the pair perpendicular to the plane, since |zA - zB| takes its maximal value of r, whereas it takes its minimum value (0) in the parallel case. The crossover arises because the A hits are more likely to come from trajectories where the particles start diffusing apart and where A initially diffuses toward the plane. Thus we see more A hits at short times than expected from independent pairs. The long-time behavior is more strongly influenced by the relative inefficiency of the plane in competing with recombination discussed above. Thus, in spite of the fact that the IRT method overestimates plane hits in total, the A plane hits that do occur tend to take place faster than predicted by the IRT model. For the same reason the B particles that hit the plane are more likely to have diffused away from the plane initially, and so the plane hits are slower than predicted by the IRT method. This effect means that the IRT simulation tends to underestimate the time of the
Simulation of Geminate Recombination
J. Phys. Chem., Vol. 100, No. 32, 1996 13565
Figure 4. Results for two particles starting parallel to the plane,10 Å apart, and at various distances from the plane. Lines represent independent reaction times simulations; symbols represent random flights simulations. (a) Fraction of pairs recombining up to a given time. (b) The probability of the A particle hitting the reactive plane up to a given time. (c) The A particle survival probability. (d) Simulated density of separation between pairs of plane hits, starting at 10 or 50 Å from the plane. Results for the B particle are identical to those for A.
second hit and therefore also the width of the distribution of the separation between the two hits. Another way of thinking about the competition in this latter case is to separate the inner and outer particles. The space available to the inner (A) particle is restricted by the presence of the outer particle. The IRT approximation does not recognize this restriction and generates, for example, a reaction time for the A particle with the plane as if the B particle were not present, i.e. including all trajectories which pass through the B particle. In reality those trajectories where A initially diffuses toward the plane are more likely to react with the plane, and those where it diffuses away from the plane are more likely to recombine. Such a competition effect will always be found in a system confined between two boundaries. Diffusion relative to the two boundaries will not in general be independent. The situation is complicated here by the fact that the outer boundary (the B particle) is mobile. We analyze this effect in a simple one-dimensional model described in the next subsection. One-Dimensional Model. To analyze the differences in the results from the RF and IRT simulations, where the pair is initially perpendicular to the reactive surface, we consider a simplified one-dimensional model, which can be solved analytically. In this model the two particles diffuse along a line from which they cannot deviate. Particles are assumed to react if they encounter one another or if either particle hits the coordinate 0. The space available to the inner (A) particle, closest to the plane, is bounded by the plane and the outer particle. The outer particle is bounded only by the inner particle, or the plane if the inner particle has reacted. The particle positions are denoted
by ξ and η, respectively, and ξ0 and η0 are the starting positions, ξ0 < η0. Exact Solution. If the particles diffuse independently and without drift, then the instantaneous positions of the two particles can be represented as a point (ξ,η) in the two-dimensional ξη plane. The first reaction must either be particle A hitting the plane, which is equivalent to the point hitting the line ξ ) 0 for the first time, or a recombination of the two particles, which is equivalent to hitting the line ξ ) η. Since the coordinate of particle A must always be less than that of particle B (otherwise they must have met and reacted), the point diffuses in a 45° wedge 0 < ξ < η of the two-dimensional ξη plane, with absorbing boundaries on the lines ξ ) 0 and ξ ) η. The exact solution of this two-dimensional diffusion problem can be simplified by the application of two transformations. The method of images can be used to impose the boundary at ξ ) 0 (plane hitting) if the two-dimensional diffusion can be solved with absorbing boundaries at ξ ) η and ξ ) -η. If the Green function for the position of the diffusing point in this 90° wedge is denoted p90(ξ0,η0;ξ,η,t), then the required probability density in the 45° wedge can be obtained by introducing an image at the point (-ξ0,η0) (equivalent to the reflection principle22). The solution for the required density is equal to p90(ξ0,η0;ξ,η,t) p90(-ξ0,η0;ξ,η,t). It remains to find the Green function for the diffusion p90(ξ0,η0;ξ,η,t). However, this can be seen to be a 45° rotation of two-dimensional Brownian motion absorbed on ξ ) 0 and η ) 0, whose solution is straightforward. By solving this system, rotating the coordinates 45° back to their original orientation and using the method of images, it is possible to calculate the
13566 J. Phys. Chem., Vol. 100, No. 32, 1996
Green and Spencer-Smith
Figure 5. Results for two particles starting orthogonal to the plane 10 Å apart. The simulations start with the center of the two particles either at 5, 10, 20, 30, 40, or 50 Å from the plane. The particle starting closer to the plane is referred to as the A particle; the other, as the B particle. In the 5 Å case the A particle is touching the plane. Lines represent independent reaction times simulations; symbols represent random flights simulations. (a) Fraction of pairs recombining at a given time. (b) A particle plane hit kinetics. (c) B particle plane hit kinetics. (d) Hit separation density.
probability flux crossing the boundaries ξ ) 0 and ξ ) η. This gives analytic solutions for the marginal probability density functions of the reaction time with the plane and the recombination time, respectively. For the reaction of particle 1 with the plane,
f(t1) )
1
[
x4πDt13
ξ0e-ξ0 /4Dt1 erf 2
( ) η0
x4Dt1
-
-η02/4Dt1
η0e
For recombination of the two particles,
f(tx) )
1
x4πD′tx3
[
( )]
erf
η0 + ξ0
x4D′tx
(η0 + ξ0)e-(η0+ξ0) /4D′tx erf 2
x4Dt1
( ) ( )]
(η0 - ξ0)e-(η0-ξ0) /4D′tx erf 2
ξ0
-
η0 - ξ0
x4D′tx
(10)
(11)
where D′ ) 2D. By integrating eq 10 over all time, the probability of the A particle ultimately hitting the plane (as opposed to recombining) can be determined:
()
ξ0 4 F∞ ) 1 - tan-1 π η0
(12)
Equations 10 and 11 do not appear to be analytically integrable over the range 0 to t in closed form. We have
therefore integrated them numerically to obtain the corresponding cumulative distributions. In addition, we can obtain the asymptotic behavior of the A particle plane hits by an expansion of eq 10 suitable for large values of t1:
()
ξ0η0(η02 - ξ02) 4 -1 ξ0 F(t) ≈ 1 - tan + ... π η0 12π(Dt )2
(13)
1
Independent Reaction Times. Here the reaction times for recombination and collision with the plane are assumed independent. This means effectively that the particles can move along the line and through each other without reacting. The joint density function is simply the product of the density functions of the A+plane hits, B+plane hits and recombination. By integrating out one variable, we obtain the marginal densities of t1 and tx.
fp(t1) )
fr(tx) )
ξ0
2
x4πDt13
η0 - ξ0
x4πD′tx3
e-ξ0 /4Dt1 erf
( ) ( ) η0 - ξ0
x4D′t1
e-(η0-ξ0) /4D′tx erf 2
ξ0
x4Dtx
(14)
(15)
Integrating eq 14 over {t1 ) 0,∞} we get the probability of the A particle ultimately hitting the plane:
( )
η0 - ξ0 2 F∞ ) 1 - tan-1 π x2ξ0
(16)
Simulation of Geminate Recombination
Figure 6. Analytic solutions of the one-dimensional model. Two particles are allowed to diffuse along a half-line bounded by the origin. They react if they encounter one another of if one encounters the origin. Results are shown for the exact solution and IRT approximation (eqs 10, 14 and 11, 15).
Once again an analytic form of the asymptotic time dependence of the A particle plane hits can be deduced:
( )
ξ0(η0 - ξ0) η0 - ξ0 2 Fp(t1) ≈ 1 - tan-1 + ... (17) π x2ξ0 2x2πDt1 There are several important things to notice about the cumulative probabilities predicted by these solutions, which are shown in Figure 6. One is that the IRT and RF methods predict different final amounts of recombination and that they tend to this value in very different ways. Another is that the plane hit curves cross in a manner similar to that seen in the threedimensional equivalents. The fit at long times between the numerical integrals and the analytic asymptotic approximations is very good, indicating that the numerical integration procedure is accurate. Inspection of eqs 13 and 17 shows that the IRT predicts a 1/t dependence at long times, whereas the leading term in the exact solution is 1/t2. Calculation of the approximate asymptotes of the integrals for a large value of t1 gives a result a few hundredths of a percent under the predicted cumulative total at t1 ) ∞. The behavior of the one-dimensional solutions bears a marked similarity to the simulated behavior of the corresponding threedimensional system. The fact that the one-dimensional solutions reproduce all the unusual features seen in the three-dimensional study leads us to believe that the effects are indeed being caused by the competition between recombination and plane hitting. The effects are clearly not identical to the three-dimensional case because in three dimensions the particles can diffuse around each other. The result is that the presence of a second particle constrains the region of space through which the first particle can freely move and remains available to react with the plane. This manifests itself in a reduction in the amount of reaction with the plane in the RF simulations. 4. Conclusions The results presented in this paper demonstrate that threebody correlations have a significant effect in these simple systems where recombination competes with scavenging by a macroscopic entity, such as a flat surface. The main effect arises because of a negative correlation between the pair separation and the distance from the nearest particle to the plane. This negative correlation means that the plane is less efficient at intercepting recombining radicals than would be expected on
J. Phys. Chem., Vol. 100, No. 32, 1996 13567 the basis of the usual independent pairs approximation, preferentially picking up pairs that diffuse apart (and are therefore less likely to recombine anyway). As a consequence, pairs of hits on the plane also tend to be further apart than predicted for independent pairs. An alternative way of looking at the competition is in terms of a shielding effect where each particle restricts the space available to the other. This effect is most pronounced in trajectories where the competition between recombination and reaction with the plane is maximized. We have deliberately chosen to investigate systems in which there is significant competition between recombination and scavenging by the plane. For most of these systems the IRT method is reasonably accurate; however, there are cases where discrepancies in modeling the probability of plane-hitting approach 15%, and for these cases the asymptotic approach to this probability also seems to be wrong. We conclude that the IRT model should be used with circumspection under conditions where such a competition is likely to be present. If accurate results are required, it is advisable to use the more computationally intensive random flights simulations. It is likely that these competition effects will be greater in spurs containing more than two radicals, and work is under way to analyze systems of this type. These conclusions have implications for the interpretation of the chemistry in realistic situations such as the radiolysis of DNA in solution. Clearly, in attempts to relate observable quantities like strand-breaks to the physical and chemical processes using the IRT method or any other method based on the independent pairs approximation, stringent tests should be carried out in order to ensure that the method is reliable under the particular conditions simulated. This study identifies one type of correlation effect that should be looked for in such systems and investigated further. Work is also in progress to produce a more quantitative description of the three-body correlation effects that are being observed in these simulations. It is hoped that this will lead to simple and usable modifications of the independent reaction times approximation, which will be able to reproduce these effects yielding a fast, reliable, and possibly analytic method for future use. Further work is also in progress extending the range of geometric shapes for the macroscopic reactive entity, increasing the number of radicals in the system and including a homogeneously distributed scavenger so that two types of competition will be present. Acknowledgment. The authors would like to thank the EPSRC for the funding of this work (GR/H 98069) and Dr. S. Pimblott, University of Notre Dame, for stimulating discussions. Appendix Derivation of Eq 9. Equation 9 is an expression for the infinitesimal covariance of the distance between the two particles and the minimum of their two z coordinates. Both of these are random variables and are here denoted R and Zm, respectively. The infinitesimal covariance of the two random variables is the coefficient of the cross term in the second derivatives that appear in the joint diffusion equation for the two quantities and is defined by the conditional expectation23
σ2(r,zm) ) lim [δt-1 E(δR δZm|r,zm)]
(A1)
δtf0
conditional on the values of R and Zm. Initially we consider the infinitesimal covariance of R with one of the two z coordinates, say ZA. δR can be expressed in terms of the interparticle vector R (defined as RB - RA) as follows:
13568 J. Phys. Chem., Vol. 100, No. 32, 1996
Green and Spencer-Smith
δR ) |R + δR| - |R| ) xR2 + 2R‚δR + δR‚δR - R (A2) The increment in the relative position vector δR is a vector of three independent identically distributed (iid) normal random numbers of mean zero and variance 2D′δt, each of which is itself the difference between two iid normal random numbers of mean zero and variance 2Dδt, representing the increments in the coordinates of each particle during δt (the two diffusion coefficients can easily be made different if desired). The higher order moments of these random increments all have expectations which are o(δt), i.e. such that lim o(δt)/δt ) 0, and therefore do not contribute to the infinitesimal covariance. For this reason we take a power expansion of eq A2 to second order in δR. 2
δR )
R‚δR |δR|2 (R‚δR) + ... + R 2R 2R3
(A3)
The expectation value of δR δZA is found by resolving the vector R into components:
[
E(δR δZA) ) E
δZA (R (δX - δXA) + Ry(δYB - δYA) + R x B
]
Rz(δZB - δZA)) + ... (A4) The higher order terms, which have been omitted, are all o(δt) and make no contribution in the limit δt f 0. Because the individual increments in the coordinates of each particle are mutually independent, the only nonzero term of eq A4 when the limit is taken is the term in δZA2, whose expectation is 2Dδt. Hence, the infinitesimal covariance of R and ZA is -2DRz/R or -2D(ZB - ZA)/R. The equivalent equation for the infinitesimal covariance of R and ZB can be obtained simply by permuting the labels A and B on the two particles. If particle A is closer to the plane than particle B (ZB > ZA), it can be seen that the interparticle separation and ZA are
negatively correlated, and if the converse is true, then R and ZB are negatively correlated. Thus, defining Zm as the smaller of ZA and ZB, we find the result of eq 9.
|zB - zA| R
σ2(r,zm) ) -2D
(9)
References and Notes (1) Schwarz, H. A. J. Phys. Chem. 1969, 73, 1928. (2) Burns, W. G.; Sims, H. E.; Goodall, J. A. B. Radiat. Phys. Chem. 1984, 23, 143. (3) Clifford, P.; Green, N. J. B.; Oldfield, M. J.; Pilling, M. J.; Pimblott, S. M. J. Chem. Soc., Faraday Trans. 2 1986, 82, 2673. (4) Green, N. J. B.; Pilling, M. J.; Pimblott, S. M. Radiat. Phys. Chem. 1989, 34, 105. (5) Pimblott, S. M.; LaVerne, J. A. Radiat. Res. 1992, 129, 265. (6) Kaplan, I. G.; Miterev, A. M.; Sukhonosov, V. Ya. Radiat. Phys. Chem. 1990, 36, 493. (7) Jonah, C. D.; Matheson, M. S.; Miller, J. R.; Hart, E. J. J. Phys. Chem 1976, 80, 1267. (8) Chernovitz, A. C.; Jonah, C. D. J. Phys. Chem. 1988, 92, 5946. (9) Buxton, G. V.; Greenstock, C. L.; Helman, W. P.; Ross, A. B. J. Phys. Chem. Ref. Data 1988, 17, 513. (10) Laverne, J. A.; Pimblott, S. M. J. Phys. Chem. 1991, 95, 3196. (11) Fessenden, R. W.; Schuler, R. H. J. Chem. Phys. 1963, 39, 2147. (12) Schuler, R. H. Radiat. Phys. Chem. 1994, 43, 417. (13) Warman, J. M.; Asmus, K.-D.; Schuler, R. H. J. Phys. Chem. 1969, 73, 931. (14) Hummel, A. J. Chem. Phys. 1968, 48, 3268. (15) Green, N. J. B.; Pimblott, S. M. Mol. Phys. 1991, 74, 811. (16) Green, N. J. B. In preparation. (17) Clifford, P.; Green, N. J. B.; Pilling, M. J.; Burns, W. G. Radiat. Phys. Chem. 1987, 30, 125. (18) Pimblott, S. M.; Pilling, M. J.; Green, N. J. B. Radiat. Phys. Chem. 1991, 37, 377. (19) Clifford, P.; Green, N. J. B.; Pilling, M. J.; Pimblott, S. M. J. Phys. Chem. 1987, 91, 4417. (20) Clifford, P.; Green, N. J. B.; Pilling, M. J.; Pimblott, S. M. J. Phys. Chem. 1989, 94, 8025. (21) Green, N. J. B. Mol. Phys. 1988, 65, 1399. (22) Feller, W. An Introduction to Probability Theory and its Applications; Wiley: New York, 1971; Vol. 2. (23) Arnold, L. Stochastic Differential Equations; Wiley: New York, 1971.
JP960487R
