MHD modeling of the molecular filament evolution

Texto integral

1. ВВЕДЕНИЕ

Современные наблюдения показывают, что межзвездные облака имеют волокнистую структуру, которая прослеживается от сверхоблаков H I до отдельных молекулярных облаков [1]. Под волокнами понимаются вытянутые структуры на картах излучения межзвездной среды. Подобные волокна могут представлять собой либо цилиндрические облака, либо плоские слои молекулярного газа, видимые с ребра [2]. Большинство протозвездных облаков, в которых происходит звездообразование, располагаются внутри волокнообразных молекулярных облаков [3]. Поэтому исследование структуры и эволюции молекулярных волокон имеет важное значение для построения теории образования звезд.

Наблюдения показывают, что характерная ширина молекулярных волокон составляет порядка 10⁻¹ пк, а длина варьируется от нескольких пк до сотен пк. Температура в волокнах лежит в диапазоне от 10 до 25 К, а плотность вещества — от 10⁴ до 10⁵ см⁻³[2].

Поляризационное картирование молекулярных облаков указывает на то, что они пронизаны крупномасштабным магнитным полем [4]. В разреженных облаках магнитное поле, как правило, направлено вдоль главной оси волокна, в плотных облаках — поперек. Согласно измерениям зеемановского расщепления линий OH и оценкам с помощью метода Чандрасекара—Ферми, установлено, что интенсивность магнитного поля волокон увеличивается с лучевой концентрацией N и лежит в диапазоне от 10^–5 Гс для разреженных облаков с N = 10¹⁹ см⁻² до 10⁻³ Гс для наиболее плотных волокон с N = 10²³ см⁻².

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

В уединенных однородных изотермических волокнах гравитационная фокусировка приводит к образованию плотных ядер на концах волокна. Данный механизм иногда называется в литературе “коллапс с конца” (англ., end–dominated collapse [6]) или “фрагментация на концах” (edge fragmentation [5]). Примером такого облака является волокно S242, на концах которого наблюдаются ядра плотностью порядка 10⁵ см⁻³и размерами 1 пк [7].

Малые продольные возмущения волокон могут привести к развитию гравитационной неустойчивости [8, 9, 10]. В случае волокон без магнитного поля, неустойчивость развивается для волн, длина которых более чем в 4 раза превышает радиус однородной части волокна. Неустойчивость приводит к формированию гравитационных перетяжек и впоследствии — ядер, распределенных вдоль волокна с характерным расстоянием между ними порядка длины волны наиболее быстрорастущей моды неустойчивости. Форма ядер, образующихся в результате гравитационной фрагментации, близка к сферической [11]. Примером объекта, в котором наблюдаются признаки гравитационной фрагментации, является волокно NGC 2024S/Orion B. Ядра в этом волокне имеют размер порядка 10⁻² пк и массы ~1M_e [12]. Профиль скорости вдоль волокна имеет периодический характер с длиной волны λ ~ 0.2 пк. Ядра смещены относительно максимумов скорости на λ/4, что явно указывает на гравитационную фрагментацию волокна. Другим примером подобного волокна является WB 673 [13].

С точки зрения интерпретации наблюдений особый интерес представляет моделирование фрагментации волокон с магнитным полем и определением условий их фрагментации. В работе Сейфрида и Уолш [14] с помощью численного кода FLASH выполнено моделирование эволюции волокон с учетом турбулентности и различных ориентаций магнитного поля. В зависимости от начальных условий, авторы выделили несколько режимов фрагментации волокна: фрагментация на концах, гравитационная фрагментация волокна с последующим образованием равномерно распределенных волокон, глобальный коллапс волокна по направлению к центру облака. Авторы показали, что градиент давления продольного магнитного поля поддерживает почти постоянную толщину волокна порядка 0.1 пк. Впоследствии этот вывод был подтвержден в МГД-расчетах Дудорова и Хайбрахманова [2].

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

В разделе 2 описана постановка задачи, основные уравнения и численный код FLASH, который используется для решения системы уравнений. В параграфе 3.1 приведены результаты расчетов эволюции волокна без магнитного поля и со слабым магнитным полем. Результаты расчетов с более сильным магнитным полем приведены в параграфе 3.2. В параграфе 3.3 описаны характеристики ядер, образующихся в расчетах с магнитным полем. В заключении суммированы результаты и сформированы выводы.

