Компьютерное моделирование диффузии взаимодействующих частиц на поверхности неупорядоченного сплава Текст научной статьи по специальности «Физика»

ФИЗИКА

Вестн. Ом. ун-та. 2012. № 4. С. 64-71.

УДК 53:621.38

И.В. Аполовникова, Г.Л. Бухбиндер

КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ДИФФУЗИИ ВЗАИМОДЕЙСТВУЮЩИХ ЧАСТИЦ НА ПОВЕРХНОСТИ НЕУПОРЯДОЧЕННОГО СПЛАВА

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

Ключевые слова: поверхностная диффузия; модель решеточного газа; Монте-Карло моделирование; диффузионный коэффициент.

Введение

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

В последние несколько десятилетий диффузия адсорбированных частиц в рамках модели решеточного газа интенсивно изучалась на однородных поверхностях. В различных приближениях было исследовано влияние межчастичных взаимодействий, а также роль фазовых переходов на диффузионный процесс [3-7, 12-17]. Однако реальные металлические поверхности, вообще говоря, неоднородны, и изучение влияния эффектов неоднородности на поверхностную диффузию и, в частности, на коэффициенты диффузии представляется достаточно важным.

Взаимодействие газов с гетерогенными металлическими поверхностями сильно чувствительно к адсорбционной энергетической топографии [18]. Влияние различных типов энергетической топографии на поведение коэффициентов диффузии изучалось в ряде работ [8-11; 18; 19]. В частности, в [10] на основе метода МК рассматривалась поверхностная диффузия невзаимодействующих адсорбированных частиц на энергетически гетерогенной металлической поверхности неупорядоченного бинарного сплава. Цель данной работы - на основе метода МК вычислить коэффициент диффузии в такой системе с учетом парных взаимодействий между диффундирующими частицами.

Модель решеточного газа

Рассмотрим простую кубическую решетку, узлы которой случайным образом занимаются металлическими атомами двух сортов А и В, форми-

© И.В. Аполовникова, Г.Л. Бухбиндер, 2012

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

Обозначим через xA концентрацию А атомов. При отсутствии вакансий в металлической матрице, концентрация В атомов равна xB = 1 - xA . В каждом междоузлии, окруженном атомами металла, может находиться только одна адсорбированная частица. Соседние междоузлия разделены потенциальным барьером. Чтобы совершить прыжок в новое положение, частица должна пройти через седловую точку потенциального барьера. Как и в [10], примем, что энергия в седловой точке еЗР фиксирована по всей решетке. Однако глубины потенциальных ям для различных типов междоузлий различны и будут зависеть от числа А и В атомов, находящихся в ближайшем окружении (рис. 1).

£5Р

_£?(0)

—£?(2)

— £-(3)

— £?(4)

Рис. 1. Адсорбционные состояния, имеющиеся в системе, и соответствующие им потенциальные барьеры в зависимости от числа А и В атомов, находящихся в ближайшем окружении междоузлия; • - А атомы, о - В атомы [10]

Пусть энергия парного взаимодействия между адатомами и ближайшими А и В атомами, соответственно, есть еА и еB

(еА ФеВ). Тогда энергию взаимодействия

адатома в г-ом междоузлии с металлической подложкой можно записать в виде:

, ПВ ) = ПАеА + Пвев , (1)

где nA,пВ - числа А и В атомов, находящихся

в ближайшем окружении /-го адсорбционного узла, для которых в случае квадратной решетки пА + пв = 4 . Таким образом, всего имеется пять энергетически различных адсорбционных положений. Отметим, что учитывается только взаимодействие между адатомом и металлическими атомами, локализованными на кристаллической поверхности. Удобно далее отсчитывать энергию от значения ев , т. е. примем, что ев = 0,

е0(ПА, ПВ ) = ПАеА =-ПАе01 (2)

где еА = — е0 < 0 и пА = 0...4 . В соответствии с

(2), высоты потенциальных барьеров равны Ае/ = £зр — . Самой глубокой потенци-

альной яме соответствует пА = 4, самой низкой - пА = 0 .

При наличии взаимодействия между адатомами полная энергия адатома в узле:

е / (ПА, ПВ ) = —ПАе0 + Е ' ЗуСу , (3)

где Зу - энергия парного взаимодействия ближайших адатомов (Зу > 0 соответствует отталкиванию, Зу < 0 - притяжению);

