Астрономический журнал, 2020, T. 97, № 12, стр. 986-996

Миграция горячих Юпитеров под действием истечения атмосферы

Е. П. Курбатов 1*, Д. В. Бисикало 1**, И. Ф. Шайхисламов 2***

1 Институт астрономии РАН
Москва, Россия

2 Институт лазерной физики СО РАН
Новосибирск, Россия

* E-mail: kurbatov@inasan.ru
** E-mail: bisikalo@inasan.ru
*** E-mail: ildars@ngs.ru

Поступила в редакцию 27.04.2020
После доработки 03.07.2020
Принята к публикации 30.07.2020

Полный текст (PDF)

Аннотация

Как известно, воздействие ионизующего излучения и гравитации звезды на “горячий Юпитер” приводит к истечению его атмосферы. В результате гравитационного воздействия со стороны планеты истекающее вещество приобретает угловой момент, после чего оно накапливается на более высокой круговой орбите вокруг звезды, формируя диск или тор. Обмен угловым моментом между тором и планетой вызывает миграцию планеты к звезде. В данной работе мы рассматриваем эффективность такого механизма миграции на примере системы HD 209458. Оказалось, что за 4.5 × × 109 лет, прошедших после испарения протопланетного диска, планета может мигрировать с орбиты $ \gtrsim {\kern 1pt} 0.67$ а.е. до ее текущей орбиты $0.045$ а.е.

1. ВВЕДЕНИЕ

