Физика плазмы, 2022, T. 48, № 5, стр. 457-462

1. ВВЕДЕНИЕ

Броуновское движение в системах взаимодействующих частиц широко распространено в природе и наблюдается, например, в биологических и коллоидных растворах, в плазме продуктов сгорания, в атмосфере Земли и т.д. [1–4]. Однако на настоящий момент аналитические модели для броуновской динамики частиц разработаны только для двух простейших случаев: невзаимодействующие частицы; и одиночная заряженная частица, движение которой ограничено потенциальным полем ловушки. Анализ этих задач не позволяет исследовать влияние числа N взаимодействующих частиц на характер их броуновского движения. Для этой цели широко используется численное моделирование.

Большинство работ, посвященных численному исследованию свойств комплексной (пылевой) плазмы, опираются на решение уравнений Ланжевена для пылевых частиц с некоторой эффективной температурой Т большей, чем температура T_n окружающего их нейтрального газа [1–3]. Данное обстоятельство связано с тем, что стохастическая энергия пылевых частиц (их “кинетическая температура” Т) в таких условиях может достигать ~0.3–10 эВ, что значительно выше T_n. Механизмы “аномального разогрева” пылевых частиц обычно связывают с временными и/или пространственными изменениями их зарядов в объеме неоднородной плазмы [1]. (Флуктуации зарядов пылевых частиц, вызванные случайной природой ионных и электронных токов, заряжающих эти частицы, присущи любым типам комплексной плазмы.)

Экспериментальный, теоретический и численный анализ теплового движения взаимодействующих пылевых частиц в протяженных и ограниченных ансамблях, формирующихся в газоразрядной плазме, представлены в работах [5–11]. Отметим, что в тлеющих разрядах (как переменного, так и постоянного токов) в центре газоразрядных камер наблюдается некоторое превышение концентрации ионов плазмы над ее электронной концентрацией [12]. Данное обстоятельство приводит к формированию эффективных ловушек для отрицательно заряженных частиц пыли [1–3]. В отличие от плазмы высокочастотного емкостного разряда, где зачастую формируются слоистые пылевые структуры, объемные облака заряженных пылевых частиц обычно наблюдаются в стратах тлеющего разряда постоянного тока, или в индукционных высокочастотных разрядах [1, 13–20].

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

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

2. ОСНОВНЫЕ ПАРАМЕТРЫ И СООТНОШЕНИЯ

Рассмотрим движение одной заряженной частицы с массой M и зарядом Q в поле гармонической ловушки c характерной частотой ω₀ = = (Qα/M)^1/2 под воздействием случайной силы F_b, которая является источником стохастической (тепловой) энергии частиц; здесь α ≡ α_x = dE_x/dx – величина градиента внешнего электрического поля E_x, ловушки в направлении оси x. Средний квадрат отклонений $\langle {{x}^{2}}(t)\rangle $ такой частицы для одной степени свободы от ее начального положения за время t и автокорреляционная функция скорости $\left\langle {{{V}_{x}}(0){{V}_{x}}(t)} \right\rangle $ могут быть представлены в форме [1, 5]

Здесь и далее T – температура частиц в энергетических единицах, ν_fr – коэффициент трения частиц, $\psi = {{(1 - 4{{\xi }^{2}})}^{{1/2}}}{\text{/}}2$, ξ = ω₀/ν_fr, а угловые скобки $\left\langle \ldots \right\rangle $ обозначают усреднение по всем отрезкам времени, равным t. (При условии 1 < 4ξ², $\psi = {{(4{{\xi }^{2}} - 1)}^{{1/2}}}{\text{/}}2$, а функции ch и sh переходят в cos и sin соответственно.) С ростом времени при t → ∞ величина $\langle {{x}^{2}}(t)\rangle \to 2T{\text{/}}(M\omega _{0}^{2})$.

Формула (2) может быть получена как независимым путем, так и из соотношения

Отметим, что соотношение (3) имеет место как для одной частицы в ловушке, так и для протяженных систем различной конфигурации [5, 6].

