Астрономический журнал, 2022, T. 99, № 4, стр. 325-333

1. ВВЕДЕНИЕ

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

Недавние наблюдения, проведенные на комплексе телескопов ALMA, показали, что распределения поверхностной плотности не соответствуют гладким степенн${\text{ы}}\prime $м законам, часто используемым теоретиками. Напротив, в масштабах от десятков до сотен а.е. повсеместно встречаются яркие кольца и темные промежутки [1–3]). На данный момент эти кольца и провалы чаще всего объясняются влиянием невидимых планет [4–8]. Однако было предложено множество альтернативных сценариев для объяснения этих особенностей. Одним из них является неустойчивость, связанная с затенением звездного излучения поверхностными неоднородностями диска.

Эта неустойчивость (в англоязычных изданиях известная как “Irradiation Instability” или “Thermal Wave Instability”) может приводить к возникновению поверхностных волн, бегущих по направлению к звезде, и кольцеобразных структур в околозвездных дисках [9–12]. Механизм развития данной физической неустойчивости следующий. Если на поверхности диска образуется небольшой выступ, то освещенная сторона выступа, обращенная к звезде, получает больше звездного света, чем невозмущенная поверхность, и сильнее нагревается. Нагретые элементы выступа прогревают своим излучением более нижние слои диска, в результате чего образовавшийся выступ приподнимается еще сильнее. Поскольку обращенная к звезде сторона выступа нагревается сильнее обратной стороны выступа, то возмущение также начинает перемещаться по направлению к звезде.

