Астрономический журнал, 2021, T. 98, № 1, стр. 66-74

1. ВВЕДЕНИЕ

Транзитная спектроскопия является мощным методом исследования экзопланет, который дает информацию о наличии и составе атмосферы. Для горячих экзопланет спектральные наблюдения также дают информацию об экзосфере и параметрах звездного ветра, вступающего с ней во взаимодействие. Так, например, измерения дуплета натрия позволили обнаружить первую атмосферу у планеты HD 209458b [1], а наблюдения линии Lyα открыли существование обширной атомарной экзосферы вокруг этого горячего юпитера [2]. Аналогичная экзосфера по наблюдениям в Lyα была обнаружена у горячих юпитеров HD 189733b [3, 4], 55 Cancri b [5] и теплых нептунов GJ-436b [6] и GJ-3470b [7].

В 2011 г. для планеты HD189733b было зарегистрировано значительное поглощение в линии водорода Hα 656.3 нм на уровне 1–2% [8, 9]. Особенностью данной линии, наряду с другими нерезонансными переходами, является ее слабая подверженность влиянию межзвездной среды, геокорональной эмиссии, а также слабое поглощение атмосферой Земли. Это делает возможным наблюдения наземными телескопами, что значительно увеличивает доступность исследования экзопланет. Наблюдения в спектре Hα были сделаны также для экзопланет KELT-9b, MASCARA-2b/KELT-20b и WASP-52b [10–13]. Отметим, что другой нерезонансной линией, представляющей интерес для изучения атмосфер и экзосфер газовых планет, является линия 1083 нм метастабильного гелия [14]. В дополнение к Lyα наблюдения в Hα дают независимые данные о параметрах верхних атмосфер и экзосфер и позволяют независимо верифицировать численные модели, разработанные для расчета планетарного течения горячих экзопланет.

В работах [15, 16] было впервые реализовано моделирование атмосферы планеты HD189733b, где было учтено большое количество различных процессов, влияющих на популяцию атомов водорода, возбужденных на второй энергетический уровень, что в свою очередь приводит к поглощению в линии Hα. Основным процессом возбуждения является поглощение излучения Lyα, которое является интенсивным практически для всех звезд, а также столкновения с энергичными электронами. Наибольшую сложность представляет расчет переноса Lyα фотонов сквозь атмосферу планеты. В работах [15, 16] использовался одномерный расчет методом Монте-Карло, а распределения плотности и температуры атмосферы по радиусу моделировались при помощи одномерной гидростатической модели. В качестве результата было получено хорошее соответствие спектра поглощения с наблюдаемыми данными. Также из работ [15, 16] следует, что радиационное возбуждение второго энергетического уровня атома водорода, по сравнению со столкновительным, может происходить гораздо глубже в атмосфере, где концентрация невозбужденных атомов водорода высока. Следует отметить, что в настоящей работе рассматриваются сразу два горячих юпитера, HD 189733b и HD 209458b. Результаты, касающиеся планеты HD 189733b, служат как для сравнения, так и для верификации используемого подхода. Сравнение с предыдущими исследованиями показывает, что в статье [15] вообще не учитывается процесс радиационного возбуждения на уровень 2p, при этом получены результаты, которые могут объяснить наблюдаемые данные. Между тем в работе [16], помимо приведенных в таблице процессов, учитывается процесс радиационного возбуждения/деактивации атомов водорода Lyα фотонами, и показано, что он также способен объяснить наблюдаемое поглощение. При сравнении работ [15] и [16] роль процесса резонансного рассеяния Lyα фотонов при поглощении в линии Hα не является прозрачной, хотя из текста следует, что учет этого процесса безусловно важен. Таким образом, численное моделирование данной проблемы требует дальнейшего развития.

