Астрономический журнал, 2021, T. 98, № 3, стр. 197-211

1. ВВЕДЕНИЕ

Образование звезд и планет – одна из ключевых проблем астрофизики и одна из наиболее бурно развивающихся областей наблюдательной астрономии. Звезды образуются в результате гравитационного коллапса компактных газопылевых дозвездных конденсаций. Существенная часть вещества дозвездной конденсации, прежде чем попасть на протозвезду, образует околозвездный диск благодаря сохранению углового момента аккрецирующего вещества. В этом диске могут образовываться планеты в результате слипания пылевых частиц и их последующего гравитационного объединения [1], либо в результате гравитационной неустойчивости самого газового диска [2, 3], либо в результате комбинации процессов формирования сгустков, их миграции и дегазации [4, 5].

Несмотря на успехи построения общей картины образования звезд и планет, детали эволюции протопланетного диска и роль различных физических процессов до сих пор окончательно не выяснены (см. обзор [6]). Во многом это связано с трудностью наблюдений протопланетных дисков (ПД). На протяжении многих лет одним из основных источников информации о протопланетных дисках являлся анализ спектрального распределения энергии излучения всего диска (см., напр., [7, 8]). Именно особенности спектра (так называемые инфракрасные избытки) легли в основу классификации молодых звездных объектов [9]. Относительно простые аналитические модели позволяют восстановить параметры дисков (массу, размер, профиль распределения плотности) из его спектра [10, 11]. Однако восстановление информации о структуре диска из его интегрального спектра – сложная обратная задача, решение которой может быть неоднозначным [12].

С появлением пространственно разрешенных изображений дисков в различных длинах волн стало возможным изучать структуру диска непосредственно. В то время как оптические и инфракрасные изображения дисков, полученные, например, с помощью телескопов Hubble [13] и VLT [14], проявляют поверхностные слои диска, радио-интерферометрические наблюдения предоставляют информацию о внутренней структуре диска (см., напр., [15]). По-настоящему революционный шаг в наблюдательном изучении ПД произошел после начала работы интерферометра ALMA (см. раздел 4), полученные в рамках проекта DSHARP [16]. В целом изображения ПД в различных частотных диапазонах продемонстрировали разнообразные морфологические особенности дисков (кольца, спирали, полости), сильные радиальные и вертикальные градиенты обилий молекул, связанные с изменением физических условий, градиенты спектрального индекса, обусловленные, вероятно, эволюцией и миграцией пыли и т.д. Все это свидетельствует об ограниченности простых моделей (в основе которых “стандартные” предположения о монотонном распределении поверхностной плотности, постоянном отношении массы пыли к массе газа и др.) для интерпретации наблюдений. Вместе с тем новые наблюдательные данные предоставляют почву для разработки детальных теоретических моделей эволюции дисков и исследования ключевых процессов в них.

В направлении разработки детальных теоретических моделей дисков и интерпретации наблюдений на их основе уже сделаны впечатляющие шаги. Например, в работе [17] трехмерная МГД-модель протопланетного диска совместно с расчетом переноса излучения была использована для визуализации дисков, эволюция которых управляется магнито-вращательной неустойчивостью. В статье [18] трехмерная гидродинамическая модель была совмещена с расчетом переноса излучения для изучения наблюдательных проявлений щелей, расчищаемых планетами в околозвездных дисках. Объединение термохимической модели диска с переносом излучения, а также детальный обзор о самосогласованном моделировании структуры и изображений дисков можно найти в статье [19].

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

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

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

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

Блок-схема, описывающая основные этапы нашего моделирования, приведена на рис. 1. В следующих разделах мы последовательно опишем эти этапы и обсудим полученные результаты.

Рис. 1.

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

2. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ ДИНАМИЧЕСКОЙ ЭВОЛЮЦИИ ДИСКА