За последнее время достигнут заметный прогресс в изучении данной неустойчивости (см. [13, 14]. В частности, в работе [14] в рамках аналитического приближения показано, что эта неустойчивость действительно имеет место, причем ключевую роль в ее развитии играет детальный учет зависимости $H{\text{/}}h$ от расстояния (где $H$ – оптическая высота диска, $h$ – характерная гидростатическая шкала). Авторы статьи [14] помимо аналитической модели представили также и полуаналитическую модель, в которой для расчета эволюции экваториальной температуры использованы результаты расчета функции нагрева с помощью кода переноса излучения RADMC-3D.

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

2. КВАЗИГИДРОСТАТИЧЕСКАЯ 1+1D МОДЕЛЬ ПРОТОПЛАНЕТНОГО ДИСКА

Для расчета эволюции диска за основу взята модель из работ [15–17]. В данной модели решается задача переноса излучения, при этом учитывается нагрев звездным и межзвездным излучением, а также диффузия теплового (инфракрасного) излучения самого диска. Совместно с расчетом температуры происходит восстановление вертикальной структуры диска. Диффузия теплового излучения в модели рассчитывается только в вертикальном ($z$) направлении. Для ее моделирования решается система моментных уравнений переноса в эддингтоновском приближении:

где $T$ – температура среды, $E$ – плотность лучистой энергии, $z$ – вертикальная координата, $\rho $ – объемная плотность, ${{c}_{V}}$ – теплоемкость газопылевой среды, $c$ – скорость света, $a$ – постоянная Стефана, ${{\kappa }_{P}}$ [см$^{2}$/г] – коэффициент поглощения, усредненный по Планку (на ед. массы газопылевой среды), ${{\kappa }_{R}}$ [см$^{2}$/г] – коэффициент поглощения, усредненный по Росселанду, $S$ [эрг/с/г] – функция нагрева (на ед. массы газопылевой среды) звездным и межзвездным излучением, $S = S{\kern 1pt} *\; + {{S}_{{{\text{b}}g}}}$.

Интенсивность ультрафиолетового излучения, необходимая для вычисления функции нагрева звездным излучением $S{\kern 1pt} *$, вычисляется для каждой ячейки путем прямого интегрирования уравнения переноса излучения от звезды до рассматриваемого элемента среды по всему диску. Эта двумерная процедура отличается от метода, описанного в работе [15], где сделано предположение о постоянстве угла между направлением на звезду и поверхностью диска, и тем самым задача там сводилась к одномерной. В модифицированном методе при интегрировании уравнения переноса вдоль луча находятся все его пересечения с границами ячеек (см. рис. 1), что используется для точного вычисления полной оптической толщины.

Принципиальным моментом при вычислении функции нагрева звездным излучением в рамках нашей 1+1D-мерной модели является учет радиального градиента плотности внутри ячейки. В нашей модели функция нагрева $s{\kern 1pt} *$ [эрг см^‒3 с^‒1] звездным излучением вычисляется следующим образом:

где $L$ – светимость звезды, $\kappa $ [см²/г] – коэффициент поглощения звездного излучения, $R$ – радиальное расстояние от звезды до центра ячейки, $\tau $ – полная оптическая толщина на луче зрения до момента вхождения луча в ячейку, $\Delta \tau = \kappa \rho {\kern 1pt} *\Delta l$ – оптическая толщина самой ячейки вдоль луча, $\Delta l$   – длина отрезка луча внутри ячейки, $\rho {\kern 1pt} * = \frac{1}{4}\left( {{{\rho }_{L}} + 2{{\rho }_{i}} + {{\rho }_{r}}} \right)$ – усредненная плотность вдоль луча. При выводе формулы (3) из формального решения уравнения переноса излучения предполагалось, что плотность вдоль луча внутри ячейки линейно изменяется от ${{\rho }_{L}}$ до ${{\rho }_{i}}$ и от ${{\rho }_{i}}$ до ${{\rho }_{R}}$ (см. схему на рис. 2). Значения ${{\rho }_{L}}$ и ${{\rho }_{R}}$, в свою очередь, вычисляются с помощью линейной интерполяции плотности между центром текущей ячейки и центрами прилегающих слева и справа ячеек:

Таким образом, при вычислении функции нагрева вместо центральной плотности в ячейке (которая в процедуре расчета вертикальной структуры диска предполагается постоянной внутри ячейки) берется усредненное определенным образом по соседним ячейкам значение. При отсутствии такого (или более точного) усреднения ячейки с бóльшими $R$ не будут влиять на ячейки с меньшими $R$, и, как следствие, будет отсутствовать исследуемый нами механизм возникновения тепловой волны, бегущей снаружи внутрь. При учете интерполяции плотности потенциальный механизм возникновения волновой неустойчивости становится следующим. Рассмотрим два соседних столбца, один из которых назовем внутренним, другой – внешним. Пусть внешний столбец нагревается, тогда плотность ${{\rho }_{{i + 1}}}$ в его верхних слоях возрастает, т.к. возрастает характерная шкала высоты диска. Повышение плотности во внешних ячейках приводит к росту усредненной плотности $\rho {\kern 1pt} *$ в прилегающих внутренних ячейках, что ведет к увеличению $s{\kern 1pt} *$. Дальнейшее развитие процесса будет зависеть от характерных тепловых времен и геометрических особенностях перехвата излучения внутренним столбцом.

Нагрев межзвездным УФ-излучением рассчитывался нами по формуле:

где ${{\kappa }_{{{\text{uv}}}}}$ – коэффициент поглощения для межзвездного излучения, $W = {{10}^{{ - 14}}}$, ${{T}_{{{\text{b}}g}}} = {{10}^{4}}$ K – дилюция и температура межзвездного излучения, $\tau $ – оптическая толщина от поверхности диска до текущего элемента объема вдоль вертикального направления, $\Delta \tau $ – оптическая толщина текущего элемента, $\sigma $ – постоянная Стефана$ - $Больцмана. Нагрев межзвездным излучением включен в модель, поскольку он может быть сопоставим с нагревом от звезды во внешних частях диска в моменты затмений звезды на возникающих неоднородностях поверхности диска.

Решение системы уравнений (1), (2) находится с использованием полностью неявной схемы, принцип построения которой аналогичен численной схеме для решения квазилинейного уравнения теплопроводности с переменными коэффициентами, описанной в книге Н.Н. Калиткина “Численные методы” [18]. При линеаризации уравнений и использовании метода Ньютона итерации могут расходиться при достаточно высоком шаге по времени. Именно такая ситуация иногда (но достаточно редко) встречается и в нашем случае. В случае расходимости итерационного процесса Ньютона мы используем процедуру дробления шага по времени. Используемый метод позволяет рассчитывать тепловую эволюцию во всех участках диска, включая оптически толстые области, в которых характерные времена процессов нагрева и охлаждения сопоставимы с динамическими временами.

Метод моделирования тепловой структуры в данной модели тесно связан с методом восстановления вертикальной структуры диска в предположении локального гидростатического равновесия, которое находится из следующего уравнения:

где ${{k}_{{\text{B}}}}$ – постоянная Больцмана, ${{\mu }_{{\text{g}}}} = 2.3$ – средний молекулярный вес, ${{m}_{{\text{H}}}}$ – атомная единица массы, $G$ – гравитационная постоянная, ${{M}_{*}}$ – масса звезды. Отметим, что при восстановлении вертикальной структуры диска мы не учитываем самогравитацию диска, что оправдано для параметров диска, используемых в данной работе. Расчет вертикальной структуры диска позволяет получить полную информацию о распределении плотности и температуры в диске. Для решения уравнения гидростатического равновесия также используется неявный метод.

Принципиальным условием эффективности работы этих методов является оптимальный выбор пространственной сетки в $z$-направлении. Пространственная сетка должна отслеживать все заранее неизвестные особенности решения (градиенты плотности и температуры) с учетом существенного ограничения на количество ячеек и больших интервалов плотности газа (до 10 порядков). Нами разработан алгоритм построения и адаптивной модификации такой сетки, основанный на приближенном быстром решении уравнения гидростатического равновесия. Данный метод восстановления вертикальной структуры диска с расчетом переноса излучения был тщательно протестирован и сопоставлен с другими методами. В стационарном режиме распределения температуры хорошо согласуются с результатами моделирования структуры диска, полученными другими авторами. В нестационарном режиме характерные времена прихода к тепловому равновесию соответствуют аналитическим оценкам. Более подробное описание данного метода можно найти в статье [15].

В модели предполагается, что единственным источником непрозрачности является пыль, причем температуры газа и пыли равны. Отношение массы пыли к массе газа предполагается постоянным по всему диску и равным 0.01. Особенностью тепловой модели является использование усредненных по Планку и Росселанду непрозрачностей, зависящих от температуры. Эти коэффициенты взяты нами из работы [16], где они подробно описаны.

Основными входными параметрами модели являются масса и светимость звезды, которые мы предполагаем солнечными, а также распределение поверхностной плотности, которое выбирается нами в виде:

где ${{\Sigma }_{0}}$ – поверхностная плотность вблизи внутренней границы ${{R}_{{AU}}}$=1 а.е., $q = - 1$ – наклон в распределении плотности, ${{R}_{0}} = 3$ а.е., $p = 8$ – параметры сглаживания распределения плотности вблизи внутренней границы диска. Внутренняя и внешняя границы диска равны 1 и 100 а.е. соответственно. Отметим, что при отсутствии сглаживания распределения поверхностной плотности вблизи внутренней границы диска внутренние ячейки перехватывают большую часть излучения звезды и сильно влияют на структуру и эволюцию диска, затрудняя анализ непосредственно самой поверхностной неустойчивости. Начальное состояние диска рассчитывается нами в предположении о постоянном угле вхождения излучения в диск в рамках тепловой модели [15]. Шаг по времени выбран постоянным и равным 0.1 года, что меньше характерных тепловых времен при используемых параметрах модели. Расчет проводится на сетке 200 радиальных $ \times $ 120 вертикальных ячеек.

Отметим, что преимуществом этой численной модели перед используемой в работе [14] аналитической моделью является рассмотрение детальной двумерной структуры диска. При этом проводится расчет нагрева диска УФ-излучением звезды с помощью прямого интегрирования переноса излучения звезды, а также моделирование диффузии собственного теплового излучения в вертикальном направлении.

3. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ

Моделирование эволюции диска в рамках описанного 1+1D-приближения показывает, что в диске самопроизвольно возникают возмущения, бегущие по направлению к звезде. Характеристики этих возмущений сильно зависят, в частности, от заданной поверхностной плотности диска. Рассмотрим результаты моделирования для диска с ${{\Sigma }_{0}} = {{10}^{2}}$ г/см³ (соответствующая масса диска $7 \times {{10}^{{ - 3}}}{{M}_{ \odot }}$). На рис. 3 приведена структура диска (распределения плотности, температуры, функции нагрева УФ-излучением и плотности энергии ИК-излучения) на момент времени 10 лет после начала эволюции диска. В окрестности 35 а.е. видно возмущение в форме излома в распределении плотности, за этим изломом (фронтом) происходит увеличение характерной высоты диска. Развертка результатов моделирования по времени показывает, что данное возмущение распространяется снаружи внутрь. За фронтом отчетливо видно повышение функции нагрева звездным излучением в верхних слоях атмосферы, что связано с менее острым углом вхождения звездного излучения в диск. Максимуму функции нагрева в окрестности возмущения соответствует и максимум температуры. Температура также повышена в экваториальных слоях диска на некотором отдалении от излома. Это связано с конечным временем прогрева внутренних слоев ИК-излучением поверхностных слоев. В окрестности возмущенного слоя повышена плотность энергии ИК-излучения, что также свидетельствует о большем прогреве данного слоя.

Временн${\text{у}}\prime $ю эволюцию диска удобно анализировать с помощью двумерной диаграммы c зависимостью распределения экваториальной температуры от времени, изображенной на левой панели рис. 4. Бегущие внутрь волны хорошо на ней выделяются, им соответствуют максимумы, перемещающиеся справа налево с увеличением времени. Волны распространяются от внешней границы диска до области $ \approx $3 a.e., которая прозрачна к УФ-излучению звезды (из-за малой плотности за счет сглаживания распределения $\Sigma (R)$). На данной диаграмме также построено распределение характерного теплового времени, вычисленного по формуле:

где ${{\tau }_{{{\text{IR}}}}} = {{\kappa }_{P}}\Sigma $ – оптическая толщина в вертикальном направлении для собственного теплового излучения, ${{T}_{{\text{m}}}}$ – экваториальная температура, ${{c}_{p}}$ – удельная теплоемкость, $\sigma $ – постоянная Стефана$ - $Больцмана. Формула (9) получена в предположении, что диск является оптически толстым к собственному тепловому излучению, и служит оценкой для характерного времени нагрева (или охлаждения) экваториальных слоев до температуры ${{T}_{{\text{m}}}}$. Периодичность возникновения волн, $ \approx {\kern 1pt} 3$ года, близка к максимальному характерному тепловому времени $ \approx {\kern 1pt} 7$ лет, т.е. возмущенные области диска действительно успевают прогреваться до экваториальных слоев.

На правой панели рис. 4 показана эволюция экваториальной температуры для модели диска с теми же параметрами, но в которой отсутствует интерполяция плотности в радиальном направлении при вычислении темпа нагрева звездным излучением, т.е. где принято $\rho {\kern 1pt} * = {{\rho }_{i}}$, см. ф-лу (3). В данном случае в диске возникают несколько горбов, перехватывающих излучение звезды и затеняющих расположенные за ними области. Медленное движение этих горбов по направлению к звезде на начальных этапах эволюции связано с процессами установления самосогласованного решения на интервале характерных тепловых времен диска. Регулярных поверхностных волн в такой модели не возникает. Этот расчет показывает важность использования высоких порядков аппроксимации при интегрировании уравнения переноса излучения в данной задаче.

На рис. 5 приведены распределения физических величин для модели с увеличенной на порядок поверхностной плотностью, ${{\Sigma }_{0}} = {{10}^{3}}$ г/см³, на момент времени 10 лет от начала эволюции. Структура данного диска неоднородная, в распределении плотности энергии ИК-излучения можно выделить два возмущения: в окрестности 10 и 40 а.е. Морфология данных возмущений в целом повторяет картину, описанную для предыдущей модели. Однако в отличие от предыдущей модели область прогрева в окрестности возмущений не достигает экваториальной плоскости. Это хорошо видно из левой верхней панели рис. 6, где показана эволюция распределения экваториальной температуры в первые 20 лет. На данной диаграмме волны не проявляются, за исключением слабых осцилляций в окрестности 3 а.е. Бегущие внутрь волны, однако, хорошо просматриваются на распределениях поверхностной температуры (левая нижняя панель рис. 6) и температуры на половине поверхностной плотности до экватора (правая верхняя панель рис. 6). На рассматриваемом временнóм интервале возмущения зарождаются в окрестности 20 а.е. и распространяются внутрь примерно за 3 года. По мере приближения к внутренней границе диска скорость распространения возмущений уменьшается. Период прохождения этих волн существенно короче характерного теплового времени ${{t}_{{{\text{therm}}}}}$, рассчитанного для всей толщи диска. Это видно из правой нижней панели рис. 6, где построено распределение характерного теплового времени для всей толщи диска. Видно, что время ${{t}_{{{\text{therm}}}}}$ составляет сотни лет для области между 3 и 20 а.е. На той же панели приведена эволюция экваториальной температуры диска в течение 1000 лет, т.е. на временах, сопоставимых с характерным тепловым временем диска. На данной диаграмме волны внутри 20 а.е. также не отождествляются, однако видны квазипериодические возмущения с периодом возникновения $ \approx {\kern 1pt} 100$ лет, распространяющиеся внутрь с внешней границы диска. Очевидно, что волны во внешней области диска отождествляются на распределении экваториальной температуры в связи с тем, что характерные тепловые времена в этой области существенно меньше, чем в более внутренних частях, т.е. диск, будучи возмущенным на поверхности в этих областях, успевает прогреться вплоть до экватора.

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

Основной целью данной работы являлась проверка возможности самопроизвольного возбуждения поверхностных тепловых волн в протопланетных дисках. Для достижения этой цели мы модифицировали разработанную нами ранее 1+1D модель диска, включив в нее более точный (двумерный) расчет нагрева диска излучением звезды. В модифицированной модели нагрев рассчитывается с помощью интегрирования уравнения переноса вдоль всего диска и учитывается радиальный градиент плотности внутри ячейки, для которой вычисляется функция нагрева. Это позволяет учесть геометрические эффекты, имеющие принципиальное значение для данной задачи. Результаты моделирования показали, что в диске действительно развивается неустойчивость, приводящая к появлению возмущений, распространяющихся снаружи внутрь. Тем самым мы подтверждаем вывод о возможности самопроизвольного возбуждения поверхностных тепловых волн в протопланетных дисках, полученный ранее другими авторами.

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

Безусловно, при температурах порядка тысяч градусов предположение нашей модели о том, что непрозрачность вещества обусловлена только пылью, не выполняется. При таких температурах становятся важным испарение пыли, т.е. основной вклад в непрозрачность начинает вносить газ. При высоких температурах значительными становятся и другие процессы, такие как диссоциация водорода. При этом можно ожидать появление в диске различных неустойчивостей и связанной с ними сложной динамики. Это отдельная тема, требующая детального изучения. Мы можем заметить, что при используемых параметрах модели такие высокие температуры получаются только в самых внутренних ($ < {\kern 1pt} 3$ а.е.) областях диска и не затрагивают более отдаленные слои, где зарождаются и начинают распространяться исследуемые нами волны.

В своей модели мы также не учитывали процессы осаждения, миграции, роста и разрушения пыли, которые играют важнейшую роль в протопланетных дисках. Их учет может существенно модифицировать или даже полностью подавить исследуемую нами неустойчивость диска, связанную с эффектами самозатенения. Отметим также, что само по себе наличие пыли приводит к своим неустойчивостям, таким как потоковая неустойчивость [19]. Потоковая неустойчивость сейчас становится очень популярной для объяснения процессов турбулизации и формирования зародышей планет. Изучение этих процессов – отдельное направление, для чего разрабатываются сложные динамические модели. В своей работе мы не рассматривали все эти процессы, отложив изучение их взаимодействия с исследуемой нами неустойчивостью на будущее.

Существенным отличием наших результатов от выводов аналитической модели [14] являются полученные нами характерные временны́е масштабы распространения возмущений. В рамках аналитической модели [14] характерное время распространения возмущений соответствует характерному тепловому времени для всей толщи диска. Этот результат закономерен, поскольку в аналитической модели высота диска определяется именно экваториальной температурой. В нашей модели времена распространения для оптически толстых дисков оказываются существенно ниже времен ${{t}_{{{\text{therm}}}}}$. Такая ситуация объясняется тем, что бегущие возмущения не затрагивают всю толщину диска, т.е. механизм возбуждения волн может работать в приповерхностных слоях.

Полученные результаты свидетельствуют о необходимости изучения данной неустойчивости в рамках более согласованной модели. Действительно, в основе использованного нами 1+1D подхода лежит несколько ключевых приближений, которые могут существенно искажать реальную картину. Такими приближениями являются: 1) гидростатическое равновесие в вертикальном направлении; 2) отсутствие диффузии теплового излучения в радиальном направлении; 3) отсутствие динамических эффектов, в том числе отсутствие динамической связи между областями диска в радиальном направлении. Влияние этих эффектов на возбуждение поверхностных тепловых волн необходимо исследовать в рамках двумерной или трехмерной гидродинамической модели с полноценным расчетом переноса излучения. Особое внимание в такой модели необходимо уделить деталям расчета функции нагрева звездным излучением.