Физика плазмы, 2022, T. 48, № 7, стр. 613-627

1. ВВЕДЕНИЕ

Последовательное изучение резонансного взаимодействия волн и частиц в плазме началось с классической работы Ландау [1], в которой было открыто затухание ленгмюровской волны в плазме с максвелловским распределением электронов по скоростям, обусловленное их резонансным взаимодействием с волной. Подход, использованный Ландау, основан на линеаризации кинетического уравнения с самосогласованным полем, при котором в членах, пропорциональных полю волны, функция распределения заменяется на ее невозмущенное значение. Этот подход первоначально был предложен Власовым [2], так что кинетическое уравнение с самосогласованным полем называют уравнением Больцмана–Власова.

Подход, предложенный Ландау для анализа затухания ленгмюровских волн, был затем обобщен на волновые моды в магнитоактивной плазме, в частности, на свистовые волны [3]. Кроме того, важнейшим результатом исследования резонансного взаимодействия волн и частиц в плазме стало установление факта, что в случае неравновесного распределения заряженных частиц плазмы по скоростям это взаимодействие может приводить к усилению волн. С этого началось изучение линейных кинетических неустойчивостей в плазме, исследованию которых посвящено большое число оригинальных работ. Ссылки на многие из них могут быть найдены в двухтомной монографии Михайловского [4], посвященной исследованию плазменных неустойчивостей в однородной и неоднородной плазме, которая сама по себе является важным вкладом в теорию плазменных неустойчивостей.

Линеаризация кинетического уравнения, указанная выше, эквивалентна использованию невозмущенных траекторий в качестве характеристик кинетического уравнения. Такой подход имеет, очевидно, ограниченную область применимости и справедлив только на временах, меньших нелинейного времени τ [5], на котором возмущение траекторий резонансных частиц под действием поля волны становится существенным. Для ленгмюровской волны нелинейное время $\tau = (m{\text{/}}e{{E}_{0}}k{{)}^{{1/2}}},$ где m и $ - e$ – масса и заряд электрона, ${{E}_{0}}$ – амплитуда ленгмюровской волны и k – ее волновой вектор. Для свистовой волны, распространяющейся вдоль внешнего магнитного поля, нелинейное время, наряду с параметрами волны, зависит также от поперечной скорости частицы ${{v}_{ \bot }}$, определяясь выражением ${{(m\omega {\text{/}}eE{{k}^{2}}{{v}_{ \bot }})}^{{1/2}}},$ где теперь E – амплитуда электрического поля свистовой волны, и ω – ее частота [6]. При выполнении условия $\gamma \tau \gg 1$, где γ – декремент волны, последняя успевает затухнуть за время $ \lesssim {\kern 1pt} \tau $, так что неравенство $\gamma \tau \gg 1$ определяет область применимости линейной теории резонансного взаимодействия частиц с монохроматической волной в однородной плазме. В противоположном случае $\gamma \tau \ll 1$, детально исследовавшемся, начиная с работы [5] (см. монографию [7], где также имеются многочисленные ссылки на другие оригинальные работы в этом направлении), изменение амплитуды волны в течение нелинейного времени мало, так что при решении уравнений движения частиц амплитуду волны можно считать постоянной. Это соответствует другому способу линеаризации кинетического уравнения, при котором поле в нелинейном члене кинетического уравнения считается заданным. При таком подходе изменение амплитуды волны возникает в следующем приближении по малому параметру $\gamma \tau \ll 1$ и находится из закона сохранения энергии в системе волна–резонансные частицы. Как показано в [5], в случае монохроматической волны и при выполнении условия $\gamma \tau \ll 1$ декремент волны выходит на ноль и, соответственно, ее амплитуда – на постоянное значение на временах, превышающих несколько нелинейных периодов τ. Это обусловлено фазовым перемешиванием частиц в резонансной области фазового пространства и эргодизацией их функции распределения [5].

Метод решения комплексного дисперсионного уравнения и определения декремента, развитый в работе [1], применим только при условии, что декремент волны много меньше ее частоты. Решение комплексного дисперсионного уравнения, которое следует из линеаризованной системы Больцмана–Власова–Максвелла, в случае, когда неравенство $\gamma \ll \omega $ не выполнено, исследовалось в ряде работ (см., например, [8]) в предположении, что невозмущенное распределение частиц описывается достаточно простой аналитической функцией, в частности, максвелловской или би-максвелловской. Для более сложных аналитических функций были предложены численные методы решения дисперсионного уравнения (см., например, [9]).

В случае, когда инкремент волны много меньше ее частоты, а зависимость частоты от волнового вектора известна, инкремент волны может быть получен из закона сохранения энергии в системе волна–резонансные частицы, с тем же результатом, который следует из решения комплексного дисперсионного уравнения [10]. При этом вычисление инкремента сводится к интегрированию по поперечной скорости некоторого выражения, зависящего от производных функции распределения, взятых при резонансной продольной скорости. Другим используемым методом для получения инкремента при условии $\gamma \ll \omega $ является вычисление анизотропии функции распределения, через которую инкремент выражается известной формулой из работы [11]. Такой метод использовался, например, в работе [12], в то время, как упомянутый перед этим метод применялся в работах [13, 14].

