ЖЭТФ, 2021, том 160, вып. 6 (12), стр. 885-897

ЭФФЕКТИВНОЕ ТРЕНИЕ И ПОДВИЖНОСТЬ ГРАФЕНОВЫХ

НАНОЧАСТИЦ (НАНОЛЕНТ И НАНОТРУБОК)

НА ПЛОСКОЙ МНОГОСЛОЙНОЙ ПОДЛОЖКЕ h-BN

А. В. Савин^*

Федеральный исследовательский центр химической физики им. Н. Н. Семенова

Российской академии наук (ФИЦ ХФ РАН)

119991, Москва, Россия

Российский экономический университет им. Г. В. Плеханова

117997, Москва, Россия

Поступила в редакцию 12 июня 2021 г.,

после переработки 16 июля 2021 г.

Принята к публикации 25 июля 2021 г.

Методом молекулярной динамики с использованием двумерной цепной модели показано, что движение

графеновых наночастиц (нанолент и нанотрубок) на плоской термализованной многослойной подложке

h-BN описывается как движение частиц в вязкой среде с постоянным коэффициентом трения. Возни-

кающее при движении эффективное трение имеет волновую природу, причиной торможения является

взаимодействие наночастицы с тепловыми изгибными колебаниями листов подложки. Величина коэф-

фициента трения монотонно увеличивается с ростом температуры и уменьшается с увеличением размера

наночастицы. Для нанолент возникающее трение разделяется на два типа: на внутреннее и краевое (на

трение внутренней поверхности наноленты и трение ее краев с поверхностью подложки). При длинах

L < 35 нм главную роль в торможении движения играет краевое, а при L > 35 нм — внутреннее трение.

Под действием постоянной продольной силы динамика наночастиц всегда выходит на режим движения с

постоянной скоростью, значение которой прямо пропорционально величине силы и обратно пропорцио-

нально значению коэффициента трения. Моделирование движения наноленты при наличии нормальной

нагрузки (давления) показывает, что рост нагрузки приводит к уменьшению внутреннего трения из-за

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

краевого трения из-за вдавливания краев наноленты в подложку. В силу этого эффект уменьшения тре-

ния при увеличении нормальной нагрузки может наблюдаться только для нанолент достаточно большого

размера (L > 250 нм), когда главную роль играет внутреннее трение.

DOI: 10.31857/S0044451021120105

слоистым материалам, которые могут демонстриро-

вать различные новые физические свойства по срав-

нению с их однородными аналогами [15-17]. Так,

1. ВВЕДЕНИЕ

было показано, что использование гетероструктур

G/h-BN позволяет получить нужные электронные

Двумерные (2D) слоистые материалы, такие как

свойства [18,19], а также существенно понизить тре-

графен (G), гексагональный нитрид бора (h-BN),

ние между слоями [20].

дисульфид молибдена (MoS₂) и вольфрама (WS₂),

вызывают большой интерес из-за своих уникаль-

Важной задачей для механических устройств

ных электронных [1-3] и механических [4-7] свойств.

нано- и микроразмеров является максимально воз-

Благодаря очень низкому трению слоев эти матери-

можное уменьшение трения [21]. Стандартные схе-

алы могут быть использованы в качестве высоко-

мы смазки на таких размерах перестают работать,

эффективной сухой смазки [8-14]. В последнее вре-

здесь необходимо от жидкой смазки перейти к су-

мя повышенное внимание уделяется гетерогенным

хой, связанной со скольжением плоских молеку-

лярных слоев. Такой подход, впервые теоретически

^* E-mail: asavin@center.chph.ras.ru

предложенный несколько десятилетий назад [22],

885

А. В. Савин

ЖЭТФ, том 160, вып. 6 (12), 2021

тов трения от размера наночастиц и температуры

подложки.

Работа построена следующим образом. В разд. 2

строится цепная 2D-модель, используемая далее для

моделирования движения графеновых наночастиц

(нанолент и нанотрубок) по плоской многослойной

подложке h-BN. В разд. 3 проводится моделирова-

ние торможения свободного движения наночастиц

в результате их взаимодействия с термализован-

ной подложкой. Определяется коэффициент возни-

Рис. 1. Прямоугольный лист графена размера 2.02 ×

кающего эффективного трения и анализируется его

× 1.88 нм² (160 атомов углерода), лежащий на плоской

природа. В разд. 4 анализируется подвижность на-

поверхности кристалла h-BN

ночастиц, моделируется их движение под действи-

ем боковой силы. В разд. 5 проводится провер-

ка выполнения для нанолент эмпирического закона

позволил достичь чрезвычайно низких коэффициен-

Амонтона - Кулона (моделируется влияние на тре-

тов трения [23-25]. Особенно сильно понизить тре-

ние нормальной нагрузки). Заключительные заме-

ние позволяет использование двумерных материа-

чания приводятся в разд. 6.

лов, таких как графен и гексагональный нитрид

бора h-BN [26]. Слоистые структуры из этих ма-

териалов могут обладать сверхскольжением слоев

2. МОДЕЛЬ

[15, 20, 27].

Использование слоистых 2D-структур требует

Для описания динамики нанолент и нанотрубок

фундаментального понимания истоков силы трения

графена на плоской поверхности многослойной под-

на атомном уровне. В работах [20, 27, 28] проводи-

ложки удобно использовать двумерную модель сис-

лось моделирование медленного протягивания слоя

темы молекулярных цепей [31-33]. Если считать,

графена, расположенного над слоем h-BN, под дей-

что наночастица (нанолента или нанотрубка графе-

ствием боковых сил. Величина трения оценивалась

на) и листы h-BN подложки лежат так, что у всех

как максимальное значение возникающей силы, пре-

направление зигзаг совпадает с осью x (см. рис. 1),

пятствующей скольжению. Фактически находилась

то двумерная цепная модель будет описывать сече-

сила трения покоя при нулевой температуре. Функ-

ние системы наночастица + многослойная подложка

ционирование механических устройств нано- и мик-

вдоль оси x (см. рис. 2). Тогда одной частице в дву-

роразмеров невозможно описать без учета тепловых

мерной модели будут соответствовать все атомы на-

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

ноленты (нанотрубки), имеющие одинаковые коор-

колебаний на трение имеет особое значение. В от-

динаты x, z.

личие от стандартных сценариев трения, где тре-

Если атомы, расположенные вдоль одной линии,

ние всегда уменьшается с повышением температуры

параллельной оси y, двигаются синхронно, меняя

[29], в двухслойной структуре G/h-BN трение может

только координаты x, z, то гамильтониан одной на-