Для анализа теплового движения заряженных пылевых частиц удобно использовать отношение функции массопереноса $D(t) = \langle {{x}^{2}}(t)\rangle {\text{/}}2t$ к коэффициенту диффузии ${{D}_{0}} = T{\text{/}}\left( {{{\nu }_{{fr}}}M} \right)$ невзаимодействующих частиц: $D(t){\text{/}}{{D}_{0}}$. Легко увидеть, что на малых временах наблюдения, ${{\nu }_{{fr}}}t \ll (1 + \xi )$, для частицы в ловушке характерен баллистический режим движения ($\langle {{x}^{2}}(t)\rangle \approx T{{t}^{2}}{\text{/}}М$, $D(t) = \langle {{x}^{2}}(t)\rangle {\text{/}}$ ${\text{/}}2t \propto t$), как и в случае невзаимодействующих частиц, см. рис. 1а. Численный анализ показывает, что максимум функции D(t) наблюдается в диапазоне времен, соответствующем ${{\omega }_{0}}t$ от 1 до 2 вне зависимости от величины ${{\nu }_{{fr}}}$ (см. рис. 1б).

Рис. 1.

Зависимость функций D(t)/D₀ от ν_fr t (а) и от ω₀t (б): 1 – баллистический режим, D(t) ∝ t; 2 – невзаимодействующие частицы [4, 5]; гармонический осциллятор (1) при различных параметрах ξ: 3 – 0.033; 4 – 0.38; 5 – 3; 6 – 15.

Для описания физических свойств исследуемых структур введем также параметр $\Gamma = $ $ = {{Q}^{2}}{{n}^{{1/3}}}{\text{/}}T \equiv {{Q}^{2}}{\text{/}}{{l}_{p}}T$, где $n = l_{p}^{{ - 3}}$ – концентрация частиц, а l_p – среднее расстояние между ними. При этом в линейном электрическом поле $n \cong 3\alpha {\text{/}}(4\pi Q)$, а ${{l}_{p}} \cong {{(4\pi Q{\text{/}}3\alpha )}^{{1/3}}}$. Тогда для оценки радиуса ограниченной структуры имеет место соотношение $R \cong {{(3N{\text{/}}4\pi n)}^{{1/3}}}$, где N – число частиц в системе. Реальный (эффективный) радиус ансамбля будет несколько увеличиваться с ростом температуры частиц.

Эволюция среднеквадратичного смещения D(t) в однородных и протяженных неидеальных системах исследовалась численно в работах [9, 11] для частиц, взаимодействующих посредством экранированного кулоновского потенциала (типа Юкавы)

Здесь r – расстояние между двумя заряженными частицами, λ – длина экранирования, κ = l_p/λ – параметр экранирования, а l_р – среднее межчастичное расстояние. (Для моделирования таких систем использовались периодические граничные условия [1, 2].)

Анализ численных исследований для протяженных трехмерных систем с экранированным кулоновским потенциалом [9–11], показывает, что при κ = l_p/λ < 6 за их равновесные характеристики отвечают два безразмерных параметра: $\Gamma {\kern 1pt} * = {{Q}^{2}}(1 + \kappa + {{\kappa }^{2}}{\text{/}}2)\exp ( - \kappa ){\text{/}}(T{{l}_{р}})$ и $\xi {\kern 1pt} * = \omega {\kern 1pt} *{\text{/}}{{\nu }_{{fr}}}$; здесь $\omega {\kern 1pt} * = Q{{[(1 + \kappa + {{\kappa }^{2}}{\text{/}}2)\exp ( - \kappa ){\text{/}}(\pi Ml_{p}^{3})]}^{{1/2}}}$. Значение данных параметров, совместно с температурой частиц, определяет пространственную корреляцию, точки фазовых переходов и процессы переноса (диффузию, вязкость) в рассматриваемых системах при условии $20 < \Gamma {\kern 1pt} * < \Gamma _{c}^{*} \approx 106$; здесь $\Gamma _{c}^{*}$ соответствует точке их кристаллизации [10, 11]. Отметим, что при κ < 1 величина $\Gamma {\kern 1pt} * \approx G \equiv {{Q}^{2}}{\text{/}}{{l}_{p}}T$, а $\omega {\kern 1pt} * = Q{\text{/}}{{(\pi Ml_{p}^{3})}^{{1/2}}}$. Расчеты показывают, что диффузионный режим движения частиц в таких структурах наблюдается для времен $t \gg 1{\text{/}}\omega {\kern 1pt} *$, когда их коэффициент диффузии D(t) → D_с = const.