Задача, рассматриваемая в настоящей работе, связана с измерениями волн и частиц на спутниках, в частности, Van Allen Probs [15] и MMS [16]. При анализе этих данных возникает необходимость вычисления инкремента волн в случае, когда функция распределения энергичных частиц определяется их дифференциальными потоками и имеет только численное представление. Подчеркнем, что неаналитическое задание функции распределения и нарушение условия $\gamma \ll \omega $ – это два различных и независимых момента, которые могут присутствовать в задаче о вычислении инкремента волны. (Далее в этой работе мы будем говорить о свистовых волнах.) Например, в случае наблюдения волн во внутренней магнитосфере на спутниках Van Allen Probes условие $\gamma \ll \omega $, как правило, выполняется, в то время как в случае волн в хвосте магнитосферы, наблюдаемых на спутниках MMS, это условие может нарушаться.

Следует отметить, что численным методам решения дисперсионного уравнения для волн в плазме по-прежнему уделяется большое внимание. Помимо упомянутых выше работ укажем также две недавние работы, посвященные развитию математических и численных подходов для решения дисперсионного уравнения в “горячей” магнитоактивной плазме [17, 18]. Если в работе [17] для решения дисперсионного уравнения при численном задании функции распределения используется интерполяция кубическими сплайнами, то в работе [18] используется новый метод гибридного аналитического продолжения. Аналитическое продолжение функции, определяющей дисперсионное уравнение, является центральным моментом и в нашей работе. Используемый нами метод аналитического продолжения основан на принципе непрерывности.

Настоящая работа исследует частный случай задачи, рассмотренной в [18], а именно, в ней развивается метод решения дисперсионного уравнения для свистовых волн, распространяющихся вдоль внешнего однородного магнитного поля, без предположения о малости инкремента, а функция распределения задается численно на основе измеряемых дифференциальных потоков частиц. Мы показываем, что в этом важном для приложений частном случае задачу удается свести к одномерной в пространстве скоростей, так что дисперсионное уравнение по форме совпадает с уравнением, полученным в работе [1], с важным отличием в том, что функция распределения задана численно. Естественно, что при невыполнении условия $\gamma \ll \omega $ требуется другой метод решения самого дисперсионного уравнения. Для того, чтобы сделать переход к одномерной задаче более наглядным, мы выводим дисперсионное уравнение из первых принципов, вместо того, чтобы упрощать достаточно громоздкие общие формулы для дисперсионного уравнения для волн в магнитоактивной плазме.

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

Будем считать, что внешнее магнитное поле ${{{\mathbf{B}}}_{0}}$ направлено вдоль оси $z$ декартовой системы координат. Тогда электромагнитное поле свистовой волны, распространяющейся вдоль внешнего магнитного поля, зависит от одной координаты $z$. Такая волна, как известно [19], является чисто поперечной и имеет круговую поляризацию.

Запишем действительное электрическое и магнитное поле волны в виде

где k – продольная компонента волнового вектора ${\mathbf{k}} = (0,0,k)$. Из закона индукции Фарадея тогда имеем

В случае, когда зависимость всех величин от координаты z дается множителем ${{e}^{{ikz}}}$, знак действительной части в уравнении (2) можно снять и получить следующие соотношения между комплексными компонентами электрического и магнитного поля волны:

Точно такие же соотношения, очевидно, имеют место и между компонентами ${\mathbf{E}}(t)$ и $\partial {\mathbf{B}}(t){\text{/}}\partial t$. Обратим внимание, что без конкретизации зависимости поля от времени мы не получаем соотношения между E и B, а только между E и $\partial {\mathbf{B}}{\text{/}}\partial t$.

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

Соотношения (3) и (4) позволяют выразить все компоненты электромагнитного поля волны через одну компоненту магнитного поля и ее производную, например, через ${{B}_{x}}$ и $\partial {{B}_{x}}{\text{/}}\partial t$.

Для вывода уравнения для поля волны будем исходить из уравнения Максвелла

С учетом (1) это уравнение можно записать в виде системы двух уравнений

Умножим теперь первое уравнение на i и сложим со вторым. Выражая все компоненты поля через ${{\widetilde B}_{x}}$ и $\partial {{\widetilde B}_{x}}{\text{/}}\partial t$ и учитывая, что для любой комплексной величины Z имеет место равенство ${\text{Re}}{\kern 1pt} ( - iZ) = {\text{Im}}{\kern 1pt} (Z)$, получим одно уравнение для комплексной величины ${{\widetilde B}_{x}}$

где учтено, что согласно (1) $\partial {{\widetilde B}_{x}}{\text{/}}\partial z = ik{{\widetilde B}_{x}}$.

Перейдем теперь к уравнению для функции распределения электронов, которая определяет плотность тока электронов, входящую в (7). Ограничимся линейным приближением по полю волны и представим функцию распределения электронов $f({\mathbf{r}},{\mathbf{v}},t)$ в виде

Поскольку поле волны зависит от одной координаты z, то же верно и для возмущения функции распределения, так что линеаризованное кинетическое уравнение имеет вид

В дальнейшем будем полагать, что невозмущенная функция распределения ${{f}_{0}}$ зависит только от продольной скорости ${{v}_{\parallel }}$ и модуля перпендикулярной скорости ${{v}_{ \bot }}$.

