Физика плазмы, 2019, T. 45, № 3, стр. 258-267

1. ВВЕДЕНИЕ

Броуновское движение в системах взаимодействующих частиц широко распространено в природе и наблюдается, например, в биологических и полимерных коллоидных растворах, в плазме продуктов сгорания, в атмосфере Земли и т.д. [1–6]. Особый интерес вызывает влияние теплового движения заряженных частиц на их динамику во внешних магнитных и электрических полях [7–12].

Экспериментальный, теоретический и численный анализ теплового движения взаимодействующих пылевых частиц в протяженных и ограниченных ансамблях, формирующихся в газоразрядной плазме без магнитного поля, представлен в работах [13–19]. Зависимости траекторий и среднеквадратичных смещений одиночной броуновской частицы (N = 1) в электростатической ловушке под воздействием магнитного поля исследовались в теоретических работах [8–10]. Влияние теплового движения ограниченного облака заряженных частиц (N ≤ 500) на их динамику в постоянных электрических и магнитных полях недавно исследовалось численно в [11, 12]. Экспериментальные исследования динамики пылевых частиц в магнитном поле представлены в работах [20–23] для газоразрядной плазмы.

Значительный рост интереса к исследованиям динамики заряженных частиц во внешних электромагнитных полях, наблюдаемый в настоящее время, по большей части связан с проблемами эффективности энергетических установок для управляемого термоядерного синтеза [24–26], а также с развитием плазменных методов переработки отработанного ядерного топлива (ОЯТ) и радиоактивных отходов (РАО) [27–31].

В настоящей работе представлены аналитические и численные исследования спектральной плотности для тепловых смещений центра масс ансамблей заряженных частиц в потенциальном поле электростатической ловушки под воздействием постоянного магнитного поля. В обычных тлеющих разрядах (как переменного, так и постоянного токов) в отсутствие магнитного поля, В = 0, в центре газоразрядных камер наблюдается некоторое превышение концентрации ионов над электронной концентрацией [32]. Данное обстоятельство приводит к формированию эффективных ловушек для отрицательно заряженных частиц (например, для частиц пыли [5, 6]). В плазме с магнитным полем, В ≠ 0, в центре газоразрядной камеры ситуация может быть обратной (за счет “замагничивания” электронов плазмы), т.е. могут существовать условия для удержания положительно заряженных частиц и/или ионов. Так, например, наличие электростатических ловушек для положительно заряженных ионов в установках по разделению компонентов ОЯТ может возникать за счет “замагничивания” электронов на осях камер разрядов зеркального (отражательного) типа, которые обычно используются для данных целей [27–31].

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

2. ОСНОВНЫЕ УРАВНЕНИЯ

Запишем уравнения движения (уравнения Ланжевена) для одной заряженной частицы с массой M и зарядом Q в постоянном электрическом поле ловушки ${\mathbf{E}} = [{{E}_{x}},{{E}_{y}},{{E}_{z}}]$ и в магнитном поле В = В_y (направленном по оси y) под воздействием случайной силы ${{{\mathbf{F}}}_{b}} = [{{F}_{{bx}}},{{F}_{{by}}},{{F}_{{bz}}}]$, которая является источником стохастической (тепловой) энергии частиц

Здесь y, x, z – смещения частицы от ее равновесного положения, ${{V}_{y}} = dy{\text{/}}dt$, ${{V}_{x}} = dx{\text{/}}dt$, ${{V}_{z}} = dz{\text{/}}dt$, ν – коэффициент трения заряженных частиц за счет их столкновений с нейтралами окружающего газа, ω_t =(Qα/M)^1/2 – характерная частота ловушки, ω_В = QВ/M – циклотронная частота вращения в магнитном поле В, α – величина градиента внешнего электрического поля E, а В – величина индукции магнитного поля.

Корни характеристического уравнения (ω_В = 0) для (1) имеют хорошо известный вид

Корни характеристического уравнения для системы (2), (3) можно записать в виде

