Расчет скорости ионизации вещества атмосферы земли протонами галактических и солнечных космических лучей

Одной из важных задач в физике космических лучей (КЛ) является исследование распространения элементарных частиц в атмосфере Земли, их влияния на окружающую среду и различные системы через последовательные процессы взаимодействий. Поэтому на протяжении многих лет остается актуальной оценка скорости ионизации, индуцированной КЛ [Bazilevskaya et al., 2008; Usoskin et al., 2009; Velinov et al., 2017]. Поскольку прогресс в вычислительных технологиях шагнул далеко вперед, эти исследования проводятся не только традиционными экспериментальными методами, но и при помощи численного моделирования. На станции КЛ «Апатиты» для этих целей, помимо уже имеющейся системы мониторинга, был разработан программный комплекс RUSCOSMICS [Maurchev, Balabin, 2016; ] на базе набора инструментов GEANT4 [Agostinelli et al., 2003; Maurchev et al., 2015]. Развиваемая методика даст возможность проводить оценку дозы радиации на разных высотах в атмосфере Земли с использованием данных стандартных, недорогих и непрерывных измерений, получаемых на наземных станциях нейтронных мониторов. Это позволит обеспечить радиационную безопасность во время полетов на пассажирских лайнерах, особенно при трансатлантических перелетах. В представленной работе путем моделирования прохождения протонов через атмосферу Земли проводится исследование взаимодействия этого вида частиц как галактических (ГКЛ), так и солнечных (СКЛ) космических лучей во время событий наземного возрастания (ground level enhancement, GLE) интенсивности КЛ GLE65 и GLE67 (см. таблицу). Получены численные характеристики интенсивности потоков различных компонент вторичных КЛ (мюоны, протоны, нейтроны, электроны, гамма-кванты), а также рассчитаны скорости ионизации для множества значений жесткости геомагнитного обрезания. Некоторые результаты использования спектров ГКЛ уже были описаны в работе [Maurchev, Balabin, 2016] и поэтому здесь приведены не будут. Верификация модели проводилась сравнением высотных профилей скорости счета, полученных во время запусков шаров-зондов и полетов на аэробусе и рассчитанных путем моделирования прохождения протонов ГКЛ через атмосферу Земли.

МЕТОДИКА

На начальном этапе необходимо определить геометрию модели, для которой будут проводиться расчеты. Делается это при помощи параметризации таким образом, чтобы полученная конфигурация была не только максимально приближена к реальным значениям физических параметров вещества атмосферы Земли, но и оптимизирована для расчетов. В основе метода лежит понятие плоской геометрии, когда столб воздуха делится на N слоев, для которых определяются масса, плотность, температура, а также процентное содержание химических элементов по модели NRLMSISE-00 [Picone et al., 2002] с их последующим усреднением. На практике было получено, что для удовлетворительного конечного результата, расходящегося с данными для верификации в пределах 15 %, оптимальное значение N=20, т. е. масса каждого слоя составляет 5 % от общей массы столба. Алгоритм программы, которая производит этот расчет по табличным данным, представлен на рис. 1. Этот процесс повторяется N– 1 раз, где N — количество слоев атмосферы, а на N расчет проводится простым усреднением оставшихся значений. Показатели процентного содержания H, He, O, N и Ar вычисляются так же, как и в случае с плотностью, по уже известным данным.

Параметры дифференциальных энергетических спектров солнечных протонов в событиях GLE65 и GLE67

GLE №	Дата	Параметры энергетических спектров
		J 0	E 0	J 1	γ
65	28.10.2003	1.2·104	0.60	1.5·104	4.4
67	02.11.2003	4.6·104	0.51	9.7·103	6.3

Таким образом, описываемый подход позволяет оптимизировать расчеты путем изменения процентного значения массы в зависимости от требуемой точности. Параметры ( T , ρ, O % , N % , He % , Ar % ), полученные для 67°34′03″ N, 33°23′36″ E (Апатиты), приведены на рис. 2, 3. Следует заметить, что для других

Рис. 1. Алгоритм расчета значений высоты и плотности для одного слоя атмосферы в модели прохождения протонов КЛ: m _layer — расчетная масса столба атмосферы из табличных значений NRLMSISE; m = M / N — масса слоя атмосферы, M — масса всего столба атмосферы; h _layer, h _nrlm — высота слоя и высота для табличного значения (100 м) соответственно