В настоящей статье приводятся результаты моделирования распределения атомов водорода, возбужденных на второй уровень. В качестве первого приближения рассматривается только радиационное возбуждение. Перенос Lyα фотонов моделируется методом Монте-Карло. Рассмотрены планеты HD189733b и HD209458b, для которых рассчитано поглощение в линии Hα и сделано сравнение с наблюдениями. Параметры атмосферы, а именно, распределения объемной плотности и температуры атомов водорода взяты из результатов моделирования трехмерным многожидкостным газодинамическим кодом. Данный код был разработан нами ранее и применялся для расчета верхней истекающей атмосферы горячих юпитеров HD209458b [17, 18], WASP-12b [19] и теплого нептуна GJ-436b [20]. Разработка кода, способного быстро и адекватно моделировать распределение возбужденного водорода, представляет большое значение не только для интерпретации наблюдательных данных, но и для учета дополнительного канала нагрева верхней атмосферы горячих экзопланет. Как было показано в работе [21] на примере горячего Юпитера KELT-9b, фотоионизация водорода с возбужденного уровня значительно расширяет доступную радиационную энергию, поскольку энергия излучения звезд в интервале <365 нм может быть на много порядков больше, чем в интервале <91.2 нм.

2. ИСХОДНЫЕ ПАРАМЕТРЫ

В качестве исходных параметров планетарных атмосфер были взяты результаты трехмерных расчетов, полученные с использованием модели [17–20]. В модели решаются газодинамические уравнения непрерывности, импульса и энергии для отдельных атомарных, молекулярных и ионных компонент плазмы. Основными процессами взаимопревращения частиц являются фото-ионизация, электронный удар, рекомбинация и перезарядка. Фотоионизация также приводит к нагреву вещества через порождаемые фотоэлектроны. Для каждой конкретной звезды используется наиболее подходящий XUV спектр в диапазоне 10–912 A. Прохождение и ослабление ионизующей радиации вычисляются в каждом спектральном интервале в соответствии с зависящим от длины волны сечением поглощения. Помимо нагрева, также рассчитывается охлаждение за счет возбуждения атомарного водорода и инфракрасного излучения молекулы H3+. Обмен импульсом и энергией между различными компонентами происходит за счет ион-атомных и Кулоновских столкновений. Верхняя атмосфера планет берется состоящей из компонент H, H+, H2, H2+, H3+, между которыми рассчитывается полный набор химических реакций. Код реализован в обобщенном виде, позволяющем включение любых других независимых компонентов, в частности наиболее распространенных – C, O, Mg, Si. Для самосогласованного 3D моделирования на масштабе всей звездной системы те же алгоритмы используются для вычисления динамики звездного ветра. Помимо планетарной атмосферы, другой граничной поверхностью является корона звезды, на которой задаются параметры корональной плазмы. При необходимости также рассчитывается сила радиационного давления на частицы со стороны звездной радиации. Модель выполняет вычисления в неинерциальной сферической системе координат, привязанной к центру планеты и вращающейся в фазе с орбитальным движением, т.е. в приливно захваченном режиме. При этом сама планета может вращаться с любым периодом. В этой системе отсчета вводятся соответствующие неинерциальные силы – обобщенный гравитационный потенциал и сила Кориолиса. Радиальный шаг сетки сильно неоднородный с наименьшим значением вблизи поверхности планеты, где необходимо разрешать сильно стратифицированную верхнюю атмосферу. Код полностью распараллелен для вычислений на кластерах.

Распределения объемной плотности атомов водорода и температуры вдоль радиуса атмосферы планеты показаны на рис. 1. По горизонтальной оси взяты единицы измерения r относительно радиуса планет. Основным фактором, определяющим верхнюю атмосферу горячих юпитеров, является ионизующее излучение, а также компонентный состав атмосферы. Для анализа влияния различных условий на поглощение в линии Hα рассмотрено несколько вариантов модельной верхней атмосферы. Трехмерные расчеты для планеты HD189733b были проведены для двух различных значений параметра XUV, который значительно влияет на ионизацию атмосферы и скорость планетарного течения. Здесь и далее параметр XUV физически соответствует потоку рентгеновского и далекого ультрафиолетового излучения звезды λ < 91.2 нм, измеряемого в единицах [эрг с^–1 см^–2] на расстоянии 1 а.е. Для HD209458b, помимо вариации параметра XUV, было также рассмотрено два различных варианта изначального водородного состава верхней атмосферы: молекулярная и атомная, тогда как для HD189733b атмосфера предполагается изначально атомарной. На рисунках молекулярный вариант отмечен значком (m). Подробный анализ параметров верхней атмосферы HD209458b, получаемый в численном моделировании, и сравнение с другими моделями выполнены в наших предыдущих работах, см., например, [17–20, 22].

