Астрономический журнал, 2021, T. 98, № 12, стр. 1043-1056

2.1. Моделирование транзитов

Транзитная кривая блеска зависит от большого числа величин, определяющих свойства звезды, планеты и ее орбиты. В качестве выходных параметров нас в первую очередь интересует размер планеты, кроме того мы определяем наклонение орбиты и, с некоторой оговоркой (см. ниже), эксцентриситет. Для получения значений этих величин прочие (входные) параметры должны быть заданы, а если это невозможно, то так или иначе определены. Мы осуществляем это перебором некоторых коэффициентов (например, задающих потемнение звездного диска к краю) в широком диапазоне возможных значений и минимизацией ${{\chi }^{2}}$ для полученных кривых блеска.

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

Координаты центра планеты в картинной плоскости определялись нами как функция от времени и кеплеровых элементов орбиты. Зная период обращения планеты по орбите вокруг материнской звезды, можно установить линейную зависимость средней аномалии от времени $M(t)$. Для этого используется уравнение Кеплера:

где $e$ – эксцентриситет, $M$ – средняя аномалия, а $E$ – эксцентрическая. Через среднюю аномалию определяется эксцентрическая. Так как уравнение Кеплера является трансцендентным, то его решение проводилось при помощи последовательных приближений:

Далее, зная эксцентрическую аномалию и ориентацию орбиты в пространстве, можно найти координаты центра планеты в картинной плоскости для каждого момента времени. В программе это выполняется за несколько шагов. Сначала рассчитываются координаты для точек эллипса, лежащего на картинной плоскости, с перицентром, совпадающим с точкой восходящего узла, и линией апсид, совпадающей с осью $x$. В таком случае координаты центра диска планеты на картинной плоскости в разные моменты времени (что неявно задается через значение эксцентрической аномалии $E$) рассчитываются как:

где $a$ – большая полуось, ${{x}_{{02}}}$ и ${{y}_{{02}}}$ – координаты проекции центра планеты на картинную плоскость в первой итерации.

Затем производится поворот эллипса на угол, равный аргументу перицентра, осуществляемый при помощи двумерной матрицы поворота. Соответственно, координаты центра диска планеты получаются равными:

Здесь $\omega $ – аргумент перицентра, ${{x}_{{01}}}$ и ${{y}_{{01}}}$ – координаты проекции центра планеты на картинную плоскость после второй итерации.

После этого производится поворот по оси $x$, отвечающий за учет наклона орбиты, что дает нам финальные координаты центра диска планеты на картинной плоскости:

где $i$ – наклон орбиты (угол между плоскостью орбиты и картинной плоскостью), ${{x}_{0}}$ и ${{y}_{0}}$ – итоговые координаты проекции центра планеты на картинную плоскость.

После определения геометрии задачи необходимо рассчитать яркость различных участков диска звезды. В нашей программе используется центрально-симметричная модель потемнения к краю, описываемая либо линейным, либо квадратичным законами [10]. Линейный закон использовался в основном для иллюстративных целей, в частности, для иллюстрации вырождения по параметрам (см. ниже раздел 4).

Квадратичный закон записывается в виде:

где $\mu = cos\gamma $, $\gamma $ – угол между лучом зрения и нормалью к участку излучающей поверхности. Так как линейный закон представляет собой частный случай квадратичного, в программе он задавался путем фиксирования коэффициента ${{u}_{2}} = 0$.

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

2.2. Расчет эксцентриситета и аргумента перицентра

Так как орбиты большинства экзопланет отличны от круговых, характеристики их вытянутости и ориентации относительно наблюдателя также необходимо учитывать при обработке кривых блеска. Большую полуось орбиты можно установить с относительно хорошей точностью, зная период обращения и массу системы. Период определяется измерениями промежутка между транзитами, масса звезды может быть получена, основываясь на ее спектральных характеристиках. Соотношение масс звезды и планеты может быть измерено на основе данных о вариации лучевой скорости звезды, если они доступны, в противном случае предположение о том, что масса звезды намного превосходит массу планеты, позволяет получить достаточно точные оценки.