Перейдем от переменных ${{v}_{x}},\;{{v}_{y}},\;{{v}_{z}}$ к переменным ${{v}_{ \bot }},\;\varphi ,\;{{v}_{\parallel }}$ согласно соотношениям:

В новых переменных с учетом соотношений (1)–(4) уравнение (8) принимает вид

Уравнение (10) имеет структуру

где $\hat {L}$ – линейный дифференциальный оператор с действительными коэффициентами. Поэтому, если мы найдем решение уравнения то ${\kern 1pt} {\text{Re}}{\kern 1pt} \{ \widetilde \Phi \} $ будет, очевидно, удовлетворять уравнению (10). Поскольку коэффициенты оператора $\hat {L}$ не зависят ни от z, ни от φ, то величина $\widetilde \Phi $ должна иметь вид

Подставляя (12) в (10), проводя дифференцирование в левой части уравнения по z и φ и сокращая затем обе части уравнения на $\exp (ikz + i\varphi )$, получим уравнение для (комплексной) величины Φ:

Линейное возмущение функции распределения, то есть действительная величина $\delta f$, выражается через Φ соотношением

где * означает комплексное сопряжение. Нетрудно видеть, что комплексный ток ${{j}_{x}} - i{{j}_{y}}$, входящий в уравнение (7) для ${{\widetilde B}_{x}}(t,z) = B(t){{e}^{{ikz}}}$, выражается через функцию Φ соотношением так что с учетом (7) и (15) уравнение для (комплексного) поля ${{B}_{x}}(t)$ принимает вид

Введем теперь функцию $F(t,{{v}_{\parallel }})$ согласно соотношению

Умножая уравнение (13) на $\pi v_{ \bot }^{2}$ и интегрируя по $d{{v}_{ \bot }}$, получим уравнение для функции $F(t,{{v}_{\parallel }})$ в виде

где а уравнение (16) для поля ${{B}_{x}}(t)$ с учетом (17) принимает вид

Уравнения (18), (20) образуют замкнутую систему уравнений для определения электромагнитного поля свистовой волны, распространяющейся вдоль внешнего магнитного поля в однородной магнитоактивной плазме, и возмущения функции распределения в линейном приближении. Входящее в эту систему поле ${{B}_{x}}$ зависит только от времени t, а функция распределения F – только от t и ${{v}_{\parallel }}$, так что эта система вполне аналогична системе Больцмана–Власова для ленгмюровской волны [2].

Следуя работе [1], будем искать решение этой системы с помощью преобразования Фурье–Лапласа с комплексной частотой ω. Для поля ${{B}_{x}}(t)$, удовлетворяющего условию ${{B}_{x}}(t) = 0$ при $t < 0$ и ${\text{|}}{{B}_{x}}(t){\text{|}} < M{{e}^{{\alpha t}}}$ при $t \geqslant 0$, где M и α – две положительные константы, это преобразование записывается в виде [20]

Интегрирование во втором выражении (21) производится в комплексной плоскости ω по прямой, параллельной действительной оси, лежащей над ней на расстоянии $\sigma > \alpha $. Аналогичное преобразование по переменной t имеет место и для возмущения функции распределения $F(t,{{v}_{\parallel }})$. (Поскольку мы обозначаем прямое и обратное преобразование одной и той же буквой, мы везде указываем аргументы соответствующих функций).

Проводя стандартные вычисления, из (18), (20) получим

Выражая величину $F(\omega ,{{v}_{\parallel }})$ из (23) и подставляя ее в (22), получим замкнутое уравнение для величины ${{B}_{x}}(\omega )$

3. ДИСПЕРСИОННОЕ УРАВНЕНИЕ

Если величина ${{B}_{x}}(\omega )$ известна, то зависимость поля от времени определяется вторым соотношением в (21), где интегрирование ведется в верхней полуплоскости ω по прямой, параллельной действительной оси. На этой прямой и выше нее Лаплас-образ поля ${{B}_{x}}(\omega )$ является аналитической функцией ω. Однако для практического вычисления этого интеграла контур интегрирования следует замкнуть в нижней полуплоскости ω. Как следует из (24), величина ${{B}_{x}}(\omega )$ равна отношению правой части уравнения (24) и величины

равенство нулю которой представляет собой дисперсионное уравнение для комплексной частоты ω как функции волнового числа k при заданной функции распределения ${{f}_{{0\parallel }}}$ и определяемой ею функции $dg{\text{/}}d{{v}_{\parallel }}$. Выражение для величины $D(\omega ;k)$ может быть переписано в виде где мы учли, что интеграл $\int_{ - \infty }^\infty {{{f}_{{0\parallel }}}d{{{v}}_{\parallel }}} $ равен невозмущенной плотности плазмы ${{n}_{0}}$ и ввели электронную плазменную частоту согласно соотношению $\omega _{p}^{2} = 4\pi {{n}_{0}}{{e}^{2}}{\text{/}}m$.

Уравнение (26) определяет величину $D(\omega ;k)$ в верхней полуплоскости ${\kern 1pt} {\text{Im}}{\kern 1pt} (\omega ) > \sigma > 0$, где $D(\omega ;k)$ является аналитической функцией. Однако корень дисперсионного уравнения