Принято считать, что основной причиной существования газовых гигантов на орбитах, близких к звезде ($0.1$ а.е. и меньше), является миграция планет в результате гравитационного взаимодействия с протопланетным диском. В зависимости от массы планеты, плотности диска и от физических условий в газе выделяют различные механизмы миграции [1, 2]. Планеты с массой порядка нескольких масс Земли подвержены миграции типа I [3, 4]: гравитационное воздействие планеты возбуждает в диске приливные волны, темп обмена орбитальным угловым моментом определяется интенсивностью гравитационного взаимодействия планеты и приливных волн. Для планеты с массой порядка массы Юпитера в газовом диске с большой плотностью, как правило, предполагают миграцию типа II [5]: угловой момент сил, действующий со стороны планеты, приводит к выметанию газа из окрестности ее орбиты, образуя щель (gap, “гэп”), темп обмена угловым моментом определяется временем заполнения щели, которое зависит от вязкости газа. Возможны также и другие механизмы миграции (см., напр., [6]. Планеты промежуточных масс могут участвовать в смешанном сценарии I + II или в миграции типа III (см. работу [7] и ссылки в ней). Если брать в расчет одновременное присутствие в газовом диске нескольких планет, число возможных сценариев миграции сильно возрастает. В случае отсутствия газового диска обмен орбитальным угловым моментом может происходить с диском планетезималей, через механизм Козаи и др. [2].

Наблюдения показывают, что время жизни газового протопланетного диска составляет ${{10}^{6}}{\kern 1pt} - {\kern 1pt} $ 107 лет [8], и это находится в согласии с теорией [9]. Характерное время миграции планетной орбиты зависит от плотности протопланетного диска и для типичных дисков происходит в более короткой шкале времени, от $ \gtrsim {\kern 1pt} {{10}^{4}}$ до $ \gtrsim {\kern 1pt} {{10}^{6}}$ лет, в зависимости от механизма миграции [1]. Это приведет к тому, что большинство планет должно падать на звезду, либо должно иметь очень низкие орбиты. Для того, чтобы разрешить противоречие с существующими взглядами на процесс образования массивных планет, а также с наблюдениями, предлагаются различные сценарии, призванные замедлить миграцию или даже обратить ее вспять [2, 7, 10]. Эта проблема, однако, еще далека от решения.

Для планеты с массой ${{M}_{{\text{p}}}} \sim {{M}_{{{\text{Jup}}}}}$ обычно полагают, что миграция орбиты происходит по типу II, т.е. время миграции в конечном итоге контролируется вязкостью диска. Если диск имеет поверхностную плотность $\Sigma $ и вязкость $\nu $, это время можно оценить как ${{t}_{{{\text{mig}}}}} \sim {{M}_{{\text{p}}}}{\text{/}}(\Sigma {{\nu }^{2}})$ [2]. Однако, если планета изначально находилась во внешней части диска, где плотность мала, характерное время изменения ее орбиты может превышать время жизни газового диска. В настоящей работе мы предлагаем механизм миграции массивной планеты в случае, когда газового протопланетного диска уже нет, а другие способы изменения орбитального углового момента неэффективны. Идея заключается в следующем.

Звезда и массивная планета образуют двойную систему. Если за время жизни газового протопланетного диска орбита планеты успевает приблизиться к звезде достаточно близко, то включается процесс газодинамического истечения верхней атмосферы под воздействием ионизующего излучения и гравитации звезды. Более того, в первые сотни миллионов лет светимость звезды в коротковолновом спектре может быть на два порядка больше, чем в последующие эпохи [11]. Характерная скорость истечения составляет $ \sim 10$ км/с, что на порядок меньше типичной орбитальной скорости, в результате планетарное вещество накапливается на орбите. В отсутствие звездного ветра это вещество должно формировать кольцо (тор) на планетной орбите. Как показали газодинамические и МГД расчеты [1216], вещество, покинувшее атмосферу, получает дополнительный импульс от ветра и, следовательно, кольцо должно формироваться на орбите выше планетной. Общие потери вещества планетой типа горячий юпитер достаточно велики [17] и образующийся под действием диссипативных процессов диск может быть массивным. Приливное взаимодействие планеты с этим диском должно приводить к уменьшению орбитального углового момента планеты и, как следствие, к миграции ее в направлении звезды.

Численное моделирование процесса формирования планетарного тора находится пока на начальном этапе. В работе [18] в 3D МГД-расчетах обнаружено формирование двух характерных режимов взаимодействия планетарного течения со звездным ветром. В случае сильного ветра планетарное вещество полностью уносится с орбиты планеты, при этом формируется истечение, похожее на хвост кометы. В случае относительно слабого звездного ветра часть планетарного вещества быстро теряет угловой момент и аккрецирует на звезду, в то время как другая часть вытягивается вдоль орбиты планеты. В работе [19] был проведен 3D гидродинамический расчет истечения планетарного вещества с учетом радиационного давления и при условии слабого звездного ветра. Расчет показал формирование тора из планетарного вещества, сосредоточенного в основном за пределами орбиты планеты, при этом радиационное давление частично остановило аккрецию вещества на звезду. В работе [20] расчет для горячего Юпитера HD 209458b, в условиях относительно слабого звездного ветра, показал формирование тора планетарного вещества с толщиной порядка $20$ радиусов планеты ${{R}_{{\text{p}}}}$ и радиальной шириной порядка $40{{R}_{{\text{p}}}}$. В работе [21] в рамках 3D МГД-моделирования была рассчитана относительно долговременная динамика накопления планетарного вещества в торе на протяжении сотен орбитальных периодов. Однако эффекты, которые мы изучаем в настоящей работе, проявляются на временах порядка десятков и сотен миллионов лет, и приведенные выше примеры численного моделирования, хотя и подтверждают предположение о формировании тора (диска), не могут быть использованы непосредственно для оценок миграции планеты.

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

2. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим систему, содержащую звезду солнечного типа массой ${{M}_{{\text{s}}}} = {{M}_{ \odot }} = 2 \times {{10}^{{33}}}$ г и планету массой ${{M}_{{\text{p}}}} = {{M}_{{{\text{Jup}}}}} \approx {{10}^{{ - 3}}}{{M}_{ \odot }} = 2 \times {{10}^{{30}}}$ г и радиусом ${{R}_{{\text{p}}}} = {{R}_{{{\text{Jup}}}}} \approx 7 \times {{10}^{9}}$ см. Планета находится на круговой орбите с большой полуосью $a$ и образует со звездой двойную систему.

Атмосфера планеты разогревается УФ-излучением звезды. Если температура в экзосфере становится порядка или больше ${{10}^{4}}$ К, начинается истечение вещества [22, 23]. Для звезды солнечного типа с возрастом 4.5 млрд. лет это будет происходить на радиусе орбиты $0.3$ а.е., для звезды возрастом 200 млн. лет – на радиусе 2 а.е. [22]. Качественно эта зависимость выглядит как $\dot {M} \propto $ $ \propto {{L}_{{{\text{XUV}}}}}{\text{/}}{{a}^{2}}$, где $a$ – радиус орбиты планеты [24], а ${{L}_{{{\text{XUV}}}}}$ – светимость звезды в экстремальном ультрафиолетовом диапазоне. Эта светимость убывает с возрастом звезды примерно как ${{L}_{{{\text{XUV}}}}} \propto {{t}^{{ - 1}}}$ [25]. В общем случае темп истечения имеет вид

(1)
$\dot {M} = {{\dot {M}}_{{{\text{ref}}}}}\mathop {\left( {\frac{t}{{4.6 \times {{{10}}^{9}}\;{\text{лет}}}}} \right)}\nolimits^{ - 1} \mathop {\left( {\frac{a}{{{{a}_{{{\text{ref}}}}}}}} \right)}\nolimits^{ - 2} ,$
где ${{\dot {M}}_{{{\text{ref}}}}}$ и ${{a}_{{{\text{ref}}}}}$ – некоторые референсные значения. Мы будем нормировать величину $\dot {M}$ на современный темп истечения атмосферы планеты HD 209458b. По анализу затменной кривой блеска в системе HD 209458 авторы работы [26] установили, что атмосфера теряет вещество со скоростью $\dot {M} = (3{\kern 1pt} - {\kern 1pt} 7) \times {{10}^{{10}}}$ г/с. Некоторые теоретические оценки приводили к более высоким значениям, порядка ${{10}^{{11}}}$ г/с [24] или ${{10}^{{12}}}$ г/с [22]. Для HD 209458b имеем ${{a}_{{{\text{ref}}}}} = 0.045$ а.е, ${{\dot {M}}_{{{\text{ref}}}}} = $ $ = (3 \times {{10}^{{10}}} - {{10}^{{12}}})$ г/с. Этой зависимости хорошо соответствует также планета в системе OGLE-TR-113 [24].

Основной поток вещества из атмосферы планеты будет выходить из окрестностей точек Лагранжа L1 и L2. При этом динамическое давление звездного ветра может запирать точку L1, либо фокусировать в направлении L2 поток, выходящий из L1 [27]. Газ, истекающий из атмосферы, спустя один оборот планеты попадет в поле действия ее гравитационной силы. В результате гравитационного рассеяния угловой момент вещества струи будет меняться. Угловой момент сил можно оценить в баллистическом приближении [2]:

(2)
${{\dot {J}}_{{{\text{gr}}}}} \sim \mathop {\left( {\frac{{{{M}_{{\text{p}}}}}}{{{{M}_{{\text{s}}}}}}} \right)}\nolimits^2 \mathop {\left( {\frac{a}{b}} \right)}\nolimits^3 \Sigma {{a}^{4}}\Omega _{{\text{p}}}^{2},$
где $b$ – радиальное расстояние от планеты до струи, $\Sigma $ – плотность струи газа, ${{\Omega }_{{\text{p}}}} = {{(G{{M}_{{\text{s}}}}{\text{/}}{{a}^{3}})}^{{1/2}}}$ – орбитальная угловая скорость планеты. Расстояние $b$, очевидно, порядка радиуса Хилла – масштаба полости Роша планеты,
(3)
${{r}_{{{\text{Hill}}}}} = \mathop {\left( {\frac{{{{M}_{{\text{p}}}}}}{{3{{M}_{{\text{s}}}}}}} \right)}\nolimits^{1/3} a \approx 0.0682a.$
В результате приливного взаимодействия за один орбитальный период вещество струи отбирает у планеты угловой момент
(4)
$\Delta J \sim {{\dot {J}}_{{{\text{gr}}}}}\Omega _{{\text{p}}}^{{ - 1}} \sim \mathop {\left( {\frac{{{{M}_{{\text{p}}}}}}{{{{M}_{{\text{s}}}}}}} \right)}\nolimits^2 \mathop {\left( {\frac{a}{b}} \right)}\nolimits^5 {{J}_{{{\text{str}}}}},$
где ${{J}_{{{\text{str}}}}} \sim \Sigma {{b}^{2}}{{\Omega }_{{\text{p}}}}{{a}^{2}}$ – угловой момент газа в струе [2].

Само по себе истечение вещества из планеты не способно существенно изменить ее угловой момент. Действительно, для того, чтобы планета потеряла весь свой угловой момент, необходимо порядка $N \equiv {{M}_{{\text{p}}}}{{a}^{2}}{{\Omega }_{{\text{p}}}}{\text{/}}\Delta J$ орбитальных периодов. Характерное время выметания вещества струи из окрестности орбиты планеты много меньше орбитального времени: ${{J}_{{{\text{str}}}}}{\text{/}}{{\dot {J}}_{{{\text{gr}}}}} \sim {{(b{\text{/}}a)}^{5}}\Omega _{{\text{p}}}^{{ - 1}} \ll \Omega _{{\text{p}}}^{{ - 1}}$. Следовательно, масса струи может накапливаться только в течение одного орбитального периода, ее можно оценить также как ${{m}_{{{\text{str}}}}} \equiv \dot {M}\Omega _{{\text{p}}}^{{ - 1}}$. За время $N\Omega _{{\text{p}}}^{{ - 1}}$ планета должна потерять массу

(5)
$N{{m}_{{{\text{str}}}}} \sim \mathop {\left( {\frac{{{{M}_{{\text{s}}}}}}{{{{M}_{{\text{p}}}}}}} \right)}\nolimits^2 \mathop {\left( {\frac{b}{a}} \right)}\nolimits^5 {{M}_{{\text{p}}}}.$
Подставляя сюда $b = {{r}_{{{\text{Hill}}}}}$, получим $N{{m}_{{{\text{str}}}}} \sim {{M}_{{\text{p}}}}$, т.е. для того, чтобы планета потеряла свой угловой момент, она должна потерять всю свою массу.

Совсем другая картина возникает, если допустить, что вещество струи, испытав гравитационное рассеяние планетой, накапливается на некоторой равновесной орбите вокруг звезды с радиусом $a + b + \Delta b$. Величину радиального смещения $\Delta b$ оценим в баллистическом приближении, с помощью (4) и формулы для кеплеровского момента: $\Delta J{\text{/}}{{J}_{{{\text{str}}}}} = {{[1 + \Delta b{\text{/}}(a + b)]}^{{1/2}}} - 1$. Полагая $b = {{r}_{{{\text{Hill}}}}}$, получим $\Delta b = 1.81a$. В действительности самопересечение струи должно приводить к диссипации части углового момента и циркуляризации ее на гораздо меньшем радиусе.

Самопересечение струи, вероятно, приведет также к возбуждению в газе турбулентности и, следовательно, турбулентной вязкости. Под действием вязкости вещество начнет “расплываться” по радиусу до тех пор, пока не установится баланс гравитационной и вязкой силы. В газовом протопланетном диске этот механизм описывает формирование щели. Угловой момент вязких сил можно оценить как [2]

(6)
${{\dot {J}}_{{{\text{visc}}}}} \sim {{\nu }_{{\text{t}}}}\Sigma {{a}^{2}}{{\Omega }_{{\text{p}}}}.$
Полагая распределение газа равновесным или медленно меняющимся и приравнивая выражения (2) и (6), получим оценку коэффициента кинематической турбулентной вязкости:
(7)
${{\nu }_{{\text{t}}}} \sim \mathop {\left( {\frac{{{{M}_{{\text{p}}}}}}{{{{M}_{{\text{s}}}}}}} \right)}\nolimits^2 \mathop {\left( {\frac{a}{{{{H}_{{{\text{gap}}}}}}}} \right)}\nolimits^3 {{a}^{2}}{{\Omega }_{{\text{p}}}},$
где ${{H}_{{{\text{gap}}}}}$ – полуширина щели, так мы переобозначили параметр $b$ в уравнении (2). Если параметризовать вязкость $\alpha $ параметром Шакуры–Сюняева [28] как ${{\nu }_{{\text{t}}}} = \alpha {{H}^{2}}{{\Omega }_{{\text{p}}}}$, то
(8)
$\alpha \sim \mathop {\left( {\frac{{{{M}_{{\text{p}}}}}}{{{{M}_{{\text{s}}}}}}} \right)}\nolimits^2 \mathop {\left( {\frac{a}{{{{H}_{{{\text{gap}}}}}}}} \right)}\nolimits^3 \mathop {\left( {\frac{a}{H}} \right)}\nolimits^2 ,$
где $H$ – вертикальная шкала распределения газа, ее можно оценить в изотермическом приближении как $H = {{c}_{{\text{s}}}}{\text{/}}{{\Omega }_{{\text{p}}}}$ или
(9)
$H = (26.7{{R}_{{{\text{Jup}}}}})\mathop {\left( {\frac{T}{{{{{10}}^{4}}\;{\text{К}}}}} \right)}\nolimits^{1/2} \mathop {\left( {\frac{a}{{0.1\;{\text{а}}.{\text{е}}.}}} \right)}\nolimits^{3/2} $
или
(10)
$\frac{H}{a} = 0.125\,\mathop {\left( {\frac{T}{{{{{10}}^{4}}\;{\text{К}}}}} \right)}\nolimits^{1/2} \mathop {\left( {\frac{a}{{0.1\;{\text{а}}.{\text{е}}.}}} \right)}\nolimits^{1/2} ,$
где $T$ – температура газа.

Трудно дать надежную оценку для полуширины щели ${{H}_{{{\text{gap}}}}}$. С одной стороны, ${{H}_{{{\text{gap}}}}}$ не может быть меньше радиуса гравитационного захвата планеты – радиуса Хилла (3). С другой стороны, она определяется величиной вязкости и структурой распределения вещества. Часто для грубых оценок полагают, что ширина щели равна вертикальной шкале распределения газа [1, 2]. Численные расчеты, проведенные в работе [29], показали, что для планеты с массой Юпитера полуширина щели на половине ее глубины составляет ${{H}_{{{\text{gap}}}}} \approx 0.4a$, что почти на порядок превышает вертикальную шкалу диска11.

Накопление и расплывание вещества должно привести к формированию тора или диска. Химический состав диска будет идентичен составу газа, истекающего из экзосферы планеты. Кроме того, вещество диска будет подвергаться воздействию такого же по величине потока экстремального УФ-излучения от звезды. Если принять, для примера, что плотность диска по порядку величины равна плотности газа на уровне экзобазы в атмосфере планеты, $ \gtrsim {\kern 1pt} {{10}^{{10}}}$ см–3, то температура будет порядка ${{10}^{4}}$ К; для плотности ${{10}^{7}}$ см–3 температура окажется примерно $7 \times {{10}^{3}}$ К [24].

В поле экстремального УФ излучения звезды диск, как и атмосфера планеты, будет оптически толстым [30]. Однако для более длинноволнового излучения, при температуре $ \gtrsim {\kern 1pt} {{10}^{4}}$ К и концентрации $N = {{10}^{{10}}}$ см–3, основной вклад в поглощение будет вносить рассеяние на свободных электронах с коэффициентом непрозрачности ${{\kappa }^{{{\text{es}}}}}$ = = 0.4 см2/г. На вертикальной шкале распределения газа оптическая толщина будет порядка ${{\kappa }^{{{\text{es}}}}}\rho H \sim 0.001$ (при условии полной ионизации, которое заведомо не выполняется при такой большой плотности), а на радиальной шкале – лишь на порядок больше.

Плотность вещества в диске будет определяться скоростью расплывания газа в радиальном направлении, т.е. величиной вязкости. Заметим, что выражение (7) дает лишь верхнюю допустимую границу значений коэффициента вязкости, при условии открытия щели, и если структура газа полагается двумерной. Последнее условие можно записать как ${{r}_{{{\text{Hill}}}}} \gtrsim H$ [1] (см. также результаты численных расчетов [29]). Для температуры $T = {{10}^{4}}$ К это не выполняется. Это означает, что при достигнутом балансе сил (2) и (6) газ может проникать в щель, приводя к эффективному уменьшению ее ширины.

Как видно, и полуширина щели, и коэффициент вязкости определены довольно плохо. Так, полагая $T = {{10}^{4}}$ К и ${{H}_{{{\text{gap}}}}} = H$, из выражения (8) будем иметь $\alpha = 0.033$, а принимая ${{H}_{{{\text{gap}}}}}{\text{/}}a = 0.4$, как в работе [29], получим $\alpha = 0.001$. Недавние наблюдения протопланетных дисков дают верхние значения $\alpha \lesssim 0.003{\kern 1pt} - {\kern 1pt} 0.007$ [31, 32], при том, что температура газа там существенно ниже ${{10}^{4}}$ К.

Исходя из приведенных выше соображений относительно ширины щели и прочего, можем оценить характерное время миграции планеты с угловым моментом ${{J}_{{\text{p}}}} = {{M}_{{\text{p}}}}{{a}^{2}}{{\Omega }_{{\text{p}}}}$:

(11)
$\frac{{{{J}_{{\text{p}}}}}}{{{{{\dot {J}}}_{{{\text{gr}}}}}}} \sim ({{10}^{9}}\;{\text{лет}})\mathop {\left( {\frac{\Sigma }{{0.01\;{\text{г/с}}{{{\text{м}}}^{2}}}}} \right)}\nolimits^{ - 1} \mathop {\left( {\frac{{{{H}_{{{\text{gap}}}}}}}{H}} \right)}\nolimits^3 \frac{a}{{0.1\;{\text{а}}.{\text{е}}.}}.$
Если моделирование турбулентности до сих пор является непосильной задачей, то согласованный расчет формирования щели вполне может быть выполнен. Сформулируем задачу. Нас интересует распределение плотности, которое образуется в результате истечения атмосферы планеты. Истекающее вещество подвержено, с одной стороны, приливному отталкиванию под действием гравитации планеты, а с другой стороны – вязкой диффузии. Аккреция вещества извне отсутствует. Требуется рассчитать эволюцию распределения плотности, угловой момент силы, действующий на вещество и эволюцию орбиты планеты в результате потери углового момента.

3. МОДЕЛЬ ДИСКА

3.1. Перенос массы

Мы будем считать распределение вещества вокруг звезды осесимметричным и геометрически тонким. Предположим также, что перенос массы полностью контролируется переносом углового момента. Течение с такими свойствами описывается моделью Прингла [33]:

(12)
$\frac{{\partial \Sigma }}{{\partial t}} + \frac{1}{r}\frac{{\partial (rF)}}{{\partial r}},$
(13)
$\frac{\partial }{{\partial t}}(\Sigma r\Omega ) + \frac{1}{{{{r}^{2}}}}\frac{{\partial ({{r}^{2}}Fr\Omega )}}{{\partial r}} = \frac{1}{{{{r}^{2}}}}\frac{{\partial ({{r}^{2}}W)}}{{\partial r}} + \Sigma \tau .$
Здесь $F$ – радиальный поток массы; $\Omega $ – распределение угловой скорости вещества, предполагается заданным; $W$ – компонент “$r\phi $” тензора турбулентных напряжений, или тензора Рейнольдса. В правую часть уравнения (13) добавлен угловой момент поверхностной силы $\Sigma \tau $, который учитывает гравитационное взаимодействие вещества и планеты.

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

(14)
$W = {{\nu }_{{\text{t}}}}\Sigma r\frac{{\partial \Omega }}{{\partial r}},$
где ${{\nu }_{{\text{t}}}}$ – турбулентная кинематическая вязкость. Если распределение угловой скорости стационарно, можно выразить поток массы из уравнения (13),
(15)
$F = \mathop {\left[ {\frac{{\partial ({{r}^{2}}\Omega )}}{{\partial r}}} \right]}\nolimits^{ - 1} \left\{ {\frac{1}{r}\frac{{\partial ({{r}^{2}}W)}}{{\partial r}} + \Sigma \tau } \right\}.$
Приливное воздействие планеты выметает газ из ближайшей окрестности ее орбиты, образуя щель. Массовая плотность момента приливной силы может быть записана в виде [3, 36]:
(16)
$\tau = \frac{{{{C}_{0}}}}{\pi }\frac{{{{G}^{2}}M_{{\text{p}}}^{2}a}}{{{{{(r - a)}}^{2}}}}\frac{{{{\Omega }_{{\text{p}}}} - \Omega }}{{{{{(a{{\Omega }_{{\text{p}}}} - r\Omega )}}^{3}}}},$
где ${{C}_{0}} = (4{\text{/}}9){{[2{{K}_{0}}(2{\text{/}}3) + {{K}_{1}}(2{\text{/}}3)]}^{2}} \approx 2.82$. В выражении для потока (15) учет приливной силы приводит к появлению неотрицательной добавки к потоку, которая, впрочем, убывает с радиусом довольно быстро, как $\Sigma {\text{/}}{{r}^{{3/2}}}$, для $r \gg a$ и кеплеровского распределения угловой скорости.

Истечение атмосферы начинается на расстоянии порядка радиуса сферы Хилла от планеты (3):

(17)
${{r}_{0}} = a + {{r}_{{{\text{Hill}}}}} = \left[ {1 + \mathop {\left( {\frac{{{{M}_{{\text{p}}}}}}{{3{{M}_{{\text{s}}}}}}} \right)}\nolimits^{1/3} } \right]a.$
На радиусе ${{r}_{0}}$ будет располагаться внутренняя граница расчетной области и там будет задан поток вещества
(18)
${{F}_{0}} = \frac{{\dot {M}}}{{2\pi {{r}_{0}}}}.$
Поскольку темп истечения атмосферы $\dot {M}$ зависит от времени как явным образом, так и через радиус орбиты планеты (см. (1)), поток массы будет нестационарным.

Будем полагать также, что угловой момент поступающего вещества равен угловому моменту газа в диске на радиусе ${{r}_{0}}$. Распределение углового момента положим всюду кеплеровским, $\Omega = \sqrt {GM{\text{/}}{{r}^{3}}} $. Поток массы в этом случае примет вид

(19)
$\begin{gathered} F = - \frac{3}{{{{r}^{{1/2}}}}}\frac{\partial }{{\partial r}}({{r}^{{1/2}}}{{\nu }_{{\text{t}}}}\Sigma ) + \\ + \;\frac{{2{{C}_{0}}}}{\pi }\mathop {\left( {\frac{{{{M}_{{\text{p}}}}}}{{{{M}_{{\text{s}}}}}}} \right)}\nolimits^2 \frac{{{{{(G{{M}_{{\text{s}}}}r)}}^{{1/2}}}a}}{{{{{(r - a)}}^{2}}}}\frac{{{{r}^{{3/2}}} - {{a}^{{3/2}}}}}{{{{{({{r}^{{1/2}}} - {{a}^{{1/2}}})}}^{3}}}}\Sigma . \\ \end{gathered} $
Приливное взаимодействие планеты и вещества тора приводит к уменьшению углового момента планеты. Если планета обращается по кеплеровской орбите, взаимодействие приведет к тому, что ее орбита будет мигрировать со скоростью
(20)
$\dot {a} = - \frac{{2a}}{{{{J}_{{\text{p}}}}}}2\pi \int\limits_0^\infty {dr\,r\Sigma \tau } ,$
где ${{J}_{{\text{p}}}} = {{M}_{{\text{p}}}}\sqrt {G{{M}_{{\text{s}}}}a} $.

Совместная эволюция планетной орбиты и тора описывается совокупностью уравнений (12) и (17)–(20). В качестве начальных условий примем нулевое распределение плотности. При численном решении этой системы могут возникнуть сложности, вызванные тем, что расчетная область ($r \geqslant {{r}_{0}}$) имеет подвижную внутреннюю границу. Для того, чтобы этого избежать, приведем уравнения к безразмерному виду с помощью преобразований $t = {{t}_{0}}s$, $r = {{r}_{0}}x$, $\Sigma = {{\Sigma }_{0}}\sigma $, $F = {{F}_{0}}f$, ${{\nu }_{{\text{t}}}} = {{\nu }_{0}}n$, $a = {{r}_{0}}\xi $. Если принять определения

(21)
$\begin{gathered} {{t}_{0}} = \frac{1}{{{{\Omega }_{0}}}},\quad {{\Sigma }_{0}} = \frac{{\dot {M}}}{{2\pi r_{0}^{2}{{\Omega }_{0}}}}, \\ {{\nu }_{0}} = r_{0}^{2}{{\Omega }_{0}},\quad {{\Omega }_{0}} = \Omega ({{r}_{0}}), \\ \end{gathered} $
то эти уравнения примут полностью безразмерный вид:
(22)
$\frac{{\partial \sigma }}{{\partial s}} + \frac{1}{x}\frac{{\partial (xf)}}{{\partial x}} + \frac{{{{t}_{0}}{{{\dot {\Sigma }}}_{0}}}}{{{{\Sigma }_{0}}}}\sigma = 0,$
(23)
$\begin{gathered} f = - \frac{3}{{{{x}^{{1/2}}}}}\frac{\partial }{{\partial x}}({{x}^{{1/2}}}n\sigma ) + \\ + \;\frac{{2{{C}_{0}}}}{\pi }\mathop {\left( {\frac{{{{M}_{{\text{p}}}}}}{{{{M}_{{\text{s}}}}}}} \right)}\nolimits^2 \frac{{\xi {{x}^{{1/2}}}}}{{{{{(x - \xi )}}^{2}}}}\frac{{{{x}^{{3/2}}} - {{\xi }^{{3/2}}}}}{{{{{({{x}^{{1/2}}} - {{\xi }^{{1/2}}})}}^{3}}}}\sigma , \\ \end{gathered} $
(24)
$\frac{{da}}{{ds}} = {{t}_{0}}a\frac{{{{M}_{{\text{p}}}}\dot {M}}}{{M_{{\text{s}}}^{2}}}\frac{{2{{C}_{0}}}}{\pi }\int\limits_1^\infty \,dx\frac{{{{\xi }^{{1/2}}}x}}{{{{{(x - \xi )}}^{2}}}}\frac{{{{x}^{{3/2}}} - {{\xi }^{{3/2}}}}}{{{{{({{x}^{{1/2}}} - {{\xi }^{{1/2}}})}}^{3}}}}\sigma ,$
(25)
$s \geqslant {{t}_{{{\text{ini}}}}}{\text{/}}{{t}_{0}},\quad x \geqslant 1,$
(26)
${{\left. \sigma \right|}_{{s = {{t}_{{{\text{ini}}}}}/{{t}_{0}}}}} \equiv 0,\quad {{\left. f \right|}_{{x = 1}}} \equiv 1,\quad {{\left. a \right|}_{{s = {{t}_{{{\text{ini}}}}}/{{t}_{0}}}}} = {{a}_{{{\text{ini}}}}}.$
Как видно, в уравнении непрерывности (22) появился дополнительный член, вызванный тем, что масштабный множитель ${{\Sigma }_{0}}$ зависит от времени посредством ${{r}_{0}}(a(t))$. Физическое решение записывается как $\Sigma (t,r) = {{\Sigma }_{0}}\sigma ({{\Omega }_{0}}t,r{\text{/}}{{r}_{0}})$. Заметим, что изменение массы планеты предполагается малым. По этой причине величины ${{M}_{{\text{p}}}}$ и $\xi = a{\text{/}}{{r}_{0}}$ принимаются константами.

Можно показать, что модель, описываемая уравнениями (22)(26), не имеет физически корректного стационарного решения. Действительно, проинтегрируем выражение для плотности потока массы (23), учитывая то, что в стационарном пределе поток не будет зависеть от радиальной координаты, или $xf \equiv 1$:

(27)
$n\sigma = \frac{2}{3}\left( {\frac{{{{X}^{{1/2}}}}}{{{{x}^{{1/2}}}}} - 1} \right),$
где
(28)
${{X}^{{1/2}}} = \frac{3}{2}n(1)\sigma (1) + 1 + \int\limits_1^x {dx{\kern 1pt} '\frac{{\omega \sigma }}{{2x{\kern 1pt} {{'}^{{1/2}}}}}} ,$
$n(1)$ и $\sigma (1)$ – значения на внутренней границе, $x = 1$; через $\omega \sigma $ обозначено второе слагаемое в правой части (23). Подынтегральная функция в правой части последнего выражения убывает довольно быстро с расстоянием, и величина $X$ стремится к некоторому постоянному значению ${{x}_{*}}$. На радиусе $x = {{x}_{*}}$ комбинация $n\sigma $ обращается в нуль и далее становится отрицательной. Таким образом, задачу о переносе массы необходимо рассматривать в нестационарной постановке.

