Физика плазмы, 2023, T. 49, № 2, стр. 154-164

Усовершенствованное квазистатическое приближение

П. В. Туев ab*, Р. И. Спицын ab, К. В. Лотов ab

a Институт ядерной физики им. Г.И. Будкера СО РАН
Новосибирск, Россия

b Новосибирский государственный университет
Новосибирск, Россия

* E-mail: p.v.tuev@inp.nsk.su

Поступила в редакцию 19.09.2022
После доработки 01.11.2022
Принята к публикации 10.11.2022

Полный текст (PDF)

Аннотация

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

Ключевые слова: плазменное кильватерное ускорение, квазистатическое приближение, моделирование, численная модель

1. ВВЕДЕНИЕ

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

Сложные явления, происходящие в плазменных кильватерных ускорителях, удается проанализировать аналитически только в простейших приближениях, поэтому основным методом теоретического исследования является численное моделирование [2]. Большинство интересных процессов происходит в небольшой области пространства, которая движется вместе с пучками. Это облегчает моделирование, так как позволяет использовать короткое движущееся окно. Однако в задаче есть сильно различающиеся временные и пространственные масштабы [3, 4], разница между которыми может достигать многих порядков, от длины волны лазерного импульса (около микрона) до полной длины ускорения (сотни метров) [57]. По этой причине не всегда возможно моделировать плазменное кильватерное ускорение с помощью численных кодов, использующих основанный на исходных уравнениях метод “частиц-в-ячейках” (PIC), и приходится использовать различные упрощенные модели [2, 8, 9].

Одной из популярных упрощенных моделей является так называемое квазистатическое приближение (QSA) [10, 11]. Оно реализовано в ряде кодов [1216], иногда в сочетании с другими упрощениями, такими как уравнение огибающей для лазерного импульса [11, 1720] или жидкостная модель для плазмы [2123]. Квазистатическое приближение опирается на тот факт, что свойства исследуемых объектов в сопутствующей системе координат изменяются гораздо медленнее, чем в лабораторной (рис. 1). Это касается и свойств плазменной волны (метка 1 на рис. 1), и плотности и скорости пучков частиц (2), и интенсивности лазерного излучения (3). Если пучки распространяются в направлении z, то замена переменных с продольной координаты z и времени t на

(1)
$s = z,\quad \xi = z - ct,$
где c – скорость света, приводит к тому, что все функции при изменении переменной s (при постоянной ξ) меняются намного медленнее, чем при изменении переменной ξ (при постоянной s), равно как и при изменении z и ct по отдельности. В то время как для обычных PIC-кодов необходимы малые шаги сетки по времени и продольной координате для разрешения быстрых изменений моделируемых величин (метка 5, рис. 1), в квазистатическом приближении требуется использовать маленький шаг только по ξ, а по s допустимы большие шаги (метка 6). Использование большого шага сетки делает QSA-коды на несколько порядков быстрее обычных PIC-кодов. Выигрыш в скорости определяется отношением длины эволюции пучка, соответствующей изменению его характерной формы, к длине волны лазерного излучения или плазменной волны, в зависимости от типа драйвера. Отметим также, что в квазистатическом приближении зависимость моделируемых величин от ξ при некотором z можно интерпретировать как их зависимость от z при фиксированном t (метка 9, рис. 1), если свойства пучка и плазмы изменяются незначительно на соответствующих пространственных и временных интервалах.

Рис. 1.

Иллюстративное изображение объектов, моделируемых в задачах плазменного кильватерного ускорения: плазменная волна (1), ускоряемый пучок (2), лазерный драйвер (3), градиент плотности плазмы (4). Расчетные сетки, используемые в обычных PIC-кодах (5), (7) и в квазистатических кодах (6), (8) для моделирования динамики пучка (5), (6) и временной эволюции плазменной волны (7), (8). Пары точек с почти одинаковым откликом плазмы (9).

Квазистатическое приближение является точным, если пучки неизменной формы распространяются в продольно однородной плазме. В этом случае в задаче появляется дополнительная симметрия: идентичные частицы плазмы, изначально расположенные в одном и том же поперечном положении ${{{\mathbf{r}}}_{ \bot }}$, но с разными z, копируют движение друг друга с некоторой временной задержкой. Симметрия приводит к еще одному преимуществу QSA: при вычислении отклика плазмы размерность задачи понижается на единицу, так как одна пространственная координата $\left( z \right)$ исчезает, сливаясь со временем t во времениподобную координату ξ. Понижение размерности позволяет использовать гораздо меньшее количество “макрочастиц” для расчета отклика плазмы. Эти “макрочастицы” представляют собой не группы реальных частиц, как в PIC-кодах, а “струи частиц”, состоящие из реальных частиц, начинающих свое движение из заданного поперечного положения с заданным начальным импульсом, но с разными z.

Третьим преимуществом QSA является его эффективность при моделировании долговременной динамики плазменной волны. В PIC-кодах для расчета состояния плазмы в следующий момент времени необходима информация о текущем состоянии соседних слоев плазмы. Поэтому моделировать временную эволюцию можно только в протяженных областях, охватывающих несколько периодов кильватерной волны ${{\lambda }_{p}}$, чтобы минимизировать влияние продольных границ (метка 7, рис. 1). В рамках квазистатического приближения для вычисления отклика плазмы по ξ (8) не требуется знать состояние соседних слоев. Модель предполагает, что эти слои копируют рассматриваемый плазменный слой с некоторой временной задержкой или опережением. В результате выигрыш в скорости вычислений в несколько раз превышает отношение ${{\lambda }_{p}}{\text{/}}\Delta z$. Благодаря этому преимуществу QSA-коды удерживают рекорд по моделированию долговременной эволюции плазменной волны [24, 25].

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

