ЖЭТФ, 2022, том 161, вып. 4, стр. 570-582

УСКОРЕНИЕ КОСМИЧЕСКИХ ЛУЧЕЙ ДО ЭНЕРГИЙ ВЫШЕ

10¹⁵ эВ ТРАНСРЕЛЯТИВИСТСКИМИ УДАРНЫМИ ВОЛНАМИ

А. М. Быков^*, С. М. Осипов^**, В. И. Романский

Физическо-технический институт им. А. Ф. Иоффе Российской академии наук

194021, Санкт-Петербург, Россия

Поступила в редакцию 8 сентября 2021 г.,

после переработки 9 ноября 2021 г.

Принята к публикации 30 ноября 2021 г.

МГД-течения с четырех-скоростями Γβ ∼ 1 (трансрелятивистские), реализующиеся в сверхновых звездах

типа Ib/c и после слияния нейтронных звезд, представляют интерес как потенциальные источники косми-

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

в астрофизических объектах с трансрелятивистскими МГД-течениями и ударными волнами. Построе-

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

ускоренных частиц на локальную структуру течения в окрестности ударной волны и механизмов сверх-

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

Показано, что максимальные энергии протонов, ускоренных трансрелятивистскими МГД-течениями, мо-

гут заметно превышать 10¹⁵ эВ и вносить вклад в наблюдаемые потоки космических лучей в окрестности

излома спектра.

Статья для специального выпуска ЖЭТФ, посвященного 100-летию А. Е. Чудакова

DOI: 10.31857/S0044451022040113

О природе источников КЛ с энергиями в диа-

EDN: DQJPRP

пазоне 1-1000 ПэВ, вносящих вклад в наблюдае-

мые потоки КЛ [13], нет общего мнения, хотя име-

1. ВВЕДЕНИЕ

ется ряд моделей [14,15]. В качестве одного из воз-

можных источников рассматривают подкласс сверх-

Совокупность наблюдений указывает на то, что

новых типа Ib/c, характеризующихся умеренно ре-

источники космических лучей (КЛ) в Галактике с

лятивистскими скоростями u_f истечения, с удар-

энергиями до излома спектра около петаэлектрон-

ной волной, распространяющейся в ветер массивной

вольта (1 ПэВ = 10¹⁵ эВ), вероятно, связаны с остат-

звезды прародителя сверхновой [16-18]. Для квази-

ками сверхновых звезд [1, 2], а источниками КЛ

стационарных МГД-течений анализ (см., в частнос-

сверхвысоких энергий выше 1000 ПэВ могут быть

ти, [3]), основанный на условии, что время ускоре-

релятивистские истечения активных ядер галактик

ния частицы должно быть меньше динамического

[3-5]. Поиск гамма-излучения источников КЛ черен-

времени течения, дает приближенную, но удобную

ковскими телескопами был предложен более 50 лет

оценку максимальной энергии E_max частицы заряда

назад [6] и в последние годы позволил детектиро-

Z, ускоренной МГД-течением с плотностью потока

вать гамма-фотоны с энергиями до 100 ТэВ от попу-

L_M магнитной энергии на единицу телесного угла:

ляции молодых галактических остатков сверхновых

√

[7-9]. Модели ускорения космических лучей удар-

φ(β_f )

L_M

E_max [эВ] ≈ Z · 10¹⁸

(1)

ными волнами в наблюдаемых остатках сверхновых

Γ_f

3 · 10⁴¹ эрг/с

различных типов [10-12] демонстрируют возмож-

Данная оценка не учитывает потерь энергии части-

ность ускорения КЛ до энергий, соответствующих

цы на излучение, что в большинстве случаев при-

излому наблюдаемого спектра КЛ.

менимо для ядер. Она предполагает выполненным

^* E-mail: byk@astro.ioffe.ru

условие вмороженности магнитного поля в течение

^** E-mail: osm2004@mail.ru

плазмы (в общем случае релятивистское) с безраз-

570

ЖЭТФ, том 161, вып. 4, 2022

Ускорение космических лучей до энергий. . .

нерелятивистскими, так и с Γ_f ≫ 1) требуют учета

мерной скоростью β_f = u_f /c, c — скорость света. В√

нелинейных эффектов обратного влияния на струк-

выражении (1) Γ_f = 1/

1 - β2f. Функция φ(β_f) для

туру течения.

ультрарелятивистского течения с Γ_f ≫ 1 стремится

√

В данной работе мы приводим результаты мо-

к единице, а в случае β_f ≪ 1 имеем φ(β_f ) ∝

β_f

делирования спектров КЛ, ускоренных трансреля-

[3]. Для ультрарелятивистских ударных волн оцен-

тивистскими ударными волнами, в рамках нелиней-

ки максимальных энергий ускоренных КЛ обсужда-

ных моделей для характерных параметров одномер-

лись в работах [19-21].

ных МГД-течений плазмы.

Соотношение (1) показывает, что для заданного

значения «магнитной светимости» L_M максималь-

ные значения энергии ускоренных частиц достига-

2. МОДЕЛИРОВАНИЕ СПЕКТРОВ

ются в трансрелятивистских МГД-течениях с β_f <∼

НЕТЕПЛОВЫХ ЧАСТИЦ В

1, и такие течения могут реализовываться в неко-

ТРАНСРЕЛЯТИВИСТСКОЙ УДАРНОЙ

∼

торых классах сверхновых [22]. Радионаблюдения

ВОЛНЕ МЕТОДОМ PiC

сверхновой SN 2009bb [18] позволили оценить ско-

В данном разделе представлены результаты рас-

рость расширения остатка за первый год со ско-

чета эффективности инжекции частиц в нетепло-

ростью, соответствующей лоренц-фактору течения

вое распределение в зависимости от угла наклона

Γ_f ≈ 1.5, а наблюдаемый эффект синхротронного

бесстолкновительной трансрелятивистской ударной

самопоглощения дает оценку магнитного поля по-

волны в электрон-протонной плазме. Моделирова-

рядка 0.03 Гс на расстоянии 0.1 пк. Высокая на-

ние структуры ударной волны выполнено с помо-

пряженность магнитного поля на расстоянии, пре-

щью PiC-кода Smilei [32]. Инициализация ударной

вышающем 10⁵ радиусов звезды, может указывать

волны выполнена в двумерной области, заполненной

на распространение ударной волны в звездном ветре

плазмой, где заданы следующие граничные условия:

с медленным убыванием азимутальной компоненты

по оси x через правую границу втекает однород-