лежит ниже прямой ${\kern 1pt} {\text{Im}}{\kern 1pt} (\omega ) = \sigma $, и для его нахождения функцию $D(\omega ;k)$ необходимо аналитически продолжить в нижнюю полуплоскость $\omega $. Прежде, чем перейти к этому вопросу, заметим следующее. Использованный нами подход, который заключается во введении вместо невозмущенной функции распределения ${{f}_{0}}({{v}_{ \bot }},{{v}_{\parallel }})$, зависящей от двух переменных ${{v}_{ \bot }},\;{{v}_{\parallel }}$, двух функций $g({{v}_{\parallel }})$ и ${{f}_{{0\parallel }}}({{v}_{\parallel }})$ (19), зависящих от одной переменной ${{v}_{\parallel }}$, сводит отыскание комплексной частоты свистовой волны к решению уравнения $D(\omega ;k) = 0$, которое аналогично дисперсионному уравнению для ленгмюровских волн, подход к решению которого был развит в работе [1]. Для аналитического продолжения интеграла по ${{v}_{\parallel }}$, определяющего функцию $D(\omega ;k)$, Ландау предложил сместить контур интегрирования по ${{v}_{\parallel }}$ в комплексную плоскость ${{v}_{\parallel }}$, считая производную функции распределения, стоящую в числителе под знаком интеграла, аналитической функцией ${{v}_{\parallel }}$. В нашем случае использование этого метода не очевидно, поскольку функция распределения задана численно и только при действительных значениях скорости. Для аналитического продолжения интеграла в нижнюю полуплоскость мы также будем считать стоящие в числителе под знаком интеграла (26) функции аналитическими, однако для аналитического продолжения интеграла в нижнюю полуплоскость ω мы используем принцип непрерывности [21]. (Отметим, что в случае неустойчивой плазмы, корень дисперсионного уравнения лежит в верхней полуплоскости ω, и для его нахождения аналитического продолжения интеграла не требуется.) Применительно к рассматриваемой нами задаче, принцип непрерывности можно сформулировать следующим образом. Пусть ${{\Psi }_{1}}(\omega )$ – аналитическая функция ω в верхней полуплоскости и непрерывная вплоть до действительной оси, и пусть ${{\Psi }_{2}}(\omega )$ – аналитическая функция ω в нижней полуплоскости, также непрерывная вплоть до действительной оси. Тогда, если ${{\Psi }_{1}}(\omega ) = {{\Psi }_{2}}(\omega )$ на действительной оси, то функция ${{\Psi }_{2}}(\omega )$ является аналитическим продолжением функции ${{\Psi }_{1}}(\omega )$ в нижнюю полуплоскость.

Обозначим числитель выражения под интегралом (26) через $Q({{v}_{\parallel }})$:

Для определенности будем считать $k > 0$. Очевидно, что интеграл

является аналитической функцией ω в верхней полуплоскости. Его значение на действительной оси задается пределом при ${\kern 1pt} {\text{Im}}{\kern 1pt} (\omega ) \to + 0$, который равен где мы представили комплексную частоту ω в виде $\omega = a + ib$ и использовали формулу Сохоцкого (см., например, [22].) Символ P перед интегралом означает главное значение соответствующего интеграла. Интеграл (29) является аналитической функцией также и в нижней полуплоскости $\omega $, однако его предел при ${\kern 1pt} {\text{Im}}{\kern 1pt} (\omega ) \to - 0$ равен

Считая $Q({{{v}}_{\parallel }})$ аналитической функцией ${{{v}}_{\parallel }}$, рассмотрим тогда функцию ${{I}_{2}}(\omega )$, заданную в нижней полуплоскости ω соотношением

Очевидно, что

т.е. на действительной оси ${{I}_{2}}(\omega ) = {{I}_{1}}(\omega ).$ Из принципа непрерывности тогда следует, что функция ${{I}_{2}}(\omega )$, определяемая соотношением (32), является аналитическим продолжением функции ${{I}_{1}}(\omega )$ в нижнюю полуплоскость.

Напишем теперь выражение для интеграла (29) в верхней полуплоскости ω, справедливое вплоть до действительной оси и пригодное для численных расчетов, когда функция $Q({{v}_{\parallel }})$ задана численно. Используя, как и выше, запись комплексной частоты в виде $\omega = a + ib$, получим

где введено обозначение

Разобьем теперь интегралы в правой части (33) на три интервала

считая величину A достаточно большой, так что функция $Q({{v}_{\parallel }})$ при ${{v}_{\parallel }} < {{v}_{R}} - A$ и ${{v}_{\parallel }} > {{v}_{R}} + A$ обращается в ноль, а величину $\delta v$ достаточно малой, так что изменением функции $Q({{v}_{\parallel }})$ на интервале $[{{v}_{R}} - \delta v,{{v}_{R}} + \delta v]$ можно пренебречь. Тогда значение первого интеграла в правой части (33) по промежутку $[{{v}_{R}} - \delta v,{{v}_{R}} + \delta v]$ обращается в ноль как интеграл от нечетной функции по симметричному промежутку, а второй интеграл по этому промежутку вычисляется и дает

В результате получим

Приведем выражение для величины $D(\omega ;k)$ (26), определяющей дисперсионное уравнение (27), в важном для магнитосферных приложений случае, когда функцию распределения электронов можно представить в виде суммы “холодной” и “горячей” компонент. Будем отмечать величины, относящиеся к этим компонентам, верхними индексами $^{{(c)}}$ и $^{{(h)}}$ соответственно. Запишем

