Oxygen and silicon β-factors of zircon estimated from first principles

Cover Page

Abstract


Zircon β-factors have been calibrated against temperature for isotopic substitutions of 18O/16O and 30Si/28Si. Calculations were performed using the density functional theory (DFT) with the “frozen phonon” approach. The deduced geometric parameters of the zircon unit cell, and the phonon frequencies calculated, agree well with the experimental data. The results are expressed by the cubic polynomials on x = 106/T(K)2: 1000lnβzrn(18O/16O) = 9.83055x – 0.19499x2 + 0.00388x3;  1000lnβzrn(30Si/28Si) = 7.89907x – 0.17978x2 + 0.00377x3.

The expressions deduced can be utilized to construct geothermometers if combined with β-factors of coexisting phases. New calibrations of quartz-zircon are given. The new values of 1000lnβzrn and the estimated isotope fractionation factors between quartz and zircon (1000lnβqtz–1000lnβzrn) deviate considerably from previously used experimental, empirical, and semi-empirical calibration of the isotopic equilibrium.


ВВЕДЕНИЕ

Циркон (ZrSiO4) встречается в породах различного химического состава и генезиса и представляет важные свидетельства условий, соответствующих процессам, проходящим на разных стадиях образования и преобразования магматических, осадочных и метаморфических комплексов. В ряде случаев кристаллохимические особенности циркона отражают первичные условия кристаллизации даже при наложении высокотемпературных преобразований, что дает возможность использовать циркон в качестве уникального трассера древнейших процессов в докембрийских комплексах (например, Valley et al., 2014). Возможность вхождения урана в кристаллическую решетку циркона при замещении циркония в октаэдрических позициях (Zr4+замещается на U4+) и низкая подвижность урана и продуктов его распада определяют наиболее распространенное применение циркона в качестве U-Pb геохронометра. Вместе с тем важнейшие перспективы цирконометрии обуславливаются возможностью сопоставления результатов датирования геологических процессов или процессов кристаллизации, с одной стороны, и определения условий процессов на основе анализа стабильных изотопов (кислорода, кремния) – с другой. Таким образом, возможна непосредственная петрологическая интерпретации результатов датирования (изотопная термохронометрия).

Необходимым условием использования изотопных отношений кислорода циркона в качестве геохимических трассеров (в геотермометрии; при определении условий взаимодействия пород с флюидами; оценках степени достижения равновесия реакций изотопного обмена) является знание величин (факторов) фракционирования используемых стабильных изотопов между цирконом и другими фазами в зависимости от температуры. Для определения факторов фракционирования 18O/16O циркона применялись различные методы: полуэмпирический метод динамики решетки (Kieffer, 1982), полуэмпирические методы «инкрементов» (Hoffbauer et al., 1994) и «модифицированных инкрементов» (Zheng, 1993), метод позиционного (site) потенциала (Smyth, 1989), эмпирический метод определения распределения изотопов в природных ассоциациях (Valley et al., 2003) и, наконец, экспериментальные калибровки (Krylov et al., 2002; Trail et al., 2009). Явные или неявные упрощения и предположения, определяющие применимость перечисленных методов, приводят к существенно разным результатам и, как правило, затрудняют оценки точности и достоверности полученных калибровок (Trail et al., 2009). «Из первых принципов» (без использования каких-либо эмпирических предположений) β-факторы циркона были рассчитаны на основе теории функционала плотности в сочетании с методом динамики решетки (Qin et al., 2016). Также на основе методов DFT были рассчитаны (Chaplot et al., 2006) термодинамические параметры циркона, плотность фононных состояний и ее проекции на колебания подрешеток кремния и кислорода, что также дает возможность оценивать β-факторы. Перспективы изотопных DFT калибровок обуславливаются как применимостью к неограниченному набору минералов, так и возможностью расчета факторов фракционирования любых пар изотопов.

Цель настоящей работы – расчет кислородных (18O/16O) и кремниевых (30Si/28Si) β-факторов циркона путем сочетания квантово-химических расчетов на основе теории функционала плотности (DFT) и метода супер-ячеек (например, Evarestov, Smirnov, 1997).

МЕТОДЫ РАСЧЕТА

Общая методология

В общепринятых обозначениях при равновесных условиях коэффициент (фактор) изотопного фракционирования между фазами A и B, βA/B равен:

1000 ln αA/B = 1000 lnβA – 1000 lnβB ≈ δA – δB, (1.1)

где β – величины так называемых приведенных отношений статистических сумм (или β-факторы), которые численно равны константе равновесия данной фазы (соединения) с одноатомным идеальным газом рассматривае­мого элемента (Bigeleisen, Mayer, 1947; Kieffer, 1982). Необходимо отметить, что значения β-факторов относятся к индивидуальным фазам, тогда как значения факторов изотопного фракционирования отражают распределение изотопов между разными фазами.

Для идеального кристалла в гармоническом приближении (т.е. когда частоты колебаний кристаллической решетки и объем ячейки постоянны при всех рассматриваемых условиях и соответствуют равновесной структуре кристалла) β-факторы можно вычислить (например, Meheut et al., 2007) с помощью соотношения:

lnβ=1Nqq1Ni=  13Natlnνq,i*νq,isinh(hνq,i/2kT)sinh(hνq,i*/2kT) , (1.2)

где vq,i – частота нормального колебания (фонона) с волновым вектором q первой зоны Бриллюэна и индексом фононной ветви i от 1 до 3Nat (Nat – количество атомов в примитивной ячейке). T – температура (градусы Кельвина), h – постоянная Планка, k – постоянная Больцмана, N – число атомов элемента, подвергающегося изотопному замещению в ячейке, Nq – количество векторов q, учитываемых при суммировании в первой зоне Бриллюэна. Надстрочный индекс * относится к более тяжелому изотопу. Выражение в квадратных скобках представляет собой вклад фононов с волновым вектором q с учетом полного изотопного замещения в примитивной ячейке, lnβq. Таким образом, lnβ представляется как средний вклад фононов lnβq по волновым векторам первой зоны Бриллюэна. В приведенном выражении β-фактор нормирован на единицу при температуре, стремящейся к бесконечности (с учетом правила Редлиха–Теллера).