Другим важным параметром задачи являются форма линии Lyα и ее интегральная интенсивность. Восстановленные из наблюдений профили Lyα излучения родительских звезд для планет HD189733b и HD209458b взяты из работы [23] и изображены на рис. 2.

3. ОПИСАНИЕ МОДЕЛИ

Чтобы получить распределение плотности возбужденных на второй уровень атомов водорода в атмосфере планеты, была создана численная модель, которая вычисляет перенос Lyα фотонов методом Монте-Карло. Задача решается в приближении изотропного частично когерентного рассеяния [24]. На планету радиусом R₀ с атмосферой, ограниченной радиусом R_a> R₀, падает излучение родительской звезды, находящейся на расстоянии L от планеты. Задача решается в системе отчета приливно захваченной планеты, обращенной к звезде одной стороной. Излучение звезды представляется как сферический фронт, состоящий из задаваемого количества N модельных фотонов, равномерно распределенных по поверхности фронта. Каждому модельному фотону поставлено в соответствие несколько характеристик: координаты в пространстве (x, y, z), направление распространения, частота фотона, выраженная как доплеровская скорость V атомов, излучивших данный модельный фотон, а также пройденный путь, что является внутренними часами фотона и определяет время существования фотона в расчетной области. Начальные координаты и направление распространения определяются в соответствии с положением фотона на фронте, частота (или доплеровская скорость) случайным образом определяется из профиля распределения интенсивности линии излучения звезды [23]. Распределение по доплеровским скоростям задается в виде массива точек в качестве начального параметра. Тогда суммарная интенсивность потока фотонов I вычисляется как интеграл от распределения по скоростям. Количество w реальных фотонов в одном модельном определяется из следующего выражения для интенсивности:

где $S = 4{{R}_{{\text{a}}}}(L - {{R}_{{\text{a}}}})\arcsin \left( {\frac{{{{R}_{{\text{a}}}}}}{{L - {{R}_{{\text{a}}}}}}} \right) \approx 4R_{{\text{a}}}^{2}$ – площадь фронта на расстоянии R_a от планеты, Δt – интервал времени усреднения.

Рассматривается процесс радиационного возбуждения и переизлучения Lyα фотонов в атмосфере планеты. Атмосфера характеризуется сферически симметричными распределениями объемной плотности нейтрального водорода n(r) и температуры T(r), которые задаются в виде массивов в качестве фиксированных исходных параметров. При этом делается допущение, что излучение звезды не влияет на изначальные распределения.

Для трассировки фотонов применяется метод Монте-Карло (см., например, [25]), в рамках которого движение модельного фотона вычисляется следующим образом. Из выражения $\tau = \ln \left( {1 - \xi } \right)$ определяется оптический путь заданного фотона ${{\tau }_{i}}$ на шаге i с использованием равномерно распределенной случайной величины $\xi $. При этом, с другой стороны, $\tau = \int {\alpha (l){\text{d}}l} $, где $\alpha (l)$ – коэффициент ослабления вдоль траектории $l$. Случайная величина $\xi $ так же, как и все остальные случайные величины данной модели, является псевдослучайным числом, генерируемым с использованием вихря Марсенна с периодом ${{2}^{{19937}}} - 1$ из стандартной библиотеки языка C⁺⁺11 (стандарт ISO/IEC 14882:2011). По преодолении фотоном определенного таким образом оптического пути происходит событие поглощения фотона на атоме водорода с дальнейшим переизлучением, при этом фотон поглощается атомом водорода, который в свою очередь переходит в состояние 2p, затем атом живет время ${{t}_{0}} \approx 1.6$ нс (время жизни уровня) и спонтанно излучает фотон. Частота излученного фотона тесно коррелирует с частотой поглощенного фотона и определяется исходя из соображений сохранения энергии: в системе отсчета атома поглощаемый и излучаемый фотоны имеют одинаковые частоты, но в “лабораторной” системе отсчета фотон все же имеет другую частоту, которая определяется углом рассеяния, а также локальной температурой среды (см. подробнее [24]).