где ${{\Psi }_{1}} = \nu (1 + {{D}_{1}}{\text{/}}\sqrt 2 ){\text{/}}2$, ${{\Psi }_{2}} = \nu (1 - {{D}_{1}}{\text{/}}\sqrt 2 ){\text{/}}2$, ${{\Omega }_{1}} = $ $ = ({{\omega }_{B}} + \nu {{D}_{2}}{\text{/}}\sqrt 2 ){\text{/}}2$, ${{\Omega }_{2}} = ({{\omega }_{B}} - \nu {{D}_{2}}{\text{/}}\sqrt 2 ){\text{/}}2$,

Средний квадрат отклонений $\left\langle {y{{{(t)}}^{2}}} \right\rangle $ частицы в направлении y от ее начального положения и автокорреляционная функция $\left\langle {y(t*)y({{t}^{*}} + t)} \right\rangle $ могут быть представлены в форме [14, 15, 34, 35]

Здесь и далее T – температура частиц в энергетических единицах, $\psi = {{(4{{\xi }^{2}} - 1)}^{{1/2}}}{\text{/}}2$, $\xi = {{\omega }_{t}}{\text{/}}\nu $, а угловые скобки, $\left\langle {} \right\rangle $, обозначают усреднение по всем отрезкам времени, равным t [1–4, 14, 15, 34, 35]. При этом при t → 0: $\left\langle {y{{{({{t}^{*}})}}^{2}}} \right\rangle = 2T{\text{/}}(М \omega _{t}^{2})$, $\left\langle {y{{{(t)}}^{2}}} \right\rangle = 0$, а при t → ∞: $\left\langle {y{{{({{t}^{*}})}}^{2}}} \right\rangle = 0$ и $\left\langle {y{{{(t)}}^{2}}} \right\rangle = 2T{\text{/}}(М \omega _{t}^{2})$.

Спектральная плотность случайного процесса задается косинусом – преобразованием Фурье для соответствующей автокорреляционной функции [36]. Таким образом, спектральная плотность, G_e(ω), для случайных смещений частицы в ловушке в направлении y (при ω_В = 0), т.е. спектральная плотность классического затухающего осциллятора, может быть записана как [35, 37]

Рассмотрим решения задачи для системы уравнений (2), (3) (при ω_В ≠ 0). Средний квадрат отклонений, $\left\langle {x{{{(t)}}^{2}}} \right\rangle \equiv \left\langle {z{{{(t)}}^{2}}} \right\rangle $, частицы для этой системы в направлениях x и z от ее положения равновесия и автокорреляционная функция, $\left\langle {x({{t}^{*}})x({{t}^{*}}\, + \,t)} \right\rangle \, \equiv \,\left\langle {z({{t}^{*}})z({{t}^{*}}\, + \,t)} \right\rangle $, при условии ${{({{\Psi }_{{1,2}}})}^{2}}\, \ll $ $ \ll {{{\text{(}}{{\Omega }_{{1,2}}})}^{2}}$ (см. Приложение) могут быть представлены как

где ${{\psi }_{1}}\, = \,{{(4\xi _{1}^{2}\, - \,1)}^{{1/2}}}{\text{/}}2$, ${{\xi }_{1}} = {{\omega }_{1}}{\text{/}}{{\nu }_{1}}$, ${{\psi }_{2}} = {{(4\xi _{2}^{2} - 1)}^{{1/2}}}{\text{/}}2$, ${{\xi }_{2}} = {{\omega }_{2}}{\text{/}}{{\nu }_{2}}$, а $\langle x{{({{t}^{*}})}^{2}}\rangle \equiv \langle z{{({{t}^{*}})}^{2}}\rangle = 2T{\text{/}}(М \omega _{t}^{2})$. Здесь ${{\nu }_{1}}\;\; = \;\;{{\Psi }_{1}}$, ${{\nu }_{2}}\;\, = \;\,{{\Psi }_{2}}$, ${{\omega }_{1}}\;\, = \;\,{{(\Psi _{1}^{2}\;\, + \;\,\Omega _{1}^{2})}^{{1/2}}}$, ${{\omega }_{2}}\;\, = $ $ = {{(\Psi _{2}^{2} + \Omega _{2}^{2})}^{{1/2}}}$, ${{B}_{1}} = \sqrt 2 (\omega _{t}^{2} - \omega _{2}^{2}){\text{/}}({{D}_{2}}{{\omega }_{B}}\nu )$, а ${{B}_{1}} = $ $ = \sqrt 2 (\omega _{t}^{2} - \omega _{2}^{2}){\text{/}}({{D}_{2}}{{\omega }_{B}}\nu )$.