Рис. 2. Графическое представление параметров модели атмосферы Земли (плотность и температура), используемых для расчета прохождения через нее протонов КЛ и исследования каскадов вторичных частиц

географических координат значения плотности будут отличаться от приведенных, пример такой зависимости представлен на сайте проекта []. Подробнее с возможностями модели NRLMSISE можно ознакомиться из описания [Picone et al., 2002].

На втором этапе определяется модельный источник протонов. Его основными параметрами являются положение в пространстве и плотность вероятности генерируемых частиц, которая соответствует дифференциальному энергетическому спектру. Исходя из анализа результатов проведенных ранее вычислений, было принято решение использовать точечную реализацию генератора, т. е. на верхней границе столба атмосферы в центральной части источник первичных частиц задается без назначения ему пространственного распределения. Такой подход обеспечивает хорошую статистическую точность при сборе информации об интенсивности потоков частиц вторичных КЛ. Для задания энергетических характеристик протонов ГКЛ (в с^–1м^–2ср^–1МэВ^–1) использовалась формула [ГОСТ 25645.122-85] .

DE ^α

( 10 - 2 E + B ) 4

где F ( E ) — поток частиц; E — энергия первичных протонов; D = 5.2; B = 6.2; α = 1.4; C = 2 , что соответствует минимуму солнечной активности). Спектр протонов СКЛ описывается при помощи модели, разработанной в ПГИ, Апатиты [Vashenyuk et al., 2008, 2011; Perez-Peraza et al., 2006, 2008] , согласно которой в спектре присутствуют две компоненты — быстрая (PC) и медленная (DC), их энергетические зависимости могут быть выражены (в с^–1м^–2ср^–1ГэВ^–1) как

J PC = J 0 exp ( –E/E 0 );

J _DC = J ₁ E^{– γ},

где J 0 , E 0 , J 1 , γ — коэффициенты, различающиеся в зависимости от номера события GLE. Для возрастаний с номерами 65 и 67 были использованы параметры, представленные в таблице. Соответствующие

Сalculating the ionization rate… дифференциальные спектры показаны на рис. 4. Видно (рис. 4), что значения для ГКЛ и СКЛ различаются вплоть до нескольких порядков для энергии до 5 ГэВ, но, начиная с 10 ГэВ, потоки уменьшаются и суммарные спектры становятся абсолютно одинаковыми.

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

РЕЗУЛЬТАТЫ

В результате моделирования прохождения протонов ГКЛ и СКЛ через атмосферу Земли были получены большие массивы, включающие в себя энергетические спектры вторичных частиц и их суммарные интенсивности потоков. На рис. 5 представлены типовые распределения в случае только ГКЛ и для GLE65, когда ГКЛ суммируются с медленной компонентой СКЛ. Хорошо видно (рис. 5), что величины потоков сильно различаются до уровня 400 г см–2 (~7.5 км над уровнем моря), после чего, в силу характера спектра первичных частиц, становятся практически одинаковыми.

Особый интерес представляют также данные о скорости ионизации вещества атмосферы Земли частицами КЛ. Модельные детекторы учитывают все виды ионизации, сохраняя конечный результат в формате высотных профилей. Профили, полученные для значений жесткости геомагнитного обрезания 0.65 и 2 ГВ протонов ГКЛ и СКЛ, показаны на рис. 6.

Помимо функций, о которых было сказано выше, модуль RUSCOSMICS для моделирования прохождения частиц КЛ через атмосферу Земли способен рассчитывать интенсивности потоков и скорости ионизации для всех географических координат. На рис. 7 представлены результаты, полученные при использовании в качестве входных параметров спектра протонов ГКЛ и значений жесткости геомагнитного обрезания, выведенных с помощью IGRF [Thébault et al., 2015], шаг сетки составляет 5°. Подробные наборы данных для высот от 1 до 80 км доступны в виде как рисунков, так и архивов на сайте проекта [].

ВЕРИФИКАЦИЯ ВЫЧИСЛЕНИЙ

Для оценки корректности проведенных вычислений в представленной работе используются различные наборы экспериментальных данных. Результаты ранее проведенных сравнений описаны в работе [Маурчев и др., 2018] . Для верификации используются данные как шаров-зондов [Stozhkov et al., 2009] , так и разработанного с целью расширения географии измерений портативного детектора ионизирующего излучения (рис. 8). При помощи этого устройства был записан профиль скорости счета

