Wigner function method for describing electromagnetic field in plasma-like media with spatial dispersion and resonant dissipation

Е. D. Gospodchikov; Господчиков Е. Д.; A. A. Balakin; Балакин А. А.; A. G. Shalashov; Шалашов А. Г.

doi:10.31857/S0367292124080037

Wigner function method for describing electromagnetic field in plasma-like media with spatial dispersion and resonant dissipation

Authors: Gospodchikov Е.D.¹, Balakin A.A.¹, Shalashov A.G.¹
Affiliations:
1. Gaponov-Grekhov Institute of Applied Physics, Russian Academy of Sciences
Issue: Vol 50, No 8 (2024)
Pages: 866-886
Section: МАГНИТНЫЕ ЛОВУШКИ
URL: https://journal-vniispk.ru/0367-2921/article/view/280611
DOI: https://doi.org/10.31857/S0367292124080037
EDN: https://elibrary.ru/OBFRTH
ID: 280611

Cite item

Full Text

Abstract
Full Text
About the authors
References
Supplementary files
Statistics

Abstract

The paper gives a systematic presentation of the Wigner function method (Weyl formalism) for modeling the propagation and absorption of electromagnetic waves in anisotropic and gyrotropic dissipative media with spatial dispersion. A general kinetic equation for the Wigner function (tensor) is formulated, and its asymptotic expansion up to the second order for smoothly inhomogeneous and weakly dissipative media is constructed. As a result, a modification of the method of the kinetic equation for rays is proposed, based on the stochastic description of rays, making it possible to increase the accuracy of numerical modeling of wave problems with strong transverse inhomogeneity of the absorption coefficient without increasing the amount of calculations. The technique developed can be used to describe the propagation, absorption, and scattering of electron-cyclotron waves in high-temperature plasma of magnetic traps for controlled fusion in cases where standard modeling methods do not provide the necessary accuracy.

Keywords

asymptotic methods in electromagnetics, wave beams, wave propagation in inhomogeneous gyrotropic media, waves in high-temperature plasma, microwave heating of plasma

Full Text

1. Введение

Практическое моделирование распространения, поглощения и рассеяния волн электронно-циклотронного (ЭЦ) диапазона в турбулентной горячей магнитоактивной плазме, удерживаемой в больших магнитных ловушках УТС, традиционно сталкивается с проблемой огромного разброса пространственных масштабов – длина волны и характерный масштаб пространственной дисперсии волн оказываются на два–четыре порядка меньше геометрических размеров установок. Это обстоятельство делает полноволновое моделирование, опирающееся на прямое решение уравнений Максвелла, дополненных уравнениями для возбуждаемого в среде тока, крайне затруднительным.

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

$λ ≪ a ≪ L,$ (1)

где $λ$ $-$ длина волны, $a$ $-$ поперечный размер волнового пучка, $L$ $-$ масштаб неоднородности среды. Второе неравенство может нарушаться, например, в областях резонансного поглощения.

Не всегда решает указанные проблемы и следующий шаг в развитии асимптотических средств моделирования – beam-tracing [3]. В этом подходе используется представление об автомодельности гауссовского распределения поля в волновом пучке, справедливое для ряда упрощенных случаев. Для таких модельных задач фиксируется поперечная структура пучка, а управляющие параметры (ширина, кривизна волнового фронта и т. д.) находятся из решения эволюционных уравнений вдоль луча. При этом предполагается, что коэффициент поглощения волн неоднороден только вдоль несущего геометрооптического луча, что хорошо отвечает традиционной схеме электронного циклотронного нагрева термоядерной плазмы с поперечным экваториальным вводом мощности [4–6]. Такой подход хоть и позволяет достаточно правдоподобно восстановить распределение поля, страдает рядом врожденных недостатков, не позволяющих использовать его в тех случаях, когда коэффициент поглощения заметно неоднороден в поперечной плоскости пучка [7]. Так происходит, например, при вертикальном вводе излучения в тороидальные установки [8, 9]. Серьезную проблему для этого подхода создают эффекты отражения и формирования каустик на трассе распространения пучка и в области поглощения. Такие эффекты важны для нагрева плазмы в аксиально-симметричных магнитных ловушках [10], идейно родственных “методу магнитного берега” [11], и для случая так называемой “тангенциальной инжекции” для безындукционной генерации тока с высокой плотностью в тороидальных установках [12, 13]. При этом в рамках формализма beam-tracing могут возникать нефизические эффекты завышенного или отрицательного поглощения; с другой стороны, реальные эффекты совместного влияния дифракции и резонансной диссипации могут выпасть из результатов моделирования [7, 13, 14].

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

· численная схема решения эволюционного уравнения становится в задачах без естественной факторизации довольно громоздкой и требовательной к вычислительным ресурсам;

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

· по этой же причине затруднено прямое сравнение квазиоптического подхода с полноволновым моделированием, круг модельных задач, для которых это возможно, очень ограничен.

Развивается также подход, основанный на разложении Габора исходного скалярного поля на набор элементарных гауссовых пучков с последующим суммированием с учетом их взаимной фазы [21]. Этот метод позволяет описывать некоторые аберрационные эффекты, выходящие за пределы применимости стандартного beam-tracing [22].

Появление вычислительных мощностей, позволяющих проводить параллельные расчеты с высокой эффективностью, стимулировали развитие еще одного подхода к нахождению асимптотического решения уравнений Максвелла, основанного на решении “лучевого кинетического уравнения” для функции Вигнера [23], известного также как формализм Вейля [24]. Формально точно вычисленная функция Вигнера позволяет восстановить любую билинейную форму реального поля. Этот метод является адаптацией хорошо известного из квантовой механики и оптики вигнеровского подхода к построению квазиклассического приближения, позволяющего сводить уравнение Шредингера в некоторых приближениях к кинетическому уравнению Больцмана [25, 26]. Идейно этот подход близок и к многолучевой геометрической оптике, в которой волновой пучок моделируется целым набором независимых [27] или взаимодействующих [28] геометрооптических лучей. Такой подход не лишен недостатков и ограничений:

· по смыслу используемых разложений он должен классифицироваться как безаберрационный [29];

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

· в практических приложениях он требует очень большого числа лучей ( $10^{6}$ и выше) для обеспечения приемлемой точности [30].

· Тем не менее, метод обладает несомненными достоинствами:

· возможностью последовательного учета пространственной и временной дисперсии и диссипации;

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

Развитие кодов на основе этого подхода происходит в последнее время довольно стремительно [31–35].

Наиболее сложным эффектом для моделирования с использованием безаберрационных кодов является неоднородная на поперечном масштабе квазиоптического волнового пучка диссипация [7, 14] и пространственная дисперсия диссипации [36]. При этом значительная неоднородность коэффициента поглощения поперек волнового пучка характерна для целого ряда схем нагрева плазмы в тороидальных и аксиально-симметричных магнитных ловушках. Это связано как с резонансным характером поглощения СВЧ-волн, так и с релятивистскими эффектами, приводящими к негладким профилям поглощения (при достаточно слабом магнитном поле поглощение отсутствует полностью). С математической точки зрения это выражается в том, что компоненты тензора диэлектрической проницаемости представляются неаналитическими функциями ([37], $§$ 16 гл. 4).

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

Структура работы следующая. В разд. 2 вводится формализм, используемый далее $-$ представление Вейля для псевдодифференциальных линейных операторов и функция Вигнера для скалярного поля. В разд. 3 вводится тензорная функция Вигнера для монохроматического электромагнитного поля и формулируется кинетическое уравнение для нее. В разд. 4 и 5 вводится скалярная функция Вигнера для электромагнитного поля и обсуждается ее связь с квазиоптическим приближением. В разд. 6 рассматривается первый асимптотический порядок в плавнонеоднородной и слабодиссипативной среде, приводящий к кинетическому уравнения Больцмана для скалярной функции Вигнера (кинетический лучевой подход). В разд. 7 рассматривается второй порядок $-$ показывается, что эффекты неоднородного поглощения приводят к переходу от уравнения Больцмана к уравнению Фоккера $-$ Планка, которое может быть решено добавлением стохастической силы в уравнения Гамильтона для лучей (диффузионное приближение). В разд. 8 обсуждаются начальные условия для кинетического уравнения. В разд. 9 выводятся нетривиальные для сред с пространственной дисперсией выражения для плотностей потока энергии и поглощаемой мощности через введенные функции Вигнера. В разд. 10 приводится пример использования развитой техники для решения задачи о распространении волнового пучка через область замагниченной плазмы с неоднородным циклотронным поглощением в плазме токамака. Пример носит демонстрационный характер и не имеет самостоятельной ценности, однако он показывает разницу между первым и вторым порядком асимптотической теории, проявляющейся из-за резонансной диссипации. В Заключении подводятся итоги работы.

2. Функция вигнера для скалярного поля

