Промежуточная динамическая функция рассеяния и эффекты памяти в решеточных флюидах Текст научной статьи по специальности «Физика»

УДК 531.19

П. Аргиракис, профессор (Университет им. Аристотеля, г. Салоники, Греция);

Л. Скарпалезос, профессор (Университет им. Аристотеля, г. Салоники, Греция);

В. С. Вихренко, доктор физико-математических наук, профессор, заведующий кафедрой высшей математики (БГТУ)

ПРОМЕЖУТОЧНАЯ ДИНАМИЧЕСКАЯ ФУНКЦИЯ РАССЕЯНИЯ И ЭФФЕКТЫ ПАМЯТИ В РЕШЕТОЧНЫХ ФЛЮИДАХ

Моделированием по динамическому методу Монте-Карло выполнено детальное исследование промежуточной функции некогерентного рассеяния решеточного флюида с межчастичным притяжением ближайших соседей при температуре несколько выше критической. Предложено аналитическое выражение, совпадающее по форме с гидродинамическим пределом промежуточной функции некогерентного рассеяния, позволяющее воспроизвести результаты компьютерного моделирования. Это выражение включает коэффициент самодиффузии, зависящий от волнового вектора и времени. Времена релаксации группируются вблизи 10, 100 и 1000 шагов алгоритма Монте-Карло.

The self-intermediate scattering function of a lattice fluid with attractive nearest neighbor interpar-ticle interactions at a temperature slightly above the critical one is investigated in detail by means of Monte-Carlo simulations. An analytical expression is suggested to reproduce the simulation data. This expression is the generalization of that for the hydrodynamic limit with the wave vector, the time dependent tracer diffusion coefficient. The relaxation times group in three well separated regions around 10, 100 and 1000 Monte-Carlo steps and show weak concentration dependence.

Введение. Наряду с широко известным применением решеточных систем для изучения фазовых переходов они интенсивно используются в качестве простейших модельных систем при исследовании реальных физических объектов: субмонослоев на твердых поверхностях, ионных кристаллов, интерколяционных соединений, электронных ловушек и многих других. Важные приложения связаны с исследованием диффузионных процессов [1, 2].

Наиболее прямым методом изучения кинетических коэффициентов является теория линейного отклика [3-5]. С другой стороны, пространственно-временные функции распределения и полученные на их основе с помощью преобразований Фурье и Лапласа промежуточная функция рассеяния и динамический структурный фактор содержат ценную информацию о кинетических свойствах системы в широком диапазоне изменения пространственных и временных масштабов. Для различных жидкостей и флюидов эти характеристики подробно исследованы [4, 5], тогда как для решеточных флюидов таких результатов существенно меньше [6, 7].

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

Структура промежуточной функции некогерентного рассеяния. Пространственно-временная функция распределения Ван-Хова может быть представлена выражением