увеличиваться с повышением температуры [30].

ноленты (нанотрубки) графена (h-BN) будет иметь

Целью настоящей работы является объяснение

вид гамильтониана цепи, расположенной в плоскос-

на атомарном уровне возникающего трения при дви-

ти xz:

жении наночастиц (нанолент и нанотрубок графена)

[

]

∑

по плоской поверхности многослойной подложки. В

H_i =

M_i( u_n, u_n) + V_i(R_n) + U_i(θ_n)

(1)

качестве подложки будет использована плоская по-

n=1

верхность гексагонального кристалла нитрида бора

h-BN (см. рис. 1). Для определения коэффициен-

Здесь индекс i = 1, если рассматривается нано-

тов возникающего эффективного трения будет про-

лента графена (G), и i = 2, если рассматривает-

моделировано торможение свободного движения на-

ся нанолента бор нитрида (BN). Двумерный век-

ночастиц вследствие их взаимодействия с термали-

тор u_n = (x_n, z_n) задает координаты n-ой части-

зованной подложкой (вследствие перекачки началь-

цы цепи. Масса частицы для цепи G совпадает с

ной кинетической энергии наночастиц в подложку).

массой атома углерода M₁ = M_C = 12m_p, а для

Будет проанализирована зависимость коэффициен-

цепи BN — со средней массой атомов бора и нит-

886

ЖЭТФ, том 160, вып. 6 (12), 2021

Эффективное трение и подвижность графеновых наночастиц. ..

cos(θ_n) = -(vn-1, v_n)/Rn-1R_n,

вектор v_n = un+1 - u_n.

Параметры потенциалов (2), (3) для цепи G опре-

делены в [31, 32] из анализа дисперсионных кривых

наноленты графена. Продольная жесткость K₁ =

√

3/2

= 1.228Å

= 405 Н/м, шаг цепи R₁ = r_CC

(r_CC = 1.418Å — длина валентной связи C-C в ли-

сте графена), энергия ϵ₁ = 3.5 эВ.

Нанолента h-BN имеет такую же структуру, как

и нанолента графена. В двумерной модели ее га-

мильтониан тоже будет иметь вид (1). Для нахож-

дение параметров цепи BN была рассмотрена нано-

лента h-BN. Взаимодействие атомов наноленты опи-

сывалось расширенным потенциалом Терсоффа для

нитрида бора [34]. Вычисления показали, что нано-

лента h-BN имеет такую же гексагональную струк-

туру, как и нанолента графена. В основном состоя-

нии длина валентной связи B-N r_BN = 1.445685Å

незначительно превышает длину валентной связи

C-C. Анализ дисперсионных кривых наноленты по-

казал, что она имеет частотный спектр изгибных

Рис. 2. Вид стационарных состояний однослойных листов

(внеплоскостных) колебаний 0 ≤ ω ≤ 747 см-1 и

графена длины L = 2.33, 4.79, 9.70 нм (число звеньев це-

спектр продольных (вутриплоскостных) колебаний

пи Nc = 20, 40, 80) и однослойных нанотрубок с индексом

0 ≤ ω ≤ 1688 см-1. Вычисления показали, что дис-

киральности (10, 0), (20, 0), (40, 0) (число звеньев цикли-

персионные кривые двумерной цепи с гамильтониа-

ческой цепи Nc = 20, 40, 80) на многослойной подложке,

образованной поверхностью кристалла h-BN (части a, b,

ном (1) наилучшим образом соответствуют диспер-

c и d, e, f). Число слоев листов h-BN K - 1 = 3, число

сионным кривым полноатомной модели наноленты

звеньев в BN-цепи N_bn = 360. Черная прямая показывает

при жесткости K2 = 480 Н/м, шаге цепи R2 =

√

положение поверхности неподвижной подложки

= r_BN

3/2 = 1.252Å и энергии ϵ₂ = 1.10 эВ (см.

рис. 3).

Отметим, что гамильтониан цепи (1) дает энер-

рида M₂

= (M_B + M_N )/2 = 12.4085m_p (m_p =

гию деформации наноленты, приходящуюся на про-

= 1.6603 · 10-27 кг — масса протона).

√

дольную полосу ширины Δy = R_i/

3, поэтому, если

Отметим, что используемая цепная 2D-модель

энергию системы цепей далее будем нормировать по

(1) не позволяет описать поперечные (вдоль оси y)

наноленте графена, то энергию деформации нано-

колебания наноленты, но позволяет описать про-

лент h-BN нужно умножить на нормирующий мно-

дольные (вдоль оси x) и изгибные (вдоль оси z) ко-

житель c = R₁/R₂ = r_CC/r_BN = 0.9808.

лебания наноленты.

Для вычисления эффективного потенциала

Потенциал

невалентного взаимодействия узлов цепей (энер-

гии взаимодействия одного атома с поперечной

V_i(R) =

K_i(R - R_i)²

(2)

линией атомов) были использованы потенциалы

описывает продольную жесткость цепи, K_i — жест-

Леннарда - Джонса

кость взаимодействия, R_i — равновесная длина свя-

[

]

зи (шаг цепи), R_n = |un+1 -u_n| — расстояние между

V (r) = ϵ₀

(r₀/r)¹² - 2(r₀/r)⁶

(4)

соседними узлами n и n + 1.

Потенциал

со значениями энергии ϵ₀ и равновесными длинами

взаимодействия r₀, приведенными в таблице.

U_i(θ) = ϵ_i[1 + cos(θ)]

(3)

Проведенные вычисления показали (см. рис. 4),

описывает изгибную жесткость цепи, θ — угол меж-

что взаимодействия узлов цепей, соответствующих

ду двумя соседними связями,

слоям бор-нитридной подложки и слою графена,

887

А. В. Савин

ЖЭТФ, том 160, вып. 6 (12), 2021

Рис. 4. Вид потенциалов взаимодействия узлов цепей

Wi(r), i = 1, 2 (кривые 1, 2). Сплошные линии дают за-

висимости, полученные численно, штриховые линии — за-

висимости, описываемые потенциалом Леннарда - Джон-

са (5)

ε₁ = 0.01511 эВ, ε₂ = 0.01433 эВ, равновесные рас-

стояния r₁ = 3.642Å, r₂ = 3.701Å.

При моделировании динамики многослойной

подложки необходимо ограничить число слоев.

Поэтому будем считать, что первый (самый ниж-

Рис. 3. Вид дисперсионных кривых наноленты h-BN (крас-

ний) слой уже взаимодействует с неподвижной

ные линии 1 соответствуют поперечным, а синие 2 — про-

плоской поверхностью кристалла (на рис.

эта

дольным колебаниям наноленты). Штриховые кривые 3 и

поверхность показана черной линией). Энергия