Однако если орбита отлична от круговой, то расстояние, на котором будет находиться планета в момент транзита, будет больше или меньше большой полуоси и, соответственно, от этого будут зависеть скорость планеты и продолжительность транзита. Кроме того, расстояние между планетой и звездой, а также скорость движения планеты по орбите, будут меняться в течение самого транзита, что приведет к тому, что кривая блеска будет иметь некоторую асимметрию относительно середины транзита. Однако при небольших значениях эксцентриситета ($e \lesssim 0.25$) вследствие того, что продолжительность транзита намного меньше периода обращения, расстояние до звезды в течение транзита и скорость планеты меняются незначительно. Поэтому мы пренебрегали асимметрией транзита, связанной с эксцентричностью орбиты. Формально можно сказать, что это соответствует предположению о том, что транзит происходит либо в точке апоцентра, либо в точке перицентра. Таким образом, аргумент перицентра не рассчитывался в нашей программе. Единственная величина, которую можно получить при точности тех наблюдений, что имелись в нашем распоряжении, это среднее расстояние от звезды до планеты в момент транзита, характеризуемое величиной $esin\omega $.

Это позволило сократить вычисления, которые в ином случае потребовались бы для определения средней точки транзита, совпадающей в данном случае с серединой транзита.

2.3. Связь радиуса планеты и наклона орбиты

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

В момент касания планетой края звездного диска для декартовых координат центра планеты можно записать следующую систему уравнений (плоскость $xy$ совпадает с картинной, ось $z$ направлена на наблюдателя).

Здесь $p$ – фокальный параметр, угол $\phi $ равен разности истинных аномалий $\Delta \nu $. Решив эту систему, получаем функцию $r(i)$,

2.4. Определение диапазона параметров, в пределах которого производится аппроксимация

Чтобы найти максимальное и минимальное значения наклона орбиты, между которыми будет осуществляться поиск наиболее оптимального решения, вычисляется медианное значение блеска для точек кривой блеска вблизи середины транзита. С учетом величины заданных погрешностей определяется диапазон возможных значений блеска в средней точке. После определяются значения наклона орбиты (и соответствующие им значения радиуса), при которых блеск в середине транзита будет иметь соответствующие величины. Для этого проводится моделирование блеска в одной только средней точке. Далее получившийся диапазон значений наклона орбиты делится на шесть частей. Для каждого из получившихся значений описанная выше программа моделирования строит модель кривой и рассчитывает величину ${{\chi }^{2}}$ как сумму квадратов разностей O–C (наблюдаемые значения минус рассчитанные), деленных на погрешности измерения для каждой отдельной точки. Затем точки справа и слева от той, которой соответствует наименьший квадрат отклонений, принимаются в качестве новых границ аппроксимации. При $N$-кратном повторении данной процедуры диапазон возможных значений наклона орбиты сужается в ${{3}^{N}}$ раз. В нашей работе хватало проведения 6–7 подобных итераций, прежде чем разница между двумя границами становилась много меньше неточностей, связанных с погрешностями фотометрии транзита.

Мы не стали использовать стандартные методы нахождения минимума числовой функции (например, метод золотого сечения или парабол) из-за особенностей функции ${{\chi }^{2}}(i)$, в частности, наличия у нее в некоторых случаях локальных минимумов.

После определения минимального значения ${{\chi }^{2}}$ вычислялись величины погрешностей по наклону орбиты из соотношения:

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

2.5. Определение прочих параметров модели

Ключевым параметром, определяемым по транзитным кривым в рамках нашей модели, является радиус планеты. Однако его значение зависит также от коэффициентов потемнения к краю и эксцентриситета. Нами проводился перебор значений этих параметров в достаточно широких реалистичных диапазонах. Кроме того, использовались два значения аргумента перицентра: $sin\omega = 1$ и $sin\omega = - 1$.

Для каждого фиксированного набора значений проводилось фитирование транзитной кривой блеска по радиусу и наклону орбиты и, соответственно, рассчитывалась величина суммы квадратов отклонений. В качестве конечного результата принимался тот набор параметров, при котором фитирование по радиусу и наклону орбиты давало наименьшую сумму квадратов отклонений.

