- Email: [email protected]
MD simulation of rotational relaxation in liquid nitrogen. Comparison with NMR and Raman data
CHEMICAL PHYSICS LETTERS
Volume 2 15, number 4
3 December 1993
MD simulation of rotational relaxation in liquid nitrogen. Comparison with NMR and Raman data S.I. Temkin ’ and W.A. Steele Department ofChemistry
152
DaveyLaboratory,Penn State University, University Park. PA 16802, USA
Received 9 August 1993; in final form 27 September 1993
A simulation of liquid nitrogen has been performed along the coexistence curve. Angular momentum, rotational energy and second-order orientational correlation times 7,, 7E, ros2were calculated. The first and the last were compared with the available data. The results for sE were compared with halfwidth data from the isotropic Raman scattering measurements. For the entire coexistence curve, the validity of stochastic perturbation theory was shown: ~0~7~~ 1, where wo is an average frequency of the Q branch. The rotational contribution to the halfividth (estimated from the extreme narrowing formula) decreases with cooling, but does not become less than 25% of the whole width. The ratio 7J7, varies from 1.7 to 1 in the region between the critical and the triple points, which is inconsistent with the Langevin description of the angular momentum relaxation.
1. Introduction The purpose of this Letter is to compare the results of a MD simulation of rotational relaxation in dense nitrogen with experiment. Starting with the work of Barojas et al. [ I], nitrogen has been the non-spherical molecule most simulated in the fluid state. Comparisons with the experimental kinetic data usually involve the orientational relaxation times or correlation functions determined from molecular spectroscopy [ 2 1. Such comparisons are not easy due to many-particle contributions to the spectroscopic time-correlation functions. In this respect NMR measurements are simpler to interpret because the spin interactions vanish sufficiently rapidly with distance to break the total magnetic coupling into a sum of single-molecule terms. However, the method gives only correlation times, i.e. the first moments of the corresponding functions [ 3-5 1. The desired NMR experiments were performed in the mid-seventies [ 61. Together with earlier results [ 71, these data provide angular momentum times rJ and orientational relaxation times of the second order rO,*along the nitrogen liquid-vapor coexistence line. For co’ Permanent address: The Institute of Chemical Kinetics and Combustion, Novosibirsk 630090, Russian Federation. Elsevier Science Publishers B.V.
existence conditions, precise measurements of the isotropic Raman line width were also reported [ 8 1. To our knowledge, this array of information has never been used before in checking simulated correlation times. Meanwhile, a pioneering experiment [ 8 ] appeared that was of great importance for the verification of the quasi-classical theory of the isotropic Q branch collapse [ 9 1. However, this theory needs the rotational energy relaxation time rE as an input parameter. Even for rarified gases, a restricted amount of such data is available from routine acoustics methods [ 10 1. In ref. [ 9 1, the Keilson-Storer model (KS) [ 11,12 ] was used to express rE in terms of the measured rJ. In this way, gas [ 13 ] and liquid [ 81 isotropic Raman data were successfully treated even though neither rJ nor rE were known in the gas phase. However, this approach [ 9 ] yielded values for the cross section aJ, which were confirmed by NMR measurements [ 14,15 1. However, the KS model is based on an impact approximation for molecular collisions and thus obtains a strongly nonadiabatic relation between rotational relaxation times: 1 G 7,/7&2. To what extent this may be considered to be an acceptable approximation can be revealed uniquely by numerical experiment. Thus, the objectives of the present publication are: - to perform a MD simulation for bulk nitrogen 285
Volume 2 15, number 4
CHEMICAL PHYSICS LETTERS
along the vapor-liquid coexistence line; - to use the available NMR and Raman data for the same conditions to compare theoretical models and simulation. It will be shown that good agreement between calculation and experiment is obtained for re,2over the whole coexistence curve. Good agreement is found for the r., values as well, but the time dependence of the correlation functions for this quantity are more erratic, so the precision of the simulated values is less. The most difficult calculation was that of TV,because of the large fluctuations in this quantity. It is also shown here that the “pure” rotational contribution to the isotropic Raman width decreases with temperature, but even at the lowest temperature it is not smaller than f of the intrinsic vibrational contribution.
2. Simulation Molecular dynamics simulations were performed for 256 molecules using cubic boundary conditions. A pair-wise, site-site potential was employed with Lennard-Jones site-site parameters equal to E/ k= 36.5 K and u= 3.29 A. The density-temperature dependence along the vapor-liquid coexistence was taken from ref. [ 16 1. The time step was chosen to be 0.0004 in dimensionless units (rn~~/e)‘/~, and thus corresponds to 1.2658 fs. Five thousand steps were run at the equilibration stage, and ten or twenty thousand for the simulation itself. Correlation functions of angular momentum and the Legendre polynomial of the second order were evaluated and integrated, thus giving the desired relaxation times. The influence of a finite integration interval was controlled by shifting the truncation point to larger times, having the number of the time points per correlation function fmed. Consequently, the statistics of the averaging becomes poorer through the decrease of the sampling number. Results are plotted in fig. 1. Solid lines represent the experimental results from ref. [ 61. In fact, the high-temperature Arrhenius part of T_,data is taken from ref. [ 71. To give an idea of the scatter, the most deviated experimental point is shown. Fig. 2 presents the temperature dependence of the rotational energy relaxation times. The convergence 286
3 December 1993
1
‘f----y7 0.1 5 16
15
14
13
12
11
10
9
8
1000/T(K) Fig. 1. Liquid nitrogen correlation times re,r and r,along the liquid-vapor coexistence; simulation versus the NMR experiment. (-) ?I, 15Ndata from refs. [6,7]; (- - -) r@,r,“N data from refs. [ 6,7]; (0) T,, MD simulation; ( + ) zar, MD simulation; (Cl ) one measurement from ref. [ 61 is exposed to demonstrate the largest scatter of the data; ( x ) to give an idea of the scatter of r, from simulation for the highest temperature, a point is drawn for the nm of ten and twenty thousand steps; ( A ) MD simulation for the boiling point [ 11, where ~,r and T, coincide; ( * ) MD simulation from ref. [ 191. Upper point is r,+,r,lower point r,.
I
I
1.6 1.4 -
+
1.2 x j
, _
I
+
i:
i
f
+ *
+ 0
0.6 -
0
0
0.4 -
16
_
0
0.8 -
0.2-o
+I
+ +
+
+
+
0 I 15
o
I 14
on0 I 13
I 12
0 0 0 I 11
I 10
I 9
8
1000/T(K) Fig. 2. Rotational energy relaxation time for liquid nitrogen. (0) rE simulated along the coexistence line; ( + ) ratio 7&, from simulation. For this variable ordinate axis scale is dimensionless; ( 0 ) MD simulation from ref. [ 221 for the boiling point.
of the integrations of these correlation functions over the time was, as before, checked by increasing the integration interval. In addition, independent runs were made for several states. A simple interpolation of the simulation points was made, yielding rE for the temperatures along the coexistence curve for which the isotropic scattering measurements of ref. [8] were made. For the calculation of the average
frequency of the Q branch a,T/T,, [171, a value for the rotation-vibration interaction constant was taken from ref. [ 181, (w,=O.O174 cm-i, and a rotational temperature of T,,,,= fi*/21k= 2.93 K. wQ=
3. Discussion
TJ
rJ,
KJ
7J
=M(t)
.
(1)
Eq. ( 1) links the correlation functions KJ(t)=
-f&(t)
,
(2)
and results in the exact property [ 20 ] TM=
JJL!f(t) df =o
w*>
7~=
( [J(f) -J(O) I*WO) > I-em CM*> * lim
(4)
Due to (2) the integration of be truncated when the simultaneously calculated zM is close to 0. The details of the method are to be found in ref. [ 2 11. The largest fluctuations appeared in the rE calculation, compared with the case, as expected for the correlations of a stochastic variable of higher power. An improvement of the statistics through a larger number of particles or/and an increase in the length of the runs seems to be the only way to proceed further. In some cases, independent runs were employed to estimate the precision of the calculation, as shown in fig. 2. Note also the excellent coincidence with the only known simulation [ 221 of rE in liquid nitrogen. Fig. 2 also shows the ratio Although only nonadiabatic channels contribute to the relaxation of the rotational energy, the angular momentum relaxation has additional channels associated with adiabatic reorientations of J. That is why is a qualitative measure of the adiabaticity of a given relaxation channel. It tends to infinity for a purely adiabatic interaction in which the relaxation involves only changes of direction of J. For the op posite case, is of the order of 1, with the exact values being dependent on a given model for the nonadiabatic relaxation (e.g. 1 for the Gordon extended J-diffusion model [ 41) . Cooling and an increase of density combine to give an increase of the efftciency of the energy relaxation so that both times are nearly equal in the vicinity of the triple point. (In ref. [ 231 it was also suggested that the adiabatic channels are suppressed in dense media. Thus the larger cross section GEobtained for nitrogen at the highest pressures was explained.) To get an idea of the uE temperature dependence we use the MD data from ref. [24] for the isothermal ( T= 300 K) compression of nitrogen. The simulation made in ref. [ 241 for a density of 600 amagat can be compared with ours for the liquid at T= 90 K and nearly the same density (597 amagat). Calculating [24] the corresponding effective cross sea tions from the decay times, we have G~(300 K) = 24.4 A2 (7,9=0.378 PS), (~J~0.285 ps). Introducing an effective cross section for the liquid phase through the gas-kinetic formula KJ(
t)
can
7J
The comparison presented in fig. 1 is a severe test of the MD simulation. Along the coexistence curve, the density changes from 42 1 amagat for 120 K to 689 amagat for 65 K. The earliest comparison of experiment with a MD calculation was made in ref. [ 7 ] for the known values re,z and rJ for 64.1 K [ 19 ] and at the boiling point [ I]. Although these simulations were performed using a different computational scheme, the values differ insignificantly from ours. Fig. 1 shows results from refs. [ I,19 ] chosen because of the proximity of their state conditions to the coexistence line. In ref. [ 11, rJ equals to re,* for T= 77 K, while in our case this happens at 75 K, for 64.1 K rS:< 2 and r$,* > 3 [ 19 1. The NMR measurements are limited by the uncertainty in the quadrupole constant that results in 6W uncertainty in re,z [ 6 1. The scatter of the data is shown in fig. 1. The differences between 7_, values simulated in independent runs are also shown. For all the simulated correlation functions have virtually no negative regions. Nevertheless, the convergence of the correlation function integration with increasing integration interval was not as stable as in the case of 70.2.In this case there is an independent check of the calculation, based on the stochastic equation defining the J(t) evolution along a given simulated trajectory via the total torque M(t), $J(t)
3 December 1993
CHEMICAL PHYSICS LETTERS
Volume 2 15, number 4
’
Taking into account ( 1 ), (3) is expressed as
(3)
rE/7J.
rE/7J
7,_?/7J
rE/zJ=
0~(300K)=32.3
A*
287
Ti=(nDiU)-’
(5)
5
we have aE(90K)=38 A2 (~~~0.445 ps), aJ(90 K) ~52.5 AZ (r.,=O.322 ps). Note that TJ~., changes slightly: from 1.38 at 90 K to 1.32 at 300 K, showing that an increase of temperature causes an increase of the efficiency of the nonadiabatic channels. Correspondingly, ratios of the cross sections are
e(90 K) =
1
56
0,(300K)
*
aA90 K)
= 1 62 * ’
a,(300K)
’
Near the critical point, rE definitely exceeds rJ, their ratio being close to 1.7. Note that the KS theory forbids this ratio to be greater than 1. The most popular Langevin description [ 25 ] of rotational relaxation (corresponding to the particular KS case of extremely correlated collisions) gives rJ=2rE and therefore is also inadequate for this case. Recently [26] it has been shown how the KS theory can be extended beyond this limitation. In particular, the idea was to incorporate negative values of the KS parameter y, while still maintaining the constraint )y ( < 1. For y values close to - 1 the average result of collision is a reversal of the rotation direction. In the classical description of a rotation, such an interaction event is a quasiadiabatical one. In the liquid phase, the model should at least be considered as competitive with other semiquantitative descriptions of the cage effect in the relaxation of dynamical variables. Using the theoretical relation between the times we have [ 12 ] TE -=ZJ
from (6) that y= -0.4. In ref. [ 91 the estimation was y= + 0.3, being obtained from the Raman data [ 81 using an earlier version of the KS theory. A recent generalization [ 271 has shown that the expression for an isotropic Q-branch contour depends on y2. Therefore, y= -0.4, is in accord with the previous result. The negative sign just reverses the inequality of and rE, making it physically proper. Certainly, the KS generalization for negative y still has disadvantages of an impact theory. Nevertheless, it seems to be a convenient model for a preliminary analysis [ 28 ] of a MD result. Analysis [ 26] shows that for 7% - 0.4 the KS model does not yield a negative tail for an angular momentum correlation function. This is also in accord with our simulation results [21]. In fig. 3 the isotropic Raman halfividth data [ 81 are represented in dimensionless variables in which one plots the pure rotational contribution estimated from the motional narrowing formula [ 91 rJ
Ac#2 =c&.
1 1+y’
(6)
rJ
zJ=
0.599
TE/rJ
(7)
Fig. 3 demonstrates that 1 ->lO. 0,
(8)
TE
This shows that the limit of fast fluctuation of the rotational energy is valid for the whole coexistence region. The rotational contribution is comparable 0.1
7
0.09
-
I 0
0.08 $0.07
At T= 120 K, two independent runs of ten thousand and twenty thousand steps were performed. Comparing the corresponding values for and rE indicated in figs. 1 and 2, we have: ps (10 kilo rJ~O.673 ps (20 kilo timesteps), timesteps), rE=0.76 ps (,lO kilo timesteps), rE=0.998 ps (20 kilo timesteps). (In fig. 1, the value re,2 for the run of 20 kilo timesteps is not shown, because it is close to the 10 K s run values: 0.182 and 0.179 ps. ) Thus, the extreme boundaries for the ratio are 1.13 and 1.67. For the greatest value 1.67, we find 288
3 December 1993
CHEMICAL PHYSICS LETTERS
Volume 2 15, number 4
0
o
-
E
0.06 -
x
0.05 -
5
;::; _ 0.02 0.01 0
I 0
lo
20 l/rE
I
30 x
40
50
60
WQ
Fig. 3. Rotational contribution to the isotropic Raman halfwidth. (0) Experimental halfwidths (hwhm) from ref. [S] in dimensionless units; (-) the perturbation theory result, calculated with the help of (7); (Cl ) MD calculated hallividth from ref. [22].
Volume 215, number 4
CHEMICALPHYSICS LETTERS
with the measured widths near the critical point, and provides nearly 20% of the width near the triple point,
in agreement with the prediction in ref. [91.At the boiling point, the MD calculation of the rotational width from ref. [22] is indicated in fig. 3. The fact that the point lies on the curve proves that (7) is an excellent analytical calculation of the Q-branch collapsed rotational structure. Within the assumption of the additivity of the dephasing and rotational contributions, $ of the width at the boiling point T= 77 K . . gives a dephasmg time T,_+230 ps, which can be compared with a very rough estimation Tdpe 150 ps from ref. [ 29 1. It is highly attractive to have Raman scattering as a way of measuring 7E in the liquid state. Whether this is possible depends on the share of other contributions to the Q-branch width and, what may be more important, on the cross correlations between all relaxation channels. Relaxational processes in nitrogen at high densities have been studied very intensively by different techniques. The earliest measurements on compressed gas [ 13 ] and liquid [ 8,3033] were complemented recently by a breakthrough series of measurements: for eight temperatures between the critical and triple points, the stimulated Raman line was recorded under pressure [ 34 1. Raman spectra in the density region 200-800 amagat at room temperature were measured [ 351 and timeresolved picosecond results were reported in ref. [ 361 for the same temperature for even higher densities. For the interpretation of these data, a variety of theoretical approaches have been put forward. The main features of the phenomenon are relatively well understood. It is evident that vibrational energy relaxation does not contribute, leaving “pure” dephasing, resonant energy transfer (which, in effect, can be reduced to vibrational phase modulation [ 371) and rotation-vibration coupling). Having started from analytical attempts to investigate the problem, the theory later incorporated computer simulations. It is not easy to solve the problem by “head-on” attack [ 381 because of the hierarchy of time scales: 10R2 ps for intramolecular vibration, about a fraction of a picosecond to several picoseconds for the relaxation of single molecule orientations and tens or hundreds of picoseconds for the averaged dephasing. The crucial simplification is due to the large relaxation rates of relative coordinates
3 December 1993
and orientations, whose stochastic time evolution affects the intermolecular coupling and dephases the intramolecular vibrations. The applicability of a stochastic perturbation theory provides a way of calculation through the “input noise” autocorrelation functions. Within the adiabatic approximation [ 22,391 this correlation function contains no vibrational coordinates. A criterion for its applicability is well known from magnetic resonance [ 3 1, &07,<<
1
(9)
where the SWis the dispersion of the vibrational frequency shift, and 7, is the correlation time of this time-dependent frequency displacement. Semiquantitative estimations indicated [ 401 that this criterion is fulfilled for both the resonance transfer and pure dephasing channels. The majority of the authors [39,40] treat only dephasing channels and consider them in the motional narrowing limit. Several articles discuss a parallel dephasing channel via the vibration-rotation coupling [ 22,411. Note that in ref. [ 4 11 the rotational contribution was calculated within the Gaussian-Markovian model for angular momentum relaxation. The model is equivalent to Langevin relaxation, therefore 7E= 1 7J. In ref. [ 4 1 ] 7J was expressed via the Enskog time for the calculation. From fig. 2 it is clear that the rotational contribution was diminished by a factor of 2 at T=77 K up to 3.4 at T= 120 IL Finally, ref. [41] questions the fast modulation limit for the long-range attractive force contribution, which results in partial inhomogeneity of the contour. In ref. [ 221 the strong cancellation between different multipole contributions into the dephasing is reported. A negative cross correlation between pure dephasing and rotational terms is assumed and the claim is made that there is no reason to calculate the rotational relaxation separately since it is either negligible or is completely compensated by subtracting the cross correlation. The controversy of the situation is enhanced, if one complements the discussion by gas phase data and their interpretation. In this case, the absence of correlation is assumed successfully for the vibrational and rotational contributions [ 3 5 1. The cross correlation between vibrational dephasing and rotation is still an intriguing and challenging problem. Its existence and negative sign is supported by a MD simulation made in ref. [22] within the 289
Volume 2 15, number 4
CHEMICAL PHYSICS LETTERS
motional narrowing limit (for one temperature). In ref. [ 4 1] the cross correlation was taken to be zero, even taking into account an inhomogeneous broadening resulting from slowly varying attraction forces. Our calculations of the pure rotational contribution to the Q-branch width for all coexistence conditions demonstrate that it is significant compared with the pure dephasing term. Its value at T=77 K is 25% of the whole width, only 5% larger than the smallest value near the triple point. At present, the MD simulation is the only way to clarify this problem.
References [ 1] J. Barojas, D. Levesque and B. Quentrec, Phys. Rev. A 7 (1973) 1092.
[ 21 A. Mueller, W.A. Steele and H. Versmold, in press. [3] AI. Abragam, Principles of nuclear magnetism (Oxford Univ. Press, Oxford, 1983). [4] R.G. Gordon, Advan. Magn. Res. 3 (1968) 1. [ 51W.A. Steele, Advan. Chem. Phys. 34 ( 1976) 1. [6] L.M. Ishol, T.A. Scott and M. Goldblatt, J. Magn. Res. 23 (1976) 313. [7] K Krynicki, E.J. Rahkamaa and J.G. Powles, Mol. Phys. 29 (1975) 539. [8] M.J. Clouter and H. Kiefte, J. Chem. Phys. 66 (1977) 1736. [9] S.I. Ternkin and A.I. Burshtein, Chem. Phys. Letters 66 (1979) 52,57. [ lo] G.J. Prangsma, A.H. Alberga and J.J.M. Beenakker, Physica 64 (1973) 278. [ 111 J. Keilson and J.E. Storer, Quart. Appl. Math. 10 (1952) 243. [ 121 A.I. Burshtein and S.I. Ternkin, Spectroscopy of molecular rotation in gases and liquids (in Russian) (Nauka, Moscow, 1982); (Cambridge Univ. Press, Cambridge), in press. [ 131 A.D. May, J.C. Stryland and G. Varghese, Can. J. Phys. 48 (1970) 2331. [ 141 N.S. Golubev, A.I. Burshtein and S.I. Ten&in, Chem. Phys. Letters 91 (1982) 139. [ 15 ] C.J. Jameson, AK. Jameson and N.C. Smith, J. Chem. Phys. 86 (1987) 6833.
290
3 December 1993
[ 161 J.S. Rowlinson, Liquids and liquid mixtures, 2nd Ed. (Butterworths, London, 1969). [ 171 S.I. Temkin and A.I. Burshtein, JETP Letters 24 ( 1976) 86. [ 181 J. Bendtsen, J. Raman Spectry. 2 (1974) 133. [ 19 ] P.S.Y. Cheung and J.G. Powles, Mol. Phys. 30 ( 1975) 92 1. [20] E.M. Purcell, Phys. Rev. 117 (1960) 828. [ 211 S.I. Temkin and W.A. Steele, to be submitted. [22] D. Levesque, J.-J. Weis and D.W. Oxtoby, J. Chem. Phys. 72 (1980) 2744. [ 231 S.I. Temkin, J.M. Thuet, L. Bonamy, J. Bonamy and D. Robert, Chem. Phys. 158 (1991) 89. [24] A.V. Storoxhevand R.M. Lynden-Bell, Chem. Phys. Letters 183 (1991) 316. 125] P.S. Hubbard, Phys. Rev. A 6 (1972) 2421. [ 261 B.M. Abdmkhmanov, A.I. Burshtein and S.I. Temkin, Chem. Phys. 143 (1990) 297. [ 27 ] S.I. Temkin and A.A. Suvemev, in Springer proceedings in physics, Vol. 63. Coherent Raman spectroscopy, eds. G. Marowsky and V.V. Smimov (Springer, Berlin, 1992) p. 49). [ 281 B.M. Abdrakhmanov and S.I. Ten&in, J. Mol. Liquid, in press. [29] SF. Fischer and A. Laubereau, Chem. Phys. Letters 35 (1975) 6. [ 301 W.R.L. Clements and B.P. Stoicheff, Appl. Phys. Letters 12 (1968) 246. [ 3 1 ] A. Laubereau, Chem. Phys. Letters 27 ( 1974) 600. [32] M. Scotto, J. Chem. Phys. 49 (1968) 5362. [33] S.A. Akhmanov, F.N. Gadjiev, N.I. Koroteev, R.Yu. Orlov and I.L. Shumai, JETP Letters 27 (1978) 243. [34] B. Oksengom, D. Fabre, B. Lavorel, R. Saint-Loup and H. Berger, J. Chem. Phys. 94 (1991) 1774. [ 351 B. Lavorel, B. Oksengom,D. Fabre, R. Saint-Loup and H. Berger, Mol. Phys. 75 (1992) 397. [36] R. Kroon, R. Sprik and A. Lagendijk, Chem. Phys. Letters 161 (1989) 137; R. Kroon, M. Baggen and A. Lagendijk, J. Chem. Phys. 91 (1989) 74. [37] R.M. Lynden-Bell, Mol. Phys. 33 (1977) 907. [38] M. Tuckerman and B.J. Beme, J. Chem. Phys. 98 (1993) 7301. [ 391 D.W. Oxtoby, D. Levesque and J.-J. Weis, J. Chem. Phys. 68 (1978) 5528. [40] W.G. Rothschild, J. Chem. Phys. 65 (1976) 2958. [41] K.S. Schweizer and D. Chandler, J. Chem. Phys. 76 (1982) 2297.
Volume 2 15, number 4
3 December 1993
MD simulation of rotational relaxation in liquid nitrogen. Comparison with NMR and Raman data S.I. Temkin ’ and W.A. Steele Department ofChemistry
152
DaveyLaboratory,Penn State University, University Park. PA 16802, USA
Received 9 August 1993; in final form 27 September 1993
A simulation of liquid nitrogen has been performed along the coexistence curve. Angular momentum, rotational energy and second-order orientational correlation times 7,, 7E, ros2were calculated. The first and the last were compared with the available data. The results for sE were compared with halfwidth data from the isotropic Raman scattering measurements. For the entire coexistence curve, the validity of stochastic perturbation theory was shown: ~0~7~~ 1, where wo is an average frequency of the Q branch. The rotational contribution to the halfividth (estimated from the extreme narrowing formula) decreases with cooling, but does not become less than 25% of the whole width. The ratio 7J7, varies from 1.7 to 1 in the region between the critical and the triple points, which is inconsistent with the Langevin description of the angular momentum relaxation.
1. Introduction The purpose of this Letter is to compare the results of a MD simulation of rotational relaxation in dense nitrogen with experiment. Starting with the work of Barojas et al. [ I], nitrogen has been the non-spherical molecule most simulated in the fluid state. Comparisons with the experimental kinetic data usually involve the orientational relaxation times or correlation functions determined from molecular spectroscopy [ 2 1. Such comparisons are not easy due to many-particle contributions to the spectroscopic time-correlation functions. In this respect NMR measurements are simpler to interpret because the spin interactions vanish sufficiently rapidly with distance to break the total magnetic coupling into a sum of single-molecule terms. However, the method gives only correlation times, i.e. the first moments of the corresponding functions [ 3-5 1. The desired NMR experiments were performed in the mid-seventies [ 61. Together with earlier results [ 71, these data provide angular momentum times rJ and orientational relaxation times of the second order rO,*along the nitrogen liquid-vapor coexistence line. For co’ Permanent address: The Institute of Chemical Kinetics and Combustion, Novosibirsk 630090, Russian Federation. Elsevier Science Publishers B.V.
existence conditions, precise measurements of the isotropic Raman line width were also reported [ 8 1. To our knowledge, this array of information has never been used before in checking simulated correlation times. Meanwhile, a pioneering experiment [ 8 ] appeared that was of great importance for the verification of the quasi-classical theory of the isotropic Q branch collapse [ 9 1. However, this theory needs the rotational energy relaxation time rE as an input parameter. Even for rarified gases, a restricted amount of such data is available from routine acoustics methods [ 10 1. In ref. [ 9 1, the Keilson-Storer model (KS) [ 11,12 ] was used to express rE in terms of the measured rJ. In this way, gas [ 13 ] and liquid [ 81 isotropic Raman data were successfully treated even though neither rJ nor rE were known in the gas phase. However, this approach [ 9 ] yielded values for the cross section aJ, which were confirmed by NMR measurements [ 14,15 1. However, the KS model is based on an impact approximation for molecular collisions and thus obtains a strongly nonadiabatic relation between rotational relaxation times: 1 G 7,/7&2. To what extent this may be considered to be an acceptable approximation can be revealed uniquely by numerical experiment. Thus, the objectives of the present publication are: - to perform a MD simulation for bulk nitrogen 285
Volume 2 15, number 4
CHEMICAL PHYSICS LETTERS
along the vapor-liquid coexistence line; - to use the available NMR and Raman data for the same conditions to compare theoretical models and simulation. It will be shown that good agreement between calculation and experiment is obtained for re,2over the whole coexistence curve. Good agreement is found for the r., values as well, but the time dependence of the correlation functions for this quantity are more erratic, so the precision of the simulated values is less. The most difficult calculation was that of TV,because of the large fluctuations in this quantity. It is also shown here that the “pure” rotational contribution to the isotropic Raman width decreases with temperature, but even at the lowest temperature it is not smaller than f of the intrinsic vibrational contribution.
2. Simulation Molecular dynamics simulations were performed for 256 molecules using cubic boundary conditions. A pair-wise, site-site potential was employed with Lennard-Jones site-site parameters equal to E/ k= 36.5 K and u= 3.29 A. The density-temperature dependence along the vapor-liquid coexistence was taken from ref. [ 16 1. The time step was chosen to be 0.0004 in dimensionless units (rn~~/e)‘/~, and thus corresponds to 1.2658 fs. Five thousand steps were run at the equilibration stage, and ten or twenty thousand for the simulation itself. Correlation functions of angular momentum and the Legendre polynomial of the second order were evaluated and integrated, thus giving the desired relaxation times. The influence of a finite integration interval was controlled by shifting the truncation point to larger times, having the number of the time points per correlation function fmed. Consequently, the statistics of the averaging becomes poorer through the decrease of the sampling number. Results are plotted in fig. 1. Solid lines represent the experimental results from ref. [ 61. In fact, the high-temperature Arrhenius part of T_,data is taken from ref. [ 71. To give an idea of the scatter, the most deviated experimental point is shown. Fig. 2 presents the temperature dependence of the rotational energy relaxation times. The convergence 286
3 December 1993
1
‘f----y7 0.1 5 16
15
14
13
12
11
10
9
8
1000/T(K) Fig. 1. Liquid nitrogen correlation times re,r and r,along the liquid-vapor coexistence; simulation versus the NMR experiment. (-) ?I, 15Ndata from refs. [6,7]; (- - -) r@,r,“N data from refs. [ 6,7]; (0) T,, MD simulation; ( + ) zar, MD simulation; (Cl ) one measurement from ref. [ 61 is exposed to demonstrate the largest scatter of the data; ( x ) to give an idea of the scatter of r, from simulation for the highest temperature, a point is drawn for the nm of ten and twenty thousand steps; ( A ) MD simulation for the boiling point [ 11, where ~,r and T, coincide; ( * ) MD simulation from ref. [ 191. Upper point is r,+,r,lower point r,.
I
I
1.6 1.4 -
+
1.2 x j
, _
I
+
i:
i
f
+ *
+ 0
0.6 -
0
0
0.4 -
16
_
0
0.8 -
0.2-o
+I
+ +
+
+
+
0 I 15
o
I 14
on0 I 13
I 12
0 0 0 I 11
I 10
I 9
8
1000/T(K) Fig. 2. Rotational energy relaxation time for liquid nitrogen. (0) rE simulated along the coexistence line; ( + ) ratio 7&, from simulation. For this variable ordinate axis scale is dimensionless; ( 0 ) MD simulation from ref. [ 221 for the boiling point.
of the integrations of these correlation functions over the time was, as before, checked by increasing the integration interval. In addition, independent runs were made for several states. A simple interpolation of the simulation points was made, yielding rE for the temperatures along the coexistence curve for which the isotropic scattering measurements of ref. [8] were made. For the calculation of the average
frequency of the Q branch a,T/T,, [171, a value for the rotation-vibration interaction constant was taken from ref. [ 181, (w,=O.O174 cm-i, and a rotational temperature of T,,,,= fi*/21k= 2.93 K. wQ=
3. Discussion
TJ
rJ,
KJ
7J
=M(t)
.
(1)
Eq. ( 1) links the correlation functions KJ(t)=
-f&(t)
,
(2)
and results in the exact property [ 20 ] TM=
JJL!f(t) df =o
w*>
7~=
( [J(f) -J(O) I*WO) > I-em CM*> * lim
(4)
Due to (2) the integration of be truncated when the simultaneously calculated zM is close to 0. The details of the method are to be found in ref. [ 2 11. The largest fluctuations appeared in the rE calculation, compared with the case, as expected for the correlations of a stochastic variable of higher power. An improvement of the statistics through a larger number of particles or/and an increase in the length of the runs seems to be the only way to proceed further. In some cases, independent runs were employed to estimate the precision of the calculation, as shown in fig. 2. Note also the excellent coincidence with the only known simulation [ 221 of rE in liquid nitrogen. Fig. 2 also shows the ratio Although only nonadiabatic channels contribute to the relaxation of the rotational energy, the angular momentum relaxation has additional channels associated with adiabatic reorientations of J. That is why is a qualitative measure of the adiabaticity of a given relaxation channel. It tends to infinity for a purely adiabatic interaction in which the relaxation involves only changes of direction of J. For the op posite case, is of the order of 1, with the exact values being dependent on a given model for the nonadiabatic relaxation (e.g. 1 for the Gordon extended J-diffusion model [ 41) . Cooling and an increase of density combine to give an increase of the efftciency of the energy relaxation so that both times are nearly equal in the vicinity of the triple point. (In ref. [ 231 it was also suggested that the adiabatic channels are suppressed in dense media. Thus the larger cross section GEobtained for nitrogen at the highest pressures was explained.) To get an idea of the uE temperature dependence we use the MD data from ref. [24] for the isothermal ( T= 300 K) compression of nitrogen. The simulation made in ref. [ 241 for a density of 600 amagat can be compared with ours for the liquid at T= 90 K and nearly the same density (597 amagat). Calculating [24] the corresponding effective cross sea tions from the decay times, we have G~(300 K) = 24.4 A2 (7,9=0.378 PS), (~J~0.285 ps). Introducing an effective cross section for the liquid phase through the gas-kinetic formula KJ(
t)
can
7J
The comparison presented in fig. 1 is a severe test of the MD simulation. Along the coexistence curve, the density changes from 42 1 amagat for 120 K to 689 amagat for 65 K. The earliest comparison of experiment with a MD calculation was made in ref. [ 7 ] for the known values re,z and rJ for 64.1 K [ 19 ] and at the boiling point [ I]. Although these simulations were performed using a different computational scheme, the values differ insignificantly from ours. Fig. 1 shows results from refs. [ I,19 ] chosen because of the proximity of their state conditions to the coexistence line. In ref. [ 11, rJ equals to re,* for T= 77 K, while in our case this happens at 75 K, for 64.1 K rS:< 2 and r$,* > 3 [ 19 1. The NMR measurements are limited by the uncertainty in the quadrupole constant that results in 6W uncertainty in re,z [ 6 1. The scatter of the data is shown in fig. 1. The differences between 7_, values simulated in independent runs are also shown. For all the simulated correlation functions have virtually no negative regions. Nevertheless, the convergence of the correlation function integration with increasing integration interval was not as stable as in the case of 70.2.In this case there is an independent check of the calculation, based on the stochastic equation defining the J(t) evolution along a given simulated trajectory via the total torque M(t), $J(t)
3 December 1993
CHEMICAL PHYSICS LETTERS
Volume 2 15, number 4
’
Taking into account ( 1 ), (3) is expressed as
(3)
rE/7J.
rE/7J
7,_?/7J
rE/zJ=
0~(300K)=32.3
A*
287
Ti=(nDiU)-’
(5)
5
we have aE(90K)=38 A2 (~~~0.445 ps), aJ(90 K) ~52.5 AZ (r.,=O.322 ps). Note that TJ~., changes slightly: from 1.38 at 90 K to 1.32 at 300 K, showing that an increase of temperature causes an increase of the efficiency of the nonadiabatic channels. Correspondingly, ratios of the cross sections are
e(90 K) =
1
56
0,(300K)
*
aA90 K)
= 1 62 * ’
a,(300K)
’
Near the critical point, rE definitely exceeds rJ, their ratio being close to 1.7. Note that the KS theory forbids this ratio to be greater than 1. The most popular Langevin description [ 25 ] of rotational relaxation (corresponding to the particular KS case of extremely correlated collisions) gives rJ=2rE and therefore is also inadequate for this case. Recently [26] it has been shown how the KS theory can be extended beyond this limitation. In particular, the idea was to incorporate negative values of the KS parameter y, while still maintaining the constraint )y ( < 1. For y values close to - 1 the average result of collision is a reversal of the rotation direction. In the classical description of a rotation, such an interaction event is a quasiadiabatical one. In the liquid phase, the model should at least be considered as competitive with other semiquantitative descriptions of the cage effect in the relaxation of dynamical variables. Using the theoretical relation between the times we have [ 12 ] TE -=ZJ
from (6) that y= -0.4. In ref. [ 91 the estimation was y= + 0.3, being obtained from the Raman data [ 81 using an earlier version of the KS theory. A recent generalization [ 271 has shown that the expression for an isotropic Q-branch contour depends on y2. Therefore, y= -0.4, is in accord with the previous result. The negative sign just reverses the inequality of and rE, making it physically proper. Certainly, the KS generalization for negative y still has disadvantages of an impact theory. Nevertheless, it seems to be a convenient model for a preliminary analysis [ 28 ] of a MD result. Analysis [ 26] shows that for 7% - 0.4 the KS model does not yield a negative tail for an angular momentum correlation function. This is also in accord with our simulation results [21]. In fig. 3 the isotropic Raman halfividth data [ 81 are represented in dimensionless variables in which one plots the pure rotational contribution estimated from the motional narrowing formula [ 91 rJ
Ac#2 =c&.
1 1+y’
(6)
rJ
zJ=
0.599
TE/rJ
(7)
Fig. 3 demonstrates that 1 ->lO. 0,
(8)
TE
This shows that the limit of fast fluctuation of the rotational energy is valid for the whole coexistence region. The rotational contribution is comparable 0.1
7
0.09
-
I 0
0.08 $0.07
At T= 120 K, two independent runs of ten thousand and twenty thousand steps were performed. Comparing the corresponding values for and rE indicated in figs. 1 and 2, we have: ps (10 kilo rJ~O.673 ps (20 kilo timesteps), timesteps), rE=0.76 ps (,lO kilo timesteps), rE=0.998 ps (20 kilo timesteps). (In fig. 1, the value re,2 for the run of 20 kilo timesteps is not shown, because it is close to the 10 K s run values: 0.182 and 0.179 ps. ) Thus, the extreme boundaries for the ratio are 1.13 and 1.67. For the greatest value 1.67, we find 288
3 December 1993
CHEMICAL PHYSICS LETTERS
Volume 2 15, number 4
0
o
-
E
0.06 -
x
0.05 -
5
;::; _ 0.02 0.01 0
I 0
lo
20 l/rE
I
30 x
40
50
60
WQ
Fig. 3. Rotational contribution to the isotropic Raman halfwidth. (0) Experimental halfwidths (hwhm) from ref. [S] in dimensionless units; (-) the perturbation theory result, calculated with the help of (7); (Cl ) MD calculated hallividth from ref. [22].
Volume 215, number 4
CHEMICALPHYSICS LETTERS
with the measured widths near the critical point, and provides nearly 20% of the width near the triple point,
in agreement with the prediction in ref. [91.At the boiling point, the MD calculation of the rotational width from ref. [22] is indicated in fig. 3. The fact that the point lies on the curve proves that (7) is an excellent analytical calculation of the Q-branch collapsed rotational structure. Within the assumption of the additivity of the dephasing and rotational contributions, $ of the width at the boiling point T= 77 K . . gives a dephasmg time T,_+230 ps, which can be compared with a very rough estimation Tdpe 150 ps from ref. [ 29 1. It is highly attractive to have Raman scattering as a way of measuring 7E in the liquid state. Whether this is possible depends on the share of other contributions to the Q-branch width and, what may be more important, on the cross correlations between all relaxation channels. Relaxational processes in nitrogen at high densities have been studied very intensively by different techniques. The earliest measurements on compressed gas [ 13 ] and liquid [ 8,3033] were complemented recently by a breakthrough series of measurements: for eight temperatures between the critical and triple points, the stimulated Raman line was recorded under pressure [ 34 1. Raman spectra in the density region 200-800 amagat at room temperature were measured [ 351 and timeresolved picosecond results were reported in ref. [ 361 for the same temperature for even higher densities. For the interpretation of these data, a variety of theoretical approaches have been put forward. The main features of the phenomenon are relatively well understood. It is evident that vibrational energy relaxation does not contribute, leaving “pure” dephasing, resonant energy transfer (which, in effect, can be reduced to vibrational phase modulation [ 371) and rotation-vibration coupling). Having started from analytical attempts to investigate the problem, the theory later incorporated computer simulations. It is not easy to solve the problem by “head-on” attack [ 381 because of the hierarchy of time scales: 10R2 ps for intramolecular vibration, about a fraction of a picosecond to several picoseconds for the relaxation of single molecule orientations and tens or hundreds of picoseconds for the averaged dephasing. The crucial simplification is due to the large relaxation rates of relative coordinates
3 December 1993
and orientations, whose stochastic time evolution affects the intermolecular coupling and dephases the intramolecular vibrations. The applicability of a stochastic perturbation theory provides a way of calculation through the “input noise” autocorrelation functions. Within the adiabatic approximation [ 22,391 this correlation function contains no vibrational coordinates. A criterion for its applicability is well known from magnetic resonance [ 3 1, &07,<<
1
(9)
where the SWis the dispersion of the vibrational frequency shift, and 7, is the correlation time of this time-dependent frequency displacement. Semiquantitative estimations indicated [ 401 that this criterion is fulfilled for both the resonance transfer and pure dephasing channels. The majority of the authors [39,40] treat only dephasing channels and consider them in the motional narrowing limit. Several articles discuss a parallel dephasing channel via the vibration-rotation coupling [ 22,411. Note that in ref. [ 4 11 the rotational contribution was calculated within the Gaussian-Markovian model for angular momentum relaxation. The model is equivalent to Langevin relaxation, therefore 7E= 1 7J. In ref. [ 4 1 ] 7J was expressed via the Enskog time for the calculation. From fig. 2 it is clear that the rotational contribution was diminished by a factor of 2 at T=77 K up to 3.4 at T= 120 IL Finally, ref. [41] questions the fast modulation limit for the long-range attractive force contribution, which results in partial inhomogeneity of the contour. In ref. [ 221 the strong cancellation between different multipole contributions into the dephasing is reported. A negative cross correlation between pure dephasing and rotational terms is assumed and the claim is made that there is no reason to calculate the rotational relaxation separately since it is either negligible or is completely compensated by subtracting the cross correlation. The controversy of the situation is enhanced, if one complements the discussion by gas phase data and their interpretation. In this case, the absence of correlation is assumed successfully for the vibrational and rotational contributions [ 3 5 1. The cross correlation between vibrational dephasing and rotation is still an intriguing and challenging problem. Its existence and negative sign is supported by a MD simulation made in ref. [22] within the 289
Volume 2 15, number 4
CHEMICAL PHYSICS LETTERS
motional narrowing limit (for one temperature). In ref. [ 4 1] the cross correlation was taken to be zero, even taking into account an inhomogeneous broadening resulting from slowly varying attraction forces. Our calculations of the pure rotational contribution to the Q-branch width for all coexistence conditions demonstrate that it is significant compared with the pure dephasing term. Its value at T=77 K is 25% of the whole width, only 5% larger than the smallest value near the triple point. At present, the MD simulation is the only way to clarify this problem.
References [ 1] J. Barojas, D. Levesque and B. Quentrec, Phys. Rev. A 7 (1973) 1092.
[ 21 A. Mueller, W.A. Steele and H. Versmold, in press. [3] AI. Abragam, Principles of nuclear magnetism (Oxford Univ. Press, Oxford, 1983). [4] R.G. Gordon, Advan. Magn. Res. 3 (1968) 1. [ 51W.A. Steele, Advan. Chem. Phys. 34 ( 1976) 1. [6] L.M. Ishol, T.A. Scott and M. Goldblatt, J. Magn. Res. 23 (1976) 313. [7] K Krynicki, E.J. Rahkamaa and J.G. Powles, Mol. Phys. 29 (1975) 539. [8] M.J. Clouter and H. Kiefte, J. Chem. Phys. 66 (1977) 1736. [9] S.I. Ternkin and A.I. Burshtein, Chem. Phys. Letters 66 (1979) 52,57. [ lo] G.J. Prangsma, A.H. Alberga and J.J.M. Beenakker, Physica 64 (1973) 278. [ 111 J. Keilson and J.E. Storer, Quart. Appl. Math. 10 (1952) 243. [ 121 A.I. Burshtein and S.I. Ternkin, Spectroscopy of molecular rotation in gases and liquids (in Russian) (Nauka, Moscow, 1982); (Cambridge Univ. Press, Cambridge), in press. [ 131 A.D. May, J.C. Stryland and G. Varghese, Can. J. Phys. 48 (1970) 2331. [ 141 N.S. Golubev, A.I. Burshtein and S.I. Ten&in, Chem. Phys. Letters 91 (1982) 139. [ 15 ] C.J. Jameson, AK. Jameson and N.C. Smith, J. Chem. Phys. 86 (1987) 6833.
290
3 December 1993
[ 161 J.S. Rowlinson, Liquids and liquid mixtures, 2nd Ed. (Butterworths, London, 1969). [ 171 S.I. Temkin and A.I. Burshtein, JETP Letters 24 ( 1976) 86. [ 181 J. Bendtsen, J. Raman Spectry. 2 (1974) 133. [ 19 ] P.S.Y. Cheung and J.G. Powles, Mol. Phys. 30 ( 1975) 92 1. [20] E.M. Purcell, Phys. Rev. 117 (1960) 828. [ 211 S.I. Temkin and W.A. Steele, to be submitted. [22] D. Levesque, J.-J. Weis and D.W. Oxtoby, J. Chem. Phys. 72 (1980) 2744. [ 231 S.I. Temkin, J.M. Thuet, L. Bonamy, J. Bonamy and D. Robert, Chem. Phys. 158 (1991) 89. [24] A.V. Storoxhevand R.M. Lynden-Bell, Chem. Phys. Letters 183 (1991) 316. 125] P.S. Hubbard, Phys. Rev. A 6 (1972) 2421. [ 261 B.M. Abdmkhmanov, A.I. Burshtein and S.I. Temkin, Chem. Phys. 143 (1990) 297. [ 27 ] S.I. Temkin and A.A. Suvemev, in Springer proceedings in physics, Vol. 63. Coherent Raman spectroscopy, eds. G. Marowsky and V.V. Smimov (Springer, Berlin, 1992) p. 49). [ 281 B.M. Abdrakhmanov and S.I. Ten&in, J. Mol. Liquid, in press. [29] SF. Fischer and A. Laubereau, Chem. Phys. Letters 35 (1975) 6. [ 301 W.R.L. Clements and B.P. Stoicheff, Appl. Phys. Letters 12 (1968) 246. [ 3 1 ] A. Laubereau, Chem. Phys. Letters 27 ( 1974) 600. [32] M. Scotto, J. Chem. Phys. 49 (1968) 5362. [33] S.A. Akhmanov, F.N. Gadjiev, N.I. Koroteev, R.Yu. Orlov and I.L. Shumai, JETP Letters 27 (1978) 243. [34] B. Oksengom, D. Fabre, B. Lavorel, R. Saint-Loup and H. Berger, J. Chem. Phys. 94 (1991) 1774. [ 351 B. Lavorel, B. Oksengom,D. Fabre, R. Saint-Loup and H. Berger, Mol. Phys. 75 (1992) 397. [36] R. Kroon, R. Sprik and A. Lagendijk, Chem. Phys. Letters 161 (1989) 137; R. Kroon, M. Baggen and A. Lagendijk, J. Chem. Phys. 91 (1989) 74. [37] R.M. Lynden-Bell, Mol. Phys. 33 (1977) 907. [38] M. Tuckerman and B.J. Beme, J. Chem. Phys. 98 (1993) 7301. [ 391 D.W. Oxtoby, D. Levesque and J.-J. Weis, J. Chem. Phys. 68 (1978) 5528. [40] W.G. Rothschild, J. Chem. Phys. 65 (1976) 2958. [41] K.S. Schweizer and D. Chandler, J. Chem. Phys. 76 (1982) 2297.