Динамическая эволюция диска моделировалась с помощью программного комплекса F-EOSAD, детальное описание которого представлено в работе [20]. Этот код является дальнейшим развитием двумерной (r, φ) модели формирования и долговременной эволюции околозвездного газопылевого диска [21]. В рамках данной модели, в частности, удается воспроизвести режим эпизодической аккреции и объяснить вспышки фуоров падением гравитационно-связанных фрагментов, формирующихся в аккреционном диске и мигрирующих к звезде. Модель FEOSAD является одной из самых детальных в мире для описания самосогласованной динамики пыли и газа с учетом большого числа ключевых физических процессов. Пыль в модели представлена в виде двух компонентов – мелкого и крупного. Мелкая пыль (с размерами менее 10^–4 см) сцеплена с газом, в то время как крупная пыль (с размерами более 10^–4 см) может испытывать дрейф относительно газа. При расчете дрейфа крупной пыли учитывается обмен импульсом с газом и обеспечивается корректный расчет миграции пыли в широком диапазоне чисел Стокса. Содержание крупной пыли в модели может увеличиваться за счет столкновений с мелкой, при этом максимальный размер крупной пыли изменяется. Рост пыли в модели ограничен так называемым фрагментационным пределом, при котором столкновения пылинок приводят не к их слипанию, а к разрушению. С помощью динамической модели FEOSAD получены результаты об эволюции пыли в протопланетном диске [22], исследовано влияние внутренней области диска на его эволюцию [23], исследована динамика гальки [24].

В наших расчетах основные параметры соответствуют базовой модели из статьи [20], за исключением тех, что перечислены ниже. Масса начального молекулярного облака выбрана равной 0.5 ${{M}_{ \odot }}$, значение фрагментационной скорости пыли уменьшено до 3 м/c. В качестве внутреннего граничного условия мы использовали модель “умной” ячейки из работы [23] с эффективностью аккреции $\xi = 0.5$. Мы также использовали обновленные коэффициенты поглощения и рассеяния для пылевого населения (см. подробнее раздел 3.2). При использованных параметрах диск формируется через 30 тыс. лет после начала коллапса облака.

Для исследования были выбраны три момента времени. На рис. 2 представлены распределения основных физических параметров диска (поверхностной плотности газа, выросшей и мелкой пыли, экваториальной температуры и максимального размера выросшей пыли) в каждый из выбранных моментов времени. Левый столбец соответствует моменту времени 81 тыс. лет, на котором проявляется молодой, фрагментирующий гравитационно-неустойчивый диск с двумя ярко выраженными спиралями. В среднем столбце рис. 2 изображена система на момент 83 тыс. лет. В данный момент происходит вспышка аккреции вследствие развития магнито-ротационной неустойчивости (МРН). В используемой гидродинамической модели вспышка МРН активизируется внутри центральной поглощающей ячейки при достижении в этой области порогового значения температуры 1500 K, т.е. температуры, достаточной для термической ионизации среды. Наконец, проэволюционировавший менее плотный осесимметричный диск, расширившийся в процессе эволюции, изображен в правом столбце. Возраст данного диска составляет 1 млн. лет. Регулярная структура и отсутствие экстремальных условий в таком диске сопоставимы с простейшими моделями, используемыми при интерпретации наблюдений.

Рис. 2.

Результаты моделирования динамической эволюции газопылевого диска с помощью гидродинамического кода FEOSAD. Распределения основных параметров диска в области 300 × 300 а.е. для различных времен эволюции: 81 тыс. лет, 83 тыс. лет и 1 млн. лет в левом, среднем и правом столбцах панелей соответственно. Верхняя строка соответствует распределению поверхностной плотности газа, 2-я строка – поверхностной плотности выросшей пыли, 3-я строка – поверхностной плотности мелкой пыли, 4-я строка – температуры в экваториальной плоскости, нижняя строка – распределению максимальных размеров пылинок. Все величины представлены в логарифмической шкале цветом.

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

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