2. Модель

2.1. Постановка задачи

В работе моделируется гравитационный коллапс цилиндрического молекулярного облака (волокна) длиной H₀ = 10 пк и радиусом r₀ = 0.2 пк. Молекулярный вес газа равен m = 2.31, температура T₀ = 10 К, концентрация n₀ = 10⁵ см⁻³. Линейная масса волокна M/L = 658M_⊚ /пк превышает критическое значение (M/L)_crit = 16.6M_⊚ /пк [9, 10], т.е. волокно является гравитационно неустойчивым. Модели коллапса с учетом переноса излучения показывают, что в диапазоне концентраций n = [10⁵,10¹¹] см⁻³тепловая энергия сжимающегося газа эффективно высвечивается и температура газа остается примерно постоянной. На основе этого вывода в настоящей работе для простоты принимается, что газ характеризуется уравнением состояния с эффективным показателем адиабаты γ = 1.001, соответствующим изотермическому сжатию. Соответствующее отношение тепловой энергии к модулю гравитационной энергии ε_T = 0.003. Скорость звука в волокне равна c_s = 0.19 км/с.

Для исследования роли магнитного поля в эволюции волокон проведены три расчета с разными значениями отношения магнитной энергии облака к его модулю гравитационной энергии: ε_m = 0 (расчет ГД), ε_m = 0.03 (расчет МГД-1), ε_m = 0.27 (расчет МГД-2). Соответствующие интенсивности магнитного поля равны: B = 0, 1.9×10⁻⁴, 6×10⁻⁴ Гс соответственно. В обоих случаях магнитное поле параллельно волокну. Волокно находится в равновесии по давлению с внешней средой с концентрацией и температурой n = 10⁴ см⁻³и T = 100 К. Характерное время свободного падения для выбранной плотности составляет t_ff ≈ 10⁵ лет.

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

$R_m=\frac{v_0{l}_0}{\eta }$                                                                                              (1)

где v₀ и l₀ — характерные скорость газа и пространственный масштаб, η — коэффициент диффузии магнитного поля. В качестве l₀ выберем радиус волокна r₀, в качестве характерной скорости — v₀ = r₀/t_ff . Диссипация магнитного потока может быть обусловлена омической диффузией (ОД) или магнитной амбиполярной диффузией (МАД). Для оценок соответствующих коэффициентов используем формулы из статьи Дудорова и Хайбрахманова [15]:

