Расчет жесткости геомагнитного обрезания с помощью трейсинга на основе метода Бунемана−Бориса
- Авторы: Кручинин П.А.1, Малахов В.В.1, Голубков В.С.1, Майоров А.Г.1
-
Учреждения:
- Московский инженерно-физический институт
- Выпуск: Том 64, № 6 (2024)
- Страницы: 750-759
- Раздел: Статьи
- URL: https://journals.rcsi.science/0016-7940/article/view/283259
- DOI: https://doi.org/10.31857/S0016794024060032
- EDN: https://elibrary.ru/QOQMPN
- ID: 283259
Цитировать
Полный текст
Аннотация
Работа включает в себя разработку метода определения жесткости геомагнитного обрезания, основанного на трассировке заряженных частиц в магнитном поле Земли по методу частица-в-ячейке, реализованного в схеме Бунемана–Бориса. Для тестирования метода были проведены расчеты жесткости геомагнитного обрезания в поле идеального диполя и в поле, заданном моделью IGRF. В первом случае полученные данные сопоставлялись с аналитическими значениями. Точность расчета в этом случае составила 3 МВ. Во втором случае была воспроизведена картина полутени в различных географических точках, за различные периоды, а также исследована устойчивость метода к малым возмущениям начальных параметров. В качестве основных результатов в работе построены и проанализированы карты жесткости геомагнитного обрезания на высотах низкоорбитальных спутников для разных направлений в пространстве, а также их вариации с 1900 по 2015 года.
Полный текст
ВВЕДЕНИЕ
Расчет жесткости геомагнитного обрезания (ЖГО) на сегодняшний момент является важной как теоретической, так и практической задачей. Целями могут быть оценка радиационного фона вблизи Земли, расчеты отклика нейтронных мониторов [Koldobskiy et al., 2019; Mishev et al., 2020; Mishev and Poluianov, 2021], моделирование механизма CRAND [Selesnick et al., 2007; Sarkar and Roy, 2022] и другие. Отдельной задачей является изучение вариаций ЖГО во время солнечных событий [Adriani, 2016; Птицына и др., 2021; Данилова и др., 2023], также важны и географические карты ЖГО в магнитоспокойное время и их долгопериодические вариации. Карты в большинстве случаев строятся лишь для вертикальной жесткости [Kress et al., 2015; Smart and Shea, 2019], тогда как для многих задач важна и более детальная картина. В настоящее время существует достаточно много методик, позволяющие решать важные практические задачи, связанные с траекторными расчетами в различных моделях магнитного поля (https://www.geomagsphere.org; https://tools.izmiran.ru/; http://cosmos.hwr.arizona.edu/Util/rigidity.php;). Большинство использует метод Рунге–Кутта для решения уравнения движения частицы в магнитном поле, в то время как существуют и другие схемы для решения этой задачи [Голубков и Майоров, 2021]. Кроме того, сама форма реализации таких методик делает крайне трудоемкой задачу массовых расчетов ЖГО (построение подробных карт, их вариаций во времени). Поэтому в качестве среды для анализа мы выбрали разрабатываемую в данный момент среду для моделирования “жизни” заряженной частицы в околоземном пространстве [https://spacephysics.mephi.ru/beta/GTsimulation/GTsimulation\_8p0.pdf], призванной воспроизвести различные эффекты в этой области и промоделировать потоки заряженных частиц. Основой этой среды является алгоритм трейсинга заряженных частиц, основанный на методе частица–в–ячейке, реализованный в схеме Бунемана–Бориса [Boris, 1970; Boris, 1971]. Данный метод позволяет решать уравнения движения в магнитном поле с сохранением кинетической энергии при длительном движении в сложных магнитных полях и, в частности, в магнитном поле Земли, благодаря чему, метод оказывается точнее и устойчивее, чем метода Рунге–Кутта, и быстрее, чем другие схемы с сохранением энергии [Mao and Wirz, 2011; Qin et al., 2013].
Целью работы является, во-первых, тестирование алгоритма трейсинга и в целом разрабатываемой программной среды, воспроизведение известных эффектов, а во-вторых, восстановление более полной картины ЖГО, в частности получение карты для разных направлений в пространстве на высотах низкоорбитальных спутников.
В качестве тестов рассчитывались значения ЖГО в заранее заданном дипольном поле и сравнивались с соответствующими аналитически вычисленными значениями; проверялась стабильность метода при применении его в реалистичном магнитном поле Земли при небольшом варьировании начальных параметров прилета частиц (широты, долготы, высоты, зенитного и азимутального угла). Расчеты проводились для протонов в условиях магнитоспокойного времени для разных географических точек, высот и направлений, которые затем сопоставлялись с предыдущими аналогичными результатами.
МЕТОДОЛОГИЯ ИССЛЕДОВАНИЯ
Для расчета ЖГО, для выбранной точки пространства и выбранного направления проводился обратный трейсинг заряженной частиц в магнитном поле до тех пор, пока не было выполнено одно из следующих условий остановок:
а) частица достигла поверхности Земли (траектория альбедо);
б) частица достигла расстояния, равного 30 радиусам Земли (галактическая траектория);
в) время жизни частицы от начала трейсинга достигло 100 секунд, при этом ни одно из предыдущих условий не выполнено.
Первый и третий тип траектории соответствует частице с жесткостью ниже пороговой и с точки зрения теории Штёрмера [Stoёrmer, 1955] такие траектории попадают в запрещенную область. Ко второму типу относится частица с жесткостью выше пороговой, что соответствует разрешенной области.
Для проверки метода использовалась модель идеального дипольного поля и модель главного поля IGRF-13 [Alken et al., 2021]. ЖГО определялась подбором жесткости трассируемой частицы в области теоретически рассчитанного значения методом последовательных приближений до точности 0.1 МВ для дипольного поля и шагом 2 МВ для модели главного поля IGRF. Значения шага подобраны исходя из наиболее оптимального соотношения скорости работы алгоритма и точности результата.
ТЕСТИРОВАНИЕ МЕТОДА
3.1. Дипольное поле
На данном этапе задавалось поле идеального диполя с магнитным моментом равным магнитному моменту Земли в 2000 году. Для определения ЖГО использовался метод последовательных приближений. Он заключается в следующем: начальные жесткости выбираются исходя из рассчитанной по формуле Штермера [D.F. Smart and M.A. Shea., 2005], одна (R1) ниже и одна (R2) выше нее на некоторое значение ΔR0, выбранные таким образом, что верхнее значение соответствует заведомо разрешенной траектории, а нижнее – запрещенной. Для этих жесткостей рассчитывается среднеарифметическое значение R0. Далее проводится бэктрейсинг с соответствующей жесткостью и определяется тип траектории. Если жесткость попадает в разрешенную область, т.е. частица покидает пределы магнитосферы, то максимальному значению жесткости присваивается текущее значение, а разница между ними устанавливается ΔR1 = R0 – R1. Если запрещенной, то минимальной жесткости присваивается текущее значение, в этом случае ΔR1 = R2 – R0. Считается среднее между новыми значениями. Процедура повторяется до тех пор, пока ΔRi не станет меньше 0.1 МВ.
В данном случае расчеты проводились для высоты 5 радиусов Земли и диапазона широт λ = [0°; 70°] с шагом в 10°; азимутальный угол ξ менялся в пределах [0°;360°] с шагом в 10°; зенитный угол ε − в пределах [0°;90°] с шагом в 10°.
В результате были получены распределения относительных и абсолютных отклонений жесткости геомагнитного обрезания от теоретического значения по зенитному и азимутальному углам. На рис. 1 для примера представлено такое распределение для широты 30°. Видно, что распределение ошибок равномерное, без выделенных направлений, но в среднем не превышает 1.5%. С учетом этого в качестве единой характеристики точности ЖГО на конкретной широте было использовано средние отклонение по всем направлениям. Они были вычислены для ряда широт в указанном диапазоне, график зависимости среднего относительного отклонения от широты приведен на рис. 2.
Рис. 1. Карта отклонений ЖГО смоделированной частицы относительно значения, полученного по формуле Штёрмера для геомагнитной широты λ = 30°. Расстояние 5 радиусов Земли от центра диполя.
Рис. 2. Зависимость от широты относительной разницы между рассчитанным и теоретическим значениями ЖГО. Расстояние: 5 радиусов Земли от центра диполя.
Сравнение с теоретическими расчетами показало хорошее совпадение в экваториальной области и в области средних широт, заметно ухудшаясь в полярных областях. В районе 70° относительная ошибка составила 28%, при этом отклонение по абсолютному значению для всех широт оказалось одинаковым и лежало в пределах 1–3 МВ.
3.2. Поле IGRF
Основной задачей тестирования в поле IGRF было воспроизведение характерной картины полутени (рис. 3), где нет четкой границы, а диапазоны разрешенных и запрещенных жесткостей чередуются. Для сопоставления были выбраны несколько точек, для которых ранее проводились расчеты Смартом и Ши [D.F. Smart and M.A. Shea.,1994; D.F. Smart and M.A. Shea.,2005; D.F. Smart and M.A. Shea.,2009].
Рис. 3. Иллюстрация полутени геомагнитного обрезания. Sioux Falls. Модель IGRF-13 (слева). 1975 год. Высота 30 км. Зенитный угол прилета частицы 5°. Справа – результат Смарта и Ши [D.F. Smart and M.A. Shea, 2009]. Белым обозначены разрешенные области, черным – запрещенные.
Сравнение проводилось для различных географических точек (Palestine, Sioux Falls, Cape Giradeau, MT. Washington, Newark, Climax), высот и направлений прилета частиц. На рис. 3 представлена картина полутени геомагнитного обрезания. Черная область является запрещенной, а белая – разрешенной. По оси абсцисс отложены азимутальные углы, каждому значению которого соответствует своя полутень. Азимутальные углы расположены так, чтобы направление прилета частицы менялось с запада на восток. Можно выделить несколько устойчивых зон: широкая запрещенная область, узкая запрещенная область у верхней границы полутени. Также выделяется область шума, где каждая тонкая линия запрещенных траекторий при варьировании начальных параметров не совпадают между собой. Для каждой точки наблюдается общее совпадение структуры расположения разрешенных и запрещенных зон, сдвинутой однако на 100–200 МВ между полученными данными и результатами Смарта и Ши.
3.3. Устойчивость полутени жесткости геомагнитного обрезания
Для дополнительной проверки полученного результата была исследована устойчивость картины полутени ЖГО при варьировании начальных параметров прилета частиц: зенитных и азимутальных углов, под которыми частица достигает заданной точки; географических координат; высоты (рис. 4).
Рис. 4. Исследование устойчивости полутени геомагнитного обрезания от начальных параметров прилета частиц.
Данные иллюстрации показывают, что варьирование параметров либо не влияет, либо приводит к плавному изменению основных особенностей полутени (кроме хаотично меняющихся запрещенных и разрешенных зон между верхней запрещенной широкой полосой и основной запрещенной областью снизу), что говорит об устойчивости метода.
КАРТА ЖЕСТКОСТИ ГЕОМАГНИТНОГО ОБРЕЗАНИЯ И ЕЕ ВАРИАЦИИ
Поскольку в поле IGRF не существует точной границы между запрещенной и разрешенной областью, то обычно используется некоторое условное значение, лежащее в пределах области полутени. Мы будем использовать в качестве ЖГО значение, вычисляемое по следующей формуле [Тясто и др., 2015; http://cosray.unibe.ch/~laurent/planetocosmics/doc/planetocosmics\_sum.pdf]:
, (1)
где Rmax – верхняя граница полутени; nallowed – число значений начальных жесткостей, при которых траектория попала в разрешенную область; ΔR – шаг по жесткости.
На рис. 5 показаны карты ЖГО для четырех зенитных углов (0, 15, 30), при фиксированном азимутальном (равным 200°) на высоте 500 км на 5 июля 2006 года (азимутальный угол считается от направления на геомагнитный север). Из карт видно, что для геомагнитного экватора с увеличением зенитного угла картина становится все более ассиметричной вдоль геомагнитной долготы с ярко выраженным минимумом в области ЮАА и максимумом в противоположном полушарии. Кроме того, увеличивается и максимальное значение ЖГО с 16 до 20 ГВ.
Рис. 5. Карта жесткости геомагнитного обрезания для: (а) вертикального направления прилета частиц; (б) зенитный угол 15°, азимутальный угол 200°; (в) зенитный угол 30°, азимутальный угол 200°. Модель IGRF-13. 2006 год. Высота 500 км.
На рис. 6 приведены карты ЖГО для той же высоты и того же времени, но с фиксированным зенитным углом (30°) и разными азимутальными углами от 0 до 350 с шагом 50°. Здесь видно, что при вращении вектора направления прилета по часовой стрелке максимум ЖГО смещается в западном направлении и совершает один оборот вокруг Земли, что хорошо согласуется с эффектом восточно-западной ассиметрии.
Рис. 6. Карты ЖГО для азимутальных углов от 0° до 350° с шагом 50°. Зенитный угол = 30°. 2006 год. Высота 500 км.
4.1. Временные вариации
Для восстановления вековых вариаций было построено 7 карт вертикальной ЖГО с 1903 по 2023 года с шагом в 20 лет на высоте 500 км. Для описания магнитного поля в каждом из временных отрезков также использовалась модель IGRF.
На рис. 7 карта относительных отклонений между значениями ЖГО в 1903 и 2023 годах. На рис. 7а можно видеть, что наибольшая величина вариаций приходится на область Южно-Атлантической магнитной аномалии. Отклонения в данном регионе составляют порядка 101% и по абсолютному значению – 3 ГВ. Это хорошо согласуется с другими результатами. Например, на рис. 7б показана зависимость ЖГО (кресты) в точке (–40°; –50°), соответствующей центру и напряженности магнитного поля (квадраты) [Gelvam et al., 2009] от времени. Видно, что зависимости хорошо коррелируют.
Рис. 7. (а) Карта относительных отклонений значений жесткости геомагнитного обрезания в 1903 и 2023 годах; (б) зависимость от времени жесткости геомагнитного обрезания (кресты) и напряженности магнитного поля (квадраты) [Gelvam et al., 2009] в области Южно-Атлантической магнитной аномалии (40° S; 50° W).
ЗАКЛЮЧЕНИЕ
В работе на примере задачи расчета ЖГО и построения карт ЖГО был протестирован метод трассировки заряженных частиц Бунемана–Бориса по схеме частица в ячейке в составе разрабатываемого инструментария для моделирования “жизни” заряженной частицы в околоземном пространстве. Было получено хорошее совпадение с аналитическими значениями для дипольного поля, точность при этом составил 3 МВ, что соответствует 1.5% в экваториальной области. Были воспроизведены основные особенности ЖГО в “реальном” поле, заданном моделью IGRF. В частности, восстановлена характерная картина полутени в нескольких географических точках, которая хорошо совпала по структуре с ранее полученными результатами. Данный метод трассировки по сравнению широко используемым методом Рунге–Кутта показал высокую точность, устойчивость и может использоваться в таких методиках как (https://www.geomagsphere.org/index.php; http://cosmos.hwr.arizona.edu/Util/rigidity.php ) особенно при необходимости массовых расчетов.
Кроме того, построены карты ЖГО и их вековые вариации в условиях спокойной магнитосферы в зависимости от разных направлений, в которых также были воспроизведены известные эффекты, что является дополнительной валидацией метода.
ФИНАНСИРОВАНИЕ
Исследование выполнено за счет гранта Российского научного фонда № 19-72-10161, (https://rscf.ru/project/19-72-10161/).
Об авторах
П. А. Кручинин
Московский инженерно-физический институт
Автор, ответственный за переписку.
Email: kruchinin_01@inbox.ru
Россия, Москва
В. В. Малахов
Московский инженерно-физический институт
Email: vvmalakhov@mephi.ru
Россия, Москва
В. С. Голубков
Московский инженерно-физический институт
Email: vlad10433@mail.ru
Россия, Москва
А. Г. Майоров
Московский инженерно-физический институт
Email: agmayorov@mephi.ru
Россия, Москва
Список литературы
- Голубков В.С., Майоров А.Г. Пакет программ для численных расчетов траектории частиц в магнитосфере Земли и его применение для обработки данных эксперимента PAMELA // Изв. РАН. Сер. физ. Т. 85. № 4. С. 512–514. 2021.
- Данилова О.А., Птицына Н.Г., Тясто М.И., Сдобнов В.E. Изменения жесткостей обрезания космических лучей во время бури 8–11 марта 2012 г. в период CAWSES-II // Солнечно-земная физика. Т. 9. № 2. 2023. https://doi.org/10.12737/szf-92202310
- Птицына Н.Г., Данилова О.А., Тясто М.И., Сдобнов В.Е. Динамика жесткости геомагнитного обрезания космических лучей и параметров магнитосферы во время разных фаз бури 20 ноября 2003 г. // Геомагнетизм и аэрономия. Т. 61. № 2. С. 160–171. 2021. https://doi.org/10.31857/S0016794021010120
- Тясто М.И., Данилова О.А., Птицына Н.Г., Сдобнов В.Е. Вариации жесткостей обрезания космических лучей во время сильной геомагнитной бури в ноябре 2004 г. // Солнечно-земная физика. Т. 1. № 2. C. 97–105. 2015. https://doi.org/10.12737/7890
- Adriani O., Barbarino G.C., Bazilevskaya G.A. et al. PAMELA’smeasurements of geomagnetic cutoff variations during the 14 December 2006 storm // Space Weather. V. 14. P. 210–220. 2016. https://doi.org/10.1002/2016SW001364
- Alken P., Thébault E., Beggan C.D. et al. International Geomagnetic Reference Field – The Thirteen Generation // Earth, Planets and Space. V. 73:49. 2021. https://doi.org/10.1186/s40623-020-01288-x
- Boris J.P. The acceleration calculation from a scalar potential // Technical report MATT-152, Princeton: Princeton Univ. 1970.
- Boris J.P. Relativistic plasma simulation-optimization of a hybrid code // Proc. 4th Conf. on Numerical Simulation of Plasmas, Washington. P. 3. 1971.
- Gelvam A., Hartmann, Pacca I.G. Time evolution of the South Atlantic Magnetic Anomaly // Anais da Academia Brasileira de Ciências. V.81. №2. P. 243–255. 2009.
- Koldobskiy S.A., Bindi V., Corti C., Kovaltsov G.A., Usoskin I.G. Validation of the Neutron Monitor Yield Function Using Data From AMS-02 Experiment, 2011–2017 // J. Geophys. Res. Space Physics. 124. P. 2367 – 2379. 2019. https://doi.org/10.1029/2018JA026340
- Kress B.T., Hudson M.K., Selesnick R.S., Mertens C.J., Engel M. Modeling geomagnetic cutoffs for space weather applications // J. Geophys. Res. Space Physics. V. 120. 5694–5702. https://doi.org/10.1002/2014JA020899
- Mao H. and Wirz R.E. // Comparison of Charged Particle Tracking Methods for Non-Uniform Magnetic Fields // 42nd AIAA Plasmadynamics and Lasers Conference. 2011.
- Mishev A.L., Koldobskiy S.A., Kovaltsov G.A., Gil A., Usoskin I.G. Updated Neutron-Monitor Yield Function: Bridging Between In Situ and Ground-Based Cosmic Ray Measurements // J. Geophys. Res. Space Physics. 125. 2020. https://doi.org/10.1029/2019JA027433
- Mishev A., Poluianov S. About the Altitude Profile of the Atmospheric Cut-Off of Cosmic Rays: New Revised Assessment // Solar Physics. V. 296:129. 2021. https://doi.org/10.1007/s11207-021-01875-5
- Qin R., Zhang S., Xiao J., Liu J., Sun Y. Tang W. // Why is Boris algorithm so good? // Physics of Plasmas V. 20.8. P. 084503. 2013. https://doi.org/10.1063/1.4818428
- Sarkar R., Roy A. Monte Carlo simulation of CRAND protons trapped at low Earth orbits // Adv. Space Res. V. 69. P. 197–208. 2022. https://doi.org/10.1016/j.asr.2021.10.006
- Selesnick R.S., Looper M.D., Mewaldt R.A. A theoretical model of the inner proton radiation belt // Space Weather. V. 5. 2007. https://doi.org/10.1029/2006SW000275.
- Smart D.F., Shea M.A. Geomagnetic cutoffs: A review for space dosimetry applications // Adv. Space Res. V. 14. № 10. P. 787–796. 1994.
- Smart D.F., Shea M.A. A review of geomagnetic cutoff rigidities for earth-orbiting spacecraft// Adv. Space Res. V. 36. P. 2012–2020. 2005.
- Smart D.F., Shea M.A. Fifty years of progress in geomagnetic cutoff rigidity determinations // Adv. Space Res. V. 44. P. 1107–1123. 2009.
- Smart D.F., Shea M.A. Vertical Geomagnetic Cutoff Rigidities for Epoch 2015 // Proceedings of science. 2019. // 36th International Cosmic Ray Conference. ICRC2019. July 24th August 1st. 2019.
- Stormer C. The Polar Aurora. Clarendon Press Oxford. 1955
Дополнительные файлы