Подчеркнем, что плотность горячей компоненты, равная $\int f_{0}^{{(h)}}({{{v}}_{ \bot }},{{{v}}_{\parallel }}){{d}^{3}}{\mathbf{v}}$, не предполагается малой по сравнению с ${{n}^{{(c)}}}$. При этом величина $\omega _{p}^{2}$, входящая в (26), пропорциональна сумме ${{n}^{{(c)}}} + {{n}^{{(h)}}}$. Из определений (19) следует, что

так что величина ${{Q}^{{(c)}}}({{{v}}_{\parallel }})$ (28), связанная с холодной компонентой распределения, равна

Тогда, после учета вклада холодной компоненты в (26), дисперсионное уравнение принимает вид

Таким образом, обсуждавшееся выше вычисление интеграла, входящего в дисперсионное уравнение, фактически необходимо лишь для учета горячей компоненты электронов плазмы. В пренебрежении горячей компонентой, ${{\omega }_{p}} = \omega _{p}^{{(c)}}$, и уравнение (40) дает известное выражение для квадрата показателя преломления свистовых волн, распространяющихся вдоль внешнего магнитного поля

Выражение (36) для интеграла, входящего в дисперсионное уравнение (26) и (40), справедливо только в верхней полуплоскости. Однако при решении дисперсионного уравнения область, где лежит его корень, заранее не известна. Поэтому необходимо иметь общее выражение для этого интеграла и, следовательно, для дисперсионного уравнения, справедливое во всей комплексной плоскости ω. Используя полученное выше аналитическое продолжение интеграла, входящего в дисперсионное уравнение, а именно, соотношение (32), запишем выражение для этого интеграла в виде, справедливом во всей комплексной плоскости и удобном для его численного вычисления:

Нетрудно видеть, что при ${\text{|}}b{\text{|}} \to 0$ представление (42) переходит в выражение, которое дается формулой Сохоцкого в предположении $b \to + 0$ и соответствует правилу обхода полюса Ландау.

Прежде чем перейти к результатам расчета инкремента, сделаем два важных замечания, относящихся к численным расчетам. Интервал интегрирования в правой части (42), вне которого по определению функция $Q({{v}_{\parallel }})$ обращается в ноль, симметричен относительно ${{v}_{R}} < 0$. В то же время, реально задаваемая численно функция $Q({{v}_{\parallel }})$ обращается в ноль вне меньшего интервала, который является частью интервала $({{v}_{R}} - A,{{v}_{R}} + A)$ и симметричен относительно некоторой близкой к нулю точки ${{v}_{0}}$. Тем не менее, разные по длине реальные промежутки интегрирования левее ${{v}_{R}} - \delta v$ и правее ${{v}_{R}} + \delta v$ должны быть разбиты на маленькие интервалы одинаковой длины, но не на равное число интервалов. Другими словами, численное интегрирование должно обеспечивать точное равенство нулю интеграла

по любому симметричному относительно ${{{v}}_{R}}$ интервалу, т.е. при произвольных ${{A}_{1}}$, $\delta {{{v}}_{1}}$, для произвольного значения b. Что касается величины $\delta {v}$, то, как указано выше, она должна быть настолько мала, чтобы изменением величины $Q({{{v}}_{\parallel }})$ при изменении ${{{v}}_{\parallel }}$ на $\delta {v}$ можно было пренебречь. Так, если функция $Q({{{v}}_{\parallel }})$ задана с шагом $\Delta {v}$, то величина $\delta {v}$ должна быть в несколько раз меньше. Естественно, что результат вычислений не должен зависеть от величины $\delta {v}$. В наших расчетах мы использовали $\delta {v} = \Delta {v}{\text{/}}2.$

Изложенный выше метод решения дисперсионного уравнения в области отрицательных значений γ использует аналитическое представление функции $Q({{{v}}_{\parallel }})$ в нижней полуплоскости. Приведем теперь метод решения дисперсионного уравнения, свободный от данного ограничения и применимый в случае численного задания функции распределения, в частности, на основе измерений дифференциальных потоков частиц. Для этого заметим, что изначально нам необходимо иметь лишь аналитическое представление самого интеграла (42) в комплексной плоскости $\omega $, а не функции $Q({{v}_{\parallel }})$, причем интеграл (42) является аналитической функцией $\omega $ во всей верхней полуплоскости вплоть до действительной оси, независимо от аналитических свойств функции $Q({{{v}}_{\parallel }})$. Поэтому мы можем вначале найти интеграл (42) на действительной оси, который равен

а затем продолжить этот интеграл в комплексную плоскость, что сводится к замене в выражении для $I(a;k)$ действительной величины a на комплексную величину ω. Заметим, что сама величина $I(a;k)$ является комплексной функцией действительной переменной a и действительного параметра k. Cамо же представление интеграла $I(a;k)$ в виде аналитической функции легко может быть получено после его численного вычисления, например, с помощью функции polyfit в программе Matlab. Естественно, что в области положительных γ оба способа решения дисперсионного уравнения дают одинаковые результаты. Мы использовали оба способа решения дисперсионного уравнения для вычисления инкрементов на основе спутниковых измерений дифференциальных потоков частиц. Результаты этих вычислений приводятся ниже.

