Астрономический журнал, 2020, T. 97, № 10, стр. 807-819

1. ВВЕДЕНИЕ

Формирование и эволюция протопланетных дисков (ПД) у молодых звезд являются одними из наиболее интригующих вопросов астрофизики. На самых ранних этапах эволюции многие молодые звездные системы проявляют признаки эпизодической аккреции; эти объекты известны как фуоры и экзоры (см., напр., обзоры [1, 2]). Есть теоретические предпосылки полагать, что фазу эпизодической аккреции проходят все ПД на начальных этапах своей эволюции, что и приводит к переменной светимости молодых звездных объектов [3]. Однако физические механизмы данной переменности до сих пор являются дискуссионными. Проблема переменности тесно связана с более общим вопросом о механизмах переноса углового момента в аккреционных дисках. В качестве механизмов, которые могут не только обеспечивать перенос углового момента в дисках, но и отвечать за нерегулярный характер аккреции, активно обсуждаются гравитационная, магниторотационная и тепловая неустойчивости (см. [2]). Так, нерегулярный характер аккреции, вызываемый гравитационной неустойчивостью, связан с формированием сгустков и их падением на звезду, что приводит к вспышкам светимости (см., напр., [4]). Вспышечный характер, обусловленный магниторотационной неустойчивостью, обусловлен положительной обратной связью этой неустойчивости со степенью ионизации вещества [5]. Тепловая неустойчивость обусловлена ростом непрозрачности газа с повышением температуры в условиях частично-ионизованного газа [6].

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

Целью нашей работы является более детальное исследование нерегулярного аккреционного режима, полученного в Статье I. Настоящая статья организована следующим образом: в разделе 2 дано краткое описание используемой модели. В разделе 3 изучено влияние параметров модели на характер аккреции в рамках простейшего описания падения вещества из оболочки, используемого в Статье I. Раздел 4 посвящен исследованию долговременной эволюции диска при более реалистичном описании аккреции вещества из оболочки. В Заключении кратко перечислены результаты данного исследования.

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

В данной работе используется модель эволюции протопланетного диска из Статьи I. Приведем здесь только ключевые элементы данной модели. Мы рассматриваем аксиально-симметричный геометрически тонкий вязкий кеплеровский диск без радиального градиента давления. Эволюция поверхностной плотности такого диска моделируется с помощью уравнения Прингла [12] с учетом аккреции вещества из оболочки на диск:

где $\Sigma (R,t)$ – поверхностная плотность, $R$ – расстояние до звезды, $t$ – время, $\nu (R,t)$ – коэффициент кинематической вязкости, $W(R,t)$ – темп притока вещества из оболочки в предположении, что удельный момент импульса оседающего вещества совпадает с таковым в диске. Для идентификации конвективно-неустойчивых областей диска одновременно с моделированием радиальной эволюции происходит восстановление вертикальной структуры диска. В приближении гидростатически-равновесного диска рассчитываются распределения плотности и температуры в полярном направлении. Основным фактором, управляющим эволюцией диска в рамках данной модели, является зависимость коэффициента вязкости $\nu (R,t)$ от радиуса. В модели принимается, что коэффициент вязкости является суммой фоновой и конвективной вязкости (идея о дополнительной вязкости в околозвездных дисках, связанной с конвекцией, более подробно обсуждается в работах [13–15]): Фоновая вязкость связана с некоторым неконвективным механизмом переноса углового момента (например, с магнитовращательной неустойчивостью [16, 17]) и обеспечивает непрерывную аккрецию газа. Эта вязкость задается в виде: где параметр $\beta = 1$ выбран из условия воспроизведения закона распределения плотности в наблюдаемых протопланетных дисках [18]. Источником возникновения конвекции в данной модели является выделение тепла за счет вязкой диссипации газа: где ${{\Gamma }_{{{\text{vis}}}}}$ – темп вязкой диссипации, приходящийся на единичную площадку диска на данном радиусе, ${{M}_{ * }}$ – масса звезды.

Конвективная вязкость ${{\nu }_{{\text{c}}}}$ отлична от нуля в конвективно-неустойчивых областях и рассчитывается в виде:

где $\gamma $ – массовая доля вещества в конвективно-неустойчивой области на данном радиусе, вычисленная в процедуре восстановления вертикальной структуры диска, $H$ – локальная высота диска, ${{V}_{{\text{c}}}}$ – характерная скорость конвекции. Рассчитанная таким образом конвективная вязкость ${{\nu }_{{\text{c}}}}$ сглаживается по радиусу на масштабе $H$, формируя ${{\tilde {\nu }}_{{\text{c}}}}$, входящую в формулу (2). В Статье I для нахождения характерной скорости конвекции ${{V}_{{\text{c}}}}$ мы использовали предположение, что вся выделившаяся в результате вязкой диссипации энергия переходит в кинетическую энергию его конвективного движения, что является верхней оценкой для этой величины. В данной работе мы дополнительно вводим коэффициент $\eta \leqslant 1$, который характеризует эффективность перехода тепловой энергии в конвективное движение: где ${{\rho }_{0}}$ – экваториальная плотность.