При этом спектральная плотность, G_em(ω), для случайных смещений частицы в направлениях x и/или z в постоянном магнитном поле имеет вид

При ω_В → 0 уравнения (10)–(12) переходят в соотношения (7)–(9) соответственно. Иллюстрация зависимости спектральной плотности G(ω) ≡ ≡ G_em(ω) от ω/ν для различных ω_В/ν представлена на рис. 1.

Следует отметить, что для диагностики параметров ограниченных и протяженных систем в лабораторных экспериментах зачастую используется анализ поведения среднего квадрата отклонений частиц от их начального положения (например, функции $\langle x{{(t)}^{2}}\rangle $ в отсутствие магнитного поля, ω_В = 0) [14, 17]. Иллюстрация зависимостей $\langle x{{(t)}^{2}}\rangle $ при различных ξ = ω_t/ν для ω_В = 0 показана на рис. 2. В случае ограниченных систем такие измерения позволяют найти величину среднеквадратичного отклонения $2T{\text{/}}(М \omega _{t}^{2})$ при t → ∞. Откуда при наличии информации о среднем квадрате скорости анализируемых частиц, T/М, легко определить значение характерной частоты ловушки ω_t. Данные о коэффициенте трения частиц ν в этом случае могут быть получены путем анализа максимума функции $\langle x{{(t)}^{2}}\rangle $ (для ξ < 3), см. рис. 2а, или путем “подгонки” численных и экспериментальных данных [14, 17, 19, 20]. Так, например, при $\xi = {{\omega }_{t}}{\text{/}}\nu > 5$ спадание пиков кривых $\langle x{{(t)}^{2}}\rangle $ определяется функцией exp(–νt) и, соответственно, отражает величину ν (см. рис. 2б). Тем не менее тепловые флуктуации в реальных системах могут сильно исказить величину восстанавливаемого коэффициента ν.

Добавим также, что такая процедура хороша только для систем, имеющих одну характерную частоту ω_t. При наличии двух, или более гармоник (характерных частот) в исследуемой системе простой анализ функций $\langle x{{(t)}^{2}}\rangle $ становится затруднительным, см. рис. 3. Тем не менее исследование спектральных характеристик анализируемой системы позволяют легко решить данную задачу (см. рис. 1, а также рисунки ниже для систем с параметрами, аналогичным представленным на рис. 3).

Следует отметить, что уравнения типа (1)–(3) могут использоваться для анализа движения центра масс любого ограниченного ансамбля частиц с попарным взаимодействием, а также для отдельной частицы в системе, состоящей из N частиц, в случае, когда влиянием межчастичного взаимодействия можно пренебречь. Подчеркнем, что для попарных (взаимных) сил межчастичного взаимодействия (F_ij = F_ji, где F_ij – величина силы, действующая со стороны j-й частицы на частицу i; j ≠ i) их сумма в ансамбле, состоящем из N частиц, равна 0. И, соответственно, для движения центра масс системы таких частиц имеют место уравнения типа (1)–(3) [38].

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

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

Численное моделирование динамики систем заряженных частиц в линейной и изотропной электростатической ловушке выполнялось методом молекулярной динамики Ланжевена. Расчеты проводились для одной частицы (N = 1) и для ансамблей, состоящих из N = 50, 250 и 500 частиц. Техника моделирования подробно описана в работах [5, 6]. Шаг интегрирования составлял от $\Delta t \cong {{(40\max [{{\omega }_{t}};\nu ])}^{{ - 1}}}$ до $\Delta t \cong {{(100\max [{{\omega }_{t}};\nu ])}^{{ - 1}}}$ в зависимости от начальных условий задачи. Время расчетов t_c после установления равновесия в моделируемых системах варьировалось от $\sim {\kern 1pt} {{10}^{3}}{\text{/}}\min [{{\omega }_{t}};\quad\nu ]$ до $\sim {{10}^{4}}{\text{/}}\min [{{\omega }_{t}};\quad\nu ]$.