поля. Оценка, по критерию Хилласа, максимальной

ный поток плазмы с постоянными электрическим

энергии КЛ, ускоренных в таком объекте, соглас-

и магнитным полями; на левой границе установле-

но авторам работы [18], достигает 4000 ПэВ. По-

на отражающая стенка. При столкновении потока

этому представляется важным построение количе-

со стенкой образуется ударная волна. В поперечном

ственных моделей ускорения КЛ трансрелятивист-

направлении поставлены периодические граничные

скими ударными волнами.

условия.

Кинетические методы прямого моделирования

Рассматривались следующие параметры втекаю-

плазмы с помощью кода PiC (particle-in-cell) [23]

щего потока плазмы: лоренц-фактор Γ = 1.5, замаг-

позволяют исследовать как структуру ударных волн

ниченность σ = B²/4πΓρ_bc² = 4 · 10-3 (B — ин-

в бесстолкновительной плазме, так и процессы

дукция магнитного поля, ρ_b — плотность плазмы в

инжекции и ускорения энергичных частиц. Для

системе покоя течения за фронтом ударной волны),

трансрелятивистских ударных волн такие исследо-

температура в безразмерных единицах k_BT/m_pc² =

вания выполнены, в частности, в работах [24-27].

= 10-4 (k_B — постоянная Больцмана, m_p — масса

Мощность современных компьютеров не позволя-

протона). Угол наклона магнитного поля к направ-

ет пока рассчитать методом PiC спектры частиц

лению распространения ударной волны, θ, варьиро-

в диапазоне, превышающем три декады по им-

вался в пределах от 0 до 90^◦. Зависимость всех рас-

пульсу, что недостаточно для моделирования спект-

сматриваемых величин от концентрации масштаби-

ров ПэВ-протонов. Для моделирования ускорения

руется через безразмерные параметры. В качестве

протонов высоких энергий используют упрощенные

параметра, в зависимости от которого будут мас-

уравнения переноса [28] или методы Монте-Карло

штабироваться результаты, удобнее выбрать плаз-

[29]. Эффективное ускорение КЛ с существенной

менную электронную частоту

конверсией энергии в спектры ускоренных частиц

√

может приводить к тому, что давление КЛ [29]

4πn_ee²

модифицирует МГД-течение, натекающее на удар-

ω_e =

Γm_e

ную волну, а амплитуды магнитных флуктуаций

растут за счет анизотропии убегающих ускоренных

где n_e — концентрация электронов в системе по-

надтепловых частиц [30, 31]. Расчеты спектров КЛ,

коя течения за фронтом ударной волны, m_e

—