Отметим, что в рамках данной модели основной рост пыли происходит на ранних стадиях эволюции. Таким образом, на моменты 81 и 83 тыс. лет максимальный размер пылинок достигает $ \sim 1$ мм (нижний ряд рис. 2). Впоследствии максимальный размер пыли определяется процессами фрагментации и миграции, что ведет к уменьшению размеров пыли со временем. Для модели 1 млн. лет максимальный размер падает до 0.3 мм. Интересной особенностью является пониженный размер пылевых частиц в центральных областях диска. Это связано с повышенной температурой там, а значит и пониженным фрагментационным барьером.

Рассмотрим более подробно тепловую структуру диска вдоль одного из радиальных направлений в диске для моментов времени 81 тыс. лет и 1 млн. лет. На левой панели рис. 3 показаны темпы нагрева поверхности диска УФ-излучением звезды ${{\Gamma }_{{{\text{UV}}}}}$, экваториальной плоскости ИК-излучением самого диска ${{\Gamma }_{{{\text{IR}}}}}$, а также темп нагрева вязкой диссипацией ${{\Gamma }_{\nu }}$. Эти темпы вычисляются следующим образом:

Рис. 3.

Радиальная тепловая структура газопылевого диска по результатам расчета динамической модели на момент 81 тыс. лет (верхняя группа панелей) и 1 млн. лет (нижний ряд). Радиальный срез соответствует углу $\varphi = \tfrac{{5\pi }}{4}$, отсчитанному от положительного направления оси X на рис. 2 против часовой стрелки. Левый столбец: распределения по радиусу вязкого нагрева (оранжевая кривая), нагрева в экваториальной плоскости ИК-излучением (синяя кривая) и нагрева поверхности УФ-излучением звезды (фиолетовая кривая). Серым и коричневым во втором столбце показаны зависимости вертикальных оптических толщин, рассчитанных по Планку и Росселанду. Третий столбец: фактическая температура в экваториальной плоскости диска (красная кривая), равновесная температура, вычисленная из условия равенства вязкого нагрева и охлаждения (оранжевая кривая), равновесная температура, вычисленная из условия равенства ИК-нагрева и охлаждения (синяя кривая). На правых панелях показаны характерные временные шкалы различных процессов нагрева: ИК-излучением (синяя кривая), вязкостью (оранжевая кривая), адиабатическим сжатием (зеленая кривая). Голубым цветом обозначено характерное время охлаждения. Вертикальными пунктирными линиями обозначено положение спирали на данном азимутальном срезе.

Для модели с возрастом 81 тыс. лет вязкий нагрев доминирует над нагревом излучением в экваториальной плоскости диска внутри 20 а.е., несмотря на то, что поток звездного излучения на поверхности диска превышает темп вязкой диссипации только в области r < 5 а.е. Такая ситуация реализуется благодаря высоким значениям оптической толщины диска к тепловому излучению внутри 20 а.е. – оптическая толщина здесь составляет сотни единиц (см. вторую слева панель на рис. 3). Выделившаяся в результате вязкой диссипации энергия задерживается в диске, повышая его температуру, тогда как высокая оптическая толщина уменьшает скорость нагрева экваториальных слоев диска поверхностными. Как следствие, экваториальная температура диска внутри 20 а.е. определяется в основном вязким нагревом, а снаружи 20 а.е. – излучением звезды (см. третью панель на рис. 3).