4 соответствуют изгибным и продольным колебаниям цеп-

взаимодействия атомов слоев с неподвижной под-

ной 2D-модели наноленты

ложкой тоже вычислялась с помощью потенциалов

Леннарда - Джонса (4).

Таблица. Значения параметров потенциала Лен-

Вычисления показали (см. рис. 5), что потенци-

нарда- Джонса (4) для различных пар взаимодей-

ал взаимодействия с неподвижной подложкой W₀(h)

ствующих атомов [35]

(зависимость энергии от расстояния h узла цепи до

плоскости подложки) с высокой точностью может

быть описан (k, l) потенциалом Леннарда - Джонса

BB NN BN CB CN

ϵ₀, МэВ

7.81

2.99

4.81

5.96

3.69

W₀(h) = ε₀[k(h₀/h)^l - l(h₀/h)^k]/(l - k),

(6)

r₀,Å

4.083

3.660

3.866

3.976

3.756

где степень l = 10, k = 3.75, энергия взаимодействия

ϵ₀ = 0.0974 эВ, равновесное расстояние h₀ = 3.49Å.

Рассмотрим K-слойные структуры, представлен-

можно с высокой точностью описать потенциалом

ные на рис. 2. Пусть первые k = 1, . . . , K - 1 сло-

Леннарда - Джонса

ев соответствуют BN-цепям (слоям кристалла h-

BN), состоящим из N_bn звеньев. Эти слои лежат

[

]

W_i(r) = ε_i

5(r_i/r)¹¹ - 11(r_i/r)⁵

/6.

(5)

на плоской твердой подложке и взаимодействуют с

ней (будем считать, что поверхность твердой под-

Здесь r — расстояние между взаимодействующи-

ложки совпадает с плоскостью z = 0). Последняя

ми узлами (индекс i = 1, если описывается взаи-

K-тая G-цепь, имеющая N_c < N_bn звеньев, соответ-

модействие узлов цепей BN, i = 2 — взаимодей-

ствует графеновой наночастице (наноленте, нано-

ствие узлов цепей BN и G). Энергии взаимодействия

трубке), лежащей на плоской деформируемой мно-

888

ЖЭТФ, том 160, вып. 6 (12), 2021

Эффективное трение и подвижность графеновых наночастиц. ..

вектор v_n,k = un+1,k -u_n,k, расстояние между узла-

ми разных цепей k₁ и k₂ rn,k1;l,k2 = |ul,k2 - un,k1 |.

В формуле потенциальной энергии первое слагае-

мое (8) задает энергию деформации BN-цепей с уче-

том энергии их взаимодействия с твердой подлож-

кой, второе слагаемое (9) — энергию невалентного

взаимодействия BN-цепей, третье слагаемое (10) —

энергию деформации G-цепи, а последнее слагаемое

(11) — энергию невалентного взаимодействия G-це-

пи с BN-цепями.

3. ТОРМОЖЕНИЕ СВОБОДНОГО

ДВИЖЕНИЯ НАНОЧАСТИЦЫ

Рис. 5. Потенциал взаимодействия W0(h) узла цепи BN

Рассмотрим систему, состоящую из K - 1 = 3

c неподвижной плоской поверхностью кристалла h-BN.

слоев BN-цепей с N_bn = 2400 звеньями. Длина такой

Сплошная линия дает зависимость, полученную численно,

трехслойной подложки L = (N_bn - 1)R₂ = 300.4 нм.

штриховая линия — зависимость, описываемую потенциа-

В качестве графеновой наночастицы, расположен-

лом Леннарда - Джонса (6)

ной на этой подложке, возьмем линейную G-цепочку

из N_c = 20, 40, 80, 160, 320 звеньев (длина соответ-

ствующей наноленты графена L_c = 2.3, 4.8, 9.7, 19.5,

гослойной (K-1-слойной) подложке. Координаты уз-

39.2 нм). Также в качестве графеновой наночасти-

лов этой системы K цепей задаются векторами

цы будут рассмотрены циклические цепочки (угле-

{u_n,k = (x_n,k, z_n,k)}Nk,Kn=1,k=1,

родные нанотрубки) из N_c = 20, 40, 80 звеньев (см.

рис. 2d,e,f). Промоделируем свободное движение на-

где N_k — число узлов в k-й цепи (N_k = N_bn при

ночастицы по такой многослойной подложке.

k = 1,...,K - 1 и N_K = N_c).

Сначала поместим наночастицу у левого края

Гамильтониан системы цепей будет иметь вид

подложки, как показано на рис. 2. Далее найдем ста-

ционарное состояние системы «наночастица + трех-

∑

^∑1

H =

cM₂( u_n,k, u_n,k)+

слойная подложка». Для этого численно методом со-

k=1 n=1

пряженных градиентов решим задачу на минимум

∑

энергии:

M₁( u_n,K, u_n,K) + E,

(7)

E → min : {u_n,k}Nk,

(12)

n=1,k=1

n=1

где число цепей K = 4, число узлов N₁ = N₂ = N₃ =

где потенциальная энергия

= N_ab = 2400, N₄ = N_c. Типичный вид стационар-

∑

ных состояний показан на рис. 2. Как видно на ри-

E =c

[V₂(R_n,k)+U₂(θ_n,k)+W₀(z_n,k)] +

(8)

сунке, нанотрубка графена в результате взаимодей-

k=1 n=1

ствия с подложкой с увеличением ее радиуса (чис-

∑

ла звеньев N_c) начинает приобретать сплющенную

W₁(rn,k1;l,k2 ) +

(9)

каплевидную форму.

k1=1 k2=k1+1 n=1 l=1

Для получения термализованного состояния сис-

∑

темы поместим ее в термостат Ланжевена. Для

[V₁(R_n,K ) + U₁(θ_n,K )] +

(10)

n=1

этого зафиксируем x-координаты краевых частиц

BN-цепей (положим скорости { x1,k ≡ 0, xNab,k ≡

∑

W₂(rn,k;l,K).

(11)

≡ 0}K-1k=1) и x-координату центральной частицы

k=1 n=1 l=1

G-цепи ( xNc/2,K ≡ 0) и численно проинтегрируем

систему уравнений Ланжевена

Здесь расстояние между соседними узлами k-й цепи

R_n,k = |un+1,k - u_n,k|, косинус угла между двумя

∂H

соседними связями

M₂ ü_n,k = -