Необходимо отметить, что Nq должно быть достаточно большим для достижения необходимой точности вычисления β-факторов (суммирование по бесконечному количеству волновых векторов обеспечит «точный» результат, но недостижимо). В настоящей работе ограничения, обусловленные конечностью набора волновых векторов при суммировании, Nq преодолеваются путем определения зависимости lnβ от Nq (Krylov, Evarestov, 2018).

Вычисления

Частоты фононов вычислены с использованием метода «замороженных фононов» (реализованного в программе CRYSTAL14, Dovesi et al., 2016) в гармоническом приближении с применением набора гауссовых полноэлектронных базисов и гибридного функционала B3LYP-D*, в котором обменный функционал дополнен локальным корреляционным функционалом с дисперсионной поправкой (Grimme, 2006). Для Si и O использовались наборы полноэлектронных базисов, опубликованных на сайте CRYSTAL (все наборы базисов и их обозначения приведены на сайте http://www.crystal.unito.it) с совокупностью базисных функций 88-1111G(d) (Porter et al., 1999) и m-6-311G(d) (Heyd, 2005), соответственно. Для внешних электронных оболочек Zr адаптирован набор базисных функций SC_HAYWSC-3111(32111dfG) (Sophia et al., 2013). Однако, вместо использованного ранее (ibid.) псевдопотенциала (в котором внутренние электронные орбитали представлены усредненным потенциалом, ECP), нами адаптирована остовная часть полноэлектронного базиса Dovesi (http://www.crystal.unito.it). Для улучшения точности вычислений (которые в программе CRYSTAL контролируются параметром TOLINTEG) выбраны значения 8 8 8 9 30, критерий сходимости энергии самосогласованного поля (SCF) установлен на уровне 10–11 Хартри. Узлы суммирования обратной решетки определяются фактором сжатия IS = 3, что соответствует шести (пространственная группа симметрии циркона I41/amd) независимым k-векторам неприводимой части зоны Бриллюэна. При расчетах матричных элементов обменно-корреляционного потенциала применена прецизионная сетка из 99 радиальных и 1454 угловых точек (параметр XXLGRID).

Процесс оптимизации геометрии решетки включал релаксацию координат ядер атомов и параметров решетки (FULLOPTG) по квази-ньютоновскому алгоритму. Сходимость в процессе оптимизации геометрии оценивалась по среднеквадратичному отклонению (RMS) и абсолютному значению наибольшей компоненты градиентов и смещений ядер для всех атомов. Для максимальных и среднеквадратичных (RMS) градиентов выбраны уточненные значения (0.0001, 0.00003 a.u); для максимальных значений смещений атомов приняты стандартные величины (0.00180, 0.00120 a.u.). Достоверность проведенных DFT вычислений подтверждается сравнением вычисленных в результате оптимизации геометрических параметров ячеек с экспериментальными. Отклонения от параметров решетки, полученных экспериментально (Mursic et al., 1992) для кристаллографических параметров примитивной ячейки, составляют 0.3% (6.592 Å), –1.3% (6.081 Å) и –0.8% (264.247 Å3) для a, c и V, соответственно.

Результаты расчетов

Для определения температурной зависимости β-факторов (теоретических калибровок β) были: а) определены частоты фононов для различных изотопологов циркона и изменения частот при изотопных замещениях; б) вычислены значения β (в общепринятых терминах 1000 lnβ) для различных дискретных значений температур (с учетом пределов устойчивости циркона); в) проведена интерполяция 1000 lnβzrn с помощью подходящего многочлена, 1000 lnβzrn(T).

Вычисления частот фононов

Как отмечено выше, вычисленные величины частот колебаний зависят от выбранных q-векторов, по которым производится суммирование в пределах зоны Бриллюэна. В рамках расчетов, реализованных в программе CRYSTAL, точность средних значений в уравнении (1.2), полученных при суммировании по q-векторам, определяется построением расширенных ячеек (супер-ячеек, supercells) посредством линейных комбинаций (целочисленных линейных преобразований) единичных векторов примитивных ячеек (Dovesi et al., 2016). Расчеты по супер-ячейкам большего объема (т.е. с большими определителями матриц преобразований, detL = Nq) в общем уменьшают ошибки определения средних значений, но ограничиваются вычислительными ресурсами. Для оценки влияния Nq на результаты расчетов в настоящей работе частоты колебаний циркона вычислены по расширенным ячейкам варьирующего объема, соответствующим Nq = 1 (исходная ячейка), Nq = 2, 4, 16 и 27. Преобразования определялись матрицами с элементами строк {1 0 0/0 1 0/0 0 1} (исходная ячейка), {0 1 1/1 0 1/1 1 0} (Nq = 2), {2 1 1/1 2 1/1 1 2}  (Nq = 4), {2 0 0/0 2 0/0 0 2} (Nq = 8), {0 2 2/2 0 2/2 2 0}  (Nq = 16), {3 0 0/0 3 0/0 0 3} (Nq = 27). Такие расширенные ячейки соответствуют суммированию дисперсий фононов по 1, 2, 4, 8, 16 и 27 k-точкам в пределах первой зоны Бриллюэна по схеме Монкхерста–Пака (Monkhorst, Pack, 1976; Evarestov, Smirnov, 1997). Независимо от объема супер-ячейки, наборы частот колебаний в центре зоны Бриллюэна (q = 0) остаются неизменными (центр супер-ячейки не сдвигается).