Механизм гибели фотонов определяется их выходом из расчетной области, превышением времени жизни наперед заданного времени t_max, а также физическими процессами. В качестве основного физического процесса выбрано явление фотоионизации со второго энергетического уровня атомов водорода за счет полного спектра излучения звезды. Вероятность такой фотоионизации была оценена в соответствии с [16] как 21.5 с^–1 для HD 189733b и 9.4 с^–1 для HD 209458b, что значительно превосходит вероятность гибели Lyα квантов за счет других механизмов. Так, например, вероятность гибели фотонов в результате столкновений с электронами в нашем случае меньше в  ~30 раз (учитывая характерное значение максимальной концентрации электронов, полученной в нашей газодинамической модели ~10⁸ см^–3), а вероятность двухфотонного перехода 2s–1s на 3–4 порядка меньше. Численно гибель фотонов реализуется ограничением количества событий, что соответствует накоплению вероятности фотоионизации с возбужденного состояния. Пройденный фотоном путь ${{l}_{i}}$ между “событиями” в нашем случае вычисляется из выражения (2).

Здесь ${{\sigma }_{{L\alpha }}}$ – сечение поглощения фотона в среде, которое вычисляется из выражения (3), $l$– пройденный фотоном путь, отсчитываемый от предыдущего события.

Сечение поглощения вычислялось как свертка профиля Фойгта с аппроксимацией, полученной в работе [26]:

Здесь ${{\sigma }_{{L\alpha D}}} = 5.9 \times {{10}^{{ - 14}}}$ см² и ${{\sigma }_{{L\alpha N}}} = 2.6 \times {{10}^{{ - 19}}}$ см² есть сечение поглощения за счет доплеровского уширения и за счет естественной ширины линии атома, соответственно, ${{x}^{2}} = {{m}_{p}}{{V}^{2}}{\text{/}}2kT$, ${{m}_{p}}$ – масса протона, $k$ – постоянная Больцмана, а явный вид функции q приведен в [26].

Атмосфера дискретизировалась с шагом dr = = 0.01 так, что уменьшение шага сетки в два раза не влияло значительным образом на результат вычислений.

Для каждого события регистрируется момент времени, а также координаты. Таким образом, строится усредненная по интервалу времени Δt динамическая картина распределения возбужденных в состояние 2p атомов водорода в атмосфере планеты, образовавшихся от одного фронта фотонов $n_{{2{\text{p}}}}^{1}$. Поскольку задача симметрична относительно оси “планета-звезда” (ось z), то расчет траекторий фотонов проводился в цилиндрических координатах с шагом сетки $\Delta r$ и $\Delta z$. Объемная плотность возбужденных атомов находится из количества событий ∆N, произошедших в объеме пространства $2\pi r\Delta r\Delta z$ за время Δt:

Полученная картина возбужденных атомов является своего рода импульсной характеристикой системы – реакцией воздействия на группу фотонов от бесконечно тонкого фронта. При этом в приближении малого влияния фотонов на атмосферу планеты система предполагается линейной. Для получения реакции на непрерывное излучение звезды необходимо взять свертку импульсной характеристики с функцией Хевисайда η:

Искомое решение получается, когда выражение выходит через некоторое время на стационарное состояние.

Поглощение в линии Hα вычислялось из распределения ${{n}_{{2{\text{p}}}}}(r,z)$ аналогично [26] из соотношения (5):

Здесь интенсивность I имеет индексы “transit” и “out”, что соответствует интенсивности линии для случаев во время транзита и вне его. При этом ${{\tau }_{{H\alpha }}}(V) = \int_{ - R}^R {{{n}_{{2{\text{p}}}}}(r,z){{\sigma }_{{H\alpha }}}(V,T)dz} $, где ${{\sigma }_{{H\alpha }}}(V,T)$ вычислялось аналогично выражению (3) с соответствующими параметрами перехода: ${{\sigma }_{{H\alpha D}}} = 4.9 \times $ $ \times \;{{10}^{{ - 13}}}$ см² и ${{\sigma }_{{H\alpha N}}} = 5.14 \times {{10}^{{ - 18}}}$ см², которые, в свою очередь, рассчитывались для усредненных значений силы осциллятора и времени жизни перехода Hα [27].

4. РЕЗУЛЬТАТЫ

