ЖЭТФ, 2022, том 162, вып. 6 (12), стр. 991-1003
© 2022
ИССЛЕДОВАНИЕ СПЕКТРАЛЬНЫХ СВОЙСТВ
ПРОСТРАНСТВЕННО-НЕОДНОРОДНОЙ СИСТЕМЫ ЧАСТИЦ
ЮКАВЫ В ПАРАБОЛИЧЕСКОМ КОНФАЙНМЕНТЕ
И. В. Вороновa,b*, В. С. Николаевa,b, А. В. Тимофеевb,a,c
a Московский физико-технический институт (Национальный исследовательский университет)
141700, Долгопрудный, Московская обл., Россия
b Объединенный институт высоких температур Российской академии наук
125412, Москва, Россия
c Национальный исследовательский университет ¾Высшая школа экономики¿
101000, Москва, Россия
Поступила в редакцию 17 июня 2022 г.,
после переработки 27 июля 2022 г.
Принята к публикации 2 августа 2022 г.
Исследуются динамические свойства системы из конечного числа одноименно заряженных частиц, взаи-
модействующих по экранированному кулоновскому потенциалу в поле параболической ловушки. Методом
анализа нормальных мод в квазигармоническом приближении получено пространственное распределе-
ние амплитуд колебаний частиц в системе: высокочастотные колебания локализуются в ее центральной
области. Это приводит к тому, что в разных областях структуры спектры колебаний частиц различны.
Для обобщения этого результата на широкий температурный диапазон проведено численное модели-
рование систем из конечного числа заряженных частиц в ловушке методом классической молекулярной
динамики. Показана пространственная неоднородность спектров тепловых колебаний частиц: происходит
смещение спектра в низкочастотную область при увеличении расстояния от рассматриваемой подсистемы
до центра структуры. Различные методы расчета амплитуды тепловых колебаний частиц как в рамках
квазигармонического приближения, так и методом молекулярной динамики подтверждают неоднород-
ность радиального профиля амплитуды: на периферии системы амплитуда тепловых колебаний частиц
имеет более высокое значение, чем в центре структуры. Показано, что в локальном приближении динами-
ческие свойства подсистем структур из конечного числа заряженных частиц в ловушке при определенных
условиях существенно отличаются от динамических свойств систем заряженных частиц близкой плот-
ности в периодических граничных условиях. Результаты настоящего исследования могут быть полезны
при описании динамических свойств лабораторных плазменно-пылевых и коллоидных систем.
DOI: 10.31857/S0044451022120197
тенциал, совпадающий с отталкивательным потен-
EDN: LFJRGM
циалом Юкавы по функциональной форме:
Q2
UCoul(rij) =
exp(-κrij),
(1)
r
ij
1. ВВЕДЕНИЕ
где Q электрический заряд частицы, rij = |ri -rj|
Система частиц, взаимодействующих по потен-
расстояние между частицами i и j, κ параметр
циалу Юкавы, используется в качестве модельной
экранирования. В силу своего отталкивательного
как для систем элементарных частиц [1], так и для
характера потенциал такой формы может быть ис-
коллоидной [2] и пылевой [3-18] плазмы. При описа-
пользован только для описания модельных ¾беско-
нии систем одноименно заряженных частиц в плаз-
нечных¿ однородных систем заряженных частиц в
ме используется экранированный кулоновский по-
периодических граничных условиях. Структурные,
динамические и транспортные свойства таких си-
* E-mail: ilya.v.voronov@gmail.com
стем детально изучены в работах [19-21]. В идеаль-
991
И. В. Воронов, В. С. Николаев, А. В. Тимофеев
ЖЭТФ, том 162, вып. 6 (12), 2022
ных условиях описанные структуры образуют кри-
рот, увеличивается от центра системы к ее краю.
сталлическую решетку из одноименно заряженных
В работах [26, 27] подчеркивается, что структурная
частиц: в трехмерном случае ОЦК или ГЦК, в дву-
неоднородность системы заряженных частиц связа-
мерном гексагональную.
на как с наличием удерживающей ловушки, так и с
В лабораторных экспериментах исследуются
экранировкой (κ = 0) потенциала взаимодействия.
структуры из конечного числа частиц с открытой
В работе [32] предлагается локальный подход к
границей. Примерами могут служить коллоидные
описанию динамических свойств пространственно-
системы, пылевая и однокомпонентная плазма.
неоднородных систем заряженных частиц в удер-
В таких системах одним из подходов для удержа-
живающей ловушке. Показано, что амплитуда теп-
ния отталкивающихся частиц от разлета является
ловых колебаний частиц увеличивается от центра
воздействие внешней электростатической ловуш-
структуры к ее краю. Также в работе [32] обсужда-
ки. Например, в плазменно-пылевых системах
ется неоднородность пространственного распределе-
конфайнмент может быть связан со следующими
ния температуры плавления: отдельные области си-
факторами: в тлеющем разряде ловушка вызвана
стемы плавятся при различных температурах. Де-
амбиполярной диффузией [22], а в высокочастот-
лается вывод, что в таких системах может наблю-
ном разряде она создается намеренно с помощью
даться сосуществование двух фаз: твердой в цен-
специального кольца [23]. В первом приближении
тре системы и жидкой на периферии. Это проис-
конфайнмент в таких системах имеет параболиче-
ходит из-за того, что от центра к периферии ам-
ский профиль [24]:
плитуда тепловых колебаний частиц растет быст-
Utrap(xi, yi, zi) = 0.5Q(αx2i + βy2i + γz2i),
(2)
рее, чем среднее межчастичное расстояние. Отно-
шение этих двух параметров определяет величину
где α, β, γ параметры ловушки. В зависимости от
параметра Линдеманна, что в первом приближении
их соотношения получаются трех-, двух- и одномер-
и указывает на неоднородный характер процесса
ные системы. Обоснованность использования функ-
плавления. Позднее в работе [33] показано, что влия-
циональной формы потенциалов (1), (2) подтвер-
ние неустойчивости связанных мод может привести
ждается теоретическими результатами работ [12,25]
к диаметрально противоположному эффекту: сна-
по исследованию пространственного распределения
чала плавится центральная часть квазидвумерной
плазмы в положительном столбе тлеющего разряда.
структуры, а затем периферия.
Результаты этих работ указывают на то, что в об-
ласти существования пылевой структуры профиль
Наряду с амплитудными свойствами исследуют-
плотности ионов и электронов плазмы в радиаль-
ся и спектральные характеристики систем из конеч-
ном направлении сглаживается, градиент плотно-
ного числа заряженных частиц в ловушке [34, 35].
сти ионов и электронов вблизи оси разрядной труб-
В работах [34, 35] для исследования спектральных
ки уменьшается. Радиальная компонента электри-
свойств систем применяется метод анализа нор-
ческого поля в этих условиях имеет линейную за-
мальных мод. Он позволяет получить информацию
висимость от радиального расстояния, что соответ-
о колебательных модах и сравнить колебательные
ствуют ловушке с параболическим профилем вбли-
процессы между системой из конечного числа час-
зи центра системы. Таким образом, модель, включа-
тиц и однородной системой в периодических гранич-
ющая в себя потенциалы (1) и (2), может использо-
ных условиях. В работе [34] проведен анализ спек-
ваться в качестве первого приближения для описа-
тров колебаний частиц в системе из конечного чис-
ния подсистемы пылевых частиц в плазме газового
ла частиц и приведено замечание, что существует
разряда.
некоторая максимальная частота колебаний частиц,
В данной работе в центре внимания находятся
ограничивающая спектр колебаний. Это ограниче-
двумерные системы из конечного числа заряжен-
ние частоты колебаний частиц в системе связано с
ных частиц в параболической ловушке, для которых
числом частиц и, соответственно, с размером систе-
α = β ≪ γ. Такие системы хорошо изучены с точ-
мы. В работе [35] показано, что колебания частиц
ки зрения структурных свойств [26-31]. Им присуща
в таких системах можно разделить на колебатель-
пространственная неоднородность, проявляющаяся
ные процессы, имеющие сходство с волнами сжатия
в монотонном убывании плотности частиц с увели-
и волнами сдвига в ¾бесконечных¿ юкавовских си-
чением радиального расстояния как n(r) = a - br2
стемах. Отметим, что в работах [34,35] анализ нор-
[30]. Следовательно, среднее межчастичное рассто-
мальных мод проведен для систем из малого числа
яние Δlocal(r) ∝ n-1/2(r) в таких системах, наобо-
заряженных частиц и не было уделено достаточно
992
ЖЭТФ, том 162, вып. 6 (12), 2022Исследование спектральных свойств пространственно-неоднородной системы. . .
внимания влиянию неоднородности пространствен-
неоднородных ¾конечных¿ и однородных ¾бесконеч-
ных свойств систем в ловушке на их динамические
ных¿ системах. Обнаружено различие в амплитудах
свойства. В настоящем же исследовании проводится
тепловых колебаний частиц в центральной квазиод-
анализ нормальных мод с целью установления вли-
нородной подсистеме конечной структуры и одно-
яния неоднородности пространственных свойств та-
родной системе в ПГУ при совпадении их плотно-
ких систем на их динамические свойства.
стей вплоть до числа частиц N = 105 в обеих систе-
мах.
Изучение влияния неоднородности простран-
ственных свойств на динамические свойства
могло бы быть полезным для описания некоторых
2. ОПИСАНИЕ МОДЕЛИ И ЛОКАЛЬНОГО
экспериментальных результатов. К примеру, в
ПОДХОДА К РАСЧЕТУ ДИНАМИЧЕСКИХ
экспериментальных работах [36-38] из плазменно-
СВОЙСТВ КОНЕЧНОЙ СИСТЕМЫ
пылевой структуры выделяется центральная
В настоящей работе с помощью численного мето-
область, свойства которой далее описываются с
да, аналогичного методу классической молекуляр-
помощью численного моделирования в периодиче-
ной динамики, моделируется двумерная система од-
ских граничных условиях. Заметим, что, несмотря
ноименно заряженных частиц, взаимодействующих
на обилие работ, в которых используется данный
по экранированному кулоновскому потенциалу (1)
подход, вопрос обоснованности его применения все
и находящихся в поле параболической электроста-
еще нельзя считать закрытым [39-41].
тической ловушки (2). Моделирование проводится
Настоящая статья построена следующим обра-
с применением пакета LAMMPS. Полная потенци-
зом. В разд. 2 подробно описываются рассматрива-
альная энергия системы записывается следующим
емая модель, методика исследования и локальный
образом:
подход к изучению динамических свойств системы
∑∑
из конечного числа частиц в ловушке. В разд. 3 ди-
Q2
E = 0.5α (x2i + y2i ) +
exp(-κrij),
(3)
намические свойства конечной системы рассматри-
rij
i=1
i=1 i>j
ваются в низкотемпературном квазигармоническом
приближении с помощью метода анализа нормаль-
где α = 0.09 ед. СГСЭ, N = 500 частиц, Q = 5000e,
ных мод. Показано, что в зависимости от частоты
κ = 100 см-1. Численное интегрирование уравне-
колебательные процессы локализуются в различных
ний движения выполняется с применением алгорит-
ма Верле с шагом по времени τ = 10-5 с. Система
областях системы. В разд. 4 исследуются динамиче-
ские свойства системы из конечного числа заряжен-
подготавливается в мягком термостате Нозе-Гувера
с температурой T
= 300 K и параметром релак-
ных частиц в ловушке с помощью метода молеку-
лярной динамики, поскольку данный метод, в отли-
сации Tdamp = 1000τ в течение 5 · 106 шагов. По-
чие от метода анализа нормальных мод, эффекти-
сле релаксации системы происходит расчет интере-
вен при большом числе частиц в системе и применим
сующих характеристик в течение 106 шагов, мяг-
в широком диапазоне температур. Наблюдается сов-
кий термостат при этом остается включенным. Рас-
падение результатов расчета динамических свойств
считанные характеристики усредняются более чем
в разд. 3 и 4. Показана неоднородность динамиче-
по нескольким идентичным статистически независи-
мым структурам. Контроль сходимости рассчиты-
ских свойств и предложены методы их расчета. В
разд. 5 проводится сравнительный анализ динами-
ваемых характеристик проводится по зависимости
усредненной величины от времени усреднения.
ческих свойств ¾конечных¿ и ¾бесконечных¿ струк-
тур. В рассмотрении участвуют набор систем из ко-
Для реализации метода анализа нормальных
мод в разд. 3 используется динамическая матрица,
нечного числа заряженных частиц в ловушке и на-
бор однородных систем заряженных частиц в перио-
полученная численно для исследуемой методом мо-
дических граничных условиях (ПГУ). Набор систем
лекулярной динамики структуры с помощью паке-
в ловушке состоит из нескольких пространственно-
та LAMMPS. Выражение для расчета динамической
неоднородных структур с совпадающим межчастич-
матрицы имеет вид
ным расстоянием в центре и различным числом час-
2E
тиц. Набор систем в ПГУ состоит из нескольких од-
Eα,β,i,j =
,
(4)
∂rα,i∂rβ,j
нородных структур с совпадающей плотностью час-
тиц и также различным числом частиц. Проводит-
где α, β = {x, y}, т.е. rα,i обозначает координату x
ся сравнение спектров тепловых колебаний частиц в
или y частицы с номером i из моделирования, а диф-
993
И. В. Воронов, В. С. Николаев, А. В. Тимофеев
ЖЭТФ, том 162, вып. 6 (12), 2022
ференцируемая величина E полная потенциаль-
ная энергия системы (3). Затем с применением соб-
ственного программного кода происходит поиск соб-
ственных чисел и векторов динамической матрицы.
Полученные характеристики усредняются по стати-
стически независимым структурам.
В рамках численного моделирования в разд. 4
конечные системы заряженных частиц в ловушке
рассматриваются с применением локального подхо-
да. В структуре выделяются области, такие что в
них плотность частиц приблизительно постоянна,
а число частиц достаточно велико для проявления
подсистемой макроскопических свойств. Целесооб-
разность такого подхода заключается в том, что во
многих классах систем из конечного числа частиц
(например, пылевая плазма, коллоидные системы)
наблюдается наличие зависимости плотности час-
тиц от радиального расстояния. В силу цилиндри-
Рис. 1. В конечной двумерной системе мысленно выде-
ческой симметрии рассматриваемой в данной рабо-
ляются кольца, в каждое из которых входит подсистема
те системы, в качестве описанных выше областей в
с постоянной плотностью и числом частиц N ≫ 1 (A).
ней можно выделить кольца, в которых число час-
В разд. 5 каждому такому кольцу конечной структуры (A)
тиц велико, и плотность частиц близка к постоян-
ставится в соответствие однородная система в периодиче-
ной. Характерное разделение исследуемой структу-
ских граничных условиях с совпадающей плотностью (B)
ры на кольца можно видеть на рис. 1 A. В насто-
ящей работе сравниваются характеристики, усред-
IV
ненные по всем частицам в указанных подсистемах.
Также в настоящей работе в разд. 5 соотносятся
0.016
свойства выделенных цветом квазиоднородных под-
III
систем конечной системы со свойствами однородной
0.014
системы в периодических граничных условиях с сов-
падающей плотностью (см. рис. 1 B). Рассматрива-
II
емая в работе структура из 500 частиц имеет ра-
0.012
диус R0 ≈ 0.15 см. Разбиение системы на кольце-
I
вые фрагменты проводится с учетом неоднородно-
0.010
сти радиального профиля плотности частиц в систе-
0.00
0.05
0.10
0.15
ме [26,27]. Принцип выбора границ разбиения пока-
r, cm
зан на рис. 2. Во фрагментах находится приблизи-
тельно следующее число частиц: I
200, II
150,
III
100, IV
50.
Рис. 2. Радиальный профиль межчастичного расстояния
local(r). Штриховыми линиями обозначено разбиение си-
стемы на кольцевые фрагменты в локальном подходе
3. АНАЛИЗ ДИНАМИЧЕСКИХ СВОЙСТВ
СИСТЕМ ЗАРЯЖЕННЫХ ЧАСТИЦ
векторов колебаний частиц в различных модах. Все-
В ЛОВУШКЕ В НИЗКОТЕМПЕРАТУРНОМ
го в произвольной системе из N частиц существует
КВАЗИГАРМОНИЧЕСКОМ
DN нормальных мод, где D = 2 размерность про-
ПРИБЛИЖЕНИИ
странства. Из DN нормальных мод ftrans мод соот-
В настоящей работе применяется метод анали-
ветствуют поступательному движению центра масс
за нормальных мод для исследования динамических
системы, frot мод соответствуют вращению системы
свойств системы заряженных частиц в ловушке в
как целого. В рассматриваемом в работе двумерном
рамках низкотемпературного квазигармонического
случае ftrans = 2, frot = 1. Таким образом, в дву-
приближения. Метод анализа нормальных мод ис-
мерной системе из N частиц наблюдается 2N - 3
пользуется для нахождения собственных частот и
нормальных колебательных мод. Поступательные и
994
ЖЭТФ, том 162, вып. 6 (12), 2022Исследование спектральных свойств пространственно-неоднородной системы. . .
Собственный вектор колебания um,p показыва-
0.012
ет направление и безразмерную амплитуду смеще-
ния данной частицы p в моде колебаний m. Длину
0.010
собственного вектора можно понимать как количе-
ственное описание того, насколько энергично части-
0.008
ца участвует в данной моде колебаний. На рис. 4A-
0.006
C положения точек соответствуют координатам час-
тиц в конечной структуре, а яркость точек отражает
0.004
величину собственного вектора колебания частиц в
данной моде колебаний. Таким образом, рис. 4A-C
0.002
изображает пространственное распределение приве-
0.000
денной к безразмерному виду амплитуды колебаний
0
25
50
75
100
125
частиц в различных модах. На рис. 4D-F показан
/2
,
Hz
радиальный профиль амплитуды собственного век-
тора колебаний частиц в различных модах, полу-
Рис. 3. Спектральная плотность энергии колебаний частиц
ченный из соответствующих графиков на рис. 4A-C.
SPD(ω) в системе из конечного числа заряженных частиц
Согласно рис. 4, при увеличении частоты моды коле-
в ловушке. Наблюдается ограничение спектра максималь-
баний колебания локализуются в области с большей
ной частотой колебаний частиц ωmax/2π ≈ 125 Гц
плотностью (в рассматриваемой системе в цен-
тре). Этот результат дополняет вывод, сделанный
ранее по спектральной плотности энергии (рис. 3):
вращательные моды исключаются из рассмотрения
колебания наивысшей доступной частоты ωmax тес-
в данной работе, поскольку именно колебательные
но связаны с минимальным расстоянием между час-
моды представляют интерес.
тицами в системе.
Нормальные частоты и направления колебаний
Из работ [26,27] известна связь среднего межча-
частиц могут быть найдены как собственные числа
стичного расстояния с расстоянием до центра струк-
и векторы динамической матрицы
E (4) [42]. При
туры:
расчете динамической матрицы используются коор-
1
Δlocal ∝ n-1/2(r) =
,
(5)
динаты всех частиц системы, т.е. локальный подход
a - br2
в данном разделе не применяется.
где a, b коэффициенты, определяемые параметра-
В результате анализа нормальных мод получены
ми системы. Следовательно, имеется тесная связь
частоты колебаний частиц в системе ωm (m но-
между максимальной частотой колебаний частицы
мер моды), которые представлены на рис. 3 в виде
в системе и расстоянием от данной частицы до цен-
спектра колебаний частиц. Спектр колебаний, полу-
тра. На рис. 5 показано, что максимальная доступ-
ченный таким образом из анализа нормальных мод,
ная частицам частота колебаний ωmax уменьшается
принято называть спектральной плотностью энер-
от центра к краю системы. Таким образом, происхо-
гии (SPD(ω)) [35]. Диапазон частот, проявляющихся
дит локализация колебаний высокочастотных мод в
в спектре колебаний частиц, может быть оценен из
центре конечной системы.
структурных характеристик системы. В кристалле
Согласно стандартной формулировке метода
могут распространяться колебательные процессы с
анализа нормальных мод, амплитуда собственного
длиной волны из диапазона от 2Δ0 до 2D0, где Δ0
вектора um,p приведена к безразмерному виду
межчастичное расстояние, D0 линейный раз-
таким образом, чтоm |um,p|2 = 1. В случае, когда
мер системы. В таком случае максимальная часто-
частицы испытывают тепловые колебания при
та колебаний частиц ωmax в спектре SPD(ω) зави-
произвольной температуре T , амплитуда этих коле-
сит от числа частиц в системе и определяется мини-
баний может быть получена путем перенормировки
мальным межчастичным расстоянием Δminlocal, о чем
суммарной длины собственных векторов um,p на
упоминается в том числе в работе [34]. Связь мак-
абсолютную температуру по шкале Кельвина [45].
симальной частоты и минимального межчастичного
В результате имеем выражение для амплитуды
расстояния может быть в первом приближении оце-
тепловых колебаний частицы p в конечной системе:
нена из закона дисперсии, который хорошо изучен
для однородной системы Юкавы как в длинновол-
∑
∑
DkBTp 1
Ap =
ũ m,p
=
u m,p
 ,