Нормальные колебательные моды циркона в центре зоны Бриллюэна (q = 0) по симметрии разделяются по неприводимым представлениям на

Г: 2A1g + A2g + A1u + 4A2u + 4B1g + B2g + B1u + 2B2u + 5Eg + 5Eu.      (2.1)

Моды A1g, B1g, B2g и Eg активны в спектре комбинационного рассеяния (Рамана), а моды A2u и Eu активны в инфракрасном спектре (моды 1Eu и 1A2u соответствуют трансляциям), остальные моды неактивны. Для оценки качества вычислений полученные значения частот колебаний сопоставлены с имеющимися экспериментальными данными при Nq= 1 (табл. 1). Средние абсолютные отклонения вычисленных частот относительно экспериментальных данных (|Δω| = Σ |ωi – ωiexp| /n; n = 22) составляют около 9.12 см–1, среднеквадратичная ошибка вычисленных частот составляет 12.59 см–1. Среднее отклонения частот с учетом знака (параметр Δω = Σ(ωi – ωiexp)/n), принятое в качестве обобщенной меры неопределенности вычисленных частот (Montanari et al., 2006), равно 0.16. Таким образом, при выбранных параметрах экспериментальные частоты воспроизводятся с хорошей точностью. Расчеты плотностей состояний фононов (PDOS) по результатам расчетов колебаний супер-ячейки (Nq = 27) также хорошо соответствуют экспериментальной PDOS (рис. 1).

Определение β-факторов

 

Рис. 1. Плотности состояний фононов PDOS. По горизонтали – энергия фононов, по вертикали – плотности состояний (произвольная шкала). Залитые кружки – данные по неупругому рассеянию нейтронов (Chaplot et al., 2006), сплошная и пунктирная линии – результаты настоящей работы (для изотопологов циркона 18O и 16O, соответственно при Nq = 27).

 

Таблица 1. Частоты колебаний для изотопомеров 16O (ω16), 18O (ω18) и 30Si (ω30), вычисленные по примитивной ячейке циркона (Nq = 1)

Симметрия моды

ω16 (см-1)

ω18 (см–1)

ω1816

ω30 (см–1)

ω3028

Данные экспериментов(1)

A1g

428.8074

404.2192

0.9427

428.8074

1.0000

431.51

A1g

992.0560

935.1706

0.9427

992.0560

1.0000

975.12(3)

A1u

444.6490

419.1525

0.9427

444.6490

1.0000

411.34

A2g

258.4119

243.5943

0.9427

258.4119

1.0000

263.74

A2u

324.5738

314.173

0.9682

324.3707

0.9994

339.56(2)

A2u

609.6329

591.0584

0.9694

600.6355

0.9852

606.53

A2u

987.5852

955.5091

0.9675

973.9879

0.9862

977.54(2)

B1g

222.9582

222.9209

0.9998

222.8004

0.9993

211.32

B1g

391.3808

375.8959

0.9613

386.8127

0.9883

388.76

B1g

647.4292

620.5426

0.9561

641.0871

0.9902

661.37

B1g

1008.2596

973.4347

0.9670

995.6188

0.9875

1009.00(3)

B1u

119.6253

112.7659

0.9427

119.6253

1.0000

143.57

B2g

257.3132

242.5587

0.9427

257.3132

1.0000

274.23

B2u

576.1301

543.0943

0.9427

576.1301

1.0000

 

B2u

964.2912

908.9978

0.9427

964.2912

1.0000

 

Eg

199.1376

189.7384

0.9540

198.0704

0.9946

204.86

Eg

226.9030

222.8685

0.9811

226.0819

0.9964

224.22

Eg

376.6658

364.8171

0.9687

374.5551

0.9944

361.34

Eg

559.8031

537.1074

0.9591

554.8781

0.9912

553.30

Eg

924.6715

890.6277

0.9633

914.1405

0.9886

926.73

Eu

287.0058

278.3264

0.9701

286.6207

0.9987

288.75

Eu

384.8811

362.9566

0.9429

384.7502

0.9997

388.76

Eu

431.3852

418.7353

0.9703

424.7942

0.9847

435.54

Eu

869.483

838.4171

0.9644

858.8815

0.9878

871.88(2)

Примечание. ω1816, ω3028 – сдвиги частот при изотопных замещениях, выраженные как отношения частот соответствующих колебательных мод.

(1)Данные по цирконам с природными отношениями изотопов: результаты экспериментов по неупругому рассеянию нейтронов (Chaplot et al., 2006), (2) (Nipko, Loong, 1997); (3)данные рамановской спектроскопии (Dawson et al., 1971).

 

β-факторы изотопного фракционирования кислорода и кремния циркона (в виде логарифмических функций 1000 lnβzrn) вычислены согласно уравнению (1.2) для 16O/18O и 28Si/30Si при температурах от 0 дo 2000°C с интервалом 20°C (в табл. 2 выборочно с шагом 100°C показаны 1000 lnβzrn для Nq = 1, 4, 8, 16 и 27). Результаты представлены посредством кубических полиномов Ax + Bx2 + Cx3, x = 106/T(K)2 (Clayton, Kieffer, 1991; Chacko et al., 2001) с помощью аппроксимации β-факторов методом наименьших квадратов в пределах рассматриваемого интервала температур. Значения 1000 lnβzrn исходной ячейки (Nq = 1) существенно отличается от 1000 lnβzrn, полученных по супер-ячейкам (например, при 0°С 1000 lnβ исходной ячейки на 4% меньше, чем по супер-ячейке Nq = 27, табл. 2). При этом наблюдается линейная зависимость между полиномиальными коэффициентами разложения 1000 lnβzrn и обратной величиной Nq–1. В пределе, при стремлении количества волновых векторов к бесконечности (или, что то же самое, при стремлении  Nq–1 к 0), полиномиальные коэффициенты стремятся к значениям A = 9.83055, B = –0.19499 и C = 0.00388 (18O/16O). Таким образом, 1000 lnβzrn при стремлении количества волновых векторов к бесконечности представляется в виде:

1000 lnβzrn(18O/16O) = 9.83055x – 0.19499x2 +  + 0.00388x3. (2.2)

Аналогично, в пределе Nq = ∞ полиномиальные коэффициенты 1000 lnβzrn(30Si/28Si) сходятся к значениям:

1000 lnβzrn(30Si/28Si) = 7.89907x – 0.17978x2 +  + 0.00377x3. (2.3)

 

Таблица 2. Вычисленные 1000 lnβzrn при разных температурах (°C) и объемах супер-ячеек и аппроксимирующие коэффициенты кубических полиномов

Nq

1

2

4

8

16

27

Chaplot et al., 2006

Qin et al., 2016

16O/18O

0

102.00

103.93

105.18

105.57

105.74

105.88

103.30

 

100

59.58

60.66

61.42

61.68

61.79

61.86

60.47

62.13

200

38.80

39.48

39.99

40.16

40.24

40.28

39.38

40.25

300

27.13

27.60

27.96

28.08

28.14

28.17

27.56

28.15

400

19.97

20.32

20.58

20.68

20.72

20.74

20.31

20.75

500

15.30

15.56

15.76

15.83

15.87

15.88

15.56

15.91

600

12.08

12.28

12.44

12.50

12.53

12.54

12.30

12.57

700

9.77

9.93

10.07

10.11

10.13

10.15

9.95

10.18

800

8.06

8.20

8.31

8.35

8.36

8.37

8.22

8.40

900

6.77

6.88

6.97

7.00

7.02

7.03

6.90

7.05

1000

5.76

5.85

5.93

5.96

5.97

5.98

5.87

6.00

1100

4.96

5.04

5.11

5.13

5.14

5.15

5.06

5.17

1200

4.31

4.38

4.44

4.46

4.47

4.48

4.40

4.50

1300

3.79

3.85

3.90

3.92

3.93

3.93

3.86

3.95

1400

3.35

3.41

3.45

3.47

3.48

3.48

3.42

3.50

1500

2.99

3.03

3.08

3.09

3.10

3.10

3.05

3.12

1600

2.68

2.72

2.76

2.77

2.78

2.78

2.73

2.79

1700

2.41

2.45

2.49

2.50

2.50

2.51

2.46

2.52

1800

2.19

2.22

2.25

2.26

2.27

2.27

2.23

2.28

1900

1.99

2.02

2.05

2.06

2.07

2.07

2.03

2.08

2000

1.82

1.85

1.88

1.89

1.89

1.89

1.86

1.90

A

9.4451

9.6013

9.7305

9.7780

9.7977

9.8098

9.6178

9.87055

B

-0.1868

-0.1873

-0.1912

-0.1932

-0.1940

-0.1944

-0.1957

-0.22965

C

0.0037

0.0037

0.0038

0.0038

0.0039

0.0039

0.0040

0.00832

28Si/30Si

0

80.82

81.75

82.06

82.45

82.57

82.64

84.18

 

100

47.63

48.24

48.45

48.72

48.80

48.84

50.76

47.88

200

31.15

31.59

31.75

31.94

31.99

32.02

33.29

31.14

300

21.83

22.15

22.28

22.41

22.45

22.48

23.39

21.84

400

16.10

16.34

16.44

16.54

16.57

16.59

17.28

16.13

500

12.34

12.53

12.60

12.69

12.71

12.72

13.27

12.38

600

9.75

9.90

9.96

10.03

10.04

10.06

10.50

9.79

700

7.89

8.02

8.06

8.12

8.13

8.14

8.50

7.93

800

6.51

6.62

6.66

6.70

6.72

6.72

7.02

6.56

900

5.47

5.55

5.59

5.63

5.64

5.64

5.90

5.51

1000

4.65

4.73

4.76

4.79

4.80

4.80

5.02

4.69

1100

4.01

4.07

4.10

4.12

4.13

4.14

4.33

4.04

1200

3.49

3.54

3.57

3.59

3.60

3.60

3.77

3.51

1300

3.06

3.11

3.13

3.15

3.16

3.16

3.31

3.09

1400

2.71

2.75

2.77

2.79

2.79

2.80

2.93

2.73

1500

2.41

2.45

2.47

2.49

2.49

2.49

2.61

2.44

1600

2.16

2.20

2.21

2.23

2.23

2.24

2.34

2.18

1700

1.95

1.98

2.00

2.01

2.01

2.02

2.11

1.97

1800

1.77

1.80

1.81

1.82

1.83

1.83

1.91

1.79

1900

1.61

1.64

1.65

1.66

1.66

1.66

1.74

1.63

2000

1.47

1.50

1.51

1.52

1.52

1.52

1.59

1.49

A

7.6405

7.7682

7.8168

7.8705

7.8861

7.8954

8.2191

7.71985

B

-0.1639

-0.1718

-0.1751

-0.1780

-0.1789

-0.1794

-0.1802

-0.20269

C

0.0033

0.0035

0.0036

0.0037

0.0037

0.0038

0.0027

0.00782

Примечание. Полиномиальные коэффициенты аппроксимированы с шагом 20°C с коэффициентом детерминации R2 ≈ 0.99; в таблице приведены выборочные 1000 lnβzrn с шагом 100°C. Расчеты по данным (Chaplot et al., 2006) проведены с применением программы В.Б. Полякова.

 

