Using the PHONON program for interpreting inelastic neutron scattering data

Nuclear Instruments and Methods in Physics Research A 354 (1995)66-72 &MElwoos IN PHYSICS

Using the PHONON program for interpreting inelastic neutron scattering data J-c. Lia,*, M. Leslieb aDepartment of Physics, University of Sal&ord, Saljord, M6 4wT, UK b Daresbury Laboratory Warrington, WA4 4AD. UK

Abstract The PHONON program was developed into a generally useful form for lattice dynamic calculations using the quasi-harmonic approximation and made available to other users interested in the interpretation of INS data. By diagonalizing the dynamical matrix, the actual values of the o(q) for each wave vector q are determined. The dispersion curves can, therefore, be obtained and plotted using NAG graphic routines on a suitable computer. The weighted phonon density of states can also be obtained by integration over all the possible values of q lying in the first Brillouin zone of the crystal. There is a wide range of two-body potentials programmed into the code which are used for the short-range repulsive parts of the potential. Ionic polarisation is modelled using the shell model. There are also a number of many-body potentials for molecular systems and lattice structures. These include angle bending three-body potentials and four-body (torsional) potentials.

1. Introduction Thermal (or subthermal) neutrons are an excellent probe for studies of vibrational dynamics for a given solid and liquid, because several remarkable properties make thermal neutrons unique: for instance, the thermal neutron energy is comparable to phonon energies and the wavelength associated with the neutron is of the same order as the interatomic distances in condensed phases. Another characteristic of this probe is that the neutron mass is of the same order as the mass of the scattering nuclei. The scattering is, therefore, sensitive to the structure and dynamics of the system. In order to obtain physical significance from the quantities measured using neutron scattering technique, the spectrum has to be reproduced based on a particular lattice dynamic treatment of the system (the molecular dynamics approach is an another way). In the lattice dynamical approach, the role of frequency distribution functions can be discussed systematically. The calculation begins with the adiabatic approximation which enables us to treat the solutions in the electronic problem as interaction potentials in the nuclear problem. By assuming harmonic forces and periodic boundary conditions, we can obtain a normal mode distribution

function of the nuclear displacements. Hence the problem then reduces to a classical system of coupled oscillators. Under these conditions the one-phonon dispersion relation can be evaluated and hence the measured scattering intensities can be rigorously reconstructed. The lattice dynamics program PHONON is made just for this purpose based on the potential model.’ The method used to calculate phonon dispersion curves from a potential model has been described in the literature [l], so only brief details will be given here. A potential model describes the potential energy of a crystal. It is assumed that the internal coordinates (basis atom positions) of the crystal are at their fully relaxed positions. In the harmonic approximation, the terms greater than second order in the expanded potential energy in the displacements of atoms in the crystal are neglected. The force constants are: PV(r)

F(r) = K.

1

(1) J

In the quasi-harmonic approximation, the cell constants are not relaxed using the potential model, although the internal coordinates are. This gives force constants which are temperature or pressure dependent. The equations of

*Corresponding author. 016%9002/95/$09.500 1995 Elsevier Science B.V. All rights reserved SSDl 0168-9002(94)01035-8

1The manual is available on request from Dr M. Leslie.

J-c. Li, kf. Leslie/Nucl.

Instr. and Meth. in Phys. Rex A 3.54 (1995) 66-72

61

motion are constructed from this model and plane wave solutions are sought. The frequencies of vibrational modes at a given q can be determined by diagonalising the dynamical matrix D(q) of the crystal:

(2) where C(q) is the matrix of the wavevector-dependent mass-weighted eigenvectors and o(q) is the wave vector-dependent eigenvalue. The dynamic matrix D(q) can be constructed using F(q) and B(q) matrices (D(q) = M-1’2B*(q)F(q)B(q)M-1’2 (M is the diagonal matrix of the atomic mass). B(q) is 3n x 3n x m dimensional matrix, where n is the number of elements in the unit cell and m is the number of force parameters. Therefore, the larger unit cell means that a larger matrix must be diagonalised at each perturbation cycle for each value of q. Hence the measured dispersion curves in the (Q, w) space of the instrument can be compared with the calculation results of the o(q). On the other hand, in many cases, measuring the dispersion relations are practically impossible. Alternatively, the frequency distribution function, weighted by the amplitudes of modes and scattering cross-section, G(w), can be measured. This quantity can also be calculated by integration of all the possible frequencies of the vibrational modes with wave vector q lying in the first Brillouin zone of the crystal,

G(W)= C Oinc C (C,S*Q)*exp(- (u~)Q’)/(~N/ww). 1

S

(3) Because the scattering intensity observed from the sample, proportional to the double differential cross-section per unit neutron energy loss, can be converted to the weighted phonon density of state, G(w) =

(4)

therefore, the measured and calculated G(w) can be compared at this point. Conventionally, the G(w) are calculated by integrating all the frequencies for each q value outside the lattice dynamics program. High-quality calculations require fine mesh points of the q-vector in the three-dimensional reciprocal lattice, resulting in the calculations being very time consuming. In the PHONON program, the G(w) can be calculated internally by defining the first BZ boundary. Equally the full set of dispersion cures in the principal directions can be plotted using NAG graphic routines. In addition, there are many available options in the program. Ionic charges may be used, in which case the sums are carried out using the Eward method. It is also possible to model ionic polarization using the shell model. This effectively introduces massless

30 20 IO 0

A

K

0

M

0

9

Fig. 1. Dispersion curves for single-crystal D20 ice Ih measured by Renker’s using IN8 at ILL [3] and fitted by his five force constants model.

particles into the equations of motion. Furthermore, three-body potentials include angle bending forces and four-body interactions involving torsions are also in the program.

2. Lattice dynamics of ice Ih Although, for many years, ice has been one of the mos‘ widely studied molecular solids, a complete understanding of its peculiar behaviour has not been achieved at the microscopic level. A number of incoherent [2] and coherent [3] inelastic neutron scattering measurements have been carried out previously, but for various reasons, and in particular because of limitations of instrument performance, little real progress has been made towards understanding the lattice dynamics of even the common form of ice, ice Ih. The dispersion curve measurements of Renker [3] on ice Ih (D20) made by time-of-flight using a chopper spectrometer were extensive but not complete in that information above 20 meV in the (00 1) direction and above 30 meV in other directions is missing (see Fig. l), presumably because of a lack of scattering intensity and of the nature of the disordering of the protons in the ice structure. This information is crucial for clarifying the lattice models proposed by Renker [3] and others [4-lo], in particular with reference to the splitting between the peaks at _ 28 meV and N 37 meV for Hz0 as observed by inelastic incoherent neutron scattering (IINS) measurements [ 111, i.e. for determining whether the splitting is due to a greater stretching force constant along the c-axis of the hexagonal ice or to some other cause. Recently, in connection with our long-term programme of studies of lattice dynamics of ice [ 121, we have made measurements on a single crystal of ice Ih (D20) on PRISMA (Coherent Collective Excitation Spectrometer) at ISIS [13], in an attempt to reach a better under-

NEUTRON SCATTERING DATA ANALYSIS

68

J-c. Li, M. Leslie/Nucl. Imtr. and Meth. in Phys. Res. A 354 (1995) 66-72

1

1.5

2

Wavevector

0

2.5

3

3.5

Fig. 2. Dispersion curves for single-crystal D20 ice Ih measured using the PRISMA at RAL in the crystal orientation of (110).

standing of the high-energy molecular optic dispersion curves in the region between 20 and 40 meV and the hindered rotational modes above 40 meV. This instrument is an inverted geometry, time-of-flight spectrometer which provides measurements of excitation spectra from single-crystal samples over a wide range of (Q, o) space concurrently (energy transfer from a few meV to a few hund:eds meV, momentum transfer from 0.1 to 10.0 A-‘) with good resolution [14]. In order to facilitate measurements of the well-defined optic branches between 20 and 40 meV, high-quality single-crystal ice Ih (DzO) were made in the Ice laboratory, University of Birmingham for this particular experiment. The technique used for growing the single-crystal sample can be found elsewhere [15]. Fig. 2 shows one of the particular measurement along the direction of (100). The results show that there are well-defined acoustic modes which agree with both ours and Renker’s model (see Fig. 1). However, for the molecular optic modes at higher energy transfers, there are no well-defined dispersion curves observed, apart from a particular branch which increases in energy from the BZ boundary (30meV) to the zone centre (36 meV). There are two flat and broad features which cross the whole BZ at -28 and -37 meV as observed in the inelastic incoherent spectrum for DzO [ 111. Similar behaviour is also seen in other crystal orientations. This indicates that the optic modes are completely localised, as there is little dispersion. Theoretical representations of the vibrational spectrum of ice Ih have been proposed by various authors [3-lo]. However, a final and unambiguous explanation of the IR, Raman and Neutron data has still not yet been achieved, mainly because of the limitation in our understanding of the unusual physical phenomena associated with the intrinsic disorder of the H positions in ice Ih. The importance of the dipole effects arising from this

disorder has been recognised in a number of papers [4,5]. However, these effects have not previously been properly included in lattice dynamical calculation for ice, partly due to the fact that there has been no accurate experimental data on the system. Renker’s calculation [S], which was based on an ordered proton model for ice Ih, used five force constants, where the O-O stretching force constant along the c-axis is about 1.7 times larger than the stretching force constant in the &plane. It predicts that the peak in the PDOS integrated over all the vibrational amplitudes along the c-axis gives only the peak at 37 meV (298 cm-‘) while the PDOS integrated over all the vibrational amplitudes in the &-plane has a peak at only 28 meV (225 cm-‘). Thus, this model is incompatible in these two orientations with our IINS data for single crystalline ice Ih which showed that the peaks at 28-37 meV do not depend on the orientation of the crystal. A number of alternative explanations have also been proposed. Faure et al. [4] suggested that there could be a partial ordering of the O-H bonds parallel to the c-axis of ice Ih due to long-range forces, without breaking Bernal-Fowler’s ice rules [16], but offered no detailed examination. Wong and Whalley [S] suggested that the H-bond stretching force constant depends on the relative orientation of the two nearest water molecules. The differing configurations of the bond asymmetries cause significant differences in the forces which determine the vibrations of molecules parallel and normal to their polar axis. The interaction will tend to destroy translational symmetry and could be important in determining the forms of the normal vibrations and their IR and Raman activity. However, these authors did not reproduce successfully the peak at 35 meV in the IR spectrum. In 1986 Marchi et al. [6] have reported results of lattice dynamic calculations on ice Ih using a periodically replicated disordered sample of 128 water molecules using an effective pair potential for the intermolecular interactions. The results obtained for the phonon density of states at the Brillouin zone (BZ) centre only, still do not agree well with the IINS data, because the optic modes at about 30 meV is not sufficiently split and, hence, the peak at 37 meV is missing. In order to explain the two optic peaks in IINS spectrum, Klug and Whalley [7] suggested that the two peaks are due to TO and LO splitting, as is seen in ionic crystals such as NaCl, but give no detailed theoretical calculation for what the neutron spectrum would look like based on this postulate. Although there were some disagreements as to whether the idea of TO and LO splitting is valid for ice, since the argument for it depends on the existence of long-range order, it has been generally accepted for the explanation of the weak shoulder at 35 meV in IR and Raman spectra. Lately it has been extended to explain the two molecular optic peaks at 28 and 37 meV seen in the IINS spectrum [17]. There are, however, a number of strong arguments against this

J-c. Li, M. Leslie/Nucl.

69

Ins@. and Meth. in Phys. Rex A 354 (1995) 66-72

Table 1 Force constants

0 (unit mdyn/A)

O-O harmonic force constant Three-body bending force constant Kz’ A type D-D interaction K,ab’ D type K,” B type K,“bZ C type

0.1782 0.0421 0.1620 0.1620 0.0000 0.0000

K

G

&L (C)

CD)

Fig. 3. The upper diagram shows a section of the Ice I network with a typical arrangement of strong and weak bonds; the strong H-bonds are indicated by the solid bonds and the weak bonds are shown as double lines: The lower diagram shows four possible orientations of molecule pairs in ice Ih. In ice Ic only the (C, strong bond) and (D, weak bond) types of molecule pairs are found.

postulation that have emerged as result of our recent measurements on a range of different ice structures using IINS, including all “recovered” high-pressure phases of ice. These studies show that the two peaks at 28 and 37 meV do not necessarily appear at the same time in the spectrum for different phases of ice [12]. In order to explain all the properties of these spectra of ice as a whole, we have to conclude that the two molecular peaks are associated with different dipole-dipole (D-D) configurations 1181 as shown in Fig. 3. We postulate that the relative intensities of the two optical bands are entirely dependent on the relative number (or ratio) of the configurations which are related to the strong and weak H-bonds in the ice structure. For instance, in ice Ic (which has an identical spectrum to ice Ih), the protons are completely disordered. Hence, statistically, there is one configuration D to every two configuration C. Therefore, in ice Ic, there will be one weak H-bond for every two strong H-bonds (ratio is l/3:2/3). They are randomly and isotropically distributed in the ice structure which gives the integrated peak area for the low/high-energy optic modes in a ratio of about 1: 2, as observed. Furthermore, the diEerence of the two H-bond force constants are considerably larger than can be explained on the basis of electrostatic effects only. A simple calculation using the relationship of w N Jk/m (where w is the vibrational frequency and k is the force constant) gives an estimated ratio of the two H-bond force con-

stants as kl : k2 = (ml : o_I~)~= 282 : 372 = 1: 1.8 (the exact values are determined by the lattice dynamic calculation). The molecular orbitals are influenced strongly by correlations in the orientation of neighbouring molecules. The actual source of the splitting of the H-bond, therefore, may lie at the quantum mechanical level [lS].

3. Calculations of phonon density of states As we have seen in the Section 2, the determination of the normal mode frequencies within the harmonic or quasi-harmonic approximation requires at the outside a knowledge of the equilibrium atomic positions, i.e. the crystal structure, the atomic masses and the interatomic potentials (or force constants). The next stage is to make any possible simplification arising from the symmetry conditions or other restrictions on the interactions that may be imposed. In the case of ice Ih, in order to represent the disordering of the proton arrangements, a large “superlattice” (cluster of water molecules) must be constructed based on the ice rules [16]. For such models, dispersion curves lose their meaning and only an overall amplitude weighted phonon density of states can be calculated. It is clear that a larger cluster may be closer to disordered ice Ih. However, this will be difficult to achieve because (a) a larger cluster means that a whole set of possible arrangements needs to be studied; (b) computational requirements will increase considerably because the dynamical matrix becomes very large. We have therefore constructed a cluster containing eight primary unit cells (32 molecules). The force constants used for this calculation are listed in Table 1. The PDOS calculation involves full integration of the first BZ using standard “root sampling method” [19]. Fig. 4 illustrates the G(o) obtained from IINS data and from the calculation. As can be seen, the main features of the calculated G(w) agree well with the experimental one; firstly, the G(w) calculated from the model reproduces the acoustic modes very well - see for example, the sharp acoustic peaks at 7 meV for G(w) as observed in the experimental data. Secondly, the dramatic improvement in these calculations is in the optic modes between 15 and

NEUTRON SCATTERING DATA ANALYSIS

70

J-c. Li, M. Leslie/Nucl. Instr. and Meth. in Phys. Rex A 354 (1995) 66-72

0 0.5 A

‘;”

0(1/A)

1.0

1 .o

K

M

0.0 I

B 3 30

0

10

20 ENERGY

30

40

50

(meV)

Fig. 4. The weighted phonon density of states G(w): (a) measured, (b) calculated.

40 meV and, not only in the peak positions, but also in the general triangular shapes of the two peaks at 28 and 37 meV. These features result directly from the randomness of the strong and weak force constants and a cluster larger than one unit cell. The right-hand diagram shows dispersion curves for ice Ih calculated using the model, which show two flat regions at 28 and 37 meV, giving the sharp higher-energy cut-offs. The strong D-D interactions push some dispersion curves above 30 meV and the large cluster for representation of the disordering breaks down the translation symmetries and the localised modes exhibit triangle shapes with high-energy cut-offs. Although the level of noise in the calculated spectrum is rather high, this could be improved by increasing the number of mesh points in the calculation (40 points along each direction have been used for this calculation) and by convoluting the result with the resolution function of the instrument used.

4. Calculated dispersion curves Because the strong and weak stretching H-bonds are randomly located in the ice crystal, it is possible that in the high-energy transfer region there are no well-defined dispersion curves. In the large “superlattice” model, dispersion curves cover the whole region from 25 to 40 meV and the density of the curves increases from the left to the right for each of the optic bands at 28 and 37 meV for Hz0 ice. To clarify this problem, we chose a few unit cells in the large “superlattice” cluster and regard these as primary unit cells, namely Cl and C2, and calculate their dispersion curves [12]. Fig. 5 illustrates the dispersion

10

Fig. 5. Calculated dispersion curves using the PHONON program. The upper curves are calculated using the Cl model and the lower curves are calculated using the C2 model. curves for the two cells. As can be seen, initially in the lower-energy transfer region ( < 20 meV), the two sets of dispersion curves are almost identical to Renker’s in Fig. 1, because the acoustic modes are insensitive to the stretching force constants. There are at least four branches close at 28 meV in Renker’s model and two higher energies which cannot produce enough intensity in the G(o) to fit the IINS experimental data [ll]. In our model, the optic branches are separated equally into two regions in the Cl model and show more branches above 30 meV in the C2 model, having very flat dispersion curves at -37 meV, thus producing a sharp cut-off. Detailed assignments of the modes using group theory were well-established by Forsind [20] for ice Ih. The six vibrational modes Ai,, Bz,, El, and Elg for the ordered Renker model, and the disordered Cl and C2 structures are illustrated in Table 2. Because there are asymmetrical forces on all atoms, the degenerate modes El, and Ezg are separated and unable to vibrate purely in either the (10 0) direction or the (1 10) direction and Al,, Bzu in the (00 1) direction as the ordered cases in Renker model as shown in Table 2. Hence assignments are rather difficult for them. For both cells of ice Ih, a numerical investigation of the molecule displacement eigenvectors using the PHONON program shows that in the case of the low-energy bands, and in particular the three acoustic modes Al, and EL,, all the

71

J-c. Li, M. Leslie/Nucl. Instr. and Meth. in Phys. Res. A 354 (1995) 66-72 Table 2 Schematic illustration of the vibrational modes for the ordered model and the disordered models

Modes

for

ordered ice Ih

0t

I OT

OT

28.1

I

t0

LO

28.1

28.1

lo

10

YO

28.1

I

10’

TO’ CO

Energy (meV)

01

36.8

I

36.8 01

Modes for

I

0t

Cl model TO’

I

10

Energy (meV)

28.2

37.8

28.2

37.8

37.6

27.5

28.4

27.2

37.8

38.0

32.5

35.1

Modes for C2 model

Energy (meV)

water molecules in the unit cell oscillate as a whole at the BZ centre. For other optic modes, atomic vibrations could have all components in the three directions, for instance, the Ei, and EZg equivalent branches as shown in Table 2. However, for the two vertical vibrational modes, Ai, and Bzur the displacements are close to the c direction, except that the amplitudes vary from one atom to another in the unit cell. This may cause a change of the dipole moments, not only to first order but also to second order, which may indicate the reason for both IR and Raman activity. Acknowledgements

The authors would like to thank the Engineering and Physical Sciences Research Council (UK) for financial support and the Rutherford Appleton Laboratory for the use of neutron facilities.

References [l] G. Dolling, in: Methods in Computational Physics, Vol 15, ed. G. Gilat (Academic Press, New York, 1976) p. 1. [Z] H.J. Prask, S.F. Trevino, J.D. Gault and K.W. Logan, J. Chem. Phys. 56 (1972) 3217. [3] B. Renker, in: Physics and Chemistry of Ice, ed. E. Whalley, S.J. Hones and L.W. Gold (University of Toroto Press, Toroto, 1973), p. 82. [4] P. Faure and A. Chosson, in: Light Scattering in Solids, ed. Balkanski, M. (Flammarion Sciences, Paris, 1971), p. 272. c51P.T.T. Wong and E. Whalley, J. Chem. Phys. 65 (1976) 829. PI M. Marchi, J.S. Tse and M.L. Klein, J. Chem. Phys. 85 (1986) 2414. II71D.D. Klug and E. Whalley, J. Glociology 85 (1978) 55. PI R.E. Shawyer and P. Dean, J. Phys. C 5 (1972) 2028. PI J.R. Scherer and R.G. Snyder, J. Chem. Phys. 67 (1977) 4794.

NEUTRON SCATTERING DATA ANALYSIS

12

J-c. Li, M. Leslie/Nucl. Instr. and Meth. in Phys. Res. A 354 (1995) 66-72

[lo] E. Whalley, in: Physics and Chemistry of Ice, eds. E. Whalley, S.J. Hones and L.W. Gold (1978), p. 13. [1 1] J-c. Li, D.K. Ross, L. Howe and T. Tomkinson, J. Physica B 156 & 157 (1989) 376. [12] J-c. Li and D.K. Ross, in: Physics and Chemistry of Ice, eds. N. Maeno and T. Hondoh (Hokkaida University press, Hokkaida, 1992), p. 27. [13] U. Steigenberger, M. Hagen, R. Caciuffo, C. Petrillo, F. Cilloco and F. Sacchetti, Nucl. Instrum. Methods B 53 (1991) 87. [14] M. Hagen and U. Steigenberger, Nucl. Instrum. Methods B 53 (1991) 87.

[15] M. Ohtomo, S. Ahmad and R.W. Whitworth, J. de Physique Colloque Cl, Suppl. au No 3, Tome 48 (1987) Cl595. [16] J.D. Bernal and R.H. Fowler, J. Chem. Phys. 1(1933) 515. [17] D.D. Klug, J.S. Tse and E. Whalley, in: Physics and Chemistry of Ice, eds. N. Maeno and T. Hondoh (Hokkaida University press, Hokkaida, 1992), p. 88. [lS] J-c. Li and D.K. Ross, Nature 365 (1993) 327. [19] G. Gilat, in: Methods in Computational Physics, Vol 15, ed. G. Gilat (Academic Press, New York, 1976), p. 317. [20] E. Forslind, Swed, Cement Concrete Inst., Roy. Inst. Technol., Stockholm, Proc., No. 21 (1954) p. 1.