3.2. Турбулентность

В модели Прингла перенос массы определяется механизмом переноса углового момента, который зависит от тензора турбулентных вязких напряжений. Если рассматривать неустойчивости в качестве причин турбулентности, то можно отметить, что диск должен быть оптически тонким в континууме, как мы отметили в п. 2. Тогда основным источником тепловой энергии будет выступать УФ излучение звезды, которое поглощается нейтральными атомами водорода. Непрозрачность газа при этом падает с ростом температуры. Как видно, в диске отсутствуют условия для вертикальной и радиальной конвективной неустойчивости (последняя необходима для реализации бароклинной неустойчивости [37]). Конечная толщина диска оставляет потенциальную возможность для развития вертикальной сдвиговой неустойчивости [38] или магниторотационной неустойчивости при наличии магнитного поля. Возбуждение в диске приливных волн достаточно большой амплитуды может привести к развитию чисто гидродинамической неустойчивости [39]. Кроме того, даже нелинейные звуковые волны во вращающемся течении приводят к развитию турбулентности [40].

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

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

(29)
$\frac{{{{a}^{2}}}}{{{{\nu }_{{\text{t}}}}}} \sim (300\;{\text{лет}})\mathop {\left( {\frac{\alpha }{{0.001}}} \right)}\nolimits^{ - 1} \mathop {\left( {\frac{a}{{0.1\;{\text{а}}.{\text{е}}.}}} \right)}\nolimits^{1/2} .$
Поскольку уравнение переноса массы не может дать стационарного решения (см. раздел 3.1), изменение профиля плотности происходит на гораздо большем масштабе времени. Кроме того, мы ожидаем, что турбулентность будет существенно дозвуковой. Все это позволяет нам рассматривать перенос турбулентности как быстрый процесс, независимый от переноса массы. Далее будем полагать распределение турбулентности стационарным.