Расчеты проводились для систем частиц с кулоновским взаимодействием в широком диапазоне их параметров неидеальности: от Γ ~ 0.1 до Γ ~ 100. Величина индукции магнитного поля В изменялась в пределах от 300 до $3 \times {{10}^{4}}$ Гс. Значение параметра ξ = ω_t/ν варьировалось от ~1 до ~70, отношение ω_В/ν – от ~0.15 до ~15.

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

Парные корреляционные функции g(l) для ансамблей из N = 500 частиц с различными параметрами Г показаны на рис. 4. В качестве нормировки величины g(l), представленной на этих рисунках, использовалось предположение однородной концентрации частиц равной n ≅ ≅ 3α/(4πQ). Легко увидеть, что первый пик функций g(l) для Г ≥ 0.1 хорошо соответствует величине ${{l}_{p}} \cong {{(4\pi Q{\text{/}}3\alpha )}^{{1/3}}}$, полученной в приближении однородной системы, а для оценки радиуса неидеальных систем может быть использовано соотношениеR ≅ (3N/4πn)^1/3.

Следует подчеркнуть, что перечисленные выше характеристики ($\left\langle {{{{(\Delta x)}}^{2}}} \right\rangle $, g(l), l_p, R) не зависели от величины индукции магнитного поля В.

Численные исследования показали, что вне зависимости от параметра неидеальности Г и числа частиц N в моделируемых ансамблях спектральные плотности смещений их центра масс хорошо соответствуют предлагаемым аналитическим соотношениям (9), (12).

Нормированные спектральные плотности ${{f}^{*}}(\omega ) = NG(\omega ){\text{/}}(2T{\text{/}}M\omega _{t}^{2}\nu )$, полученные для различных параметров задачи (ω_t/ν, ω_В/ν, Г, N), а также аналитические решения (9), (12) показаны на рис. 5–7. (Вычисления спектральной плотности проводились на основе численных расчетов смещений x(t), y(t) и z(t) при помощи процедуры “N-D fast Fourier transform” в пакете прикладных программ MATLAB.)

Следует отметить, что для всех представленных случаев величина ${{({{\Psi }_{{1,2}}})}^{2}} \ll {{({{\Omega }_{{1,2}}})}^{2}}$, за исключением случая ω_В/ν ≅ 1.5, ω_t/ν ≅ 1 (см. рис. 6б), когда ${{({{\Psi }_{{1,2}}})}^{2}} \cong 0.25{{({{\Omega }_{{1,2}}})}^{2}}$. Тем не менее и для этого случая получены хорошие совпадения между численными и аналитическими результатами.

Численные исследования также показали, что отклонения формы спектральных распределений для отдельных частиц ансамбля от аппроксимирующих функций (9), (12) наблюдаются при Γ > < 0.1. С ростом величины Γ характерная частота спектра смещается в сторону меньших частот относительно частоты колебаний центра масс системы как для смещений частиц вдоль магнитного поля (ω_В = 0), так и для смещений в плоскости ортогональной магнитному полю (ω_В/ν ≠ 0). Иллюстрация поведения нормированных спектральных плотностей $f{\text{*}}(\omega {\text{/}}\nu )$ для отдельных частиц ансамбля с ростом Γ для случая ω_В = 0 представлена на рис. 8. Следует отметить, что простая подгонка численных данных, полученных для отдельных частиц таких систем (при Γ > 0.1), аналитическими соотношениями (9), (12) не позволяет получить физически обоснованных результатов. Тем не менее можно предположить, что данное обстоятельство связано с уменьшением коэффициента диффузии частиц с ростом параметра Γ, которое, в свою очередь, происходит за счет роста эффективной диссипации в сильно коррелируемых системах [14, 15, 35].