Выбранный радиальный срез пересекает спиральную волну плотности в окрестности 40 а.е., где наблюдается небольшой скачок экваториальной температуры. Как видно из темпа вязкой диссипации, распределения оптической толщины и профиля равновесной температуры (вычисленной из условия равенства вязкого нагрева и темпа охлаждения), повышение экваториальной температуры здесь по крайней мере частично связано с вязким нагревом в условиях большой оптической толщины. Стоит отметить, что некоторый вклад в повышенную температуру внутри спирали может вносить также нагрев, связанный с адиабатическим сжатием этой области. Для демонстрации этого на правой панели рис. 3 приведены характерные времена процессов нагрева–охлаждения, а также характерное динамическое время (время распространения звука в вертикальном направлении) в зависимости от расстояния до звезды. Характерные времена нагрева и охлаждения рассчитаны как $t = \tfrac{E}{{2G}}$, где E – полная тепловая энергия, полученная из данных о температуре и плотности газа, а G – один из источников нагрева, описанных уравнениями (2), (3) или охлаждения за счет высвечивания ИК-излучения (см. [26], формула (5)):

В модели с возрастом 1 млн. лет оптическая толщина и температура внутри 40 а.е. уменьшились примерно на порядок величины по сравнению с моделью для 81 тыс. лет. Как видно из левой панели рис. 3, вязкий нагрев сопоставим с нагревом ИК-излучением от поверхностных слоев только внутри 3 а.е. Отметим также, что нагрев ИК-излучением от поверхностных слоев диска превышает поток звездного излучения в области 30–100 а.е., что связано с дополнительным вкладом фонового излучения при расчете ${{\Gamma }_{{{\text{IR}}}}}$ в модели FEOSAD.

3. ФОРМИРОВАНИЕ МОДЕЛИ ДЛЯ РАСЧЕТА ПЕРЕНОСА ИЗЛУЧЕНИЯ

Для моделирования переноса излучения и построения синтетических изображений околозвездного диска мы используем трехмерный код RADMC-3D, разработанный К. Дуллемондом (С. Dullemond). Этот код основан на методе Монте-Карло и позволяет рассчитывать перенос излучения на пыли с учетом процессов поглощения, рассеяния и теплового переизлучения. Этот комплекс активно используется для моделирования протопланетных дисков, околозвездных оболочек, молекулярных облаков и доступен для свободного использования¹¹.

3.2. Формирование коэффициентов поглощения и рассеяния

В исходной модели FEOSAD [20] использовались усредненные по частоте непрозрачности из статьи [28]. Однако эти непрозрачности рассчитаны для фиксированного распределения пыли по размерам. В то же время пыль в модели F-EOSAD эволюционирует – в каждой расчетной ячейке может быть свой максимальный размер пылинок. Чтобы привести в соответствие динамическую и тепловую модели пыли, мы перешли к новым коэффициентам непрозрачности. Для этого с помощью теории Ми вначале были рассчитаны спектральные коэффициенты поглощения и рассеяния для разных значений максимального размера пылинок. При этом распределение пылинок по размерам бралось степенны́м, $n(a) \propto {{a}^{{ - 3.5}}}$, c   минимальным радиусом пылинок ${{a}_{{min}}} = 5\; \times $ × 10^‒7 см в соответствии с динамической моделью пыли. Пылинки предполагались полностью силикатными. Зависящие от частоты и максимального размера пылинок коэффициенты поглощения и рассеяния приведены на рис. 4. Отметим, что полученные спектральные коэффициенты поглощения существенно зависят от максимального размера пыли, что подчеркивает необходимость их введения.

Рис. 4.

Верхний ряд: зависимости усредненных по Планку (слева) и Росселанду (справа) непрозрачностей от температуры и размеров пылинок (используется в динамической модели FEOSAD). Нижний ряд: зависимости коэффициентов поглощения (слева) и рассеяния (справа) от частоты и размеров пылинок (используется для расчета переноса излучения с помощью RADMC-3D).

Далее на основе полученных спектральных коэффициентов поглощения и рассеяния были вычислены усредненные по Планку и по Росселанду непрозрачности. Зависимости планковской ${{\kappa }_{P}}(T$, a_max) и росселандовской непрозрачностей от температуры и максимального размера пылинок также приведены на рис. 4. Как и следовало ожидать, двумерная морфология распределения ${{\kappa }_{P}}(T,{{a}_{{max}}})$ повторяет морфологию распределения ${{\kappa }_{\nu }}({{a}_{{max}}})$. При формировании усредненных коэффициентов непрозрачности мы, в отличие от работы [28], не учитывали возможность испарения пылинок при высоких температурах и не использовали газовые непрозрачности.