штрих у суммы означает, что суммирование ведется по ближайшим соседям г-го узла и [ 1, если узел / занят,

' [0, если узел / свободен.

Следует отметить, что взаимодействие между адатомам будет эффективно приводить к изменению высот потенциальных барьеров, которые в этом случае равны

Ае; =еэР — е, =еэР + ПАе0 — Е 'ЗуС!. (4)

Наконец, гамильтониан системы запишем в виде

(5)

где е{ дается равенством (2), З у = З0 и символ < у > обозначает суммирование по ближайшим соседям.

Поверхностные диффузионные коэффициенты

При экспериментальном изучении поверхностной диффузии измеряется так называемый химический диффузионный коэффициент, определенный формулой Кубо-Грина [2]:

Б =

(6)

где <(ЗМ)2 > - средний квадрат флуктуаций числа диффундирующих частиц N, Щ) = Еу /(1) - полный поток частиц и

V/ (^) - скорость / -й частицы в момент t.

Выражение (6) может быть представлено как [2]

Б = д(М1 квТ )

д 1пв З

(7)

здесь ц- химический потенциал, кВ - постоянная Больцмана, Т - температура, в -покрытие, БЗ - прыжковый диффузионный коэффициент:

Б, = 1/т

1 1

4М N

ЕаГ

(8)

где Аг. (t) - смещение / -й частицы в момент t из начального положения. Входящий в (7) термодинамический фактор может быть преобразован к виду [2]:

д(ц/кБТ ) Г((<^ У) ^

д 1п в N

(9)

Статистика адсорбционных положений

Прежде чем перейти к исследованию диффузионного процесса, были смоделированы две взаимопроникающие кристаллические подрешетки размером Ь х Ь с периодическими граничными условиями. Первая представляет собой кристаллическую решетку металлических атомов. Случайная конфигурация А и В атомов была создана в соответствии с их концентрациями хА и хВ . Междоузлия этой решетки формируют узлы подрешетки адатомных положений с разными энергиями.

Далее, случайным образом с заданным

покрытием в распределялись вЬ частиц по N = Ь узлам. После создания начальной конфигурации система приводилась в состояние термодинамического равновесия, используя метод Метрополиса для канонического ансамбля [1], при заданном значении е0 /квТ .

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

впА = КА /(NnA )тах ,

где Nп - число занятых междоузлий сорта пА, а (Nп )тах - общее число междоузлий данного типа при заданном значении хА. Вычисления были проведены для Ь = 64 , хА = 0,1, | З0 | = е0 в интервале температур от 100 К до 1500 К при е0 = 0,05 эВ [10]. Предполагается, что система достигнет равновесного значения, когда ее энергия начнет флуктуировать около среднего значения. Для рассматриваемых значений безразмерных параметров е0 /кВТ и размера решетки система приходит в равновесное состояние в течение приблизительно 5000 шагов Монте-Карло в случае отсутствия взаимодействия и в течение приблизительно 10000 шагов Монте-Карло с учетом взаимодействия между адсорбированными частицами.

а)

б)

Рис. 2. Зависимость впот покрытия в

при отталкивающем взаимодействии между адсорбированными частицами;

■ - в4, • - в3, ▲ - в2, ▼ — вх, ♦ — в0

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

Равновесные адсорбционные конфигурации для разных типов взаимодействий приведены на рис. 3 для некоторых значений параметров. Из рисунков видно, что притягивающее взаимодействие приводит к образованию кластеров (в сравнении с газом невзаимодействующих частиц). В случае отталкивания в области вблизи в = 0,5 и низких температурах возникают упорядоченные расположения частиц типа структуры с(2 х 2).

2

/=1

Рис. 3. Равновесное распределение адсорбированных частиц: а) при отсутствии взаимодействия, в = 0,5 и ео/квТ = 4,82;

б) притяжение, в = 0,5 и ео/квТ = 4,82;

в) отталкивание, в = 0,4 и £о/квТ = 4,82;

г) отталкивание, при в = 0,5 и ео/квТ = 9,62;

■ - занятые междоузлия, □ - свободные междоузлия

Термодинамический фактор

Вычисление термодинамического фактора (9), входящего в определение химического диффузионного коэффициента, было проведено в рамках метода МК для большого канонического ансамбля. Для этого моделируется случайная конфигурация адсорбированных частиц с произвольным покрытием, задаются две безразмерные величины е/квТ , ЦквТ и система приводится в состояние термодинамического равновесия. Затем для заданного значения ц/кВТ определяется среднее покрытие в = < N >/1} . Вычисленные таким образом зависимости ц/кВТ от в (адсорбционные изотермы) приведены на рис. 4 для некоторых значений температуры е0/кВТ и разных типов взаимодействий.