В заключение данного раздела следует отметить, что измерения спектральных плотностей смещений для центра масс ограниченной системы частиц позволяет восстанавливать характерные параметры исследуемых ансамблей (ν, ω_t, ω_В) даже при наличии значительных тепловых флуктуаций в анализируемой системе.

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

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

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

Кроме того, решение уравнений движения (1)–(3) при заданной частоте трения ν зависит только от относительных значений параметров ω_B/ν, ω_t/ν. Таким образом, полученные результаты справедливы для частиц любой массы и зарядов (например, для случая пылевой плазмы, для ионов сепарируемого вещества ОЯТ и т.д.).

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

Данная работа была поддержана Российским научным фондом (грант № 14-29-00231).

ПРИЛОЖЕНИЕ

Любое решение F(t) для задачи (2), (3) можно представить в виде суперпозиции

При этом, для функции $\left\langle {x{{{(t)}}^{2}}} \right\rangle $ величина ${{С }_{0}} = {\text{ }}2T{\text{/}}(М \omega _{t}^{2})$, для $\left\langle {x({{t}^{*}})x({{t}^{*}} + t)} \right\rangle --{{С }_{0}} = 0$. Для поиска коэффициентов С_i ($i = 1,2,3,4$) указанных функций используются начальные условия задачи: $F\left( 0 \right) = {{Т }_{0}}$, $dF\left( 0 \right){\text{/}}dt = 0$, ${{d}^{2}}F\left( 0 \right){\text{/}}d{{t}^{2}} = {{Т }_{2}}$, ${{d}^{3}}F\left( 0 \right){\text{/}}d{{t}^{3}} = 0$, где Т₀ = 0 для $\left\langle {x{{{(t)}}^{2}}} \right\rangle $, ${{T}_{0}} = $ $ = {\text{ }}2T{\text{/}}(М \omega _{t}^{2})$ для $\left\langle {x({{t}^{*}})x({{t}^{*}} + t)} \right\rangle $, Т₂ = 2T/М для $\left\langle {x{{{(t)}}^{2}}} \right\rangle $ и Т₂ = –2T/М для $\left\langle {x({{t}^{*}})x({{t}^{*}} + t)} \right\rangle $. Тогда для коэффициентов С_i можно записать следующую систему уравнений:

Здесь и далее знак “–” для функции $\left\langle {x{{{(t)}}^{2}}} \right\rangle $, знак “+” для $\left\langle {x({{t}^{*}})x({{t}^{*}} + t)} \right\rangle $,

Используя значения λ_i (5а), (5б), в предположении ${{D}_{1}}^{2} \ll {{D}_{2}}{{\omega }_{В }}{\text{/}}(\nu \sqrt 2 )$ и ${{\nu }^{2}} - {{D}_{1}}^{2}{\text{/}}2 \ll 2{{\omega }_{{1,2}}}$ (где ${{\omega }_{1}} = {{(\Psi _{1}^{2} + \Omega _{1}^{2})}^{{1/2}}}$, а ${{\omega }_{2}} = {{(\Psi _{2}^{2} + \Omega _{2}^{2})}^{{1/2}}}$) формулы (6)–(9) можно представить в виде

Здесь ${{B}_{1}} = \sqrt 2 (\omega _{t}^{2} - \omega _{2}^{2}){\text{/}}({{D}_{2}}{{\omega }_{В }}\nu )$, ${{B}_{2}} = \sqrt 2 (\omega _{1}^{2} - \omega _{t}^{2}){\text{/}}$ ${\text{/}}({{D}_{2}}{{\omega }_{В }}\nu )$.

Отметим, что упомянутые выше условия ($D_{1}^{2} \ll {{D}_{2}}{{\omega }_{В }}{\text{/}}(\nu \sqrt 2 )$ и ${{\nu }^{2}} - D_{1}^{2}{\text{/}}2 \ll 2{{\omega }_{{1,2}}}$) выполняются, когда ${{({{\Psi }_{{1,2}}})}^{2}} \ll {{({{\Omega }_{{1,2}}})}^{2}}$.