Астрономический журнал, 2022, T. 99, № 3, стр. 189-210

2.1. Структура и состав пылинок

В начальный момент времени пылинки являются двухслойными шариками: ядро с плотностью ${{\rho }_{{\text{c}}}}$ и радиусом ${{a}_{{\text{c}}}}$ покрыто мантией с плотностью ${{\rho }_{{\text{m}}}}$ и радиусом ${{a}_{{\text{m}}}}$ (измеряемым от центра пылинки), причем отношение ${{a}_{{\text{c}}}}{\text{/}}{{a}_{{\text{m}}}} \equiv k$ не зависит от размера пылинки. Параметр $k$ вычисляется через начальное отношение масс мантии и ядра ${{m}_{{\text{m}}}}{\text{/}}{{m}_{{\text{c}}}}$, а также через их плотности:

Массы мантии и ядра ограничены распространенностью элементов в моделируемых облаках. Поскольку она неизвестна, используется распространенность элементов в современной солнечной системе [40]. Считается, что в ядре доминирует форстерит (Mg₂SiO₄), а в мантии – вода. Поэтому масса ядра ограничена распространенностью Mg (0.08%) и Si (0.09%) по массе (относительно водорода), а масса мантии ограничена распространенностью O (0.86%). Полагая, что основная часть этих атомов сосредоточена в пыли, и учитывая, что масса кислорода в форстерите в 1.2 раза меньше массы магния и кремния, получаем оценку ${{m}_{{\text{m}}}}{\text{/}}{{m}_{{\text{c}}}} = (0.86{\kern 1pt} - {\kern 1pt} 0.17{\text{/}}1.2){\text{/}}(0.17 + 0.17{\text{/}}1.2) = 2.3$. Эта оценка согласуется с типичным отношением масс льдов и силикатов в молодых звездных объектах [41, п. 2.6].

Предполагается, что и форстерит, и водяной лед имеют пористость 20%, поэтому ${{\rho }_{{\text{c}}}} = 2.6$ г/см³, ${{\rho }_{{\text{m}}}} = 0.8$ г/см³. Тогда из (2) следует $k = 0.5$. Поскольку мантии и ядра однородны по составу, то их испаряющаяся поверхность остается сферической. После исчезновения мантий пылинки становятся однослойными шариками.

2.2. Эволюция распределения пылинок по радиусам

Распределение пылинок по радиусам $a$ в момент времени $t$ представлено плотностью $f(a,t)$ – числом пылинок в единичном интервале радиусов в единице массы вещества (размерность $[f] = {{{\text{г}}}^{{ - 1}}}\;{\text{с}}{{{\text{м}}}^{{ - 1}}}$). Такой выбор $f(a,t)$ позволяет легко учесть сжимаемость среды. Начальная функция $f(a,t)$ считается степеннóй:

где $x \equiv a{\text{/}}{{a}_{{01}}}$ – безразмерный радиус пылинок, ${{a}_{{01}}}$ – начальный радиус наименьших пылинок, ${{f}_{{01}}} = f({{a}_{{01}}},0)$ – нормировочный коэффициент. Во всех моделях $q = - 3.5$. Вначале $x \in [1,A]$, где $A = 100$ – начальное отношение максимального и минимального радиусов.

Для простоты модели предполагается, что 1) коагуляция и фрагментация пылинок уравновешивают друг друга и не влияют на $f(a,t)$; 2) к начальному моменту времени почти вся вода уже сконденсировалась на мантиях, поэтому ее дальнейшей конденсацией можно пренебречь; 3) почти весь кремний уже находится в форстерите, поэтому рост ядер пылинок после потери мантий пренебрежимо мал; 4) пылинки разрушаются только путем отделения атомов и молекул от внешнего слоя: сначала при сублимации воды, затем при термическом разложении форстерита. Для краткости оба эти процесса называются испарением. Испарение через поры не рассматривается. Также предполагается, что 5) температуры газа и пыли одинаковы, что приемлемо в моделях коллапса протозвездных облаков при $n{{ > 10}^{5}}$ см^–3 (см., напр., [42]), кроме случаев, когда пыль проникает внутрь ультракомпактных HII зон возле массивных протозвезд [43]. Считается, что 6) температура среды растет достаточно плавно, чтобы мантии пылинок успевали испариться при $T \approx 130$ К до начала испарения ядер при $T \approx 1150$ К. Это условие может не выполняться на скачках температуры с $ \sim {\kern 1pt} 100$ К до ${{10}^{3}}{\kern 1pt} - {\kern 1pt} {{10}^{4}}$ К, например, на фронтах ионизации водорода возле массивных протозвезд [43] и на ударных волнах, ограничивающих протозвезды (см. модели с массой 2 и 14 ${{M}_{ \odot }}$ в [22]).

Начальные радиусы ядер наименьших пылинок приравниваются к минимальному радиусу стабильных наночастиц ${{a}_{{\min }}}$, при котором испарение вещества наночастиц еще отличается от их взрывного разрушения, т.е. связи атомов и молекул внутри наночастиц сохраняются. Из этого предположения следует, что ${{a}_{{01}}} = {{a}_{{\min }}}{\text{/}}k$. Считается, что ${{a}_{{\min }}} = 5$ нм, но при более детальном подходе необходимо задавать ${{a}_{{\min }}}$, исходя из состава и структуры ядер (см., напр., [35]).

Темп потери массы пылинки, состоящей из сферических слоев, $\dot {m} = - 4\pi {{a}^{2}}w$, где $w > 0$ – темп потери массы с единицы площади поверхности слоя, являющегося внешним в данный момент. С другой стороны, $\dot {m} = 4\pi {{a}^{2}}\dot {a}\tilde {\rho }$, где $\dot {a}$ – скорость изменения радиуса, $\tilde {\rho }$ – плотность внешнего слоя. Поэтому скорость изменения радиуса пылинки при испарении вещества с ее поверхности не зависит от радиуса явно и определяется только характеристиками внешнего слоя [36]:

Однако $\dot {a}$ зависит от $a$ неявно: и $w$, и $\tilde {\rho }$ меняются скачком, когда текущий радиус пылинки проходит через внутренний радиус слоя, т.е. слой исчезает. В рассматриваемой двухслойной модели пылинок при исчезновении мантии $w$ резко уменьшается, поскольку ядро является менее летучим, а $\tilde {\rho }$ резко увеличивается, поскольку ядро является более плотным, чем мантия.