- ΓM₂ ˙u_n,k + Ξ_n,k,

c ∂u_n,k

(13)

cos(θ_n,k) = -(vn-1,k, v_n,k)/Rn-1,kR_n,k,

k = 1,...,K - 1, n = 1,...,N_bn,

889

ЖЭТФ, вып. 6 (12)

А. В. Савин

ЖЭТФ, том 160, вып. 6 (12), 2021

∂H

M₁ ü_n,K = -

- ΓM₁ ˙u_n,K + Ξ_n,K ,

∂u_n,K

(14)

n = 1,...,N_c,

где Γ = 1/t_r — коэффициент трения (время релакса-

ции t_r = 0.4 пс), Ξ_n,k = (ξn,k,1, ξn,k,2) — двумерный

вектор нормально распределенных случайных сил,

нормализованных условиями

〈ξn,k1,i(t1)ξl,k2,j (t2)〉 = 2MΓkBT δnlδk1k2 δij δ(t2-t1)

— температура термостата, kB — постоянная

Больцмана, масса M = M₂/c при k < K и M = M₁

при k = K).

Возьмем начальные условия для системы урав-

нений движения Ланжевена (13), (14) соответствую-

щими основному стационарному состоянию системы

{u_n,k(0) = u0n,k,

u_n,k(0) = 0}Nk,Kn=1,k=1

и численно проинтегрируем систему в течение вре-

мени t₀ = 20t_r. За это время система придет в пол-

ное равновесие с термостатом, и мы получим терма-

лизованное состояние системы

{w_n,k = u_n,k(t₀), v_n,k = u_n,k(t₀)}Nk,Kn=1,k=1.

Рис. 6. Зависимость центра тяжести графеновой наночас-

тицы xc от времени t при температуре подложки T =

Для моделирования свободного движения нано-

= 300 K. Тонкие кривые показывают 9 траекторий дви-

частицы по термализованной многослойной подлож-

жения для различных независимых начальных реализаций

ке оставим фиксированными только x координаты

термализованного состояния системы. Сплошная толстая

краевых узлов BN-цепей, уберем взаимодействие си-

(черная) линия дает траекторию среднего значения центра

стемы с термостатом, а всем атомам наночастицы

тяжести xc, полученную по 256 независимым реализациям

термализованного состояния. Рисунок а показывает дина-

сообщим дополнительную начальную скорость v₀ >

мику для линейной цепи (наноленты), рис. б — для цик-

> 0, направленную вдоль оси x. Для этого численно

лической цепи (нанотрубки) с Nc = 20 звеньями. Штри-

проинтегрируем систему уравнений движения

ховые (красные) линии показывают зависимость (19) при

значении коэффициента эффективного трения γ = 8.3 и

∂H

M₂ü_n,k = -

8.6 нс-1

c ∂u_n,k

(15)

k = 1,...,K - 1, n = 1,...,N_bn,

подложки. Характер движения графеновых наноча-

∂H

стиц (нанолент и нанотрубок) показан на рис. 6.

M₁ ü_n,K = -

∂u_n,K

(16)

Как видно на рисунке, вид траектории x_c(t) за-

n = 1,...,N_c,

висит от реализации начального термализованного

состояния системы, но взаимодействие с подлож-

с начальными условиями

кой всегда приводит к торможению направленно-

u_n,k(0) = w_n,k,

u_n,k(0) = v_n,k,

го движения частицы (начальная скорость движе-

(17)

n = 1,...,N_bn, k = 1,...,K - 1,

ния v₀ = 500 м/с). Если траектории усреднить по

всем независимым реализациям термализованного

u_n,K(0) = w_n,K,

u_n,K(0) = v_n,K + v₀e_x,

состояния системы (взять среднее значение x_c(t) =

(18)

= 〈x_c(t)〉), то динамика наночастицы с хорошей точ-

n = 1,...,N_c,

ностью описывается как движение частицы в вязкой

где вектор e_x = (1, 0).

среде:

Проследим за движением центра тяжести нано-

частицы x_c = (x1,K + x2,K + . . . + xNc,K )/Nc вдоль

x_c(t) = x_c(0) + v₀ [1 - exp(-γt)]/γ,

(19)

890

ЖЭТФ, том 160, вып. 6 (12), 2021

Эффективное трение и подвижность графеновых наночастиц. ..

Заметим, что полученный существенный рост

трения с увеличением температуры парадоксаль-

ным образом отличается от стандартных сценариев

трения, где трение всегда уменьшается с увеличе-

нием температуры [29]. Отличие связано с тем, что

в рассматриваемой системе трение имеет волновую

природу (обусловлено взаимодействием с тепловы-

ми колебаниями). Увеличение силы трения слоев

G/h-BN с увеличением температуры также получе-

но в работе [30].

С увеличением размера наночастицы коэффици-

ент эффективного трения уменьшается. Так, при

T = 300 K для циклической G-цепи (для нанотруб-

ки) при числе звеньев N_c = 20 коэффициент трения

для ее движения по подложке γ = 8.6 нс-1, при N_c =

Рис. 7. Зависимость коэффициента трения γ от темпе-

= 40 коэффициент γ = 4.6 нс-1, а при N_c = 80 уже

ратуры T для линейной (наноленты) и циклической цепи

γ = 1.9 нс-1.

(нанотрубки) из Nc = 20 звеньев (соответственно кривые

Торможение наночастицы происходит из-за ее

1 и 2)

взаимодействия с атомами многослойной подлож-

ки. Для того чтобы определить распределение сил

торможения вдоль движущейся наноленты графе-

c коэффициентом вязкого «трения» γ > 0 (обрат-

на, найдем средние значения продольных сил F_x,n,

ная величина γ-1 соответствует времени, за кото-

действующих со стороны подложки на узлы G-цепи.

рое начальная скорость наночастицы уменьшается

Вектор силы

в e раз). Так, при N_c = 20 и температуре T = 300 K

∑

коэффициент эффективного трения для движения

∂W₂(rl,k;n,K)

F_n = (F_x,n, F_z,n) =

по подложке наноленты и нанотрубки γ = 8.3 и

∂u_n,K

k=1 l=1

8.6 нс-1, см. рис. 6.

Численное моделирование динамики наночастиц

расстояние rl,k;n,K = |u_n,K - u_l,k|.

показало, что ее вязкое торможение (19) происхо-

Для нахождения средних значений сил при чис-

дит при всех размерах наночастиц и температурах

ленном интегрировании системы уравнений движе-

подложки 0 < T ≤ 600K, но значение коэффициен-

ния (15), (16) с начальными условиями (17), (18)

