Расчет теплоемкости и среднеквадратичных смещений по фононным спектрам для кристаллов с ОЦК и ГЦК решеткой Текст научной статьи по специальности «Физика»

УДК 534.1

РАСЧЕТ ТЕПЛОЕМКОСТИ И СРЕДНЕКВАДРАТИЧНЫХ СМЕЩЕНИЙ ПО ФОНОННЫМ СПЕКТРАМ ДЛЯ КРИСТАЛЛОВ С ОЦК И ГЦК РЕШЕТКОЙ

В.Е. Холодовский, И.О. Мачихина, Е.А. Кульченков

Исследуется динамика моноатомных кубических кристаллических решеток под воздействием межатомных сил, имеющих ван-дер-ваальсовскую природу. В адиабатическом приближении получена и решена система уравнений, описывающая колебания моноатомных ОЦК и ГЦК решеток, реализован принцип длинных волн, что позволило выразить силовые константы динамической модели через упругие константы рассматриваемого вещества и произвести расчеты фононных спектров, температурных зависимостей теплоемкости и среднеквадратичных смещений для некоторых кристаллов с ОЦК и ГЦК решеткой.

Ключевые слова: динамическая модель, диполь, кристаллическая решетка, упругие константы, фононный спектр, теплоемкость, среднеквадратичное смещение.

Введение

Как известно, динамические процессы, происходящие в веществе, так или иначе определяются тем, каким образом взаимодействуют между собой отдельные атомы.

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

Настоящая работа представляет собой продолжение исследований, начатых в [1] и [2], где была построена динамическая модель для ОЦК и ГЦК кристаллических решеток, использующая силы межатомного взаимодействия, имеющие ван-дер-ваальсовскую природу, и произведены расчеты дисперсионных кривых и фононных спектров для ряда элементов 1-5 групп таблицы Д.И. Менделеева без каких бы то ни было подгоночных параметров.

Согласно принципам, сформулированным в этих работах, атом кристалла рассматривается как структуризованный объект, состоящий из ионного остова и электронов на внешних оболочках. Считается, что остов колеблется как единое целое, а колебания электронов на внешних оболочках сводятся к колебаниям их центра заряда. Положение центра заряда внешних электронных оболочек (в.э.о.) атома А определяется взаимным расположением его остова и остовов его соседей, находящихся на первой и второй координационных сферах, и не обязано совпадать с положением остова атома А . В результате внутри атома наводится дипольный момент, действующий с некоторой силой на его остов. На остов атома действует сила, вызванная излучением диполей остальных атомов решетки, которая, изменяя положение центра заряда в.э.о., наводит в атоме А дополнительный дипольный момент, частично экранирующий эту силу. Атом А, представляющий собой динамический диполь, излучает электромагнитную энергию, которую можно рассматривать как работу силы реакции на его излучение. Как показано в [3], в первом приближении

сила реакции пропорциональна плечу диполя. В состоянии термодинамического равновесия на любом временном промежутке средняя энергия, поглощаемая атомом, совпадает со средней энергией, излучаемой им. Данное условие будет выполнено, если считать, что внешняя, частично экранированная, кулоновская сила уравновешивается силой реакции. Тогда движение остова атома будет происходить лишь под действием силы внутреннего диполя, наведенного соседними атомами и имеющего квантово-механическую природу своего возникновения.

В настоящей работе получено и решено уравнение динамики для ГЦК решетки с учетом фактора её плотной упаковки и малого порядка симметрии относительно оси вращения в направлении [110]. В отличие от работ [1, 2], здесь считается, что перекрытие электронных оболочек возможно лишь между атомами, лежащими друг относительно друга на первой координационной сфере, при этом степени перекрытия электронных оболочек при тангенциальных перемещениях атомов в направлении [100] и [110] отличаются друг от друга. С использованием принципа длинных волн определены силовые константы динамической модели и произведены расчеты фонон-ных спектров, а также температурных зависимостей теплоемкости и среднеквадратичных смещений. Используя результаты работ [1, 2], также рассчитаны температурные зависимости теплоемкости и среднеквадратичных смещений для некоторых элементов с ОЦК решеткой.