Термическое разложение в отличие от сублимации относится к химическим реакциям, однако темп потери массы при этих процессах может быть рассчитан единым образом с помощью уравнения Герца-Кнудсена, использующего давление насыщенного пара $\tilde {p}$. При температуре интенсивного испарения оно много больше давления ненасыщенного пара, созданного пылинками. Поэтому можно пренебречь последним и записать уравнение Герца-Кнудсена в виде (см., напр., [44])

где $\tilde {v}$ – средний модуль скорости атомов или молекул, образующих пар. Для паров H₂O известна простая, но достаточно точная аппроксимация экспериментальных данных [45]: где температура $T$ выражена в градусах Кельвина. При термическом разложении форстерита основной вклад в давление паров дают молекулы SiO и атомы Mg [36]. Следуя [39], рассчитаем темп потери массы через давление паров SiO, у которых [36]

Скорость изменения радиуса пылинки $\dot {a}$ есть скорость ее движения в одномерном пространстве размеров. Когда мантия исчезает, $\dot {a}$ скачком уменьшается, поскольку при характерной температуре испарения воды 130 К функция (4) близка к нулю. Значит, эволюция пылинок, потерявших мантию, приостанавливается до тех пор, пока при $T \approx 1150$ К испарение форстерита не станет заметным. Поскольку начальное отношение радиусов ядра и мантии $k$ одинаково для всех пылинок, толщина мантии у мелких пылинок существенно меньше, чем у крупных (см. рис. 1). Поэтому мелкие пылинки теряют мантию раньше крупных, и делают остановку в пространстве размеров тоже раньше.

Монотонная зависимость начального радиуса ядра от начального радиуса мантии исключает перемешивание пылинок в пространстве размеров в моменты остановок и позволяет описать эволюцию пыли полуаналитическими методами, а именно выразить $f(a,t)$ через элементарные функции и их интегралы по времени.

Пока $\dot {a}$ не зависит от $a$, функция $f(a,t)$ при испарении вещества пылинок смещается влево без деформаций. Смещенная функция (3) имеет вид

Здесь $\epsilon $ – относительное смещение пылинок в пространстве размеров, $x - \epsilon $ – начальная безразмерная координата пылинок, оказавшихся в момент времени $t$ в точке $x$. Потеряв мантию, пылинка останавливается в пространстве размеров, но следующая за ней более крупная пылинка продолжает движение. Поэтому на этом участке функция $f(a,t)$ деформируется – сжимается и возрастает.

Предположения о двухслойности пылинок и плавном росте их температуры позволяют описать эволюцию пыли в рамках трех этапов. На рис. 1 схематично показана эволюция пылинок трех размеров и стрелками указаны этапы эволюции. На этапе 1 мантии испаряются у всех пылинок, поэтому в пространстве размеров они движутся все, и справедливо решение (5). Этот этап заканчивается, когда останавливаются наименьшие пылинки: при $\epsilon = k - 1$. На этапе 2 пылинки без мантий стоят, а пылинки с мантиями движутся. Этот этап заканчивается, когда останавливаются наибольшие пылинки: при $\epsilon = A(k - 1)$. На этапе 3 снова движутся все пылинки: их ядра испаряются и, достигнув радиуса ${{a}_{{{\text{min}}}}}$, исчезают. Следовательно, число пылинок в единице массы вещества на этапе 3 уменьшается.

Распределение пылинок по радиусам на этапе 2 состоит из двух частей: левую часть создают пылинки без мантий, а правую часть – пылинки с мантиями. Граница между ними – точка остановки c координатой ${{a}_{{\text{s}}}}$ – постоянно смещается вверх по течению, подобно краю пробки на дороге. Найти закон движения этой границы можно из следующих соображений. Пылинки с начальным радиусом ${{a}_{0}} + d{{a}_{0}}$ облачены в мантию на время $d{{t}_{{\text{s}}}}$ дольше, чем пылинки с начальным радиусом ${{a}_{0}}$, и останавливаются на $d{{a}_{{\text{s}}}}$ правее. Поскольку вначале их разделяет расстояние $d{{a}_{0}}$, то $d{{a}_{{\text{s}}}} = d{{a}_{0}} + \dot {a}d{{t}_{{\text{s}}}}$. С другой стороны, $d{{a}_{0}} = d{{a}_{{\text{s}}}}{\text{/}}k$. Поэтому $d{{a}_{{\text{s}}}} = k\dot {a}d{{t}_{{\text{s}}}}{\text{/}}(k - 1)$. Тогда, приравнивая $d{{t}_{{\text{s}}}}$ к $dt$, находим зависимость безразмерной координаты точки остановки от времени