Коэффициент кинематической турбулентной вязкости можно записать следующим образом:

(30)
${{\nu }_{{\text{t}}}} = \frac{{{{\ell }^{2}}}}{\tau },$
где $\ell $ – пространственный масштаб корреляции турбулентности (он связан с размером наиболее крупного вихря); $\tau $ – время корреляции. Если принять классическую параметризацию вязкости по модели Шакуры и Сюняева [28], ${{\nu }_{{\text{t}}}} = \alpha {{H}^{2}}\Omega $, то безразмерный коэффициент $\alpha $ можно записать в виде
(31)
$\alpha = \mathop {\left( {\frac{\ell }{H}} \right)}\nolimits^2 \frac{1}{{\Omega \tau }}.$
В астрофизических дисках азимутальное течение является сверхзвуковым, с характерным временем переноса в азимутальном направлении порядка кеплеровского времени ${{\Omega }^{{ - 1}}}$, оно же вертикальное звуковое время. Распространение турбулентности в азимутальном направлении должно происходить в этой шкале времени. Следовательно, именно локальное кеплеровское время нужно принять в качестве характерного времени корреляции турбулентности [41]. Другими словами, в выражении (31) необходимо положить $\Omega \tau \sim 1$.

Если предполагать, что единственным источником турбулентности в диске может быть только неустойчивость Кельвина–Гельмгольца на его внутренней границе (см. выше), то течение будет эквивалентно так называемому свободному турбулентному течению в плоской струе, за тем исключением, что в нашем случае имеется радиальный градиент азимутальной скорости. Лабораторные эксперименты по течениям такого типа обнаруживают, что пространственный масштаб корреляции пропорционален толщине турбулентной области [42], которая в нашем случае, очевидно, будет порядка толщины диска. На основании этого положим $\ell \propto H$, следовательно, $\alpha \approx {\text{const}}$. При этом радиальное распределение коэффициента вязкости будет определяться зависимостью от радиуса величины ${{H}^{2}}\Omega \propto T{\text{/}}\Omega $, где $T$ – температура газа.