При реконструкции вертикальной структуры диска наряду с вязкой диссипацией газа учитывается нагрев излучением центрального объекта. При этом его светимость складывается из фотосферной светимости звезды (принимаем ее равной солнечной) и аккреционной светимости, рассчитываемой по формуле:

где ${{R}_{*}}$ – радиус звезды, ${{\dot {M}}_{*}}$ – темп аккреции из диска на звезду.

При решении уравнения (1) мы используем фиксированные значения поверхностной плотности на внутренней (0.2 а.е.) и внешней (200 а.е.) границах модельного диска. Значения плотности на границах выбраны сравнительно малыми, что реализует свободное вытекание вещества из рассматриваемой области. Как отмечалось в Cтатье I, более сложные граничные условия требуют отдельного изучения. Исследованию важности внутреннего граничного условия при моделировании протопланетных дисков посвящена, в частности, работа [19].

В данной работе мы исследуем влияние четырех параметров модели на проявление эпизодического характера аккреции: 1) темпа аккреции из оболочки на диск $\dot {M}$; 2) области питания, на которую аккрецирует вещество из оболочки ${{R}_{{{\text{ring}}}}}$; 3) величины фоновой вязкости ${{\nu }_{0}}$; 4) коэффициента эффективности конвекции $\eta $. В Статье I задавался постоянный приток газа из оболочки на диск в кольцо между 10 и 20 а.е. с темпом аккреции ${{10}^{{ - 7}}}\,{{M}_{ \odot }}$/год, а также были использованы коэффициенты ${{\nu }_{0}} = {{10}^{{15}}}$ см²/c и $\eta = 1$. Эту модель в дальнейшем мы будем называть базовой. В данной работе мы рассмотрели модели с темпами аккреции, отличающимися от базового на два порядка в обе стороны, что соответствует разбросу наблюдаемых в ПД темпов аккреции [20, 21]. Наряду с зоной питания 10–20 а.е. мы также исследуем случай, когда газ аккрецирует на диск в кольцо между 1 и 2 а.е., что более характерно для начальных этапов эволюции ПД.

Другим важным параметром модели является фоновая вязкость, которая предопределяет время жизни диска и его массу. В базовой модели мы использовали значение ${{\nu }_{0}} = {{10}^{{15}}}$ см²/c, которое соответствует относительно высокой турбулентной вязкости. Согласно оценкам из Статьи I, соответствующий этому значению альфа-параметр турбулентной вязкости [22, 23] равен $\alpha = 0.1$. В представленной статье мы дополнительно рассмотрели модель с фоновой вязкостью на порядок меньше базовой, что ближе к наблюдаемым оценкам для проэволюционировавших ПД. Коэффициент эффективности конвекции $\eta = 1$, принятый в базовой модели, заведомо переоценивает переход тепловой энергии в конвективную, поскольку часть энергии должна переноситься излучением. Для ответа на вопрос, не исчезнут ли вспышки, вызываемые конвекцией, если существенная часть выделившейся энергии переносится излучением, мы рассмотрели модель с $\eta = 0.1$. В табл. 1 перечислены параметры исследуемых моделей, а в разделе 3 описаны результаты моделирования.

В дополнение к этим моделям мы рассмотрели модель для изучения долговременной эволюции диска с более реалистичными темпом аккреции и зоной питания, параметры которых меняются со временем. Соответствующая последнему случаю функция $W(R,t)$ и результаты описаны в разделе 4.

3. РЕЖИМЫ ДИСКОВОЙ АККРЕЦИИ

В данном разделе представлены результаты для моделей с постоянным (по времени и пространству) притоком вещества в диск из оболочки. Соответствующий темп притока вещества внутри области питания имеет вид:

где $\dot {M}$ – интегральный темп притока из оболочки на диск, ${{R}_{{{\text{ring}},1}}}$ и ${{R}_{{{\text{ring}},2}}}$ – внутренний и внешний радиусы области питания. Мы проводили моделирование до начала возникновения вспышек, если таковые возникали, и исследовали их характеристики. При этом мы не интересовались дальнейшей эволюцией диска. В этих расчетах масса центральной звезды (${{M}_{*}} = 1\;{{M}_{ \odot }}$), ее радиус (${{R}_{*}} = 1\;{{R}_{ \odot }}$) и ее светимость (${{L}_{*}} = 1\;{{L}_{ \odot }}$) считались постоянными. В то же время учитывалась переменность интегральной светимости центрального объекта, связанная с аккрецией вещества из диска на звезду (см. формулу (7)).