В данной работе представлено усовершенствованное квазистатическое приближение (AQSA), которое свободно от ограничений обычного квазистатического приближения и сохраняет его основные преимущества – скорость расчетов и пониженную размерность. В новом подходе учитывается обмен информацией между соседними слоями плазмы, и поэтому корректно моделируются продольные неоднородности плазмы, движение появляющихся в плазме быстрых частиц и волны с ненулевой групповой скоростью, если эти явления разрешимы на выбранной сетке моделирования. В разделе 2 сформулирована физическая модель и обсуждается ее численная реализация на примере двумерной геометрии. В разделе 3 результаты моделирования сравниваются с аналитическими решениями, а также с результатами моделирования другими QSA- и PIC-кодами. Здесь же продемонстрировано, что новая модель воспроизводит эффекты, отсутствующие в QSA-кодах. Раздел 4 посвящен обсуждению возможных направлений развития и области применимости новой модели.

2. РАСШИРЕНИЕ КВАЗИСТАТИЧЕСКОЙ МОДЕЛИ

2.1. Физическая модель

Уравнения Максвелла, записанные в виде

(2)
$\nabla \times {\mathbf{E}} = - \frac{1}{c}\frac{{\partial {\mathbf{B}}}}{{\partial t}},\quad \nabla \cdot {\mathbf{E}} = 4\pi \rho ,$
(3)
$\nabla \times {\mathbf{B}} = \frac{{4\pi }}{c}{\mathbf{j}} + \frac{1}{c}\frac{{\partial {\mathbf{E}}}}{{\partial t}},\quad \nabla \cdot {\mathbf{B}} = 0,$
где E и B – электрическое и магнитное поля, а j и ρ – плотности тока и заряда, можно объединить:

(4)
$\Delta {\mathbf{E}} - \frac{1}{{{{c}^{2}}}}\frac{{{{\partial }^{2}}{\mathbf{E}}}}{{{{\partial }^{2}}t}} = 4\pi \nabla \rho + \frac{{4\pi }}{{{{c}^{2}}}}\frac{{\partial {\mathbf{j}}}}{{\partial t}}.$

При переходе к переменным $(s,\xi )$ производные изменяются следующим образом:

(5)
$\frac{\partial }{{\partial z}} = \frac{\partial }{{\partial s}} + \frac{\partial }{{\partial \xi }},\quad \frac{\partial }{{\partial t}} = - c\frac{\partial }{{\partial \xi }}.$

Если все величины медленно зависят от $s$, то соответствующая производная является малым параметром: ${{\partial }_{s}} \ll {{\partial }_{\xi }}$. В традиционном квазистатическом приближении всеми малыми членами, содержащими ${{\partial }_{s}}$, пренебрегают. Мы же удержим члены первого порядка малости:

(6)
$\left( {{{\Delta }_{ \bot }} + 2\frac{{{{\partial }^{2}}}}{{\partial s\partial \xi }}} \right){\mathbf{E}} = 4\pi \hat {\nabla }\rho - \frac{{4\pi }}{c}\frac{{\partial {\mathbf{j}}}}{{\partial \xi }},$
где $\hat {\nabla } = \left( {{{\nabla }_{ \bot }},\;{{\partial }_{s}} + {{\partial }_{\xi }}} \right)$, а нижний индекс $ \bot $ обозначает двумерные (поперечные) векторы и операторы. Вторая производная $\partial _{{ss}}^{2}$, которой пренебрегли, отвечает за распространение излучения назад [26] и важна только для низкочастотных возмущений, которые распространяются намного медленнее скорости света [17, 27]. Аналогично, для магнитного поля получаем

(7)
$\left( {{{\Delta }_{ \bot }} + 2\frac{{{{\partial }^{2}}}}{{\partial s\partial \xi }}} \right){\mathbf{B}} = - \frac{{4\pi }}{c}\left[ {\hat {\nabla } \times {\mathbf{j}}} \right].$

Смешанные производные в уравнениях (6) и (7) отвечают за распространение свободного излучения.

Кроме уравнений электромагнитного поля необходимо также модифицировать уравнения движения макрочастиц плазмы, т.е. законы изменения их параметров: поперечного положения ${{{\mathbf{r}}}_{ \bot }}$, трех компонент импульса p и заряда q. Эти величины должны быть вычислены как функции от ξ при фиксированном значении s (рис. 2). Продольная координата макрочастицы не является параметром, так как она объединяется со временем t и переходит из параметров в аргументы. Масса m пропорциональна заряду q и не требует отдельного рассмотрения. В лабораторной системе отсчета параметры частиц плазмы изменяются в соответствии с уравнениями

(8)
$\frac{{d{{{\mathbf{r}}}_{ \bot }}}}{{dt}} = {{{\mathbf{v}}}_{ \bot }},\quad \frac{{dq}}{{dt}} = 0,$
(9)
$\frac{{d{\mathbf{p}}}}{{dt}} = q\left( {{\mathbf{E}} + \frac{1}{c}\left[ {{\mathbf{v}} \times {\mathbf{B}}} \right]} \right),$
где ${\mathbf{v}} = {\mathbf{p}}c{\text{/}}\sqrt {{{p}^{2}} + {{m}^{2}}{{c}^{2}}} $ – скорость частицы. Если обозначить произвольный параметр через $\chi (s,\xi )$, а правые части соответствующих уравнений (8), (9) через F, то в квазистатических переменных (1) получится
(10)
$\begin{gathered} \frac{{d\chi }}{{dt}} = {{\left. {\frac{{\partial \chi }}{{\partial s}}} \right|}_{\xi }}\frac{{ds}}{{dt}} + {{\left. {\frac{{\partial \chi }}{{\partial \xi }}} \right|}_{s}}\frac{{d\xi }}{{dt}} = \\ = {{\left. {\frac{{\partial \chi }}{{\partial s}}} \right|}_{\xi }}{{{v}}_{z}} + {{\left. {\frac{{\partial \chi }}{{\partial \xi }}} \right|}_{s}}({{{v}}_{z}} - c) = F, \\ \end{gathered} $
откуда

(11)
${{\left. {\frac{{\partial \chi }}{{\partial \xi }}} \right|}_{s}} = \frac{1}{{{{{v}}_{z}} - c}}\left( {F - {{{v}}_{z}}{{{\left. {\frac{{\partial \chi }}{{\partial s}}} \right|}}_{\xi }}} \right).$