Деформацию распределения пылинок по радиусам в точке остановки на этапе 2 можно найти из следующих соображений. Пылинки с начальными радиусами ${{a}_{0}}$ и ${{a}_{0}} + d{{a}_{0}}$ после испарения мантий имеют радиусы ${{a}_{{\text{s}}}} = k{{a}_{0}}$ и $a_{{\text{s}}}^{'} = k({{a}_{0}} + d{{a}_{0}})$. Поэтому все пылинки из интервала с начальной шириной $d{{a}_{0}}$ оказываются в более узком интервале с шириной $d{{a}_{{\text{s}}}} = kd{{a}_{0}}$. В терминах лагранжева описания в пространстве размеров это означает, что лагранжев элемент с начальной шириной $d{{a}_{0}}$ сжимается до ширины $kd{{a}_{0}}$. Поскольку число пылинок в этом элементе сохраняется, то лагранжева плотность распределения возрастает в $1{\text{/}}k$ раз.

На этапе 3 все сжатое распределение $f(a,t)$ смещается влево без деформаций и обрезается на границе $a = {{a}_{{\min }}}$. Относительное смещение пылинок в пространстве размеров на этом этапе

где ${{t}_{{\text{c}}}}$ – время начала испарения ядер.

Введя переменные ${{x}_{*}} \equiv x{\text{/}}k$ и $\chi \equiv {{x}_{{\text{s}}}}{\text{/}}k = $ ${{x}_{{\text{s}}}} - \epsilon = \epsilon {\text{/}}(k - 1)$, запишем распределение пылинок по радиусам на всех трех этапах эволюции:

где первая строка относится к этапу 1, вторая и третья строки – к этапу 2, четвертая – к этапу 3.

Нормировочный коэффициент ${{f}_{{01}}}$ вычисляется из условия, что интеграл от начальной массы пылинок по всему распределению (3) должен быть равен начальному отношению массы пыли к массе вещества ${{Y}_{{\text{d}}}}$:

где ${{m}_{{01}}} = (4{\text{/}}3)\pi a_{{01}}^{3}{{\rho }_{{\text{g}}}}$ – начальная масса наименьших пылинок, ${{\rho }_{{\text{g}}}} = {{\rho }_{{\text{m}}}} + {{k}^{3}}{{\rho }_{{{\text{cm}}}}}$ – начальная плотность материала пылинок, усредненная по их объему, ${{\rho }_{{{\text{cm}}}}} \equiv {{\rho }_{{\text{c}}}} - {{\rho }_{{\text{m}}}}$.

На рис. 2 приведены примеры распределения пылинок по радиусам (6) при нескольких выбранных значениях $\chi $ и $\delta $ (без привязки ко времени). В логарифмическом масштабе прямыми линиями являются графики $f(a,t)$ на старте и на этапе 2 слева от точки остановки (пылинки без мантий). В остальных случаях графики $f(a,t)$ являются выпуклыми кривыми. Видно, что в точке остановки функция $f(a,t)$ терпит разрыв. Постепенно он смещается вправо – вверх по течению частиц в пространстве размеров. Поскольку радиусы пылинок только уменьшаются, все графики $f(a,t)$ сдвигаются влево. Поскольку $k = 0.5$, то $\max [f(a,t)] = 2\max [f(a,0)]$ (этап 2, левая граница: $a = {{a}_{{\min }}}$).

2.3. Интегральные характеристики пыли

В качестве интегральных характеристик пыли рассматриваются радиус $\langle a\rangle $, площадь поперечного сечения $\langle s\rangle $ и масса пылинки $\langle m\rangle $, усредненные по распределению радиусов пылинок; ${{Y}_{{{\text{ev}}}}}$ – отношение испарившейся массы пыли к массе вещества, ${{X}_{{\text{g}}}}$ – отношение концентраций пылинок и частиц газа.

Введем интеграл

Тогда из (6) следуют выражения для средней величины радиуса пылинок Средняя площадь поперечного сечения $\langle s\rangle = \pi \langle {{a}^{2}}\rangle $. Средняя масса пылинки вычисляется более сложно из-за различия плотностей мантий и ядер: Отношение испарившейся массы пыли к массе вещества:

Пусть ${{N}_{{\text{g}}}} = {{N}_{{\text{g}}}}(t) = \int \,f(a,t)da$ – число пылинок в единице массы вещества. Тогда отношение концентраций пылинок и частиц газа ${{X}_{{\text{g}}}} = {{N}_{{\text{g}}}}\mu $, где $\mu = 2.3$ а.е. м. – средняя масса частиц газа. На первом и втором этапах эволюции пылинок величина ${{N}_{{\text{g}}}}$ не меняется, но на этапе 3 монотонно уменьшается:

Выражения (7)–(10) протестированы на модели свободного коллапса однородного шарообразного облака из состояния покоя. Движение любой частицы на краю такого облака есть свободное падение в поле с потенциалом $ - GM{\text{/}}R$, где $M$ – масса облака, $R = R(t)$ – его радиус, $G$ – гравитационная постоянная. Поэтому скорость частицы $dR{\text{/}}dt = - \sqrt {2GM({{R}^{{ - 1}}} - R_{0}^{{ - 1}})} $, где ${{R}_{0}} = R(0)$. Переменные $t$ и $R$ в этом уравнении разделяются и

Свободный коллапс заканчивается при $R = 0$ и $t = {{t}_{{{\text{ff}}}}}$, где здесь ${{t}_{{{\text{ff}}}}}$ – время свободного падения, ${{\rho }_{0}}$ – начальная плотность облака. После выражения $R$ через безразмерную плотность $\psi \equiv \rho {\text{/}}{{\rho }_{0}}$ решение (11) записывается в виде Температура считается кусочно-степеннóй функцией плотности (см. раздел 3.1). Зависимость интегральных характеристик пыли от температуры при ${{Y}_{{\text{d}}}} = 0.01$ показана на рис. 3. Испарение водяного льда становится заметным при $T \approx 111$ К и заканчивается при $T \approx 136$ К. Это приводит к почти скачкообразному уменьшению $\langle a\rangle $, $\langle s\rangle $ и $\langle m\rangle $ в 2, 4 и 3.8 раза соответственно. У однослойных пылинок масса уменьшилась бы в 8 раз. На этом этапе $\max [{{Y}_{{{\text{ev}}}}}] = 0.0074$, ${{X}_{{\text{g}}}} = {\text{const}}$. Испарение ядер становится заметным при $T \approx 1150$ К и заканчивается при $T \approx 1490$ К. На этом этапе $\max [{{Y}_{{{\text{ev}}}}}] = 0.01$, а ${{X}_{{\text{g}}}}$ плавно уменьшается до нуля. Значения $\langle a\rangle $, $\langle s\rangle $ и $\langle m\rangle $ вблизи $T \approx 1455$ К достигают максимумов, что вызвано уменьшением количества мелких пылинок. Затем эти величины резко спадают до нуля, поскольку исчезают все пылинки. Максимумы у $\langle a\rangle $, $\langle s\rangle $ и $\langle m\rangle $ соответственно в 5, 26 и 106 раз выше начальных значений, однако это не проявляется во взаимодействии пыли с частицами газа, поскольку такое повышение компенсируется спадом ${{X}_{{\text{g}}}}$ в 11 тысяч раз.

3.1. Модель коллапса протозвездного облака

Гравитационный коллапс протозвездного облака моделируется численно в приближении двухкомпонентной магнитной газодинамики. Межзвездная плазма рассматривается как сплошная среда, которая является смесью 1) электрически нейтральных частиц (нейтралов) и 2) электропроводящего компонента, состоящего из электронов, ионов и заряженных пылинок. Исследуются стадии коллапса до температуры существенной тепловой ионизации водорода (10 тыс. К), поэтому считается, что плотность и тепловое давление среды определяются только концентрацией нейтралов ${{n}_{{\text{n}}}}$: $\rho = \mu {{n}_{{\text{n}}}}$, $P = {{n}_{{\text{n}}}}{{k}_{{\text{b}}}}T$, где ${{k}_{{\text{b}}}}$ – постоянная Больцмана. Скорость среды ${\vec {v}}$ приравнивается к скорости нейтралов. Поскольку в численном методе используется лагранжева разностная сетка (движущаяся со скоростью среды), то уравнение движения единицы объема среды записывается через субстанциональную производную $d{\text{/}}dt$:

где $\vec {g}$ – ускорение свободного падения. Используется система единиц СГСЭ.