Из рисунка видно ступенькообразное

поведение ц/квТ при низких температурах е0/кВТ = 9,62 . Причина этого может быть в следующем. При низких температурах с ростом в идет поэтапное заполнение междоузлий вначале с самыми глубокими потенциальными ямами, а затем заполняются междоузлия с менее глубокими ямами. Пока идет заполнение потенциальных ям одинаковой глубины, химический потенциал меняется плавно, так как внесенные в систему новые частицы попадают в адсорбционные состояния с одинаковой энергией. Как раз

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

Термодинамический фактор

д(ц/кВТ)/Э^в) может быть найден в результате численного дифференцирования изотерм, представленных на рис. 4, либо

вычисления отношения < N >/<SN2 > .

ц/квТ

Рис. 4. Зависимость безразмерного химического потенциала у/квТ от покрытия в для Ха = 0,1 и Ео/квТ = 9,62 (во вставке показаны кривые для £о/квТ = 0,6);

■ - невзаимодействующая система,

• - отталкивающее взаимодействие, ▲ - притяжение

На рис. 5 изображены зависимости

< N >/< 8N2 > от покрытия в для разных типов взаимодействий. Качественно все кривые ведут себя подобным образом. Однако можно отметить, что в сравнении с невзаимодействующей системой отталкивающее взаимодействие приводит к увеличению термодинамического фактора, в то время как притяжение - к его уменьшению.

Из рисунка видно, что термодинамический фактор при низких температурах имеет характерные пики. Эти пики соответствуют ступеням в поведении ц/кВТ (рис. 4). С повышением температуры происходит значительное сглаживание пиков.

Рис. 5. Зависимость термодинамического фактора <М>/фМ¥> от покрытия в при ха = 0,1 и £о/квТ = 9,62 (во вставке показаны кривые для £о/квТ= 0,6);

■ - невзаимодействующая система,

• - отталкивающее взаимодействие, ▲ - притяжение

Вычисление диффузионных коэффициентов

Для вычисления коэффициентов диффузии используются выражения (7) и (8). Рассмотрим вначале вычисление прыжкового коэффициента диффузии БЗ .

Поверхностная диффузия моделируется в рамках метода МК, после того как система была приведена в статистическое равновесие при заданной температуре Т и покрытии в . Элементарным актом процесса является прыжок адатома в ближайший свободный узел. Если время жизни в г -м узле есть т, то вероятность того, что в течение промежутка времени Аt < т адатом совершит прыжок в соседнее вакантное положение, определяется отношением [2]:

Р =—,

(10)

где время жизни в потенциальной яме глубины Ае есть

квТ

(11)

и (01 - частота колебаний частицы в потенциальной яме, а Ае дается равенством (4).

Чтобы уменьшить время вычислений, в качестве единицы времени, соответствующей одному МК-шагу, возьмем интервал Аt из (10), равный времени жизни частицы в узле с минимальной глубиной потенциальной ямы:

\

(12)

Примем далее, что частоты 0 для всех ям приближенно равны [6]. Тогда, с учетом (11) и (12) выражение для прыжковой вероятности можно записать в виде:

Р=

■((Ае) ■ ) (

IV •) тт/ = ехр

'(Ае)

Ае-(Ае)

г У г / П

КТ

(13)

При таком выборе Аt частицы с единичной вероятностью за один МК-шаг будут преодолевать барьеры с наименьшей высотой. Отметим также, что вероятности (13) сохраняют принцип детального равновесия [6].

В рамках метода МК алгоритм моделирования диффузионного процесса состоит в следующем. Для /-ой частицы случайным образом выбирается один из соседних с ней узлов. Если узел занят, то переходим к следующей частице. Если свободен, то частице предоставляется возможность совершить туда прыжок. Для этого генерируется случайное число г, равномерно распределенное на интервале [0,1]. Если г < р, то частица занимает новое положение, в противном случае она остается на месте. При обходе таким образом всех адатомов будет совершен один МК-шаг, соответствующий времени Аt из (12).

Проделав т шагов Монте-Карло за время t = т (в единицах Аt), находим:

ЕАг/

= |Аг|2 = Ах2 + Ау2

(14)

Ах = ЕАх/, ау = Еау/

/=1

/=1

где Ах, =(Г+ — Г- )а , Ау =(Г+- Г- )а , ]±х и !

- число шагов каждой частицы соответственно вдоль осей - и у в положительном и отрицательном направлениях, а - длина прыжка.