ускоренных сильными ударными МГД-волнами (как

масса электрона, e — элементарный заряд. Чис-

571

А. М. Быков, С. М. Осипов, В. И. Романский

ЖЭТФ, том 161, вып. 4, 2022

Расчеты показывают, что для высокой замагни-

ченности плазмы распределение протонов за фрон-

том трансрелятивистской ударной волны имеют вы-

раженную нетепловую компоненту для углов накло-

на магнитного поля θ < 40^◦. Детальное микроскопи-

ческое моделирование распределения частиц кине-

тическими методами PiC возможно лишь в ограни-

ченном интервале энергий. Для изучения возмож-

ности ускорения частиц высоких энергий можно до-

полнить модель макроскопическим расчетом с ис-

пользованием метода Монте-Карло, который пред-

полагает использование параметризации процессов

переноса заряженных частиц.

3. РЕЛЯТИВИСТСКАЯ МОДЕЛЬ

МОНТЕ-КАРЛО

Моделирование процесса ускорения частиц

бесстолкновительными ударными волнами мето-

дом Монте-Карло предполагает задание закона

Рис. 1. (В цвете онлайн) Функции распределения протонов

f(E) в системе покоя течения за фронтом ударной волны

рассеяния частиц плазмы МГД-флуктуациями

при разных значениях угла наклона внешнего магнитного

(рассеивающими центрами). В отсутствие куло-

поля к нормали ударной волны, θ. Расчет выполнен ме-

новских столкновений длина свободного пробега

тодом PiC для ударной волны с лоренц-фактором нате-

частиц параметризуется как функция гирорадиуса

кающего потока Γ = 1.5 и параметром замагниченности

частицы в магнитном поле [33]. Метод позволяет

σ = 4 · 10-3

учитывать эффекты обратного влияния давления

ускоренных частиц на структуру сверхзвукового

течения и ударной волны. В модели рассчитывается

полная функция распределения частиц плазмы с

ленные параметры моделирования; шаг по вре-

произвольной степенью анизотропии, что позволяет

мени Δt

= 0.09ω-1e; шаг пространственной сет-

исследовать развитие плазменных неустойчивостей

ки Δx

= 0.2c/ω_e; полное время моделирования

и согласованно рассчитывать длины свободного

45000ω-1e; количество ячеек в области симуляции

пробега частиц и структуру течения с учетом

N_x = 256000, N_y

= 2000. Для экономии вычис-

сильного роста магнитных флуктуаций, вызывае-

лительных ресурсов масса электрона увеличена по

мого токами и градиентом давления ускоренных

сравнению с реальной и равна m_e = m_p/100.

частиц [31].

Эффективность ускорения частиц зависит от уг-

Ниже мы приведем результаты моделирова-

ла наклона магнитного поля. Для ускорения необ-

ния методом Монте-Карло локальной структуры

ходимо, чтобы частицы могли уходить от фронта

трансрелятивистской квазипродольной бесстолк-

ударной волны вверх по течению, двигаясь вдоль

новительной ударной волны (с углами наклона

силовых линий магнитного поля. Для этого угол

меньше критических углов, обсуждавшихся в

наклона магнитного поля в системе отсчета на-

разд. 2) и спектров ускоренных космических лу-

текающего потока должен удовлетворять условию

чей для параметров, соответствующих остатку

cosθ > v_sh/c, где v_sh — скорость ударной волны. При

релятивистской сверхновой.

углах, больших критического, ускорение практиче-

Стационарная и плоскопараллельная модель

ски исчезает. При меньших углах оно нарастает при

Монте-Карло описывает ускорение частиц квази-

увеличении угла от нуля до критического значения.

продольными бесстолкновительными ударными

Графики спектров ускоренных частиц за фронтом

волнами и усиление магнитных полей за счет плаз-

ударной волны на расстояниях 2000c/ω_e от фронта

менных неустойчивостей в предфронте ударной

при различных углах наклона поля приведены на

волны. Ускорение частиц бесстолкновительными

рис. 1.

ударными волнами происходит по механизму Ферми

572

ЖЭТФ, том 161, вып. 4, 2022

Ускорение космических лучей до энергий. . .

первого порядка при их рассеянии на флуктуациях

правлен от фронта ударной волны в далекий пред-

магнитных полей и многократном пересечении

фронт. Таким образом, усиленные волны движутся

фронта ударной волны. В процессе ускорения

в определенном направлении, и мы вводим скорость

значительная доля потока энергии, натекающего

движения рассеивающих центров v_scat (x) относи-

на ударную волну, трансформируется в энергию

тельно фоновой плазмы. Предполагается, что уско-

распределения ускоренных частиц. Поскольку

ренные частицы рассеиваются упруго и изотропно

давление ускоренных частиц в предфронте ударной

в системе покоя, движущейся относительно фрон-

волны сравнимо с потоком энергии, натекаю-

та ударной волны со скоростью u (x) + v_scat (x), а

щим на ударную волну, его градиент приводит

фоновые частицы — в системе покоя, движущейся

к модификации течения в предфронте ударной

относительно фронта ударной волны со скоростью

волны.

u (x), где u (x) — скорость фоновой плазмы.

В модели Монте-Карло мы разделяем частицы

В численной модели Монте-Карло на основе ите-

на фоновые и ускоренные. Ускоренными считают-

рационного процесса достигается выполнение зако-

ся частицы, которые хотя бы раз пересекли фронт

нов сохранения потока импульса и энергии вблизи

ударной волны в направлении, обратном течению

ударной волны. Закон сохранения потока числа час-

плазмы. Ускоренные частицы описываются во всей

тиц выполняется в данной модели автоматически. В

расчетной области в виде частиц и распространя-

процессе итераций подбираются профиль скорости

ются на основе метода Монте-Карло малоуглового

фоновой плазмы u (x) в предфронте ударной волны

рассеяния [34-38]. Фоновая компонента плазмы в

и полное сжатие ударной волной R_tot.

большей части предфронта описывается на основе

гидродинамических уравнений. Вблизи и за фрон-

3.1. Длина свободного пробега частиц

том ударной волны она описывается, как и уско-

В численной модели Монте-Карло частицы при

ренная компонента, в виде частиц, что позволяет не

распространении рассеиваются упруго и изотроп-

вводить дополнительного параметра, определяюще-

но в системах покоя рассеивающих центров. Меж-

го инжекцию частиц в процесс ускорения. Частица

ду рассеяниями частицы движутся прямолинейно

распространяется до тех пор, пока не уйдет из рас-

и равномерно. Расстояние, которое проходит части-

четной области далеко за фронт ударной волны ли-

ца между рассеяниями, пропорционально длине сво-

бо не пересечет границу свободного ухода в далеком

бодного пробега частицы в системе покоя рассеива-

предфронте ударной волны при x = x_FEB. Грани-

ющих центров. Длина свободного пробега является

ца свободного ухода соответствует размерам облас-

функцией импульса частицы и значения спектраль-

ти ускорения частиц в реальном астрофизическом

ной плотности энергии турбулентности W (x, k) в

объекте и определяет максимальный импульс уско-

точке нахождения частицы (k — модуль волнового

ренных частиц.

числа волн, из которых состоят магнитные флукту-

Функция распределения ускоренных частиц име-

ации).

ет значительную анизотропию в предфронте удар-

В данной работе мы используем следующее вы-

ной волны. Анизотропия функции распределения

ражение для длины свободного пробега частиц:

приводит к развитию плазменных неустойчивостей,

(

)-1

которые усиливают магнитные поля в предфрон-

λ(x, p) =

(2)

те. Таким образом, задача об ускорении частиц бес-

λ_B,st (x, p)

λ_ss (x, p)

столкновительной ударной волной является сущест-

венно нелинейной, так как ускорение частиц опреде-

где p — импульс частицы, бомовская длина пробега

ляется флуктуациями магнитного поля, а усиление

в стохастическом крупномасштабном поле

магнитных полей — анизотропией функции распре-

λ_B,st (x, p) =

(3)

деления ускоренных частиц.

eB_ls,st (x, k_res)

В данной работе мы рассматриваем неустойчи-

крупномасштабное турбулентное поле

вости, вызванные током ускоренных частиц относи-

(

тельно фоновой плазмы. В расчеты включены бел-

)

) ∫k

ловская нерезонанская [39] и резонансная неустой-

)

B_ls,st (x, k) =

√4π W (x, k) dk,

(4)

чивости. При развитии резонансой неустойчивости

kmin

усиливаются волны, распространяющиеся только в

направлении тока ускоренных частиц в системе по-

k_min в (4) и k_max в (5) — соответственно минималь-

коя фоновых частиц. Данный ток в предфронте на-

ное и максимальное волновые числа численной сет-

573

А. М. Быков, С. М. Осипов, В. И. Романский

ЖЭТФ, том 161, вып. 4, 2022

ки, длина пробега относительно мелкомасштабных

импульса, уносимый ускоренными частицами через

мод

границу свободного ухода (Q^ESCpx

= F^crpx (x_FEB)),

⎛

⎞

-1

F^crpx (x) — поток импульса ускоренных частиц.

∫

(

Закон сохранения потока энергии имеет вид

(x, k)

λ_ss (x, p) =

⎝ W

dk^⎠

(5)

πe

[

]

kres

γ² (x) β (x)

m_pc²n (x) + Φ_th (x) + Φ_w (x)

(частицы в расчетах предполагаются протонами).

+ Fcren (x) = Fen0 + Q^ESCen,

(10)

Резонансное волновое число k_res определяется

где Fen0 — поток энергии в далеком невозмущен-

следующим соотношением:

ном предфронте, Q^ESCen — поток энергии, уносимый

k_respc

ускоренными частицами через границу свободного

= 1,

(6)

eB_ls (x, k_res)

ухода (Q^ESCen = F^cren (x_FEB)), F^cren

(x) — поток энер-

гии ускоренных частиц.

где крупномасштабное магнитное поле

Отметим что, в нашей геометрии QESCen < 0,

(

)

Q^ESCpx > 0.

) ∫k

)

B_ls (x, k) =

√4π W (x, k) dk + B²⁰,

(7)

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

начальные приближенные профили u (x), v_scat

(x)

kmin

в предфронте ударной волны, величина W (x, k),

B₀ — постоянное продольное магнитное поле.

определяющая длину пробега частиц, и полное