та эффективного трения γ зависит от температуры,

найдем распределение вдоль G-цепи продольных

размера и типа наночастицы. Зависимость γ от тем-

сил взаимодействия с подложкой:

пературы показана на рис. 7. Как видно на рисунке,

∫

коэффициент трения монотонно растет с увеличени-

F_x,n(v) =

F_x,n(t)dt

ем температуры. При T ≤ 300K уменьшение темпе-

t₁

ратуры приводит к практически линейному умень-

шению эффективного трения γ до значений, близ-

при начальной скорости движения v = 500 м/с.

ких к нулю. Это позволяет заключить, что причи-

Возьмем время усреднения t₁ = 100 пс (как вид-

ной торможения наночастиц является их взаимо-

но на рис. 6, на начальном интервале времени [0, t₁]

действие с тепловыми колебаниями. Главную роль

наночастица всегда сохраняет направленное движе-

в этом взаимодействии играют изгибные колебания

ние). Найдем также средние значения этих сил для

цепей — если уменьшить их амплитуду, то значи-

неподвижной частицы

F_x,n(0) (значения сил при на-

тельно уменьшается и эффективное трение. Так, ес-

чальной скорости движения v = 0). Тогда можно

ли увеличить в десять раз энергию деформации ва-

определить силы торможения, препятствующие на-

лентных углов ϵ_i, i = 1, 2, увеличивая изгибную

правленному движению как

жесткость цепей, то коэффициент трения для цепи

из N_c = 20 звеньев при T = 300K уменьшится от

Fn,b(v) =

F_x,n(v) -

F_x,n(0).

значения γ = 8.3 до значения 3.5 нс-1 для линейной

и от значения γ = 8.6 до 3.2 нс-1 для циклической

Распределение сил торможения вдоль G-цепи

цепи.

(вдоль наноленты графена) показано на рис. 8. Как

891

А. В. Савин

ЖЭТФ, том 160, вып. 6 (12), 2021

Рис. 9. Зависимость коэффициента трения γ от числа

Рис. 8. Распределение вдоль G-цепи сил торможения дви-

звеньев линейной G-цепи Nc (от длины наноленты) при

жения (сил трения) F_b при числе звеньев цепи Nc = 40, 80,

температуре T = 100, 300, 500 K (кривые 1, 2, 3). Точ-

160 (кривые 1, 2, 3). Начальная скорость движения цепи

ки — значения, найденные численно, линии — зависимо-

v = 500 м/с, температура подложки T = 300 K

сти γ = γ1 + 2γ2/Nc со значениями коэффициентов γ1, γ2

соответственно 0.14, 22.6; 0.54, 77.5; 0.71, 152 нс-1

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

подложкой оказывают на края наноленты — на кра-

ях силы торможения значительно выше сил тормо-

внутреннего трения γ₁ = 0.14, коэффициент краево-

жения во внутренней части цепи. Моделирование

го трения γ₂ = 22.6 нс-1, при T = 300 K коэффи-

показало, что силы торможения прямо пропорцио-

циенты γ₁ = 0.54, γ₂ = 77.5 нс-1, при T = 500 K

нальны значению начальной скорости наночастицы

коэффициенты γ₁ = 0.71, γ₂ = 152 нс-1. Как сле-

v. Все это позволяет говорить о наличии двух сил

дует из этих цифр, коэффициент краевого всегда

трения — краевого трения, действующего на конце-

значительно превышает коэффициент внутреннего

вые атомы и значительно более слабого внутреннего

трения (γ₂ ≫ γ₁). Оба вида трения будут давать

трения, действующего на все остальные атомы цепи.

одинаковый вклад в торможение движения цепи при

Наличие краевого трения позволяет объяснить

числе звеньев

N_c = 2γ₂/γ₁. При меньшем числе зве-

зависимость коэффициента трения γ от длины цепи.

ньев главный вклад в торможение будет давать кра-

Зависимость γ от числа звеньев линейной цепи N_c

евое, а при большем — внутреннее трение. При T =

представлена на рис. 9. Как видно на рисунке, с уве-

= 100 K это характерное число узлов

N_c = 323,

N_c

личением длины цепи коэффициент эффективного

при T = 300 K число

= 287, при T = 500 K

трения монотонно уменьшается при всех значениях

N_c = 428. Таким образом, при размерах наночастиц

температуры.

L < 35 нм основной вклад в торможение ее движе-

Предположим, что при движении цепи по под-

ния по многослойной подложке будет давать кра-

ложке со скоростью v на каждый ее узел действует

евое трение (трение краев наночастицы с подлож-

одинаковая сила трения F₁ = -M₁γ₁v (сила тре-

кой), а при L > 35 нм — внутреннее трение (трение

ния соприкасающихся поверхностей) и сила F₂ =

внутренней поверхности наночастицы с подложкой).

= -M₁γ₂v, действующая только на краевые узлы

Отметим, что полученные оценки не зависят от

(сила краевого трения), тогда общий коэффициент

числа краевых узлов, на которые действует краевое

трения для цепи из N_c звеньев

трение. Если краевое трение действует на N_e конеч-

ных узлов с каждого края, то формула (20) примет

γ(N_c) = γ₁ + 2γ₂/N_c.

(20)

вид

Полученные из численного моделирования динами-

ки значения коэффициента трения γ для цепей с

γ(N_c) = γ₁ + 2N_eγ′2/N_c,

(21)

числом звеньев N_c = 20, 40, 80, 160, 320 с хоро-

шей точностью подтверждают зависимость (20) (см.

рис. 9). При температуре T = 100 K коэффициент

где N_eγ′2 = γ₂.

892

ЖЭТФ, том 160, вып. 6 (12), 2021

Эффективное трение и подвижность графеновых наночастиц. ..

4. ПОДВИЖНОСТЬ НАНОЧАСТИЦ

Если частица движется в вязкой среде под дей-

ствием внешней силы F_x, то ее динамика описыва-

ется уравнением движения

Mx_c = -γM ˙x_c + F_x,

(22)

где M

= N_cM₁ — масса частицы. Из уравнения

(22) следует, что со временем частица всегда будет

двигаться с установившейся постоянной скоростью

v_F = F_x/γM. Здесь подвижность частицы, отноше-

ние средней скорости ее движения к силе,

μ = v_F/F_x = 1/γN_cM₁.

(23)

Найдем подвижность наночастицы μ из прямо-

го моделирования ее движения по термализованной

подложке под действием постоянной силы F_x > 0 на

узел G-цепи с номером n = N_c/2. Здесь динамика