Далее система адатомов несколько раз приводилась в состояние статистического равновесия при заданной температуре и покрытии с последующим диффузионным пробегом из т шагов МК, описанных выше, и вычислялось суммарное смещение

^г|2 согласно (14). Затем находилось среднее:

|А,

2

ЕС|Аг| к

< | Аг |2>= к=

К

где всего было произведено К выборок и индекс к относится к к -й выборке. Наконец, прыжковый диффузионный коэффициент находится из равенства:

1

Бз =■

•< | Аг |2> .

4 Nm

Результаты и обсуждение

Результаты моделирования диффузионного процесса представлены на рис. 6, 7. для К = 10 и т = 1-103.

2

г=1

На рис. 6 изображены зависимости

прыжкового диффузионного коэффициента БЗ от покрытия в для различных типов взаимодействий.

а)

б)

в)

Рис. 6. Зависимость прыжкового диффузионного коэффициента DJот покрытия в при ха=0,1 при различных температурах;

■ - невзаимодействующая система,

• - отталкивающее взаимодействие, ▲ - притяжение

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

При низких температурах (рис. 6 а, б) и малых покрытиях частицы занимают преимущественно узлы с глубокими потенциальными ямами (см. рис. 2). Вероятность прыжка из такой ямы мала, и частицы

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

При высоких температурах (рис. 6 в) заполнение потенциальных ям разной глубины происходит более равномерно (рис. 2 б). В этом случае подвижность частиц полностью определяется количеством вакантных мест и с ростом покрытия БЗ монотонно убывает. Для невзаимодействующей системы описанное выше поведение было выявлено в работе [10]. Как видно из рис. 6 качественно такое же поведение прыжкового диффузионного коэффициента сохраняется и при включении взаимодействия.

Что касается роли парных взаимодействий, то, как интуитивно и ожидалось, в сравнении с невзаимодействующей системой притяжение приводит к уменьшению прыжкового диффузионного коэффициента БЗ невзаимодействующих частиц, а отталкивание - к его увеличению во всей области изменения температуры (рис. 6). Следует отметить, что наибольшие изменения БЗ в зависимости от типа взаимодействия имеют место при низких температурах (рис. 6 а, б). Качественно все три кривые ведут себя подобным образом.

Фактором, увеличивающим подвижность частиц, в случае отталкивания, может являться образование при низких температурах, вблизи значения в = 0,5 упорядоченной фазы с(2 х 2) (при высоких температурах здесь наблюдается образования кластеров (рис. 3б)).

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

Зависимость химического коэффициента диффузии Б от покрытия представлена на рис. 7. Как следует из равенств (7) и

(9), его поведение полностью определяется

двумя факторами БЗ и < N >/<SN2 > . Общая тенденция сохраняется - притяжение уменьшает диффузионный коэффициент, а отталкивание приводит к его увеличению в сравнении с невзаимодействующей системой. При высоких температурах Б практи-

чески постоянен во всей области изменения в (рис. 7 в). Падение БЗ при высоких покрытиях компенсируется ростом термодинамического фактора.

6)

в)

Рис. 7. Зависимость химического коэффициента диффузии D от покрытия в при ха=0,1 при различных температурах;

■ - невзаимодействующая система,

• - отталкивающее взаимодействие, ▲ - притяжение

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

Заключение

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

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

Методом МК были вычислены прыжковый БЗ и химический Б диффузионные коэффициенты.

Энергетическая неоднородность поверхности в основном сказывается при низких температурах и покрытиях. Это связано с тем, что в этой области параметров адатомы в основном занимают самые глубокие потенциальные ямы, что ограничивает подвижность частиц. С ростом покрытия адатомы начинают занимать более мелкие потенциальные ямы, что приводит к монотонному росту БЗ на начальном этапе (рис. 6 б). Последующее поведение определяется количеством вакантных узлов.

При высоких температурах потенциальные ямы разной глубины занимаются равномерно, и подвижность частиц ограничена вакантными узлами, количество которых монотонно убывает с ростом покрытия. Это обстоятельство определяет и поведение БЗ (рис. 6 в).

Химический диффузионный коэффициент определяется произведением БГ и

< N >/<SN2 > (рис. 7). При больших покрытиях падение БЗ компенсируется ростом термодинамического фактора.

Описанное поведение характерно для всех типов взаимодействий.