При описании диффузии магнитного поля в коллапсирующих протозвездных облаках можно пренебречь трением проводящих компонентов плазмы друг с другом [46]. Поэтому предполагается, что компонент сорта $\alpha $, состоящий из частиц с массой ${{m}_{\alpha }}$, зарядом $ \pm e$ и концентрацией ${{n}_{\alpha }}$, имеет проводимость ${{\sigma }_{\alpha }} = {{e}^{2}}n_{\alpha }^{2}{\text{/}}{{R}_{{\alpha {\text{n}}}}}$, где ${{R}_{{\alpha {\text{n}}}}} = {{\mu }_{{\alpha {\text{n}}}}}{{n}_{\alpha }}{{n}_{{\text{n}}}}{{\langle sv\rangle }_{{\alpha {\text{n}}}}}$ – коэффициент трения с нейтралами, ${{\mu }_{{\alpha {\text{n}}}}} = \mu {{m}_{\alpha }}{\text{/}}(\mu + {{m}_{\alpha }})$ – приведенная масса, ${{\langle sv\rangle }_{{\alpha {\text{n}}}}}$ – коэффициент замедления. Для ионов ${{\langle sv\rangle }_{{{\text{in}}}}} = 2 \times {{10}^{{ - 9}}}$ см³ с^–1 (см. [47, п. 2.2]). Для заряженных пылинок ${{\langle sv\rangle }_{{{{{\text{g}}}_{ + }}{\text{n}}}}} = {{\langle sv\rangle }_{{{{{\text{g}}}_{ - }}{\text{n}}}}} = \langle s\rangle {{v}_{{\text{n}}}}$, где ${{v}_{{\text{n}}}} = \sqrt {8{{k}_{{\text{b}}}}T{\text{/}}\pi \mu } $ – средняя тепловая скорость нейтралов. Коэффициент замедления электронов задается кусочно-степеннóй аппроксимацией табличных данных из [48, с. 393],

где ${{E}_{{\text{e}}}} = 6.25 \times {{10}^{{11}}}{{m}_{{\text{e}}}}v_{{\text{e}}}^{2}{\text{/}}2$ – средняя тепловая энергия электронов в эВ, ${{v}_{{\text{e}}}} = \sqrt {8{{k}_{{\text{b}}}}T{\text{/}}\pi {{m}_{{\text{e}}}}} $ – их средняя тепловая скорость, ${{m}_{{\text{e}}}}$ – масса электрона.

Проводимость плазмы $\sigma = \sum \,{{\sigma }_{\alpha }}$, т.е. циклотронное вращение заряженных частиц не учитывается. Коэффициент трения проводящего компонента с нейтралами ${{R}_{{{\text{in}}}}} = \sum\nolimits_\alpha \,{{R}_{{\alpha {\text{n}}}}}$. Магнитная амбиполярная диффузия (дрейф плазмы) – медленное совместное движение всех заряженных частиц относительно нейтралов под действием электромагнитной силы – считается стационарной (см. [47, п. 13.3]), т.е. полагается, что электромагнитная сила уравновешена силой трения

где ${{\vec {v}}_{{\text{i}}}}$ – скорость проводящего компонента. Уравнение (13) является алгебраическим уравнением движения этого компонента.

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

где ${{\eta }_{{\text{O}}}} = {{c}^{2}}{\text{/}}4\pi \sigma $ – коэффициент омической диффузии, $c$ – скорость света в вакууме. Выразив ${{{\vec {v}}}_{{\text{i}}}}$ из (13) и подставив в (14), можно получить уравнение индукции, в котором магнитная амбиполярная диффузия порождает диффузию магнитного поля с коэффициентом ${{\eta }_{{\text{A}}}} = {{B}^{2}}{\text{/}}4\pi {{R}_{{{\text{in}}}}}$. Однако в численной модели эта подстановка не используется и решается уравнение (14).

Турбулентная диффузия магнитного поля не рассматривается, поскольку магнитная амбиполярная диффузия подавляет МГД турбулентность еще вначале коллапса протозвездного облака [25, 26].

В численной модели вращение облака не учитывается, поле скорости и поле плотности считаются сферически симметричными, а магнитное поле – полоидальным. В сферических координатах $(r,\theta ,\phi )$ компоненты магнитного поля представляются в виде ${{B}_{r}} = {{H}_{r}}(r,t)\cos \theta $, ${{B}_{\theta }} = - {{H}_{\theta }}(r,t)\sin \theta $. Это позволяет усреднить по полусферам и электромагнитную силу в (12), (13), и уравнение индукции (14), чтобы получить сферически симметричную модель коллапса [49]. В такой модели $g = GM{\text{/}}{{r}^{2}}$, где $M = M(r)$ – масса внутри сферы с радиусом $r$. На лагранжевой сетке масса $dM$ сферического слоя с объемом $dV$ постоянна, поэтому $\rho = dM{\text{/}}dV$, а $dV$ вычисляется через координаты узлов сетки.

Температура задается кусочно-степеннóй функцией плотности, аппроксимирующей результаты работы [50]. Показатель степени меняется в интервале (0.1, 0.66) и зависит не от плотности, а от самой температуры, поскольку с ее ростом возбуждаются колебательные и вращательные степени свободы молекул водорода, а затем происходит их диссоциация.

3.2. Модель ионизации

Для корректного расчета диффузии магнитного поля в протозвездных облаках необходимо знать концентрации электронов ${{n}_{{\text{e}}}}$, ионов ${{n}_{{\text{i}}}}$ и пылинок ${{n}_{{{{{\text{g}}}_{ - }}}}}$ и ${{n}_{{{{{\text{g}}}_{ + }}}}}$, имеющих заряд $ - e$ и $ + e$ соответственно (см. [38, 46]). Уравнения кинетики ионизации разделены на две системы: первая описывает ионизацию среды космическими лучами и радионуклидами, вторая – ее ионизацию собственным тепловым излучением. В обеих системах учитывается рекомбинация ионов на пылинках.

Считается, что при $T < 575$ K в газовой фазе есть только один вид атомов металлов – Mg, а остальные десорбируются из ядер пылинок при $T > 575$ K – после испарения мало летучей органики [51]. Предполагается, что эти атомы металлов равномерно распределены в верхней трети ядра (по радиусу), и в интервале температур от 575 до 1300 K занимаемая ими область линейно уменьшается по радиусу с ростом $T$. Распространенность таких атомов в газовой фазе пропорциональна изменению занимаемого ими объема пылинок и для каждого элемента ограничена его распространенностью в современной солнечной системе [40]. Такая модель десорбции металлов из пылинок реалистична, но строго не обоснована. Она лишь призвана описать плавный впрыск металлов в газовую фазу, не повышая сложность модели испарения пыли.