В то время как при моделировании динамической эволюции диска используются усредненные по частоте непрозрачности (зависящие от температуры), для моделирования переноса излучения с помощью RADMC-3D необходимо использовать исходные спектральные коэффициенты поглощения и рассеяния. В каждой ячейке расчетной сетки для RADMC-3D в зависимости от максимального размера пылинок в ней задаются свои распределения ${{\kappa }_{\nu }}$ и ${{\sigma }_{\nu }}$. Таким образом, как для моделировании эволюции диска, так и для моделирования переноса излучения мы используем согласованные модели пыли.

4. РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ ПЕРЕНОСА ИЗЛУЧЕНИЯ

Наше моделирование показало, что внедрение внутреннего источника нагрева при расчете температуры диска с помощью RADMC-3D значительно увеличивает расчетное время (от нескольких десятков до сотен раз по сравнению с нулевым внутренним источником). Особенно сильным это увеличение оказалось для моделей диска на ранней фазе эволюции диска, составив до нескольких месяцев непрерывного счета, что является практически неприемлемым. В частности, с нулевым внутренним источником расчет на 72-ядерной вычислительной станции составляет 24–36 ч. В случае включения внутреннего источника длительность расчета при том же количестве фотонов увеличивается до 150–200 дней. Очевидно, это связано с высокими оптическими толщинами диска в этих моделях. В результате нам удалось промоделировать тепловую структуру диска с полным набором фотонов $N = {{10}^{9}}$ без учета вязкого нагрева для всех моделей, а также с учетом вязкого нагрева для модели с возрастом 1 млн. лет (в силу меньшей плотности проэволюционировавшего диска и, как следствие, меньших оптических толщин). Для моделей на ранней стадии эволюции были проведены расчеты с сокращенным набором фотонов $N = {{10}^{7}}$. Это позволило сделать выводы о влиянии вязкого нагрева на распределение температуры. Однако при малом числе фотонов флуктуации температуры сопоставимы со значениями самой температуры, поэтому для построения изображений и спектральных распределений энергии мы будем использовать модели без вязкого нагрева. Далее мы последовательно проанализируем результаты моделирования температуры и синтетических изображений.

4.1. Тепловая структура

На рис. 5, 7 представлены результаты расчета тепловой структуры с помощью RADMC-3D для выбранных моделей диска. Для модели без вязкого нагрева с возрастом 81 тыс. лет экваториальная температура во всем диске систематически ниже той, что получена в FEOSAD. В расчете RADMC-3D она достигает своих минимальных значений в окрестности 30 а.е., а затем происходит некоторое увеличение по направлению к внешним областям. Причиной небольшого увеличения экваториальной температуры во внешней области диска служит то, что эта область становится оптически тонкой к УФ-излучению, в результате прямое и рассеянное УФ-излучение непосредственно достигает экваториальной плоскости диска и определяет температуру в ней. При учете вязкого нагрева в RADMC-3D температура внутри 30 а.е. существенно поднимается, приближаясь к результатам гидродинамической модели. За пределами 30 а.е. на распределении температуры в модели с вязким нагревом виден сильный шум, связанный с недостаточным количеством фотонов. Самым существенным отличием распределений, полученных в RADMC-3D, от динамической модели является то, что спирали в них имеют более низкую температуру, чем внешние области диска, в то время как в динамической модели температура внутри спиралей повышена относительно фона (см. рис. 5). Возможной причиной этому является то, что при моделировании переноса излучения с помощью RADMC-3D мы не учли нагрев, связанный с работой сил давления (адиабатический нагрев). Как показал анализ в разделе 2, адиабатический нагрев может быть важным фактором. Другой причиной расхождения результатов может быть то, что учет трехмерной структуры диска и использование непрозрачностей, зависящих от частоты в RADMC-3D, существенно упрощают выход теплового излучения из спиралей, тем самым приводя к меньшей равновесной температуре в них.