Следует подчеркнуть, что представленные уравнения не учитывают влияние ангармоничности, для учета которой необходимо использовать уравнение, отличающееся от уравнения Бигелейзена–Майер (1.2). Источниками ошибок 1000 lnβ могут быть отклонения сдвигов частот при изотопных замещениях вследствие ангармоничности колебаний, неполного учета электронных корреляций, использования конечных наборов базисных функций для атомов (например, Blanchard et al., 2009; Schauble, 2011; Widanagamage et al., 2014). Порядок возможных отклонений β-факторов от вычисленных значений можно оценить путем распространения ошибок Δω (Blanchard et al., 2009). В частности, зависимости полиномиальных коэффициентов от вариаций частот Δω циркона при выбранном гамильтониане передаются выражениями:

A = 9.83055 + 0.04869×Δω (R2 = 0.97),

B = –0.19499 – 0.00215×Δω (R2 = 0.93),

C = 0.00388 + 0.00006×Δω (R2 = 0.92).

Таким образом,

(∂1000 lnβzrn)/(∂Δω) = –0.04869×Ax + + 0.00215Bx2 + 0.00006Cx3,

и при погрешности вычисленных частот колебаний Δω = 0.16 (табл. 1) неопределенность 1000 lnβzrn составляет по порядку величины 0.1×106/T2. Например, при температуре 0°C возможные отклонения 1000 lnβzrn оцениваются в ≈1.3 ‰, при 500°С – ≈0.16‰, при 1000°С – ≈0.06 ‰.

 

Рис. 2. Сравнение калибровок изотопного фракционирования 18O/16O между кварцем и цирконом.

(а) Зависимость факторов изотопного фракционирования Δ = 1000 lnβqtz–1000 lnβzrn от температуры, x = 106/T2(K). В качестве 1000 lnβqtz приняты значения: 1000 lnβqtz = 12.55277x – 0.41976x2 + 0.01979x3 (Qin et al., 2016).

(б) Отклонения температур, полученных с применением различных калибровок изотопного равновесия кварц–циркон, от результатов настоящей работы.

DFT – βzrn по соотношению (2.2), Δ = 2.72222x – 0.122477x2 + 0.01591x3; Zhe – методом «модифицированных инкрементов» (Zheng, 1993), Δ = 0.72x + 4.26x1/2 – 1.79; VBP – «природная» калибровка (Valley et al., 2003), Δ = 2.64x, Kie – расчеты на основе совокупности полуэмпирических правил для определения сдвигов частот при изотопных замещениях (Kieffer, 1982), Δ = 3.1x; QWWH – ab-initio расчеты методом «линейного отклика» с использованием псевдопотенциалов в виде суперпозиции плоских волн (Qin et al., 2016), Δ = 2.68222x – 0.19011x2 + 0.01147x3. Экспериментальная калибровка Δ = 2.33x (Trail et al., 2009) в ограниченном интервале 700–1000°C в масштабе рис. 2a совпадает с DFT и показана на рис. 2б (TBWS).

Полученные значения β-факторов циркона близки значениям, полученным с применением методов динамики решетки и расчетам по данным неупругого рассеяния нейтронов (табл. 2, см. также ниже).

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Сравнение факторов изотопного фракционирования 18O/16O между цирконом и кварцем, определенных различными методами

Для использования β-факторов циркона в геотермометрии, они должны быть скомбинированы с β-факторами других минералов. Можно предположить, что более точные факторы фракционирования получаются при близких параметрах расчетов, включая одинаковые методы DFT-вычислений и согласованные наборы базисов для одинаковых атомов. Для сравнения с калибровками изотопного фракционирования кислорода циркона, полученного эмпирическими («природные» калибровки, Valley et al., 2003), полуэмпирическими (калибровки методом инкрементов, Zheng, 1993; Hoffbauer et al., 1994; калибровками с использованием закономерностей изменений частот при изотопных замещениях, Kieffer, 1982) и экспериментальными калибровками (Krylov et al., 2002; Trail et al., 2009) в качестве единой референтной фазы, чаще всего, используется кварц. Относительно кварца определяются факторы изотопного фракционирования кварц–циркон, которые согласно соотношению (1.1) равны 1000 lnβqtz–1000 lnβzrn. На рис. 2а показаны факторы фракционирования между кварцем и цирконом, определенные разными методами в зависимости от температуры. На рис. 2б демонстрируются возможные отклонения температур, определенных с помощью различных калибровок относительно рассчитанных в настоящей работе. Значительные (более ста градусов) расхождения могут обуславливаться, например, отсутствием равновесия и неточностью определения температур в природных минеральных ассоциациях («природные» эмпирические калибровки); влиянием состава растворов и условий при экспериментальных калибровках; отсутствием критериев достоверности и точности применения эмпирических закономерностей в полуэмпирических калибровках (Horita, Clayton, 2007; Trail et al., 2009).

Квазигармоническое приближение и оценка влияние давления на изотопное фракционирование циркона

 

Таблица 3. Изменение частот колебаний (16O и 18O) при изменении объема ячейки циркона

T,°C

25–500

500–1000

1000–1100

1100–1750

Симметрия моды

∂ωi16/∂V

∂ωi18/∂V

∂ωi16/∂V

∂ωi18/∂V

∂ωi16/∂V

∂ωi18/∂V

∂ωi16/∂V

∂ωi18/∂V

A1g

-1.526

-1.438

-0.302

-0.284

-8.152

-7.685

-0.938

-0.885

A1g

-5.159

-4.863

-4.559

-4.298

-11.336

-10.686

-4.570

-4.308

A1u

-0.266

-0.250

-0.461

-0.435

-0.514

-0.485

-0.120

-0.113

A2g

1.178

1.110

0.841

0.793

1.129

1.065

0.689

0.650

A2u

-2.139

-2.066

-2.091

-2.048

-6.531

-6.264

-1.470

-1.423

A2u

-2.013

-2.001

-2.033

-2.025

-3.587

-3.530

-1.682

-1.682

A2u

-5.171

-4.914

-5.142

-4.883

-9.263

-8.831

-4.314

-4.088

B1g

-1.602

-1.593

-1.795