Поскольку моделируется коллапс облаков с начальной концентрацией частиц $n \simeq {{10}^{5}}{\kern 1pt} - {\kern 1pt} {{10}^{6}}$ см^–3, то пренебрегается ионизацией среды межзвездным ультрафиолетом, и ион C⁺ не учитывается, хотя он значим во внешних слоях таких облаков (см. [52]). При наличии мелких (5–10 нм) пылинок на начальной стадии коллапса с $T \approx 10$ K доминирует ион ${\text{H}}_{3}^{ + }$ [53], а в их отсутствие – ион HCO⁺ [52, 53]. На более теплой стадии коллапса начинается сублимация летучих веществ из мантий пылинок, поэтому доминируют ионы HCO⁺ и H₃O⁺, независимо от наличия мелких пылинок [53]. Учитывая мелкую пыль, мы рассматриваем только один молекулярный ион – ${\text{H}}_{3}^{ + }$. Пренебрежение ионами HCO⁺ и H₃O⁺ влияет на остаточное поле незначительно, поскольку в зонах сильной диффузии магнитного поля проводимость определяется именно мелкой пылью [53].

Предполагается, что ионы Mg⁺ и других металлов образуются в процессе диссоциативной перезарядки иона ${\text{H}}_{3}^{ + }$ с атомами, а также при их фотоионизации тепловым излучением среды на поздних стадиях коллапса. Ион ${\text{H}}_{3}^{ + }$ возникает в результате быстрых реакций, начинающихся с ионизации водорода и гелия космическими лучами и радионуклидами, поэтому темп производства H$_{3}^{ + }$ в единице объема

где $\zeta $ – вероятность ионизации космическими лучами одной частицы за секунду на границе облака с радиусом $R$, $\tau $ – показатель ослабления темпа ионизации в сферическом слое с радиусом $r$, ${{\chi }_{\zeta }} \approx 100$ г см^–2 – столбцовая плотность, при которой темп ионизации ослабевает в $e$ раз. Измерения потока космических лучей космическими аппаратами Pioneer и Voyager приводят к оценке [54] $\zeta \approx 3 \times {{10}^{{ - 17}}}$ с^–1. Выражение (15) также учитывает ионизацию среды двумя радионуклидами [55]: $^{{60}}$Fe, ионизирующий одну частицу с вероятностью ${{\zeta }_{1}} = 8 \times {{10}^{{ - 20}}}$ с^–1 и обладающий периодом полураспада ${{t}_{1}} = 1.5 \times {{10}^{6}}$ лет, а также ²⁶Al, у которого, соответственно, ${{\zeta }_{2}} = 6 \times $ $ \times {{10}^{{ - 19}}}$ с^–1 и ${{t}_{2}} = 7 \times {{10}^{5}}$ лет.

Первая система уравнений кинетики ионизации основана на модели Кунца и Мусковиаса [56] и дополнена учетом рекомбинации ионов на пылинках при $T > 575$ K. Эта система состоит из 4 уравнений для расчета ${{n}_{{{{{\text{g}}}_{ - }}}}}$ и ${{n}_{{{{{\text{g}}}_{ + }}}}}$, а также для расчета концентраций ${{n}_{{{{{\text{a}}}_{ + }}}}}$ и ${{n}_{{{{{\text{m}}}_{ + }}}}}$ ионов Mg⁺ и ${\text{H}}_{3}^{ + }$:

где ${{n}_{{{{{\text{a}}}_{0}}}}} = {{n}_{{\text{a}}}} - {{n}_{{{{{\text{a}}}_{ + }}}}}$ – концентрация нейтральных атомов Mg в газовой фазе, ${{n}_{{\text{a}}}}$ – полная концентрация Mg в газовой фазе, ${{n}_{{{{{\text{g}}}_{0}}}}} = {{n}_{{\text{g}}}} - {{n}_{{{{{\text{g}}}_{ + }}}}} - {{n}_{{{{{\text{g}}}_{ - }}}}}$ – концентрация нейтральных пылинок, ${{n}_{{\text{g}}}}$ – полная концентрация пылинок. Индекс $\beta $ пробегает два значения: a (атом), m (молекула). Обозначения констант реакций расшифрованы в табл. 1. Коэффициент ${{C}_{{{{{\text{a}}}_{ + }}{\text{e}}}}}$ аппроксимируется соотношением для водородоподобного иона (см. [47, п. 5.1]).

Вторая система уравнений кинетики ионизации состоит из 10 уравнений для расчета концентраций однозарядных ионов 10 химических элементов: H, Li, Na, Mg, Al, Si, K, Ca, Rb, Cs. В этой системе уравнений в расчет приняты только те атомы Mg, которые десорбированы пылью при $T > 575$ K. Каждое из уравнений второй системы решается только при $T > 575$ K и имеет вид

где индекс $\beta $ означает химический элемент, есть вероятность фотоионизации атома сорта $\beta $ за секунду (см. [47, п. 5.2]) собственным тепловым излучением среды, ${{u}_{\nu }} = 8\pi h{{(\nu {\text{/}}c)}^{3}}{{({{e}^{{h\nu /{{k}_{{\text{b}}}}T}}} - 1)}^{{ - 1}}}$ – его плотность энергии, $\nu $ – частота, ${{\nu }_{\beta }} = {{I}_{\beta }}{\text{/}}h$ – пороговая частота ионизации, ${{I}_{\beta }}$ – энергия ионизации, ${{s}_{\beta }}(\nu )$ – сечение фотоионизации, $h$ – постоянная Планка, ${{C}_{{{{\beta }_{0}}{{{\text{m}}}_{ + }}}}}$ – константа диссоциативной перезарядка H$_{3}^{ + }$ с атомом, ${{C}_{{{{\beta }_{ + }}{\text{e}}}}}$ – константа фоторекомбинации (см. табл. 1).

Поскольку обилие металлов мало, считается, что их фотоионизация не влияет на поле излучения. Зависимости ${{s}_{\beta }}(\nu )$ для разных атомов сильно различаются, однако при исследовании магнитной диффузии достаточно знать пороговые значения ${{s}_{\beta }}({{\nu }_{\beta }})$. Дело в том, что в области эффективной магнитной диффузии $T < 5000$ К и $h\nu {\text{/}}{{k}_{{\text{b}}}}T \gg 1$, поэтому с ростом частоты функция ${{u}_{\nu }}$ настолько быстро убывает, что для всех зависимостей ${{s}_{\beta }}(\nu )$ справедлива следующая аппроксимация интеграла (21) (см. [58]):