узлов G-цепи будет задаваться системой уравнений

движения

∂H

M₁ ü_n,K = -

+ F_xδn,Nc/2ex,

Рис. 10. Зависимость центра тяжести графеновой наночас-

∂u_n,K

(24)

тицы xc из Nc = 20 звеньев от времени t при температуре

n = 1,...,N_c.

подложки T = 300 K и внешней силе Fx = 0.001 эВ/Å.

Тонкие кривые показывают 9 траекторий движения для

Таким образом, для моделирования вынужденного

различных независимых начальных реализаций термали-

движения наночастицы нужно численно проинте-

зованного состояния системы. Сплошная толстая (черная)

грировать систему уравнений движения (15), (24) с

линия дает траекторию среднего значения центра тяже-

начальными условиями (17), (18), где скорость на-

сти xc, полученную по 256 независимым реализациям тер-

чального движения v₀ = F_x/γN_cM₁.

мализованного состояния системы. Рисунок а показывает

Численное моделирование движения показало,

динамику для линейной цепи (наноленты), рис. б — для

что под действием постоянной силы F_x > 0 нано-

циклической цепи (нанотрубки). Скорости установившего-

ся движения частиц vF = 542 и 530 м/с

частица всегда (при любом значении температуры

T и числа звеньев N_c) выходит на режим движе-

ния вдоль подложки с постоянной скоростью v_F > 0

Значение постоянной скорости движения нано-

(см. рис. 10). Конечно, вид траектории движения

частицы v_F практически линейно зависит от вели-

центра тяжести наночастицы x_c(t) зависит от вы-

чины силы F_x. Поэтому мы всегда можем найти по-

бора начального термализованного состояния систе-

движность наночастицы μ = v_F /F_x. Зависимость

мы {w_n,k, v_n,k}Nj,Kn=1,k=1 — могут происходить неболь-

подвижности наночастицы от температуры подлож-

шие флуктуации скорости движения частицы. Но

ки показана на рис. 11. Как видно на рисунке, по-

если усреднить по всем независимым реализациям

движность монотонно убывает с ростом температу-

начального термализованного состояния, то движе-

ры, а значение подвижности хорошо приближается

ние центра тяжести будет происходить с постоянной

формулой (23).

скоростью

Таким образом, полученные значения подвиж-

ности подтверждают, что движение наночастиц на

x_c(t) = 〈x_c(t)〉 = x_c(0) + v_F t.

(25)

многослойной термализованной подложке носит ха-

Таким образом, движение наночастицы под дейст-

рактер движения частицы в вязкой среде, вязкость

вием постоянной силы по многослойной подложке

которой сильно зависит от значения температуры.

происходит с постоянным значением скорости, а при

С увеличением температуры возрастает сопротив-

выходе частицы на край подложки происходит ее от-

ление среды движению наночастицы. Возникающее

ражение от края (см. рис. 10).

эффективное трение имеет волновую природу, его

893

А. В. Савин

ЖЭТФ, том 160, вып. 6 (12), 2021

Рис. 11. Зависимость подвижности μ линейной и цикли-

Рис. 12. Распределение сил торможения движения (сил

ческой G-цепи (наноленты и нанотрубки) из Nc = 20 зве-

трения) F_b вдоль G-цепи из Nc = 80 звеньев при силе дав-

ньев от температуры подложки T. Сплошные кривые 1

ления Fz = 0 (кривая 1) и Fz = 0.03 эВ/Å. Начальная

и 2 дают значения, полученные из прямого моделирова-

скорость движения цепи v = 500 м/с, температура под-

ния движения наночастицы под действием внешней силы,

ложки T = 300 K

штриховые кривые 3 и 4 — значения, полученные из фор-

мулы (23). Кривые 1 и 3 дают зависимости для наноленты,

2 и 4 — для нанотрубки

мающей ее к подложке вертикальной силы нужно

численно проинтегрировать систему уравнений дви-

жения (15), (26) с начальными условиями (17), (18).

значение обусловлено в первую очередь взаимодей-

Численное моделирование динамики показало,

ствием наночастицы с тепловыми поперечными ко-

что добавление вертикальной силы, прижимающей

лебаниями слоев подложки. Чем сильнее эти коле-

G-цепь к подложке, приводит к возрастанию крае-

бания (чем выше температура), тем сильнее это вза-

вых и уменьшению внутренних сил торможения (см.

имодействие и тем сильнее происходит отток кине-

рис. 12). Так, при начальной скорости движения v =

тической энергии наночастицы в подложку.

= 500 м/с, температуре T = 300 K для цепи из N_c =

= 80 звеньев в отсутствие давления (при силе F_z =

= 0) сила внутреннего торможения

5. ВЛИЯНИЕ НА ТРЕНИЕ СИЛЫ

ДАВЛЕНИЯ НАНОЧАСТИЦЫ НА

∑

ПОДЛОЖКУ

Fb,1 =

F_n,b = 0.259 · 10-3 эВ/Å,

n=11

Согласно эмпирическому закону Амонтона - Ку-

лона сила трения скольжения одного тела по по-

сила краевого торможения

верхности другого тела (опоры) пропорциональна

∑

силе давления тела на опору (чем больше давле-

Fb,2 =

F_n,b+

F_n,b = 0.848 · 10-3 эВ/Å.

ние, тем больше сила трения). Проверим этот за-

n=1

n=Nc-9

кон для наночастиц, находящихся на многослойной

плоской подложке, выполняющей роль опоры. Для

При давлении F_z

= 0.03 эВ/Å сила внутренне-

этого рассмотрим движение линейной G-цепи (на-

го торможения уменьшается Fb,1 = 0.249 эВ/Å, а

ноленты) в присутствии силы F_z, вертикально при-

сила краевого торможения увеличивается Fb,2

жимающей ее к подложке. В этом случае динамика

= 1.068 эВ/Å.

звеньев G-цепи будет описываться системой уравне-

С увеличением силы давления F_z коэффициент

ний движения:

эффективного трения наночастицы с подложкой γ

монотонно увеличивается (см. рис. 13). Полученные

∂H

M₁ ü_n,K = -

- F_ze_z, n = 1,...,N_c,

(26)

зависимости коэффициента трения от числа звеньев

∂u_n,K

цепи N_c (длины наноленты) при разных значениях

где вектор e_z = (0, 1). Таким образом, для моделиро-

силы давления с хорошей точностью аппроксимиру-

вания движения наночастицы при наличии прижи-

ются формулой (20). Поэтому при каждом значении

894

ЖЭТФ, том 160, вып. 6 (12), 2021