Рис. 5.

Результаты моделирования тепловой структуры молодого спокойного диска на 81 тыс. лет (левая панель), диска со вспышкой на 83 тыс. лет (центральная панель) и проэволюционировавшего диска на 1 млн. лет (правая панель). Каждая панель состоит из распределений плотности (верхняя строка), температуры FEOSAD (вторая строка), температуры RADMC-3D с нулевым внутренним источником (третья строка), температуры RADMC-3D с включенным внутренним нагревом (нижняя строка). В левом столбце панели показаны экваториальные распределения, в правом столбце – меридиональный срез вдоль угла $\varphi = \tfrac{\pi }{4}$, отсчитываемого от положительного направления оси Х.

Рис. 6.

Распределения экваториальной температуры в исходной динамической модели (синие профили) и после расчета с помощью комплекса RADMC-3D: оранжевые профили соответствуют модели без вязкого нагрева, зеленые – с вязким нагревом. Левая, центральная и правая панели соответствуют возрастам диска 81 тыс. лет, 83 тыс. лет и 1 млн. лет.

Рис. 7.

Синтетические изображения газопылевого диска на длинах волн 1.5 мкм (левая пара панелей), 300 мкм (центральная пара) и 1.3 мм (правая пара) для времени 81 тыс. лет. Каждая пара представляет собой вид с углом наклона к центральной оси 0° (слева) и 60° (справа). Данные об интенсивности приведены цветом в логарифмической шкале [Ян/пиксель]. Область, охватываемая изображением, соответствует линейному масштабу 150 × 150 а.е.

В модели со вспышкой (83 тыс. лет) распределения экваториальной температуры в динамической модели и по результатам RADMC-3D отличаются не так сильно. В этом случае определяющим источником нагрева во всем диске является нагрев звездным излучением, поэтому различия, связанные с внутренними источниками нагрева и межзвездным нагревом, проявляются слабее. На распределениях экваториальной температуры, полученных и в FEOSAD, и в RADMC-3D, видна пониженная температура внутри спиралей, однако в результатах RADMC-3D холодные спирали более выражены. Включение внутреннего нагрева не приводит к значительным изменениям в общей картине. Температурное распределение повторяет распределение в случае с нулевым внутренним источником, за исключением некоторых пиковых отклонений, которые являются следствием шума при малом количестве фотонов.

Для проэволюционировавшего диска (1 млн. лет) экваториальная температура по результатам RADMC-3D систематически ниже той, что в исходной динамической модели, как и в случае модели на момент 81 тыс. лет. Учет вязкого нагрева в данной модели выполнен с задействованием полного количества фотонов, а значит и наиболее согласован по сравнению с другими моделями. Как видно, включение внутреннего источника немного увеличивает температуру внутри 20 а.е., что указывает на необходимость рассмотрения этого процесса даже на этой фазе эволюции диска. Однако даже с учетом вязкого нагрева температура по результатам RADMC-3D систематически ниже динамической. Похожие различия в распределениях экваториальной температуры были также отмечены в работе [26] в рамках (2+1)-мерной модели с нестационарным методом переноса излучения в вертикальном направлении. Мы считаем, что систематически более высокая во всем диске температура в гидродинамической модели связана с приближенным характером функций нагрева и охлаждения, используемых в FEOSAD (см. [26], формулы (5), (6)).

5. ВОССТАНОВЛЕНИЕ ПАРАМЕТРОВ ДИСКА НА ОСНОВЕ РАССЧИТАННЫХ СИНТЕТИЧЕСКИХ ИЗОБРАЖЕНИЙ