(6)
новом, так и в коротковолновом пределе [43, 44].
M ωm
m
m
995
И. В. Воронов, В. С. Николаев, А. В. Тимофеев
ЖЭТФ, том 162, вып. 6 (12), 2022
100
50
3
2
0.4
0.5
0.6
0.7
0.8
0.9
1.0
r/R0
Рис. 5. Сверху показан нормированный на радиус системы
R0 радиальный профиль максимальной частоты колеба-
ний частиц в системе ωmax частоты, которой ограничен
спектр на рис. 3. Снизу показан нормированный таким же
образом радиальный профиль амплитуды тепловых коле-
баний частиц в системе, рассчитанной методом анализа
нормальных мод с помощью формулы (6). Штриховые ли-
нии проведены для наглядности
лученный из анализа нормальных мод, показан на
рис. 5: амплитуда тепловых колебаний в квазигар-
моническом приближении увеличивается от центра
к краю системы. Этот результат согласуется с ра-
ботой [32], где амплитуда тепловых колебаний заря-
женных частиц в ловушке рассчитана с применени-
ем метода классической молекулярной динамики.
Рис. 4. Распределение длин собственных векторов коле-
баний всех частиц в конечной системе в трех различных
4. ИССЛЕДОВАНИЕ ДИНАМИЧЕСКИХ
модах колебаний m, соответствующих частотам колебаний
СВОЙСТВ СИСТЕМ ЗАРЯЖЕННЫХ
ωm/2π = {3, 101, 125} Гц. На рис. A-C показано распреде-
ЧАСТИЦ В ЛОВУШКЕ В РАМКАХ
ление амплитуды собственного вектора каждой частицы по
ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ
структуре. Чем ярче изображена точка, тем больше дли-
на собственного вектора колебания соответствующей час-
В численном моделировании, аналогичном ме-
тицы в данной моде. На рис. D-F показана зависимость
тоду классической молекулярной динамики, систе-
длины собственного вектора частицы от радиального рас-
ма из конечного числа заряженных частиц рассмат-
стояния, соответствующая графикам A-C для данной мо-
ривается в рамках описанного в разд. 2 локально-
ды колебаний. Расстояния нормированы на радиус систе-
мы R0
го подхода. Численное моделирование позволяет ис-
следовать систему из значительно большего числа
частиц и в более широком диапазоне параметров,
где D размерность пространства, kB постоянная
чем метод анализа нормальных мод, примененный
Больцмана, M масса частицы, а Tp ≡ T, посколь-
в разд. 3. В каждом выделенном кольце конечной
ку температура в системе постоянна. Суммирование
системы (см. рис. 1) рассчитывается автокорреля-
ведется по модам m = {ftrans + frot + 1, ..., DN}. Та-
ционная функция скорости (АКФС):
ким образом, из суммирования исключаются посту-
〈vi(0) · vi(t)〉i
пательные и вращательные моды, присутствующие
VACF(t) =
,
(7)
〈(vi(0))2i
среди решений задачи на собственные значения ди-
намической матрицы (4). Радиальный профиль ам-
где vi(t)
скорость частицы i в момент време-
плитуды тепловых колебаний частиц в системе, по-
ни t, а усреднение ведется по частицам в равно-
996
ЖЭТФ, том 162, вып. 6 (12), 2022Исследование спектральных свойств пространственно-неоднородной системы. . .
весном ансамбле. Затем в каждом кольце вычис-
ляется плотность колебательных состояний частиц
VDOS(ω) путем применения преобразования Фурье
к АКФС [46]:
+∞
1
VDOS(ω) =
VACF(t) exp(-iωt) dt.
(8)
-∞
Сравнение плотности колебательных состояний
VDOS(ω) со спектральной плотностью энергии
SPD(ω) при различных значениях температу-
ры можно видеть на рис. 6. Частоты на рис. 6
нормированы на среднее значение
¾плазменно-
пылевой¿ частоты ωpd
=
Q2/(M〈Δlocal3), где
Q электрический заряд частицы, M ее мас-
са, а 〈Δlocal〉 = 1/
N/S =
πR20/N
среднее
межчастичное расстояние в системе, S площадь
структуры. Спектры колебаний частиц VDOS(ω)
и SPD(ω) имеют схожий вид: наблюдаются два
характерных пика. При этом наилучшее соответ-
ствие VDOS и SPD наблюдается в низкочастотной
Рис. 6. Сравнение спектров колебаний всех частиц в ко-
области спектра. Отметим, что при обоих значе-
нечной системе SPD(ω) (рис. 3) и VDOS(ω) (формула (8))
ниях температур спектры SPD трудно различимы.
при двух значениях температуры: → T = 100 K (A) и 300 K
Спектры VDOS, напротив, при разных температу-
(B). Частоты нормированы на ωpd. Наблюдается схожесть
рах различаются значительно. Это связано с тем,
положения характерных пиков, в особенности низкочастот-
ного. Однако можно видеть различное поведение спектров
что метод анализа нормальных мод не выходит за
в высокочастотной части. Видимая непрерывность спек-
рамки квазигармонического приближения, а ре-
тров связана только со сглаживанием исходных дискрет-
зультаты его применения определяются силовыми
ных спектров для наглядности
константами и положениями частиц, т.е. зависят от
структурных характеристик. Спектр VDOS, в свою
ваются приблизительно на одной и той же частоте
очередь, получен напрямую из расчета динамики
ωmax, которая может быть оценена как
частиц, поэтому позволяет более точно описать вид
спектра колебаний частиц в широком диапазоне
ωmax ≈ CωlocalE,
(9)
температур.
7
На рис.
показаны спектры колебаний час-
где величина C не зависит от радиального расстоя-
тиц VDOS(ω) в различных кольцах конечной си-
ния.
стемы. Частоты нормированы на среднее значе-
Радиальный профиль максимальной частоты
ние ¾плазменно-пылевой¿ частоты ωpd. Наблюда-
ωmax в различных безразмерных единицах пред-
ется ограничение максимальной частоты спектра
ставлен на рис. 9. Отметим, что величина C из
колебаний частиц величиной ωmax. Значение ωmax
выражения (9) приблизительно равна 1.5 в рас-
уменьшается с увеличением межчастичного рассто-
сматриваемой в работе структуре. Величина
ωE
яния.
среднее значение эйнштейновской частоты, рас-
Спектры колебаний частиц VDOS(ω) в различ-
считанное для среднего значения межчастичного
ных кольцах конечной системы могут быть нор-
расстояния 〈Δlocal〉 c использованием выражения
мированы на локальную величину эйнштейновской
из работы [47]. Причины ограничения максималь-
частоты ωlocalE, которая зависит от локального зна-
ной частоты спектра обсуждаются в разд.
3
чения межчастичного расстояния Δlocal. Для нор-
величина максимальной частоты ωmax тесно свя-
мировки используются результаты расчета эйн-
зана с межчастичным расстоянием в кольце и,
штейновской частоты для однородных юкавовских
соответственно, с радиальным расстоянием. Таким
систем из работы [47]. На рис. 8 показаны спек-
образом, неоднородность спектра колебаний частиц,
тры колебаний частиц в различных кольцах конеч-
определяющего динамические свойства системы,
ной системы с описанной нормировкой на величи-
напрямую связана с неоднородностью структурных
ну ωlocalE. При такой нормировке спектры ограничи-
характеристик.
997
И. В. Воронов, В. С. Николаев, А. В. Тимофеев
ЖЭТФ, том 162, вып. 6 (12), 2022
Рис.
7. (В цвете онлайн) Спектр колебаний частиц
Рис.
8. (В цвете онлайн) Спектр колебаний частиц
VDOS(ω) в областях конечной системы (синим цветом):
VDOS(ω) в областях конечной системы (синим цветом):
центральной (A), промежуточной (B) и периферийной (C).
центральной (A), промежуточной (B) и периферийной (C).
Частоты нормированы на среднее значение ¾плазменно-
Частоты нормированы на локальное значение эйнштей-
пылевой¿ частоты ωpd. Спектры ограничены величиной
новской частоты ωlocalE. Спектры ограничены величиной
максимальной частоты ωmaxpd, которая уменьшается с
максимальной частоты ωmaxEocal, которая не зависит от
увеличением межчастичного расстояния. Оранжевой ли-
радиального расстояния. Оранжевой линии соответствует
нии соответствует сглаженный для наглядности спектр
сглаженный для наглядности спектр
Как было отмечено выше, спектр колебаний час-
тиц на рис. 6 состоит из двух пиков. Это можно свя-
зать с тем, что в структуре происходят колебатель-
ные процессы двух принципиально разных типов,
которые можно сравнить с волнами сдвига (низко-
частотный пик) и сжатия (высокочастотный пик) в
¾бесконечных¿ системах [35]. Сглаженные спектры
на рис. 7 наглядно показывают, что наличие двух
пиков в спектре выделенной подсистемы наблюда-
ется при любом разбиении системы на кольца. Со-
хранение двух пиков в спектре колебаний позволяет
сделать вывод о том, что на качественном уровне
спектральные характеристики можно сопоставлять
с бесконечным юкавовским кристаллом не только в
конечной структуре целиком, но и в ее отдельных
областях.
Рис. 9. (В цвете онлайн) Приведенный к безразмерно-
Амплитуда тепловых колебаний частиц может
му виду радиальный профиль наивысшей доступной час-
быть рассчитана по спектру колебаний частиц с при-
тоты колебаний частиц ωmax (см. рис. 5). Оранжевыми
менением низкотемпературного приближения тео-
квадратами (■) обозначена величина ωmaxEocal. Синими
рии динамики кристаллической решетки в каждом
кругами (•) обозначена величина ωmaxE. Фиолетовыми
отдельном кольце конечной системы:
звездами (⋆) обозначена величина ωmaxp
d
. Зелеными
ромбами (♦) обозначена величина ωmaxpd. Штриховые
DkBT
1
линии проведены для наглядности
〈u2〉 =
,
(10)
M ω2
998
ЖЭТФ, том 162, вып. 6 (12), 2022Исследование спектральных свойств пространственно-неоднородной системы. . .
числа частиц в ловушке, но и однородные системы
в периодических граничных условиях.
Предлагаемый локальный подход к описанию ко-
нечных систем позволяет сравнить с точки зрения
динамических свойств различные области конечной
системы с соответствующими им по плотности одно-
родными системами в ПГУ. Структура из конечно-
го числа частиц мысленно разбивается на участки
с близкой к постоянной плотностью частиц. Каж-
дая из полученных таким образом квазиоднородных
подсистем сравнивается с однородной структурой в
ПГУ (см. рис. 1). Это позволяет сделать выводы о
границах применимости подхода к описанию конеч-
ных структур, который заключается в переносе хо-
рошо изученных свойств систем в ПГУ на близкую
Рис. 10. Зависимости амплитуды тепловых колебаний час-
к однородной внутреннюю область структур из ко-
тиц в конечной системе от безразмерного параметра экра-
нечного числа частиц.
нирования κ∆local. Круги (•) соответствуют RMSD (12),
квадраты (■) соответствуют теоретической формуле (10),
треугольники (▴) обозначают результаты анализа нор-
5.1. Сравнение спектров колебаний частиц
мальных мод (6). Штриховая линия соответствует аппрок-
симации 〈u2(r)〉 ∝ ∆local(r) exp (κ∆local(r)) из работы [32]
На рис. 11 показан результат применения ло-
кального подхода к сравнению спектров колебаний
где используется величина 〈1/ω2〉, полученная по
частиц в конечной системе в ловушке со спектрами
спектру колебаний согласно следующему выраже-
соответствующих по плотности однородных систем
нию [48, 49]:
в ПГУ.
Особенности спектров колебаний частиц в конеч-
1
VDOS(ω)ω-2
=
(11)
ной системе уже подробно описаны в разд. 3 и 4. В
ω2
VDOS(ω) dω
данном разделе мы отметим особенности спектров,
проявляющиеся при сравнении систем из конечного
В некоторых случаях может быть более удоб-
числа частиц и однородных бесконечных систем в
но провести расчет амплитуды тепловых колебаний
ПГУ.
как среднеквадратического отклонения частиц от
Как в конечных структурах, так и в ¾бесконеч-
их положений равновесия:
ных¿ системах спектр колебаний частиц ограничи-
〈ui〉 =
〈(ri(t) - 〈rit)2t.
(12)
вается предельной частотой, причем величины ωmax
имеют близкие друг к другу значения в обеих си-
Результаты расчетов амплитуды тепловых ко-
стемах при одинаковом межчастичном расстоянии.
лебаний заряженных частиц в ловушке, проведен-
Этот результат еще раз подчеркивает вывод о том,
ных по формулам (6), (10) и (12), представлены на
что максимальная частота колебаний частиц опре-
рис. 10. Показана неоднородность радиального про-
деляется минимальным межчастичным расстояни-
филя амплитуды тепловых колебаний. Наблюдается
ем в системе.
согласие предложенных методов расчета амплитуды
тепловых колебаний как между собой, так и с тео-
ретической оценкой из работы [32]:
5.2. Сравнение амплитуд колебаний частиц
〈u2(r)〉 ∝ Δlocal(r) exp (κΔlocal(r)).
(13)
В настоящем разделе рассматриваются (а) цен-
тральная часть системы из конечного числа частиц
и (б) система в ПГУ при той же плотности час-
тиц. Двумерная постановка задачи требует исследо-
5. СРАВНЕНИЕ ДИНАМИЧЕСКИХ
вания влияния числа частиц на величину амплиту-
СВОЙСТВ КОНЕЧНЫХ И БЕСКОНЕЧНЫХ
ды тепловых колебаний частиц, поскольку в двумер-
СИСТЕМ ЗАРЯЖЕННЫХ ЧАСТИЦ
ных системах наблюдается логарифмическая расхо-
С помощью приведенных в разд. 2 и 4 методов
димость амплитуды тепловых колебаний по числу
можно исследовать не только системы из конечного
частиц [50-52]. Также отметим, что существенным
999
И. В. Воронов, В. С. Николаев, А. В. Тимофеев
ЖЭТФ, том 162, вып. 6 (12), 2022
Рис.
12.
Зависимость параметра Линдеманна
〈u2(r)〉/∆local от полного числа частиц в системе.
Каждая точка соответствует отдельной системе. Закра-
шенные символы (•) соответствуют центральной области
конечных систем в ловушке, светлые (□) однородным
структурам в ПГУ
этом в центре конечной системы наблюдается гекса-
гональная решетка, аналогичная решетке в струк-
туре в периодических граничных условиях. В ра-
боте [53] сравнение конечной системы с системой в
ПГУ проводится не только с точки зрения струк-
турных свойств, но и с точки зрения упругих ха-
рактеристик, в частности, модуля упругости. Вы-
двинуто предположение о том, что при числе час-
тиц большем 3200 внутренняя область системы из
Рис. 11. (В цвете онлайн) Спектры колебаний частиц в
конечного числа частиц может описываться с помо-
выделенных подсистемах конечной системы (соответству-
щью теоретических представлений, существующих
ют буквам A, C) и соответствующих по плотности одно-
родных системах (соответствуют буквам B, D). Частоты
для систем в ПГУ. Этот результат получен для зна-
нормированы на локальное значение эйнштейновской час-
чений параметра экранировки потенциала из диа-
тоты ωlocalE. Спектры ограничены величиной максималь-
пазона 0.2 < κΔlocal < 0.5. В работе [35] отмечено
ной частоты ωmaxEocal. Спектр выделенной подсистемы
сходство спектров колебаний частиц в конечной си-
с постоянной плотностью ограничивается приблизительно
стеме и в системе в ПГУ. Показано наличие в конеч-
на той же частоте, что и соответствующая однородная си-
ной системе колебательных процессов, аналогичных
стема в ПГУ. Оранжевой линии соответствует сглаженный
волнам сжатия и сдвига в ¾бесконечной¿ системе.
для наглядности спектр
Тем не менее полученные в приведенных работах ре-
зультаты не позволяют разрешить вопрос о возмож-
отличием конечной системы от системы в ПГУ яв-
ности переноса динамических свойств систем в ПГУ
ляется отсутствие в последней удерживающей ло-
на конечные структуры в связи с ограниченностью
вушки. В данном разделе проводится сравнение ам-
рассмотренного диапазона числа частиц N ≲ 103.
плитуд тепловых колебаний частиц в близкой к од-
В приведенных работах также не обсуждается во-
нородной центральной области системы в ловушке
прос обеспечения неизменного межчастичного рас-
и в системе в ПГУ.
стояния в центре конечной системы. Этот вопрос
Вопрос сравнения свойств конечной системы в
имеет принципиальное значение, поскольку при уве-
личении числа частиц в конечной системе при по-
ловушке и системы в ПГУ поднимался ранее [35,53,
54]. В работе [54] отмечено, что конечноразмерные
стоянных Q, κ, α межчастичное расстояние в цен-
тральной области структуры монотонно уменьшает-
эффекты в двумерной системе становятся незначи-
тельными с точки зрения влияния на структурные
ся, что показано в работах [26, 30]. Таким образом,
при увеличении числа частиц требуется уменьшать
свойства при числе частиц в системе N > 50. При
1000
ЖЭТФ, том 162, вып. 6 (12), 2022Исследование спектральных свойств пространственно-неоднородной системы. . .
параметр ловушки α для обеспечения неизменности
нечного числа одноименно заряженных частиц в па-
межчастичного расстояния в центральной области
раболической электростатической ловушке методом
конечной структуры. При этом влияние ловушки
классической молекулярной динамики. В локальном
на систему из конечного числа частиц уменьшает-
подходе рассмотрены спектры и амплитуды тепло-
ся, и ее динамические свойства постепенно сходятся
вых колебаний частиц. Показано, что от центра к
к свойствам системы в ПГУ.
краю структуры происходит смещение спектра ко-
На рис. 12 представлена зависимость параметра
лебаний частиц в низкочастотную область, что свя-
Линдеманна
〈u2〉/Δlocal от полного числа частиц в
зано с увеличением межчастичного расстояния в си-
системе. Величина параметра Линдеманна в конеч-
стеме от центра к краю. Неоднородность радиально-
ной системе рассчитывается для центральной обла-
го профиля амплитуды тепловых колебаний частиц
сти, а в периодических граничных условиях для
продемонстрирована как методом анализа нормаль-
всей системы. Как можно видеть из графика, ампли-
ных мод, так и с помощью расчетов методом моле-
туда тепловых колебаний частиц в системе в ПГУ
кулярной динамики. Результаты расчетов величины
логарифмически растет с ростом числа частиц. Этот
амплитуды тепловых колебаний частиц различны-
эффект подтверждается результатами работ [50-52].
ми методами согласуются друг с другом при низких
В то же время амплитуда тепловых колебаний час-
температурах.
тиц в центре конечной системы монотонно умень-
С помощью численного моделирования с приме-
шается с ростом полного числа частиц в структу-
нением локального подхода сравнивается централь-
ре. При значении числа частиц N ≈ 105 величина
ная подсистема конечной структуры из заряженных
амплитуды тепловых колебаний частиц в конечной
частиц в ловушке и соответствующая ей по плотно-
системе сравнивается с величиной в системе в ПГУ
сти ¾бесконечная¿ однородная система заряженных
в пределах ошибки расчета. Отметим, что структу-
частиц в периодических граничных условиях. Даже
ры из 105 частиц редко наблюдаются в лаборатор-
при условии совпадения межчастичных расстояний
ных экспериментах [3]. Таким образом, перенос хо-
при числе частиц меньшем 105 наблюдается разли-
рошо изученных свойств однородных систем в ПГУ
чие динамических свойств центральной подсистемы
на структуры из конечного числа частиц может при-
конечной структуры и соответствующей однород-
водить к заметным систематическим ошибкам при
ной системы в рассматриваемых условиях. Показа-
описании пылевой подсистемы в плазме в условиях
но, что при малом числе частиц амплитуды тепло-
лабораторного эксперимента.
вых колебаний частиц существенно различны в цен-
тре конечной системы в ловушке и в совпадающей с
ней по плотности однородной системе в периодиче-
6. ЗАКЛЮЧЕНИЕ
ских граничных условиях. Таким образом, перенос
динамических свойств однородных систем в перио-
В настоящей работе исследованы динамические
дических граничных условиях на системы из конеч-
свойства двумерной системы из конечного числа
ного и недостаточно большого числа частиц может
взаимодействующих по экранированному кулонов-
приводить к ошибкам как при разработке теорети-
скому потенциалу одноименно заряженных частиц
ческих моделей, так и при обработке эксперимен-
в ловушке: спектры, частоты и амплитуды тепло-
тальных результатов.
вых колебаний частиц. Система рассматривается
Финансирование. Работа В. С. Николаева бы-
при сильной экранировке (κΔ ≳ 1).
ла поддержана грантом Фонда развития теорети-
Методом анализа нормальных мод показано про-
ческой физики и математики ¾БАЗИС¿. Работа
странственное распределение амплитуд колебаний
А. В. Тимофеева поддержана в рамках Программы
частиц в различных модах. Сделан вывод о том, что
фундаментальных исследований НИУ ВШЭ.
в низкочастотные колебания вовлечены все частицы
системы, в то время как высокочастотные колебания
локализуются в области с наибольшей плотностью
частиц. В рассматриваемой двумерной системе в па-
ЛИТЕРАТУРА
раболической ловушке этой области соответствует
центральная часть структуры.
1. H. Yukawa, Proc. Phys.-Math. Soc. Jpn. 17,
48
(1935).
Для обобщения результатов анализа нормаль-
ных мод на более широкий температурный диапазон
проведено численное моделирование системы из ко-
2. P. Pieranski, Contemp. Phys. 24, 25 (1983).
1001
13
ЖЭТФ, вып. 6 (12)
И. В. Воронов, В. С. Николаев, А. В. Тимофеев
ЖЭТФ, том 162, вып. 6 (12), 2022
3.
V. E. Fortov and G. Morfill, Complex and Dusty
22.
D. N. Polyakov, V. V. Shumova, and L. M. Vasilyak,
Plasmas: From Laboratory to Space, CRC Press
Plasma Sources Sci. Technol. 26, 08LT01 (2017).
(2009).
23.
T. S. Ramazanov, K. N. Dzhumagulova, A. N.
4.
С. Н. Антипов, Э. И. Асиновский, А. В. Кирил-
Jumabekov, and M. K. Dosbolayev, Phys. Plasmas
лин, С. А. Майоров, В. В. Марковен, О. Ф. Петров,
15, 053704 (2008).
В. Е. Фортов, ЖЭТФ 133, 948 (2008).
24.
U. Konopka, G. E. Morfill, and L. Ratke, Phys. Rev.
5.
D. N. Polyakov, V. V. Shumova, and L. M. Vasilyak,
Lett. 84, 891 (2000).
Surf. Eng. Appl. Electrochem. 51, 143 (2015).
25.
D. N. Polyakov, V. V. Shumova, L. M. Vasilyak, and
V. E. Fortov, Phys. Lett. A 375, 300 (2011).
6.
V. N. Tsytovich, N. G. Gusein-Zade, and A. M.
Ignatov, Plasma Phys. Rep. 43, 981 (2017).
26.
C. Henning, H. Baumgartner, A. Piel, P. Ludwig, V.
Golubnichiy, M. Bonitz, and D. Block, Phys. Rev. E
7.
О. С. Ваулина, Э. А. Саметов, ЖЭТФ 154, 407
74, 056403 (2006).
(2018).
27.
C. Henning, P. Ludwig, A. Filinov, A. Piel, and M.
8.
В. Е. Фортов, В. С. Филинов, А. П. Нефедов,
Bonitz, Phys. Rev. E 76, 036404 (2007).
О. Ф. Петров, А. А. Самарян, А. М. Липаев,
ЖЭТФ 111, 889 (1997).
28.
A. M. Ignatov, Plasma Phys. Rep. 46, 410 (2020).
9.
В. Ю. Карасев, А. Ю. Иванов, Е. С. Дзлиева,
29.
N. P. Kryuchkov, Phys. Rev. Lett. 121, 075003
А. И. Эйхвальд, ЖЭТФ 133, 460 (2008).
(2018).
10.
М. М. Васильев, Л. Г. Дьячков, С. Н. Антипов,
30.
H. Totsuji, C. Totsuji, and K. Tsuruta, Phys. Rev. E
О. Ф. Петров, В. Е. Фортов, Письма в ЖЭТФ 86,
64, 066402 (2001).
414 (2007).
31.
M.G. Hariprasad, P. Bandyopadhyay, G. Arora, and
11.
О. С. Ваулина, Е. А. Лисин, А. В. Гавриков,
A. Sen, Phys. Plasmas 25, 123704 (2018).
. Ф. Петров, В. Е. Фортов, ЖЭТФ 137, 751 (2010).
32.
V. S. Nikolaev and A. V. Timofeev, Phys. Plasmas
26, 073701 (2019).
12.
G. I. Sukhinin, A. V. Fedoseev, S. N. Antipov,
O. F. Petrov, and V. E. Fortov, Phys. Rev. E: Stat.,
33.
V. S. Nikolaev and A. V. Timofeev, Phys. Plasmas
Nonlinear, Soft Matter Phys. 87, 013101 (2013).
28, 033704 (2021).
13.
D. Samsonov, S. K. Zhdanov, R. A. Quinn,
34.
V. A. Schweigert and F. M. Peeters, Phys. Rev. B 51,
S. I. Popel, and G. E. Morfill, Phys. Rev. Lett. 92,
7700 (1995).
255004 (2004).
35.
A. Melzer, Phys. Rev. E 7, 016411 (2003).
14.
A. M. Ignatov, Plasma Phys. Rep. 31, 46 (2005).
36.
V. Nosenko and J. Goree, Phys. Rev. Lett. 93, 155004
15.
А. В. Филиппов, В. В. Решетняк, А. Н. Старостин,
(2004).
И. М. Ткаченко, В. Е. Фортов, Письма в ЖЭТФ
37.
Z. Donkó, P. Hartmann, and J. Goree, Mod. Phys.
110, 658 (2019).
Lett. B 21, 1357 (2007).
16.
Р. Е. Болтнев, ЖЭТФ 153, 679 (2018).
38.
P. Hartmann, M. C. Sándor, A. Kovács, and Z.
Donkó, Phys. Rev. E 84, 016404 (2011).
17.
O. S. Vaulina and X. G. Koss, Phys. Rev. E 92,
042155 (2015).
39.
H. C. Lee, D. Y. Chen, and B. Rosenstein, Phys. Rev.
E 56, 4596 (1997).
18.
O. C. Ваулина, И. И. Лисина, Е. А. Лисин, ЖЭТФ
148, 819 (2015).
40.
H. C. Lee and B. Rosenstein, Phys. Rev. E 55, 7805
(1997).
19.
M. O. Robbins, K. Kremer, and G. S. Grest, J. Chem.
Phys. 88, 3286 (1988).
41.
Б. А. Клумов, УФН 180, 1095 (2010).
20.
S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, J.
42.
A. Melzer, M. Klindworth, and A. Piel, Phys. Rev.
Chem. Phys. 105, 7641 (1996).
Lett. 87, 115002 (2001).
21.
S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin,
43.
L. Bonsall and A. A. Maradudin, Phys. Rev. B 15,
Phys. Rev. E 56, 4671 (1997).
1959 (1977).
1002
ЖЭТФ, том 162, вып. 6 (12), 2022Исследование спектральных свойств пространственно-неоднородной системы. . .
44. S. H. X. W. S. Nunomura, J. Goree, S. Hu, X. Wang,
49. S. A. Khrapak, Phys. Rev. Res. 2, 012040 (2020).
and A. Bhattacharjee, Phys. Rev. E 65, 066402
(2002).
50. Л. Д. Ландау, Е.М. Лифшиц, Теоретическая фи-
зика, т. V, Статистическая физика, ч. 1., Наука,
45. Normal Mode Analysis: Theory and Applications to
Москва (1976).
Biological and Chemical Systems, ed. by Q. Cui and
I. Bahar, CRC press (2005).
51. S. Toxværd, Phys. Rev. Lett. 51, 1971 (1983).
46. P. Vashishta, R. K. Kalia, A. Nakano, and J. P. Rino,
J. Appl. Phys. (Melville, NY, U. S.) 101, 103515
52. N. D. Mermin, Phys. Rev. 176, 250 (1968).
(2007).
47. S. Khrapak and B. Klumov, Phys. Plasmas 25,
53. T. E. Sheridan, Phys. Plasmas 14, 032108 (2007).
033706 (2018).
48. D. A. Young and B. J. Alder, J. Chem. Phys. 60,
54. Yu. Ivanov and A. Melzer, Phys. Rev. B 79, 036402
1254 (1974).
(2009).
1003
13*