Эффективное трение и подвижность графеновых наночастиц. ..

γ₁ монотонно уменьшается с увеличением силы дав-

ления F_z .

Полученные зависимости γ₁(F_z), γ₂(F_z) позволя-

ют заключить, что при размерах наночастицы, при

которых главную роль в торможении ее движения

по подложке играет краевое трение, закон Амонто-

на - Кулона выполняется, а при размерах, при кото-

рых главную роль играет внутреннее трение, — не

выполняется (здесь увеличение нагрузки будет уже

приводить к уменьшению общего трения). Таким об-

разом, для нанолент графена на многослойной под-

ложке h-BN эмпирический закон Амонтана-Кулона

будет выполняться при длинах L < 250 нм и пере-

станет выполняться при L > 250 нм.

Заметим, что понижение трения между листами

графена и h-BN при увеличении давления, т. е. пря-

Рис. 13. Зависимость коэффициента трения γ от числа

мое нарушение закона Амонтана - Кулона, было об-

звеньев линейной G-цепи Nc (от дины наноленты) при си-

наружено в работе [30] при моделировании протас-

ле давления Fz = -0.01, 0, 0.01, 0.03, 0.05 эВ/Å (кри-

кивания листа графена по листу нитрида бора (при

вые 1, 2, 3, 4, 5). Точки — значения, найденные числен-

увеличении давления сила трения может уменьшит-

но, линии — аппроксимирующие эти значения зависимости

ся на 30 %). При протаскивании бесконечных листов

γ = γ1 + 2γ2/Nc. Температура подложки T = 300 K

краевое трение отсутствует, поэтому этот резуль-

тат согласуется с понижением коэффициента внут-

реннего трения γ₁ при повышении внешней нагруз-

ки. Здесь имеет место эффект «отрицательного ко-

эффициента трения (ОКТ)» — уменьшение трения

при увеличении нормальной нагрузки [36, 37]. Это

парадоксальное поведение связано с тем, что нор-

мальная нагрузка приводит к уменьшению ампли-

туды изгибных тепловых колебаний слоев (в силу

этого уменьшается коэффициент внутреннего тре-

ния γ₁). Но нормальная нагрузка на конечный лист

графена приводит к вдавливанию краев листа в под-

ложку. Из-за этого краевое сопротивление (коэффи-

циент краевого трения γ₂) увеличивается пропор-

ционально увеличению нагрузки. Поэтому эффект

OKT будет наблюдаться только для листов доста-

точно больших размеров, когда трение краев листа

будет оказывать более слабое влияние, чем трение с

Рис. 14. Зависимость коэффициентов краевого γ2 и внут-

подложкой внутренней поверхности листа.

реннего γ1 трения (рис. а,б) от силы давления на нанолен-

ту Fz при температуре подложки T = 300 K

6. ЗАКЛЮЧЕНИЕ

силы F_z мы можем определить значение коэффици-

Моделирование свободного движения графено-

ентов внутреннего γ₁ и краевого γ₂ трения. Зависи-

вых наночастиц (нанолент и нанотрубок) показало,

мости этих коэффициентов от силы давления пред-

что их движение на плоской многослойной подлож-

ставлены на рис. 14. Как видно на рисунке, вели-

ке h-BN можно описать как движение частицы в

чина коэффициента краевого трения (трения краев

вязкой среде с постоянным коэффициентом трения

наночастицы с подложкой) γ₂ линейно растет, а ве-

(коэффициентом релаксации скорости) γ. Торможе-

личина коэффициента внутреннего трения (трения

ние движения происходит из-за перекачки кинети-

внутренней поверхности наночастицы с подложкой)

ческой энергии наночаститцы в подложку. Величи-

895

А. В. Савин

ЖЭТФ, том 160, вып. 6 (12), 2021

на коэффициента эффективного трения наночасти-

та трения» (уменьшение трения при увеличении

цы с подложкой монотонно увеличивается с увели-

нормальной нагрузки) может наблюдаться только

чением температуры T . При T ≤ 300K уменьшение

для нанолент достаточно большого размера (L >

температуры приводит к линейному уменьшению γ

> 250 нм), когда главную роль играет внутреннее

до очень низких значений. Такое поведение объяс-

трение.

няется волновой природой эффективного трения —

причиной торможения движения наночастиц явля-

Финансирование. Научно-исследовательская

ется их взаимодействие с тепловыми (в первую оче-

работа выполнена за счет субсидии, выделен-

редь изгибными) колебаниями листов подложки.

ной Федеральным исследовательским центром

С увеличением размера наночастицы коэффици-

химической физики Российской академии наук

ент эффективного трения монотонно уменьшается

на выполнение государственного задания (тема

при всех значениях температуры. Анализ распреде-

0082-2019-0005). Вычислительные ресурсы предо-

ления вдоль движущейся наноленты сил торможе-

ставлены межведомственным суперкомпьютерным

ния показывает, что трение разделяется на два ти-

центром Российской академии наук.

па: на внутреннее и краевое (на трение внутренней

поверхности наноленты и ее краев с поверхностью

ЛИТЕРАТУРА

подложки). Такое разделение сил трения объясняет

K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,

зависимость значения коэффициента от длины на-

Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and

ноленты. Коэффициент краевого трения γ₂ всегда

A. A. Firsov, Science 306(5696), 666 (2004).

значительно превышает коэффициент внутреннего

трения γ₁. Общий коэффициент γ = γ₁ +2γ₂a/L, где

A. H. Castro Neto, F. Guinea, N. M. R. Peres,

K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys.

a — продольный шаг наноленты, L — ее длина. При

81, 109 (2009).

длинах L < 35 нм главную роль в торможении дви-

жения наноленты играет краевое, а при L > 35 нм —

E. Koren, I. Leven, E. Lörtscher, A. Knoll, O. Hod,

внутреннее трение.

and U. Duerig, Nature Nanotech. 11, 752 (2016).

Моделирование подвижности наночастиц пока-

J. C. Meyer, A. K. Geim, M. Katsnelson, K. Novose-

зало, что под действием постоянной продольной си-

lov, T. Booth, and S. Roth, Nature 446, 60 (2007).

лы их динамика всегда выходит на режим движе-

ния с постоянной скоростью, значение которой пря-

C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science

мо пропорционально значению силы и обратно про-

321, 385 (2008).

порционально значению коэффициента трения γ.

A. Falin, Q. Cai, E. J. G. Santos, D. Scullion, D. Qian,

Это подтверждает, что динамика графеновых на-

R. Zhang, Z. Yang, S. Huang, K. Watanabe, T. Ta-