Для оценки массы протопланетных дисков широко используется связь массы диска с потоком излучения в (суб)мм-диапазоне [30, формула (2)]:

Более детальный метод восстановления параметров диска может быть основан на анализе спектрального распределения энергии (SED) от всего диска. Для этого мы использовали программный комплекс для подгонки спектров протопланетных дисков “SED-Fitter” [31]. Этот комплекс используется для анализа наблюдательных данных, (см., напр., [32]). В основе комплекса лежит феноменологическая модель диска с 14 свободными параметрами, среди которых характеристики звезды (масса, радиус, температура), параметры оболочки (темп падения вещества из оболочки на диск, внешний радиус) и параметры самого диска (масса, границы, темп аккреции, геометрические характеристики). Программа осуществляет поиск наиболее близких спектров на основе предрассчитанной сетки моделей, насчитывающей ~20 тыс. моделей и охватывающей широкий набор звездных объектов (массой от 0.1 до 50 ${{M}_{ \odot }}$) на различных этапах их эволюции – от самых ранних (с возрастом порядка нескольких тысяч лет) до 10 млн. лет. Для каждой из моделей введено по 10 углов наклона к оси наблюдения. Таким образом, общее число спектров соответствует ~200 тыс. В наших расчетах мы не исследовали влияние угла наклона на точность восстановления параметров и поэтому предполагаем, что ось диска направлена к наблюдателю. Отметим, что на ранних фазах эволюции диск должен быть окружен протяженной оболочкой, которая затеняет диск даже в полярном направлении. При моделировании спектрального распределения энергии с помощью RADMC-3D этот эффект нами не учитывался.

Подгонка спектрального распределения энергии производится на определенных длинах волн. Для этого в SED-Fitter используется система фильтров. Фильтры могут быть аналогичны используемым в телескопах при реальных наблюдениях, а также монохроматическими. В работе задействован набор монохроматических фильтров для длин волн, согласующихся с набором частот, использованных в расчетах RADMC-3D. Тем не менее такое согласование несколько сокращает количество точек, участвующих в подгонке, поскольку наборы длин волн, использованные в RADMC-3D и доступные в базе моделей SED-Fitter, различаются. В итоге из 150 частот, использованных в расчете переноса излучения и модельных диаграмм, при подгонке были использованы 70.

В ходе подгонки алгоритмом SED-Fitter было отобрано от 17 до 48 (в зависимости от модели) близких спектральных распределений энергии для каждого исходного спектра. Результаты подгонки представлены на рис. 8. Как и следовало ожидать, поток излучения от диска во время вспышки ощутимо выше во всем спектре. В модели проэволюционировавшего диска (t = 1 млн. лет), в среднем, поток ниже, чем для молодого спокойного диска; в нем также заметно некоторое уплощение максимума в сторону длинноволновой области по сравнению с молодым диском. Последнее связано с увеличением размеров диска: диск перехватывает больше излучения от звезды, однако, температура внешних областей меньше, что выражается в появлении дополнительного низкочастотного излучения.

Рис. 8.

Синтетические спектральные распределения энергии и результаты их подгонки с помощью комплекса “SED-fitter”. На левой панели отображены данные, соответствующие модели с возрастом 81 тыс. лет, на центральной – 83 тыс. лет, на правой – 1 млн. лет. Оранжевые точки соответствуют рассчитанному в RADMC-3D распределению энергии, серыми профилями показаны близкие распределения имеющихся в базе данных моделей. Сплошная черная линия – наиболее близкое распределение из имеющихся в базе. Максимальный уровень потока в модели 81 тыс. лет обозначен штриховой синей линией.