4. РАСЧЕТ ИНКРЕМЕНТА СВИСТОВЫХ ВОЛН ПО ДАННЫМ СПУТНИКОВ VAN ALLEN PROBE-A И MMS

Применим развитый выше подход для вычисления инкремента свистовых волн на основе измерений плотности холодной плазмы и распределения энергичных частиц, выполненных на спутнике Van Allen Probe-A 12 декабря 2017 г. в 11:25:34 UT. В это время спутник находился на L‑оболочке $L = 3.21$ вблизи экватора (магнитная широта составляла –0.28°). Плотность холодной плазмы равнялась 1383 см^–3. Надежные измерения потоков энергичных электронов были в пяти энергетических каналах 33, 54, 80, 108 и 143 кэВ, на девяти питч-углах 24.5, 40.9, 57.3, 73.6, 90.0, 106.4, 122.7, 139.1 и 155.5 градуса. Эти измерения проводились прибором MagEIS (The Magnetic Electron Ion Spectrometer) [23], который выполнял измерения потоков энергичных частиц с временным разрешением ~10 с. Измеренные потоки позволяют получить функцию распределения энергичных электронов, используя соотношение [13]

где ${{f}_{{{\kern 1pt} {\text{СГС}}}}}$ – функция распределения в системе единиц СГС, а поток частиц J и энергия W в правой части соотношения (44) выражены в практических единицах 1/(см² ⋅ с ⋅ ср ⋅ кэВ) и кэВ соответственно. Переменные, от которых зависит функция распределения, никак не влияют на соотношение (44): функция распределения всегда считается нормированной так, что $\int f{{d}^{3}}{\mathbf{v}} = n$, где $d{\mathbf{v}}$ – элемент объема в пространстве скоростей, а $n$ – плотность частиц, и все величины выражены в системе СГС. При этом сама функция распределения может быть выражена в любых переменных. В частности, для использования приведенных выше соотношений, нам необходима функция распределения в переменных $({{v}_{\parallel }},{{v}_{ \bot }})$, для получения которой в соотношении (44) необходимо выразить энергию и питч-угол, от которых зависит поток J, через переменные $({{v}_{\parallel }},{{v}_{ \bot }})$. Зная функцию распределения в переменных $({{v}_{\parallel }},{{v}_{ \bot }})$, мы можем рассчитать функции ${{f}_{{0\parallel }}}({{v}_{\parallel }})$ и $g({{v}_{\parallel }})$ (19), входящие в дисперсионное уравнение. Эти функции для указанных выше момента времени и положения спутника Van Allen Probe-A показаны на рис. 1.

В обсуждаемых данных электронная плотность ${{N}_{e}}$ измеряется отдельно, независимо от потоков энергичных электронов, и намного превышает плотность энергичных частиц. В этом случае частоту и инкремент волны следует определять из соотношения (40), где ${{\omega }_{p}}$ – плазменная частота с учетом горячей компоненты плазмы, причем, как указано выше, $\omega _{p}^{2} \simeq {{(\omega _{p}^{{(c)}})}^{2}}$.

Соотношение (40) определяет действительную и мнимую часть частоты как функцию волнового числа k. Однако в рассматриваемом случае, когда дисперсия волн определяется, в основном, холодной плазмой, действительная часть частоты и волновое число являются взаимно однозначными функциями, так что инкремент волны может быть выражен как функция частоты. Эта зависимость, полученная из численного решения уравнения (40), показана на рис. 2 сплошной линией.

В случае, когда дисперсия волны определяется плотной холодной компонентой плазмы, а ее инкремент – энергичными частицами малой плотности, которые не влияют на дисперсионные свойства волн, ее линейный инкремент определяется выражением (см., например, [24])

Выражение (45) предполагает, что показатель преломления волны $kc{\text{/}}\omega $ является известной величиной, зависящей от частоты волны, а также определенное выражение для плотности энергии волны, также определяемое дисперсионным соотношением. Зависимость инкремента γ от частоты волны ω, рассчитанная по формуле (45) с использованием той же самой измеренной функции распределения, что и выше, показана на рис. 2 штриховой линией.

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

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

Мы приведем результаты решения дисперсионного уравнения (26) на основе данных по потокам частиц, полученным на одном из спутников MMS [16] 19 июня 2017 г., в 03:57:38.625, когда спутник находился в хвосте магнитосферы на расстоянии 17 радиусов Земли. Согласно спутниковым измерениям электронная циклотронная частота ${{f}_{c}}$ в области наблюдения составляла ~0.515 кГц. Плотность частиц, рассчитанная по измеренной функции распределения горячих электронов, составляла 0.0887 см^–3, что соответствует плазменной частоте 2.6736 кГц Отметим, что существенное превышение плазменной частоты электронной гирочастоты вовсе не подразумевает возможность использования дисперсионного соотношения для свистов (41), справедливого в холодной плазме.