Уравнение (11) описывает состояние той частицы плазмы в “струе частиц”, которая в данный момент имеет продольную координату s. Эта частица может первоначально (до прихода драйвера) находиться на некотором расстоянии от этого сечения. Таким образом, уравнение (11) описывает параметры различных физических частиц при различных ξ. Если пренебречь производной ${{\partial }_{s}}$, то уравнение (11) сводится к обычным уравнениям квазистатического приближения [28].

Для замыкания системы уравнения (6), (7) и (11) надо дополнить уравнениями для заряда и тока. Они такие же, как в квазистатическом приближении:

(12)
${\mathbf{j}} = A\sum\limits_i \frac{{{{q}_{i}}{{{\mathbf{v}}}_{i}}}}{{c - {{{v}}_{{z,i}}}}} + {{{\mathbf{j}}}_{b}},$
(13)
$\rho = A\sum\limits_i \frac{{{{q}_{i}}}}{{c - {{{v}}_{{z,i}}}}} + {{\rho }_{b}},$
где ${{{\mathbf{j}}}_{b}}$ и ${{\rho }_{b}}$ – плотности тока и заряда драйвера, A – нормировочный коэффициент, а суммирование ведется по “струям частиц”, пересекающим ячейку, в которой вычисляются j и ρ. Знаменатели в уравнениях (12) и (13) появляются потому, что вклад “струи частиц” в плотность и ток зависит от скорости макрочастиц в сопутствующем окне моделирования.

Уравнения для частиц пучка – такие же, как в обычной квазистатической модели (см., например, [28]):

(14)
$\frac{{d{{{\mathbf{r}}}_{{b \bot }}}}}{{ds}} = \frac{{{{{\mathbf{v}}}_{{b \bot }}}}}{c},\quad \frac{{d{{\xi }_{b}}}}{{ds}} = \frac{{{{{v}}_{{bz}}}}}{c} - 1,$
(15)
$\frac{{d{{{\mathbf{p}}}_{b}}}}{{ds}} = \frac{{{{q}_{b}}}}{c}{\mathbf{E}} + \frac{{{{q}_{b}}}}{{{{c}^{2}}}}\left[ {{{{\mathbf{v}}}_{b}} \times {\mathbf{B}}} \right],$
(16)
${{{\mathbf{v}}}_{b}} = \frac{{{{{\mathbf{p}}}_{b}}c}}{{\sqrt {m_{b}^{2}{{c}^{2}} + p_{b}^{2}} }},$
где $({{{\mathbf{r}}}_{{b \bot }}},{{\xi }_{b}})$ – координаты макрочастицы пучка в окне моделирования, а ${{m}_{b}}$, ${{q}_{b}}$ и ${{{\mathbf{p}}}_{b}}$ – ее масса, заряд и импульс. Модель можно дополнить уравнением огибающей для лазерного импульса и пондеромоторной силой, действующей на частицы плазмы [11], но в данной работе это не пригодится.

2.2. Численная реализация

Новая модель реализована на базе существующего квазистатического 2D3V кода LCODE [29]. Представленная численная схема не оптимизирована для достижения наилучшей производительности и имеет целью показать, что новая модель в принципе может быть включена в квазистатический код. Код может работать в декартовой или осесимметричной системе координат, и модель проверялась в обеих. Далее выписаны уравнения только для декартовой системы координат, их модификация для цилиндрической геометрии аналогична работе [30]. В рассмотренных примерах три компоненты поля тождественно равны нулю (${{E}_{\varphi }} = {{B}_{r}} = {{B}_{z}} = 0$ в цилиндрической и ${{E}_{y}} = {{B}_{x}} = {{B}_{z}} = 0$ в декартовой системах координат), поэтому уравнения для этих компонент опускаем.

Отклик плазмы в LCODE рассчитывается послойно в направлении убывания координаты ξ (справа налево на рис. 2). Поскольку для расчета поля ${{E}_{x}}$ в уравнении (6) требуются производные тока по ξ, используется следующая схема предиктор-корректор:

• рассчитываем движение макрочастиц плазмы из слоя a в слой b под действием полей слоя a;

• вычисляем ток и плотность заряда в слое b;

• вычисляем все поля в слое b;

• рассчитываем движение частиц плазмы из слоя a в слой b под действием средних полей слоев a и b;

• снова вычисляем ток и плотность заряда в слое b;

• снова вычисляем все поля в слое b;

• в третий раз рассчитываем движение частицы плазмы из слоя a в слой b под действием средних полей.

Чтобы обеспечить устойчивость алгоритма, методом конечных разностей вместо уравнения (6) решается следующее уравнение для поперечного электрического поля:

(17)
$\left( {\frac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}} + 2\frac{{{{\partial }^{2}}}}{{\partial s\partial \xi }}} \right){{E}_{x}} - {{E}_{x}} = 4\pi \frac{{\partial \rho }}{{\partial x}} - \frac{{4\pi }}{c}\frac{{\partial {{j}_{x}}}}{{\partial \xi }} - {{\tilde {E}}_{x}},$
где ${{\tilde {E}}_{x}}$ – некоторое предсказание для ${{E}_{x}}$. В конечно-разностном представлении это уравнение принимает вид
$\begin{gathered} \frac{1}{2}\left( {\frac{{E_{{x,n,l + 1}}^{m} - 2E_{{x,n,l}}^{m} + E_{{x,n,l - 1}}^{m}}}{{\Delta {{x}^{2}}}} + } \right. \\ \, + \left. {\frac{{E_{{x,n,l + 1}}^{{m - 1}} - 2E_{{x,n,l}}^{{m - 1}} + E_{{x,n,l - 1}}^{{m - 1}}}}{{\Delta {{x}^{2}}}}} \right) + \\ \end{gathered} $
(18)
$\, + \frac{2}{{\Delta s}}\left( {\frac{{ - 3E_{{x,n,l}}^{m} + 4E_{{x,n - 1,l}}^{m} - E_{{x,n - 2,l}}^{m}}}{{2\Delta \xi }} - } \right.$
$\begin{gathered} \left. {\, - \frac{{ - 3E_{{x,n,l}}^{{m - 1}} + 4E_{{x,n - 1,l}}^{{m - 1}} - E_{{x,n - 2,l}}^{{m - 1}}}}{{2\Delta \xi }}} \right) - \frac{{E_{{x,n,l}}^{m} + E_{{x,n,l}}^{{m - 1}}}}{2} = \\ \, = 4\pi \left( {\frac{{\rho _{{n,l + 1}}^{{m - \frac{1}{2}}} - \rho _{{n,l - 1}}^{{m - \frac{1}{2}}}}}{{2\Delta x}} + \frac{1}{c}\frac{{j_{{x,n,l}}^{{m - \frac{1}{2}}} - j_{{x,n - 1,l}}^{{m - \frac{1}{2}}}}}{{\Delta \xi }}} \right) - \tilde {E}_{{x,n,l}}^{{m - \frac{1}{2}}}, \\ \end{gathered} $
где индексы l, m, n обозначают слои сетки по x, s, ξ, соответственно (рис. 3), $\Delta x$, $\Delta s$, $\Delta \xi $ – шаги сетки, а полуцелые индексы обозначают полусумму значений в близлежащих точках сетки, например,