4. ЭВОЛЮЦИЯ ПД ПРИ ПЕРЕМЕННОМ ПРИТОКЕ ВЕЩЕСТВА ИЗ ОБОЛОЧКИ

Для исследования долговременной эволюции диска необходимо задать реалистичную функцию $W(R,t)$, описывающую темп притока вещества из оболочки. Для вычисления данной функции мы используем приближение о сохранении локального углового момента вещества аккрецирующей оболочки – остатка родительского протозвездного облака. В рамках данного приближения элемент объема, первоначально находившийся на расстоянии $l$ от полярной оси, падает на так называемый центробежный радиус ${{R}_{c}}$:

на котором его угловая скорость становится равной кеплеровской. В данном выражении $\Omega $ – начальная угловая скорость рассматриваемого элемента, ${{M}_{*}}$ – текущая масса звезды. Таким образом, в принятой модели предполагается, что элементы облака последовательно оседают на диск с локальной кеплеровской скоростью, причем каждому элементу облака соответствует свой радиус осаждения, вычисляемый из условия сохранения момента импульса элемента.

Если предположить, что исходное протозвездное облако является сферически-симметричным и вращается твердотельно, то функция $W(R,t)$ имеет вид [24]:

где $\dot {M}(t)$ – текущий полный темп аккреции из оболочки на диск, ${{R}_{c}}(t)$ – граница области питания, т.е. центробежный радиус для аккрецируемого элемента из экваториальной плоскости. Для задания функций $\dot {M}(t)$ и ${{R}_{c}}(t)$ могут использоваться различные подходы [24]. Мы берем эти функции, аппроксимируя и экстраполируя результаты численного моделирования коллапса облака и последующей аккреции оболочки на звезду из работы [25]. Это можно сделать, поскольку в данной работе использовался лагранжевый метод, в рамках которого прослеживается эволюция отдельных элементов. При этом была использована угловая скорость начального облака $\Omega = {{10}^{{ - 14}}}$ c^–1, что является характерным значением для ядер молекулярных облаков [26]. На рис. 5 приведены используемые нами функции $\dot {M}(t)$, ${{R}_{c}}(t)$, а также вид функции $W(R,t)$.

Темп аккреции из оболочки на диск в интервале 0.2–0.7 млн. лет хорошо аппроксимируется экспоненциальной функцией. В связи с отсутствием данных для дальнейшей эволюции оболочки мы используем эту аппроксимацию и для больших времен. Центробежный радиус увеличивается со временем, достигая значения ≈180 а.е. на 0.8 млн. лет. Для больших времен используется постоянное значение 180 а.е. Из формы функции $W(R,t)$ видно, что падающее на диск вещество неравномерно заполняет область внутри области питания диска. В частности, максимум $W(R,t)$ вблизи нуля связан с падением вещества на диск из околополярных областей оболочки.

Поскольку в рамках данной модели исследуется долговременная эволюция диска, необходимо учитывать, что масса звезды растет за счет притока вещества из диска. Мы учитываем это, приняв, что начальная масса звезды составляет 0.3 ${{M}_{ \odot }}$ и увеличиваем ее в соответствии с аккрецированной массой. Одновременно с массой должны меняться также радиус звезды и ее фотосферная светимость, однако в целях упрощения мы пренебрегаем этим, полагая радиус и светимость звезды равными солнечным. Как и в модели с постоянным притоком, аккреционная светимость центрального объекта переменна и вычисляется по формуле (7). Отметим, что именно аккреционная светимость вносит больший вклад в светимость центрального объекта на ранних этапах эволюции диска. Значения коэффициентов ${{\nu }_{0}} = {{10}^{{15}}}$ см²/c и $\eta = 1$ были взяты из базовой модели и от времени не зависели.