С (г, <) = ^ Е8(г + г (0) - г(Г ^ =

= С8(г, 0 + Ой(т, 0, (1)

где г - радиус-вектор; * - время; 8 - дельта-функция Дирака; угловые скобки обозначают усреднение по каноническому ансамблю из п частиц.

Суммирование выполняется по положениям всех частиц системы, определяемым радиус-векторами гг- и г;-. В решеточном флюиде каждая частица может занимать только один узел. Занятие узла несколькими частицами запрещено. Некогерентная часть функции распределения О. (г, *) включает положения одной и той же частицы ( = у) при двух различных моментах времени (0 и *), в то время как когерентная часть Gd (г, *) содержит члены с 1 Фу.

В гидродинамическом пределе, т. е. при больших расстояниях и временах, некогерентная часть функции распределения подчиняется уравнению диффузии

^Щг11 = '^(г,*)), (2)

от

где - тензор коэффициентов диффузии меченых атомов (коэффициент самодиффузии). Для решеток кубической симметрии тензор коэффициентов диффузии меченых атомов сводится к коэффициенту диффузии

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

представлены в виде С8(0, 0) = 1 и Gs(r Ф 0, 0) = 0, из соотношения (2) следует

Fs(k, t) = exp(-Dsк2t).

(3)

Коэффициент самодиффузии Д в гидродинамическом пределе не зависит от волнового вектора и времени. В этом случае ПФНР является гуссианом в к-пространстве с уменьшающейся с течением времени полушириной. Соотношение (3) позволяет вычислить коэффициент самодиффузии через ПФНР:

Ds =-

lnFs(k, t) k2t '

(4)

Поскольку динамика решеточного флюида описывается основным кинетическим уравнением, решение вида (3) и выражение (4) не могут быть использованы в области больших волновых векторов и малых времен. Для частиц на плоской квадратной решетке при рассмотрении некоррелированных прыжков в ближайшие соседние узлы [6] так же, как и в решеточном флюиде при пренебрежении влиянием эффектов памяти [8], к в соотношении (4) следует заменить функцией

П(к) = -

• 2 kxa . 2 kya sin2 —-— + sin2 2 2

(5)

Все иные вклады эффектов памяти могут быть включены в ^(к, 0 посредством зависимости от времени и волнового вектора коэффициента диффузии меченых атомов:

Fs(k, t) = exp[-Ds(k, t)n(k)t].

(6)

Таким образом, соотношение (4) может быть переписано как

Ds(k, t) = -

ln Fs(k, t) n(k)t '

(7)

Разумеется, при малых к соотношения (4) и (7) совпадают.

В соответствии с соотношениями (3) или (6) логарифм ПФНР экспоненциально уменьшается с течением времени. Эффекты памяти могут существенно изменить зависимость коэффициента диффузии от времени. Более того, в общем случае коэффициент диффузии, определяемый соотношением (7), зависит также от волнового вектора. В дальнейшем выражение (7) будет использовано для определения коэффициента диффузии как функции ^ и к.

Процедура моделирования. Для моделирования по динамическому методу Монте-Карло использована плоская квадратная решетка (ЬхЬ = 50x50 решеточных узлов) в сочетании с периодическими граничными условиями.

Первоначальная конфигурация соответствовала случайному распределению частиц по решетке с концентрацией с. В этом случае волновой вектор может быть введен как

к = Щ (к + т\), Ьа

m = 0,1, 2,...,L

h + m Ф 0.

(8)

Здесь 1 и ] - единичные векторы, направленные вдоль осей х и у соответственно.

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

p = exp

( - zJ >

KT

(9)

где У - параметр взаимодействия, связанный с критической температурой Тс соотношением У / квТС = 1,76; кв - постоянная Больцмана.

Один шаг алгоритма Монте-Карло (МКШ) определяется как п = схЬхЬ попыток совершить прыжок случайно выбранной частицей. Некоторое начальное число шагов алгоритма (~104 МКШ) отводится на переход системы в равновесное состояние, после чего в течение 2000 МКШ на каждом шаге вычисляется ПФНР. Отслеживая положение каждой частицы (х, у) в каждый момент времени, находим функцию К;(к, 0 для всех волновых векторов кх = 2кк /Ьа и ку = 2пт / Ьа, согласно соотношению

Fs(k, t) =

1 n (2п N

=n Xcos I —(h(-,t- - ,0)+m( y,t- y^-,o))

(10)

где к, т - две целочисленные независимые переменные, изменяющиеся от 0 до Ь / 2 (отрицательные значения этих величин не рассматриваются по соображениям симметрии). Координаты частиц определяются в единицах а. Полученные результаты усреднялись по 20 000 независимым траекториям.

Моделирование было выполнено для температуры Т = 1,2ТС, незначительно выше критической, и для концентраций с в диапазоне от 0,3 до 0,7 с шагом 0,05. Постоянная решетки а принималась равной 1.

Анализ результатов моделирования. Типичные результаты моделирования для коэффициента самодиффузии представлены на рис. 1. Он резко уменьшается в области малых к, но при приближении к к границе первой зоны Бриллю-эна скорость уменьшения снижается. Касательная

к рассматриваемой зависимости является горизонтальной при к = 0 и к = кт, где кт - максимальное значение волнового вектора (кт = п / а), а постоянная решетки а принимается равной 1.

0,009-

сО

0,006-

0,003

0,000

10 15 20 25

к

■ 1 ■ I I I 1 I 1 ■ ■ 1 1 ■ 1 ■ 1 ■ ■ 1

1,5 ' д А А Л -

3

1,0 ; 2 ....... _

)к(

сц : .•' 1

0,5 - 4 . ■

0,0 ................. . . . ...... ■

10

б

15

20

25

Рис. 1. Зависимость коэффициента самодиффузии О. (а) и функции В (б) от модуля волнового вектора к, направленного вдоль ребра решетки, при Т = 1,2 Тс и при различных временах моделирования: 1 - 20 МКШ; 2 - 200 МКШ; 3 - 2000 МКШ

Для аппроксимации данных моделирования разумно предложить следующее соотношение:

О. (к, *) = Д0(*)ехр[-В2(к, *)],

(11)

где В(0, *) = 0 и первая производная от В(к, *) по волновому вектору на границе первой зоны Бриллюэна также равна нулю.

Для нахождения аппроксимирующей функции для коэффициента диффузии вычислим В(к, *) = л1О.(к, *)/ О.0 (*). Типичный результат представлен на рис. 1, из которого следует, что эта функция экспоненциально возрастает, стремясь к некоторому предельному значению. Для обеспечения равенства нулю первой производной функции на границе первой зоны Бриллюэна используем слагаемое, пропорциональное волновому числу к, так что окончательное выражение запишем в виде

В(к, *) = Ь(X)Вт(*), Вт(*) = В(хт, *),

== Л. х=Ко, Хт = «У)'

Ь( х) = Ь0[1 - ехр(-х) - х ехр(-хт)], 1

(12)

Ь =

1 - (1 + хт)еХР(-хт ) '

где х - безразмерная величина; к - целочисленная переменная, изменяющаяся от 0 до 25; параметр к определяется в единицах (2п / Ьа). Ь(к, *) = Ь(х) - некоторая универсальная функция параметра х, зависящая от времени только через параметр к, Ь(0) = 0 и Ь(хт) = 1. Эта функция хорошо воспроизводит представленную на рис. 1 зависимость В(к, *) от волнового числа.

Предэкспоненциальный фактор в соотношении (12) определяется предельным значением

оА* )=кт о.(к, *).

к т0

(13)

Поскольку прямое вычисление при к = 0 невозможно, соответствующие результаты при к = ±1 и ±2 для определения Ок0(*) были аппроксимированы параболической зависимостью.

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

\

О.0 (*) = О„ + Ох ехр I — + О2 ехр

IТ1 )

( -* >

V Т2 )

. (14)

Две другие функции аппроксимируются трех-экспоненциальными зависимостями.

Все три времени релаксации оказываются разделенными на три разнесенных интервала: от 5 до 18, от 50 до 100 и от 500 до 1000 МКШ. Для каждого из рассматриваемых параметров разделение на интервалы даже более четко выражено и два последующие времени релаксации различаются на порядок величины.

Зависимость времен релаксации от концентрации является слабой. В отличие от времен релаксации предэкспоненциальные факторы показывают достаточно сильную зависимость от волнового числа.

Вклады релаксационных процессов на временах порядка десятков и сотен МКШ являются приблизительно равными, в то время как масштабный фактор сильно снижается за время порядка 5 МКШ. Общее снижение коэффициента

а

самодиффузии и функции Вт составляет 30 и 50% соответственно по сравнению с их длин-новолнововыми пределами. Уменьшение параметра к является существенно большим. Он уменьшается от 17 до 60 раз в зависимости от концентрации.

Глубина релаксации может быть оценена по отношению разницы между начальным и конечным (при ^ ^ те) значениями параметра к его конечному значению. Глубина релаксации коэффициента самодиффузии при к = 0 изменяется от 0,3 до 1,0 при росте концентрации от 0,3 до 0,7. Очень сильная релаксация наблюдается для масштабного параметра к. Он увеличивается приблизительно от 14 до 52 и незначительно снижается от 10 до 8 с ростом концентрации для волнового вектора, ориентированного вдоль ребра решетки и по диагонали соответственно. Таким образом, интенсивность релаксации масштабного параметра является сильно анизотропной в пространстве волновых векторов. Релаксация параметра Вт не зависит от концентрации и приводит к ее снижению приблизительно в 1,5 раза.

На рис. 2 гидродинамический предел коэффициента диффузии меченых атомов (при к = 0 и ^ ^ те) сравнивается с результатами моделирования по методу Монте-Карло, полученными посредством анализа среднего квадрата смещения в соответствии с выражением

((Ar)2) = 4DtJ.

(15)

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

Рис. 2. Зависимость коэффициента самодиффузии от концентрации с при температуре Т = 1,2 Тс

Значения, вычисленные в соответствии с выражением (15), систематически немного ниже по

сравнению с аналогичными величинами, наиденными на основе анализа ПФНР. Это связано с тем, что средний квадрат смещения частицы определялся по траектории в 10 000 МКШ, в то время как данные по ПФНР были получены по 2000 МКШ. Если уменьшить число шагов алгоритма Монте-Карло при вычислении среднего квадрата смещения до 2000 МКШ, различия в результатах двух подходов при больших концентрациях не превышают 2%.

Заключение. Результаты моделирования по методу Монте-Карло ПФНР решеточного флюида с межчастичным притяжением ближайших соседей представлены аналитическим выражением, которое включает коэффициент самодиффузии, зависящий от волнового вектора и времени. Параметры этого выражения аппроксимированы экспоненциальными функциями времени, что позволило выявить характерные масштабы времени, на которых протекают релаксационные процессы. Для коэффициента самодиффузии в длинноволновом пределе характерны два слабо зависящих от концентрации основных времени релаксации около 80 и 800 МКШ. Для двух других параметров три времена релаксации также разнесены примерно на порядок величины.

Литература

1. Gomer, R. Diffusion of adsorbates on metal surfaces / R. Gomer // Rep. Progr. Phys. - 1990. -Vol. 53. - P. 917-1002.

2. Allnatt, A. R. Atomic transport in solids / A. R. Allnatt, A. B. Lidiard. - Cambridge: University Press, 2004. - 600 p.

3. Zwanzig, R. Nonequilibium statistical mechanics / R. Zwanzig. - Oxford: Oxford University Press, 2001. - 222 p.

4. Boon, J. P. Molecular hydrodynamics / J. P. Boon, S. Yip. - New York: McGraw-Hill International, 1980. - 417 p.

5. Hansen, J.-P. Theory of simple liquids / J.-P. Hansen, I. R. McDonald. - London: Academic Press, 1986. - 418 p.

6. Kutner, R. Diffusion in concentrated lattice gases: Intermediate incoherent dynamical scattering function for tagged particles on a square lattice / R. Kutner, K. W. Kehr // Phys. Rev. B. - 1990. -Vol. 41. - P. 2784-2793.

7. Feng, E. H. A diagrammatic kinetic theory for a lattice model of a liquid: II. Comparison of theory and simulation results / E. H. Feng, H. C. Andersen // J. Chem. Phys. - 2004. - Vol. 121. -P. 3598-4005.

8. Statistical-mechanical description of diffusion in interacting lattice gases / G. S. Bokun [et al.] // Physica A. - 2001. - Vol. 296. - P. 83-105.

Поступила 02.03.2011