Воспользуемся простым модельным представлением для безразмерного коэффициента вязкости, который будем применять при решении уравнений (22)–(26). Положим, что полутолщина зависит от радиуса степенным образом, тогда безразмерный коэффициент вязкости для модели (22)–(26) примет вид

(32)
$n = \frac{{{{\nu }_{{\text{t}}}}}}{{r_{0}^{2}{{\Omega }_{0}}}} = \alpha {{h}^{2}}{{x}^{\beta }},$
где обозначено $h = H({{r}_{0}}){\text{/}}{{r}_{0}}$. В рамках этого представления можно описать и те случаи, когда турбулентность порождается и в объеме газа (см. выше).

4. ЭВОЛЮЦИЯ ТОРА И ОРБИТЫ ПЛАНЕТЫ

Модель переноса массы (22)–(26) с законом вязкости (32) зависит от нескольких параметров: массы звезды и планеты, темпа потери массы планетой $\dot {M}$, температуры газа на внутреннем радиусе (от нее зависит относительная полутолщина диска $h$), параметра турбулентности $\alpha $ на внутреннем радиусе диска и показателя степени $\beta $ в распределении коэффициента вязкости по радиусу.

Для расчета модели была введена пространственная сетка и использована стандартная дискретизация дифференциального оператора. Плотность и коэффициент вязкости задавались в центрах ячеек с координатами ${{x}_{{i + 1/2}}}$, потоки – в узлах ячеек с координатами ${{x}_{i}}$. В итоге была получена система ОДУ, для решения которой применялся неявный метод33. Численная схема выглядит следующим образом44:

(33)
$\frac{{d{{\sigma }_{{i + 1/2}}}}}{{ds}} + \frac{1}{{{{x}_{{i + 1/2}}}}}\frac{{{{x}_{{i + 1}}}{{f}_{{i + 1}}} - {{x}_{i}}{{f}_{i}}}}{{{{\Delta }_{{i + 1/2}}}}} = 0,$
(34)
$\begin{gathered} {{f}_{i}} = - \frac{3}{{x_{i}^{{1/2}}{{\Delta }_{i}}}}(x_{{i + 1/2}}^{{1/2}}{{n}_{{i + 1/2}}}{{\sigma }_{{i + 1/2}}} - x_{{i - 1/2}}^{{1/2}}{{n}_{{i - 1/2}}}{{\sigma }_{{i - 1/2}}}) + \\ \, + {{\omega }_{i}}\frac{{{{\sigma }_{{i - 1/2}}} + {{\sigma }_{{i + 1/2}}}}}{2}, \\ \end{gathered} $
(35)
${{\omega }_{i}} = \frac{{2{{C}_{0}}}}{\pi }\mathop {\left( {\frac{{{{M}_{{\text{p}}}}}}{{{{M}_{{\text{s}}}}}}} \right)}\nolimits^2 \frac{{\xi x_{i}^{{1/2}}}}{{{{{({{x}_{i}} - \xi )}}^{2}}}}\frac{{x_{i}^{{3/2}} - {{\xi }^{{3/2}}}}}{{{{{(x_{i}^{{1/2}} - {{\xi }^{{1/2}}})}}^{3}}}},$
(36)
${{x}_{{i + 1/2}}} = ({{x}_{i}} + {{x}_{{i + 1}}}){\text{/}}2,$
(37)
${{\Delta }_{{i + 1/2}}} = {{x}_{{i + 1}}} - {{x}_{i}},\quad {{\Delta }_{i}} = ({{\Delta }_{{i - 1/2}}} + {{\Delta }_{{i + 1/2}}}){\text{/}}2.$
Начальные и граничные условия:
(38)
$\sigma {{{\text{|}}}_{{s = 0}}} = 0,$
(39)
${{x}_{0}} = 1,\quad {{f}_{0}} = 1,$
(40)
${{x}_{N}}{{f}_{N}} = {{x}_{{N - 1}}}{{f}_{{N - 1}}}.$
Во всех расчетах использовалась сетка на $N = 1000$ узлов с логарифмическим шагом. Для того, чтобы уменьшить влияние граничных условий на внешней границе расчетной области, координата последнего узла сетки была установлена равной ${{x}_{N}} = ({{10}^{4}}\;{\text{а}}.{\text{е}}.)/{{r}_{0}}({{t}_{{{\text{ini}}}}})$. Начальный шаг по времени, для численного интегрирования, соответствовал кеплеровскому орбитальному периоду планеты.

Профиль полутолщины диска, $H(r)$, задается балансом нагрева газа со стороны звезды и турбулентной вязкости, и радиативного охлаждения. В стандартном оптически толстом $\alpha $-диске, если непрозрачность в континууме определяется электронным рассеянием, температура от радиуса зависит как ${{r}^{{ - 9/10}}}$ [45], чему будет соответствовать показатель наклона $\beta = 3{\text{/}}5$ в коэффициенте вязкости. Этот случай вряд ли реализуется в нашей задаче, поскольку при температуре порядка ${{10}^{4}}$ К и концентрации ${{10}^{{10}}}$ см–3 газ не будет ионизован полностью, а значит, как было отмечено в п. 2, диск будет оптически тонким. С другой стороны, если источником турбулентности в торе является только область его внутренней границы, то коэффициент вязкости будет пропорционален отношению $T{\text{/}}\Omega $, что в случае изотермического распределения даст $\beta = 3{\text{/}}2$. В действительности температура газа, вероятнее всего, падает с радиусом, поэтому профиль распределения вязкости будет более пологим. Для полноты мы рассмотрим три случая: $\beta = 3{\text{/}}2$, $\beta = 0$ (постоянный по радиусу коэффициент вязкости) и промежуточный случай, $\beta = 1$. Также возьмем два значения параметра турбулентности $\alpha $: ${{10}^{{ - 3}}}$ (условно – слабая турбулентность) и ${{10}^{{ - 2}}}$ (сильная турбулентность), которые соответствуют наблюдательным оценкам.