По результатам подгонки были определены все 14 параметров системы, среди которых мы выделили шесть наиболее интересующих нас: масса диска, темп аккреции на звезду, внешняя и внутренняя границы диска, полная светимость звезды и возраст системы. Восстановленные и исходные параметры системы представлены на рис. 9. Видно, что массы дисков, полученные в ходе подгонки SED-Fitter, систематически ниже, чем исходные значения, особенно в моделях без вспышки (t = 81 тыс. лет и t = 1 млн. лет). Также в моделях без вспышки значения темпов аккреции на несколько порядков величин выше исходных. В то же время для модели со вспышкой (t = 83 тыс. лет) SED-Fitter хорошо воспроизводит исходные массу диска и темп аккреции. Систематическую недооценку масс дисков мы получили не только при подгонке спектров, но и при расчете масс с помощью выражения (11). Результаты расчетов масс по формуле (11) нанесены на верхний ряд панелей зелеными штриховыми линиями и отличаются от реальных данных в меньшую сторону. Как уже было отмечено выше, это, очевидно, связано с ограничениями формулы (11), в частности, с предположением об оптически тонком диске.

Рис. 9.

Исходные параметры модели диска и восстановленные параметры с помощью комплекса “SED-fitter”. Левый, средний и правый столбцы соответствуют моделям 81 тыс. лет, 83 тыс. лет и 1 млн. лет. Первый ряд показывает соответствие параметров в плоскости “темп аккреции–масса диска”; второй ряд: “внешняя граница диска–внутренняя граница диска”; третий ряд: “светимость звезды–возраст системы”; четвертый ряд: “масса диска–внешняя граница диска”. Оранжевым цветом обозначены модели с наиболее близкими спектральными распределениями энергии, полученные во время процесса подгонки. Синие звездочки соответствуют данным, самосогласованно полученным в ходе гидродинамического моделирования. Зеленые штриховые линии в верхнем ряду соответствуют массам, рассчитанным согласно выражению (11). Серые точки – полный набор моделей в использованном каталоге SED-Fitter.

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

Светимость звезды вполне успешно восстановлена для всех моделей диска (см. третий ряд рис. 9). Несмотря на небольшие отклонения, основная часть наилучших моделей сконцентрирована достаточно близко к реальным показателям светимости. О возрасте системы, напротив, достаточно трудно судить, исходя из данных, полученных в ходе подгонки спектров.

Нижний ряд рис. 9 иллюстрирует результаты восстановления параметров в пространстве: “масса диска–размер диска”. Эти параметры можно считать важнейшими наблюдательными характеристиками протозвездных дисков. Видно, что в рамках комплекса SED-Fitter не удается надежно восстановить оба параметра одновременно.

На основе полученных результатов можно утверждать, что комплекс SED-Fitter позволил относительно неплохо восстановить только часть исходных параметров модельного протопланетного диска. Основными причинами расхождений между входными и восстановленными параметрами диска являются, очевидно, различия между моделями диска, используемыми для генерации спектра и для его восстановления. Кроме того, как было отмечено в работе [31], задача о восстановлении параметров диска является вырожденной, т.е. различные особенности протопланетных систем могут проявляться на диаграммах спектрального распределения энергии одинаково. Все это свидетельствует о необходимости использования пространственно разрешенных наблюдений дисков в различных спектральных диапазонах, анализ которых должен осуществляться на базе физически согласованных моделей диска.

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

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

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

• Продемонстрированы различия в тепловой структуре диска, полученной в гидродинамическом моделировании и в вычислениях трехмерным кодом переноса излучения RADMC-3D. Показан вклад различных механизмов нагрева в диске. В молодых дисках важным источником нагрева во внутренних оптически толстых областях диска является вязкий нагрев. Другим существенным источником нагрева в неоднородных структурах диска может являться адиабатический нагрев сжатием (ударные волны). Однако учет внутренних источников нагрева в расчете переноса излучения с помощью RADMC-3D приводит с многократному повышению вычислительного времени. Это ведет к необходимости разработки и использованию специализированных алгоритмов расчета переноса излучения.

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