$\eta =\left\{\begin{array}{l}480x^{-1}{T}^{1/2}{\mathrm{см}}^2/\mathrm{с} \left(\mathrm{ОД}\right)\\ \frac{B^2}{4\pi x{\rho }^2{\mathrm{\eta }}_{\mathrm{i}n}} \left(\mathrm{МФД}\right)\end{array}\right.$                                                                (2)

где x — степень ионизации, T — температура газа, B — магнитная индукция, ρ — плотность газа,  — коэффициент взаимодействия ионов массой m_i = 30m_H и нейтралов массой m_n = μm_H, где ⟨σv⟩_in = 2×10⁻⁹ см³/с, m_H — масса атома водорода. В рассматриваемом диапазоне плотностей степень ионизации можно оценить из баланса ионизацией космическими лучами со скоростью x и лучистыми рекомбинациями: , где a_r = 6.21×10⁻¹¹T ^–½ см³/с — коэффициент лучистых рекомбинаций, n — концентрация газа [16].

Используя типичные параметры для межзвездной среды, получим для случая омической диффузии:

$R_m^{ОД}=1.7·{10}^{14}\left(\frac{r_0}{0.2\mathrm{пк}}\right)\times {\left(\frac{\xi }{{10}^{-17}{с}^{-1}}\right)}^{1/2}{\left(\frac{T_0}{10K}\right)}^{-1/4}$                                         (3)

для амбиполярной диффузии:

$R_m^{\mathrm{МАД}}=6.6\times {10}^3{\left(\frac{r_0}{0.2\mathrm{пк}}\right)}^2\left(\frac{B_0}{1.9\times {10}^{-4}\mathrm{Гс}}\right){\left(\frac{\xi }{{10}^{-17}{с}^{-1}}\right)}^{1/2}{\left(\frac{n_0}{{10}^5{\mathrm{см}}^{-3}}\right)}^2{\left(\frac{T_0}{10K}\right)}^{1/4}.$ (4)

Формулы (3) и (4) показывают, что при выбранных параметрах R_m≫1, то есть магнитное поле вморожено в газ.

2.2. Основные уравнения и методы решения

Исследуем эволюцию молекулярного волокна с помощью системы уравнений идеальной МГД:

$\frac{\partial \mathrm{\rho }}{\partial t}+\nabla \left(\mathrm{\rho }\mathbf{v}\right)=0$,                                                                                            (5)

$\frac{\partial \mathrm{\rho }}{\partial t}+\left(\mathbf{v}\nabla \right)\mathbf{v}=-\frac{1}{\mathrm{\rho }}\nabla P\mathit{-}\nabla \mathrm{\Phi }-\frac{1}{4\mathrm{\pi \rho }}\mathbf{B}\times \left(\nabla \times \mathbf{B}\right)$,                                               (6)

$\frac{\partial \mathbf{B}}{\partial t}=\nabla \times \left(\mathbf{v}+\mathbf{B}\right)$                                                                                          (7)

$\begin{array}{l}\frac{\partial }{\partial t}=\left[\rho \left(\epsilon +\frac{v^2}{2}+\mathrm{\Phi }\right)+\frac{B^2}{8\pi }\right]=\\ -\nabla \left[\rho v\left(\epsilon +\frac{v^2}{2}+\frac{P}{\rho }+\mathrm{\Phi }\right)+\frac{1}{4\pi }\mathbf{B}\mathbf{\times }\left(v\mathrm{\times }B\right)\right],\end{array}$                                                  (8)

$∆\mathrm{\Phi } = 4\mathrm{\pi G\rho },$                                                                                             (9)

$P = \left(\gamma -1\right)\rho \epsilon ,$                                                                                           (10)

где r, v и P — плотность, вектор скорости и давление газа, Ф — гравитационный потенциал, B — магнитная индукция, ε — внутренняя энергия газа, G — гравитационная постоянная, g — показатель адиабаты.

Для моделирования эволюции волокна используем численный код FLASH 4, в котором реализована технология адаптивно–встраиваемых сеток (AMR) [17]. В коде уравнения идеальной МГД (5)–(8) решаются с помощью схемы MUSCL годуновского типа [18]. Рассматривается трехмерная постановка задачи в декартовых координатах. Ось  соответствует оси симметрии волокна. Размеры расчетной области в направлениях x × y × ɀ составляют 1.93 × 1.93 × 12.9 пк³, используется 7 уровней вложенности AMR-сетки. Размеры самой большой ячейки в направлениях x × y × ɀ равны 0.24 × 0.24 × 1.61 пк³, размеры самой малой — 0.0037 × 0.0037 × 0.025 пк³. Эффективное разрешение сетки на 7-м уровне вложенности составляет 512 × 512 × 512. Гравитационный потенциал находится с помощью древесного алгоритма Барнса—Хата [19].

3. Результаты расчетов

3.1. Общая картина эволюции волокна

 

Рис. 1. Распределение плотности волокна в плоскости x – ɀ в ГД-расчете в моменты времени t = 0 (а), 0.8t_ff (б), 1t_ff (в).

 

На рис. 1 приведены срезы плотностей в плоскости x - ɀ для ГД-расчета в моменты времени t = 0, 0.8t_ff , 1t_ff . Расчеты показывают, что в случае без магнитного поля волокно свободно сжимается вдоль радиуса и к моменту времени t = 1t_ff плотность в центре волокна вырастает на 3 порядка, толщина волокна вдоль радиуса составляет 0.004 пк. Фрагментации волокна не происходит.

 

Рис. 2. Распределение плотности (цветовая заливка) и линий магнитного поля (черные линии со стрелками) волокна в плоскости x – ɀ в расчете МГД-1 в моменты времени t = 0 (а), 1t_ff (б), 1.28t_ff (в).

 

На рис. 2 приведены срезы плотностей в плоскости x - ɀ для расчета МГД-1 в моменты времени t = 0, 1t_ff , 1.28t_ff . В расчете МГД-1 к моменту времени t = 1t_ff волокно сжимается до радиуса r ≈ 0.1 пк и плотности n ≈ 1.77×10⁸ см⁻³. К этому моменту времени эффективный показатель адиабаты увеличивается от 1 до 2, сжатие вдоль r останавливается градиентом магнитного давления, и в дальнейшем облако осциллирует вдоль радиуса.

На концах облака к моменту времени t = 1.28t_ff выделяются ядра плотностью n = 1.7×10⁸ см⁻³. Свойства ядер обсуждаются далее в разделе “Характеристики образующихся ядер”.

3.2. Влияние магнитного поля на эволюцию магнитного поля

 

Рис. 3. Распределение плотности и линий магнитного поля волокна в плоскости x – ɀ в расчете МГД-2 в моменты времени t = 0 (а), 1t_ff (б), 1.28t_ff (в), 1.9t_ff (г).

 

На рис. 3 приведены срезы плотностей в плоскости x – ɀ для расчета МГД-2 в моменты времени t = 0.1t_ff , 1.28t_ff , 1.9t_ff . В расчете МГД-2 картина коллапса качественно похожа на результаты для варианта МГД-1. Волокно сжимается до радиуса r ≈ 0.3 пк и концентрации n ≈ 2×10⁵ см⁻³к моменту времени 1t_ff , после чего сжатие останавливается и начинаются осцилляции волокна вдоль радиуса. Плотность образующихся на торцах волокна ядер к моменту времени 1.9t_ff составляет n ≈ 2×10⁷ см⁻³.

 

Рис. 4. Панель (а): профили плотности волокна вдоль оси ɀ для МГД-расчетов в моменты времени t = 0, 1.28t_ff и 1.9t_ff . Панель (б): профили скорости vɀ вдоль оси ɀ для МГД-расчетов в моменты времени t = 1.28t_ff и t = 1.9t_ff .

 

На рис. 4 показана эволюция профилей плотности и скорости вдоль волокна в МГД-расчетах. Расположение торцов волокна в момент времени t = 0 соответствует координатам ɀ_L = 1.6 пк и ɀ_R = 11.3 пк. В расчете МГД-1 к моменту времени t = 1.28t_ff на концах волокна (ɀ_L = 1.9 пк и ɀ_R = 11 пк) наблюдаются пики плотностей n = 1.7×10⁸ см⁻³и скоростей v_ɀ = 3.6 км/с. Пики соответствует ядрам, образующимся на концах волокна. Скорость левого ядра положительна, правого — отрицательна. Для данных значений скоростей число Маха составляет M = 6, то есть ядра двигаются со сверхзвуковыми скоростями навстречу друг к другу.

 

Рис. 5. Распределение плотности (цветовая заливка), поля скоростей (зеленые стрелки) и линий магнитного поля (черные линии со стрелками) в области формирования ядер в расчете МГД-1 в моменты времени t = 1t_ff (а) и 1.28t_ff (б).

 

В расчете МГД-2 плотности и скорости ядер, расположенных в момент времени t = 1.9t_ff на координатах ɀ_L = 2.15 пк и ɀ_R = 10.75 пк, равны соответственно n = 2×10⁸ см⁻³и | v_ɀ | = 5 км/с, а соответствующее число Маха равно M = 8.3.

3.3 Характеристики образующихся ядер

На рис. 5 приведено распределение плотности, линий магнитного поля и поля скоростей для расчета МГД-1 в области формирования ядер в моменты времени t = 1 и 1.28t_ff в плоскости x – ɀ. В силу симметрии, для примера показаны только ядра, находящиеся на “левом” торце волокна. На рис. 6 приведено аналогичное распределение для расчета МГД-2 в моменты времени t = 1.28t_ff и 1.9t_ff в плоскости x – ɀ.

 

Таблица 1. Характеристики ядер в МГД-расчетах

Расчет	Время, tff	dr , пк	dɀ, пк	n, см⁻³	M, M⊚	vɀ, км/с
МГД-1	1.28	0.0075	0.025	1.7×10⁸	12.5	3.6
МГД-2	1.28	0.06	0.025	5.3×10⁵	2.5	2.7
	1.9	0.03	0.025	2×10⁷	23.6	5.3

 

Рис. 6. Распределение плотности (цветовая заливка), поля скоростей (зеленые стрелки) и линий магнитного поля (голубые изолинии) в области формирования ядер в МГД-2 расчете в моменты времени t = 1t_ff (а), 1.28t_ff (б), 1.9t_ff (в).

 

В табл. 1 приведены следующие характеристики ядер в МГД-расчетах: размеры вдоль радиуса волокна r и оси ɀ (столбцы 3 и 4), концентрация n (столбец 5), масса M (столбец 6) и скорость v_ɀ (столбец 7). Таблица показывает, что в расчете с большей начальной интенсивностью магнитного поля образуются ядра большего радиуса d_r и меньшей плотности. Это связано с тем, что в случае большей начальной интенсивности магнитного поля усиливается влияние градиента магнитного давления на динамику волокна. Размеры ядер вдоль оси ɀ не зависят от интенсивности начального магнитного поля, так как продольное магнитное поле не препятствует сжатию вдоль волокна. Масса ядер растет в процессе их движения и достигает 12.5M_⊚ в расчете МГД-1 и 23M_⊚ в расчете МГД-2.

4. Заключение

В работе исследовано влияние магнитного поля на эволюцию молекулярных волокон и на характеристики ядер, формирующихся в волокнах в результате фрагментации на концах. Для этого проведено численное МГД-моделирование гравитационного коллапса цилиндрического молекулярного волокна с различными значениями продольного магнитного поля.

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

В МГД-расчетах градиент магнитного давления препятствует коллапсу и приводит к затухающим колебаниям волокна вдоль радиуса. В ходе эволюции на концах волокна образуются плотные сгустки (ядра), которые передвигаются к центру облака со сверхзвуковыми скоростями в диапазоне от 3 до 5 км/с. В облаках с большей начальной интенсивностью образуются ядра больших размеров и меньшей плотности, поскольку с увеличением интенсивности увеличивается влияние градиента давления магнитного поля вдоль радиуса. Масса ядер увеличивается в процессе эволюции волокна и лежит в диапазоне » H(10–20)M_e.

Проведенные расчеты указывают, что коллапс с конца является естественным результатом эволюции волокна с продольным магнитным полем. Можно предположить, что волокна, в которых уплотнения наблюдаются только на концах (пр., [7]), поддерживаются от гравитационной фрагментации продольным магнитным полем. Дополнительную поддержку против гравитации может оказывать турбулентность в волокне [14, 20].

Дальнейшее развитие представленной модели предполагает моделирование эволюции волокон с учетом их неоднородности, а также вращения и/или внутренней турбулентности. Особый интерес представляет исследование гравитационной фрагментации волокна с продольным магнитным полем в рамках механизма гравитационных перетяжек Чандрасекара—Ферми [8].

ФИНАНСИРОВАНИЕ

Работа выполнена при финансовой поддержке Фонда перспективных научных исследований ЧелГУ (проект 2023/7). Расчеты проводились на вычислительном кластере Челябинского государственного университета. С. А. Хайбрахманов выражает благодарность за финансирование в рамках гос. задания по проекту FEUz-2020-0038.

КОНФЛИКТ ИНТЕРЕСОВ

Авторы данной работы заявляют, что у них нет конфликта интересов.

Sobre autores

I. Sultanov

Chelyabinsk State University

Autor responsável pela correspondência
Email: syltahof@yandex.ru
Rússia, Chelyabinsk

S. Khaibrakhmanov

Saint Petersburg State University; Chelyabinsk State University; Ural Federal University

Email: syltahof@yandex.ru
Rússia, Saint Petersburg; Chelyabinsk; Yekaterinburg

Bibliografia

Arquivos suplementares