-1.788

-3.668

-3.654

-1.368

-1.365

B1g

-1.880

-1.853

-1.733

-1.698

-4.580

-4.542

-1.468

-1.443

B1g

-1.746

-1.631

-1.765

-1.673

-3.614

-3.265

-1.401

-1.294

B1g

-4.654

-4.483

-4.741

-4.553

-8.487

-8.208

-3.697

-3.574

B1u

4.711

4.441

3.781

3.564

8.902

8.391

2.820

2.659

B2g

1.440

1.358

0.950

0.895

3.082

2.905

0.667

0.629

B2u

-2.222

-2.095

-1.799

-1.696

-5.562

-5.243

-1.552

-1.463

B2u

-4.205

-3.964

-4.312

-4.065

-8.459

-7.974

-3.836

-3.616

Eg

0.454

0.486

0.629

0.736

0.428

0.522

0.389

0.503

Eg

-0.839

-0.898

-0.598

-0.755

-1.834

-1.979

-0.372

-0.521

Eg

-4.083

-3.907

-4.142

-3.974

-7.796

-7.452

-3.681

-3.530

Eg

-0.618

-0.690

-0.566

-0.636

-2.429

-2.505

-0.642

-0.717

Eg

-4.713

-4.454

-4.751

-4.485

-8.295

-7.760

-4.085

-3.846

Eu

-1.586

-1.400

-2.125

-1.922

-3.042

-2.757

-1.946

-1.798

Eu

-1.512

-1.462

-1.391

-1.351

-4.559

-4.311

-1.049

-0.993

Eu

-0.721

-0.931

-0.913

-1.095

-0.172

-0.592

-0.634

-0.787

Eu

-5.038

-4.753

-4.846

-4.570

-8.858

-8.321

-4.059

-3.830

Примечание. В таблице представлены результаты, соответствующие примитивной ячейке, Nq = 1.

 

В рамках гармонического приближения (HA) частоты колебаний не зависят от температуры и объема ячейки. Известные данные свидетельствуют об отсутствии существенного влияния ангармоничности на сдвиги частот при изотопных замещениях большинства минералов при умеренных температурах и давлениях (например, Polyakov, 1998; Schauble et al., 2004). Наши вычисления при различных объемах ячейки циркона, соответствующих тепловому расширению при 25, 500, 1000, 1100 и 1750°C (Mursie et al., 1992), демонстрирую сдвиги частот каждой колебательной моды в зависимости от объема, ∂ωi/∂V и эти данные использованы для оценки влияния ангармоничности на β-факторы циркона (в табл. 3 представлены результаты вычислений при Nq = 1). В рамках квазигармонического приближения (только когда частоты явно не зависят от T) для заданной температуры справедливы соотношения ωi = ω0i + (∂ωi/∂T)P(T – T0) (Gillet et al., 1989) и (∂ωi/dT)P = (∂ωi/dV)T × αV × V, где αV – объемный коэффициент теплового расширения. С поправкой на тепловое расширение до 1750°C (верхняя граница температуры обусловлена условиями эксперимента, Finch, Hanchar, 2003) в рамках квазигармонического приближения (QHA):

1000 lnβzrn (18O/16O: QHA) = 10.16937x – 0.19810x2 +  + 0.00321x3 (0 < T,°C < 1750). (3.1)

Эффект давления на β выражается (Polyakov, Kharlashina, 1994; Polyakov, 1998):

lnβPT=  γTKTlnβTV  , (3.2)

где γ – модальный параметр Грюнайзена, KT – изотермический модуль объемной упругости; (∂lnβ/∂T)V определяется дифференцированием соотношения (3.1). В качестве оценки KT принято значение 225 ГПa (Ozkan, 2008). Вместо модального параметра Грюнайзена γ для оценки зависимости β-фактора от давления принято значение экспериментально полученного теплового параметра Грюнайзена (γ = 1.14, Ozkan, Jamieson, 1978). Таким образом, влияние давления на β-факторы изотопного фракционирования кислорода циркона можно оценить посредством соотношений:

(∂(1000 lnβ)/∂P)T (кбар) ≈ 0.00865 × 106/T(K)2;

1000 lnβP – 1000 lnβ0 ≈ 0.00865 × 106/T(K)2 × P. (3.3)

Поскольку модальный параметр γ обычно меньше теплового параметра Грюнайзена (до 25%, см. Hofmeister, Mao, 2002), то приведенные оценки зависимости β-факторов от давления можно считать максимальными. При высоких давлениях около 100 кбар циркон становится неустойчивым и переходит в высокобарическую шеелитоподобную модификацию I41/a (реидит), поэтому, судя по соотношениям (3.3), эффект давления на изотопное фракционирование кислорода циркона при температурах более 750°С не превышает 0.15–0.20‰, а при T > 1000°С не превышает 0.05–0.10‰.

ВЫВОДЫ

Впервые методом «замороженных фононов» в рамках теории функционала плотности (DFT) определены β-факторы (18O/16O и 30Si/28Si) циркона в гармоническом и квазигармоническом приближениях с поправками на ограничения суммирования по волновым векторам и соответствующие факторы изотопного равновесия кварц–циркон для изотопов кислорода. β-факторы рассчитаны при дискретных температурах от 0 до 2000°С с шагом 20°С и интерполированы кубическими полиномами.

В рамках гармонического приближения, без учета изменения колебательного спектра и геометрических параметров кристаллической решетки с температурой:

1000 lnβzrn(18O/16O) = 9.83055x – 0.19499x2 + 0.00388x3,

1000 lnβzrn(30Si/28Si) = 7.89907x – 0.17978x2 + 0.00377x3.

В рамках квазигармонического приближения (QHA):

1000 lnβzrn(18O/16O: QHA) = 10.16937x – 0.19810x2 + 0.00321x3 (0 < T,°C < 1750).