где $\varepsilon = {{I}_{\beta }}{\text{/}}{{k}_{{\text{b}}}}T$. Пороговые значения ${{s}_{\beta }}({{\nu }_{\beta }})$ берутся из таблиц [48, 59].

Концентрация электронов находится из условия электронейтральности,

где сумма берется по всем ионам, существующим при данной температуре.

3.3. Начальные и граничные условия, численные методы

Начальная температура среды ${{T}_{0}}$ постоянна по всему облаку, а распределение плотности вещества и модуля магнитной индукции вычисляется из условия гидростатического равновесия с постоянным отношением магнитного и теплового давлений. Это отношение вычисляется через концентрацию частиц ${{n}_{{{\text{c0}}}}}$ и магнитную индукцию ${{B}_{{{\text{c0}}}}}$ в центре облака (табл. 2). Среда вначале покоится. Поскольку магнитное натяжение в расчете равновесия не учитывается, то начальная равнодействующая сила в каждом слое облака отлична от нуля. Однако она существенно слабее силы тяготения, поэтому скорость коллапса меньше, чем в случае изначально однородного облака.

Начальные характеристики протозвездных облаков (табл. 2) вычисляются на основе наблюдательных данных о столбцовой плотности, радиусе облака и магнитной индукции вдоль луча зрения [10]. Последняя умножается на $\sqrt 2 $, чтобы получить ${{B}_{{{\text{c0}}}}}$. Значение ${{n}_{{{\text{c0}}}}}$ подбирается так, чтобы радиус и столбцовая плотность гидростатически равновесного облака совпали с наблюдаемыми значениями. После этого вычисляется масса облака. Начальные отношения тепловой и магнитной энергий облака к модулю его гравитационной энергии, ${{\varepsilon }_{{\text{t}}}}$ и ${{\varepsilon }_{{\text{m}}}}$, в сумме $ < {\kern 1pt} 1$, поэтому устойчивое гидростатическое равновесие облака недостижимо, и оно, иногда после некоторых колебаний, испытывает гравитационный коллапс.

Граничные условия в центре облака: скорость меняет знак, плотность и магнитное поле непрерывны. На внешней границе облака: скорость среды вычисляется из условия постоянства массы и плотности внешнего сферического слоя, производная скорости проводящего компонента непрерывна, магнитная индукция непрерывна.

Задача с начальными и граничными условиями для одномерных уравнений движения и индукции решается численно методом Лакса-Вендроффа с добавленной искусственной вязкостью. Этот метод достаточно прост в реализации, но, будучи явным, требует слишком малых шагов по времени при сильной диффузии магнитного поля. Основа численного кода разработана Дудоровым и Сазоновым [49] и существенно доработана авторами для включения современных моделей ионизации и испарения пыли.

Для исследования влияния пыли на эволюцию магнитного поля протозвездных облаков рассматривается два варианта моделей, отличающихся параметрами пыли и начальным обилием Mg в газовой фазе ${{X}_{{\text{a}}}}(0)$: 1) мелкая пыль с нормальной массовой долей ${{Y}_{{\text{d}}}}$, приводящая к низкому значению ${{X}_{{\text{a}}}}(0)$, 2) крупная пыль с низкой массовой долей, сохраняющая высокое значение ${{X}_{{\text{a}}}}(0)$ (табл. 3). Такая антикорреляция ${{Y}_{{\text{d}}}}$ и ${{X}_{{\text{a}}}}(0)$ описывает на качественном уровне так называемое вымерзание атомов на пылинках, что означает преобладание их адсорбции над десорбцией.

Для прояснения роли пыли рассматривается также нереальный, но поучительный вариант 1а без испарения пыли, имеющий те же параметры, что и вариант 1. Десорбция металлов из пылинок при $T > 575$ K в варианте 1а учитывается, чтобы не занижать степень ионизации слишком сильно.

Уравнения (16)–(20) решаются в каждом сферическом слое одномерной модели облака методом Гира 2 порядка точности. Возникающие при этом системы нелинейных уравнений решаются методом Ньютона, а системы линейных уравнений – методом Гаусса.

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

3.4.1. Концентрация заряженных частиц. Рисунок 4 показывает зависимости температуры и относительной концентрации заряженных частиц от концентрации нейтральных частиц в центре облака в трех вариантах модели W3 (main). Пылинки являются доминирующими носителями заряда при $2 \times {{10}^{8}} < {{n}_{{\text{n}}}} < 7 \times {{10}^{{13}}}$ см^–3 в варианте 1 и при $2 \times {{10}^{8}} < {{n}_{{\text{n}}}} < 3 \times {{10}^{{15}}}$ см^–3 в варианте 1а. В этих случаях ион ${\text{H}}_{3}^{ + }$ доминирует только при ${{n}_{{\text{n}}}}{{ < 10}^{8}}$ см^–3. В варианте 2 почти все время доминируют ионы металлов и электроны; лишь при $T \approx 8000$ K их сменяют протоны с электронами; концентрация иона H$_{3}^{ + }$ всюду уступает ионам металлов на порядок величины и более.

В варианте 1 испарение ледяных мантий при $T \approx 120$ K вызывает повышение концентраций электронов в 2 раза, а ионов Mg⁺ в 54 раза, поскольку из-за уменьшения площади поперечного сечения пылинок (в 4 раза) происходит снижение темпа рекомбинации на пыли у доминирующего иона H$_{3}^{ + }$. При $T > 575$ K относительная концентрация иона H$_{3}^{ + }$ быстро снижается из-за диссоциативной перезарядки с атомами металлов, десорбированных из пылинок. При $T \approx 1350$ K пылинки исчезают совсем. Тогда концентрация электронов возрастает на 7 порядков, а ионов Mg⁺ – на 5 порядков величины.

В варианте 1а резких изменений концентраций всего два: при $T > 575$ K из-за десорбции металлов быстро падает концентрация H$_{3}^{ + }$ и возрастает концентрация Mg⁺, а при $T > 1300$ K быстро растет относительная концентрация электронов из-за фотоионизации металлов. Резкий спад относительной концентрации Mg⁺ при $T \approx 1500{\kern 1pt} - {\kern 1pt} 1800$ K в вариантах 1 и 2 вызван усилением рекомбинации из-за быстрого роста концентрации электронов при фотоионизации щелочных металлов. При более высоких температурах начинается интенсивная фотоионизация самих атомов Mg.