Рис. 4. Дифференциальные энергетические спектры протонов солнечных КЛ в диапазоне от 100 МэВ до 100 ГэВ для GLE65 ( а ) и GLE67 ( б )

Рис. 3. Процентные соотношения элементов, входящих в состав атмосферы Земли, полученные путем усреднения. Видно, что эти параметры слабо зависят от высоты, оставаясь практически неизменными до 60 км над уровнем моря

во время полета на аэробусе по маршруту Мурманск — Апатиты. Предварительно восемь счетчиков СТС-6 были заменены на более компактные четыре счетчика СТС-5, что позволило значительно уменьшить габариты оборудования. Автономное питание осуществлялось от Li-Ion-аккумулятора 5 А. Регистрируемый счетчиками Гейгера общий поток частиц будет включать вклады заряженных частиц и ~1 % гамма-квантов:

J общ = J p + J _e+e- + J _ц+ц- + 0.01 J Y ,            (4)

ee μμ где Jp — суммарный поток протонов; Je+e– — суммарный поток электронов и позитронов; Jµ+µ– — суммарный поток мюонов; Jγ — суммарный поток

Рис. 5. Энергетические спектры электронов вторичных ГКЛ ( а ) и ГКЛ+СКЛ для DC-компоненты ( б ) на различных глубинах атмосферы (в скобках указана эквивалентная высота над уровнем моря)

Расчет скорости ионизации вещества атмосферы Земли… гамма-квантов. На рис. 9 приводится сравнение экспериментальных данных, полученных при запуске шаров-зондов (Апатиты) для нескольких периодов, соответствующих минимуму солнечной активности, а также профилей с мобильного детектора, с результатами, полученными моделированием. Уже на текущем этапе наблюдается хорошее соответствие; некоторая разница может объясняться отсутствием в модели учета частиц электронов ядер гелия и элементов с Z > 2. Для профилей, полученных при измерениях на самолете, также сказывается различие периодов, для которых эти данные получены. В настоящее время осуществляется работа над устранением указанных недостатков модуля RUSCOSMICS.

ЗАКЛЮЧЕНИЕ

Были проведены расчеты прохождения протонов КЛ через атмосферу Земли с применением внесенных в модуль RUSCOSMICS улучшений. Благодаря оптимизации, достигнутой путем усреднения

Сalculating the ionization rate… физических параметров при описании геометрии модели и, соответственно, уменьшения количества логических объемов, удалось получить хорошее соотношение между скоростью вычислений и корректностью результатов, которая определяется сравнением с наборами экспериментальных данных. При этом авторы не останавливаются на использовании только одного вида детекторов. В настоящее время ведутся работы по реализации возможности моделирования с использованием параллельных вычислений. Полученные данные могут как иметь прикладное значение (например, при оценке эквивалентной дозы), так и использоваться в фундаментальных исследованиях (например, при детальном изучении характеристик каскадов). Все изменения релизов, их описания и архивы полученных данных находятся в свободном доступе на сайте проекта [].

Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 18-3200626.

Рис. 6. Высотные профили скорости ионизации, полученные как результат моделирования прохождения через атмосферу Земли протонов первичных КЛ во время GLE65 ( а , в ) и GLE67 ( б , г ). Значение жесткости геомагнитного обрезания 0.65 ГВ ( а , б ) и 2 ГВ ( в , г )

Рис. 7. Значения скорости ионизации, полученные в результате моделирования прохождения первичных ГКЛ через атмосферу Земли и образования в ней каскадов вторичных КЛ с использованием глобальной геометрии (шаг сетки 5°). Данные, представленные на графиках, соответствуют высотам 1 км ( а ) и 10 км ( б ) над уровнем моря

Рис. 8. Блок-схема мобильного детектора заряженной компоненты КЛ, разработанного и введенного в эксплуатацию

на станции КЛ «Апатиты»

Рис. 9. Сравнение высотных профилей скорости счета, полученных при запусках шаров-зондов для разных периодов времени ( а ) и измерениях на самолете ( а, б ), с результатами моделирования прохождения протонов ГКЛ для двух пунктов: Апатиты ( а ), Москва ( б )