сжатие R_tot. Далее проводится распространение

3.2. Нахождение профилей течения и

частиц с целью получения их функции распределе-

плотности энергии турбулентного

ния в области расчета. За фронтом ударной волны

магнитного поля

все частицы изотропны в системе покоя течения,

соответственно, можно определить показатель

Запишем релятивистские законы сохранения

адиабаты всего газа за фронтом ударной волны,

вблизи ударной волны в стационарном случае.

Γp2, на основе полученной после распространения

Все уравнения записаны в системе покоя ударной

функции распределения всех частиц. Уходящие

волны. Закон сохранения потока частиц

потоки определяются, как

γ (x)β (x)n (x) = Fn0,

(8)

Q^ESCen = F^cren (x_FEB), Q^ESCpx = F^crpx (x_FEB).

где γ (x) — лоренц-фактор течения, β (x) = u (x)/c —

В наших расчетах отрицательные значения коорди-

отношение скорости течения фоновой плазмы к ско-

наты x соответствуют предфронту, положительные

рости света, n (x) — концентрация фоновой плазмы,

значения — области за фронтом ударной волны, ко-

Fn0 — поток частиц в далеком предфронте.

торый находится в точке x = 0, R_tot = β₀/β₂ (индекс

Законы сохранения потока импульса представим

«0» означает величины в далеком невозмущенном

в виде

предфронте, индекс «2» — за фронтом ударной вол-

[

]

ны).

γ² (x) β² (x)

m_pc²n (x) + Φ_th (x) + Φ_w (x)

Для нахождения полного сжатия ударной вол-

+ Pth (x) + P_w (x) + Fcrpx (x) = Fpx0 + Q^ESCpx,

(9)

ной запишем уравнения сохранения потоков (8), (9)

и (10) для значения координаты далеко за фронтом

P_th (x) — давление фоновой плазмы (в области пред-

ударной волны и в далеком предфронте:

фронта, где фоновые частицы описываются гидро-

динамическими уравнениями и предполагается, что

γ₂β₂n₂ = γ₀β₀n₀,

(11)

функция распределения фоновых частиц изотропна

в системе покоя течения), P_w (x) — давление турбу-

[

]

лентности,

γ²²β²²

m_pc²n₂ + Φp2 + Φw2

+Pp2 +Pw2 =

[

]

∫

=γ²⁰β²⁰

m_pc²n₀ + Φth0 + Φw0

Γ_th (x) P_th (x)

P_w (x) =

W (x, k) dk, Φ_th (x) =

Γ_th (x) - 1

+Pth0 +Pw0 +Q^ESCpx,

(12)

(k)

[

]

Γ_th

— показатель адиабаты фоновой плазмы,

γ²²β₂

m_pc²n₂ + Φp2 + Φw2

Φ_w (x) = 3P_w (x), Fpx0 — поток импульса в дале-

[

]

=γ²⁰β₀

m_pc²n₀ + Φth0 + Φw0

+Q^ESCen,

(13)

ком невозмущенном предфронте, QESCpx — поток

574

ЖЭТФ, том 161, вып. 4, 2022

Ускорение космических лучей до энергий. . .

где Pp2 — давление всех частиц за фронтом ударной

где G (x, k) — показатель роста турбулентной энер-

волны, Φp2 = Γp2Pp2/(Γp2 - 1). Мы определяем те-

гии в системе покоя фоновой плазмы за счет плаз-

кущее значение на данной итерации величин Q^ESCen ,

менных неустойчивостей, L (x, k) — диссипация тур-

Q^ESCpx и Γp2 из функции распределения, полученной

булентной энергии, Π(x, k) — поток турбулентной

после распространения частиц (Pw2 также вычисля-

энергии по спектру (турбулентный каскад),

ется на основе приведенного ниже уравнения), по-

√

C^∗k11/2

∂

⁽W (x,k)⁾

этому в уравнениях (11)-(13), остаются только три

Π(x, k) = -

√

W (x, k)

(16)

неизвестных: R_tot, Pp2 и n₂. Решая эти уравнения,

ρ (x)

∂k

k²

находим новое значение R_tot, а для следующей ите-

где ρ (x) — плотность фоновой плазмы,

рации R_tot находится усреднением нового значения

и старого.

C^∗ =

C-3/2Colm,

(17)

Затем находим новый профиль скорости тече-

ния в предфронте по формуле, построенной на ос-

C_Colm

— колмогоровская константа, в расчетах

нове (9),

CColm = 1.6. Выражение (16) получено на осно-

γ_new (x)β_new (x) = γ_old (x) β_old (x) +

ве работы [40]. В нашу работу включены показате-

ли роста токовых неустойчивостей — нерезонансной

(x)

Fpx0 + Q^ESCpx - F^crpx (x) - F^thpx (x) - F^wpx

(14)

белловской и резонансной, аналогично работе [31].

γ₀β₀m_pc²n₀

Ток ускоренных частиц, определяющий показатели

где в правой части выражения стоят величины, по-

роста, вычисляется в системе покоя рассеивающих

лученные после распространения частиц, величина

центров после распространения частиц.

γ_old (x) β_old (x) определяется профилем течения на

Дифференцируя уравнения (9) и (10) по x, ис-

предыдущей итерации,

пользуя закон сохранения (8) и исключая из урав-

F^thpx (x) = γ² (x)β² (x) Φ_th (x) + P_th (x)

нений слагаемое, пропорциональное n (x), получим

следующее соотношение:

)

F^wpx (x) = γ² (x) β² (x) Φ_w (x) + P_w (x)

⁽ dΦ_th (x)

dΦ_w (x)

dF^cren (x)

β (x)

— соответственно потоки импульса фоновой плазмы

и турбулентности.

d (γ (x) β (x))

Подбор профиля течения на основе выражения

(Φ_th (x) + Φ_w (x)) =

(14) в области, где фоновые частицы описывают-

γ (x)

ся в виде частиц, дает хорошие результаты в слу-

)

⁽ dP_th (x)

dP_w (x)

dF^crpx (x)

чае нерелятивистского движения при вычислении

= β (x)

(18)

в этой области потока F^thpx (x) на основе функции

распределения фоновых частиц. В релятивистском

Если считать, что каждая из компонент чувствует

случае, как показано в работе [29], поток импульса

только изменение γ (x) β (x), т. е. изменение потоков

в итерационном процессе сходится к большему зна-

энергии импульса происходит адиабатически из-за

чению, чем в далеком предфронте, при использова-

изменения скорости потока, то мы можем разделить