Зависимость относительной концентрации протонов и ионов металлов от концентрации нейтральных частиц и температуры в центре облака в варианте 1 модели W3 (main) показана на рис. 5. С началом десорбции металлов из пылинок при $T > 575$ K основным конкурентом Mg⁺ становится Si⁺, поскольку Mg и Si близки по распространенности и кинетическим коэффициентам. Концентрация остальных ионов металлов при $T < 1100$ K также соответствует распространенности. Тепловая ионизация сначала заметна у металлов с наименьшей энергией ионизации (Cs, Rb, K при $T \approx 1000$, 1050, 1150 K соответственно), но выйти в лидеры им не удается из-за низкой распространенности. При $T < 1400$ K за счет высокой распространенности доминируют Mg⁺ и Si⁺, а затем их сменяют K⁺, Na⁺ и Al⁺. Таким образом, в используемой модели можно пренебречь ионизацией Ca, Li, Rb, Cs.

В моделях объектов NGC 2024 и DR 21 OH1 относительные концентрации заряженных частиц зависят от плотности и температуры почти так же, как в модели W3 (main), поскольку и температуры, и ослабление потока космических лучей во всех моделях различаются не сильно.

3.4.2. Коэффициенты диффузии и магнитное число Рейнольдса. Диффузия магнитного поля сильна, если магнитное число Рейнольдса ${{{\text{R}}}_{{\text{m}}}} < 1$. Выбрав радиальную координату $r$ в качестве характерного размера течения, положим ${{{\text{R}}}_{{\text{m}}}} = vr{\text{/}}({{\eta }_{{\text{A}}}} + {{\eta }_{{\text{O}}}})$, где, напомним, ${{\eta }_{{\text{A}}}}$ и ${{\eta }_{{\text{O}}}}$ – коэффициенты амбиполярной и омической диффузии. Магнитная амбиполярная диффузия максимальна там, где минимальна концентрация массивных заряженных частиц с большим сечением рассеяния на нейтралах, т.е. ионов и пылинок. Омическая диффузия магнитного поля максимальна там, где минимальна относительная концентрация легких заряженных частиц, т.е. электронов и ионов. Поэтому зоны эффективной амбиполярной и омической диффузий обычно не совпадают, а испарение пыли существенно влияет на их соотношение.

Рассмотрим финальные зависимости ${{\eta }_{{\text{A}}}}(r)$ и ${{\eta }_{{\text{O}}}}(r)$ в трех вариантах модели W3 (main) (рис. 6). Эти зависимости немонотонны и имеют пики и провалы в зонах испарения мантий пылинок (30–80 а.е.) и ядер пылинок (3–6 а.е.). Во внешних частях облака ($r > 100$ а.е.) во всех вариантах ${{\eta }_{{\text{A}}}} > {{\eta }_{{\text{O}}}}$. Если же пыли мало и она крупная (вариант 2), то ${{\eta }_{{\text{A}}}} > {{\eta }_{{\text{O}}}}$ при $r > 0.3$ а.е., т.е. амбиполярная диффузия доминирует в средних и внешних слоях облака. В варианте 1 глобальные максимумы ${{\eta }_{{\text{A}}}}$ и ${{\eta }_{{\text{O}}}}$ одинаковы по порядку величины, причем в варианте 1а это совпадение сохраняется. Зона высоких значений ${{\eta }_{{\text{A}}}}$ наиболее широка в варианте 2, поскольку сцепление проводящего компонента с нейтралами ослаблено из-за дефицита заряженных пылинок. Зона высоких значений ${{\eta }_{{\text{O}}}}$ наиболее широка в варианте 1а, поскольку концентрации электронов и ионов снижены из-за неослабевающей рекомбинации на пылинках.

Зависимость магнитного числа Рейнольдса от радиальной координаты в трех вариантах модели W3 (main) на момент завершения расчетов также представлена на рис. 6. Диффузия магнитного поля считается сильной при ${{R}_{{\text{m}}}} < 1$, умеренной при $1 < {{R}_{{\text{m}}}} < 10$ и слабой при ${{R}_{{\text{m}}}} > 10$. Зона сильной диффузии имеет ширину 3 а.е. в варианте 1 и ширину 22 а.е. в варианте 1а, но отсутствует в варианте 2. Зон умеренной диффузии может быть несколько: две в вариантах 1 и 1а, но три в варианте 2. Пренебрежение испарением пыли расширило зону сильной диффузии в 7 раз и умеренной диффузии в 2 раза.

Поведение коэффициентов диффузии и магнитного числа Рейнольдса в моделях NGC 2024 и DR 21 OH1 на качественном уровне такое же, как в модели W3 (main), однако положения и величины экстремумов отличаются. В табл. 4 сравниваются основные показатели моделей: минимальная относительная концентрация электронов в центре облака и соответствующая ей температура $T{\kern 1pt} '$, а также магнитная индукция в центре облака при плотности фотосфер молодых звезд. Кроме того, в табл. 4 указаны границы зоны сильной диффузии магнитного поля (мертвой зоны), в которой должны быть существенно снижены турбулентная вязкость и радиальная диффузия дециметровых твердых тел [60]. Такие тела могут собираться в планетезимали на границах мертвой зоны из-за неустойчивости волн Россби [61], а также внутри этой зоны из-за потоковой неустойчивости [60].

Таблица 4 показывает, что минимальные относительные концентрации электронов в центре облака в разных моделях и вариантах отличаются на 5.5 порядков величины, но значения $T{\kern 1pt} '$ различаются всего на $260$ К, поскольку температуры полного испарения пыли и существенной тепловой ионизации металлов тоже близки. В моделях Мершана и др. [38], учитывающих и не учитывающих испарение пыли, минимальные относительные концентрации электронов отличаются на порядок величины, но расположение минимума также почти не меняется. Таким образом, пренебрежение испарением пыли в моделях коллапса протозвездных облаков существенно усиливает диффузию магнитного поля, но сохраняет характерное расстояние от центра до мертвой зоны, 1–10 а.е.