На рис. 3 приведены распределения объемной плотности возбужденных на второй уровень атомов водорода (n_2p) для обеих планет, полученные с использованием описанной выше Монте-Карло модели. На рис. 4 приведены профили концентрации вдоль оси симметрии z при r = 1 и поперек оси при z = 1. Как следует из выражения (2), оптический путь линейно зависит от концентрации невозбужденных атомов водорода, что и наблюдается в наших расчетах – профили n_2p для всех случаев прямо коррелируют с профилями невозбужденных атомов.

На рис. 5 представлено распределение поглощения в центре линии Hα вдоль прицельного параметра b – расстояния от центра планеты перпендикулярно лучу зрения. Для планеты HD189733b поглощение происходит в слое порядка полутора радиусов планеты, что соответствует n_2p в пределах от ~10⁰ до ~10². При этом интенсивность XUV оказывает незначительное влияние на поглощение, так как не сильно влияет на профиль концентрации невозбужденных атомов в данном слое, как и на температуру (см. рис. 1). Это связано с тем, что увеличение ионизующего излучения увеличивает локальную фотоионизацию атомов водорода и нагрев газа, но это компенсируется увеличением скорости течения.

Для планеты HD209458b видна большая разница для случаев изначально молекулярной и атомарной атмосфер, что, очевидно, обусловлено различиями в концентрациях невозбужденных атомов. Различие в профилях атомарной и молекулярной атмосферы (рис. 1) связано в свою очередь с тем, что молекулярная атмосфера имеет в два раза меньшую шкалу высот. В атомарном случае поглощение, как и для HD189733b, происходит в слое порядка 1.5 радиусов планеты, что соответствует n_2p в пределах от ~10⁰ до ~10¹ см^–3. В случае изначально молекулярной атмосферы наблюдается значительное уменьшение поглощающего слоя радиусов планеты, при этом n_2p находится в пределах ~10^–1 до ~10⁰ см^–3.

На рис. 6 изображены спектры пропускания (полученные с использованием формулы (5)) в линии Hα и их сравнение с наблюдательными данными. В работах [8, 9] для планеты HD189733b в центре линии наблюдалось поглощение на уровне 1.5–2%, что в 3–4 раза выше поглощения ~0.5%, полученного с использованием данной модели. Как показано в [16], значительный вклад в Hα поглощение также дают Lya-фотоны, которые рождаются в атмосфере планеты под действием внутренних процессов.

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

Таким образом, разработана Монте-Карло модель, вычисляющая перенос Lyα фотонов через верхнюю атмосферу горячих экзопланет. На основе разработанной модели для горячих юпитеров HD189733b и HD209458b получены распределения объемной плотности атомов водорода, возбужденных на второй энергетический уровень. Полученные распределения позволяют вычислить спектр поглощения данных планет в линии Hα. Полученные результаты свидетельствуют о значительном вкладе звездного излучения в поглощение на данной линии. Интенсивность XUV незначительно сказывается на поглощении, так как слабо влияет на концентрацию атомов в поглощающем слое.

Для планет HD189733b и HD209458b с атомарной атмосферой расчет показал, что поглощение в центре линии порядка ~0.5% с шириной ~40 км/с. При этом поглощение происходит в слое ~1.5 радиусов планеты, где концентрация возбужденных на второй уровень атомов водорода больше 10⁰ см^–3.

Для случая HD209458b с молекулярной атмосферой поглощение оказалось менее 0.1%. Отметим, что для данной планеты интерпретация имеющихся наблюдательных данных в линии Hα не очевидна [8] и скорее свидетельствует о низком уровне поглощения (<1%).

По результатам работы можно сделать заключение, что транзитное поглощение горячих юпитеров в линии Hα значительно зависит от Lyα излучения звезды, но не может быть полностью им объяснено. Поглощение образуется в пределах полости Роша планет, где скорость истечения еще невелика. Поэтому ширина линии поглощения хорошо соответствует температурному максимуму термосферы (1–1.5) × 10⁴ К. Таким образом, измерения в Hα потенциально позволяют зарегистрировать раздувание верхней атмосферы горячих юпитеров за счет поглощения ионизующей радиации родительских звезд. Основной вывод заключается в том, что значительный вклад в поглощение могут вносить внутриатмосферные процессы, приводящие к рождению Lyα фотонов, в первую очередь электронное возбуждение в области максимума температуры. Это предполагается изучить в следующей работе.