Начальный момент времени в модели, ${{t}_{{{\text{ini}}}}}$, соответствовал возрасту звезды ${{10}^{7}}$ лет, т.е. расчет начинался после того, как исчез газовый протопланетный диск. Начальный радиус орбиты планеты подбирался в каждом расчете так, чтобы к настоящему моменту (возраст звезды $4.6\;\; \times $ ×  109 лет) планета мигрировала на орбиту 0.045 а.е. Во всех расчетах температура газа предполагалась равной 104 К. Параметры всех моделей приведены в табл. 1.

Таблица 1.

Параметры моделей

$\alpha $ $\beta $ ${{a}_{{{\text{ini}}}}}$, а.е.
0.001 0 0.4106
0.001 1 0.6690
0.001 1.5 0.6739
0.01 0 0.5338
0.01 1 0.6734
0.01 1.5 0.6739

Примечание. $\alpha $ – интенсивность турбулентности, $\beta $ – показатель наклона профиля турбулентной вязкости, ${{a}_{{{\text{ini}}}}}$ – начальный радиус орбиты планеты.

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

Рис. 1.

Эволюция профиля поверхностной плотности на моменты времени $1$, $10$, $100$, $1000$ и $10{\kern 1pt} {\kern 1pt} 000$ лет. Профили показаны в логарифмической шкале (верхний рисунок) и в линейной шкале (нижний рисунок). Планета расположена на орбите $0.6739$ а.е. (нижний левый угол рисунка).

Результаты расчетов на финальный момент времени, для всех комбинаций параметров, приведены на рис. 2. Как видно, в окрестности планеты образуется щель, за которой поверхностная плотность возрастает до максимального значения и далее монотонно убывает. Таким образом, распределение вещества имеет вид тора. Ширина тора уменьшается с ростом $\alpha $ (с усилением турбулентности) и с увеличением показателя наклона $\beta $. В случае $\beta = 0$ ширина тора, если оценивать ее на половине высоты распределения плотности, составляет 1 а.е для слабой турбулентности и $0.5$ а.е. для сильной турбулентности. С ростом показателя наклона ширина тора уменьшается. Так, при $\alpha = 0.1$ и $\beta = 1.5$ тор имеет ширину примерно $0.1$ а.е. Максимальное значение поверхностной плотности для сильной турбулентности примерно в пять раз ниже, чем в случае слабой турбулентности. Ширина щели также уменьшается с ростом турбулентной вязкости. Это связано, по-видимому, с тем, что больший коэффициент вязкости приводит к более сильной зависимости потока массы от градиента плотности. По этой причине внешняя область тора подвержена более сильному оттоку вещества, а щель закрывается быстрее.

Рис. 2.

Профили поверхностной плотности в торе на момент, когда орбита планеты мигрировала к радиусу $0.045$ а.е. Все расчеты начинались с возраста звезды ${{t}_{{{\text{ini}}}}} = {{10}^{7}}$ лет. Начальные орбиты планеты, в зависимости от параметров распределения вязкости, $\alpha $ и $\beta $, приведены в табл. 1.

Любопытно, что начальные радиусы орбиты планеты практически совпадают в тех расчетах, где $\beta \geqslant 1$, как для слабой, так и для сильной турбулентности. Это значит, что во всех этих случаях, несмотря на различие в распределении поверхностной плотности, темп обмена угловым моментом был одинаковым.

5. ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Ключевой элемент модели, представленной в данной работе, это предположение о том, что вещество, истекающее из атмосферы горячего юпитера, способно накапливаться на орбите вокруг звезды, вблизи орбиты планеты. В разделе 2 было показано, что одного лишь истечения атмосферы не достаточно для отвода углового момента – вызванная истечением “реактивная сила” оказывается слишком мала. Если же допустить, что газ может взаимодействовать с планетой долгое время, обмен угловым моментом будет заметным.

Накоплению вещества в торе способствуют два фактора. Во-первых, темп истечения растет по мере того, как планета приближается к звезде (см. (1) и [24]), а скорость миграции (24) пропорциональна темпу истечения. Во-вторых, избыточный угловой момент газа, который был приобретен путем гравитационного воздействия планеты, диссипирует под действием турбулентной вязкости, и газ стремится заполнить щель. С другой стороны, турбулентность приводит к расплыванию диска и уменьшению его плотности. Сильнее всего на профиль плотности оказывает влияние интенсивность турбулентности $\alpha $ и показатель наклона закона вязкости $\beta $. Для реалистичных значений этих параметров ($\alpha = 0.01{\kern 1pt} - {\kern 1pt} 0.001$, $\beta \geqslant 1$) оказалось, что за время с момента исчезновения протопланетного диска планета успевает мигрировать с орбиты примерно $0.67$ а.е. до орбиты $0.045$ а.е.

В нашей модели мы учли далеко не все факторы, которые могут повлиять на скорость миграции планеты. Например, на расстоянии от Солнца $a < 10\;{{R}_{ \odot }} \approx 0.046$ а.е. находится зона ускорения солнечного ветра излучением звезды [46]. Когда орбита горячего юпитера станет много меньше этого радиуса, звездный ветер может остановить движение внутренней границы тора, тем самым эффективно увеличивая щель и снижая темп обмена угловым моментом.

На низких орбитах вокруг звезды, где орбитальный период планеты становится меньше периода обращения звезды, их приливное взаимодействие приводит к диссипации орбитального углового момента [47], а следовательно, к более быстрой миграции. Для звезд типа T Tau период вращения составляет несколько дней [48], т.е. радиус орбиты должен быть меньше $ \sim {\kern 1pt} 0.03$ а.е. Период вращения современного Солнца составляет $24.5$ дня, это радиус коротации примерно $0.16$ а.е. Однако лишь на расстоянии около $0.03$ а.е. характерное время приливной миграции становится меньше миллиарда лет для планеты с массой Юпитера. Кроме того, приливные эффекты могут напрямую влиять на скорость истечения атмосферы планеты. В работе [24] было показано, что на расстоянии $0.015$ а.е. от звезды темп истечения может возрастать на порядок величины по сравнению с трендом (1).

Согласно общепринятым представлениям, по прошествии не более ${{10}^{6}}{\kern 1pt} - {\kern 1pt} {{10}^{7}}$ лет в газовом протопланетном диске формируются планетезимали [1]. Хотя газ, покидающий внешние слои атмосферы планеты, несколько обеднен тяжелыми элементами, по сравнению с солнечным составом [24], характерные времена, возникающие в нашей модели, достаточно велики для того, чтобы проявилась роль пыли как стока для газа, а также роль планетезималей как источника динамических возмущений.

Фотоиспарение считается одной из причин исчезновения газовых протопланетных дисков [49]. В этом случае вещество испытывает на себе действие экстремального ультрафиолета от звезды, находящейся на стадии классической Т Тельца. Поток ионизующего излучения при этом часто оценивают как ${{L}_{{{\text{XUV}}}}} \sim ({{10}^{{ - 3}}}{\kern 1pt} - {\kern 1pt} {{10}^{{ - 2}}}){{L}_{ \odot }}$. Для современного Солнца соответствующий поток на три-четыре порядка меньше, но при этом мы в нашей постановке задачи располагаем гораздо большим временем для этого эффекта. Заметим, что согласно гидродинамическим расчетам фотоиспарение следует учитывать на расстояниях от звезды не ближе $2$ а.е. [1]. В наших расчетах для $\beta \geqslant 1$ основная масса вещества в торе находится внутри этого радиуса, примерно на $0.1{\kern 1pt} - {\kern 1pt} 1$ а.е., поэтому для реалистичных параметров модели эффекта фотоиспарения может не быть.