Что касается роли взаимодействия, то притяжение уменьшает диффузионные коэффициенты БЗ и Б невзаимодействующей системы. Наоборот, отталкивание увеличивает их. Это сильнее проявляется при низких температурах. Фактором, ограничивающим подвижность частиц в случае притяжения, может являться образование кластеров (см. рис. 3б). В случае отталкивания причиной может являться образование при в = 0,5 упорядоченной с(2 х 2) фазы. Диффузионные прыжки могут нарушать периодическую структуру, а отталкивающее взаимодействие между адатомами будет стремиться ее восстановить. Это может привести к повышению подвижности всей системы и к увеличению диффузионных коэффициентов.

ЛИТЕРАТУРА

[1] Жданов В. П., Замараев К. И. Модель решеточного газа для описания хемосорбции на по-

верхности металлов // УФН. 1986. Т. 149. С. 635-669.

[2] Gomer R. Diffusion of adsorbates on metal surface // Rep. Progr. Phys. 1990. Vol. 53. P. 917-1002.

[3] Тарасенко A. А., Чумак А. А. Влияние отталкивания между адсорбированными атомами на их диффузию на поверхности кристалла // ФТТ. 1982. Т. 24. С. 2972-2976.

[4] Тарасенко A. А, Чумак А. А. Исследование диффузии в двумерной модели решеточного газа с сильным взаимодействием методом ре-нормгруппы // Поверхность. Физика, химия, механика. 1989. Т. 11. С. 98-105.

[5] Tringides M., Gomer R. A Monte-Carlo study of oxygen diffusion on the (110) plane of tungsten // Surf. Sci. 1984. Vol. 145. P. 121-144.

[6] Uebing C., Gomer R. A Monte-Carlo study of surface diffusion coefficients in the presence of adsorbate-adsorbate interactions. I. Repulsive interaction // J. Chem. Phys. 1991. Vol. 95. P. 7626-7635.

[7] Uebing C., Gomer R. Surface diffusion in presence of phase transitions. Monte Carlo studies of a simple lattice gas model // Surf. Sci.

1995. Vol. 331-333. P. 930-936.

[8] Ramirez-Pastor A. J., Nazzarro M. S., Riccardo J. L., Zgrablich G. Dimer physisorption on heterogeneous substrates // Surf. Sci. 1995. Vol. 341. P. 249-261.

[9] Uebing C., Pereyra V., Zgrablich G. Diffusion of interacting lattice gases on heterogeneous surface with simple topographies // Surf. Sci.

1996. Vol. 366. P. 185-192.

[10] Nieto F., Uebing C. Diffusion of adsorbates on random alloy surfaces // Eur. Phys. J. B. 1998. Vol. 1. P. 523-531.

[11] Bulnes F., Nieto F., Pereyra V., Zgrablich G, Uebing C. Energetic topography effects on surface diffusion // Langmuir. 1999. Vol. 15. P. 5990-5996.

[12] Tarasenko A. A., Nieto F., Jastrabik L., Uebing C. Adatom diffusion on square lattice. Theoretical and numerical studies // Eur. Phys. J. 2000. Vol. D12. P. 311-322.

[13] Tarasenko A. A., Nieto F., Jastrabik L., Uebing C. Collective surface diffusion of repulsively interacting particles on a triangular lattice: Comparison of real-space renormalization group and Monte Carlo approaches // Phys. Rev. B. 2001. Vol. 64. P. 075413-1-16.

[14] Tarasenko A. A., Jastrabik L., Muller T. Modeling diffusion on heterogeneous lattices: Derivation of general analytical expressions and verification for a two-dimensional square lattice // Phys. Rev. B. 2007. Vol. 75. P. 085401-11.

[15] Pawlikiwski K., Pekalski A. Surface diffusion of charged particles: Monte Carlo study // Phys. Rev. B. 2009. Vol. 79. P. 045419-1-10.

[16] Pinto O. A., Ramirez-Pastor A. J., Nieto F. Adsorption thermodynamics of a lattice-gas model with non-additive lateral interactions on triangular and honeycomb lattice // Physica. A. 2010. Vol. 389. P. 3456-3464.

[17] Garcia G. D., Sanchez-Varreti F. O., Bulnes F., Ramirez-Pastor A. J. Monte Carlo study of binary mixtures adsorbed on square lattice // Surf. Sci. 2012. Vol. 606. P. 83-90.

[18] Ricardo J. L., Chade M. A., Pereqra V. D., Zgrablich G. Adsorption and surface diffusion on generalized heterogeneous surface // Langmuir. 1992. Vol. 8. P. 1518-1531.

[19] Zgrablich G., at. al. Molecular processes on heterogeneous solid surface // Langmuir. 1996. Vol. 12. P. 129-138.