Начнем с введения формализма Вейля для псевдодифференциальных операторов. Линейный оператор, действующий на скалярное комплексное поле $ψ (r')$ , можно представить через интегральное ядро, возможно содержащее обобщенную дельта-функцию Дирака и ее производные, как

$\hat{A} ψ = \int A (r, r') ψ (r') d r' .$ (2)

Здесь $d r' = d x' d y' d z'$ обозначает элемент объема (а не приращение вектора); в таком же смысле мы будем понимать дифференциалы от векторных величин во всех остальных интегральных соотношениях, используемых в статье. Эрмитово-сопряженный оператор ${\hat{A}}^{+}$ при стандартном определении скалярного произведения отвечает ядру $A^{*} (r', r)$ , где звездочкой обозначена операция комплексного сопряжения.

Будем называть вейлевским представлением оператора следующую функцию [26, 38]:

$A (R, K) = \int A (R + ξ / 2, R - ξ / 2) exp (- i K ξ) d ξ .$ (3)

Здесь и далее каллиграфическое написание относится к ядрам операторов в форме (2), а аналогичное курсивное – к их вейлевскому представлению (3). Обратное выражение для ядра линейного оператора через его вейлевское представление есть

$A (r, r') = \int A (\frac{r + r'}{2}, K) exp (i K (r - r')) \frac{d K}{{(2 π)}^{3}} .$ (4)

Таким образом возникает взаимно-однозначное соответствие между операторами (2) и их вейлевскими представлениями (3).

Вычисляя комплексное сопряжение от формулы (3), получаем

$A^{*} (R, K) = \int {[A (R + ξ / 2, R - ξ / 2) exp (- i K ξ)]}^{*} d ξ =$

$= \int A^{*} (R - ξ / 2, R + ξ / 2) exp (- i K ξ) d ξ =$

$= \int A^{+} (R + ξ / 2, R - ξ / 2) exp (- i K ξ) d ξ .$

То есть комплексное сопряжение вейлевского представления линейного оператора порождает эрмитово-сопряженный оператор ${\hat{A}}^{+}$ . Из этого следует, что эрмитову оператору отвечает действительное вейлевское представление, а антиэрмитовому оператору – чисто мнимое представление.

Рассмотрим, например, оператор умножения на функцию координат

$\hat{Q} [ψ] \equiv q (r) \cdot ψ = \int ψ (r^{'}) q (r) δ (r^{'} - r) d r^{'} .$

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

$Q = \int q (R + ξ / 2) δ (ξ) e^{- i K ξ} d ξ = q (R) .$

Другой пример. Операторы дифференцирования

${\hat{D}}_{_{0}} [ψ] \equiv \frac{\partial ψ}{\partial r}, {\hat{D}}_{_{1}} [ψ] = \frac{\partial}{\partial r} r ψ, {\hat{D}}_{_{2}} [ψ] = r \frac{\partial}{\partial r} ψ$

имеют вейлевские представления

$D_{0} = \int e^{- i K ξ} \nabla δ (ξ) d ξ = i K,$

$D_{1,2} = \int (R \mp ξ / 2) e^{- i K ξ} \nabla δ (ξ) d ξ = i K R \mp \frac{1}{2} .$

При этом действительное представление $K R$ будет отвечать эрмитовому оператору $({\hat{D}}_{1} + {\hat{D}}_{2}) / 2 i$ .

Вейлевское представление $(A B)$ суперпозиции $\hat{A} \hat{B}$ линейных операторов $\hat{A}$ и $\hat{B}$ имеет вид

$(A B) = \int A (R + \frac{r_{1}}{2}, K + k_{1}) B (R + \frac{r_{2}}{2}, K + k_{2}) \times$

$\times exp (i k_{2} r_{1} - i k_{1} r_{2}) \frac{d r_{1} d k_{1} d r_{2} d k_{2}}{{(2 π)}^{6}} .$

Разложим подынтегральные функции в ряды по $r_{1}$ , $k_{1}$ , $r_{2}$ , $k_{2}$ и выполним интегрирование по координатам $r_{1}$ и $r_{2}$ . В результате получим

$(A B) = \sum_{m, n, k, l = 0}^{\infty} \frac{1}{m! n! k! l!} \int {(\frac{1}{2 i} \frac{\partial_{A}}{\partial R} \frac{\partial}{\partial k_{2}})}^{m} {(k_{1} \frac{\partial_{A}}{\partial K})}^{n}$

$\times {(\frac{i}{2} \frac{\partial_{B}}{\partial R} \frac{\partial}{\partial k_{1}})}^{k} {(k_{2} \frac{\partial_{B}}{\partial K})}^{l} A B δ (k_{1}) δ (k_{2}) d k_{1} d k_{2} .$

Нижний индекс в частных производных определяет “направление” дифференцирования: дифференциальный оператор, записанный через $\partial_{A}$ , применяется только к функции $A (R, K)$ , а дифференциальный оператор $\partial_{B}$ применяется только к $B (R, K)$ . Слагаемые в этой сумме отличны от ноля только при $m = l$ и $n = k$ . После этого мы получаем ответ в виде ряда, который может быть формально записан через экспоненту от оператора:

$(A B) = \sum_{m, n = 0}^{\infty} \frac{1}{m! n!} {(\frac{i}{2} \frac{\partial_{A}}{\partial R} \frac{\partial_{B}}{\partial K})}^{m} {(\frac{1}{2 i} \frac{\partial_{B}}{\partial R} \frac{\partial_{A}}{\partial K})}^{n} A B =$

$= exp (\frac{i}{2} \frac{\partial_{A}}{\partial R} \frac{\partial_{B}}{\partial K} - \frac{i}{2} \frac{\partial_{B}}{\partial R} \frac{\partial_{A}}{\partial K}) A (R, K) B (R, K) .$ (5)

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

Для произвольного комплексного поля $ψ (r)$ можно определить функцию Вигнера [25]

$W (R, K) = \frac{1}{{(2 π)}^{3}} \int ψ (R + \frac{ξ}{2}) ψ^{*} (R - \frac{ξ}{2}) e^{- i K ξ} d ξ .$ (6)

Функция Вигнера принимает действительные значения. Кроме того, она имеет следующие свойства:

${|ψ (r)|}^{2} = \int W (r, K) d K,$

${|ψ_{k} (k)|}^{2} = \frac{1}{{(2 π)}^{3}} \int W (R, k) d R,$

$ψ (r) ψ^{*} (r') = \int W (\frac{r + r'}{2}, K) e^{i K (r - r')} d K,$ (7)

$⟨ψ A ψ⟩ \equiv \iint ψ (r) ψ^{*} (r') A (r, r') d r d r' =$

$= \iint W (R, K) A (R, K) d R d K .$

Здесь $ψ_{k}$ есть Фурье-образ $ψ$ . Эти свойства позволяют в квантовой механике использовать функцию Вигнера как аналог функции распределения в фазовом пространстве [25, 26]. Заметим, что таким образом определенная функция хотя и является действительной, не обязана быть знакоопределенной, что может нарушать некоторые интуитивные представления, связанные с обычной функцией распределения.

Рассмотрим представление решения произвольного дифференциального или интегро-дифференциального линейного уравнения через функцию Вигнера. Пусть имеется уравнение

$\hat{D} [ψ] = 0,$ (8)

где $\hat{D}$ $-$ линейный оператор с интегральным ядром $D (r, r')$ . Это уравнение равносильно операторному уравнению

$\hat{D} \hat{W} = 0, W (r, r') = \frac{1}{{(2 π)}^{3}} ψ (r) ψ^{*} (r'),$ (9)

где $\hat{W}$ $-$ линейный оператор с интегральным ядром $W (r, r')$ . Необходимость условия (9) для выполнения (8) следует из очевидного соотношения

$ψ^{*} \hat{D} [ψ] = 0 \Rightarrow \int D (r, r^{'}) ψ (r') ψ^{*} (r^{″}) d r^{'} d r^{″} = 0.$

Его достаточность следует из соотношения

$\hat{W} [ψ] = μ ψ, μ = \frac{1}{{(2 π)}^{3}} \int {|ψ|}^{2} d r^{'} > 0$

и линейности оператора $\hat{D}$ . Таким образом, уравнения (8) и (9) полностью эквивалентны.

Вейлевское представление $W (R, K)$ оператора $\hat{W}$ в точности совпадает с функцией Вигнера (6) для поля $ψ (r)$ . Воспользовавшись формулой (5) для суперпозиции линейных операторов, уравнение (9) можно переписать в вейлевском представлении как

$exp (\frac{i}{2} \frac{\partial_{D}}{\partial R} \frac{\partial_{W}}{\partial K} - \frac{i}{2} \frac{\partial_{W}}{\partial R} \frac{\partial_{D}}{\partial K}) D (R, K) W (R, K) = 0.$ (10)

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

3. Тензор вигнера для монохроматического электромагнитного поля

Поскольку функция Вигнера квадратична по полю, то комплексное векторное поле порождает тензорную функцию Вигнера. В частности, для распределения комплексной амплитуды монохроматического электрического поля $E$ можно ввести тензорную функцию Вигнера (тензор Вигнера)

$\begin{matrix} W_{α β} (R, K) = \frac{1}{{(2 π)}^{3}} \int E_{α} (R + \frac{1}{2} ξ) \cdot \\ \cdot E_{β}^{*} (R - \frac{1}{2} ξ) e^{- i K ξ} d ξ . \end{matrix}$ (11)

Греческими буквами $α$ , $β$ , ... мы будем обозначать индексы, нумерующие декартовы компоненты электрического поля. Эти индексы во всех последующих выкладках пробегают значения x, y и z.

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

${|E (r)|}^{2} = \sum_{α} \int W_{α α} (r, K) d K,$

$E_{α} (r) E_{β}^{*} (r^{'}) = \int W_{α β} (\frac{r + r^{'}}{2}, K) e^{i K (r - r^{'})} d K .$ (12)

Рассмотрим систему волновых уравнений на комплексные амплитуды электрического поля (три декартовы компоненты $E_{α}$ ), следующую из уравнений Максвелла для монохроматического поля с частотой $ω$ . В анизотропной линейной среде с тензором диэлектрической проницаемости $ε_{α β} (ω)$ эта система определяется как

$\sum_{β} {\hat{D}}_{α β} [E_{β}] = 0,^{} {\hat{D}}_{α β} = \frac{\partial^{2}}{\partial x_{α} \partial x_{β}} - δ_{α β} Δ - k_{0}^{2} {\hat{ε}}_{α β} .$

(например, [37], $§$ 3 гл. 2). Здесь D $-$ оператор Лапласа, $k_{0} = ω / c$ $-$ размерная постоянная, имеющая смысл волнового числа в вакууме для монохроматического волнового поля. В соответствии с результатами предыдущего раздела, эта система уравнений эквивалента матричной системе уравнений

$\sum_{β} {\hat{D}}_{α β} {\hat{W}}_{β γ} = 0,$

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

$\begin{matrix} exp (\frac{i}{2} \frac{\partial_{D}}{\partial R} \frac{\partial_{W}}{\partial K} - \frac{i}{2} \frac{\partial_{W}}{\partial R} \frac{\partial_{D}}{\partial K}) \\ \sum_{β} D_{α β} (R, K) W_{β γ} (R, K) = 0, \end{matrix}$ (13)

где

$D_{α β} (R, K) = - K_{α} K_{β} + δ_{α β} K^{2} - k_{0}^{2} ε_{α β} (R, K) .$ (14)

Получившееся матричное кинетическое уравнение, аналогичное (10), позволяет строить решения для тензорной функции Вигнера с любой необходимой точностью. Заметим, что в средах с пространственной дисперсией тензор диэлектрической проницаемости должен задаваться нелокальным интегральным оператором вида (2) с ядром $ε_{α β} (r, r^{'})$ ; в средах без пространственной дисперсии действие тензора диэлектрической проницаемости сводится к умножению на числовую матрицу $ε_{α β}$ .

4. Нулевой порядок. Скалярная функция Вигнера для электромагнитного поля

В плавнонеоднородной среде, в которой масштабы неоднородности среды и характерные размеры распределения интенсивности комплексной амплитуды поля много больше длины волны, $λ / L ≪ 1$ , возможно построение упрощенной асимптотической теории на базе кинетического уравнения (13). При этом следует учесть, что каждая операция $\partial^{2} / \partial R \partial K$ в операторной экспоненте приводит к повышению порядка малости по параметру геометрической оптики $λ / L$ .

В нулевом порядке теории возмущений операторная экспонента в (13) заменяется тождественным оператором, что дает

$\sum_{β} D_{α β} (R, K) W_{β γ} (R, K) = 0.$ (15)

Это матричное тождество, которое должно выполняться покомпонентно для всех $α$ и $γ$ . Для решения этого уравнения нам потребуется использовать дополнительные физические соображения, ограничивающие вид волнового оператора. Во-первых, заметим, что для существования нетривиального решения для тензора Вигнера необходимо, чтобы одно из собственных значений матрицы волнового оператора $D_{α β}$ принимало нулевое значение при определенных R и K. Обозначим это собственное значение H(R,K). Во-вторых, будем считать, что в рассматриваемой задаче отсутствует поляризационное вырождение, то есть для всех R и K выделенному собственному значению Н отвечает ровно один собственный вектор $e (R, K)$ с компонентами $e_{α}$

$\sum_{β} D_{α β} e_{β} = H e_{α} .$ (16)

Физически это означает, что мы не рассматриваем изотропные среды и исключаем специальный случай линейного взаимодействия нормальных волн, которое реализуется в окрестности точки поляризационного вырождения уравнений Максвелла в анизотропных и гиротропных средах с глобально снятым поляризационным вырождением [39]. Вектор $e (R, K)$ можно трактовать как вектор поляризации электрического поля для нормальной моды среды, а выполняющееся на нетривиальных решениях условие $H (R, K) = 0$ $-$ как соответствующее геометрооптическое дисперсионное соотношение, с учетом сделанных замечаний полностью эквивалентное условию разрешимости волнового уравнения $det (D_{α β}) = 0$ . Из (16) следует очевидное равенство

$H (R, K) = \sum_{α, β} e_{α}^{*} e_{β} D_{α β},$ (17)

справедливое, если вектор поляризации нормирован на единицу, $e e^{*} = 1$ .

С формально математической точки зрения приведенные выше условия означают, что размерность ядра волнового оператора $D_{α β}$ не может быть больше единицы и в точности равна ей в решениях (15). Поэтому образ оператора $W_{β γ}$ , отвечающего решению (15), должен состоять из одного (с точностью до растяжений) элемента $e (R, K)$ . Другими словами, $W_{β γ}$ должен проектировать любой вектор на вектор поляризации $e$ :

$\sum_{γ} W_{β γ} x_{γ} = μ e_{β} (\forall x_{γ}) .$

Эрмитов оператор, удовлетворяющий этому свойству, имеет общий вид

$W_{β γ} (R, K) = e_{β} (R, K) e_{γ}^{*} (R, K) W (R, K) .$ (18)

Чтобы показать это достаточно перейти в базис, первым из векторов которого является вектор поляризации $e$ . В этом базисе матрица $W_{β γ}$ имеет единственный ненулевой элемент $W_{11} = W$ .

Заметим, что тензорное представление (18) описывает электрическое поле вида

$E (r) = \hat{e} (r) u (r),$ (19)

где $\hat{e}$ есть оператор поляризации, однозначно определяемый по своему вейлевскому представлению $e (R, K)$ , а $u (r)$ есть скалярная комплексная амплитуда волнового поля. Функция $W (R, K)$ в точности соответствует функции Вигнера, отвечающей распределению скалярной комплексной амплитуды $u (r)$ в соответствии с (6), поэтому будем называть эту величину скалярной функцией Вигнера. Эта функция принимает только действительные значения.

Напомним, что формула (18) определяет решение уравнения (15) только если $e$ отвечает ядру волнового оператора, поэтому

$W (R, K) = \tilde{W} (R, K) δ (H (R, K)),$ (20)

где $\tilde{W} (R, K)$ $-$ произвольная действительная функция, $δ (\dots)$ $-$ дельта-функция Дирака, отвечающаа за выполнение дисперсионного соотношения. Таким образом, мы построили точное решение кинетического уравнения для тензора Вигнера в нулевом порядке теории возмущений.

5. Связь с квазиоптическим приближением

С методической точки зрения полезно представить матричные уравнения (13) как одно уравнение на скалярную функцию Вигнера, заданную в виде в виде (18)

$exp (\frac{i}{2} \frac{\partial_{H}}{\partial R} \frac{\partial_{W}}{\partial K} - \frac{i}{2} \frac{\partial_{W}}{\partial R} \frac{\partial_{H}}{\partial K}) H (R, K) W (R, K) = 0.$ (21)

Уравнениям нулевого порядка (15) соответствует

$H (R, K) W (R, K) = 0.$

Эквивалентное (21) уравнение на электрическое поле, действующее в исходном (не вигнеровском) пространстве, есть

$\hat{H} [u] = 0,$ (22)

где $\hat{H}$ $-$ оператор, заданный вейлевским представлением (17), которое в свою очередь определяется как собственное число вейлевского представления исходного волнового оператора $D_{α β}$ . Оператор $\hat{H}$ будем называть оператором Гамильтона уравнения на скалярную амплитуду. Это определение несколько условно, поскольку наш оператор не обязательно эрмитов. Как обычно, эрмитова часть оператора $\hat{H}$ отвечает за “обратимый” диэлектрический отклик и, как мы увидим ниже, порождает гамильнонову динамику в консервативной среде, а антиэрмитова часть отвечает за необратимые процессы $-$ диссипацию электромагнитного поля или его усиление в неравновесной среде.

Представление поля в плавнонеоднородной среде в виде невзаимодействующих мод вида (19) с фиксированными поляризациями и неизвестными скалярными амплитудами, определяемые уравнением (22), соответствует идеям аберрационного квазиоптического приближения, развитого в работах [15–19, 40]. Обсуждая связь нашего подхода с последовательной квазиоптической теорией, необходимо отметить два обстоятельства.

Во-первых, уравнение (22), которое в цитируемых выше работах рассматривается как исходное, вводится хотя и естественно, но в значительной степени феноменологически. Строгий вывод этого уравнения в рамках квазиоптики возможен лишь в простейших частных случаях. В неоднородных средах с пространственной дисперсией и диссипацией остается определенный произвол в определении как волнового гамильтониана $\hat{H}$ , так и скалярной амплитуды $u$ (ее связи с вектором электрического поля) даже при условии, что все физические критерии применимости квазиоптического приближения выполнены. Как было показано выше, метод функций Вигнера открывает возможность формального определения этих величин. Подчеркнем, что связь с квазиоптикой следует прежде всего из соотношения (16) и следующего из него “точного” дисперсионного соотношения $H (R, K) = 0$ . С другой стороны, вектора поляризации нормальных волн в стандартном геометрооптическом приближении определяются как собственные векторы матрицы $D_{α β}$ , отвечающие нулевому собственному значению. Отсюда можно сделать нетривиальное и отсутствующее в цитируемых работах по квазиоптике утверждение о том, что вейлевское представление оператора $\hat{e}$ , фигурирующего в (19), на самом деле является вектором поляризации одной из нормальных волн в геометрооптическом приближении. В анизотропной или гиротропной среде без пространственной дисперсии каждому направлению $K$ отвечает ровно два решения дисперсионного уравнения $-$ обыкновенная и необыкновенная электромагнитная волна. В среде с пространственной дисперсией общее число нормальных мод определяется количеством веток дисперсионного уравнения, при этом, если исключить линейное взаимодействие, каждому отдельному корню дисперсионного уравнения $H = 0$ отвечает отдельный вектор поляризации (отдельная мода).

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

$\hat{H} [\bar{u} e^{\int i k_{τ} d τ}]= 0 \Rightarrow \frac{\partial \bar{u}}{\partial τ} \approx i {\hat{H}}_{QO}[\bar{u}] .$ (23)

При этом возникают естественные ограничения на максимальную кривизну опорного луча, которые могут ограничить применение метода, например, при расчете распространения квазистатических сильнозамедленных мод в теплой плазме. Кроме того, в процессе укорочения волнового гамильтониана возникают дополнительные формальные сложности с определением антиэрмитовой части ${\hat{H}}_{QO}$ , ответственной за резонансную диссипацию, которые приходится преодолевать феноменологическими приемами [17–19]. Как будет показано ниже, асимптотическая теория, основанная на кинетическом уравнении для функции Вигнера, позволят обойти эти трудности, при этом опираясь только на лучевую оптику и уравнение переноса интенсивности вдоль луча. Однако за это приходится “расплачиваться” тем, что этот метод дает информацию только о квадратичных по электрическому полю усредненных величинах.

6. Первый порядок. Кинетический лучевой подход

Мы пришли к квазиоптическому уравнению для скалярной амплитуды, опираясь на нулевое приближение в строгом кинетическом уравнении для тензора Вигнера. Следующее приближение мы будем строить так, чтобы сохранить квазиоптическое представление поля, которое уже показало свою применимость для широкого круга задач. С физической точки зрения это означает, что мы по прежнему пренебрегаем линейным взаимодействием мод, а также рассеянием во встречную волну для заданной моды $e$ . C формально-математической точки зрения это означает, что мы “навязываем” диагональную форму для тензора Вигнера (18). То есть мы будем учитывать только поправки к диагональным членам, возникающие при разложении матричного кинетического уравнения (13) по малому параметру $λ / L$ , при этом потребуется последовательно избавляться не только от старших порядков, но и от всех членов, приводящих к возникновению внедиагональных компонент $W_{i j}^{0}$ и связи разных диагональных элементов $W_{i i}$ .

Для этого будем сразу искать решение уравнений (13) в виде (18). Получим скалярное кинетическое уравнение (21), которое будем рассматривать как исходное уравнение для построения асимптотической теории для скалярной функции Вигнера $W$ .

Разложим операторную экспоненту в (21) в ряд Тейлора до первого порядка по $λ / L$ . Кроме того, дополнительно мы будем считать среду слабодиссипативной, т.е. полагать малой мнимую часть представления Вейля для оператора Гамильтона (17)

$|H^{″}| \underset{˜}{<} (λ / L) |H^{'}| .$ (24)

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

$H^{'} W = 0,$ (25)

$\frac{\partial H^{'}}{\partial K} \frac{\partial W}{\partial R} - \frac{\partial H}{\partial R} \frac{\partial W}{\partial K} = 2 H^{″} W .$ (26)

Первое уравнения соответствует приближению нулевого порядка, в котором отброшены диссипативные члены. Его решение аналогично (20):

$W = \tilde{W} (R, K) δ (H^{'} (R, K)),$ (27)

где $\tilde{W} (R, K)$ $-$ произвольная действительная функция, которая должна уточняться из второго уравнения, а выполняющееся на решениях (25) условие

$H^{'} (R, K) = 0$ (28)

воспроизводит геометрооптическое дисперсионное соотношение в слабодиссипативной среде.

Уравнение Больцмана (26) решается методом характеристик. Рассмотрим функции R(t) и K() $R (τ)$ и $K (τ)$ , удовлетворяющие уравнениям Гамильтона

$\frac{d R}{d τ} = \frac{\partial H^{'}}{\partial K}, \frac{d K}{d τ} = - \frac{\partial H^{'}}{\partial R} .$ (29)

Волновой вектор K здесь играет роль канонического импульса, сопряженного координате R. Для величины

$W (τ) = \int W (R, K) δ (R - R (τ)) δ (K - K (τ)) d R d K$

получаем уравнение

$\frac{d W}{d τ} = 2 H^{″} W .$ (30)

Его решение вдоль характеристики имеет вид

$W (τ) = W_{0} exp (\int_{0}^{τ} 2 H^{″} (τ^{'}) d τ^{'}) .$ (31)

$W_{0}$ определяется начальными условиями при $τ = 0$ , которые задаются на некоторой поверхности Г₀в шестимерном фазовом пространстве. Области с $H^{″} < 0$ отвечают за диссипацию поля, области с $H^{″} > 0$ отвечают за усиление в неравновесных средах.

Чтобы восстановить значение $W (R, K) = W (τ^{*})$ в некоторой заданной точке фазового пространства, достаточно посчитать “лучевое расстояние” $τ^{*} (R, K)$ от этой точки до поверхности Г₀. Для этого нужно просчитать решения уравнений (29) с начальными условиями $(R, K)$ в обратном направлении вплоть до момента, когда характеристика достигнет Г₀, т.е. надо найти решение $R (τ), K (τ)$ такое, что

$R (0) = R, K (0) = K, (R (- τ^{*}), K (- τ^{*})) \in$ Г₀.

Заметим, что в практических расчетах обычно требуется вычислить сразу все поле $W (R, K)$ или его интегральные характеристики, поэтому обратная трассировка лучей не эффективна. Вместо этого можно выпустить набор лучей с поверхности начальных условий Г₀и некоторым образом фиксировать прохождение лучей с учетом набираемого веса $W (τ)$ через разбитое на ячейки фазовое пространство.

Уравнения характеристик сохраняют гамильтониан $H'$ . Поэтому, уравнение (30) и его решение применимы как к $W$ , так и к введенной в (27) функции $\tilde{W}$ (при этом меняется смысл и размерность $W_{0}$ ). Первый случай реализуется, если для начальных условий выполнено дисперсионное соотношение (28), при этом $W (τ)$ определяет решение обоих уравнений (25), (26). Однако с технической точки зрения нам будет удобнее оставить $δ (H')$ в явном виде, поэтому решение системы уравнений (25), (26) представим в виде

$W (R, K) = W (τ^{*} (R, K)) δ (H' (R, K)),$ (32)

где $W (τ)$ определяется (31). Способ задания $W_{0}$ мы обсудим в разд. 8.

Покажем, что полученное решение не зависит от нормировки гамильтониана. Рассмотрим новую функцию Гамильтона $\tilde{H} = C H$ , полученную умножением исходной на произвольную положительную функцию $C (R, K)$ . Очевидно, такое умножение не влияет на выполнение условия (28), а в уравнениях на характеристики (29) может быть скомпенсировано переопределением эволюционной координаты $d \tilde{τ} = C d τ$ . Поэтому ${\tilde{H}}^{″} d \tilde{τ} = H^{″} d τ$ , то есть аргумент экспоненты в (31) является инвариантом нашего преобразования, а ослабление или усиление поля вдоль луча не зависит от нормировки $H$ .

Уравнения (29) в совокупности с решением (32) принято называть кинетическим лучевым подходом. Отметим отличия описываемого подхода от метода уравнения эйконала в геометрической оптике [41]. Во-первых, уравнение эйконала действует в трехмерном пространстве r, а уравнения на функцию Вигнера в шестимерном пространстве $(R, K)$ . В отличие от трехмерного пространства, в шестимерном пространстве для стационарного гамильтониана фазовые траектории не пересекаются. Поэтому в кинетическом подходе нет проблемы “схлопывания силовой трубки” и связанной с ней неоднозначности решения. В рамках метода уравнения эйконала эта проблема также может решаться переходом к другой трехмерной проекции шестимерного фазового пространства [42, 43]. Во-вторых, не возникает необходимости интегрирования фазы вдоль луча. В-третьих, уравнение переноса сводится к тривиальному, но в определенном смысле строгому виду (30) в отличие от приближения геометрической оптики, в котором уравнение переноса излучения вдоль луча либо феноменологическое, либо должно обосновываться достаточно сложными внешними способами [44].

7. Второй порядок. Диффузионное приближение

Учет следующего порядка по малому параметру $λ / L$ в уравнении (21) приводит к поправке диффузионного типа. Вместо уравнения (26) мы получаем

$\frac{\partial H^{'}}{\partial P} \frac{\partial W}{\partial Q} = 2 H^{″} W + \sum_{m, n} F_{m n} \frac{\partial^{2} W}{\partial Q_{m} \partial Q_{n}},$ (33)

$F_{m n} = - \frac{1}{4} \frac{\partial^{2} H^{″}}{\partial P_{m} \partial P_{n}} .$

Здесь для сокращения записи введены 6-компонентные векторы

$Q = (R_{x}, R_{y}, R_{z}, K_{x}, K_{y}, K_{z}),$

$P = (K_{x}, K_{y}, K_{z}, - R_{x}, - R_{y}, - R_{z}) .$

Левая часть уравнения (33) совпадает с (26). Член со вторыми производными функции $W$ в правой части описывает диффузию в импульсном пространстве, если поглощение неоднородно, и в координатном пространстве, если присутствует пространственная дисперсия поглощения.

Для физической корректности уравнения (33) в диссипативных средах наряду с знакоопределенностью мнимой части гамильтониана (H˝ < 0) необходимо также потребовать положительной определенности квадратичной формы $F_{m n}$ . Следует также отметить, что малость поглощения в условии (24) не всегда гарантирует малость пространственных производных $H''$ в $F_{m n}$ , которые при резонансном характере поглощения могут быть весьма значительными [45–47]. Кроме того, на длинных трассах распространения эффект диффузии функции Вигнера будет “накапливаться”, что может приводить к качественно новым эффектам. Особенно ярко это проявляется в случаях распространения вдоль резонансных поверхностей [7, 13, 14].

Формально появление диффузионных членов лишает кинетический лучевой подход главного численного преимущества $-$ возможности независимого расчета характеристик и эволюции функции Вигнера вдоль них. В этом случае численное решение уравнения (33) становится едва ли не сложнее, чем прямое полноволновое моделирование. Однако если рассматривать уравнение (33) как аналог кинетического уравнения на некоторую функцию распределения, то члены, отвечающие за учет неоднородности поглощения, могут интерпретироваться как переход от уравнения Больцмана к уравнению Фоккера $-$ Планка. Поэтому напрашивающимся эффективным численным методом решения обозначенной выше проблемы является переход к стохастическому описанию динамики отдельных лучей, приводящей к диффузии после усреднения по статистическому ансамблю [48, 49].

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

$\frac{d Q}{d τ} = \frac{\partial H^{'}}{\partial P} + G, ⟨ G_{m} ⟩ = 0, ⟨ G_{m} G_{n} ⟩ = F_{m n},$

где $G$ $-$ стохастическая сила в 6-ти мерном пространстве, обладающая нулевым средним и тензором дисперсии $F_{m n}$ . Эффективный способ генерации псевдослучайной последовательности, обладающей указанными свойствами, изложен в Приложении 1. Выделяя координаты и импульсы по отдельности, получаем лучевые уравнения в привычной форме Гамильтона

$\frac{d R}{d τ} = \frac{\partial H^{'}}{\partial K} + G_{R}, \frac{d K}{d τ} = - \frac{\partial H^{'}}{\partial R} + G_{K},$ (34)

но с добавочным ланжевеновскими источниками G_R и G_K, представляющими, соответственно, первые три и последние три компоненты случайного вектора $G$ . Если интенсивность $W$ вдоль луча меняется по закону (31), то в приближении Фоккера $-$ Планка функция распределения для такого семейства лучей

$W (R, K) = ⟨W_{0} exp (\int_{0}^{τ^{*} (R, K)} 2 H^{″} (τ^{'}) d τ^{'})⟩$

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

Здесь необходимо вспомнить, что кроме диффузионного кинетического уравнения, нам необходимо обеспечить выполнение уравнения (25) и следующего из него дисперсионного соотношения (28). Для этого стохастическая сила G должна обеспечивать сохранение гамильтониана, что возможно только если диффузия происходит в пространстве меньшей размерности. Возможность понизить размерность пространства, в котором происходит диффузия, связана с тем, что функция Вигнера на “соседних” лучах может меняться заметно вследствие накопления эффектов диффузии на трассе, в то время как неоднородность функции Вигнера вдоль луча определяется только локальной диссипацией. Поэтому вдоль луча функция Вигнера меняется плавно. Наиболее простой способ учесть подобные соображения $-$ это полностью запретить стохастическую диффузию вдоль луча. Этого можно добиться модификацией случайной силы в уравнениях (34):

${\tilde{G}}_{R} = G_{R} - (G_{R}, l) l,$

${\tilde{G}}_{K} = G_{K} - ((G_{K}, l) + {|\frac{\partial H^{'}}{\partial K}|}^{- 1} ({\tilde{G}}_{R}, \frac{\partial H^{'}}{\partial R})) l,$

где $l = (\partial H^{'} / \partial K) / |\partial H^{'} / \partial K|$ $-$ единичный вектор вдоль траектории. При таком выборе силы стохастически возмущаются только величины, поперечные к лучу, изменение продольного импульса при этом подчиняется закону сохранения гамильтониана.

Формально такой переход может быть строго (но очень громоздко) осуществлен на основе теоремы о понижении порядка автономной гамильтоновой системы ([50], $§$ 141). Отметим также, что предложенное модельное описание “продольной” части диффузии в импульсном пространстве приводит к отказу от учета возможного отражения луча от области неоднородности (в частности, от резонанса среды), то есть перескоку с одной ветки решения уравнения $H^{'} = 0$ на другую, отвечающую волне той же поляризации, но распространяющейся в обратном направлении [46, 47].

8. Начальные условия для кинетического уравнения

Для того, чтобы воспользоваться полученными решениями в практических расчетах, необходимо задать значение скалярной функции Вигнера при $τ = 0$ для всего пятимерного набора характеристик Г₀ (поверхность в шестимерном пространстве). Поскольку лучевые уравнения (29), (30) и (34) не содержат переменную $τ$ явным образом, мы имеем право выбирать поверхность Г₀как нам удобно.

Пусть распределение комплексной амплитуды поля задано на плоскости $z = 0$ , причем волна распространяется из полупространства $z < 0$ в полупространство $z > 0$ . Будем отсчитывать момент $τ = 0$ для каждой характеристики от момента, когда она пересекает плоскость $z = 0$ . Координаты и волновыве вектора (импульсы) на этой плоскости задаются двумерными (поперечными) векторами $R_{⊥}$ и $K_{⊥}$ , которые являются свободными параметрами при выборе начальных условий для лучевых уравнений. Воспользовавшись дисперсионным уравнением (28), можно определить начальный импульс $K_{z 0} (R_{⊥}, K_{⊥})$ при $z = 0$ и всех возможных значениях начальных поперечных координат и импульсов. Таким образом мы задаем полный набор начальных координат и импульсов в шестимерном пространстве, необходимый для расчета характеристик в соответствии с лучевыми уравнениями (29) или (34).

Очевидно, что начальная функция Вигнера $W_{0} (R_{⊥}, K_{⊥})$ в (31) также зависит только от четырех свободных параметров – компонент поперечных координат и импульсов. При этом истинная функция Вигнера в шестимерном фазовом пространстве на плоскости $z = 0$ должна зависеть от пятого параметра $-$ продольного импульса $K_{z}$ . Эта зависимость определяется дисперсионным соотношением, которое выражено в виде множителя $δ (H')$ в (32). В итоге функция Вигнера при $z = 0$ представляется в виде

$\begin{matrix} W^{(z = 0)} (R_{⊥}, K) = {|\frac{\partial H^{'}}{\partial K_{z}}|}^{- 1} \\ δ (K_{z} - K_{z 0} (R_{⊥}, K_{⊥})) W_{0} (R_{⊥}, K_{⊥}) . \end{matrix}$ (35)

Для того, чтобы определить $W_{0}$ , воспользуемся соотношением (7) для заданного при $z = 0$ распределения комплексной амплитуды $E_{0} = e u$ или, эквивалентно, $u = e^{*} E_{0}$ :

$u (r_{⊥}) u^{*} (r'_{⊥}) =$

$= \int W^{(z = 0)} (\frac{r_{⊥} + r'_{⊥}}{2}, K) e^{i K_{⊥} (r_{⊥} - r'_{⊥})} d K =$

$= \int {|\frac{\partial H'}{\partial K_{z}}|}^{- 1} W_{0} (\frac{r_{⊥} + r'_{⊥}}{2}, K_{⊥}) e^{i K_{⊥} (r_{⊥} - r'_{⊥})} d K_{⊥} .$

Выполняя обратное преобразование Фурье, получаем

$W_{0} (R_{⊥}, K_{⊥}) = {|\frac{\partial H^{'} (R_{⊥}, K_{⊥})}{\partial K_{z}}|}_{z = 0} \times$ (36)

$\times \int u (R_{⊥} + \frac{ξ_{⊥}}{2}) u^{*} (R_{⊥} - \frac{ξ_{⊥}}{2}) e^{- i K_{⊥} ξ_{⊥}} \frac{d ξ_{⊥}}{{(2 π)}^{2}} .$

Итак, мы пришли к физически понятному результату – шестимерная функция Вигнера для полного поля (35) определяется через четырехмерную начальную функцию Вигнера (36), вычисленную по заданному двумерному распределению комплексной амплитуды на границе области моделирования.

В качестве важного практического примера рассмотрим волновое поле с обобщенным гауссовским распределением

$\begin{matrix} u (x, y) = \frac{1}{\sqrt{π a_{x} a_{y}}} exp (- \frac{x^{2}}{2 a_{x}^{2}} - \frac{y^{2}}{2 a_{y}^{2}} - \frac{i k_{0} x^{2}}{2 R_{x}} - \\ - \frac{i k_{0} y^{2}}{2 R_{x}} - \frac{i k_{0} x y}{R_{x y}} + i k_{x} x + i k_{y} y), \end{matrix}$

заданное в вакууме,

$H' = K_{x}^{2} + K_{y}^{2} + K_{z}^{2} - k_{0}^{2} .$

Отвечающая этим условиям начальная функция Вигнера (35) есть

$W^{(z = 0)} (X, Y, K_{x}, K_{y}, K_{z}) =$

$= \frac{1}{π^{2}} δ (K_{z} - \sqrt{k_{0}^{2} - K_{x}^{2} - K_{y}^{2}}) \times$

$\times exp [- \frac{X^{2}}{a_{x}^{2}} - a_{x}^{2} {(K_{x} - k_{x} + \frac{k_{0} X}{R_{x}} + \frac{k_{0} Y}{R_{x y}})}^{2} -$

$- \frac{Y^{2}}{a_{y}^{2}} - a_{y}^{2} {(K_{y} - k_{y} + \frac{k_{0} Y}{R_{y}} + \frac{k_{0} X}{R_{x y}})}^{2}] .$

9. Плотность потока энергии

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

Рассмотрим сначала случай среды без пространственной дисперсии, в которой средняя плотность потока энергии равна вектору Пойнтинга, см., например, формулу (9.5) § 9 гл. 2 [37]:

$S_{0} = \frac{c}{8 π} Re [E \times B^{*}] .$

Учитывая связь между электрическим и магнитным полем, $B^{*} = i k_{0}^{- 1} [\nabla \times E^{*}]$ , следующую из уравнений Максвелла для $exp (- i ω t)$ -процессов, получаем

$S_{0 γ} = \frac{c}{8 π k_{0}} \underset{α}{} Re (i E_{α} \frac{\partial E_{α}^{*}}{\partial x_{γ}} - i E_{α} \frac{\partial E_{γ}^{*}}{\partial x_{α}}),$

где $S_{0 γ}$ – проекция вектора Пойнтинга на декартову ось $γ$ . Воспользовавшись (12), можно выразить квадратичные по полю величины через тензор Вигнера:

$\begin{matrix} S_{0 γ} = \frac{c}{16 π k_{0}} \underset{α}{} \int \{[2 K_{γ} W_{α α} - K_{α} (W_{α γ} + W_{γ α})] - \\ - i \frac{\partial}{\partial x_{α}} (W_{γ α} - W_{α γ})\} d K . \end{matrix}$ (37)

Подынтегральную часть в квадратных скобках можно интерпретировать как суммирование по вкладам от всех плоских волн. Оставшееся слагаемое с производными от тензора Вигнера описывает дополнительный поток энергии, возникающий, например, в пучке циркулярно поляризованных волн из-за неоднородного распределения момента импульса поля или за счет линейного взаимодействия нормальных мод. Этот вклад важен для интерпретации тонких эффектов, таких как смещение Федорова–Имберта [51]. Однако в квазиоптическом приближении члены с производными от тензора Вигнера появляются в $S_{0}$ в следующем порядке и ими можно пренебречь. С этой же точностью тензор Вигнера $W_{γ α}$ определяется по скалярной функции Вигнера W с помощью соотношения (18).

Известно, что в среде с пространственной дисперсией плотность потока энергии $S = S_{0} + T$ отличается от вектора Пойнтинга ([37], $§$ 9 гл. 2). Связанный с пространственной дисперсией дополнительный член $T$ в плотности потока энергии довольно изящным способом можно связать с тензорной функцией Вигнера. Этот вывод для случая плавнонеоднородной среды приведен в Приложении 2. Поправка к плотности потока энергии может быть выражена следующим образом:

$T_{γ} = - \frac{ω}{16 π} \sum_{α, β} \int \frac{\partial {ε^{'}}_{α β}}{\partial K_{γ}} W_{β α} d K,$ (38)

где ${ε^{'}}_{α β} (R, K)$ – эрмитова часть тензора вейлевского представления диэлектрического оператора . Комбинируя (37), в котором отбросили последнее слагаемое, и (38), получаем плотность потока энергии в виде

$S_{γ} = \frac{c}{16 π k_{0}} \sum_{α, β} \int \{[2 K_{γ} W_{α α} - K_{α} (W_{α γ} + W_{γ α})] -$

$- k_{0}^{2} \frac{\partial ε'_{α β}}{\partial K_{γ}} W_{β α}\} d K = \frac{c}{16 π k_{0}} \sum_{α, β} \int \frac{\partial {D^{'}}_{α β}}{\partial K_{γ}} W_{β α} d K .$

Выражение в фигурных скобках удалось упростить, используя соотношение (14) для эрмитовой части $D'_{α β}$ вейлевского представления волнового оператора. Подставим под интеграл тензорную функцию Вигнера в виде (18), воспользуемся дисперсионным соотношением (28) и следующим из него соотношением $\sum_{β} D'_{α β} e_{β} = 0$ $ε_{α β}$ , получим

$\underset{}{\sum_{α, β}} \frac{\partial {D^{'}}_{α β}}{\partial K_{γ}} e_{β} e_{α}^{*} W = \sum_{α, β} \frac{\partial}{\partial K_{γ}} (e_{α}^{*} e_{β} D'_{α β}) W = \frac{\partial H'}{\partial K_{γ}} W,$

где $H^{'}$ – действительная часть вейлевского представления оператора Гамильтона (17). В итоге получаем

$S_{γ} (r) = \frac{c}{16 π k_{0}} \int \frac{\partial H' (r, K)}{\partial K_{γ}} W (r, K) d K .$ (39)

Таким образом, мы приходим к универсальному выражению, не зависящему от наличия или отсутствия пространственной дисперсии среды. Вычисление декартовой компоненты вектора плотности потока энергии сводится к вычислению свертки скалярной функции Вигнера с производной от действительной части гамильтониана по соответствующей компоненте волнового вектора. Очевидно, это равенство можно трактовать и как векторное, связывающее S и $\partial H' / \partial K$ .

Подставим теперь в выражение (39) полученное ранее решение для функции Вигнера. Прежде чем переходить к общему случаю, рассмотрим поток энергии через плоскость z = 0, в которой задаются начальные условия. Функция Вигнера в этом случае определяется выражением (35), а нормальная компонента плотности потока энергии есть

$S_{z} (R_{⊥}) = \int w d K_{⊥}, w = \frac{c}{16 π k_{0}} W_{0} (R_{⊥}, K_{⊥}) .$ (40)

При получении этой формулы производная от гамильтониана сократилась, при этом мы учли, что $\partial H' / \partial K_{z} > 0$ . Физически последнее условие означает, что все характеристики (лучи) выпускаются в “положительном” направлении оси $z$ . Полный поток энергии через плоскость $z = 0$ равен

$P_{0} = \int w d R_{⊥} d K_{⊥} .$

К общему случаю можно перейти с помощью наглядных физических соображений. Для этого заметим, что величина $w (R_{⊥}, K)$ имеет смысл доли полной мощности, приходящейся на один луч с начальными координатами $R_{⊥}$ и импульсами 0 $K_{⊥}$ . С другой стороны, диссипация этой мощности может быть описана как затухание вдоль луча в соответствии с уравнением (30) и его решением (31). Поэтому полный поток энергии, протекающий через некоторую произвольную ориентированную поверхность Г с нормалью $n$ , равен

$P = \int s i g n (\frac{\partial H^{'}}{\partial K_{n}}) \exp (\int_{0}^{τ^{*}} 2 H^{''} d τ) w d R_{⊥} d K_{⊥} .$ (41)

Интеграл берется по всем начальным поперечным координатам и импульсам, для которых характеристика пересекает поверхность Г, функция начальных условий $τ^{*} (R_{⊥}, K_{⊥})$ определяет значение эволюционной переменной, при котором луч пересекает поверхность Г, а экспонента определяет коэффициент ослабления интенсивности для каждого луча. В общем случае лучи могут пересекать поверхность с разных сторон – это учитывает дополнительный множитель $sign (\partial H^{'} / \partial K_{n})$ , который должен вычисляться в точке на поверхности; здесь $K_{n}$ – компонента импульса вдоль нормали к поверхности.

К такому же результату можно прийти и более формальным и математически строгим путем. Для этого рассмотрим поток энергии через элементарную площадку dR_Г = dR₁dR₂ поверхности Г с нормалью n. В рассматриваемой точке на поверхности разложим импульс на продольную и поперечную составляющие по отношению к нормали,

K = K_nn + K_Г,

где К_Г = (К₁, К₂) есть двумерный волновой вектор, отвечающий локальным координатам R_Г = (R₁, R₂). Поток энергии через элементарную площадку определяется из (39) как

$S_{n} d R_{Г} = \frac{c}{16 π k_{0}} \int \frac{\partial H^{'}}{\partial K_{n}} W (R_{Г} K) d K_{n} d K_{Г} d R_{Г} .$

Подставим сюда наше общее решение (32) для функции Вигнера. Интеграл по $K_{n}$ можно взять за счет $δ (H^{'})$ , при этом производные гамильтониана сокращаются (остается только знак), а оставшаяся часть функции Вигнера выражается через формулу (31). В итоге c учетом определения w (40) получаем

$S_{n} d R_{Г} = \int s i g n (\frac{\partial H^{'}}{\partial K_{n}}) \exp (\int_{0}^{τ *} 2 H^{''} d τ) w d K_{Г} d R_{Г}^{} .$ (42)

Здесь τ* по прежнему есть значение эволюционной переменной, при котором луч пересекает поверхность.

Неудобство выражения (42) связано с тем, что интегрирование в нем подразумевается по конечным точкам пути R_ГK_Г, в то время как удельный вес w и точка остановки $τ^{*}$ лучей определяются начальными условиями $R_{⊥}, K_{⊥}$ на плоскости $z = 0$ . Чтобы избежать необходимости решать обратную задачу вдоль характеристик, можно воспользоваться теорией интегральных инвариантов гамильтоновых систем ([50], § 122), позволяющей перейти к интегрированию по начальным условиям. Для этого воспользуемся теоремой о понижении порядка автономной системы уравнений Гамильтона, которая утверждает, что поскольку величина $τ$ явно не входит в гамильтониан $H$ , то эволюция величин R_Ги K_Г будет также описываться с помощью некоторой системы уравнений Гамильтона^[1]. Для этой новой системы воспользуемся теоремой Лиувилля, утверждающей о сохранении фазового объема при эволюции во времени для любой гамильтоновой системы. В нашем случае это сохранение четырехмерного фазового объема

$d R_{Г} d К_{Г} = d R_{⊥} d К_{⊥} .$

В результате в (42) можно перейти к интегрированию по начальным поперечным координатам и импульсам, после чего это выражение становится полностью эквивалентным (41).

Аналогичным образом может быть построена плотность мощности, поглощаемой в объеме $d V$ в окрестности точки $R$ . Для этого достаточно заметить, что эта мощность определяется правой частью (30) и заменить интегрирование по лучам, которые пересекают поверхность Г, на интегрирование по лучам, которые приходят в точку $R$ . В результате можно получить

$\frac{d P (r)}{d V} = - \int \int_{0}^{\infty} 2 H^{''} (τ) \exp (\int_{0}^{τ} 2 H^{''} (τ^{'}) d τ^{'}) \times$

$\times δ (r - R (τ, R_{⊥}, K_{⊥})) w d τ d R_{⊥} d K_{⊥} .$ (43)

На практике численное вычисление интегралов такого вида удобно реализовать в виде подсчета гистограммы – фазовое пространство разбивается на ячейки, выпускается набор лучей с разными начальными условиями так, чтобы покрыть всю интересующую нас область фазового пространства, и подсчитывается количество попаданий лучей в каждую ячейку фазового пространства с учетом набранного веса. В данном случае весом выступает плотность мощности вдоль луча $d W / d τ$ .

Указанные рассуждения сохраняют силу и для стохастических уравнений Гамильтона (34), если стохастическая сила сохраняет гамильтониан и все манипуляции с фазовыми объемами проводятся для одной “замороженной” реализации случайной силы G. Для этой реализации необходимо посчитать интегралы (41) и (43) по всем начальным условиям, и только потом переходить к следующей реализации для набора статистики.

В дополнение к сравнению метода функций Вигнера и метода эволюционного квазиоптического уравнения, приведенному в разделе 5, отметим, что первый метод значительно более экономичен к памяти компьютера. Действительно, при вычислении интегралов вида (41) или (43) методом гистограмм размер необходимой памяти определяется только размерностью левой части, у нас нет необходимости запоминать координаты лучей в шестимерном пространстве. Высокая размерность задачи проявляется только в пространстве начальных условий (четырехмерном), что сказывается на времени вычислений, но не на требованиях к памяти. При вычислении гистограмм начальные условия можно перебирать один за другим последовательно в любом порядке, кроме того, эта задача естественным образом распараллеливается. С другой стороны, при реализации квазиоптического расчета, основанного на уравнении (23), на каждом шаге по эволюционной координате t необходимо вычислять псевдодифференциальный оператор ${\hat{H}}_{QO}$ . В общем случае эта операция в четырехмерном пространстве, отвечающем двум координатам поперек опорного луча и сопряженным им каноническим импульсам. Поэтому, если для хранения $\bar{u} (τ, r_{⊥})$ для заданного t использовать двумерный массив с числом элементов $n \times n$ , то для вычисления ${\hat{H}}_{QO} [\bar{u}]$ потребуется хранить в памяти четырехмерный массив с числом элементов $O (n^{4})$ .

10. Пример

В качестве иллюстрации применения описанной техники расчета рассмотрим схему ЭЦ нагрева плазмы в токамаке Т-15МД на третьей гармонике гирочастоты, предложенную в [9]. В этой схеме предполагается вертикальный ввод излучения с частотой 140 ГГц и поляризацией, отвечающей необыкновенной волне, через верхний порт токамака. Магнитное поле на оси токамака при этом составляет 1.67 Тл, что соответствует ЭЦ резонансу на третьей гармонике. В статье [9] были проведены довольно грубые оценки на основе геометрической оптики, показавшие, что при характерных ожидаемых параметрах разряда в токамаке Т-15МД и углах ввода изучения в районе 3 $-$ 10° к вертикали можно добиться поглощения более 99% от вводимой мощности СВЧ-излучения. Приведенный ниже анализ приводит к более скромным оценкам поглощаемой мощности.

Для расчета ЭЦ поглощения методом кинетического уравнения использовался лучевой гамильтониан в слаборелятивистской плазме [12]. Для каждого из рассмотренных сценариев был проведен расчет:

· пространственного распределения интенсивности поля, ${|u (r)|}^{2} = \int W (r, K) d K;$

· профиля поглощения СВЧ-мощности, определяемого аналогично (43), $Π (ρ) = P_{0}^{- 1} d P / d ρ$ ,

где

$\frac{d P}{d ρ} = - \int \int_{0}^{\infty} 2 H^{″} (τ) exp (\int_{0}^{τ} 2 H^{″} (τ^{'}) d τ^{'}) \times$

$\times δ (ρ - ρ (R (τ, R_{⊥}, K_{⊥}))) w d τ d R_{⊥} d K_{⊥};$

· доли поглощенной мощности, $Q = \int Π (ρ) d ρ$ .

Для нахождения функции Вигнера в шестимерном фазовом пространстве были использованы 810 000 лучевых траекторий, отвечающей сетке $30 \times 30 \times 30 \times 30$ в черырехмерном пространстве начальных условий. Расчеты были выполнены двумя способами $-$ с помощью кинетического уравнения (33), учитывающего диффузионные поправки во втором приближении, и с помощью уравнения первого порядка (26) без учета диффузионных членов.

Для моделирования были выбраны следующие параметры на магнитной оси системы: электронная температура 3 кэВ, плотность плазмы $4 \cdot 10^{13}$ см^–3, магнитное поле 1.67 Тл. Используемые в расчетах профили температуры и плотности плазмы приведены на рис. 1. Метки магнитной поверхности вычислялись как квадратный корень из нормированного полоидального магнитного потока $ψ$ :

$ρ (r) = \sqrt{(ψ - ψ_{axis}) / (ψ_{sep} - ψ_{axis})} .$

Рис. 1. Профили электронной температуры (а) и плотности плазмы (б), использованные в моделировании.

Магнитная конфигурация была аналогична [9]. Начальный СВЧ пучок соответствовал системе ввода излучения токамака Т-15МД. Его полуширина в области перетяжки в вакууме по уровню интенсивности $1 / e$ составляла 2.2 см, на входе она составляла 3.12 см при радиусе кривизны фазового фронта в 284 см. Точка ввода пучка (координаты центра пучка) задавалась на высоте 142 см от экваториальной плоскости.

Были рассчитаны два сценария ЭЦ-нагрева, в обоих случаях вводимый СВЧ-пучок в вакууме нацеливался на магнитную ось системы. В первом сценарии угол ввода пучка составлял 5° к вертикали, во втором $-$ 10° к вертикали, при этом мы смещали горизонтальную координату точки ввода излучения.

Результаты моделирования приведены на рис. 2. Эффективность поглощения для ввода СВЧ-пучка под углом 5° к вертикали составляет 32%, при этом энерговклад локализован в окрестности магнитной оси. Моделирование показало важность учета диффузионных членов в кинетическом уравнении $-$ эффективность поглощения, рассчитанная без учета этих членов, составила 21%. Отметим, что эффективность поглощения в данном случае в основном ограничивается отражением от области резонанса. Для ввода под углом 10° к вертикали поглощение становится более эффективным, несмотря на то, что из-за рефракции точка резонанса для пучка смещается примерно на 60 см ниже магнитной оси системы. Доля поглощенной мощности при этом составляет 88%. Расчет без эффектов лучевой диффузии показал завышенную эффективность поглощения на уровне 93%

Рис. 2. Моделирование распространения и поглощения квазиоптических волновых пучков в плазме токамака: сценарий ЭЦ-нагрева на третьей гармонике с вводом излучения под углом 5° к вертикали (верхний ряд) и 10° к вертикали (нижний ряд). Приведены распределения интенсивности волнового поля | u | ²в полоидальном сечении токамака, рассчитанные без учета диффузионных членов в первом порядке теории возмущений (а, г) и с учетом диффузионных членов во втором порядке теории возмущений (б, д). Сплошной красной линией выделена область поглощения СВЧ-мощности, определенная по уровню max (dP / dV) / 2 ; пунктирная линия обозначает границу плазмы. На графиках (в, е) приведены профили поглощенной мощности П(ρ) с учетом (2) и без учета (1) диффузионных членов.

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

11. Заключение

В данной работе приведено систематическое изложение метода кинетического уравнения для функции Вигнера для моделированию распространения и поглощения электромагнитных волн в анизотропных и гиротропных диссипативных средах с пространственной дисперсией. Для случая плавнонеоднородных и слабодиссипативных сред, наиболее востребованного для практических приложений, построена асимптотическая теория. В первом порядке эта теория приводит к кинетическому уравнению Больцмана, решения которого можно строить в виде суперпозиции характеристик $-$ невзаимодействующих лучей в шестимерном фазовом пространстве. Кинетическое уравнение для лучей (формализм Вейля) было известно и раньше, однако в нашей работе мы предложили физическое обоснование ввода скалярной функции Вигнера для описания векторной волновой задачи. Во втором порядке кинетическое уравнение Больцмана переходит к уравнению типа Фоккера $-$ Планка, причем диффузионные члены появляются за счет неоднородного поглощения. Это означает, что между лучами возникает “диффузия” интенсивности, которая, в свою очередь, может быть эффективно смоделирована с помощью стохастических членов в уравнениях Гамильтона для отдельных лучей. Это позволяет учесть диффузионные эффекты непосредственно в кинетическом лучевом уравнении без заметного усложнения численной схемы его решения.

Авторы пришли к данному методу при попытке дать формальное обоснование последовательному (аберрационному) квазиоптическому подходу и, в более широком смысле, процедуре восстановления диэлектрического отклика неоднородной диспергирующей среды по его геометрооптическому пределу. Метод кинетического уравнения позволил продвинуться в данном направлении. Кроме того, он обладает некоторыми дополнительными преимуществами по сравнению с квазиоптическим подходом, основанным на укорачивании волнового уравнения на комплексную амплитуду геометрооптической моды, такими как более прозрачный учет неоднородной диссипации, отсутствие ограничений на кривизну лучей, очень экономное использование памяти при численном решении. В работе этот вопрос не рассматривается, но представляется вполне очевидным, что дальнейшее обобщение для учета линейного взаимодействия мод в рамках кинетического подхода проводится более просто и строго, чем в рамках обычного квазиоптического подхода. Действительно, уравнения (13) эквивалентны уравнениям Максвелла, поэтому описывают полное векторное поле; переход к взаимодействующим модам соответствует диагонализации $W_{β γ}$ в этих уравнениях, которая может быть выполнена приближенно. Тем не менее оба подхода, кинетический и квазиоптический, активно развиваются, какой из них лучше подходит для решения актуальных задач электродинамики высокотемпературной плазмы, покажет время.

Приведенный пример расчета сценария ЭЦ нагрева в большом токамаке показал, что неоднородность поглощения поперек квазиоптического волнового пучка может приводить к систематической ошибке при использовании кинетического уравнения первого порядка (для невзаимодействующих лучей) по сравнению с результатами использования кинетического уравнения второго порядка (диффузионное приближение). В этом предварительном примере мы хотели проиллюстрировать новые эффекты, связанные с неоднородным поглощением, которые не описываются обычным кинетическим уравнением для лучей первого порядка, как, впрочем, и остальными стандартными методами, упомянутыми во введении. Более детальное описание созданного для этой демонстрации нового кода, его верификация и результаты моделирования сценариев ЭЦ-нагрева плазмы в перспективных установках будут опубликованы отдельно. Заметим только, что вклад, связанный с поперечной неоднородностью поглощения, может оказаться заметным для тангенциальной или вертикальной инжекции волновых пучков в тороидальную ловушку [7, 13, 14], а также поперечного ввода излучения в аксиально-симметричную магнитную ловушку [11].

Работа выполнена в рамках государственного задания ИПФ РАН (тема FFUF-2023-0002) и ФТИ РАН (по договору между ФТИ РАН и ИПФ РАН $№$ 44-6-24 от 29.07.2024).

Приложение 1. Генерация случайного вектора с заданной ковариацией компонент

Для получения случайных величин $G_{i}$ с заданным тензором дисперсии можно адаптировать метод Якоби диагонализации квадратичной формы [52]. При этом для реализации Монте-Карло алгоритма нам не важна конкретная функция распределения случайных величин, нам нужны только первые два момента и ограниченность высших моментов. Приведенная ниже процедура приводит набор из $N$ случайных независимых величин $θ_{i}$ , распределенных по какому-то, не обязательно нормальному, закону с нулевым средним и единичной дисперсией,

$⟨ θ_{i} ⟩ = 0, ⟨ θ_{i} θ ⟩ = δ_{i j},$

к набору из $N$ случайных зависимых величин $G_{i}$ с нулевым средним и положительно определенным тензором дисперсии

$⟨ G_{i} ⟩ = 0, ⟨ G_{i} G ⟩ = F_{i j} .$

Процедура следующая:

$G_{1} = α_{11} θ_{1},$

$α_{11} = \sqrt{F_{11}},$

$G_{2} = α_{21} θ_{1} + α_{22} θ_{2},$

$α_{21} = F_{21} / α_{11},$

$α_{22} = \sqrt{F_{22} - α_{21}^{2}},$

$G_{3} = α_{31} θ_{1} + α_{32} θ_{2} + α_{33} θ_{3},$

$α_{31} = F_{31} / α_{11},$

$α_{32} = (F_{32} - α_{31} α_{21}) / α_{22},$

$α_{33} = \sqrt{F_{33} - α_{31}^{2} - α_{32}^{2}},$

$G_{4} = α_{41} θ_{1} + α_{42} θ_{2} + α_{43} θ_{3} + α_{44} θ_{4},$

$α_{41} = F_{41} / α_{11},$

$α_{42} = (F_{42} - α_{41} α_{21}) / α_{22},$

$α_{43} = (F_{43} - α_{41} α_{31} - α_{42} α_{32}) / α_{33},$

$α_{44} = \sqrt{F_{44} - α_{41}^{2} - α_{42}^{2} - α_{43}^{2}}$

$\dots$

При этом положительная определенность матрицы $F_{m n}$ гарантирует нам действительность всех коэффициентов, в частности $α_{i i}^{2} > 0$ .

Приложение 2. Представление дополнительного потока энергии в диспергирующей среде через тензор Вигнера

Средняя плотность мощности джоулевых потерь для монохроматического поля в произвольной среде равна

$P (r) = \frac{1}{2} Re \sum_{α, β} E_{α}^{*} σ_{α β} E_{β} = \frac{1}{2} Re P (r, r),$

где $σ_{α β}$ $-$ линейный оператор тензора проводимости среды с ядром $σ_{α β} (r, r')$ , а

$P (r, r') = \frac{1}{2} \sum_{α, β} \int σ_{α β} (r, r'') E_{α}^{*} (r') E_{β} (r'') d r'',$

можно рассматривать как интегральное ядро некоторого формального оператора $P = \sum σ_{α β} W_{β α} .$ Этот оператор является суперпозицией оператора проводимости и оператор Вигнера для электромагнитного поля. Поэтому его вейлевское представление определяется с помощью (5) как

$P (R, K) = exp (\frac{i L}{2}) \sum_{α, β} σ_{α β} (R, K) W_{β α} (R, K),$ (44)

где $W_{α β}$ $-$ тензорная функция Вигнера, $σ_{α β}$ $-$ вейлевское представление тензора проводимости,

$L = \frac{\partial_{σ}}{\partial R} \frac{\partial_{W}}{\partial K} - \frac{\partial_{W}}{\partial R} \frac{\partial_{σ}}{\partial K} = \frac{\partial}{\partial K} \frac{\partial_{σ}}{\partial R} - \frac{\partial}{\partial R} \frac{\partial_{σ}}{\partial K} .$

Восстанавливая ядро оператора по его вейлевскому представлению в соответствии с формулой (4), находим мощность джоулевых потерь

$P (r) = \frac{1}{2} Re \int P (r, K) d K .$

В плавнонеоднородной среде мы можем разложить экспоненту в (44) ряд до первого порядка. Далее, разобьем матрицу вейлевского представления оператора проводимости на эрмитову и антиэрмитову части, $σ_{α β} = σ'_{α β} + i σ {''}_{α β},$ и учтем, что $W_{α β}$ $-$ эрмитова матрица. В результате получим

$P = \sum_{α, β} \int \{\frac{1}{2} σ'_{α β} W_{β α} - \frac{1}{4} D σ {''}_{α β} W_{β α}\} d K .$

Первое слагаемое представляет собой стандартное выражение для плотности мощности джоулевых потерь (истинное поглощение Q). Второе слагаемое появляется только в средах с пространственной дисперсией. При его вычислении часть $\partial / \partial K$ оператора $L$ не дает вклад в интеграл по $K$ , а часть $\partial / \partial R$ всегда представляется в виде дивергенции от некоторого вектора. Поэтому

$P = Q + div T, T = \sum_{α, β} \int \frac{1}{4} \frac{\partial σ {''}_{α β}}{\partial K} W_{β α} d K .$

Очевидно, последняя часть должна быть включена в выражение для плотности потока энергии. Выразив оператор проводимости через тензор диэлектрической проницаемости, $σ {''}_{α β} = - ω ε'_{α β} / 4 π$ , получим окончательную форму:

$T (r) = - \frac{ω}{16 π} \sum_{α, β} \int \frac{\partial ε'_{α β} (r, K)}{\partial K} W_{β α} (r, K) d K,$

Эта формула обобщает известное выражение для дополнительного потока энергии в однородной среде с пространственной дисперсией, смотри, например, формулу (9.18) в $§$ 9 гл. 2 [37]. Для того, что бы получить формулу (9.18), надо предположить, что пространственный спектр поля достаточно узкий. Тогда $\partial ε'_{α β} / \partial K_{γ}$ можно вынести из под интегрирования и получить

$T \approx - \frac{ω}{16 π} \sum_{α, β} \frac{\partial ε'_{α β}}{\partial K} E_{α}^{*} E_{β} .$

Наша техника также естественно обобщается для учета следующих поправок, проявляющихся в не слишком плавно неоднородных средах.

^[1]Понять это утверждение можно на следующем примере. Из закона сохранения гамильтониана исходной автономной системы можно выразить импульс $K_{z} (. . .)$ как функцию всех остальных переменных. Выберем z в качестве эволюционной координаты. Тогда эволюция ортогональных координат x, y и сопряженных им импульсов может рассматриваться как новая неавтономная система с функцией Гамильтона $K_{z} (. . .)$ с $τ$ на месте z. Фазовый объем будет сохраняться на поверхностях z = const. Если гамильтониан укороченной системы умножить на произвольную положительную функцию координат и импульсов, то получим сохранение фазового объема на системе поверхностей, не конгруэнтных z = const.

About the authors

Е. D. Gospodchikov

Gaponov-Grekhov Institute of Applied Physics, Russian Academy of Sciences

Author for correspondence.
Email: egos@ipfran.ru
Russian Federation, Nizhny Novgorod

A. A. Balakin

Gaponov-Grekhov Institute of Applied Physics, Russian Academy of Sciences

Email: egos@ipfran.ru
Russian Federation, Nizhny Novgorod

A. G. Shalashov

Gaponov-Grekhov Institute of Applied Physics, Russian Academy of Sciences

Email: egos@ipfran.ru
Russian Federation, Nizhny Novgorod

References

Аликаев В.В., Литвак А.Г., Суворов Е.В., Фрайман А.А. // Высокочастотный нагрев плазмы. Материалы всесоюзного совещания. Горький: ИПФ АН, 1983.
Prater R., Farina D., Gribov Yu., Harvey R.W., Ram A.K., Lin-Liu Y.-R., Poli E., Smirnov A.P., Volpe F., Westerhof E., Zvonkov A. and the ITPA Steady State Operation Topical Group // Nuclear Fusion. 2008. V. 48. P. 035006.
Pereverzev G.V. // Physics of Plasmas. 1998. V. 5. P. 3529.
Litvak A.G., Permitin G.V., Suvorov E.V., Fraiman A.A. // Nuclear Fusion. 1977. V. 17. P. 659.
Romé M., Erckmann V., Gasparino U., Karulin N. // Plasma Phys. Control. Fusion. 1998. V. 40. P. 511.
Wagner D., Stober J., Leuterer F., Monaco F., M ü ller S., M ü nich M., Rapson C.J., Reich M., Schubert M., Sch ü tz H., Treutterer W., Zohm H., Thumm M., Scherer T., Meier A., Gantenbein G., Jelonnek J., Kasparek W., Lechte C., Plaum B., Goodman T., Litvak A.G., Denisov G.G., Chirkov A., Zapevalov V., Malygin V., Popov L.G., Nichiporenko V.O., Myasnikov V.E., Tai E.M., Solyanova E.A., Malygin S.A. and ASDEX Upgrade Team // J. Infrared Millimeter Teraherz Waves. 2016. V. 37. P. 45.
Балакин А.А., Балакина М.А., Пермитин Г.В., Смирнов А.И. // Физика плазмы. 2008. Т. 34. С. 533.
Porte L., Coda S., Alberti S., Arnoux G., Blanchard P., Bortolon A., Fasoli A., Goodman T.P., Klimanov Y., Martin Y., Maslov M., Scarabosio A., Weisen H. // Nuclear Fusion. 2007. V. 47. P. 952.
Кирнева Н.А., Борщеговский А.А., Куянов А.Ю., Пименов И.С., Рой И.Н. // ВАНТ Сер. Термоядерный синтез. 2021. Т. 44. С. 24.
Shalashov A.G., Gospodchikov E.D., Smolyakova O.B., Bagryansky P.A., Malygin V.I., Thumm M. // Phys. Plasmas. 2012. V. 19. P. 052503.
Сковородин Д.И., Черноштанов И.С., Амиров В.Х., Астрелин В.Т., Багрянский П.А., Беклемишев А.Д., Бурдаков А.В., Горбовский А.И., Котельников И.А., Магоммедов Э.М., Полосаткин С.В., Поступаев В.В., Приходько В.В., Савкин В.Я., Солдаткина Е.И., Соломахин А.Л., Сорокин А.В., Судников А.В., Христо М.С., Шиянков С.В., Яковлев Д.В., Щербаков В.И. // Физика плазмы. 2023. Т. 49. С. 831.
Балакина М.А., Смолякова О.Б., Токман М.Д. // Физика плазмы. 2003. Т. 29. С. 60.
Балакин А.А., Балакина М.А., Шашалов А.Г. // Физика плазмы. 2007. Т. 33. С. 724.
Шалашов А.Г., Господчиков Е.Д. // УФН. 2022. Т. 192. С. 1399.
Балакин А.А., Балакина М.А., Пермитин Г.В., Смирнов А.И. // Физика плазмы. 2007. Т. 33. С. 337.
Balakin A.A., Balakina M.A., Permitin G.V., Smirnov A.I. // J. Phys. D: Appl. Phys. 2007. V. 40. P. 4285.
Балакин А. А., Господчиков Е. Д., Шалашов А. Г. Письма ЖЭТФ. 2016. Т. 104. С. 701.
Shalashov A.G., Balakin A.A., Gospodchikov E.D., Khusainov T.A. // Phys. Plasmas. 2016. V. 23. P. 112504.
Шалашов А. Г., Балакин А. А., Хусаинов Т. А., Господчиков Е. Д., Соломахин А. Л. // ЖЭТФ. 2017. Т. 151. С. 379.
Maj O., Balakin A.A., Poli E. // Plasma Phys. Controlled Fusion. 2010. V. 52. P. 085006.
Tereshchenko M., Castejón F., Cappa A. // Plasma Phys. Controlled Fusion. 2013. V. 55. P. 115011.
Castejón F., Cappa A., Tereshchenko M., Fernandez A. // Nuclear Fusion. 2008. V. 48. P. 075011.
McDonald S.W. // Phys. Rev. A. 1991. V. 43. P. 4484.
McDonald S.W. // Phys. Rep. 1988. V. 158. P. 337.
Wigner E. // Phys. Rev. 1932. V. 40. P. 749.
Zachos C., Fairlie D., Curtright T. Quantum mechanics in phase space: an overview with selected papers. London: World Scientific, 2005.
Marushchenko N.B., Turkin Y., Maassberg H. // Computer Physics Communications. 2014. V. 185. P. 165.
Farina D.A // Fusion Sci. Technol. 2007. V. 52. P. 154.
Балакин А.А. // Изв. вузов: Радиофизика. 2012. Т. 55. С. 617.
Köhn A., Guidi L., Holzhauer E., Maj O., Poli E., Snicker A., Weber H. // Plasma Phys. Control. Fusion. 2018. V. 60. P. 075006.
Guidi L., Maj O., Weber H., K ö hn A., Snicker A., Poli E. // J. Phys: Confer. Ser. 2016. V. 775. P. 012005.
Snicker A., Poli E., Maj O., Guidi L., K ö hn A., Weber H., Conway G., Henderson M., Saibene G. // Nuclear Fusion. 2017. V. 58. P. 016002.
Weber H., Maj O., Poli E. // J. Phys.: Confer. Ser. 2018. V. 1125. P. 012022.
Weber H., Maj O., Poli E. // J. Computational Electronics. 2021. V. 20. P. 2199.
Weber H., Maj O., Poli E. // EPJ Web Confere. EDP Sci. 2023. V. 277. P. 01003.
Poli E., Bottino A., Korger D., Maj O., Palermo F., Weber H. // New J. Phys. 2024. V. 26. P. 013016.
Brambilla M. Kinetic theory of plasma waves: homogeneous plasmas. Oxford University Press, 1998.
Weyl H. // Zeitschrift f ü r Physik. 1927. V. 46. P. 1.
Шалашов А.Г., Господчиков Е.Д. // УФН. 2012. Т. 182. С. 157.
Балакин А.А. // Изв. высших учебных заведений. Радиофизика. 2012. Т. 55. С. 521.
Кравцов Ю.А., Орлов Ю.И. Геометрическая оптика неоднородных сред. М.: Наука, 1980.
Ziolkowski R.W., Deschamps G.A. // Radio Science. 1984. V. 19. P. 1001.
Терещенко М.А. // Физика плазмы. 2017. Т. 43. С. 9.
Кляцкин В.И. Метод погружения в теории распространения волн. М.: Наука, 1986.
Bornatici M., Engelmann F., Maroli C., Petrillo V. // Plasma Phys. 1981. V. 23. P. 89.
Звонков А.В. // Физика плазмы. 1983. Т. 9. С. 547.
Господчиков Е.Д., Чувакин П.А., Шалашов А.Г. // Физика плазмы. 2023. Т. 49. С. 953.
Иванов М.Ф., Швец В.Ф. // ДАН СССР. 1978. Т. 238. С. 1324.
Doering C.R. // Lectures in Complex Systems / Ed. L. Nadel, D.I. Stein. Boca Raton: CRC Press, 2018.
Whittaker E.T. A treatise on the analytical dynamics of particles and rigid bodies. CUP Archive, 1964.
Федоров Ф.И. // ДАН СССР. 1955. Т. 105. С. 465.
Березин И.С., Жидков Н.П. Методы вычислений. Т. II. М.: Физматлит, 1959.

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Electron temperature (a) and plasma density (b) profiles used in the simulation.

Download (89KB)

Indexing metadata

3. Fig. 2. Simulation of propagation and absorption of quasi-optical wave beams in tokamak plasma: scenario of EC heating at the third harmonic with radiation input at an angle of 5  to the vertical (upper row) and 10  to the vertical (lower row). The wave field intensity distributions | u | 2 in the poloidal cross-section of the tokamak are shown, calculated without taking into account the diffusion terms in the first order of perturbation theory (a, d) and with taking into account the diffusion terms in the second order of perturbation theory (b, d). The solid red line highlights the region of microwave power absorption, determined by the level max (dP / dV) / 2 ; the dotted line denotes the plasma boundary. The graphs (c, e) show the absorbed power profiles П(ρ) with (2) and without (1) taking into account the diffusion terms.

Download (245KB)

Indexing metadata

Username
Password
Remember me

Forgot password?	Register

Username
Password
Remember me

Forgot password?	Register

Vol 51, No 3 (2025)

Vol 51, No 3 (2025)

Wigner function method for describing electromagnetic field in plasma-like media with spatial dispersion and resonant dissipation

Full Text

Abstract

Keywords

Full Text

1. Введение

2. Функция вигнера для скалярного поля

3. Тензор вигнера для монохроматического электромагнитного поля

4. Нулевой порядок. Скалярная функция Вигнера для электромагнитного поля

5. Связь с квазиоптическим приближением

6. Первый порядок. Кинетический лучевой подход

7. Второй порядок. Диффузионное приближение

8. Начальные условия для кинетического уравнения

9. Плотность потока энергии

10. Пример

11. Заключение

Приложение 1. Генерация случайного вектора с заданной ковариацией компонент

Приложение 2. Представление дополнительного потока энергии в диспергирующей среде через тензор Вигнера

About the authors

Е. D. Gospodchikov

A. A. Balakin

A. G. Shalashov

References

Supplementary files