Измерения содержат надежные данные по дифференциальным потокам частиц в диапазоне энергий от 96 эВ до 21 кэВ, равномерно расположенных на логарифмической шкале (21 канал), в диапазоне питч-углов от 5° до 175° с шагом 10°. Измерения дифференциальных потоков частиц в так называемой “burst” моде, которые мы используем, проводились с временным разрешением ~30 мс. Такой режим измерений включался на короткое время при пересечении спутником плазменных границ и неоднородностей, таких как магнитопауза, диполизационные фронты, быстрые плазменные потоки и др. Что касается данных по плотности и, соответственно, плазменной частоте, то эти данные точно соответствуют приведенным потокам частиц и, по-видимому, не учитывают электроны с энергиями ниже 96 эВ. Поэтому плотность холодных частиц остается неопределенным параметром. Поскольку нашей целью является не вычисление инкремента для конкретного момента времени, а развитие общего метода решения дисперсионного уравнения для случая, когда ранее использовавшиеся методы неприменимы, мы приведем эти решения для нескольких значений плотности холодных частиц. Будем обозначать через ${{n}_{h}}$ и ${{n}_{c}}$ плотность горячих частиц, рассчитанную по экспериментально измеряемым потокам, и неизвестную плотность холодных частиц соответственно.

Пересчитав функцию распределения в переменных $({{v}_{\parallel }},{{v}_{ \bot }})$, можно рассчитать для случая ${{n}_{c}} = 0$ функции ${{f}_{{0\parallel }}}({{v}_{\parallel }})$, $g({{v}_{\parallel }})$ (19) и $Q({{v}_{{||}}},k)$ (28), входящие в дисперсионное уравнение. Функции ${{f}_{{0\parallel }}}({{v}_{\parallel }})$, $g({{v}_{\parallel }})$ и функция $Q({{v}_{{||}}},k)$ для $kc{\text{/}}{{\omega }_{c}} = 7.7$ и указанных выше момента времени и положения спутника MMS показаны на рис. 3.

Результаты вычислений дисперсионных характеристик и декрементов волн для случая ${{n}_{c}} = 0$ представлены на рис. 4 и 5. На рис. 4 показана зависимость действительной части частоты (левая панель) и декремента волн (правая панель) от волнового числа, которое всегда считается действительным. Мы видим, что одному и тому же значению волнового числа отвечают две различные частоты, что означает появление в этом случае новой волновой моды. Кроме того, на одной и той же непрерывной ветви колебаний, то есть в рамках одной моды, возникает немонотонная зависимость частоты от волнового числа, так что одна и та же частота может отвечать двум различным значениям волнового вектора. Поэтому, при ${{n}_{c}} = 0$ инкремент не может быть выражен как однозначная функция частоты даже в рамках одной волновой моды, хотя и частота и декремент являются при этом однозначными функциями волнового числа, что является признаком волновой моды. Из правой панели ясно видно, что высокочастотная волновая мода является сильно затухающей, также как и более низкочастотная ветвь в области достаточно больших значений волнового числа и заведомо в области, где частота волны понижается с ростом волнового числа. (Отметим, что частоты и декременты, отвечающие первой и второй моде (показанных на рисунке сплошной и штриховой линиями соответственно), и лежащие на непрерывных кривых, однозначно определяются соответствующими интервалами волновых чисел). Мы видим, что дисперсия и затухание волн в горячей плазме существенно отличаются от того, что имеет место в холодной плазме с малой примесью энергичных частиц (см. рис. 6). Это связано с тем, что в горячей плазме, в отличие от холодной, как действительная, так и мнимая часть частоты определяются трансцендентным дисперсионным уравнением вида (27). Появление новых ветвей колебаний в горячей плазме хорошо известно. В частности, в [25] указано, что в горячей плазме “появляется ряд новых ветвей колебаний магнитоактивной плазмы, слабо и сильно затухающих, длинно- и коротковолновых, причем все они отсутствуют в холодной плазме”.

На рис. 5 показана зависимость инкремента от частоты для указанного выше момента времени, где ясно видна неоднозначная зависимость $\gamma (\omega )$ в рамках одной волновой моды, связанная с немонотонной зависимостью $\omega (k)$. (Заметим, что необходимость монотонной зависимости частоты от волнового вектора в рамках волновой моды не является обязательной и не следует ни из каких общих соображений.) В предыдущие и последующие моменты времени (7 измерений, начиная с 38.505 с по 42.000 с) результаты несколько различаются количественно, но качественная картина сохраняет представленный на рисунках вид. Следовательно, небольшие изменения функции распределения, в том числе связанные с ошибками измерений, не приводят к качественным изменениям результатов.

Известно, что локальные особенности функции распределения могут влиять на устойчивость или неустойчивость плазмы. Неоднозначность же зависимости $\gamma (\omega )$ не связана с какими-либо особенностями функции распределения горячих частиц, что видно из приведенного рис. 3, а только с самим фактом отсутствия холодной компоненты плазмы. Это подтверждается тем, что при добавлении в распределение холодных частиц, немонотонная зависимость $\omega (k)$ исчезает при той же самой функции распределения горячих частиц. Если же холодные частицы отсутствуют, то как действительная, так и мнимая часть величины D из (26), задающей дисперсионное уравнение (27), определяются, в основном, входящим в D интегралом (43). В этом случае привычные представления о неустойчивом диапазоне и особенностях функции распределения, которые определяют этот диапазон, теряют свою силу, поскольку дисперсионное соотношение определяется не функцией распределения на действительной оси ${{v}_{{||}}}$, а аналитическим продолжением соответствующего интеграла в комплексную плоскость. В случае же достаточно большого ${{n}_{c}}$ и малых инкрементов неустойчивый диапазон в точности соответствует тем волновым векторам и соответствующим частотам, для которых резонансная скорость лежит в диапазоне скоростей, где величина Q положительна. Это связано с тем, что в этом случае мнимая часть величины D пропорциональна значению $Q({{v}_{R}})$.

