Collective excitations in liquid methanol: A comparison of molecular, lattice-dynamics, and neutron-scattering results

J. Chem. Phys. 96, 7696 (1992); https://doi.org/10.1063/1.462370 96, 7696 © 1992 American Institute of Physics. Collective excitations in liquid methanol: A comparison of molecular, latticedynamics, and neutron-scattering results Cite as: J. Chem. Phys. 96, 7696 (1992); https://doi.org/10.1063/1.462370 Submitted: 03 July 1991 . Accepted: 14 February 1992 . Published Online: 31 August 1998 J. Alonso, F. J. Bermejo, M. García-Hernández, J. L. Martínez, W. S. Howells, and A. Criado ARTICLES YOU MAY BE INTERESTED IN Molecular-dynamics simulations of dimethylsulfoxide-methanol mixtures The Journal of Chemical Physics 123, 154507 (2005); https://doi.org/10.1063/1.2085052 Raman Spectral Studies of Water Structure The Journal of Chemical Physics 40, 3249 (1964); https://doi.org/10.1063/1.1724992 Investigation of the temperature dependence of dielectric relaxation in liquid water by THz reflection spectroscopy and molecular dynamics simulation The Journal of Chemical Physics 107, 5319 (1997); https://doi.org/10.1063/1.474242 Collective excitations in liquid methanol: A comparison of molecular, latticedynamics, and neutron-scattering results J. Alonsoa) Rutherford Appleton Laboratory, Chilton, Didcot, Oxon OXII OQX, United Kingdom F. J. Bermejo Instituto de Estructura de la Materia, Serrano 123, E-28006 Madrid, Spain M. Garcfa-Hernandeza) Rutherford Appleton Laboratory: Chilton, Didcot, Oxon OXII OQX, United Kingdom J. L. Martinez Institut Laue Langevin, 156X, F-38042 Grenoble Cedex, France W.S. Howells Rutherford Appleton Laboratory, Chilton, Didcot, Oxon OXII OQx, United Kingdom A. Criado Departamento de F'isica Materia Condensada, Universidad de Sevilla, Apartado 1065, E-4080 Sevilla, Spain (Received 3 July 1991; accepted 14 February 1992) The collective dynamics of liquid methanol-d4 is studied by means of molecular-dynamics simulation. The model potential is validated by means of lattice energy calculations and shows a very good agreement with the experimentally obtained crystal structure. Center-of-mass density and momentum fluctuations are investigated in the (Q,w) region which is also accessible to inelastic neutron-scattering (INS) techniques. A simple viscoelastic model previously used for the analysis of INS data is tested against the dynamic structure factor computed from the simulation. A direct comparison with the INS results themselves is also ~ade and qualitative agreement is found. Also, a tentative assignment of the peaks appearing m the current-current correlations is made on the basis of lattice-dynamics calculations for the polycrystalline low-temperature a phase. I. INTRODUCTION The collective dynamical response of strongly associated liquids has attracted a substantial effort from theoretical, computer simulation, and experimental approaches (see, for instance, Ref. 1). Contrary to what was observed in molecular liquids composed of particles interacting via nearly pure Lennard-Jones potentials such as liquid nitrogen/ the presence of strong intermolecular correlations induced by the hydrogen bond was expected to render these systems amenable to experimental measurements using inelastic neutron scattering. As a matter of fact, the possibility of observing short-wavelength collective excitations in the kinematic range accessible to thermal neutrons was pointed out some time ago by Impey, Madden, and McDonald3 on the basis of molecular-dynamics (MD) simulation results for liquid water. A subsequent experiment on inelastic neutron scattering (INS)4 opened up a debate concerning the nature of the observed excitation. On the other hand, the interpretation of the experimental results was put into question by the simulation results of Wojcik and Clementi,5 but a number of papers on a simplified version of the projection operator technique, as well as computer simulations, seemed to reinforce the current interpretation of the neutron-scattering results.6 A study of the coherent inelastic response of liquid a) Pennanent address: Instituto de Estructurade 1a Materia, Serrano 119, E- 28006 Madrid, Spain. methanol was performed later7 and provided additional support for the thesis which assigned an acoustic nature to the observed modes, in opposition to those postulating a "fast-kinetic-mode" nature of the excitation. The reasons for the choice of this material were its relative simplicity compared to water and the fact that its low melting point enabled the study of the collective effects with a very small contribution from single-particle (translation and rotation) modes.8 The main difficulty arising from the analysis of INS spectra was the need to develop a simple analytic model useful for data analysis purposes which encompasses most of the relevant collective dynamical parameters. A relatively simple expression was found7 which, based upon a viscoelastic treatment,9 was able to reproduce the most prominent spectral features, and therefore was used to infer the dispersion behavior from the observed inelastic intensities. The present study has several aims. First of all, it is intended to obtain precise information regarding the collective dynamics of this system (a preliminary account of mostly single-particle results has already been given 10) which is suppossed to be simpler than water since most of the hydrogen-bond network will be constituted by long winding chains 10 instead of the complicated arrangements found in the liquid and solid phases of water. On the other hand, the present exercise will serve to test the validity of the approximation used to analyze the neutron data by performing a similar analysis of the simulated dynamical structure factors. Finally, it will try to shed some light into the origin of 7696 J. Chern. Phys. 96 (10), 15 May 1992 0021-9606/92/107696-14$06.00 @ 1992 American Institute of PhySics Alonso et a/.: Collective excitations in liquid methanol 7697 collective modes other than the normal sound whose physical origin is still an open controversy.4,5 The outline of the paper is as follows. Some relevant formulas concerning density and particle current fluctuations are reviewed in the next section. In Sec. III we summarize the main characteristics of the molecular-dynamics simulations as well as the way in which the intermolecular potential has been validated, and we introduce some results on mostly static properties. Section IV is devoted to the study of the spectral features of the autocorrelation functions of density and particle current fluctuations. A comparison with the relevant functions obtained from a latticedynamical model for the polycrystal will also be made and the origin of high-frequency modes which are shown to be common for the liquid and polycrystalline phases will be discussed. The modeling of S( Q,liJ) and the comparison with inelastic neutron-scattering results are the subject of Sec. V. Finally, our conclusions and summary are presented in the last section. Details about the lattice-dynamics calculations are given in the Appendix. II. THEORETICAL BACKGROUND Dynamic processes in molecular liquids have a greater complexity than in monatomic liquids due to the existence of additional rotational and vibrational degrees of freedom. As a first approximation, it is customary to analyze the dynamics of molecular centers of mass in terms of the formulas developed for monatomic liquids. I 1,12 In this study we will be mainly interested in the collective behavior as characterized by the fluctuations in the density and particle currents. As usual, this will be done through their time-autocorrelation functions (ACFs), F(Q,t) =J... <PQ(t+s)p_Q(s», N J/(Q,t) = ~ ([q'jQ(t+s)] [q'LQ(s)]>, J/(Q,t) =J... ([b'jQ(t+s)][b'LQ(s)]>, (1) N known as the intermediate scattering function and correlation functions of the longitudinal and transverse components of the particle current, respectively. In the preceding equations N is the number of molecules; the average is taken over initial conditions s and wave vectors Q of magnitude Q; q and b are unit vectors parallel and perpendicular, respectively, to the wave vector Q; finally, the spatial Fourier components of the microscopic number density and particle current are defined by N pQ(t) = L exp[ -iQ·Ra(t)], a=1 N jQ (I) = L Va (t) exp[ -iQ·Ra (t)], a=l (2) where Ra (t) and Va (t) are the position and velocity vectors of the center of mass of molecule a at time t. The power spectrum of the intermediate scattering function is known as the dynamic structure factor and is given by I l~ S(Q,liJ) = - dt F(Q,t) cos(liJt) , 1T 0 (3) with similar expressions for the spectra of the current autocorrelation functions (CACFs) J[(Q,t) and Jt (Q,t). From theoretical and computational standpoints it is advantageous to study the short-time behavior of the intermediate scattering function, which provides information about the frequency moments of the dynamic structure factor. Expanding F( Q,t) in a Taylor series, we obtain F(Q,t) = Scm (Q) I - liJ~ -+ liJ~liJ7 - - '" , ( t2 t4 ) 2! 4! (4) where (5) and liJi are the reduced second and fourth frequency moments of S(Q,liJ), respectively. There are analogous expressions for the longitudinal and transverse components of the particle CACF; and their reduced second moments, liJ7(Q) and liJ; (Q), are related to the Q-dependent elastic constants el1 (Q) and G~ (Q), respectively. 11 From the definitions of Eqs. (I) and (2) and exploiting the stationary property of time ACFs it follows that F( Q,t) and J[ ( Q,t) are not independent. \3 In fact, their power spectra are related through II liJ 2S(Q,liJ) = Q2J[(Q,liJ). (6) III. COMPUTATIONAL DETAILS AND STATIC PROPERTIES A. The potential model A substantial number of attempts have been registered where model potentials were developed in order to reproduce thermal and structural properties of the liquid phase as well as single-particle time-correlations functions and the dynamics of bond breaking and forming in the hydrogenbond (HB) network. 14 - 18 However, to the extent of our knowledge no serious attempt to investigate the collective dynamics in methanol from a molecular-dynamics point of view has been made so far. For this kind of computation very long MD runs are needed in order to provide acceptable statistics; therefore, flexible molecular models become impractical since they usually require time steps at least I order of magnitude smaller than the rigid models. Among the latter ones, two have been favored by most researchers. 14,15 As the reliability of both models is comparable,15,19 the choice between them becomes a matter of preference. For this work we have chosen the model proposed by Haughney, Ferrario, and Mc- Donald ls which is a site-site model, each site-site interaction being a Lennard-Jones plus a Coulombic term. B. Lattice energy minimization The adopted model for the intermolecular interactions has been also applied to the low-temperature crystal phase 20 of methanol-d4. The eqUilibrium crystal structure corresponding to this potential model has been obtained as the minimum-energy configuration. For this purpose, a New- J. Chem. Phys., Vol. 96, No.1 0, 15 May 1992 7698 Alonso et al.: Collective excitations in liquid methanol ton-Raphson minimization process with the program 21 WMIN has been carried out starting at the experimental crystal structure, using the cell parameters and molecular rotations and translations as variables. The Ewald method 22 has been used to deal with the Coulombic sums, and the cutoff distance of 12 A has been adopted for the Lennard-Jones interactions. The changes occurring in the process by the lattice parameters a, b, and care 2.9%,3.8%, and 1.6%, respectively, whereas the maximum shift in the atomic coordinates is 0.3 A. This order of discrepancy between experimental and calculated crystal structures is usual in crystal packing calculations using atom-atom potential models; therefore, we can assume that the proposed model reproduces fairly well the experimental crystal structure. C. The MD algorithm Simulations of a 256 CD 3 OD molecule system were carried out using the computer program described in Ref. 10, which was written following the hints given by Allen and Tildesley.23 The program computes the trajectory of the system subjected to cubic periodic boundary conditions using Newtonian classical mechanics (NVE-P ensemble). The molecules were treated as rigid bodies composed of six mass points (modeling the atoms in the molecule) and three interaction sites located at the oxygen, carbon, and hydroxylic deuterium positions. The Cartesian equations of motion were integrated using the velocity version of the Verlet algorithm 23 with a time step of 10 - 14 s, and the RA TILE algorithm 24 was used to implement the holonomic constraints required to keep all intramolecular distances fixed. The molecular geometry and parameters defining the interaction potential were those proposed by Haughney, Ferrario, and McDonald. 15 Only two minor modifications on the potential model were introduced. First, the masses of the hydrogen atoms have been substituted by that of deuterium in order to simulate the fully deuterated compound (which is the one studied by neutron scattering); and, second, the long-range interactions were truncated by using a switch function 23 ,25 to tum off smoothly all the interactions between pairs of molecules whose centers-of-mass separation was greater than the cutoff distance. Following the introduction of the switch function the interaction energy between two molecules, a and/3, becomes Uap(r;a,r jp ) = S(R ~p)U~p(r;a,rjp), (7) where RaP is the distance between centers of mass for both molecules, and Sex) is the unique fifth-order polynomial satisfying {I for x<R L Sex) = 2 o for x>R u, (8) and having continuous first and second derivatives at the end points of the interva1. 25 In the present calculations we took RL = 12.25 A and Ru = 12.65 A for both runs. Note that due to the large cutoff radius and the electroneutrality of methanol molecules the leading term in the truncated potential corresponds to the dipole-dipole interaction and behaves as R;;/. The function U~(r;a,rjp) is the potential function referred to as model HI by Haughney, Ferrario, and McDonald I 5 and has the form of a sum of site-site interactions. D. Accuracy of the results The use of periodic boundary conditions imposes some restrictions upon the length and time scales of the phenomena that can be studied with conventional molecular dynamics. From a dynamical point of view the quantity ofinterest is the recurrence time (defined as 'Tree = LboJc, where c is the velocity of sound) ,26 which is about 1.4 ps for the simulations reported here. In order to ensure that no measurable artifacts were introduced due to recurrence effects, the F( Q,t) for the lowest Q value analyzed (0.25 A - I) was examined in detail and no spurious features that could be attributed to recurrence effects were found for times smaller than 4 ps. Besides, correlation functions computed from finite time-length MD trajectories are subject to numerical and statistical errors;27,23 therefore, the measurable ACF is given by fMD(t) =f(t) + €(t), wheref(t) is the true ACF and €(t) represents the error term. Taking this into account, we see that the spectrum of a MD ACF is given by 1 L"" fMD(W) = - dtfMD(t)W(t) cos(wt) 1T 0 = [few) + €(w)] ® W(w), (9) where w(t) is some window function and W(w) its Fourier transform. Therefore, the attainable spectrum is a convolution of the true spectrum with the window function plus a noise term. Besides, due to the possible existence of recurrence effects as mentioned earlier, the window should decay to zero for times of the order of 'Tree to avoid the contamination of the spectra by recurrence effects. However, this worsens the resolution in the frequency domain and at low Q, where the spectral features ofinterest are narrow and located at low frequency, it seems that there is no choice but to use a wide window though it could result in an increased spectral noise level. The intermediate scattering function and particle current autocorrelation functions have been computed using the fast Fourier transform (FFT) method 23 for both thermodynamic states, and for many values of Q spanning from 0.25 to 5 A - I. In all cases an averaging over the whole set of allowed wave vectors (compatible with the periodic boundary conditions) was performed. In order to obtain an estimation of the quality of our results the 200 K run was subdivided into eight subruns of20 ps length. Autocorrelation functions were computed for these subruns and later averaged and used to estimate the standard deviation of the mean. For Q = 0.25 A -I (that with only three independent wave vectors could be considered a unfavorable case) the error in F(Q,t) was around 10% during the first 5 ps. Fourier transforming the eight subruns separately we found the error in S( Q,w )/S( Q) to be again 10% and constant in the frequency range studied. As the number of allowed wave vectors increases with Q the errors diminish accordingly down to 3 % for Q> 2 A - I. The fast decay of current autocorrelations functions with time J. Chern. Phys., Vol. 96, No. 10, 15 May 1992 Alonso et al.: Collective excitations in liquid methanol 7699 enables us to employ the error analysis of Zwanzig and Ailawadi,27 finding a maximum error of ± 0.05 for the normalized CACFs at low Q which reduces to ± 0.02 as Q increases. E. Some results on static properties Two of the three thermodynamic states studied in Ref. 10 are reanalyzed in this work, with special emphasis on the collective, time-autocorrelation functions. The average values of simple thermodynamic properties were computed along with the trajectories, and are shown together with densities and run lengths in Table I. The run lengths were chosen long enough so that we can average over the longest period fluctuations of the thermodynamics properties observed in the system. The static structure factor as well as the partial paircorrelation functions have been analyzed in a previous paper. 10 The reduced second and fourth moments of S( Q,{i) , {i)o, and (i)[o have been presented in Fig. 1 together with their respective ideal-gas limits. Our results are very similar to those of Wojcik and Clementi for waterS though the frequencies involved are somewhat smaller for methanol. The second moment shows a strong dependence on the static structure factor Q < 2.5 A -I as expected from Eq. (5), displaying a pronounced dip at Qp [i.e., the position of the main maximum of Scm (Q) ], and remaining close to its ideal value from there on. The fourth moment does not approach its ideal limit so quickly though it oscillates around its self-val- ue from relatively low Q. Besides, the dimensionless quotient of (i)/ and its single-particle counterpart, namely (i)r'lf(Q) = (n~ + 3Q2kB TIM) 112, (10) exhibits little or no dependence on the thermodynamic state, n~ being the "Einstein frequency" (see Table I). On the other hand, the oscillations around the self-value decay very slowly, being still noticeable at Q~ 5 A -I though they lose amplitude for Q> 2.5 A-I, which demonstrates that collective effects are important down to length scales of just a few angstroms (i.e., nearest-neighbor distances). Results for the second frequency moment of the spectra of the transverse CACFs also show independence of the thermodynamic state (at least for the two states studied here) when expressed in dimensionless form: {i)t (Q)/{i),;,,'f( Q), where (i)~elf(Q) = (n~ + Q2kBTIM) 112. (11) Collective effects seem to be important only in the region Q < I A - I; for larger momentum transfers the CACF behaves mostly like its single-particle counterpart. The infinite-frequency macroscopic shear modulus, Goo' can be estimated from the low-Q values of (i) t (Q). Approximating Goo by Goo (Q = 0.25 A - I) we get values (Table I) that are about half those found for liquid water.s IV. DENSITY AND CURRENT FLUCTUATIONS A. Density fluctuations The dynamic structure factor S( Q,{i) has been plotted in Figs. 2 and 3 as a function of the angular frequency (i) for the smallest accessible values of the momentum transfer. At 200 K the two lowest-Q spectra (0.25 and 0.35 A - I) exhibit distinctive Brillouin side peaks at frequencies (i) B of 0.44 and 0.56 X 10 13 radls, respectively. Then the peak becomes overdamped as the momentum transfer is increased, though it is still noticeable as a weak shoulder in the spectra for momentum transfers below 0.6 A - I. Another interesting feature of S( Q,{i) is the increase of the intensity in the frequency region around 2X 1013 rad/s, which at Q = 0.86 A -I manifests itself clearly as a shoulder in the spectrum. The presence of this shoulder suggests the existence of a second excitation in the liquid coexisting with the lower-frequency mode. At 300 K the overall situation is very similar, although the Brillouin peaks already appear as overdamped at the lowest value of the momentum transfer accessible to our MD experiment (see Fig. 3). From Fig. 2 it is clear that the peak maxima whenever visible show a noticeable dispersion (at least up to 0.4 A-I). TABLE I. Average values of simple thermodynamic properties for simulated liquid methanol·d4. Run numbers are as in Ref. 10. Quantity Units Run 1 Run 3 Density kglm 3 943.0 885.0 L"' n nm 2.5333 2.5875 Elapsed time ps 177.00 100.10 Temperature" K 202.3 ± 0.3 297.7 ± 1.0 Ulnter b kl/mol -39.88 -34.79 Pressure" MPa 52.0 105.3 n.3 ps -2 463.7 402.2 Scm (0) 0.029 ± 0.009 0.038 ± 0.003 XT 10- '0 Pa-' 6.5 ± 1.9 6.3 ± 0.5 cT mls 1290 ± 190 1340 ±60 CII (Q = 0.25 A. -') 10· J/m 3 15.77 13.16 G~ (Q=0.25A.··') 10· J/m 3 6.72 5.96 "The standard deviation of the average temperature has been estimated by taking into account the statistical inefficiency of the data (Ref. 23). "No long·range corrections have been applied. J. Chern. Phys., Vol. 96, No.1 0, 15 May 1992 7700 Alonso et al. : Collective excitations in liquid methanol 2.5- ..--... rn 2.0- --- '"Cl ro I-< M 1.5- .-t 0 ,........ '--" --- 1.0- :3 0.5- c» * * * • wz(Q) • ..... ....... 4>+' ~ .•.•• 4> #, + ~ ... ..-" ... ++ .. , + .. , .' ...... . ....... -........ e ......... . ..... - ... . ::::::::::~:::::::::~ ........... . .' * • .. * * • • • * • - .... ~ o 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 Q / A-I FIG. 1. "Dispersion relations" for the reduced second and fourth moments ofthe dynamic structure factor, S( Q,(() at 200 K. Their ideal-gas (high-Q) limits, which are given by Q(kBT 1M) 1/2 and Q(3k BT 1M) 112, respectively, are shown as straight dotted lines. The velocity of propagation for this excitation can be estimated from the frequency of the side maxima of S( Q,w), W B (Q), according t0 28 (12) which in the hydrodynamic limit reduces to the adiabatic sound velocity. Substituting the above-mentioned frequencies, we obtain sound velocities of 1760 and 1600 m/s for Q = 0.25 and 0.35 A-I, respectively. Our simulation results seem to agree reasonably well with recent measurements of the hypersonic velocity.29 In particular, the adiabatic sound velocities reported are 1475 and 1080 m/s for liquid methanol-d4 at 200 and 300 K, respectively. The higher value obtained from our 200 K simulation can be interpreted as a positive sound dispersion in the simulated system. B. Longitudinal current fluctuations In order to investigate thoroughly the dispersion of the Brillouin mode and to clarify whether a second mode exists, we have undertaken the calculation of the spectrum oflongitudinal current fluctuations which is related to S( Q,w) through Eq. (6). Some well-known benefits of studying the spectra of the current ACF are (i) J1 (Q,t) is oscillatory and decays faster than F( Q,t); therefore, the spectra are less affected by truncation or windowing effects; (ii) there is no central peak and consequently no overlapping with the quasielastic component; (iii) the spectra always exhibit, at least, a pair of peaks (Stokes and anti-Stokes); and (iv) the high-frequency part of the spectra is enhanced, allowing the study of weak high-frequency excitations. Constant-Q sections of the spectrum of the longitudinal current fluctuations are shown in Fig. 4. At first sight, the most striking feature is undoubtedly the appearance of a second peak at frequencies around 2.25 X 1013 rad/s, which is the most intense in the interval 0.70 < Q < 1.25 A-I. The Q dependence of the frequency of these two maxima has been plotted in Fig. 5. The low-frequency mode is acoustic in nature and it is related to the peak appearing in S( Q,w) at low Q, though it now appears at somewhat higher frequencies (see Fig. 5) due to the w2 factor in Eq. (6); its "dispersion relation" shows no significant qualitative difference from what has been found in simpler liquids. 28 Namely, it is strongly affected by structural effects at low Q [showing a pronounced dip in the neighborhood of Qp ], but beyond 2.5 A-I kinetic effects become predominant, as can be seen from the comparison with the ideal-gas-limiting value w:::(Q) = Q(2k B T/M) 112. The second mode exhibits a more complex behavior. To begin with, it cannot be unambiguously distinguished as a peak in the whole Q range examined in this study, but usually a shoulder can be observed at frequencies around 2 X 10 13 rad/s in all the spectra not showing a peak. Second, because of the strong overlap between the two modes it is difficult to ascertain whether this higher-frequency mode shows any measurable dispersion. In addition, we found in our previous single-particle study 10 that the center-of-mass velocity autocorrelation function (VACF) has a secondary peak at 2.04 X 1013 rad/s which agrees approximately with the position observed here. The intensity of both modes exhibits important oscillations at low Q that are progressively damped J. Chem. Phys., Vol. 96, No.1 0, 15 May 1992 Alonso et al.: Collective excitations in liquid methanol 7701 0.50 A-I 0.55 A-I 0.5 1.0 1.5 2.0 W / ( 10 13 rad/s ) 0.61 A-I 0.70 A-I 0.74 A-I o78A-I 0.82 A-I 0.86 A-I 2 3 w / ( 1013 rad/s ) FIG. 2. Centers-of-mass dynamic structure factor for liquid methanol at 200 K. S( Q,tiJ) is shown as a function of the angular frequency, tiJ, for the lowest II values of the momentum transfer. Note the different frequency ranges used in both frames. The ordinate units are (10 12 rad/s) -1. M o ...... x o 0.5 0.24 A-I 0.34 A-I 0.42 A-I 0.49 A-I 0.54 A-I 1.0 1.5 2.0 w / ( 10 13 rad/s ) FIG. 3. Centers-of-mass dynamic structure factor for liquid methanol at room temperature. S( Q,tiJ) is shown as a function of the angular frequency, fIJ, for the lowest five values of the momentum transfer. The ordinate units are (10 12 rad/s) - I. 0'0.15 II ..., ~ ...;- _0.10 ~ 2: ...;- 0.05 (a) 0.08 0.07 0'0.06 II do.0 5 ...;- 0.04 - '3 dO.0 3 ...;- 0.02 0.01 (b) 0.25 A-I 0.35 A-I 0.43 A-I 0.55 A-I 0 2 3 w / ( 1013 rad/s ) 0.61 A-I 0.70 A-I 0.74 A-I 0.78 A-I 0.86 A-I ° 234 5 W / ( 1013 rad/s ) FIG. 4. Constant-Q sections of the spectra of the normalized longitudinal current autocorrelation at 200 K as a function of the frequency. Note the different frequency ranges used in both frames. The ordinate units are (10 12 rad/s) -1. as Q increases; however, they oscillate in opposite directions. The intensity of the high-frequency mode is found to be roughly proportional to OJ m (Q) while the acoustic mode follows a trend proportional to its inverse, which is in agreement with its asymptotic behavior at high Q (i.e., ideal gas). A new sound velocity c, (Q) can be defined similarly to cB(Q) but using the frequencies from the maxima of the CACF spectra instead of those of S( Q,OJ). The high- and low-frequency limits of the sound velocity are denoted by COO (Q) and CO (Q), respectively, and are given by COO (Q) = OJ,CQ)/Q, CO CQ) = fYcT(Q), C 13) where r is the quotient of the specific heats and cT(Q) = OJ o (Q)/Q is a wave-vector-dependent generalization of the isothermal sound velocity. Finally, in Fig. 6 we have collected all our data concerning sound velocities. Note that we have plotted CT (Q) instead of the low-frequency J. Chem. Phys., Vol. 96, No. 10,15 May 1992 7702 Alonso et a/.: Collective excitations in liquid methanol "' - 2.~ 11 2.0- .... .,., '0 ..... 1.~ - 0' 1.0- 1 o.~ I~I \ \ }~'!\i+rdl~~H tHt tit· o Q j A-I FIG. 5. Frequency of the maxima of the longitudinal current spectra for methanol at 200 K as a function of the momentum transfer. Open circles represent the acoustic mode while asterisks give the location of the second excitation. The maxima observed in S( Q,Cll) at low Q have been indicated as squares. Finally, the straight dotted line represents the ideal-gas asymptotic limit. Error bars have been estimated as the half width at half maximum (HWHM) of the window function used for Fourier transforming the current ACF. To aid comparison the figure has been plotted on the same scale as Fig. I. sound velocity due to our lack of reliable information on r for the simulated system. The disagreement between C B (Q), C I (Q), and C T (Q) at low Q is a clear symptom of the strong sound dispersion [note that in order to reconcile the values ofc B (Q) and Co (Q) we would have to assumer~2]. At 200 Keven C B (Q) and C I (Q) do not coincide, indicating that the Brillouin peaks are broad and overlap substantially with the central peak. On the other hand, the low-Q behavior of c, (Q) experiences an important change at 300 K, reaching a maximum value and, possibly, beginning to converge towards Co (Q) just beyond the lowest Q we are able to investigate in our simulation (0.25 A-I). 4. \ 300 K 3. 200 K ---. en 3. - S ~ 2.~ -2. C; '(:)1. 1. o. o 0.5 1.0 1.5 2.0 0 0.5 1.0 1.5 2.0 QjA-I Q/A-l FIG. 6. Sound velocities cT(Q), c~ (Q), c,(Q), and cB (Q) as a function of the momentum transfer (see text for definitions). The open circles correspond to c,(Q), the squares to cB(Q), and the solid lines to cT(Q) and c ~ (Q), as indicated. c. Transverse current fluctuations We have undertaken the calculation of transverse CACF in order to get some insight into a part of the collective dynamics not easily amenable to experimentation. Experimental techniques such as neutron scattering do not couple easily to transverse modes (unless some kind of diffuse "umklapp" process is postulated). In fact, to date molecular dynamics is the most reliable source of information about transverse excitations in dense fluids for the kinematic range of our interest. The frequency of the maximum of the transverse CACF spectra has been plotted in Fig. 7 as a function of Q. A fit of the low-Q data of the "dispersion relation" to a straight line gives an estimate of the velocity of propagation which is 940 ± 20 mls at 200 K. Corresponding spectra and dispersion relation at 300 K look qualitatively the same. However, at this temperature there is a substantial overlap of the Stokes and anti-Stokes peaks for the lowest Q which distorts to some extent the dispersion relation; for example, at the lowest Q thevalueofJ, (Q,liJ = 0) is about 1/2 of the value at the maximum, while at 200 K it is 1/20 of it. In addition, the velocity of propagation has now reduced to 750 ± 25 m/s. The lifetime of the excitation 30 can be estimated from the reciprocal ofthe half width at half maximum ofthe peaks. At 200 K we find lifetimes of 1.06, 0.71, 0.42, and 0.35 ps at Q = 0.25, 0.35, 0.43, and 0.50 A - I, respectively, while at 300 K the lifetimes are 0.64,0.49,0.35, and 0.30 ps for similar values of momentum transfer. An extrapolation of the dispersion relation to lower Q, at 300 K, suggests that the peaks at ± liJ m (Q) will coalesce into a single one for Q < 0.1 A-I, but there is no indication of that happening for any finite Qat 200 K. Undoubtedly, the reason why methanol seems to be able to support freely propagating shear waves with wavelengths as large as 50-100 A is the existence of a very stable hydrogen-bond network, which is itself connected to the strong directionality of the hydrogen-bond interaction. ---.0. en - ~ ~ o. o ..... cr '1 o. o 0.5 1.0 1.5 2.0 Q / A-I FIG. 7. Frequency of the maximum of the transverse current autocorrelation spectra for methanol at 200 K as a function of the momentum transfer. Error bars as in Fig. 5. J. Chem. Phys., Vol. 96, No.1 0, 15 May 1992 Alonso ef al.: Collective excitations in liquid methanol 7703 For momentum transfers greater than 1.0 A-I the transverse CACF behaves mostly as a single-particle coordinate. Finally, to end the analysis of the spectra of the transverse current we just mention that a conspicuous secondary maxima can be observed at the same frequency as the third (very weak) maxima in the VACFspectra, 10 i.e., 3.57X 1013 rad/s, in almost all the spectra with Q between 0.5 and 1.5 A-I, though it is extremely weak. In the low-Q region the spectra contain some ripples that preclude the observation of weak signals as this secondary peak. D. lattice-dynamics calculations In order to explore the physical origin of the modes appearing in the simulation as well as to have a reference to compare with, we have undertaken a lattice-dynamics (LD) calculation for a polycrystal sample following the procedure described in the Appendix. From the calculated dispersion relations of the 24 modes corresponding to rotations and translations of the four molecules within the unit cell, a quantity directly comparable with the current correlation functions has been derived. Several plots corresponding to different values of the momentum transfer are shown in Fig. 8. From inspection of the figure the following comments are in order: (a) The manifold of acoustic peaks appears as a well-defined entity up to Q values of about 0.5 A, showing noticeable dispersion; and (b) two series of peaks of a mixed (rotational and translational) character are clearly apparent from 0.5 A onwards, centered at about 2.2 and 3.0X 1013 rad/s. Such frequencies can be easily correlated with prominent features appearing in the calculated vibrational density of states (DOS) which is displayed in Fig. 9. A close inspection of the figure where a comparison of atom-atom DOS is made for the crystal and liquid phases reveals the following features: • The purely acoustic modes are apparent in the polycrystal as a Debye contribution (Le., with a dependence upon (2) up to frequencies of 5 X 10 12 rad/s. The first two intense peaks centered around lOX 10 12 rad/s are of mixed (rotational and translational) character as has been evidenced from the analysis of the mode eigenvectors. • The polycrystalline DOS up to w = 15X 10 12 rad/s becomes in the liquid phase a broad, structureless rotational contribution centered at about the same frequency as in the solid. The sound mode (translational or center-of-mass) contribution appears in the liquid as a broad distribution centered at about 5 X 10 12 rad/s, which will give rise at low-Q values to a finite-frequency response in the S( Q,w). A noticeable amount of translation-rotation coupling is apparent, which evidences the fact that no pure translational or rotational modes are to be expected to occur in a molecular material within the kinematic range accessible to MD or INS studies. • A series of intense peaks centered at 22, 30, and 38 X 10 12 rad/s in the polycrystal appear as a broad structure between w = 15x 10 12 and w = 30X 10 12 rad/s in the liquid in a frequency region showing also a noticeable amount of rotationtranslation coupling. The intense peak which appears in the solid at about 50 X 10 12 rad/s also shows up in the liquid, although the modes contributing to it are now softened, which leads to a reduction in the peak frequency to about 37 X 10 12 rad/s. From the analysis of some of the mode eigenvectors the components of the mode polarizations can be obtained. As an example, it was found that the two low-frequency peaks correspond to modes with an average character of Tx(y,z) = 0.20 (0.23, 0.36) and RX(Y,z) = 0.74 (0.47, 0.11) for the low-frequency peak, and Tx(y,z) = 0.21 (0.22,0.67) and Rx(y,z) = 0.24 (0.62,0.14) for the higher-frequency one. Here Tx(y,z) and Rx(y,z) denote the averaged amplitUdes of the translational and rotational molecular degrees of freedom expressed in the principal axis of the inertia molecular frame (in increasing inertia moment order). Therefore, it seems clear that the comparison of the current-current correlations as well as the DOS for both polycrystal and liquid phases reveals that most of the observed features in the fluid phase can be assigned by taking the LD results as a reference. Further discussion on this will be deferred to the Discussion section. V.INELASTIC NEUTRON SCATTERING AND MODELING OF S{Q,ro) Recently, Bermejo et aU have measured the inelastic neutron-scattering spectra of liquid methanol-d4, I( Q,w), which is related to the atom-atom dynamic structure factors no U t no no I(Q,w) ex: I ~S~lf(Q,W) + I I bibjSi/Q,w) , i=1 41T i=lj=1 (14) where na is the number of atoms per molecule, (J inc and bi are the incoherent scattering cross section and the coherent scattering length for neutrons of atom type i, and S ij ( Q,w) is the atom-atom dynamic structure factor for atoms i and j (the superscript self makes reference to its self-part). However rigorous, the preceding equation would not be very useful as a model due to its large number offree parameters and a simplified one based in a well-known viscoelastic approximation for simple liquids has been proposed7 I model ( Q,w) ex: Scm (Q) exp( -p,Q 2) X {Rqe (Q,W) + P6 (pi - P6 )7 } [W7(W2 -PF)] 2 + (W 2 -P6)2 ® O(w), (15) where Scm (Q) denotes the structure factor for molecular centers, the exponential is a Debye-Waller term, Rqe (Q,w) represents the central quasi elastic component, the second term inside the curly brackets is the usual viscoelastic scattering law for simple liquids,9 O(w) is the experimental resolution function; and ® denotes a convolution operation. To keep the number of free parameters to a minimum the Maxwellian relaxation time 7 is estimated using the prescription due to Lovesey,9 7 -1=2[(PF-P6)!1r]1I2. (16) As usual, the Q dependence of the reduced moments from J. Chem. Phys., Vol. 96, No.1 0, 15 May 1992