нии выражения (14) вблизи фронта ударной волны.

уравнение (18) на отдельные адиабатические урав-

Следуя работе [29], мы размываем профиль тече-

нения для компонент:

ния u (x) вблизи фронта ударной волны. Новый про-

филь скорости затем усредняется со старым. Ниже

dΦ_th (x)

d(γ (x) β (x))

β (x)

Φ_th (x) =

в данном разделе мы описываем уравнения, на ос-

γ (x)

нове которых вычисляются величины, входящие в

выражения для F^thpx (x) и F^wpx (x).

dP_th (x)

= β (x)

(19)

Спектр энергии турбулентности, определяющий

поток турбулентности F^wpx (x), находим на основе ре-

шения уравнения

dΦ_w (x)

d (γ (x) β (x))

β (x)

Φ_w (x) =

∂W (x,k)

3 ∂ (γ (x)u(x))

γ (x)

γ (x) u (x)

W (x, k) +

∂x

∂Π(x, k)

dP_w (x)

= G(x,k)W (x,k) - L(x,k),

(15)

= β (x)

(20)

∂k

575

А. М. Быков, С. М. Осипов, В. И. Романский

ЖЭТФ, том 161, вып. 4, 2022

dF^cren (x)

dF^crpx (x)

скоростью u (x)+v_scat (x) относительно ударной вол-

= β (x)

(21)

ны (а фоновых частиц — в системе, движущейся со

Уравнение (20) при этом равно уравнению (15), про-

скоростью u (x)) в предположении, что при мало-

интегрированному по k в отсутствие усиления и дис-

угловом рассеянии (при диффузии в нерелятивист-

сипации. При наличии усиления и диссипации по-

ском случае) уравнение (21) выполняется при про-

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

извольном β (x). Тогда уравнение (21) примет вид

аналогичный уравнению (20):

dF^cren (x)

dF^crpx (x)

= [u (x) + v_scat (x)]

(26)

dΦ_w (x)

d (γ (x) β (x))

β (x)

Φ_w (x) =

γ (x)

∫

dP_w (x)

Усиление волн, связанное с анизотропией уско-

= β (x)

G(x, k)W (x, k) dk -

cγ (x)

ренных частиц, уменьшается при приближении

(k)

плазмы к фронту ударной волны. Вдали от фронта

ударной волны, в области быстрого роста волн,

L (x) ,

(22)

cγ (x)

скорость рассеивающих центров v_scat

= v_ampl

определяет инкремент роста G(x, k) . В области

где диссипационное слагаемое

слабого усиления волн скорость рассеивающих

∫

центров определяется альфвеновской скоростью в

L (x) = L (x, k) dk.

эффективном поле длинноволновых флуктуаций,

(k)

v_scat = v_A,eff . Соответственно, в предфронте удар-

ной волны до точки, где фоновая плазма начинает

В настоящей работе мы используем следующие вы-

описываться в виде частиц, величину v_scat (x)

ражения для диссипации турбулентной энергии:

определяем как наименьшее значение (скорости от-

рицательные) из v_ampl (x) и v_A,eff (x), определенных

L(x, k) = v_Γ (x)

W (x, k) ,

(23)

k_th

как

∫

B_ls (x, k_th)

v_Γ (x) =

√

(24)

G(x, k)W (x, k) dk

4πρ (x)

v_ampl

(x) = -(k)

(27)

√

γ (x) dF^crpx (x)/dx

k_thc

k_BT_th (x)

= 1,

(25)

eB_ls (x, k_th)

где T_th (x)

— температура фоновой плазмы

B_eff

v_A,eff (x) = -

(28)

(P_th (x)

= n(x)k_BTth (x)). В данной модели

^√4πρ,

предполагается, что за фронтом ударной волны

Π(x, k) = 0 и L (x, k) = 0.

√

∫

Далее мы предполагаем, что диссипация турбу-

B_eff (x) =

4π

W (x, k) dk + B²⁰.

(29)

лентной энергии приводит к нагреву фоновой плаз-

(k)

мы. Увеличение потока энергии при усилении тур-

булентных движений и магнитных полей происхо-

В области значительного роста магнитных флук-

дит за счет уменьшения потока энергии ускоренных

туаций имеем |v_scat (x)|

≥ |v_ampl (x)|, благодаря

частиц. Таким образом, для выполнения суммарно-

чему увеличение потока энергии турбулентности

го уравнения (18) мы должны добавить к уравнени-

за счет плазменных неустойчивостей в уравнении

ям для компонент (19) и (21) дополнительные сла-

(22) компенсируется уменьшением потока энергии

гаемые, сокращающие неадиабатические слагаемые

ускоренных частиц в уравнении (26). После точ-

из (22).

ки инжекции фоновых частиц в предфронте име-

Уравнение (21) должно быть дополнено так, что-

ем v_scat (x) = v_ampl (x). За фронтом ударной волны

бы выполнялся закон сохранения энергии при вы-

v_scat (x) = 0. В данной модели мы предполагаем, что

полнении закона сохранения импульса, т. е. уравне-

часть потока энергии ускоренных частиц в области,

ния (18). Этого можно достичь введением скорости

где |v_ampl (x)| < |v_A,eff (x)|, идет на нагрев фоновой

рассеивающих центров, если в уравнении (21) заме-

плазмы. Таким образом, уравнение, определяющее

нить u (x) на u (x) + v_scat (x) и выполнять рассеяние

изменение потока энергии фоновой плазмы, имеет

ускоренных частиц в системе покоя, движущейся со

вид

576

ЖЭТФ, том 161, вып. 4, 2022

Ускорение космических лучей до энергий. . .

dΦ_th (x)

d (γ (x) β (x))

На границе свободного ухода спектр энергии тур-

β (x)

Φ_th (x) =

γ (x)

булентности задан следующей формулой:

dP_th (x)

v_diss (x) d

(x)

= β (x)

W (x_FEB , k) =

)(

)

((

)_ζ

(

)_ζ

-1

L (x) ,

(30)

k_min

k_max

cγ (x)

=A_norm

)

(33)