Возникает естественный вопрос о влиянии на результат ошибок измерений. В связи с этим заметим, что входящие в дисперсионное уравнение функции ${{f}_{{0||}}}$ и g являются взвешенными интегралами по поперечной скорости ${{v}_{ \bot }}$ от измеряемой функции распределения. Поскольку из измерений мы знаем функцию распределения в переменных энергия (W)-питч-угол (α), то интегрирование в (19) по ${{v}_{ \bot }}$ фактически соответствует некоторому интегралу по кривой на плоскости $(W,\alpha )$ и, таким образом, уже включает некоторое усреднение, что видно из графиков функций ${{f}_{{0||}}}$ и g, приведенных на рис. 3. Следовательно, небольшие неточности измерений существенно сглаживаются и мало влияют на конечный результат. Это же подтверждают результаты расчётов в моменты времени, ближайшие к моменту времени выбранного события. Подчеркнем, что речь идет о моментах времени в интервале десятых долей секунды. Отметим также, что в расчетах мы не используем усреднение потоков во времени.

С ростом ${{n}_{c}}$ ситуация меняется. Вторая мода исчезает, и действительная часть частоты и волновое число становятся взаимно однозначными функциями, так что инкремент теперь может быть выражен как однозначная функция частоты. Для функции распределения в момент времени, для которого приведены вычисления инкремента, а именно, 19 июня 2017 г., 03:57:38.625, немонотонная зависимость $\omega (k)$ и связанная с ней неоднозначность зависимости $\gamma (\omega )$ исчезает при ${{n}_{c}}{\text{/}}{{n}_{h}} = 0.7$. Близкие значения относительной плотности холодной плазмы, при которых вторая мода исчезает, а функция $\omega (k)$ становится монотонно возрастающей функцией и, соответственно, $\gamma (\omega )$ – однозначной функцией, получаются и для других моментов времени. Результаты расчетов для ${{n}_{c}} = 3{{n}_{h}}$ и ${{n}_{c}} = 5{{n}_{h}}$ показаны на рис. 6, соответственно сплошной линией и штриховой линией. Мы видим, что в этом случае появляется неустойчивый диапазон частот, причем с ростом плотности холодной компоненты максимум инкремента смещается в сторону низких частот, а его абсолютная величина несколько понижается. Следует подчеркнуть, что в случае, когда плотность холодной компоненты плазмы на несколько порядков превышает плотность горячей компоненты, что имеет место во внутренней магнитосфере, максимум относительного инкремента (декремента) оказывается на несколько порядков меньше, чем для случая, когда холодная компонента плазмы отсутствует, или ее плотность соизмерима с плотностью горячей компоненты, что характерно для хвоста магнитосферы.

Отметим, что в случае, когда плотность холодной плазмы существенно превышает плотность горячей компоненты, а инкремент волны много меньше ее частоты, зависимость величины инкремента от плотности холодной плазмы ${{n}_{c}}$ хорошо известна. Для широкого класса функций распределения инкремент пропорционален экспоненте

(${{v}_{{{{T}_{\parallel }}}}}$ – продольная тепловая скорость), а предэкспоненциальный множитель пропорционален $1{\text{/}}n_{c}^{{3/2}}$ (см., например, формулу (3.2.6) из [7]). С ростом плотности холодной плазмы экспонента возрастает, в то время как предэкспоненциальный множитель убывает, так что окончательный результат зависит от конкуренции этих двух факторов. При ${v}_{R}^{2} \lesssim {v}_{{{{T}_{\parallel }}}}^{2}$, что характерно для горячей плазмы, основным фактором является уменьшение предэкспоненциального множителя, что ведет к уменьшению инкремента с ростом плотности холодной плазмы.

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

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

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

Использование развитого подхода проиллюстрировано в работе на двух примерах вычисления инкремента свистовых волн на основе спутниковых измерений внешнего магнитного поля и параметров плазмы, включая дифференциальные потоки электронов различных энергий и питч-углов. Первый пример относится к измерениям спутника Van Allen Probe-A в экваториальной области внутренней магнитосферы на $L = 3.21.$ В этом случае наш подход дает результаты, аналогичные тем, которые получаются при вычислении инкремента в предположении “холодной” дисперсии и малого инкремента, когда применима формула (45). Второй пример относится к вычислению инкремента на основе измерений одного из спутников MMS в хвосте магнитосферы, на расстоянии порядка $17{{R}_{e}}.$ Как показывает расчет, в этом случае результат существенно зависит от плотности “холодной” плазмы, которую не удается определить по имеющимся данным. Это было бы серьезным препятствием для достижения цели работы, если бы она состояла в вычислении инкремента в данном конкретном случае. Однако целью работы и ее основным результатом является развитие метода вычисления инкремента волны при произвольном соотношении между плотностями холодной и энергичной компонент плазмы, когда используемый обычно подход неприменим.

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

Работа выполнена при поддержке гранта РНФ 22-22-00135.