3.4.3. Эволюция магнитного поля. В используемой сферически симметричной модели коллапса [49] сплющивание облака вдоль магнитного поля не учитывается, поэтому зависимость магнитной индукции от плотности ближе к $B \propto {{\rho }^{{2/3}}}$, чем в многомерных моделях. Рассмотрим, например, зависимость магнитной индукции ${{B}_{{\text{c}}}}$ от концентрации частиц ${{n}_{{\text{c}}}}$ в центре облака в модели DR 21 OH1 (рис. 7). В нашей модели амбиполярная диффузия в центре облака неэффективна, поэтому с высокой точностью ${{B}_{{\text{c}}}} \propto n_{{\text{c}}}^{{2/3}}$ при ${{n}_{{\text{c}}}}{{ < 10}^{{12}}}$ см^–3 и ${{n}_{{\text{c}}}}{{ > 10}^{{15}}}$ см^–3, однако ${{B}_{{\text{c}}}} \propto n_{{\text{c}}}^{{0.45}}$ в условиях сильной омической диффузии при ${{10}^{{13}}} < {{n}_{{\text{c}}}}{{ < 10}^{{14}}}$ см^–3. Пренебрежение испарением пыли (вариант 1а) продлевает интервал сильной омической диффузии до ${{n}_{{\text{c}}}} \simeq {{10}^{{15}}}$ см^–3, что снижает магнитную индукцию на полпорядка величины. И наоборот, в случае дефицита пыли (вариант 2) обилие электронов достаточно велико, чтобы омическая диффузия была несущественной и лишь чуть понизила ${{B}_{{\text{c}}}}$. В итоге в диапазоне ${{n}_{{\text{c}}}}{{ = 10}^{{16}}}{\kern 1pt} - {\kern 1pt} {{10}^{{18}}}$ см^–3, соответствующем основаниям фотосфер молодых звезд, ${{B}_{{\text{c}}}} = 0.7{\kern 1pt} - {\kern 1pt} 14.5$ кГс в варианте 1, $0.2{\kern 1pt} - {\kern 1pt} 5.2$ кГс в варианте 1а, и $1.4{\kern 1pt} - {\kern 1pt} 28.3$ кГс в варианте 2. Таким образом, магнитное поле при фотосферных плотностях ослабло в 3 раза из-за пренебрежения испарением пыли и выросло в 2 раза из-за дефицита пыли.

Значения магнитной индукции ${{B}_{{\text{c}}}}$ в центре облака при плотности фотосфер молодых звезд для всех моделей и вариантов собраны в последнем столбце табл. 4. Видно, что ${{B}_{{\text{c}}}}$ варьируется от десятков гаусс до десятков килогаусс и коррелирует с начальным отношением магнитной энергии облака к модулю его гравитационной энергии. В модели W3 (main) поле занижается в 2.6 раза в отсутствие испарения пыли и завышается в 1.5 раза из-за дефицита пыли. В модели NGC 2024, соответственно, в 5 раз и в 4 раза.

Полагая варианты 1 в трех моделях наиболее корректными, получим при плотности фотосфер ${{B}_{{\text{c}}}} = 0.2{\kern 1pt} - {\kern 1pt} 49.7$ кГс, что на полтора порядка величины больше наблюдаемых магнитных полей молодых звезд. Эту нестыковку можно устранить с помощью многомерных моделей коллапса протозвездных облаков, если начинать расчет с достаточно сильным магнитным полем, которое вызывает значительное уплощение облака вдоль поля. Например, Замоздра и Жилкин [37] показали, что переход от одномерных к двумерным моделям привел к снижению ${{B}_{{\text{c}}}}$ на 2 порядка величины в модели W3 (main), на 1.5 порядка величины в модели DR 21 OH1 и на порядок величины в модели NGC 2024, еще до прохождения минимума степени ионизации. В других работах такое сравнение одномерных и двумерных моделей коллапса конкретных протозвездных облаков не выполнялось. Однако более медленный, чем (1), рост магнитного поля с плотностью на начальных стадиях коллапса показывают все многомерные модели, в которых магнитное поле достаточно сильно, чтобы замедлить поперечное сжатие (см., напр., [21, 29–32]).

Поскольку в работах по магнитометрии звезд и околозвездных дисков обычно указывается не плотность вещества, а геометрические характеристики, используем их тоже для проверки теории. На рис. 8 представлены зависимости магнитной индукции от радиальной координаты $B(r)$ в варианте 1 трех моделей, а также в классических звездах типа T Tau (CTTS) [14], аккреционном диске FU Ori [15], оболочках дисков TW Hya [17] и AS 209 [18]. Также показана оценка магнитной индукции в протосолнечной туманности на расстоянии 2 а.е. от Солнца по остаточной намагниченности хондр в LL хондрите Semarkona [16]. Как и в работе [37], ближе всего к точке FU Ori проходит кривая $B(r)$ в модели NGC 2024 с самым слабым начальным полем и самой сильной диффузией, но для прохода через точки CTTS этой кривой необходимо спуститься еще на 1.5 порядка величины. Кривым $B(r)$ в моделях DR 21 OH1 и W3 (main) необходимо спуститься еще ниже: на 2 и 2.5 порядка величины соответственно.

Вероятно, эти различия теории и наблюдений можно почти полностью устранить, перейдя от одномерных к многомерным моделям. Если это снизит остаточное поле на 1–2 порядка величины, как в работе [37], то к данным FU Ori и Semarkona ближе окажется модель DR 21 OH1, а к данным TW Hya и AS 209 – модель NGC 2024.

Заметим, что в [37] не учитывалось химическое распыление графита, поэтому считалось, что он исчезает при $T = 2300$ К. Это приводило к более медленному росту магнитного поля при сжатии среды на расстоянии около 1 а.е. от центра облака. Поэтому на меньших расстояниях кривые $B(r)$ проходили ниже примерно на 50%, чем на рис. 8.

Обратим внимание, что кривые $B(r)$ на рис. 8 имеют локальные минимумы в мертвой зоне ($r \simeq 1{\kern 1pt} - {\kern 1pt} 10$ а.е.), что указывает на формирование магнитного уплотнения (локального максимума) на внешней границе этой зоны. В моделях без испарения пыли это уплотнение более выражено. Вероятно, в трехмерных моделях, учитывающих вращение облака, магнитное уплотнение превращается в магнитную трубку, обнаруженную Вастером и др. [31]. Однако мы предполагаем, что при учете испарения пыли эта трубка не должна быть слишком заметной. Тем не менее такая структура интересна и требует дальнейшего изучения. В отличие от “магнитной стены” из двумерной модели Тассиса и Мусковиаса [62] она создает немонотонную зависимость интенсивности магнитного поля от радиальной координаты.

Наконец отметим, что возможность ослабления магнитного поля в центре облака с ростом плотности [31] наши расчеты не подтверждают, даже в вариантах без испарения пыли. Вероятно, этот эффект вызван ошибками реализации метода сглаженных частиц (SPH) в модели [31], оставшимися от ее предыдущего варианта, в котором не учитывалась диффузия магнитного поля [21]. К такому же выводу приводит сравнение с конкурирующими моделями [30, 32]: в них магнитная индукция в центре облака монотонно растет с плотностью, хотя модель [30] основана на SPH, а модель [32] использует сеточные методы.