Также нами проводилось фитирование с использованием фиксированных коэффициентов потемнения к краю. Для их получения мы использовали апплет, представленный на сайте университета Огайо²². Он позволяет получить коэффициенты потемнения к краю для вводимых пользователем параметров звезды (таких как эффективная температура, металличность и ускорение свободного падения на поверхности звезды). Работа апплета основана на результатах статьи [11]. В ней с использованием цепочек Маркова проводилась интерполяция коэффициентов из работы [10], в которой, в свою очередь, исходные коэффициенты получались путем моделирования звездных атмосфер. К сожалению, теоретические модели пока недостаточно точны. В данной работе мы приводим результаты как для фиксированных параметров потемнения к краю, так и для случая перебора (см. ниже разделы 4–5).

2.6. Teстирование программы на синтетических кривых и оценка производительности

Для тестирования программы была сгенерирована искусственная кривая блеска транзита экзопланеты. За основу мы взяли параметры экзопланеты WASP-97 b [12]: радиус звезды был задан равным $R = 1.12{{R}_{ \odot }}$ [13], коэффициенты потемнения к краю были взяты из [11]. Интервал по времени между точками был выбран равным одной минуте, таким образом вся кривая блеска с учетом участков до и после прохождения планеты по диску включала в себя 198 точек. К смоделированной кривой была добавлена шумовая дорожка шириной 0.003 от блеска звезды вне транзита. После этого для данной кривой блеска была проведена аппроксимация с помощью программы, представленной в данной статье. Полученные результаты приведены в табл. 1, а сама кривая блеска – на рис. 1. Приведенная в таблице (и далее в других таблицах) погрешность для определяемых параметров соответствует одному стандартному отклонению.

Также мы провели исследование того, каким образом на результатах отразится изменение входных параметров, таких как радиус звезды и коэффициенты потемнения к краю. Для этого нами было проведено варьирование радиуса на $0.05{{R}_{ \odot }}$ и коэффициентов потемнения к краю, соответствующих изменению температуры звезды на 200 K, что приблизительно соответствует погрешностям, приведенным в работе [13]. Полученные результаты представлены в табл. 2. Кривые блеска, соответствующие вариациям параметров, практически не отличимы от основной кривой, поэтому на рис. 1 мы их не приводим.

Сравнение данных в табл. 1 и 2 показывает, что программа хорошо определяет такие параметры планеты, как радиус и наклонение орбиты даже при недостаточно точно известном радиусе и коэффициентах потемнения к краю (хотя при существенных одновременных неопределенностях в радиусе и температуре точное значение, заложенное при синтезе кривой, перестает попадать в доверительные интервалы, соответствующие одному стандартному отклонению).

Общее время выполнения программы для данной кривой на использованном нами процессоре Intel Core i7 1.6 ГГц составило 7531 мс.

4.1. Результаты обработки транзитов

Фотометрические данные обрабатывались программой, алгоритм которой описан выше в разделе 2. В приведенных ниже таблицах приведены параметры, для которых сумма квадратов отклонений оказалась минимальной. В столбце “Наши расчеты” приведены результаты расчетов, в которых варьировались как коэффициенты потемнения к краю, так и параметр $esin\omega $. В столбце “Наши расчеты, фиксированные коэффициенты” приведены результаты расчетов, в которых в программе перебиралось только значение эксцентриситета, в то время как коэффициенты потемнения к краю принимались равными значениям, рассчитанным для данной звезды с помощью апплета Джейсона Истмана, упоминавшегося в подразделе 2.5 (см. [11]). Параметры TrES-3b, которые мы не варьировали (такие как радиус звезды, большая полуось орбиты, период обращения), были взяты из работы [18]. Для нахождения моментов начала и конца транзита, необходимых для расчета зависимости $r(i)$, мы использовали программу фитирования, представленную на сайте ETD [16] и использующую аппроксимацию кривой блеска аналитическими функциями из работы [19].