((

)_α

)((

)_ζ

где v_diss (x)

|v_A,eff (x)| - |vampl (x)|, если

k_en

k_ds

|v_ampl (x)|

|v_A,eff (x)|, а в противоположном

пределе v_diss (x) = 0.

где A_norm — нормировочная константа, ζ₀ = 2, α =

Подставляя в уравнение (30) выражение для Φ_th,

= 11/3, ζ₁

= 1, ζ₂

= 6, k_ds

= 10⁴/rg0, rg0

получаем следующее уравнение для давления фоно-

= m_pcu₀/eB₀ (в области k_en ≪ k ≪ k_ds показатель

вой плазмы:

спектра близок к -5/3). В расчетах

d Pth (x)

k_en

2π

γ (x) β (x)

k_min =

k_en =

L_en = 3 · 10¹⁹ см,

dx Γ_th (x) - 1

10⁶

d (γ (x) β (x)) Γ_th (x) P_th (x)

L_en — энергонесущий масштаб турбулентности. Вы-

Γ_th (x) - 1

ражение (33) устроено таким образом, чтобы как

v_diss (x) d

(x)

W (x_FEB , k), так и, соответственно, поток турбу-

= γ (x)

L (x) .

(31)

лентной энергии (16) при x = x_FEB точно обра-

щались в нуль на границах численной сетки k_min и

Решение уравнения

(31) определяет величину

k_max. Уравнение (15) является уравнением диффу-

F^thpx (x) в выражении (14).

зионного типа, и его решение находится на основе

неявной численной схемы с нулевыми граничными

условиями по k для W (x, k).

3.3. Значения физических параметров

Расстояния в численном коде измеряются в еди-

модели

ницах rg0 = m_pcu₀/eB₀. Отметим, в частности, что

Во всех расчетах данной работы лоренц-факто-

xF EB = -6.42 · 10⁵rg0 для моделей A и AC.

ры ударных волн γ₀ = 1.5, концентрация фоновой

плазмы в далеком предфронте n₀ = 4 · 10-3 см-3,

3.4. Результаты расчетов

граница свободного ухода находится в точке x_FEB =

= -5·10¹⁷ см. В основных расчетах постоянное про-

С помощью рис. 2 мы иллюстрируем сходи-

дольное магнитное поле B₀ = 3·10-6 Гс (обозначим:

мость нашей численной итерационной модели

A — модель без турбулентного каскада и диссипа-

Монте-Карло. На рис. 2 представлены расчеты для

ции турбулентной энергии и AC — модель c вклю-

немодифицированного случая и для модели A. В

чением турбулентного каскада и диссипации турбу-

немодифицированном расчете скорость течения

лентной энергии). Также выполнены расчеты с B₀ =

постоянна в предфронте и за фронтом ударной

= 3·10-5 Гс (обозначим: B — модель без турбулент-

волны и испытывает скачок при переходе через

ного каскада и диссипации турбулентной энергии и

фронт с полным сжатием R_tot = 4. Также в дан-

BC — модель c включением турбулентного каскада и

ном расчете отсутствуют диссипация и усиление

диссипации турбулентной энергии), B₀ = 3 · 10-4 Гс

магнитных полей, а v_scat (x) = 0. Как видно на

(обозначим: C — модель без турбулентного каскада

рис. 2б,в, потоки импульса и энергии вблизи удар-

и диссипации турбулентной энергии и CC — модель

ной волны в случае немодифицированного расчета

c включением турбулентного каскада и диссипации

значительно отличаются от их значений в далеком

турбулентной энергии). В системе также есть слу-

предфронте. Это означает, что в таком случае зако-

чайная компонента магнитного поля, а спектр энер-

ны сохранения не выполняются. В расчете модели A

гии турбулентности нормирован следующим обра-

достигнуто равенство потоков импульса и энергии

зом:

с их значениями в далеком предфронте с помощью

∫

самосогласованного подбора профилей течения и

B²

W (x_FEB, k) dk =

(32)

спектров магнитного поля на основе итерационной

4π

схемы. Таким образом, достигнуто выполнение

kmin

577

А. М. Быков, С. М. Осипов, В. И. Романский

ЖЭТФ, том 161, вып. 4, 2022

Рис. 2. (В цвете онлайн) Зависимости от расстояния до

фронта ударной волны а) величины γ (x)β (x), определяю-

щей скорость течения, б) полного потока импульса Fpx (x)

и в) полного потока энергии Fen (x). Сплошные синие ли-

нии соответствуют немодифицированному расчету, штри-

ховые зеленые — модели A

Рис. 4. (В цвете онлайн) Функции распределения частиц

по импульсам в системе покоя ударной волны. Жирные

линии соответствуют функциям распределения частиц за

фронтом ударной волны, тонкие линии — функциям рас-

пределения частиц на границе свободного ухода. Соответ-

ствие кривых моделям определено в легенде

и vscat (x) определяется значением vampl (x). Ближе

к фронту ударной волны и до точки инжекции фо-

новых частиц усиление за счет плазменных неустой-

(x) определяется

чивостей становится слабее и v_scat

Рис. 3. (В цвете онлайн) Зависимости скорости тече-

значением v_A,eff (x). В данной области пространства

ния u (x) — штриховые кривые, и скорости рассеиваю-

происходит дополнительный нагрев фоновой плаз-

щих центров относительно системы покоя фоновой плаз-

мы из-за наличия v_diss (x) в уравнении (31).

мы, vscat (x) — штрихпунктирные линии, от расстояния до

На рис. 4 приведены функции распределения

фронта ударной волны. Синие линии соответствуют моде-

ли A, зеленые — модели AC

частиц для всех моделей в системе покоя ударной

волны, за фронтом ударной волны и на границе

свободного ухода. Видно, что с ростом амплитуды

законов сохранения в трансрелятивистском течении

спектра турбулентного магнитного поля в далеком

вблизи ударной волны.

предфронте происходит увеличение максимального

На рис. 3 приведены скорость течения фоновой

импульса, который приобретают ускоряющиеся час-

плазмы u (x) и скорость рассеивающих центров от-

тицы. Спектры частиц для одинаковых значений

носительно фоновой плазмы, v_scat (x), для моделей

амплитуды турбулентного магнитного поля в дале-

A и AC. Как видно на рис. 3, на большом расстоя-

ком предфронте слабо различаются. Во всех расче-

нии от фронта ударной волны в сторону предфрон-

тах частицы могут быть ускорены до энергий выше

та величина v_scat (x) имеет максимум. На данных

1 ПэВ.

расстояниях происходит значительное усиление маг-

На рис. 5 приведены спектры турбулентной энер-

нитных полей за счет плазменных неустойчивостей

гии в различных точках вблизи ударной волны для

578

ЖЭТФ, том 161, вып. 4, 2022

Ускорение космических лучей до энергий. . .

ростом магнитного поля в далеком предфронте уси-

ление магнитных полей за счет плазменных неустой-

