Calculation on sulphur isotope fractionation between sphalerite and galena using lattice dynamics

Earth and Planetary Science Letters, 28 (1975) 172 - 180

© Elsevier Scientific Publishing Company, Amsterdam

Printed in The Netherlands

kZ3

CALCULATION OF SULPHUR ISOTOPE FRACTIONATION BETWEEN SPHALERITE AND GALENA USING LATTICE DYNAMICS* M.M. ELCOMBE Australian Atomic Energy Commission, Sutherland, N.S. W. (Australiaj

and J.R. HULSTON Institute o f Nuclear Sciences, D.S.1.R., Lower Hurt {New Zealandj

Received November 1972 Revised versions received February 1975 and September 23, 1975

Increasing use is being made of sulphide minerals in isotope geothermometry. Sulphur isotope fractionation factors for 34S exchange between sphalerite (ZnS) and galena (PbS) have been calculated for various temperatures between 0 and 1000°C. The reduced partition function ratios have been calculated using a "shell" model for the forces derived from inelastic neutron scattering studies of the lattice dynamics of sphalerite and galena. A new formalism of the current theory has been developed which enables the accuracy of the calculation to be determined. The Zn34S-pb34S equilibrium constant obtained by the shell model calculations is 1.0060 to 100°C, 1.0031 at 250°C and 1.0005 at 1000°C, in agreement with experimental determinations.

1. Introduction

Salomons [ 13] (FeS2 - P b S -S); and Czamanske and Rye [14] ( Z n S - P b S ) . To use this technique o f geothermometry effectively it is necessary to establish temperature scales for the isotope equilibrium exchange reactions between coexisting mineral pairs. Bottinga [15] has reported the calculation o f the 13 C fractionation factors for the CO2 diamond and d i a m o n d - g r a p h i t e systems using lattice dynamical force models of diamond and graphite; but no experimental data were available for comparison. Three independent sets of experimental measurements o f the 34S fractionation factor of the Z n S - P b S pair are available [9,12,14]; and calculation of this fractionation factor should show whether the full lattice dynamical treatment represents an improvement over the earlier Debye-Einstein model calculations. To date no error estimates of calculated fractionation factors have been published. Using a recent result of Elcombe [16], an alternative formalism is presented (section 2) which enables an error estimate to

The use of the fractionation of stable light isotopes to estimate the temperature of formation of minerals in geochemical systems, first suggested by Urey [1,2], has now become well established. Although the application of sulphur isotope fractionation to geothermometry was pointed out by Tudge and Thode [3] and Sakai [4], its potential for sulphides was not fully realised until comparatively recently when a number of workers [5 8] noted that in naturally occurring coexisting sulphides the 34S content consistently increased in the order galena -- sphalerite - pyrite. Experimental measurements on various sulphide systems have now been reported: Grootenboer and Schwarcz [9] (ZnS P b S - F e S E - S ) ; Kajiwara et al. [10] (FeSz Z n S - P b S ) ; Rye and Czamanske [1 1] (ZnS PbS); Kajiwara and Krouse [12] (FeS: C u F e S 2 - P b S - Z n S ) ;

* INS Contribution No. 484. 172

173

be made for the first time. The 34S fractionation factor of the sphalerite-galena pair, including an estimate of error, has been calculated using a shell model [21] for the interatomic forces in the crystal (section 3). Considerable discussion is given to the numerical and theoretical errors occurring in calculations of reduced partition function ratios and fractionation factors (section 4).

and m* are the isotopic masses of the sulphur atoms. The relation between the isotopic frequency ratios and the isotopic mass ratio for a crystal:

2. Theory

P = ln(Q*/Q),,,

The 34S exchange between sphalerite and galena can be represented by the equation: Zn34S + Pb3’S + Zn32S t Pb34S

(1)

The fractionation factor for 34S exchange between sphalerite and galena, defined as: (Y= (3”S/32S)ZnS/(34S/32S)PbS is identical to the equilibrium which is given by: a=

(2) constant

of eq. 1 [ 151

@*/Q)z,sl(Q*/Q)pbs

(3)

is the partition function ratio where @*/Q)z,,s~PI,s) of 34S sphalerite (galena) to 32S sphalerite (galena). For solids, only the partition function of the vibrational modes need be considered [ 171. In the harmonic approximation for a finite crystal:

In Qtib = k ‘g

ln {exp(-iui)/

[1 - exP(-ui)]

1

(4)

where Qvib = partition function per unit cell, i labels the normal mode of vibration, u = hv/ks T, v is the frequency of the normal mode (all frequencies in this paper are considered as zero-order frequencies), h is Planck’s constant, k, is Boltzmann’s constant, T is the absolute temperature, IV is the number of unit cells in the crystal, and r is the number of atoms per unit cell. Both ZnS and PbS have face-centred cubic structures with one formula unit per primitive unit cell (i.e. r = 2). The calculations per primitive unit cell are therefore equivalent to the molar equations (eqs. l-3). For convenience, the logarithm of the reduced partition function ratio fi, normally used in calculations for gaseous compounds, and given by P = ln(Q*/Q)red = ln(Q*/Q) - $ln(m*/m), has been calculated, where m

for a single atom substitution k [ 161, may be used to show that /3 for a crystal is also given by:

= ln(Q*/Q)

+h‘gln($/vJ

(5)

By using different approximations several formulae are available for calculating In Qvib (eq. 4). One way is first to calculate the frequency distribution g(v) and then to calculate the partition function using: *m

=knGl gnIn Cexp(-$,)l [1 ln Qvib

exp(-u,)]

> (6)

where g,, = ~‘2” g(v)du is the number of normal modes in the interval u,, to v~+~ and is normalised by means of the equation Jim g(:(v)dv = 3Nr. v, is the maximum frequency and ~1, is the number of frequency intervals into which g(v) has been divided. In this case only a knowledge of g(y) is needed to calculate the partition functions and hence determine cr (eq. 3). The accuracy of this method depends on the number of independent wavevectors (i.e. those in the irreducible volume of the first Brillouin zone) used in determiningg(u), and on the frequency interval -u,) used. (%+1 However, the frequency distribution need not be determined explicitly. Replacing the number of unit cells in the crystal (A’) by a sum over a uniform distribution of wavevectors (4) in the first Brillouin zone, and using the symmetry of the lattice to restrict calculations to the irreducible volume of the first Brillouin zone, gives:

ln Qvt, =

C4 dq $ ln{exp(-fuiq)/j[l

- exp(-rr@)j}/q’l, (7)

where d, is the degeneracy of the wavevector q and

174

uiq is the frequency of the ith normal mode with wavevector q. In this case the two steps of the first method become one. The accuracy of eq. 7 will depend only on the number of independent wavevectors used. For both methods a model of the crystal is required, which is assumed not to change with isotopic substitution. The calculations are done twice, once for the normal material and once for the isotopic material; In Qv~ and In Q*b are determined separately and the difference taken to obtain the partition function ratio. This is a lengthy procedure, the final answer, being the difference of two large numbers, is subject to rounding errors, and the error in the answer is very difficult to assess• So far eq. 7 has little advantage over eq. 6 but it can be differentiated with respect to Uiq, whereas eq. 6 cannot. Coupling differentiation with the result of Elcombe [ l 6], which relates the isotopic frequency shifts to the normal mode eigenvectors, it is possible to obtain expressions for ln(Q*/Q) and ln(Q*/Q)re d directly• For a small mass shift, Am = rn*--m:

*2

2

~iq -- Uiq 2

Uiq

*

mk --ink z for all q " ?~1~¢ " eiqk

The eigenvectors are normalised by:

ka

iqkc~ ejqkc~

and ebk (=e~qk .eiqk)is the squared eigenvector magnitude of atom k in the ith normal mode of wavevector q. For brevity in the remainder of the paper e2iqk is referred to as the eigenvector. Rewriting this and expressing Uiq in terms ofu~/ give:

dUiq _ Uiq . . . . z dm 2m* "eiqk

3r( q

ln(Q*/Q) = 3r q

dq

"=

+ exp(Uiq ) - 1

dram A

dq (9)

Using eq. 5, the approximation ln(u*/u) ~- Au/u (for small Au) and changing the variable to u give:

fi = ln(Q* /Q)red =

dq

q

.=

exp(uiq ) -- 1

~

(12)

The introduction of e~qk is no problem. In our computer program to determine the normal mode frequencies of the crystal the eigenvectors can readily be determined simultaneously, with only a slight increase in computing time. In fact, the overall computing time is reduced, as the partition flmction ratio depends on the frequencies and eigenvectors of only one composition. Rounding errors, from taking the difference of two large numbers, have also been eliminated. An approximate form ofeq. 12, valid at high temperatures, is very useful for geothermal calculations. For large T, Uiq is small and a series expansion of eq. 12 gives:

ln(Q*/Q)= ~ q

l)

"=

(8)

In a harmonic crystal the normal modes are independent so that cross terms in the derivative are zero. This leads to:

~

(11)

Substituting eq. 11 in eq. 9 gives:

• eiqk ~ In(Q*/Q) - d(ln Q) Am = ~ 3(ln Q) du/ drn / 3uj - ~ . z D n

= x, y, z) 0 otherwise_]

A

dq

"=

1

~_--~f. + ~

...

dq (10)

which is more familiar in the literature [18]. From Elcombe [16] for a single atom substitution (k) and to first order in Am~m:

• eiqk -2~m*/ q dq

(l 3)

which is valid for 0 < Uiq < 2 sinh -1 (2u~/), (uiq < 4.35464) and converges rapidly for u ~ / < 2. The first term cancels exactly with the 1/Uiq part of

175 eq. 10, thus ensuring that f3 -+ 0 as T-+ ~ : the second term gives the 1/T z dependence at high T:

= ln(Q*/Q)red 3/"

Am

h2

~_1 dq

24m* /~?T 2

riqeiqk

dq

(14)

=

3. Calculation of To make use of any of the above formulae, some model of the crystalline frequencies or the interatomic forces is required. One model which has been used with reasonable success is the Debye-Einstein model [19,20]. This model is very simple, and the reason it gives reasonable answers is clear from eq. 14. ~3is a sum over the first Brillouin zone of u2 e z and is therefore not sensitive to details of the model. However, if the model frequencies are too large on average then ~3 will be too large. If as well as large frequencies the model predicts small eigenvectors then a good value may be obtained, albeit with a large error. In order to reduce the error, as well as obtain a satisfactory calculation of/3, it is necessary to use a more

sophisticated model which reproduces the frequencies and eigenvectors more accurately. An interatomic force model based on recent lattice dynamical theory has been used (see for example Cochran [21 ]). For both materials the parameters of a "shell" model have been obtained by fitting the calculated frequencies to the observed room temperature dispersion curves in the [001], [111] and [110] directions, obtained by inelastic neutron scattering. Also included in the fitting programs were the elastic constants and the high-frequency dielectric constant. The models have then been used to calculate the normal mode frequencies and eigenvectors throughout the Brillouin zone. For lead sulphide, model IV of Elcombe [22] has been used. For zinc sulphide a new model has been fitted to the data of Bergsma [23], and the parameters are given in Table 1. The 13values of both ZnS and PbS were calculated using e q. 12, with 3142 independent wavevectors. The results are shown in Fig. 1, and the corresponding 1000 In a values in Fig. 2. Also included in Fig. 2 are the experimental values from several sources [9,12, 14]. The agreement is quite good, but the calculated value is ( 5 - 2 0 % ) higher than experiment. For comparison, the results of earlier Debye-Einstein model calculations for ZnS and PbS by Sakai [20] and

TABLE 1 Shell model parameters used in the calculations Name

Short-range force parameters

Ionic charge Ionic polarizability (Pb or Zn) Mechanical polarizability (Pb or Zn) Ionic polarizability (S) Mechanical polarizability (S) Lattice parameter * Cochran et al. [26] (Appendix).

Parameter *

Units

PbS

ZnS

A (12) B(12) A(ll) B(ll) A(22) B(22) aS /3T(12) = fiT(21) t3S Z c~l/V dl c~2]V d2 a(A)

e2/v e2/v e2/v e2/v e2/v e2/v

29.27 - 2.61 - 2.55 - 0.068 2.57 - 0.220 -21.5 3.82 1.0 1.625 0.0804 0.057 0.1882 1.008 5.953

49.80 -8.08 -4.64 0.47 2.60 -0.01 1.0 1.0 1.0 2.53 0.061 -0.628 0.1139 1.454 5.406

e 10-24em 3 e 10-24cm 3 e 10-Scm

176 1500 lOCK) 800 700

T (K)

600

4. Errors in ot

400

soo

T(°C)

j

t Zn5

Three sources of error in the calculation have been investigated in an attempt to explain the discrepancy. These are: (1) errors arising from the models, (2) computational errors and errors associated with the approximations used in deriving the formulae, and (3) errors arising from anharmonicity.

J ~5

" J ' J f

i

j-'/"

//

f/J/~

~

NO"ITz

t PbS

4. 1. Error in the models

_j.f11

(g-2)

Fig. 1. The calculated logarithms of reduced partition function ratios of lead sulphide and zinc sulphide. Full line = eq. 7 lattice dynamical model. Dot-dash line = Debye-Einstein model [81, dash line = Debye-Einstein model [201.

Hulston (in Groves et al. [8]) are also shown in Figs. 1 and 2. There appears to be a numerical error in Sakai's PbS calculation, and the figures show values recalculated using his data.

ISOO

9) 8

)ooo soo ,oo

,oo~ ~

66o

~oo

T(K)

6oo

soo

~6o

~

T(°C)

400

,~o

~6o

!

8

o

I

2

3

4

5

6

7

Fig. 2. The calculated and experimental 34S exchange fractionation factor (1000 In a) for the sphalerite-galena pair. Line symbols as for Fig. 1. Experimental points: circles, Kajiwara and Krouse [12]; squares, Grootenboer and Schwartz [9]; triangles, Czamaiaske and Rye [14].

It has not been feasible to calculate this error before because the errors in g(u) or du~/ (eqs. 9 and 10) could not be explicitly identified. Some check on the accuracy ofg(u) can be made by using it to calculate the specific heat of the crystal, which can then be compared with experimental values. This check involves only the frequencies and is very insensitive to the model used. ttowever, the replacement of du~4 by eigenvectors (eq. 11) leading to eq. 14 shows that the eigenvectors are as important a factor in the value of /3 as the frequency squared. To estimate the error in the calculation it is therefore essential that the error in the eigenvectors be estimated as well as the error in the frequencies. Since 13depends on the square of the frequency, its value depends mainly on the optic frequencies. The acoustic contribution is only ~ I 0 - 1 5 % so the error in the acoustic modes has been neglected. For ZnS, using the shell model parameters given in Table 1, the agreement between observed and calculated frequencies is - ( 0 . 3 4 + 0.34)% for the longitudinal optic branch and +(0.38 + 0.39)% for the transverse optic branch. For PbS the corresponding figures are --(1.2 + 0.7)% and (0 + 0.8)%. The accuracy of the eigenvectors is difficult to determine because they have not been measured experimentally. The fact that the models fit the frequencies very well does not imply that the eigenvectors are equally good. It is therefore necessary to compare some other experimental data, which are dependent on the eigenvectors, with their values as calculated by the models, in order to determine the reliability of the eigenvectors as calculated by the models. The most suitable quantity to use is the atomic thermal parameter (B value), given by:

177

2kBT~_ldq~(1 Bk - m k q • e~kq "1 / ~ u~/q

d

+

1 ) exp(u~q)-- 1 Uiq

q

which differs from In(Q*/Q) only by constant terms and 1/u~ (see eq. 12). Experimentally the ratio of the B values is more reliable than the individual values, so it has been used instead. This also overcomes any difficulties associated with anharmonicity since, to a first approximation, the fractional anharmonicity is the same for both atoms. For ZnS there are two sets of experimental room temperature data. Cooper et al. [24] obtained Bzn/Bs = (1.214 + 0.014) using neutron data and Barnea and Post (private communication) obtain 1.238 using X-ray data. No error is available for the latter value but it agrees well with the former. An experimental value of (1.214 -+ 0.024) has been assumed. The model calculates Bzn/B s = 1.298 which is (6.5 -+ 1.8)% higher than the experimental value. This error is converted to an error in the eigenvector by using the eigenvector normalisation which implies that ~"k mkBk = constant, as well as the assumption that B is dominated by the acoustic modes so that B k is proportional to e~ (acoustic). This assumption is only valid for zinc in ZnS. The expression obtained is:

(Ae3 e z ]Zn

1 R

(1 +(R +zXR)mzn/ms) '

(R = Bzn/Bs) Using the figures given above, the average over the zone of the acoustic eigenvector for zinc calculated by the model is (1.8 -+ 0.5)% higher than the experimental value inferred from the B values. Because of the eigenvector normalisation the identical error occurs in the average optic eigenvector for sulphur, and therefore in the dominant terms of the 13calculation. This connection between the B values and/3 has also been checked by numerical calculation. At one time we had a second force model for ZnS which gave an equally good fit to the frequencies and elastic constants. Also the specific heat calculated for each of the two models agreed to within 3~% at all temperatures.

However, the calculated ratio of the B values was 0.921 ( 2 0 - 3 0 % lower than both experiment and the present calculation). Using the two calculated values of the B ratio and the above relations predicts that fi for the second model should be 10.1% lower than for the present model. This is in excellent agreement with the calculated fall of 10.0% (at 100°C, 1000/3 (ZnS) is 8.149 for the second model and 9.053 for the present model). The agreement also confirms that the relations used are correct and that the assumption B k proportional to e 2 (acoustic) is valid for Zn in ZnS. Unfortunately the B values of PbS have not been measured, but an upper limit in the eigenvector uncertainty may be estimated. Because of the large mass difference e s2 (optic) is high (0.86 at q -- 0) and because of the greater symmetry of PbS e S2 (optic) -- l at the [111] zone boundary point [q = 27r/a (0.5, 0.5, 0.5)]. The eigenvectors are unlikely to go much outside this range, giving a maximum uncertainty of

-+71%. Combining the frequency and eigenvector errors for ZnS yields an error o f - ( 1 . 8 -+ 0.8)%; that is, the calculated value of/3 for ZnS is estimated to be 1.8% too high, with an uncertainty of 0.8%. For PbS the frequency errors are negligible compared with the eigenvector error so that the uncertainty in/3 for PbS is Because/3 for ZnS is approximately 3t3 for PbS the error in In a comes mainly from ZnS. The final value is - ( 2 . 7 -+ 5)%; that is, In a(ZnS-PbS) as calculated by the models is estimated to be 2.7% high with an uncertainty of 5%. As was stated earlier the calculation of 1000 In a is 5 - 2 0 % larger than the observations. Of this discrepancy, 2.770 comes from the systematic error in the model for ZnS. The agreement between theory and experiment is therefore satisfactory for all but the most recent data [14].

4.2. Errors from approximations and computation In the previous section we have shown that known errors in the models of ZnS and PbS give rise to an uncertainty of -+5% in/3. Before a final error in the calculation can be given errors arising from the numerical calculations and the approximations used in obtaining the formulae must also be estimated. The numerical accuracies of the different calculation methods are compared in Table 2 which gives

178 TABLE 2 Numerical errors in ~ from the different methods of calculation Temperature 373°K (100°C)

1273°K (1000°C)

9.053

ZnS 1000 ln(Q*/Q)red Method Eq. 7 Eq. 12 32S Eq. 12 34S Eq. 6 ~v = 0.01 THz

0.791 set A 0 +0.003 +0.003 -0.006 + 0.037

set B +0.003 +0.003 +0.003 -0.152

set A 0 -0.004 +0.007 -0.008

PbS 1000 ln(Q*/Q)red Method Eq. 7 Eq. 12 3eS Eq. 12 34S Eq. 6 d~v= 0.01 THz

0.266 set A 0 -0.003 0.003 +0.059 + 0.053

set B -0,001 -0,003 -0.003 +0.677

set A 0 -0.003 -0.002 +0.058

set B +0.004 -0.001 +0.010 -0.155 3.043 set B +0.001 -0.002 0 +0.685

Set A = 3142 independent wavevectors, set B = 48 independent wavevectors.

1000 ln(Q*/Q)red at 100°C and 1000°C and the differences between these values and those calculated using various formulae. For this purpose t3 was calculated using eq. 7, with 3142 independent wavevectors. This formula is the most accurate for any particular set of wavevectors provided care is taken to eliminate rounding errors. The errors introduced by using a smaller set of wavevectors (48) is very small ( - 0 . 0 0 1 to +0.004) because/3 is an average over the whole Brillouin zone. In other words, eq. 7 gives a converged result with a small set of independent wavevectors. The errors from the approximations used to obtain eq. 12, namely the differential form, the independence of the normal modes and the neglect of terms containing (Am~m) 2 and higher powers, are slightly larger ( - 0 . 0 0 4 to +0.010). These errors are insignificant compared with the errors from the models, thus justifying the use of the approximations, at least for sulphur substitution in which Am~m* = I/17. The 34S results have been included to demonstrate that normal and isotopic values can be interchanged throughout eq. 12. On the other hand eq. 6 introduces large errors which are largely independent of temperature and which can be positive or negative. The only extra factor involved in eq. 6 is the frequency interval into which g(u) is divided. An interval of 0.01 THz

(1/3 cm -1) was used in the calculations and so the large errors occurring indicate that this interval is too big. This is particularly so for PbS in which the frequency shifts on isotopic substitution are quite small (~0.10 THz). The errors are also much worse for the small wavevector set in contrast to eqs. 7 and 12. For PbS, to obtain a converged result using eq. 6 would require a frequency interval of less than 0.002 THz and a large wavevector set, and this would involve more calculation than the exact eq. 7. Eq. 6 is therefore not suitable for an accurate calculation of ~3, whereas eq. 12, with a small wavevector set gives a numerically accurate answer for less than 1% of the calculation effort. Numerical errors do not therefore contribute to the uncertainty in the calculated value of/3.

4. 3. Errors from anharmonicity Underlying all that has been stated so far is the assumption that the harmonic approximation is valid for the temperature region of the calculations. If this is not so then an estimate of error should be made to allow for this. This is much more difficult than the numerical errors, but two indicators of tile anharmonic content are available and have therefore been used. The first indicator is the difference between the experimental and calculated B values (absolute values as

179 opposed to the ratio used earlier). The difference is ~8% at room temperature despite the fact that the models fit the frequencies at room temperature and therefore include some anharmonicity. Using both eigenvector nornralisation and normal mode population probabilities an error in/3 of this magnitude is possible at ~500°C, and the calculation would be too low. On the other hand as the temperature increases, the normal mode frequencies decrease, mainly because of thermal expansion which gives rise to a reduction in the force constants. In KBr for example the optic frequencies decrease by ~ ( 4 - 8 ) % between 90°K and 400°K (OD/2 to 2(9D) [25]. Assuming similar shifts (with respect to ®D) for ZnS and PbS would give calculated/3 values at ~500°C which are respectively ~12% high and ~15% high, leading to a calculated value of In c~ ~10% high at 500°C. The eigenvectors are not expected to change nmch with temperature. It is therefore possible for the calculations to be out by ~10% either way depending on the relative magnitude of these two anharmonic contributions to the error. Coupling this with the 5% uncertainty from the previous section it is unlikely that the calculation is out by more than 15-20%. The agreement between theory and experiment is therefore very satisfactory.

5. Conclusion This paper has presented an alternative formula for calculating reduced partition function ratios of crystals. The new formula, in particular the high-temperature form, explains how the reduced partition function ratio depends on the normal mode frequencies and eigenvectors of the crystal. It has therefore been possible to estimate the accuracy of the calculation, by determining the accuracy with which the models predict other experimental quantities which depend in a similar way on the frequencies and eigenvectors. It has also been shown that to calculate a numerically accurate value offl, it is n o t necessary to calculate an accurate frequency distribution for the crystal as has been suggested [15]. An exact calculation for a few wavevectors is much more accurate than an approximate calculation at a large number of wavevectors. It is also a lot quicker. The 34S fractionation factor for Z n S - P b S calcu-

lated using the best lattice dynamical models available, is only reliable to -+5% at room temperature, becoming worse as T increases; but it agrees satisfactoriJy with the experimental data.

References 1 H.C. Urey and L.J. Greiff, Isotopic exchange equilibria, J. Am. Chem. Soc. 57(1935) 321 327. 2 lt.C. Urey, The thermodynamic properties of isotopic substances, J. Chem. Soc. (1947) 562-581. 3 A.P. Tudge and H.G. Thode, Thermodynamic properties of isotopic compounds of sulphur, Can. J. Rcs. 28B (1950) 567-578. 4 H. Sakai, Fractionation of sulfur isotopes in nature, Geochim. Cosmochim. Acta 12 (1957) 150-169. 5 S. Gavelin, A. Parwell and R. Ryhage, Sulfur isotope fractionation in sulfide mineralisation, Econ. Geol. 55 (1960) 510-530. 6 T. Tatsumi, Sulfur isotopic fractionation between coexisting sulfide minerals from some Japanese ore deposits, Econ. Geol. 60 (1965) 1645-1659. 7 R.L. Stanton and T.A. Rafter, Sulphur isotope ratios in co-existing galena and sphalerite from Broken Hill, New South Wales, Econ. Geol. 62 (1967) 1088-1091. 8 D.I. Groves, M. Solomon and T.A. Rafter, Sulfur isotope fractionation and fluid inclusion studies at the Rex ttill Mine, Tasmania, Econ. Geol. 65 (1970) 459-469. 9 J. Grootenboer and H.P. Schwarcz, Experimentally determined sulfur isotope fractionations between sulfide minerals, Earth Planet. Sci. Lett. 7 (1969) 162 166. 10 Y. Kajiwara, H.R. Krouse and A. Sasaki, Experimental study of sulfur isotope fractionation between coexisting sulfide minerals, Earth Planet. Sci. Lett. 7 (1969) 271 277. 11 R.O. Rye and G.K. Czamanske, Experimental determination of sphalerite galena sulfur isotope fractionation and application to ores at Providencia, Mexico, Geol. Soc. Am. Abstr. 7(1969) 195 196. 12 Y. Kajiwara and H.R. Krouse, Sulfur isotope partitioning in metallic sulfide systems Can. J. Earth Sci. 8 (1971) 1397-1408. 13 W. Salomons, Isotope fractionation between galena-pyrite and pyrite-sulphur, Earth Planet. Sci. Lett. 11 (1971) 236 238. 14 G.K. Czamanske and R.O. Rye, Experimentally determined sulfur isotope fractionations between sphalerite and galena in the temperature range 600° to 275°C, Econ. Geol. 69 (1974) 17-25. 15 Y. Bottinga, Carbon isotope fractionation between graphite, diamond and carbon dioxide, Earth Planet. Sci. Lett. 5 (1969) 301-307. 16 M.M. Elcombe, Determination of atomic eigenvector magnitudes by isotopic substitution, J. Phys. C: Solid State Phys. 7 (1974) L202-205.

180 17 J.E. Mayer and M.G. Mayer, Statistical Mechanics (John Wiley and Sons, New York, 1940) 239-240. 18 J. Bigeleisen and M.G. Mayer, Calculation of equilibrium constants for isotopic exchange reactions, J. Chem. Phys. 15 (1947) 261-267. 19 J.M. McCrea, On the isotopic chemistry of carbonates and the paleotemperature scale, J. Chem. Phys. 18 (1950) 849 857. 20 H. Sakai, Isotopic properties of sulfur compounds in hydrothermal processes, Geochem. J. 2 (1968) 29 49. 21 W. Cochran, Crit. Rev. Solid State Sci. 2 (1971) 1 44. 22 M.M. Elcombe, The crystal dynamics of lead sulphide, Proc. R. Soc. Lond., Set. A, 300 (1967) 210-217.

23 J. Bergsma, Lattice dynamics of magnesium stannide and zinc blende, an experimental study by means of inelastic neutron scattering, Reactor Centrum Nederland Rep. RCN-121 (1970). Briefly reported as: Lattice dynamics of zinc blende, Phys. Lett. 32A(1970) 324 325. 24 M.J. Cooper, K.D. Rouse and H. l:uess, A neutron-diffraction study of ZnS and ZnTe, Acta Crystallogr. A29 (1973) 49-56. 25 A.D.B. Woods, B.N. Brockhouse and R.A. Cowley, Lattice dynamics of alkali halide crystals II. Experimental studies o f K B r a n d NaI, Phys. Rev. 131 (1963) 1025 1029. 26 W. Cochran, R.A. Cowley, G. Dolling and M.M. Elcombe, The crystal dynamics of lead telluride, Proc. R. Soc. Lond., Set. A, 932 (1966) 433-451.