Мы использовали три кривые блеска с сайта ETD. Первая из них была получена Марком Бреттоном в 2015 г. Середина транзита соответствует HJD = 2457 208.42831. Наблюдения проводились на 82-см телескопе Observatoire des Baronnies Provencales. Вторая была получена Пере Гуэрра в 2018 г. на 40-см телескопе Observatori Astronómic Albanyа. Середина транзита: HJD = = 2458 258.59999. Третья кривая была получена Соццетти в 2007 г. и была использована в работе [18]. Результаты обработки этих трех кривых блеска представлены в табл. 4–6 и на рис. 2–4.

Результаты обработки транзита TrES-3b, наблюдения которого проводились в КГО, и полученная для него модель кривой блеска представлены в табл. 7 и на рис. 5. Кроме того, в табл. 7 приведены для сравнения параметры планеты, полученные нами при помощи сетевой версии программы EXOFAST [11]³³, а также результаты работы [18], полученные на основе нескольких кривых блеска в фильтрах $g$, $V$, $B$, $z$, $i$ с использованием метода аппроксимации, описанного в [20] (сравнение с результатами работы [18] уместно для фиксированных коэффициентов потемнения к краю). Как видно из приведенных данных, параметры, полученные с помощью нашей программы, хорошо соотносятся с опубликованными оценками [18] (а также с другими работами, см. каталоги экзопланет), в то время как E-XOFAST дает заметно большее значение для радиуса планеты. Чтобы проверить, не связано ли это с величинами коэффициентов в законе потемнения к краю, мы провели расчеты, используя те же значения, что и в EXOFAST, но расхождение не исчезло. То же самое можно сказать о различии в используемом радиусе звезды. Большой радиус планеты получается в наших расчетах при использовании значения $esin\omega $, соответствующего результатам EXOFAST, однако при этом угол наклона орбиты становится немного (примерно на один градус) меньше.

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

Как указано выше, в результате анализа полученных результатов мы пришли к выводу, что разумно использовать фиксированные значения коэффициентов потемнения к краю, основанные на данных о свойствах звезды. Если коэффициенты потемнения к краю являлись свободными параметрами, то наши расчеты в ряде случаев (см. также ниже раздел 5) приводили к нефизичным оценкам их величин. В связи с этим мы провели дополнительное исследование возможного вырождения параметров потемнения к краю и параметра $esin\omega $. Этот вопрос рассматривается в следующем подразделе 4.2.

4.2. Распределение значений квадрата отклонений по параметрам эксцентриситета и потемнения к краю

При обработке данных наблюдений TrES-3b было обнаружено, что при незначительных изменениях начальных параметров мы могли получать значительно отличавшиеся между собой величины $esin\omega $ и коэффициентов потемнения к краю. При этом значения суммы квадратов отклонений отличались друг от друга незначительно (см. ниже). Это не позволяло нам определенно сказать, какой из полученных наборов величин являлся наиболее оптимальным решением задачи аппроксимации. Тем не менее для данных результатов наблюдалась определенная закономерность: чем сильнее было падение яркости при приближении к краю звезды, тем более вытянутой была орбита планеты в конечной модели. Зависимость коэффициента потемнения к краю при линейном законе (${{u}_{2}} = 0$) от параметра $esin\omega $ представлена на рис. 6.

Также мы провели исследование аналогичной зависимости для квадратичного закона потемнения к краю. Для этого мы брали коэффициенты потемнения к краю из апплета Истмана и др. [11] для случая, когда звезда имеет такую же металличность и ускорение свободного падения на поверхности звезды, как и TrES-3, но другую эффективную температуру. Ее мы в свою очередь варьировали на 1200 K в большую и меньшую сторону от реальной температуры TrES-3, примерно равной 5700 K [12]. Полученные результаты представлены на рис. 7. Так как при прочих равных параметрах большей температуре звезды соответствует более слабый эффект потемнения к краю, можно говорить о том, что здесь также наблюдается аналогичная зависимость между потемнением к краю и расстоянием до звезды в момент транзита. При этом величина ${{\chi }^{2}}$ при увеличении температуры менялась монотонно c 2.4866 при 4500 K до 2.4779 при 6900 K.