Принимая во внимание ограниченный (конечный) размер структур, исследуемых в настоящей работе, наблюдать диффузионный режим движения частиц возможно только в диапазоне времен ${{t}_{R}} \gg t \gg 1{\text{/}}{{\omega }_{0}}$, где величина t_R определяется размером кластера R, температурой частиц T и характерными частотами (ω₀, ν_fr):

Результаты моделирования динамики ограниченных ансамблей заряженных частиц в электростатической ловушке представлены в следующем разделе. Расчеты проводились в широком диапазоне параметров ξ ~ 0.15–15, соответствующих условиям наблюдения пылевых структур в плазме газовых разрядов [1–3]: диаметр частиц ~1–10 мкм; l_p ~ 500–1000 мкм; Т ~ 0.03–10 эВ; давление буферного газа ~0.1–0.01 Торр.

3. РЕЗУЛЬТАТЫ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ И ИХ ОБСУЖДЕНИЕ

Численное моделирование динамики систем заряженных частиц в линейной и изотропной электростатической ловушке (α ≡ α_x = α_y = α_z) выполнялось методом молекулярной динамики Ланжевена. Расчеты проводились для ансамблей, состоящих из N = 50, 150, 450 и 900 частиц. Техника моделирования подробно описана в работах [1, 2]. Шаг интегрирования составлял от Δt ≅ ≅ (40max[ω₀; ν_fr])^–1 до Δt ≅ (100max[ω₀; ν_fr])^–1 в зависимости от начальных условий задачи. Время расчетов t_c после установления равновесия в моделируемых системах варьировалось от ~10³/min[ω₀; ν_fr] до ~10⁴/min[ω₀; ν_fr].

Расчеты проводились для систем частиц с кулоновским взаимодействием в широком диапазоне параметров неидеальности: от Γ ~ 1 до Γ ~ ~ 300. Частота ловушки ω₀ = (Qα/M)^1/2 изменялась от 15 до 45 с^–1, что обеспечивало среднее расстояние между частицами l_p ~ 500–1000 мкм. Параметр ξ = ω₀/ν_fr варьировался за счет коэффициента трения частиц ν_fr в пределах от 0.15 до 15.

Во всех рассмотренных случаях температура частиц не отличалась от заданной, а их функции распределения по скоростям соответствовали распределению Максвелла. При t → ∞ среднеквадратичные смещения центра масс исследуемых ансамблей частиц от их начального положения по всем степеням свободы были равны $\langle {{x}^{2}}(t)\rangle \cong \langle {{y}^{2}}(t)\rangle \cong \langle {{z}^{2}}(t)\rangle $ ≈ $T{\text{/}}(NМ\omega _{0}^{2})$.

Парные корреляционные функции g(l) для ансамблей из N = 450 частиц с различными параметрами Г показаны на рис. 2. Легко увидеть, что первый пик функций g(l) для Г ≥ 0.1 хорошо соответствует величине l_p ≅ (4πQ/3α)^1/3, полученной в приближении однородной системы. Кроме того, численное моделирование показало, что для оценки радиуса неидеальных систем при Г ≥ 1 может быть использовано соотношение R ≅ ≅ (3N/4πn)^1/3, которое имеет место для однородных систем.

Рис. 2.

Парная корреляционная функция g(l/l_p) для ансамблей из N = 450 частиц с различными параметрами Г: 1 – 100; 2 – 10; 3 – 1. Здесь l_p ≅ (4πQ/3α)^1/3.

При высоких температурах (см. рис. 3), а также при малых размерах кластеров (N ≤ 50) корректно идентифицировать диффузионный режим движения частиц весьма затруднительно. Зависимость отношения $\langle {{x}^{2}}(t)\rangle {\text{/}}{{R}^{2}}$ для системы из N = = 150 частиц с различными температурами и параметрами ξ показана на рис. 3. Отметим, что на начальных этапах наблюдения (при $t < \omega _{0}^{{ - 1}}$) режим движения частиц был близок к режиму гармонического осциллятора (1). С ростом времени (при t > t_R) кривые $\langle {{x}^{2}}(t)\rangle {\text{/}}{{R}^{2}}$ стремились к своему постоянному значению примерно равному 2/3.

Рис. 3.