Список литературы

  1. P. J. Armitage, Astrophysics of Planet Formation (Cambridge University Press, 2009).

  2. S. H. Lubow and S. Ida, in Exoplanets, edited by S. Seager (Univ. Arizona Press, 2009).

  3. P. Goldreich and S. Tremaine, Astrophys. J. 241, 425 (1980).

  4. P. Artymowicz, Astrophys. J. 419, 155 (1993).

  5. D. N. C. Lin and J. Papaloizou, Astrophys. J. 309, 846 (1986).

  6. S. Nayakshin and G. Lodato, Monthly Not. Roy. Astron. Soc. 426, 70 (2012).

  7. F. S. Masset and J. C. B. Papaloizou, Astrophys. J. 588, 494 (2003).

  8. K. E. Haisch, Jr., E. A. Lada, and C. J. Lada, Astrophys. J. 553, L153 (2001).

  9. P. J. Armitage, arXiv:1509.06382 astro-ph.SR. (2015).

  10. E. Podlewska-Gaca, J. C. B. Papaloizou, and E. Szuszkiewicz, Monthly Not. Roy. Astron. Soc. 421, 1736 (2012).

  11. L. Tu, C. P. Johnstone, M. Güdel, and H. Lammer, Astron. and Astrophys. 577, id. L3 (2015).

  12. D. Bisikalo, P. Kaygorodov, D. Ionov, V. Shematovich, H. Lammer, and L. Fossati, Astrophys. J. 764, id. 19 (2013).

  13. A. A. Cherenkov, D. V. Bisikalo, and A. G. Kosovichev, Monthly Not. Roy. Astron. Soc. 475, 605 (2018).

  14. A. A. Cherenkov, I. F. Shaikhislamov, D. V. Bisikalo, V. I. Shematovich, L. Fossati, and C. Möstl, Astron. Rep. 63, 94 (2019).

  15. A. G. Zhilkin and D. V. Bisikalo, Astron. Rep. 63, 550 (2019).

  16. A. G. Zhilkin, D. V. Bisikalo, and P. V. Kaygorodov, Astron. Rep., 64, 159 (2020).

  17. A. A. Cherenkov, D. V. Bisikalo, and P. V. Kaigorodov, Astron. Rep. 58, 679 (2014).

  18. T. Matsakos, A. Uribe, and A. Königl, Astron. and Astrophys. 578, id. A6 (2015).

  19. M. L. Khodachenko, I. F. Shaikhislamov, H. Lammer, A. G. Berezutsky, I. B. Miroshnichenko, M. S. Rumen-skikh, K. G. Kislyakova, and N. K. Dwivedi, Astrophys. J. 885, id. 67 (2019).

  20. I. F. Shaikhislamov, M. L. Khodachenko, H. Lammer, A. G. Berezutsky, I. B. Miroshnichenko, and M. S. Rumenskikh, Monthly Not. Roy. Astron. Soc. 491, 3435 (2020).

  21. A. Debrecht, J. Carroll-Nellenback, A. Frank, E. G. Black-man, L. Fossati, J. McCann, and R. Murray-Clay, Monthly Not. Roy. Astron. Soc. 493, 1292 (2020).

  22. H. Lammer, F. Selsis, I. Ribas, E. F. Guinan, S. J. Bauer, and W. W. Weiss, Astrophys. J. 598, L121 (2003).

  23. S. H. Gross, J. Atmospher. Sci. 29, 214 (1972).

  24. A. García Munoz, Planet. Space Sci. 55, 1426 (2007).

  25. K. J. Zahnle and J. C. G. Walker, Rev. Geophys. and Space Phys. 20, 280 (1982).

  26. A. Vidal-Madjar, A. Lecavelier des Etangs, J. M. Désert, G. E. Ballester, R. Ferlet, G. Hébrard, and M. Mayor, Astrophys. J. 676, L57 (2008).

  27. Д. В. Бисикало, П. В. Кайгородов, Д. Е. Ионов, В. И. Ше-матович, Астрон. журн. 90, 779 (2013).

  28. N. I. Shakura and R. A. Sunyaev, Astron. and Astrophys. 24, 337 (1973).

  29. M. R. Bate, S. H. Lubow, G. I. Ogilvie, and K. A. Miller, Monthly Not. Roy. Astron. Soc. 341, 213 (2003).

  30. T. T. Koskinen, R. V. Yelle, P. Lavvas, and N. K. Lewis, Astrophys. J. 723, 116 (2010).

  31. K. M. Flaherty, A. M. Hughes, S. C. Rose, J. B. Simon, et al., 843, id. 150 (2017).

  32. K. M. Flaherty, A. M. Hughes, R. Teague, J. B. Simon, S. M. Andrews, and D. J. Wilner, Astrophys. J. 856, id. 117 (2018).

  33. J. E. Pringle, Ann. Rev. Astron. Astrophys. 19, 137 (1981).

  34. А. С. Монин, А. М. Яглом, Статистическая гидромеханика. Механика турбулентности. Часть 1 (М.: Наука, 1965).

  35. А. С. Монин, А. М. Яглом, Статистическая гидромеханика. Механика турбулентности. Часть 2 (М.: Наука, 1967).

  36. J. C. B. Papaloizou and C. Terquem, Rep. Prog. Phys. 69(1), 119 (2006).

  37. G. Lesur and J. C. B. Papaloizou, Astron. and Astrophys. 513, id. A60 (2010).

  38. S. Kumar and C. S. Coleman, Monthly Not. Roy. Astron. Soc. 260, 323 (1993).

  39. Е. П. Курбатов, Д. В. Бисикало, П. В. Кайгородов, Успехи физ. наук 184, 851 (2014).

  40. Е. П. Курбатов, Д. В. Бисикало, Астрон. журн. 94, 477 (2017).

  41. Y. Nakao and S. Kato, Publ. Astron. Soc. Japan 47, 451 (1995).

  42. W. Rodi and D. B. Spalding, Waerme und Stoffuebertragung 3(2), 85 (1970).

  43. L. F. Shampine and M. W. Reichelt, Siam J. Sci. Compute 18, 1 (1997).

  44. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, et al., Nature Methods 17, 261 (2020).

  45. N. Shakura, Accretion Flows in Astrophysics (Springer, 2018).

  46. G. L. Withbroe, Astrophys. J. 325, 442 (1988).

  47. B. Jackson, R. Greenberg, and R. Barnes, Astrophys. J. 678, 1396 (2008).

  48. C. J. Clarke and J. Bouvier, Monthly Not. Roy. Astron. Soc. 319, 457 (2000).

  49. R. D. Alexander, C. J. Clarke, and J. E. Pringle, Monthly Not. Roy. Astron. Soc. 369, 216 (2006).

Дополнительные материалы отсутствуют.