чивостей слабеет. Значения полного магнитного по-

ля за фронтом ударной волны близки друг к другу

для всех моделей. В расчетах C и СС магнитное поле

за фронтом ударной волны даже ниже, чем в осталь-

ных моделях. Смягчение спектра ускоренных час-

тиц при высоких энергиях (см. рис. 4) в моделях A и

AC связано с максимумом значения v_scat (см. рис. 3),

обусловленного развитием плазменных неустойчи-

востей. Наличие области пространства, где v_scat (x)

по модулю сравнивается с u (x), приводит к смягче-

нию спектра ускоренных частиц, достигающих дан-

ной области пространства. Поскольку, как отмечено

выше, усиление менее существенно в остальных мо-

делях, значение по модулю v_scat (x) меньше, и влия-

ние на спектры ускоренных частиц менее выражено.

Выражение для потока частиц в системе покоя

рассеивающих центров имеет вид

∫1

J_cr(x, p) = 2π

vμf^crsc (x, p, μ) dμ,

(34)

-1

Рис. 5. (В цвете онлайн) Зависимости спектра турбулент-

ной энергии W (x, k) от волнового числа k при различных

где f^crsc (x, p, μ) — функция распределения ускорен-

положениях относительно фронта ударной волны: сплош-

ных частиц в системе покоя рассеивающих центров,

ные линии — x = xF EB, штриховые линии — x = 0.5xF EB ,

μ — косинус угла между направлением импульса

штрихпунктирные линии — за фронтом ударной волны.

частицы и нормалью к фронту ударной волны. По-

Тонкие синие линии соответствуют модели A, жирные зе-

ток частиц, пересекающих границу свободного ухо-

леные линии — модели AC

да в системе покоя рассеивающих центров, приведен

на рис. 7. Вблизи границы свободного ухода система

рассеивающих центров практически совпадает с си-

стемой покоя фоновой плазмы. Таким образом, ухо-

дящий поток частиц на рис. 7 соответствует потоку,

уходящему в межзвездную среду в ее системе покоя.

Распределение частиц по импульсам на рис. 7 сдви-

нуто к большим импульсам по сравнению с функ-

цией распределения частиц на границе свободного

ухода в системе покоя ударной волны (см. рис. 4)

из-за релятивистского относительного движения си-

стем покоя.

Максимальный импульс ускоренных частиц в

Рис. 6. (В цвете онлайн) Профили полного магнитного по-

данной модели определяется соотношением между

ля. На данном рисунке rg0 отвечает значению постоянного

длиной свободного пробега частиц и размером пред-

магнитного поля B0 = 3·10-6 Гс. Обозначение кривых для

различных моделей соответствуют рис. 4

фронта |x_FEB |. Длина свободного пробега ускорен-

ных частиц с наибольшей энергий здесь близка к

так называемому бомовскому пределу (с линейной

моделей A и AC. Пики на спектрах рис. 5 связа-

зависимостью пробега от импульса частицы) в зна-

ны с развитием плазменных неустойчивостей. Об-

чительной части далекого предфронта, где эффект

щий подъем кривых на рис. 5 обусловлен адиабати-

усиления магнитного поля мал. Частицы ускоряют-

ческим усилением магнитного поля.

ся до тех пор, пока не достигают границы свобод-

На рис. 6 приведены профили полного магнит-

ного ухода в далеком предфронте, из-за роста дли-

ного поля для всех моделей. Как видно на рис. 6, с

ны свободного пробега с импульсом. Ударная вол-

579

А. М. Быков, С. М. Осипов, В. И. Романский

ЖЭТФ, том 161, вып. 4, 2022

чений сверхновых типа Ib/c, выше 10³⁶ эрг/с. С уче-

том весьма высокой эффективности передачи энер-

гии ударной волны частицам с энергией порядка

петаэлектронвольта (см. рис. 4 и 7) указанной вы-

ше мощности источников достаточно для форми-

рования высокоэнергичной компоненты галактиче-

ских КЛ непосредственно за изломом в спектре. Для

обеспечения наблюдаемой слабой анизотропии КЛ

в модели с редкими ускорителями необходимо удер-

жание КЛ в протяженном гало Галактики с разме-

рами в несколько килопарсек на промежутках вре-

мени много больших, чем период между вспышками

сверхновых типа Ib/c.

Длина пробега ускоренных частиц для ионов

в рассмотренной модели зависит от их магнитной

жесткости — отношения энергии частицы к ее за-

ряду. Соответственно, максимальные энергии уско-

ренных частиц будут пропорциональны их зарядо-

вому числу. Химический состав вблизи максималь-

Рис. 7. Поток уходящих частиц в системе покоя рассеи-

ных энергий частиц, ускоренных релятивистскими

вающих центров. Соответствие кривых моделям такое же,

сверхновыми, будет смещаться в сторону тяжелых

как на рис. 4

ядер.

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

на с четырех-скоростью Γβ ∼ 1 не должна затор-

мозиться вплоть до радиуса R ≈ 3 · 10¹⁸ см, что

имеет место, если среда, в которую распространяет-

В работе исследованы возможности ускорения

протонов и ядер КЛ до энергий выше 10¹⁵ эВ в

ся ударная волна, достаточно разреженная. Это, в

условиях трансрелятивистских истечений сверхно-

частности, имеет место для релятивистских сверх-

вых типа Ib/c. Методами прямого кинетического мо-

новых типа Ib/c, происходящих от массивных звезд

делирования PiC (particle-in-cell) показано, что ин-

типа Вольфа - Райе с быстрым звездным ветром, ха-

жекция частиц в процессе ускорения частиц тран-

рактеризующимся скоростью потери массы порядка

срелятивистской ударной волной с Γ_f ≈ 1.5 и фор-

10-6M_⊙/год [41]. В расчетах мы предполагаем, что

мирование нетеплового распределения частиц имеет

значение |x_FEB | может составлять 1/5 от текущего

место для углов наклона поля к нормали к фронту

радиуса ударной волны и, таким образом, выбрали

ударной волны, не превышающих 40^◦.

величину |x_FEB| = 5 · 10¹⁷ см, что демонстрирует

возможность ускорения протонов до энергий выше

Для определения максимальных энергий прото-

петаэлектронвольта в такой системе.

нов, ускоренных трансрелятивистской ударной вол-

По результатам обзора Ликской обсерватории

ной, выполнено нелинейное моделирование ускоре-

сверхновые типа Ib/c составляют примерно 4 % от

ния КЛ методом Монте-Карло. Модель учитывает