(19)
$\tilde {E}_{{x,n,l}}^{{m - \frac{1}{2}}} = \frac{{\tilde {E}_{{x,n,l}}^{m} + \tilde {E}_{{x,n,l}}^{{m - 1}}}}{2}.$
Рис. 2.

Расчет отклика плазмы в квазистатическом приближении.

При вычислении поля в слое с индексом n в первый раз (предиктор) $\tilde {E}_{{x,n,l}}^{m} = E_{{x,n - 1,l}}^{m}$, иначе (корректор) $\tilde {E}_{{x,n,l}}^{m} = E_{{x,n - \frac{1}{2},l}}^{m}$, где поле $E_{{x,n,l}}^{m}$ найдено на шаге предиктора. Выбор конечно-разностной аппроксимации производных поля по ξ связан с принципом причинности. Возмущение в окне моделирования, движущемся со скоростью света, не может распространяться слева направо. Все возмущения в системе вызваны пучками и распространяются вместе с ними или отстают от них. Это означает, что последующие слои по ξ не могут влиять на предыдущие, поэтому следует использовать только правые разностные схемы для производных по ξ. Благодаря несимметричным разностным схемам возможно конвейерное распараллеливание расчетов [13]. Противопотоковая (upwind) схема второго порядка для производной поля по ξ обеспечивает устойчивость численной схемы при вычислении смешанных производных [19, 26]. Для поперечного лапласиана используется устойчивая схема Кранка–Николсона [31], которая хорошо зарекомендовала себя при решении уравнения огибающей лазерного импульса в моделировании плазменного кильватерного ускорения. В результате получается неявная численная схема, которую при необходимости можно обобщить на трехмерный случай с помощью неявного метода переменных направлений [32, 33].

Чтобы избежать вычисления производной $\partial {{j}_{z}}{\text{/}}\partial \xi $ в продольной компоненте уравнения (6), можно воспользоваться уравнением непрерывности

(20)
$\frac{{\partial \rho }}{{\partial t}} + \nabla \cdot {\mathbf{j}} = 0$
и получить
(21)
$\left( {\frac{{{{\partial }^{2}}}}{{\partial {{x}^{2}}}} + 2\frac{{{{\partial }^{2}}}}{{\partial s\partial \xi }}} \right){{E}_{z}} = \frac{{4\pi }}{c}\frac{{\partial {{j}_{x}}}}{{\partial x}} + 4\pi \frac{\partial }{{\partial s}}\left( {\rho + \frac{{{{j}_{z}}}}{c}} \right),$
что в конечных разностях дает

$\begin{gathered} \frac{{E_{{z,n,l + 1}}^{m} - 2E_{{z,n,l}}^{m} + E_{{z,n,l - 1}}^{m}}}{{\Delta {{x}^{2}}}} + \\ \, + \frac{2}{{\Delta s}}\left( {\frac{{ - 3E_{{z,n,l}}^{m} + 4E_{{z,n - 1,l}}^{m} - E_{{z,n - 2,l}}^{m}}}{{2\Delta \xi }} - } \right. \\ \end{gathered} $
(22)
$\, - \left. {\frac{{ - 3E_{{z,n,l}}^{{m - 1}} + 4E_{{z,n - 1,l}}^{{m - 1}} - E_{{z,n - 2,l}}^{{m - 1}}}}{{2\Delta \xi }}} \right) = $
$\begin{gathered} \, = 4\pi \left( {\frac{1}{c}\frac{{j_{{x,n,l + 1}}^{m} - j_{{x,n,l - 1}}^{m}}}{{2\Delta x}} + } \right. \\ \, + \left. {\frac{{\rho _{{n,l}}^{m} - \rho _{{n,l}}^{{m - 1}}}}{{\Delta s}} + \frac{1}{c}\frac{{j_{{z,n,l}}^{m} - j_{{z,n,l}}^{{m - 1}}}}{{\Delta s}}} \right). \\ \end{gathered} $

Системы уравнений (18) и (22) решаются методом прогонки.

Уравнение для магнитного поля в двумерной геометрии можно упростить. Разность z-компонент первого из уравнений (3) и второго из уравнений (2) дает

(23)
$\frac{\partial }{{\partial x}}({{B}_{y}} - {{E}_{x}}) = \frac{{4\pi }}{c}{{j}_{z}} - 4\pi \rho + \frac{{\partial {{E}_{z}}}}{{\partial s}},$
откуда

(24)
$\begin{gathered} \frac{{B_{{y,n,l}}^{m} - B_{{y,n,l - 1}}^{m}}}{{\Delta x}} - \frac{{E_{{x,n,l}}^{m} - E_{{x,n,l - 1}}^{m}}}{{\Delta x}} = \frac{{4\pi }}{c}j_{{z,n,l - \frac{1}{2}}}^{m} - \\ \, - 4\pi \rho _{{z,n,l - \frac{1}{2}}}^{m} + \frac{{E_{{z,n,l - \frac{1}{2}}}^{m} - E_{{z,n,l - \frac{1}{2}}}^{{m - 1}}}}{{\Delta s}}. \\ \end{gathered} $