Функция F(t) = 〈x²(t)〉/R² для системы из N = = 150 частиц с параметрами: 1 – ξ ≅ 15, Γ ≈ 3 (Т = = 10 эВ); 2 – ξ ≅ 0.47, Γ ≈ 3 (Т = 10 эВ); 3 – ξ ≅ 15, Γ ≈ 10 (Т = 3 эВ); а также для гармонического осциллятора (1): 4 – ξ ≅ 0.47, Т  = 10 эВ. Штриховая линия – $\langle {{x}^{2}}(t)\rangle {\text{/}}{{R}^{2}} = {\text{ }}2{\text{/}}3$.

Следует отменить, что при любом количестве и температуре частиц поведение функции D(t) для времен наблюдения $t < \omega _{0}^{{ - 1}}$ соответствовало аналитическим решениям для одной частицы (1), см. рис. 4а. При этом выраженный диффузионный режим D(t) ≅ D_b = const их движения наблюдался только в диапазоне времен ${{t}_{R}} \gg t \gg \omega _{0}^{{ - 1}}$, см. рис. 4б.

Рис. 4.

Функция D(t)/D₀ для логарифмической шкалы по времени (а) и обычной (б) при ξ = 0.47 для: 1 – гармонического осциллятора (1); и численные расчеты для N = 450 при разных параметрах Γ: 2 – 30; 3 – 60; 4 – 100; 5 – 150.

Нормированные коэффициенты $D_{b}^{*} = {{D}_{b}}({{\nu }_{{fr}}} + $ $ + \;{{\omega }_{0}})М{\text{/}}T$ в зависимости от параметра Γ, полученные в наших расчетах для ограниченных систем и усредненные для разных значений ξ, приведены на рис. 5. Там же для сравнения приведена зависимость $D_{c}^{*} = {{D}_{с}}({{\nu }_{{fr}}} + \omega {\kern 1pt} *)М{\text{/}}T$ от параметра неидеальности Γ* для протяженных трехмерных структур частиц, взаимодействующих с экранированным кулоновским потенциалом [10]. Легко заметить, что обе зависимости, $D_{b}^{*}(\Gamma )$ и $D_{c}^{*}(\Gamma )$, аналогичны в пределах ошибки численных расчетов по крайней мере до Γ* ≤ 100. Похожее поведение D*(Γ*) было зафиксировано для протяженных квазидвумерных структур частиц, взаимодействующих с экранированным кулоновским потенциалом (систем Юкавы) [5, 9].

Рис. 5.

Величина $D_{c}^{*}$ в зависимости от параметра Γ* для протяженных трехмерных систем Юкавы (сплошная линия) [9, 10]. Символами обозначены результаты наших численных расчетов $D_{b}^{*}$: (⚫) – для ξ < 0.5; (○) – для ξ > 0.5. Показаны абсолютные погрешности расчетов (±0.02) для различных N и ξ.

В заключение подчеркнем, что при любом количестве частиц, N, в ограниченном ансамбле, измерение функции D(t) для времен наблюдения $t < \omega _{0}^{{ - 1}}$ совместно с данными о среднем расстоянии l_p между частицами (например, путем регистрации g(l)) позволяют определить такие параметры пылевых частиц, как их температуру Т, характерные частоты (ν_fr, ω₀), заряд Q, а также градиент α поля ловушки. Поиск этих параметров может опираться на процедуру “подгонки” численных и аналитических данных. (Для проверки и/или коррекции результатов определения параметров можно использовать прямые измерения смещений $\langle {{x}^{2}}(t)\rangle $, и автокорреляционной функции скоростей $\left\langle {{{V}_{x}}\left( 0 \right){{V}_{x}}(t)} \right\rangle $ частиц.)

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

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

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

Рассмотрены условия лабораторных экспериментов необходимые для определения характерных параметров пылевых частиц и поля ловушки.

Представлен анализ зависимости формы парной корреляционной функции g(l) и размеров кластерных систем от температуры и количества частиц. Численное моделирование показало, что для анализа положения первого пика функций g(l) и оценки радиуса неидеальных систем при Г ≥ 1 могут быть использованы аналитические соотношения, полученные в приближении однородной системы.

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

Исследование выполнено за счет гранта Российского научного фонда № 22-22-00899, https://rscf.ru/project/22-22-00899/.