1. Уравнение динамики ГЦК решетки

Рассмотрим моноатомную кристаллическую решетку и обозначим через ц - массу остова каждого ее атома, через д - его заряд и пусть р = д2 /4пе0.

Пусть Л - какое-нибудь множество индексов, с помощью которого можно занумеровать все атомы решетки. Для каждого £ е Л обозначим через А^ соответствующий атом решетки, через

РС - узел, являющийся положением равновесия атома А^, а через и^ - смещение остова атома А| из положения равновесия в некоторый момент времени t. Обозначим, далее, через Б1 (С) -множество индексов из Л , нумерующих атомы решетки, находящиеся на /-той координационной сфере атома А^ . Пусть А^ - атом, соседний с атомом А^ . Перемещение остовов атомов А^ и А^

относительно друг друга вызывает изменение степени перекрытия орбиталей их в. э. о., что приводит к возникновению у этих атомов соответствующих дипольных моментов.

Обозначим через е ^ единичный направляющий вектор вектора РСРС , а через = Щ> - щ - вектор относительного перемещения остовов атомов А^ и А^. Пусть г ^ = е ^(е ^, ю) -радиальная, а г ^ = - г^ - тангенциальная составляющие вектора = , где в скобках обозначено скалярное произведение векторов е ^ и ю ^. В случае ГЦК решетки тангенциальная составляющая Гц > может быть разложена на два слагаемых т^< и >, где вектор т^< направлен вдоль одной из координатных осей, а вектор г2 > ортогонален векторам г ^ ^, и г^ >. В этом случае плечо дипольного момента р^,, наведенного в атоме А^ со стороны атома А^ -, лежащего на его первой координационной сфере, может быть определено формулой

' = к1г' + Ки' + K2tт1с' , (1)

где к1г, к1t, к21 - числовые параметры, постоянные для данного кристалла.

Плечо р полного дипольного момента, наведенного в атоме А со стороны всех его соседей, вычисляется путем суммирования по всем соседним атомам:

р с = X р с'. (2)

С 'С )

В формуле (2) удобно от суммирования по координационной сфере перейти к суммированию по полусфере. Для этого обозначим через Аатом, соседний к А^, противоположно расположенный по отношению к атому А . Тогда формула (2) может быть записана так:

Р, = 2 (Р, + Р, ), (3)

где S1(,) - какая-нибудь полусфера координационной сферы S1(,).

С учетом сказанного выше и работы [2], в состоянии термодинамического равновесия уравнение движения остова атома принимает вид

в

¡и, =--р, , (4)

а

где а - поляризуемость атома.

Будем считать, что рассматриваемый кристалл имеет форму куба, содержащего п3 элементарных кубических ячеек, и обозначим через a параметр решетки. Положим N = {1,2,...,2п} и зададим в пространстве систему кристаллографических координат Oxyz с единичными направляющими векторами ех, е , еz координатных осей так, чтобы положение каждого узла P = Pijk решетки могло быть задано по формуле

ОРук = |0'е x + je у + ке 2), (5)

где /, у,к е N - некоторый набор чисел. Обозначим через Л подмножество в N3, образованное всеми такими наборами (/, у,к), для которых формула (5) определяет узел решетки. Тогда для ГЦК решетки

Л = {(/, у,к) е N3 | сумма i + у + к нечетна }. Положение равновесия каждого атома А,>, находящегося на первой координационной сфере

атома А,, лежит в некоторой плоскости (100), проходящей через узел Р,. Тем самым, множество индексов 51(,) распадается на три составляющих Б1х (,), £1 у (,), Б1г (,), где, например, при е Б1х (,) узел Р, лежит в плоскости, ортогональной оси Ох. Пусть е Б1х (,), а х, у , - проекции вектора на координатные оси. Тогда т,, = хех, а формула (1) принимает вид