Расчеты “из первых принципов” являются одним из важных методов оценки изотопных факторов фракционирования. Они подтверждаются сравнением с известными экспериментальными данными по неупругому рассеянию нейтронов и сравнением с результатами расчетов, использующих другие методы (например, методы динамики решетки). В качестве критериев достоверности вычислений и для оценки погрешностей результатов можно также использовать сравнение вычисленных кристаллографических параметров фаз.

Благодарности. Автор выражает признательность В.Б. Полякову (ИЭМ РАН) за большой труд по рецензированию статьи и ценные рекомендации. Расчеты β-факторов по данным неупругого рассеяния нейтронов выполнены с применением программы рецензента.

Источники финансирования. Работа выполнена при поддержке РФФИ (грант  № 19-05-00175). Вычислительные ресурсы предоставлены Ресурсным Центром «Вычислительный центр СПбГУ» (http://cc.spbu.ru).

D. P. Krylov

Institute of Precambrian Geology and Geochronology, Russian Academy of Sciences

Author for correspondence.
Email: dkrylov@dk1899.spb.edu

Russian Federation, St. Peterburg

  1. Bigeleisen J., Mayer M.G. Calculation of equilibrium constants for isotopic exchange reactions // J. Chem. Phys. 1947. V. 15. № 5. P. 261-267.
  2. Blanchard M., Poitrasson F., Méheut M., et al. Iron isotope fractionation between pyrite (FeS2), hematite (Fe2O3) and siderite (FeCO3): A first-principles density functional theory study // Geochim. Cosmochim. Acta. 2009. V. 73. № 21. P. 6565-6578.
  3. CRYSTAL14 : A program for the ab initio investigation of crystalline solids // Int. J. Quantum Chem. 2014. V. 114. № 19. P. 1287-1317.
  4. Clayton R.N., Kieffer S.W. Oxygen isotopic thermometer calibrations // Geochem. Soc. Spec. Publ. 1991. V. 3. P. 3-10.
  5. Chacko T., Cole D.R., Horita J. Equilibrium оxygen, hydrogen and carbon isotope fractionation factors applicable to geologic systems // Rev. Mineral. Geochem. 2001. V. 43. № 1. P. 1-81.
  6. Chaplot S.L., Pintschovius L., Choudhury N., Mittal R. Phonon dispersion relations, phase transitions, and thermodynamic properties of ZrSiO4: Inelastic neutron scattering experiments, shell model, and first-principles calculations // Phys. Rev. B. 2006. V. 73. № 9. P. 094308.
  7. Dawson P., Hargreave M.M., Wilkinson G.R. The vibrational spectrum of zircon (ZrSiO4) // J. Phys. C: Solid St. Phys. 1971. V. 4. № 2. P. 240-256.
  8. Dovesi R., Orlando R., Erba A., et al. CRYSTAL14 User’s Manual. 2016. 382 p.
  9. Evarestov R.A., Smirnov V.P. Symmetrical transformation of basic translation vectors in the supercell model of imperfect crystals and in the theory of special points of the Brillouin zone // J. Phys. Condens. Matter. 1997. V. 9. № 14. P. 3023-3031.
  10. Finch R.J., Hanchar J.M. Structure and chemistry of zircon and zircon-group minerals // Rev. Mineral. Geochem. 2003. V. 53. № 1. P. 1-25.
  11. Grimme S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction // J. Comput. Chem. 2006. V. 27. № 15. P. 1787-1799.
  12. Heyd J., Peralta J. E., Scuseria G. E., Martin R. L. Energy Band Gaps and Lattice Parameters Evaluated with the Heyd-Scuseria-Ernzerhof Screened Hybrid Functional // J. Chem. Phys. 2005. V. 123. № 17. P. 174101.
  13. Hoffbauer R., Hoernes S., Fiorentini E. Oxygen isotope thermometry based on a refined increment method and its applications to granulite-grade rocks from Sri Lanka // Precambrian Res. 1994. V. 66. № 1-4. P. 199-220.
  14. Hofmeister A.M., Mao Ho-kwang. Redefinition of the mode Grüneisen parameter for polyatomic substances and thermodynamic implications // PNAS. 2002. V. 99. № 2. P. 559-564.
  15. Horita J., Clayton R.N. Comment on the studies of oxygen isotope fractionation between calcium carbonates and water at low temperatures by Zhou and Zheng (2003; 2005) // Geochim. Cosmochim. Acta. 2007. V. 71. № 12. P. 3131-3135.
  16. Kieffer S.W. Thermodynamics and lattice vibrations of minerals: 5. Applications to phase equilibria, isotopic fractionation, and high-pressure thermodynamic properties // Rev. Geophys. 1982. V. 20. № 4. P. 827-849.
  17. Krylov D.P., Evarestov R.A. Ab-initio (DFT) calculations of corundum (α-Al2O3) oxygen isotope fractionation // Eur. J. Mineral. 2018. in press.
  18. Krylov D.P., Zagnitko V.N., Hoernes S., et al. Oxygen isotope fractionation between zircon and water: experimental determination and comparison with quartz-zircon calibrations // Eur. J. Mineral. 2002. V. 14. № 4. P. 849-853.
  19. Monkhorst H., Pack J. Special points for Brillouin zone integrations // Phys. Rev. B. 1976. V. 13. № 12. P. 5188-5192.
  20. Meheut M., Lazzeri M., Balan E., Mauri F. Equilibrium isotopic fractionation in the kaolinite, quartz, water system: Prediction from first-principles density-functional theory // Geochim. Cosmochim. Acta. 2007. V. 71. № 13. P. 3170-3181.
  21. Montanari B., Civalleri B., Zicovich-Wilson C.M., Dovesi R. Influence of the exchange-correlation functional in all-electron calculations of the vibrational frequencies of corundum (α-Al2O3) // Int. J. Quantum Chem. 2006. V. 106. № 7. P. 1703-1714.
  22. Mursi Z., Vogt T., Boysen H., Frey F. Single-crystal neutron diffraction study of metamict zircon up to 2000 K // J. Appl. Cryst. 1992. V. 25. № 4. P. 519-523.
  23. Nipko J.C., Loong C.-K. Inelastic neutron scattering from zircon // Physica B. Condensed Matter. 1997. V. 241-243. P. 415-417.
  24. Ozkan H., Jamieson J.C. Pressure dependence of the elastic constants of nonmetamict zircon // Phys. Chem. Mineral. 1978. V. 2. № 3. P. 215-224.
  25. Ozkan H. Correlations of the temperature and pressure dependencies of the elastic constants of zircon // J. European Ceramic Soc. 2008. V. 28. № 16. P. 3091-3095.
  26. Polyakov V.B. On anharmonic and pressure corrections to the equilibrium isotopic constants for minerals // Geochim. Cosmochim. Acta. 1998. V. 62. № 18. P. 3077-3085.
  27. Polyakov V.B., Kharlashina N.N. Effect of pressure on equilibrium isotopic fractionation // Geochim. Cosmochim. Acta. 1994. V. 58. № 21. P. 4739-4750.
  28. Porter A.R., Towler M.D. Needs R.J. Muonium as a hydrogen analogue in silicon and germanium; quantum effects and hyperfine parameters // Phys. Rev. B. 1999. V. 60. № 19. P. 13534-13546.
  29. Qin T., Wu F., Wu Z., Huang F. First-principles calculations of equilibrium fractionation of O and Si isotopes in quartz, albite, anorthite, and zircon // Contrib. Mineral. Petrol. 2016. V. 171. № 11. P. 1-14.
  30. Schauble E.A. First-principles estimates of equilibrium magnesium isotope fractionation in silicate, oxide, carbonate and hexaaquamagnesium(2+) crystals // Geochim. Cosmochim. Acta. 2011. V. 75. № 3. P. 844-869.
  31. Schauble E., Rossman G.R., Taylor H.P. Theoretical estimates of equilibrium chromium-isotope fractionations // Chem. Geol. 2004. V. 205. № 1-2. P. 99-114.
  32. Smyth J.R. Electrostatic characterization of oxygen sites in minerals // Geochim. Cosmochim. Acta. 1989. V. 53. № 5. P. 1101-1110.
  33. Sophia G., Baranek P., Sarrazin C., et al. First-principles study of the mechanisms of the pressure-induced dielectric anomalies in ferroelectric perovskites // Phase Transitions: A Multinational Journal. 2013. V. 86. № 11. P. 1069-1084.
  34. Trail D., Bindeman I.N., Watson E.B., Schmitt A.K. Experimental calibration of oxygen isotope fractionation between quartz and zircon // Geochim. Cosmochim. Acta. 2009. V. 73. № 23. P. 7110-7126.
  35. Valley J.W., Bindeman I.N., Peck W.H. Empirical calibration of oxygen isotope fractionation in zircon // Geochim. Cosmochim. Acta. 2003. V. 67. № 17. P. 3257-3266.
  36. Valley J.W., Cavosie A.J., Ushikubo T., et al. Hadean age for a post-magma-ocean zircon confirmed by atom-probe tomography // Nature Geoscience. 2014. V. 7. № 3. P. 219-223.
  37. Widanagamage I.H., Schauble E.A., Scher H.D., Griffith E.M. Stable strontium isotope fractionation in synthetic barite // Geochim. Cosmochim. Acta. 2014. V. 147. № 1. P. 58-75.
  38. Zheng Y.-F. Calculation of oxygen isotope fractionation in anhydrous silicate minerals // Geochim. Cosmochim. Acta. 1993. V. 57. № 5. P. 1079-1091.

Supplementary files

Supplementary Files Action
1. Fig. 1. Phonon density states of PDOS. Horizontal - phonon energy, vertical - density of states (arbitrary scale). Filled circles are the data on inelastic neutron scattering (Chaplot et al., 2006), the solid and dashed lines are the results of this work (for isotopologues of zircon 18O and 16O, respectively, with Nq = 27). View (66KB) Indexing metadata
2. Fig. 2. Comparison of calibrations of isotope fractionation of 18O / 16O between quartz and zircon. (a) Dependence of isotopic fractionation factors Δ = 1000 lnβqtz – 1000 lnβzrn on temperature, x = 106 / T2 (K). As 1000 lnβqtz, the values ​​taken are: 1000 lnβqtz = 12.55277x - 0.41976x2 + 0.01979x3 (Qin et al., 2016). (b) Deviations of temperatures obtained using various calibrations of quartz – zircon isotopic equilibrium from the results of this work. DFT - βzrn by the ratio (2.2), Δ = 2.72222x - 0.122477x2 + 0.01591x3; Zhe - using the “modified increments” method (Zheng, 1993), Δ = 0.72x + 4.26x1 / 2 - 1.79; VBP - “natural” calibration (Valley et al., 2003), Δ = 2.64x, Kie - calculations based on a set of semi-empirical rules for determining frequency shifts during isotopic substitutions (Kieffer, 1982), Δ = 3.1x; QWWH - ab-initio calculations by the method of "linear response" using pseudopotentials in the form of a superposition of plane waves (Qin et al., 2016), Δ = 2.68222x - 0.19011x2 + 0.01147x3. Experimental calibration Δ = 2.33x (Trail et al., 2009) in a limited range of 700–1000 ° C in scale. 2a coincides with the DFT and is shown in Fig. 2b (TBWS). View (319KB) Indexing metadata

Views

Abstract - 27

PDF (Russian) - 29

Cited-By


PlumX

Refbacks

  • There are currently no refbacks.

Copyright (c) 2019 Russian academy of sciences