ночастиц (нанолент и нанотрубок) на многослойной

niguchi, M. R. Barnett, Y. Chen, R. S. Ruoff, and

термализованной подложке выглядит как движения

L. H. Li, Nat. Commun. 8, 15815 (2017).

частицы в вязкой среде, вязкость которой увеличи-

E. Han, J. Yu, E. Annevelink, J. Son, D. A. Kang,

вается с ростом температуры. Возникающее эффек-

K. Watanabe, T. Taniguchi, E. Ertekin, P. Y. Huang,

тивное трение имеет волновую природу, его значе-

and A. M. van der Zande, Nature Materials 19, 305

ние обусловлено в первую очередь взаимодействи-

(2020).

ем наночастицы с тепловыми поперечными колеба-

P. E. Sheehan and C. M. Lieber, Science 272, 1158

ниями слоев подложки. Чем сильнее эти колебания

(1996).

(чем выше температура), тем сильнее это взаимо-

действие и тем сильнее происходит торможение дви-

M. Dienwiebel, G. S. Verhoeven, N. Pradeep,

жения частицы.

J. W. Frenken, J. A. Heimberg, and H. W. Zand-

Моделирование движения наноленты при на-

bergen, Phys. Rev. Lett. 92, 126101 (2004).

личии нормальной (прижимающей к подложке)

10.

C. Lee, Q. Li, W. Kalb, X.-Z. Liu, H. Berger,

нагрузки показало, что рост нагрузки приводит к

R. W. Carpick, and J. Hone, Science 328, 76 (2010).

уменьшению внутреннего трения из-за уменьшения

11.

S. Cahangirov, C. Ataca, M. Topsakal, H. Sahin, and

под нанолентой амплитуды тепловых изгибных

S. Ciraci, Phys. Rev. Lett. 108, 126103 (2012).

колебаний слоев подложки и к увеличению крае-

вого трения из-за вдавливания краев наноленты в

12.

Z. Liu, J. Yang, F. Grey, J. Z. Liu, Y. Liu, Y. Wang,

подложку. Поэтому замеченный ранее в работах

Y. Yang, Y. Cheng, and Q. Zheng, Phys. Rev. Lett.

[30, 36, 37] эффект «отрицательного коэффициен-

108, 205503 (2012).

896

ЖЭТФ, том 160, вып. 6 (12), 2021

Эффективное трение и подвижность графеновых наночастиц. ..

13.

J. Yang, Z. Liu, F. Grey, Z. Xu, X. Li, Y. Liu, M. Ur-

25.

M. Z. Baykara, M. R. Vazirisereshk, and A. Martini,

bakh, Y. Cheng, and Q. Zheng, Phys. Rev. Lett. 110,

Appl. Phys. Rev. 5, 041102 (2018).

255504 (2013).

26.

D. Berman, A. Erdemir, and A. V. Sumant, ACS

Nano 12, 2122 (2018).

14.

E. Koren, E. Lörtscher, C. Rawlings, A. W. Knoll,

and U. Duerig, Science 348, 679 (2015).

27.

Y. Song, D. Mandelli, O. Hod, M. Urbakh, M. Ma,

and Q. Zheng, Nat. Mater. 17, 894 (2018).

15.

I. Leven, D. Krepel, O. Shemesh, and O. Hod,

J. Phys. Chem. Lett. 4, 115 (2013).

28.

W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod,

Nano Lett. 18(9), 6009 (2018).

16.

A. Geim and I. Grigorieva, Nature 499, 419 (2013).

29.

A. Vanossi, N. Manini, M. Urbakh, S. Zapperi, and

17.

K. S. Novoselov, A. Mishchenko, A. Carvalho, and

E. Tosatti, Rev. Mod. Phys. 85, 529 (2013).

A. H. Castro Neto, Science 353(6298), 461 (2016).

30.

D. Mandelli, W. Ouyang, O. Hod, and M. Urbakh,

18.

C. R. Woods, L. Britnell, A. Eckmann, R. S. Ma,

Phys. Rev. Lett. 122, 076102 (2019).

J. C. Lu, H. M. Guo, X. Lin, G. L. Yu, Y. Cao,

31.

A. V. Savin, E. A. Korznikova, and S. V. Dmitriev,

R. V. Gorbachev, A. V. Kretinin, J. Park, L. A. Po-

Phys. Rev. B 92, 035412 (2015).

nomarenko, M. I. Katsnelson, Y. N. Gornostyrev,

K. Watanabe, T. Taniguchi, C. Casiraghi, H. J. Gao,

32.

А. В. Савин, Е. А. Корзникова, С. В. Дмитриев,

A. K. Geim, and K. S. Novoselov, Nature Phys. 10,

ФТТ 57(11), 2278 (2015).

451 (2014).

33.

A. V. Savin, E. A. Korznikova, and S. V. Dmitriev,

19.

G. J. Slotman, M. M. van Wijk, P. L. Zhao, A. Faso-

Phys. Rev. B 99, 235411 (2019).

lino, M. I. Katsnelson, and S. Yuan, Phys. Rev. Lett.

115, 186801 (2015).

34.

J. H. Los, J. M. H. Kroes, K. Albe, R. M. Gordillo,

M. I. Katsnelson, and A. Fasolino, Phys. Rev. B 96,

20.

D. Mandelli, I. Leven, O. Hod, and M. Urbakh, Sci.

184108 (2017).

Rep. 7(1), 10851 (2017).

35.

A. K. Rappé, C. J. Casewit, K. S. Colwell, W. A. God-

21.

J. A. Williams and H. R. Le, J. Phys. D: Appl. Phys.

dard III, and W. M. Skiff, J. Amer. Chem, Soc. 114,

39, R201 (2006).

10024 (1992).

22.

K. Shinjo and M. Hirano, Surface Science 283, 473

36.

S.-W. Liu, H.-P. Wang, Q. Xu, T.-B. Ma, G. Yu,

(1993).

C. Zhang, D. Geng, Z. Yu, S. Zhang, W. Wang,

Y.-Z. Hu, H. Wang, and J. Luo, Nat. Commun. 8,

23.

O. Hod, E. Meyer, Q. Zheng, and M. Urbakh, Nature

14029 (2017).

(London) 563, 485 (2018).

37.

Y. Liu, A. Song, Z. Xu, R. Zong, J. Zhang, W. Yang,

24.

J. M. Martin and A. Erdemir, Phys. Today 71, No. 4,

R. Wang, Y. Hu, J. Luo, and T. Ma, ACS Nano 12,

40 (2018).

7638 (2018).

897