Р, = Кг Гц,' + КУЩ, ех +К2< (Щ ',,у еу + Щ^ - Гц,' ) = (К1г -К2< У,' + ех +К2< (Щ %у еу + Щ ех ). (6)

В других случаях аналогичные формулы получаются путем круговой перестановки индексов.

Подставляя (1) в (4), полагая а1г = вк1г /а , а1, = вк11 /а , а21 = вкъ/а и занося знак «-» под знак суммы, приходим к уравнению

= 2 [(°1 - Xе,,Щ,' >е, + ',хех + (щ,,',уеу + гег )] + ,' еSlx (,)

+ 2 [(СТ1г - Xе,,',' >е,,' + ,уеу + (щ,,',хех + ¿ег )] + (7)

,' е51у (,)

+ 2 [(°1 - Xе,,',>ей ' + ',zez + (щ,,',хех + ' ,уеу )]

Пусть , = (/, у, к) еЛ и = (г'', у', к' ) е £1(,). Положим = i' - i, е^> = у ' - у, е]Л! = к' - к, р = -Л. Тогда вектор е,, указывающий направление от узла р к узлу Р,', равен

е,, = (еН' ех + е уу' еу + £кк* ег

Выражая проекцию уравнения (7) на ось Ох и используя суммирование по полусферам, приходим к уравнению:

ли,,х = 2 ,х + ) + 1 2 {(СТ1г + °2г + х ) + (ст1г - )еii' екк > + )} +

+ 1 2 {(СТ1Г + СТ2г ',х + х ) + (а1г - СТ2г )еii 'е$' ,у + у )}.

2, 'еSlz (,) ' '

(8)

Проекции уравнения (7) на другие координатные оси получаются аналогично. ^04 Вестник ЮУрГУ, № 9, 2010

Принцип длинных волн

Рассмотрим уравнение (8), считая, что оно описывает колебания в длинноволновом диапазоне. В этом случае колебания решетки могут быть заданы функцией и(г, t) = щ ), которая плавно

меняется при изменении аргумента г = ОРС на величину порядка межатомного расстояния. Пусть Р ' - один из узлов, соседних к р . Положим Дг = Р С Р С,. Тогда щ, = и(г + Дг, t), ю с ' = и(г + Дг, t) - и(г, t). Для вычисления приращений проекций функции и(г, {) на координатные оси при переходе в соседний узел решетки можно воспользоваться формулой Тейлора до дифференциалов второго порядка. Считая момент времени t постоянным и опуская его в обозначениях, получим

Щ С ',х = их (г + Дг) - их (г) = Уых Дг + 2 < Дг,(У 'ых )Дг > , (9)

ЩС',х + ,х = <Дг,(У 'их)Дг> . (10)

где Уих - градиент функции их, а (V'их) - матрица, строки которой соответственно равны

удих Vдux

дх ду дг

Для произвольного узла решетки, положение которого определяется вектором г , вектор Дг смещения в положение соседнего узла равен

Дг = Дгй' = 2(£и'ех + еи'еу + £кк'ег) . Подставляя Дг в выражение (Дг,^' их )Дг>, приходим к равенству

.. , а2 2 д2их 2 д2их 2 д2и^ д2их

(ДГ (V их)Дг> = — {4 + 4' + 4' + ' еи' ^ +

4 дх ду дг дхду

д2их д2их

+ 2е Н'екк' ~ ~ + 2еи'екк' ,,,,}.

дхдг дудг

Подставляя равенства (10), (11) и подобные им в (8), приходим к уравнению

д2их а2 _ ч д2их „ ч„д2их д2их „ _ ч„д V д2иг

(11)

»Т1Г = V {2(а1г + ^ )—т + (^ 1г + 2°и + ^ +"ТТ") + 2(^1г - ^ + —^)}. (12)

дг 4 дх ду дг дхду дхдг

Согласно принципу длинных волн [4], уравнение (12) должно переходить в классическое уравнение распространения упругих волн деформаций в кристаллах. В проекциях на ось Ох для ГЦК решетки это уравнение имеет вид [5]:

,д Ч = а 2( а01 д2их + аС44 ( д 2их + д 2их ) + а(С12 + С44) ( + )

' дt2 ( 4 дх2 4 ( ду2 дг2 4 (дхду дхдг

где Сц, С12, С44 - упругие константы рассматриваемого вещества.

Сравнивая уравнение (12) с (13), приходим к соотношениям между коэффициентами динамической модели и упругими константами для ГЦК решетки:

°1г = "4(С11 + С12 + С 44) , = ^4 (2С44 - С11),а2г = -4(С11 - С12 - С 44 ) . (14)

Расчет фононных спектров и теплоемкости

Согласно общепринятым представлениям будем искать решение уравнения (4) динамики решетки в виде бегущих волн, заданных формулой

Щ = и С ( г, t) £ ,

где f = (i, j,k) еЛ , r = 2(iex + jey + kez)- радиус-вектор узла решетки, u^ = sin(Kr-mt) (или

u^ = cos(Kr - mt)), g = gxex + gyey + gzez - единичный вектор, указывающий направление поля-2п

ризации волны, а K = — (kxex + к ey + kzez) - волновой вектор. na

Пусть р - какой-либо узел решетки, rf - его радиус-вектор, а Pf и Pf - два соседних про-

a

тивоположно расположенных к нему узла. Пусть далее Ar^, = Pé P^ = = ^ Sex + s j< ey + skez) - вектор смещения из узла р в узел Pf. Тогда

2 KAra,

+ = (sin(K(r£ +Ar^,) - mt) + sin(K(r£ -Ar^,) - mt) - 2sin(Kr£ - mt)}g = -4sin —u^, (15) Проецируя равенство (15), например, на ось Ox, получим

+ , • 2 KArff л ■ 2 П(е» 'kx + Sjjky + Skkkz ) ',x + = -4sin uf gx = -4uf gx sin -2П- ' (16)

С другой стороны, также справедливо равенство

u f, x =-®2 uf gx . (17)

С учетом (16), (17) уравнение (8) принимает вид

2 . ^ .2 n(eii'kx + Sjj,ky + £kk'kz )

^m gx = 4 X °1tgx sin2- JJ У 1

>x ■ ^ w 1í&x 1 ,,

f'eSlx (f) 2n

, ^ r/ . , . . 2 n(siíkx +s jj'ky + skk ' kz) (i8) + 2 X{(ff> + CT2í )gx + (°1r -^2t )en, S kk g z }Sin ---+ , (18)

4'eS1y (f) 2n

+ 2 X {(CT1r + CT2í )gx + -^2í )sii' Zjj'gy )}sin

f eSiz (f)

. 2 n(sii'kx + Sjj'ky + Skk'kz )

Sin -2n-

Рассматривая проекции уравнения (7) на другие координатные оси, приходим к системе из трех линейных уравнений относительно неизвестных gx, gy gz, которая после соответствующих преобразований приводится к виду

cxgx - ^у - bygz = К - Х)gx

- ^х + ^у - bxgz = (^0 - Х)gy , (19)

- bygx - bxgy + = (^0 - Х)gz ,

где сх + су + cz = 0, а X = ¡лет2. Полученная система линейных уравнений является однородной и имеет симметрическую матрицу. Следовательно, ее собственные числа у = ст0 - X действительные, а собственные векторы g = gxех + gyеу + gzег, отвечающие различным собственным числам, ортогональны. Для нахождения собственных чисел матрицы системы (19) необходимо решить характеристическое уравнение

У3 + (СхСу + CyCz + CxCz - ь2х - Ь2у - Ь)у + Ь2хСх + Ь2уСу + ь1сг - СхСуСz + 2bxbybz = 0. (20)

Решая уравнение (20), мы можем найти зависимость частоты бегущей волны от набора чисел кх,ку,кг, определяющих величину и направление волнового вектора К . При этом наборы чисел

кх, ку, ^ должны выбираться так, чтобы были учтены все бегущие волны и ни одна из них не повторялась бы дважды. В случае ГЦК решетки можно считать, что кх =-п + 1,...,п,

ky = 0,...,n -1, тогда как kz = -

" n -1" n

_ 2 _ _ 2 _

В результате можем найти фононный спектр кри-

сталла.

Для вычисления фононного спектра удобно ввести сокращенные обозначения. Положим

nk nk

c(k) = cos —, s(k) = sin— . Тогда, с учетом формул (14), приходим к следующим выражениям n n

для коэффициентов уравнение (20):

cx = у (C11 - С44 )(2c(kz )c(ky ) - c(kx )c(ky ) - c(kx )c(kz )) ,

bx = a(C12 + C44 )s(ky )s(kz ) , (21)

a

= - (Cii + 2C44)(3 - c{kx)c{ky) - c(ky )c{kz) - c(kx)c(kz)),

а остальные коэффициенты вычисляются путем перестановки индексов x, y, z.

Выражения для коэффициентов cx, bx, a0 в случае ОЦК решетки приведены в [2]. Положим

- 3Р = cxcy + cycz + cxcz - bX - К - bZ , - q = bXcx + Кcy + b^cz - cxcycz + 2bxbybz ,

тогда уравнение (20) запишется в виде

Y3 - 3ру - q = 0, (22)

а его решения выразятся формулой

/— ю + 2тп Ym = 2yjр cos-3-, m = 0,1,2,

где ю = arccos —(или что то же ю = — - arctg , q =).

^ ^ 2 J4 Р3 - q2

2 -

Таким образом, частоты фононного спектра выражаются формулой

Ю = л/С^О^Тт)^ ■ (23)

Это позволило построить температурные зависимости теплоемкости и среднеквадратичных смещений атомов для ряда веществ с кубической решеткой согласно формулам:

дЕЛ 3Nah2 VT ehv7kTv2g(v)dv

Cv [dTJv kT2 o0 (ehv7kT -1)2 , (24)

Vm

г 2-1 h | g(v)dv (25)

[M ] = J ( hv / kT , (25)

4п / 0 v(e -1)

где N а - постоянная Авогадро, И - постоянная Планка, V = ю / 2п - циклическая частота моды колебаний, £(у)Су - вероятность обнаружения моды колебаний с частотой V в интервале [ V, V + Су ]■ Если частоту фононного спектра измерять в терагерцах, то формулы (24), (25) приводится к виду:

_ = 57439,5 е48т /Тт2£,(т)Ст

4 = Т2 0 (е48т/Т -1)2 , ( )

^ = (27)

т 0 т(е -1)

1 12

где т = 1^12 V, £1 (т) = 10 £(V), т - атомная масса.

На рис. 1-10 приведены фононные спектры для А1 и Си при 4 К, температурные зависимости теплоемкости и среднеквадратичных смещений для этих же элементов и для N и V при температуре 78 К, исходные данные для расчета которых взяты из работы [2]. Сравнение полученных кривых с экспериментальными данными из [6], [7], как это видно из приведенных рисунков, показывает хорошее соответствие теоретических кривых экспериментальным данным (экспериментальные данные нанесены точками).

8 ,0 0 2 4 6

V, ТГц V ТГц

Рис 1 ■ Фононный спектр А1 Рис. 2. Фононный спектр Си

Рис. 3. Температурная зависимость Рис. 4. Температурная зависимость

теплоемкости для N8 теплоемкости для V

Т , к

Т , к

Рис. 5. Температурная зависимость теплоемкости для Си

Рис. 6. Температурная зависимость теплоемкости для А!

С

о

о

3!

о

[u2],

10-22 м

T, K

Рис. 7. Среднеквадратичное смещение Na

[u2], 10-22 м

100

200

300

T, K

Рис. 8. Среднеквадратичное смещение V

[u2],

10-

2,4 1,8 1,2 0,6 0

100

200

300

T, K

Рис. 9. Среднеквадратичное смещение Al

[u2],

10-22 м21 " 0,8

0,6

0,4

0,2 0

100

200

T, K

Рис. 10. Среднеквадратичное смещение Cu

300

о

0

3

22 2 м

0

0

Статья выполнена при поддержке программы ФА по образованию «Развитие научного потенциала высшей школы» (грант РНП 2.1.1.7071).

Литература

1. Холодовский, В.Е. Дисперсионные соотношения для кубических кристаллических решеток в модели диполь- дипольных взаимодействий / В.Е. Холодовский, И.О. Мачихина, Е.А. Кульченков //Вестник ЮУрГУ. Серия «Математика, физика, химия».- 2009. - Вып.12.-№10(143).- С.92-99.

2. Холодовский, В.Е. Принцип длинных волн и фононные спектры кубических кристаллических решеток / В.Е. Холодовский, И.О. Мачихина // Вестник ЮУрГУ. Серия «Математика. Механика. Физика». - 2009. - Вып. 1. - № 22. - С. 109-116.

3. Холодовский, В.Е. Поток энергии и сила реакции на излучение подвижного диполя / В.Е. Холодовский, И.О. Сергеева // Вестник БГУ. Серия «Естественные и точные науки». - 2005. - Вып. 12.- № 4(273).- С. 266-268.

4. Борн, М. Динамическая теория кристаллических решеток / М. Борн, К. Хуан. - М.: ИЛ, 1958.-488 с.

5. Киттель, Ч. Введение в физику твердого тела / Ч. Китель. - М.: Наука, 1978. - 792 с.

6. Hultgren, R. Selected values of the thermodynamic prorerties of the elements / R. Hultgren, P. Desai // American society for metals, Metals Park, Ohio 44073.

7. Справочник. Свойства элементов. Часть 1. Физические свойства / под ред. Г.В. Самсоно-ва - М.: Металлургия, 1976. - 600 с.

Поступила в редакцию 24 декабря 2009 г.

CALCULATION OF HEAT CAPACITY AND MEAN-SQUARE DEVIATIONS ON PHONON SPECTRA FOR CRYSTALS WITH BCC AND FCC LATTICES

The paper considers the dynamics of monatomic cubic crystal lattices under the influence of the interatomic forces of Van der Waals nature. Using the adiabatic approximation, the authors received and solved a system of equations describing vibrations of monatomic bcc and fcc lattices, realized the longwave principle, which allowed to express the force constants of the dynamical model in terms of elastic constants of a given substance and to calculate phonon spectra and temperature dependences of heat capacity and mean-square deviations for some crystals with bcc and fcc lattices.

Keywords: dynamic model, dipole, crystal lattice, elastic constants, phonon spectrum, heat capacity, mean-square deviation.

Kholodovsky Vladimir Evgenievich is Cand.Sc. (Physics and Mathematics), Associate Professor of the Mathematical Analysis Department, I.G.Petrovsky Bryansk State University.

Холодовский Владимир Евгеньевич - кандидат физико-математических наук, доцент кафедры математического анализа, Брянский государственный университет им. академика И.Г. Петровского.

e-mail: tfbgubry@mail.ru

Machikhina Inna Olegovna is Post-Graduate Student at the Theoretical Physics Department , I.G.Petrovsky Bryansk State University.

Мачихина Инна Олеговна - аспирант кафедры теоретической физики, Брянский государственный университет им. академика И.Г. Петровского.

e-mail: ingibordit@yandex.ru

Kulchenkov Evgeny Aleksandrovich is senior lecturer of General Physics Department, I.G.Petrovsky Bryansk State University.

Кульченков Евгений Александрович - старший преподаватель кафедры общей физики, Брянский государственный университет им. академика И.Г. Петровского.

e-mail: evgeniy2000@mail.ru