Рассмотрим результаты расчета эволюции диска для данной модели. На рис. 6 приведены зависимость темпа аккреции из диска на звезду от времени, а также изменение массы диска со временем. Заполненная область в распределении темпа аккреции в интервале 0.17–3.7 млн. лет указывает на вспышечный режим аккреции – на данном масштабе рисунка многочисленные вспышки сливаются в единую сплошную полосу. После 3.7 млн. лет вспышки прекращаются и темп аккреции плавно уменьшается со временем. Из сравнения темпа аккреции на звезду с темпом притока вещества из оболочки (штриховая линия на левой панели рис. 6) можно сделать вывод о том, что важным процессом является накопление массы диском. В первый миллион лет в диске накапливается значительная масса (см. правую панель рис. 6) и его дальнейшая эволюция определяется перераспределением этой массы, в то время как приток вещества из оболочки становится пренебрежимо мал. Как максимальный, так и минимальный темпы аккреции во вспышечной фазе плавно уменьшаются в течение $ \approx {\kern 1pt} 2$ млн. лет, после чего остаются практически постоянными до 3.7 млн. лет. Отметим также, что изменение темпа аккреции в спокойной фазе ($t > 3.7$ млн. лет) хорошо согласуется с аналитической зависимостью $\dot {M} \propto {{t}^{{ - 5/4}}}$, которая соответствует диску с распределением вязкости $\nu \propto R$ (см. [27, формула (6)]).

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

В качестве примера проанализируем формирование вспышки для интервала времени в окрестности 1.5 млн. лет. На рис. 8 представлены поверхностная плотность диска и суммарный коэффициент вязкости $\nu (R)$ для трех близких моментов времени. В условно начальный момент времени вся внутренняя область вплоть до 30 а.е. является конвективной, что видно по высокому коэффициенту вязкости (правая панель рис. 8). С течением времени протяженность конвективной области уменьшается, т.е. ее граница двигается по направлению к звезде, достигая 2 а.е. на момент времени 227 лет. На момент времени 291 год внешняя граница этой конвективной зоны доходит до радиуса 0.35 а.е., и вскоре достигнет внутренней границы диска. На тот же момент времени видно, что образовалась новая конвективная зона внутри 0.7–5 а.е. Эта новая конвективная зона расширяется в обе стороны и впоследствии охватит всю внутреннюю зону вплоть до 30 а.е. Таким образом, новая конвективная фаза в диске начинает развиваться до того, как окончится предыдущая. Короткий пространственный интервал между границами конвективных зон (на рис. 8 интервал между 0.35 и 0.7 а.е.) и обеспечивает узкий, глубокий минимум перед максимумом аккреции.

Представленные выше результаты позволяют утверждать, что конвекция может являться важным фактором, обеспечивающим нерегулярный характер аккреции из диска на звезду. На рис. 9 (верхняя панель) показана долговременная эволюция радиального распределения поверхностной плотности и отмечены конвективно-неустойчивые области. Видно, что в процессе эволюции диска протяженность конвективно-неустойчивой области уменьшается с нескольких десятков до нескольких астрономических единиц, при этом фаза эпизодической аккреции длится менее 4 млн. лет. Эти результаты носят качественный характер, однако наша модель обладает рядом серьезных ограничений, отмеченных в Статье I. Снятие этих ограничений может существенно усложнить картину эволюции диска. Одним из таких ограничений является то, что в модели не учитываются испарение пыли и процессы диссоциации/ионизации газа, происходящие при высоких температурах. На средней панели рис. 9 показана эволюция распределения экваториальной температуры и отмечены области с температурой более 1500 K, при которых становятся важными процессы испарения пыли. Очевидно, что эти области сконцентрированы во внутренних частях диска и более выражены на начальных временах. Видно, что они частично перекрывают области конвективной неустойчивости, что, безусловно, должно сказаться на картине эволюции диска. В то же время область конвективной неустойчивости шире в пространстве и во времени, что оставляет выводы рассмотренной модели актуальными.

Другим ограничением модели является то, что в ней не учитывается самогравитация диска (уравнение (1) справедливо для кеплеровского диска). Между тем, как видно из рис. 6, масса диска на ранних этапах эволюции сопоставима с массой звезды. На нижней панели рис. 9 приведено распределение параметра Тоомре [28, 29]: $Q = \tfrac{{{{c}_{s}}\Omega }}{{\pi G\Sigma }}$, где ${{c}_{s}}$ – скорость звука, $\Omega $ – кеплеровская скорость. Низкие значения этого параметра ($Q < 1$) указывают на гравитационно-неустойчивые области. Видно, что эти области возникают на начальных временах эволюции ($t < 0.5$ млн. лет) во внешних частях диска ($R > 50$ а.е.). Возникновение данных областей также должно влиять на эволюцию диска – в нем должны возникать рукава, фрагменты, взаимодействие которых друг с другом и с диском ведет к сложной динамике (см., напр., [30]). Таким образом, на начальных временах эволюции диска конвекция может быть осложнена другими, возможно, более интенсивными управляющими процессами.

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

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

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

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

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

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

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

• Диск быстро становится конвективно-неустойчивым и остается таковым на протяжении почти 4 млн. лет. При этом неустойчивость захватывает область в несколько десятков астрономических единиц и далее постепенно сокращается.

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

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

БЛАГОДАРНОСТИ

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