Чтобы установить полную картину, мы построили модели для всех возможных значений расстояния от планеты до звезды в момент транзита и линейного закона потемнения к краю по двум транзитам TrES-3b. Первый из них наблюдался нами в КГО, а второй взят из базы ETD (кривая блеска, полученная Марком Бреттоном и описанная выше).

Коэффициент ${{u}_{1}}$ в линейном законе потемнения варьировался от 0 до 1 с шагом 0.02. Параметр $esin\omega $ варьировался в дипазоне от 0 до 0.2 с шагом 0.005. Выбор такого диапазона был сделан нами исходя из работ других исследователей данной планетной системы [21, 12 ]. Для каждой пары значений потемнения к краю и расстояния до звезды программа выполняла аппроксимацию по наклону орбиты и радиусу планеты, как это описано в разделе 2, и определяла наиболее оптимальные значения данных параметров в каждом случае и соответствующее им значение суммы квадратов отклонений ${{\chi }^{2}}$. В итоге мы получили двумерное распределение ${{\chi }^{2}}$ по коэффициенту потемнения к краю и значению $esin\omega $, которое представлено на рис. 8 и 9.

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

Другой факт, обнаруженный нами при построении данного распределения, заключается в вырождении по параметрам потемнения к краю и эксцентриситета. Величина ${{\chi }^{2}}$ по мере движения от круговых орбит с низким потемнением к краю к вытянутым орбитам со значительным потемнением к краю менялась слабо. Поэтому использование лишь метода наименьших квадратов для определения того, в какой из областей вырождения следует искать истинные значения параметров транзита, может оказаться недостаточным.

5.1. HAT-P-19b

HAT-P-19b – это горячий супер-юпитер, вращающийся вокруг звезды Главной последовательности и открытый в 2010 г. в рамках проекта Hungarian Automated Telescope Network [22].

Результаты обработки кривой блеска представлены в табл. 8 и на рис. 10. Величины, полученные нами как для случая с фиксированными коэффициентами потемнения к краю, так и при их вариации в широком диапазоне, показывают, что полученные значения параметров (как для наклона орбиты, так и для радиуса планеты) совпадают в пределах погрешностей как друг с другом, так и с результатами, полученными в работе [22]. Однако коэффициенты потемнения к краю, для которых были получены наилучшие результаты аппроксимации при варьировании данных коэффициентов в широком диапазоне, значительно отличаются от тех, что были получены исходя из теоретических моделей в работе [11].

5.2. KOI-196b

KOI-196b (другое название Kepler-41b) – горячий супер-юпитер, вращающийся вокруг желтого карлика спектрального класса G2. Эта планета была открыта с помощью спутника “Kepler” в 2011 г. Результаты обработки кривой блеска представлены в табл. 9 и на рис. 11. В данном случае полученные нами значения радиуса оказались значительно меньше, чем значения, полученные для данной планетной системы другими исследователями, в частности [12]. В то же время угол наклона орбиты по нашим расчетам получился значительно выше, что объясняет совпадение модели и исходной кривой по глубине. Чем больше угол наклона орбиты, тем меньше прицельный параметр и, соответственно, планета проходит во время транзита через более яркие области ближе к центру звезды. Из-за этого достигается та же глубина транзита при меньшем радиусе планеты.

Разница в результатах, возможно, связана с тем, что кривая блеска транзита данной планеты имеет значительно больший уровень шума по сравнению с остальными кривыми блеска, полученными в КГО. Также сказываются короткие времена наблюдения в моменты до и после транзита, что, возможно, привело к недостаточно точному определению уровня изначальной яркости звезды по данным наших измерений.

5.3. WASP-60b

WASP-60b – это горячий супер-юпитер, вращающийся вокруг звезды Главной последовательности спектрального класса F9. Планета была открыта в 2011 г. в рамках проекта SuperWASP.

Результаты обработки кривой блеска представлены в табл. 10 и на рис. 12. Для данной планеты результаты наших расчетов оказались совпадающими в пределах ошибок с результатами, полученными в статье [23]. В данном случае кривая блеска имела высокий относительно других использованных нами транзитов уровень шума, но в то же время она имела достаточно большое число точек, что позволило получить близкие к другим исследованиям результаты аппроксимации.