совокупности сверхновых всех типов [42]. Реляти-

процессы усиления магнитных флуктуаций неусто-

вистские сверхновые встречаются с частотой около

чивостями анизотропного распределения ускорен-

0.7 % от сверхновых типа Ib/c [18,41], что дает (весь-

ных частиц и модификацию трансрелятивистского

ма грубую) оценку их частоты в Галактике порядка

течения плазмы давлением ускоренных частиц, рас-

одного события в 50 000 лет. Столь редкие собы-

пространяющихся вверх по течению. Показана воз-

тия не могут вносить существенный вклад в уско-

можность ускорения протонов и ядер до энергий,

рение основной массы галактических КЛ с энерги-

заметно превышающих петаэлектронвольт в тран-

ями от гигаэлектронвольт до 100 ТэВ, однако они

срелятивистских течениях сверхновых звезд. Мак-

могут быть важными источниками КЛ с энергиями

симальные энергии ускоренных протонов зависят

в диапазоне петаэлектронвольт. Средняя мощность,

от амплитуды флуктуирующего магнитного поля в

вкладываемая в Галактику от релятивистских исте-

околозвездной среде (см. рис. 4).

580

ЖЭТФ, том 161, вып. 4, 2022

Ускорение космических лучей до энергий. . .

Сверхновые звезды с трансрелятивистскими те-

14.

А. М. Быков, УФН 188, 894 (2018) [A. M. Bykov,

чениями, составляющие, вероятно, относительно

Phys.-Usp. 61, 805 (2018)].

редкую популяцию в Галактике, тем не менее могут

15.

A. M. Bykov, D. C. Ellison, A. Marcowith, and

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

S. M. Osipov, Space Sci. Rev. 214, 41 (2018).

КЛ с энергиями выше области излома.

16.

X.-Y. Wang, S. Razzaque, P. Mészáros, and

Z.-G. Dai, Phys. Rev. D 76, 083009 (2007).

Благодарности. Авторы признательны рецен-

зенту статьи за конструктивные замечания.

17.

R. Budnik, B. Katz, A. MacFadyen et al., Astrophys.

J. 673, 928 (2008).

Финансирование. Работа поддержана Российс-

18.

S. Chakraborti, A. Ray, A. M. Soderberg et al.,

ким научным фондом (грант № 21-72-20020). Резуль-

Nature Comm. 2, 175 (2011).

таты работы получены с использованием вычис-

лительных ресурсов суперкомпьютерных центров

19.

M. Lemoine and G. Pelletier, MNRAS 402, 321

МСЦ РАН и Санкт-Петербургского политехниче-

(2010).

ского университета Петра Великого (scc.spbstu.ru).

20.

B. Reville and A. Bell, MNRAS 439, 2050 (2014).

21.

A. R. Bell, A. T.Araudo, K. M., J. H. Matthews, and

K. M. Blundell, MNRAS 471, 2364 (2018).

ЛИТЕРАТУРА

22.

R. A. Chevalier, Astrophys. J. 499, 810 (1998).

В. Л. Гинзбург, С. И. Сыроватский, Происхожде-

ние космических лучей, Изд-во АН СССР, Москва

23.

M. Pohl, M. Hoshino, and J. Niemiec, Progr. Part.

(1963).

Nucl. Phys. 111, 103751 (2020).

A. M. Hillas, J. Phys. G 31, 95 (2005).

24.

V. I. Romansky, A. M. Bykov, and S. M. Osipov, J.

Phys.: Conf. Ser. 1038, 012022 (2018).

M. Lemoine and E. Waxman, JCAP 11, 009 (2009).

25.

V. I. Romansky, A. M. Bykov, and S. M. Osipov, J.

К. В. Птицына, С. В. Троицкий, УФН 180, 723

Phys.: Conf. Ser. 1040, 022005 (2019).

(2010) [K. V. Ptitsyna and S. V. Troitsky, Phys.-Usp.

53, 691 (2010)].

26.

P. Crumley, D. Caprioli, S. Markoff, and A. Spitkov-

sky, MNRAS 485, 5105 (2019).

J. H. Matthews, A. R. Bell, K. M. Blundell, and

A. T. Araudo, MNRAS 482, 4303 (2019).

27.

A. Ligorini, J. Niemiec, O. Kobzar et al., MNRAS

502, 5065 (2021).

Г. Т. Зацепин, А. Е. Чудаков, ЖЭТФ 41, 655

(1961) [G. T. Zatsepin and A. E. Chudakov, Sov.

28.

Y. Nagar and U. Keshet, MNRAS 501, 329 (2021).

Phys. JETP 14, 469 (1962)].

29.

D. C. Ellison, D. C.Warren, and A. M. Bykov,

F. Aharonian, A. G. Akhperjanian, A. R. Bazer-Bachi

Astrophys. J. 776, 46 (2013).

et al., Astron. Astrophys. 449, 223 (2006).

30.

A. M. Bykov, A. Brandenburg, M. A. Malkov, and

F. Acero, M. Lemoine-Goumard, M. Renaud et al.,

S. M. Osipov, Space Sci. Rev. 178, 201 (2013).

Astron. Astrophys. 580, 74 (2015).

31.

A. M. Bykov, D. C. Ellison, S. M. Osipov, and

S. Funk, Ann. Rev. Nucl. Part. Sci. 65, 245 (2015).

A. E. Vladimirov, Astrophys. J. 789, 137 (2014).

10.

Е. Г. Бережко, В. К. Елшин, Л. Т. Ксенофонтов,

32.

J. Derouillat, A. Beck, F. Perez et al., Comput. Phys.

ЖЭТФ 109, 3 (1996) [E. G. Berezhko, V. K. Elshin,

Comm. 222, 351 (2018).

and L. T. Ksenofontov, JETP 82, 1 (1996)].

11.

D. C. Ellison, P. Slane, D. J. Patnaude et al.,

33.

D. C. Ellison and D. C.Eichler, Astrophys. J. 286,

Astrophys. J. 744, 39 (2012).

691 (1984).

12.

V. Ptuskin, V. Zirakashvili, and Eun-Suk Seo,

34.

F. C. Jones and D. C. Ellison, Space Sci. Rev. 58,

Astrophys. J. 763, 47 (2013).

259 (1991).

13.

N. M. Budnev, A. Chiavassa, O. A. Gress et al.,

35.

D. C. Ellison, M. G. Baring, and F. C. Jones,

Astropart. Phys. 117, 102406 (2020).

Astrophys. J. 473, 1029 (1996).

581