Statistical Investigation of Dissociation Cross-Sections for Diatoms" BY JAMES KECK Avco-Everett Research Laboratory, 2385 Revere Beach Parkway, Everett, Mass., U.S.A. Received 1st February, 1962 An efficient statistical method for investigating the mechanism and rate of simple chemical reactions has been developed. The procedure involves the random selection of representative systems from the middle of the reaction zone, followed by numerical integration of the equations of motion in both directions to determine the complete course of the reaction. The approach has the important advantage over the usual one of sampling outside the reaction zone that the fraction of systems which react is tremendously increased. The method has been used to investigate the dissociation and recombination of H2,02 and I2 in collisions with argon at temperatures of 0.01 and 0.1 the dissociation temperature D/k.A total of 2400 systems were followed, approximately half of which resulted in a reaction. The cross-section for dissociation as a function of the internal energy H12 of the diatom was found to be proportional to [l +(B-H12)/kT]-3.5, where E is the height of the rotational barrier for the diatom. In a previously published paper,l the author proposed a theory of chemical reaction rates by means of which it is possible to obtain a rigorous upper bound to the reaction rate which can be systematically lowered by a variational technique. Although similar theories were proposed more than 20 years ago by Wigner 2 and Horiuto,3 they apparently failed to attract the attention of theoretical chemists and have been lost in the literature until recently.In these theories, the reacting system is represented by a point in an appropriate classical phase space and an upper bound to the reaction rate is obtained by calculating the flux of trajectories crossing a " trial " surface dividing the reactants from the products. The bomd may then be lowered by adjusting the parameters which define the " trial " surface to produce the minimuiii crossing rate. The reitson this procedure leads only to an upper bound for the reaction rate is that, although in principle there exists a class of surfaces which are crossed only once by any one trajectory, in practice the " trial " surfaces chosen will almost never belong to this class and may be crossed many times by single trajectories.Thus, in the calculations of the reaction rate we will have included some trajectories which do not react at all, as well as others which react only after multiple traversals of the " trial " surface. The present work was initially undertaken to investigate the effect of multiple traversals of the " trial " surface on the calculation of the rate of dissociation and recombination of diatoms in collisions with noble gases carried out in ref. (1). The method involves the random sampling of phase-space trajectories crossing a '' trial " surface tangent to the top of the rotational barrier for the diatom, followed by numerical integration of the classical equations of motion in both directions in time * This work has been supported jointly by Headquarters, Ballistic Systems Division, Air Force Systems Command, United States Air Force, under Contract #AF04(694)-33 and Advanced Research Projects Agency, monitored by the Army Rocket and Guided Missile Agency, Army Ordnance Missile Command, United States Army, under Contract #DA-19-020-ORD-5476. 173174 DISSOCIATION OF DIATOMS to determine the complete trajectories. Given the complete trajectories we can deter- mine not only whether reaction occurred and the number of crossings of the " trial " surface involved, but also the dependence of the reaction probability on the initial configuration of the system.This latter possibility, which is made practical by the large fraction of trajectories which react, opens an extremely interesting field for investigation, the first results of which are the determinations of the dissociation cross-sections as a function of vibrational energy presented in fig. 6 of this paper.MATHEMATICAL FORMULATION Following the methods of ref. ( 1 ) an upper bound Ro to the equilibrium rate of a chemical reaction A+B may be obtained in the form where do is an element of the hypersurface S which divides the reactants A and products B, n is the unit normal to do, v is the generalized velocity of a point in phase space and Pe is the equilibrium density of representative points. To obtain a corresponding expression for the true rate we proceed as follows. Let R(i,j) be the partial rate associated with trajectories which make i traversals of S in the direction of n and j traversals in the opposite direction.By definition, where S(i,j) is that part of S crossed by trajectories in the class (i,j). Since the trajec- tories are continuous curves in phase space and S cannot have any " holes ", R(i,j) = 0 unless j = i,i&l. For a reaction to occur, the number of traversals of S which a trajectory makes in the direction of n must exceed the number it makes in the opposite direction. Hence, j = i- 1 and the true reaction rate R = X(i-')R(i,i-f), (3) i where the factor i-1 has been introduced to correct for the fact that a single trajectory crossing S(i,j) i times contributes i times to R(i,j). In general, it will not be possible to calculate R(i,j) because the surfaces S(i,j) are not known. However, we can obtain a statistical estimate R(i,j) by numerical methods. To do this, we select a random sample of trajectories crossing S with a weight function proportional to the flux pe(v.n). We then follow these trajectories through the reaction in both directions from S to determine their histories. The rate R(i,j) = Ro lim N(i,j)/No , No+ co (4) where NO is the number of trajectories followed and N(i,j) is the number falling in the class (i,j). In practice, the success of this procedure for investigating reaction rates depends on the fraction N(i,i-l)/No of reacting trajectories being fairly large. It is here that the present approach gains an important advantage over similar ones which have been made in the past. By starting the calculation in the middle of the reaction zone rather than outside it as is the usual procedure, the fraction of reacting trajectories is tremendously increased and what was a prohibitive calculation in terms of machine- time becomes quite modest.J . KECK 175 CALCULATIONS OF REACTION PROBABILITY We have applied the method outlined above to the problem of estimating the equilibrium rates for the reactions 2H+A+H2 +A 2 0 +A+02 + A 21 -!- A+I2 +A for values of kT/D equal to 0.01 and 0.1.I I 1 I I I I I I I - 10 -5 0 5 10 vet FIG. 1.-Typical trajectory histories for collisions involving 2H and A at kT/D = 0.01. Shown plotted in dimensionless form are the separation r12 of the H atoms, the separation r3 of A from the centre of the 2H and the difference ( B - H I ~ ) between the instantaneous height of the rotational barrier and the energy of the 2H.ve = 1.3 x 1014 sec-1 is the vibrational frequency of H2 in the ground state, /3 = 1.93 x 108 cm-1 is the decay constant appearing in the Morse potential, z2 is the most probable position of the rotational barrier, and re is the equilibrium separation of H2. The calculations were carried out in the centre-of-mass system for the three particles. As in ref. (1) the interaction potential for the three-atom system was assumed to be the sum of the interaction potentials for the individual pairs of atoms. The Morse potential 4 VM(r) = D([l-exp (-pr+gre)]2-1) was used to represent the interaction between the atoms of the diatom and the Mason- Vanderslice potential 5 K,,&) = DMv exp ( - D , AI76 DISSOCIATION OF DIATOMS was used to represent the interaction between an atom of the diatom and the argon atom.For simplicity, the Van der Waals attraction included in the calculations of ref. (1) was neglected in the present work. The trajectories were started on a surface SB tangent to the toy of the rotational barrier for the diatom. The distribution functions for the initial co-ordinates and momenta were obtained by separating the eleven dimensional integral obtained from eqn. (1) as a product of integrals following the methods of ref. (1). The six coupled equations of motion were written in Cartesian co-ordinates and integrated on an IBh4 7090 computer using a 4th order Runge-Kzrtta method. The trajectories were -4 -5 -1.0 -0.5 0 0 - 5 1'0 vet FIG. 2.-Typical trajectory plots for collisions involving 21 and A at kT/D = 0.1.Shown plotted in dimensionless form are the separation r12 of the atoms, the separation rg of the A atom from the centre of the 21, and the difference (B-Ii12) between the energy of the 21 and the instantaneous height of the rotational barrier. ve = 6.5 x 1012 sec-1 is the vibrational frequency of I2 in the ground state, /3 = 1.86 x 108 cm-1 is the decay constant appearing in the Morse potential, z2 is the most probable position of the rotational barrier and re is the equilibrium separation of 12. followed in both directions from the initial point until the potential of interaction with the argon atom fell below kT/100. As a check on the interval size we required that changes in the total energy E and angular momentum L be limited to values AESkT/100 and AL S L/lOO.The former condition was the more difficult to satisfy. For each case, 400 trajectories were followed making a total of 2400 trajectories and the total computer time involved was about 6 h. All of the 12 co-ordinates and momenta giving the relative positions and velocities of the three particles were printed out at the beginning and end of each trajectory. In addition, a tabulation was madeFIG. 1.-Section of paramagnetic resonance spectrum of 10-5 Mc in cubic ZnS. [To ,face page 176.J . KECK 177 of the internal energy H12 of the diatom, the height B of the rotational barrier for the diatom, the separation r12 of the atoms in the diatom, and the separation r3 of the argon atom from the centre of the diatom. Examples of a few trajectory histories taken from this tabulation are shown in fig. 1 and 2.The cases chosen, 2H-t-A at kT/D = 0.01, and 21+A at kT/D = 0.1, represent the two extremes in the ratio of the collision time to the vibrational period of the molecule. From the plots, one can form a general impression of the character of the particle motions, the duration of the collisions, and the excursions in the energy (B-Hl2). Note that (B-Hl2) TABLE DISTRIBUTION OF TRAJECTORIES WITH RESPECT TO CLASS OF REACTION AND NUMBER OF TRAVERSALS OF THE TRIAL SURFACE. THE PARAMETER 8 = kT/D N(b I i,i I b) N(dE/dt > 0) 2H+A 94 28 60 26 74 8 91 19 119 14 108 16 19 1 26 2 35 4 36 11 1 2 3 0 3 2 0 2 1 0 0 0 1 total 122 89 84 110 135 125 20 28 39 48 total l b class N(fl i:i-l 1 b) N(dE/dt> 0) *o 1 -1 341 55 321 74 20+A 146 33 128 31 81 4 77 14 94 8 100 6 15 0 13 1 12 5 28 2 1 7 & 4 5 3 2 0 0 0 0 0 0 0 0 0 400 400 total 181 159 85 91 102 106 15 14 17 30 total l c class N(fI i,i-1 [ b) N ( f I i,i I f> N(b I i,i I b) N(dE/dt> 0) -01 -1 348 50 346 54 2I+A 256 7 282 4 48 0 44 0 69 1 39 1 0 0 5 0 16 3 25 0 1 2 2 0 0 0 0 0 0 0 0 0 0 0 400 400 total 263 286 48 44 70 40 0 5 19 25 total -0 1 -1 389 11 0 400 395 5 0 400178 DISSOCIATION OF DIATOMS approaches a constant value at the beginning and end of each trajectory indicating that the interaction with the third particle has indeed become negligible.Also note that the vibrational period for diatoms crossing the energy surface (B-Hl2) = 0 is substantially larger than the ground-state vibrational period v, 1.In table I, the results have been catalogued according to the type of transition involved and the number of traversals of the energy surface (B-Hl2) = 0 in each direction. We have expanded the notation N(i,j) of the preceding section to indicate the initial and final state of the diatom. Thus, N(f1 i,i- 1 I b) denotes the number of cases in which a trajectory went from a free state f to a bound state b with i cros- sings of (B-Hl2) = 0 in the directionfjb and i- 1 crossings in the direction b+J Note that there were a small number of cases denoted N(d&/dt>O for which the initial traversal of the trial surface was in the wrong direction. This arose from an approxi- mation made in separating the integral in eqn. (1) in which only the interaction between the third body and its nearest neighbour was taken into account.These cases are to be treated as non-reacting. FIG. 3.-Ratio of number N of reacting trajectories to total number No followed as a function of the ratio P12/(P3+P12), where ~ 1 2 is the reduced mass of the diatom and p3 is the reduced mass for the collision. The upper scale shows the ratio of the mass rn1 of an atom of the diatom to the mass m3 of the third body. A number of interesting features are immediately apparent in table 1. First, and most welcome, is the large number of cases in which a reaction occurred. This implies that our original estimates of the rates based on the variational theory were not far wrong. Secondly, the number of trajectories which cross the trial surface more than once in each direction decreases rapidly as the number of crossings in- creases.Thirdly, within the statistical accuracy of the data there is very little indica- tion of a temperature dependence in the results. Finally, there is a relatively weak dependence of the reaction probability on the mass of the recombining atoms. In fig. 3 the fraction of reacting trajectories NIN, = N,lx(i-l)N(f I i,i-1 I b) I ( 5 ) is plotted as a function of pl2/(pl2+p3), where p12 is the reduced mass of the diatomJ, KECK I79 and ,u3 is the reduced mass of the argon and diatom. The data can be represented by the empirical equation Although the present results give no indication of the dependence of the reaction rate on the mass of the " third body ", there is some evidence in the comparison of theory and experiment made in table 1 of ref.(1) which suggests that when the '' third body " is lighter than the reacting atoms the recombination probability is near unity. Eqn. (6) is in accord with this suggestion. DISSOCIATION CROSS-SECTION As a result of the relatively large fraction of reactions occurring, it is possible to obtain information on the reaction probability as a function of conditions in the initial and final states of the system. Our first effort in this direction has been to investigate the number of reactions as a function of the energy of the bound state of the diatom. The results are shown in table 2 where we have tabulated the number of cases AN in which the ratio E = ( B - H 1 2 ) / k T falls between the values given in the first column.When correlated in this manner, the distribution exhibits very little mass or temperature dependence. TABLE 2.-DISTRIBUTION OF REACTING TRAJECTORIES WITH RESPECT TO THE PARAMETER E = (B-H12)/kT. THE PARAMETER 8 = kT/D case & I @ 0.0-0.25 0.25-0.50 0.50- 1 '0 1 -0-2-0 2.0-4.0 4.0- co 2H+A -01 0.1 42 47 23 11 16 5 8 9 8 2 11 0 20+A 2I+A -0 1 0- 1 *01 0.1 75 67 143 145 27 22 44 50 30 23 30 51 16 20 23 25 9 10 13 11 6 1 6 2 N 108 74 163 143 259 284 The data in tabIe 2 may be used to estimate the cross-section O(E) for the dis- To do this, sociation of a diatom as a function of the parameter E = (B-H12)/kT. we simply equate the differential reaction rate given by the present theory to the corresponding expression which defines the cross- section, to obtain where EX,] and [A] are respectively the concentrations of diatoms and argon atoms, Ke is the equilibrium constant for the reaction X2 -t AeZX + A, kr is the recombination180 DISSOCIATION OF DIATOMS rate constant given by eqn.(45) in ref. (l), v3 is the relative velocity of Xz and A and 83 = p3v3/2kT. For a rotating Morse oscillator, the factor where v is the frequency of the rotating oscillator in the state (B-&)/kT = E, and 22 is the most probable position of the rotational maximum which may be obtained from curves in the appendix of ref. (1). FIG. 4.-Curve of 7 = (v/ve)2D/kT as a function of the parameter E. The curve is approximately valid for 1< @re< 9 and 0.01< kTID< 0.1. For large E, 7 is asymptotic to E.Combining eqn. (6), (9) and (10) we obtain where q = (v/ve)2D/kT is shown plotted in fig. 4 as a function of E and 00 is a cross- section of the order of 10-15 cm2 which may be expressed in terms of parameters defined in ref. (1) as The quantity [F’-+F;-]/4 which must be obtained by numerical methods is plotted in fig. 5 for the cases of interest here. The ratio a(e)/oo computed from eqn. (11) using the data in table 2 is shown as a function of the parameter ( 1 + ~ ) in fig. 6. Within the statistical accuracy of the data, all the points can be fitted by the empirical equation (1 3) Whether or not this particular functional form has any basic significance is not clear to the author. However, it does seem somewhat remarkable that the results which a(E)/a, = 1.5(1 + E ) - ~ ” .FIG.3. [To face page 180,FIG. 4.J . KECK 181 13'' 0 = (kT/D) FIG. 5.--curVe of the function (F$-+F/-)/4 required for computing dissociation cross-section from eqn. (12). I I ! I I L L - ~~ I 2 3 4 5 (1 +El FIG. 6.-Calculated dissociation cross-sections as a function of 1 +E. The reference cross-section a0 is defined in eqn. (12) and E = (B-H12)/kT where B is the height of the rotational maximum and H12 is the internal energy of the molecule. The scatter of the points about the line is com- patible with the statistical errors involved.182 DISSOCIATION OF DIATOMS span so wide a range of conditions can be correlated in so simple a manner. It should be noted that (1 + E ) is the mean energy transferred to the diatom divided by kT.DISCUSSION One of the fundamental reasons for our interest in the dissociation cross-section G(E) is that it is independent of the assumption of equilibrium in the vibrational degree of freedom of the diatom. This is important for practical applications since it is virtually certain that during the process of dissociation or recombination in any experimental situation the vibrational states near the dissociation limit will be out of equilibrium. Unfortunately, direct observation of O(E) will probably be extremely difficult. It is most likely that O(E) will be observable only through its influence on the overall dissociation rate. In this connection a theory of coupled vibration and dissociation is being developed by some of the author’s colleagues and will be published shortly.With respect to the integrated reaction rate, the present work produces little change in the comparison between theory and experiment made in ref. (1). The theory predicts the experimental rates at room temperature well but overestimates the rates at high temperatures by a factor of about 5. The mass dependence of the reaction probability uncovered in the present work should slightly improve the cor- relation of the experiments shown in fig. 15 of ref. (1). It should also make possible correlation of recent data on the dissociation of Hz. Since the theory developed here has the potential of giving considerable detail about the mechanics of simple chemical reactions, a restatement of the assumptions which have been made is well worth-while.The most fundamental of these is, of course, that classical mechanics is valid. The main argument here is based on a comparison of the reduced wavelength 2 = fi/(2MkT)* for the heavy particles with the characteristic distance a0 w 0.5 A over which appreciable changes in force occur. The ratio where me is the electron mass, R is the Rydberg energy, and A is the atomic mass. It can be seen that even for hydrogen at room temperature 2/ao-+, while for all heavier molecules and higher temperatures it is substantially less. This implies that it should be possible to describe the particles by reasonably well-defined wave packets which will follow the classical orbits. The second assumption of importance is that the forces between the particles can be uniquely related to their positions and velocities. This requires that electronic transitions do not occur during the course of the collision and is related to the validity of the Born-Oppenheimer separation in chemistry. The final assumption is that the potential can be described as the sum of the interactions between pairs. While this affects the numerical results, it is not basic to the theory which can be used with any potential. In fact, the coupling of theory and experiment may well be useful for obtaining information about the nature of many-body interactions. 3/ao = (m,/M)*(R/kT)* = (86/AT)*, The author wishes to acknowledge the important contributions made to this work by Benjamin Woznik who blocked out the original IBM programme, Kendrick Ownby who made it work, and Anna Manzi who is currently assisting with the calculations. 1 Keck, J. Chem. Physics, 1960, 32, 1035. 3 Horiuto, Bull. Chem. SOC. Japan, 1938, 13, 210. 4 Herzkrg, Spectra of Diatomic Molecules @. Van Nostrand, New York, 1950). 3 Mason and Vandersiice, J. Chem. Physics, 1958,28,432. 2 Wigner, J. Chem. Physics, 1937,5, 720.