Это уравнение позволяет вычислить магнитное поле с точностью до константы. В цилиндрической геометрии ${{B}_{\varphi }} = {{E}_{r}} = 0$ на оси, что и определяет константу. В декартовой геометрии для этого используется $x$-компонента первого из уравнений (3)

(25)
$\frac{\partial }{{\partial \xi }}({{B}_{y}} - {{E}_{x}}) = - \frac{{4\pi }}{c}{{j}_{x}} - \frac{{\partial {{B}_{y}}}}{{\partial s}}.$

Достаточно найти константу только для одного значения поперечного индекса l. Это значение $\left( {l{\kern 1pt} '} \right)$ можно выбрать на периферии окна моделирования, где токи и поля малы, и производной ${{\partial }_{s}}$ можно пренебречь. В конечно-разностной аппроксимации

(26)
$\frac{{B_{{y,n,l{\kern 1pt} '}}^{m} - B_{{y,n - 1,l{\kern 1pt} '}}^{m}}}{{\Delta \xi }} = \frac{{E_{{x,n,l{\kern 1pt} '}}^{m} - E_{{x,n - 1,l{\kern 1pt} '}}^{m}}}{{\Delta \xi }} + \frac{{4\pi }}{c}j_{{x,n - \frac{1}{2},l{\kern 1pt} '}}^{m},$
откуда находится $B_{{y,n,l{\kern 1pt} '}}^{m}$, если известны остальные слагаемые.

Для расчета параметров частиц используется неявная схема, которая для уравнения (11) имеет вид

(27)
$\frac{{\chi _{n}^{m} - \chi _{{n - 1}}^{m}}}{{\Delta \xi }} = \frac{1}{{c - {v}_{{z,n - \frac{1}{2}}}^{m}}}\left( {F_{{n - \frac{1}{2}}}^{m} - {v}_{{z,n - \frac{1}{2}}}^{{m - \frac{1}{2}}}\frac{{\chi _{n}^{m} - \chi _{n}^{{m - 1}}}}{{\Delta s}}} \right).$

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

Алгоритм вычисления эволюции пучка остается неизменным. Уравнения движения решаются с помощью модифицированного метода Эйлера (метод средней точки). Поля, действующие на макрочастицу, линейно интерполируются к предсказанному местоположению макрочастицы на половине временного шага.

3. ПРОВЕРКА НОВОЙ МОДЕЛИ

3.1. Распространение лазерного импульса в плазме

Моделирование распространения короткого двумерного лазерного импульса через плазму низкой плотности ${{n}_{0}}$ в линейном режиме покажет, как возмущения с ненулевой групповой скоростью описываются новой моделью.

Групповая скорость ${{{v}}_{g}}$ лазерного импульса отличается от скорости света c не только из-за воздействия плазмы, но и из-за конечного поперечного размера импульса. Адаптируя теорию [34] к плоскому двумерному случаю, можно получить значение групповой скорости на оси вблизи точки фокуса:

(28)
$\frac{{c - {{{v}}_{g}}}}{c} = \frac{1}{{x_{0}^{2}k_{0}^{2}}} + \frac{{k_{p}^{2}}}{{2k_{0}^{2}}},$
где ${{x}_{0}}$ – ширина импульса в перетяжке, ${{k}_{0}}$ – волновое число импульса, ${{k}_{p}} = \sqrt {4\pi {{n}_{0}}{{e}^{2}}{\text{/}}\left( {{{m}_{e}}{{c}^{2}}} \right)} $ – волновое число плазмы, e – элементарный заряд, ${{m}_{e}}$ – масса электрона. Поперечный размер лазерного импульса ${{x}_{s}}$ растет согласно параксиальной теории [34]:
(29)
${{x}_{s}} = {{x}_{0}}\sqrt {1 + \frac{{{{s}^{2}}}}{{s_{R}^{2}}}} ,$
где ${{s}_{R}} = {{k}_{0}}x_{0}^{2}{\text{/}}2$ – длина Рэлея, а расстояние s измеряется от точки перетяжки.

В предположении, что лазерный импульс линейно поляризован, его поперечное электрическое поле при $s = 0$ задаем в следующем виде:

(30)
${{E}_{x}} = {{E}_{m}}\sin \left( {{{k}_{0}}\xi } \right)\exp \left( { - \frac{{{{x}^{2}}}}{{x_{0}^{2}}} - \frac{{{{{(\xi - {{\xi }_{c}})}}^{2}}}}{{{{L}^{2}}}}} \right),$
где ${{\lambda }_{0}} = 2\pi {\text{/}}{{k}_{0}} = 810$ нм, ${{x}_{0}} = 13$ мкм, $L = 9$ мкм, ${{E}_{m}} = 1.26$ МВ/м (рис. 4). Поперечное магнитное поле ${{B}_{y}} = {{E}_{x}}$. Продольное электрическое поле ${{E}_{z}}$ удовлетворяет уравнению $\nabla \cdot {\mathbf{E}} = 0$. Максимум интенсивности поля изначально расположен в точке ${{\xi }_{c}} = - 30$ мкм, которая является центром окна моделирования $(x,\xi )$ размером 180 мкм × × 60 мкм. Шаги сетки $\Delta \xi = {{\lambda }_{0}}{\text{/}}40$, $\Delta x = {{\lambda }_{0}}{\text{/}}4$, $\Delta s = 25{{\lambda }_{0}}$.

Рис. 3.

Расчетная сетка.

Рис. 4.

Моделируемый лазерный импульс при $s = 0$ (a) и $s = 2{{s}_{R}}$ (б). Для наглядности показана только часть окна моделирования с различными вертикальными и горизонтальными масштабами.

В отличие от других QSA-кодов [11, 18, 19], распространение этого импульса описывается плазменным солвером. Результаты моделирования количественно воспроизводят теоретически предсказанную скорость волнового пакета (рис. 5) и его поперечное расплывание (рис. 6). Чтобы найти групповую скорость из результатов моделирования, огибающая импульса на оси аппроксимировалась гауссовской функцией, и измерялась средняя скорость максимума за время $0.13{{s}_{R}}{\text{/}}c$, начиная от точки фокуса. Ширина импульса измерялась в соответствии с формулой

(31)
${{x}_{s}} = 2\sqrt {\frac{{\left\langle {{{x}^{2}}E_{x}^{2}} \right\rangle }}{{\left\langle {E_{x}^{2}} \right\rangle }}} ,$
где угловые скобки означают усреднение по окну моделирования. Важно отметить, что шаг продольной сетки $\Delta s \gg {{\lambda }_{0}}$, так что AQSA сохраняет главное преимущество QSA – быстродействие.

Рис. 5.

Групповая скорость лазерного импульса в плазме различной плотности.

3.2. Плазма с продольным градиентом плотности

Этот тест показывает, как новая модель описывает сильнонелинейную плазменную волну в присутствии градиентов плотности. Короткий осесимметричный электронный сгусток возбуждает волну в режиме каверны (рис. 7). Электроны плазмы в задней части каверны приобретают большую продольную скорость (близкую к c), улетают далеко вперед от своего первоначального местоположения и могут передавать информацию из одного слоя плазмы в другой. В продольно однородной плазме этот эффект не влияет на точность QSA-модели благодаря одинаковым свойствам плазмы в разных сечениях, но при наличии градиента плотности он может быть значительным.

Рис. 6.

Расширение импульса в вакууме.

Электронный пучок имеет следующие параметры: среднеквадратичные поперечный и продольный размеры 3 мкм и 10 мкм соответственно, релятивистский фактор 4 × 104, пиковый ток 11 кА. Такой пучок создает в плазме каверну, которая хорошо подходит для демонстрации возможностей новой модели и в то же время не имеет особенностей со слишком малым пространственным масштабом, которые трудно разрешить на сетке. Продольный профиль плотности плазмы содержит область линейного роста от 1.75 × × 1017 см–3 при $s = 50$ мкм до 2.5 × 1017 см–3 при $s = 220$ мкм, что соответствует изменению толщины плазменного скин-слоя $k_{p}^{{ - 1}}$ от 12.7 мкм до 10.6 мкм. До и после этой области плотность плазмы постоянна (рис. 8, верхний ряд). Моделирование проводится в цилиндрических координатах $(r,\xi )$ или $(r,z)$.

Рис. 7.

Суммарная электронная плотность n пучка и плазмы в случае, когда электронный пучок распространяется в однородной плазме (при s = 20 мкм) до области градиента плотности. Зеленым прямоугольником отмечена область, показанная на рис. 8. Моделирование выполнено с помощью FBPIC (верхняя половина) и модели AQSA на основе LCODE (нижняя половина).

Мы сравниваем расчеты, выполненные с помощью QSA-кода LCODE, осесимметричного PIC-кода FBPIC [35], и AQSA модели, реализованной на основе LCODE. Размер окна моделирования составляет 250 мкм по r и 170 мкм по ξ. В моделировании QSA и AQSA шаги сетки составляют $\Delta r = \Delta \xi = 0.1{\kern 1pt} $ мкм, $\Delta s = 10{\kern 1pt} $ мкм. В моделировании FBPIC $\Delta r = \Delta z = \Delta (ct) = 0.1{\kern 1pt} $ мкм. Во всех расчетах использовалось 12 макрочастиц плазмы на ячейку, а пучок описывался 2 × × 105 макрочастицами. Для сравнения квазистатических и PIC-расчетов их результаты нужно представить в одних и тех же координатах, и здесь это – квазистатические координаты $(r,\xi )$. Зависимости от ξ имеют смысл временных зависимостей при заданной координате z. В однородной плазме эти зависимости также можно интерпретировать как портреты пучка в некоторый момент времени, что неверно для случая плазмы с продольным градиентом плотности. Стандартная диагностика в FBPIC не позволяет выводить временные зависимости при фиксированном z, поэтому состояние плазмы и поля сохраняются каждый пятый временной шаг (для упрощения работы с большими объемами данных), если интересующая нас координата ${{s}_{0}}$ находится в окне моделирования. Затем для каждого сохраненного состояния выбираются значения в пяти ближайших к ${{s}_{0}}$ сечениях и рассматриваются как пять последовательных точек искомой временной зависимости при $z = {{s}_{0}}$.

Все три подхода согласуются в областях однородной плотности плазмы (рис. 7 и 8, левая и правая колонки). В области градиента плотности классический QSA-код дает иной отклик плазмы, чем AQSA и PIC (средний столбец на рис. 8). В отличие от QSA, AQSA воспроизводит удлинение каверны (метка 1 на рис. 8) и появление протяженной области высокой плотности вблизи оси (метка 2), наблюдаемые при PIC моделировании. Оба эффекта вызваны продольно движущимися электронами плазмы, которые приходят из более ранних областей, где плотность плазмы ниже, а длина волны больше. Причина различной формы пика поля в AQSA и FBPIC (метка 3, рис. 8) еще не до конца ясна. Амплитуда пика зависит от разброса электронов по поперечной скорости [28], поэтому это различие может появиться из-за разного численного нагрева в двух кодах. Форма “хвостовой” волны [36], наблюдаемой после области градиента плотности в QSA (метка 4, рис. 8), также отличается от AQSA или FBPIC. Эта волна (видимая как гребень плотности) образуется высокоэнергичными электронами плазмы, ускоренными в области градиента. Поскольку плотность плазмы там ниже, эти электроны появляются при больших ${\text{|}}\xi {\text{|}}$, и фронт “хвостовой” волны изгибается, что видно в моделировании AQSA и FBPIC, но не в QSA.

Рис. 8.

Сравнение расчетов с использованием моделей QSA, FBPIC и AQSA: плазменная каверна в однородной плазме (левая и правая колонки) и в плазме с градиентом плотности (средняя колонка). В первом (верхнем) ряду показано расположение рассматриваемых сечений (оранжевые точки) на профиле плотности плазмы. Во втором ряду показано продольное электрическое поле ${{E}_{z}}(\xi )$, рассчитанное с помощью рассмотренных моделей, а в остальных рядах – карты электронной плотности, полученные с помощью QSA, FBPIC и AQSA, соответственно. Конечная плотность плазмы ${{n}_{0}} = 2.5 \times {{10}^{{17}}}$ см–3 определяет единицу электрического поля ${{E}_{0}} = \sqrt {4\pi {{n}_{0}}{{m}_{e}}{{c}^{2}}} $. Синие вертикальные линии на картах плотности показывают продольное положение минимума поля в моделировании AQSA для облегчения сравнения. Цифры в кружках указывают на обсуждаемые в тексте особенности.

В заключение следует отметить, что модель AQSA воспроизводит особенности, возникающие из-за продольной неоднородности плазмы, аналогично PIC-коду FBPIC, но в 80 раз быстрее (около 560 часов процессорного времени для F-BPIC против 7 для LCODE), что близко к теоретической оценке $\Delta s{\text{/}}\Delta \xi = 100$.

4. ОБСУЖДЕНИЕ

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

Изложенная модель является лишь первым шагом к более полному и быстрому моделированию кильватерного ускорения при помощи квазистатического приближения. Уравнения и численная схема могут быть улучшены путем учета производной $\partial _{{ss}}^{2}$. Полные уравнения для полей по форме похожи на уравнение эволюции огибающей лазерного импульса, поэтому опыт, полученный при моделировании последней [26, 27, 37, 38], поможет построить эффективную численную схему. Очевидно, следующим важным логическим шагом в развитии модели должно стать включение в нее эффектов захвата. Захват электронов плазмы кильватерной волной с последующим их ускорением широко используется для формирования витнесса [3941]. Изложенная в данной работе модель пока не включает этот процесс. Технически это происходит потому, что захват частиц сопровождается уменьшением пространственного масштаба изменения параметров частиц до уровня, который невозможно разрешить. Именно по этой причине мы протестировали новую модель на примере плазмы с положительным градиентом плотности. Отрицательные градиенты могут привести к захвату частиц. Однако с частицами, которые захватываются волной, можно обращаться особым образом, аналогично частицам пучка. Этот подход дает качественно правильные результаты в рамках традиционного QSA-подхода [42, 43], и должен работать еще лучше с AQSA.

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

В то время как с помощью квазистатического приближения можно очень эффективно моделировать долговременную эволюцию плазменной волны, нет уверенности, что с AQSA это удастся делать так же хорошо. Различие траекторий и параметров отдельных частиц плазмы при разных s растет с увеличением длины окна моделирования, определяющей длительность эволюции волны. Чтобы корректно описать это различие, шаг сетки $\Delta s$ необходимо уменьшить, что обесценивает основное преимущество квазистатического приближения.

Авторы благодарны И.А. Шалимовой и А.А. Горну за полезные обсуждения. Численные эксперименты проведены на вычислительном кластере “Академик В.М. Матросов” [44].

Авторы заявляют, что у них нет конфликта интересов.

Список литературы

  1. Albert F., Couprie M.E., Debus A., Downer M.C., Faure J., Flacco A., Gizzi L.A., Grismayer T., Huebl A., Joshi C., Labat M., Leemans W.P., Maier A.R., Mangles S.P.D., Mason P., Mathieu F., Muggli P., Nishiuchi M., Oster-hoff J., Rajeev P.P., Schramm U., Schreiber J., Tho-mas A.G.R., Vay J.-L., Vranic M., Zeil K. // New J. Phys. 2021. V. 23. P. 031101. https://doi.org/10.1088/1367-2630/abcc62

  2. Vay J.-L., Lehe R. // Rev. Accelerator Science Technology. 2016. V. 9. P. 165. https://doi.org/10.1142/S1793626816300085

  3. Lotov K.V. // Nuclear Instr. Methods A. 1998. V. 410. P. 461. https://doi.org/10.1016/S0168-9002(98)00178-8

  4. Burdakov A.V., Kudryavtsev A.M., Logatchov P.V., Lo-tov K.V., Petrenko A.V., Skrinsky A.N. // Plasma Phys. Rep. 2005. V. 31. P. 292. [Бурдаков А.В., Кудряв-цев А.М., Логачев П.В., Лотов К.В., Петренко А.В., Скринский А.Н. // Физика плазмы, 2005, Т. 31, C. 327–335.]https://doi.org/10.1134/1.1904145

  5. Schroeder C.B., Esarey E., Geddes C.G.R., Benedetti C., Leemans W.P. // Phys. Rev. ST Accel. Beams. 2010. V. 13. P. 101301. https://doi.org/10.1103/PhysRevSTAB.13.101301

  6. Nakajima K., Deng A., Zhang X., Shen B., Liu J., Li R., Xu Z., Ostermayr T., Petrovics S., Klier C., Iqbal K., Ruhl H., Tajima T. // Phys. Rev. ST Accel. Beams. 2011. V. 14. P. 091301. https://doi.org/10.1103/PhysRevSTAB.14.091301

  7. Schroeder C.B., Esarey E., Leemans W.P. // Phys. Rev. ST Accel. Beams, 2012. V. 15. P. 051301. https://doi.org/10.1103/PhysRevSTAB.15.051301

  8. Vay J.-L. // Phys. Rev. Lett. 2007. V. 98. P. 130405. https://doi.org/10.1103/PhysRevLett.98.130405

  9. Vay J.-L., Geddes C.G.R., Cormier-Michel E., Gro-te D.P. // J. Computational Phys. 2011. V. 230. P. 5908. https://doi.org/10.1016/j.jcp.2011.04.003

  10. Sprangle P., Esarey E., Ting A. // Phys. Rev. Lett. 1990. V. 64. P. 2011. https://doi.org/10.1103/PhysRevLett.64.2011

  11. Mora P., Antonsen T.M. // Phys. Plasmas. 1997. V. 4. P. 217. https://doi.org/10.1063/1.872134

  12. Jain N., Palastro J., Antonsen T.M., Mori W.B., An W. // Phys. Plasmas, 2015. V. 22. P. 023103. https://doi.org/10.1063/1.4907159

  13. Sosedkin A.P., Lotov K.V. // Nuclear Instr. Methods A. 2016. V. 829. P. 350. https://doi.org/10.1016/j.nima.2015.12.032

  14. An W., Decyk V.K., Mori W.B., Antonsen Jr. T.M. // J. Computational Phys. 2013. V. 250. P. 165. https://doi.org/10.1016/j.jcp.2013.05.020

  15. Mehrling T., Benedetti C., Schroeder C.B., Osterhoff J. // Plasma Phys. Control. Fusion, 2014. V. 56. P. 084012. https://doi.org/10.1088/0741-3335/56/8/084012

  16. Pukhov A., Farmer J.P. // Phys. Rev. Lett. 2018. V. 121. P. 264801. https://doi.org/10.1103/PhysRevLett.121.264801

  17. Zhu W., Palastro J.P., Antonsen T.M. // Phys. Plasmas, 2012. V. 19. P. 033105. https://doi.org/10.1063/1.3691837

  18. Huang C., Decyk V.K., Ren C., Zhou M., Lu W., Mo-ri W.B., Cooley J.H., Antonsen Jr.T.M., Katsouleas T. // J. Computational Phys. 2006. V. 217. P. 658. https://doi.org/10.1016/j.jcp.2006.01.039

  19. Спицын Р.И. Магистерская дисс. Новосибирский государственный университет, 2016. https://www.inp.nsk.su/~dep_plasma/dip/Spitsyn_m.pdf.

  20. Terzani D., Benedetti C., Schroeder C.B., Esarey E. // Phys. Plasmas. 2021. V. 28. P. 063105. https://doi.org/10.1063/5.0050580

  21. Sprangle P., Esarey E., Krall J., Joyce G., Phys. Rev. Lett., 1992. V. 69. P. 2200. https://doi.org/10.1103/PhysRevLett.69.2200

  22. Esarey E., Sprangle P., Krall J., Ting A., Joyce G. // Phys. Fluids B. 1993. V. 5. P. 2690. https://doi.org/10.1063/1.860707

  23. Lotov K.V. // Phys. Plasmas. 1998. V. 5. P. 785. https://doi.org/10.1063/1.872765

  24. Zgadzaj R., Silva T., Khudyakov V.K., Sosedkin A., Al-len J., Gessner S., Li Z., Litos M., Vieira J., Lotov K.V., Hogan M.J., Yakimenko V., Downer M.C. // Nature Comm. 2020. V. 11. P. 4753. https://doi.org/10.1038/s41467-020-18490-w

  25. Khudiakov V.K., Lotov K.V., Downer M.C. // Plasma Phys. Control. Fusion. 2022. V. 64. P. 045003. https://doi.org/10.1088/1361-6587/ac4523

  26. Benedetti C., Schroeder C.B., Geddes C.G.R., Esarey E., Leemans W.P. // Plasma Phys. Control. Fusion. 2018. V. 60. P. 014002. https://doi.org/10.1088/1361-6587/aa8977

  27. Zhu W., Palastro J.P., Antonsen T.M. // Phys. Plasmas. 2013. V. 20. P. 073103. https://doi.org/10.1063/1.4813245

  28. Lotov K.V. // Phys. Rev. ST Accel. Beams. 2003. V. 6. P. 061301. https://doi.org/10.1103/PhysRevSTAB.6.061301

  29. https://lcode.info/.

  30. See the LCODE manual at https://lcode.info/site-files/manual.pdf.

  31. Crank J., Nicolson P. // Mathematical Proceed. Cambridge Philosophical Soc. 1947. V. 43. P. 50. https://doi.org/10.1017/S0305004100023197

  32. Peaceman D.W., Rachford H.H. // J. Soc. Industrial Applied Math. 1955. V. 3. P. 28. https://doi.org/10.1137/0103003

  33. Douglas J. // J. Soc. Industrial Applied Math. 1955. V. 3. P. 42. https://doi.org/10.1137/0103004

  34. Esarey E., Leemans W.P. // Phys. Rev. E. 1999. V. 59. P. 1082. https://doi.org/10.1103/PhysRevE.59.1082

  35. Lehe R., Kirchen M., Andriyash I.A., Godfrey B.B., Vay J.-L. // Computer Phys. Communications. 2016. V. 203. P. 66. https://doi.org/10.1016/j.cpc.2016.02.007

  36. Luo J., Chen M., Zhang G.-B., Yuan T., Yu J.-Y., Shen Z.-C., Yu L.-L., Weng S.-M., Schroeder C. B., Esa-rey E. // Phys. Plasmas. 2016. V. 23. P. 103112. https://doi.org/10.1063/1.4966047

  37. Massimo F., Beck A., Derouillat J., Grech M., Lobet M., Perez F., Zemzemi I., Specka A. // Plasma Phys. Control. Fusion. 2019. V. 61. P. 124001. https://doi.org/10.1088/1361-6587/ab49cf

  38. Terzani D., Londrillo P. // Computer Phys. Communicat. 2019. V. 242. P. 49. https://doi.org/10.1016/j.cpc.2019.04.007

  39. Pukhov A., Meyer-ter-Vehn J. // Appl. Phys. B, 2002. V. 74. P. 355. https://doi.org/10.1007/s003400200795

  40. Malka V. // Phys. Plasmas, 2012. V. 19. P. 055501. https://doi.org/10.1063/1.3695389

  41. Esarey E., Schroeder C.B., Leemans W.P. // Rev. Mod. Phys. 2009. V. 81. P. 1229. https://doi.org/10.1103/RevModPhys.81.1229

  42. Morshed S., Antonsen T.M., Palastro J.P. // Phys. Plasmas, 2010. V. 17. P. 063106. https://doi.org/10.1063/1.3432685

  43. Tuev P.V., Lotov K.V. Proc. 47th EPS Conference on Plasma Phys. 2021. P. 2.2004. http://ocs.ciemat.es/EPS2021PAP/pdf/P2.2004.pdf.

  44. Irkutsk Supercomputer Center of SB RAS (available at: http://ocs.ciemat.es/EPS2021PAP/pdf/P2.2004.pdf).

Дополнительные материалы отсутствуют.