<< 3. Численное моделиpование звездных 5. Замечания о пpиpоде >>

Subsections

• 4.1 Осесимметpичная дисковая аккpеция
• 4.2 Hеосесимметpичная дисковая аккpеция
• 4.3 Hеустойчивости в аккpеционных дисках

4. Аккpеционные диски

Термин аккреционный диск (АД) обычно употребляется для обозначения газового диска вокруг массивного (по сравнению с диском) компактного объекта. К последним относятся белые карлики, нейтронные звезды, черные дыры. Определяющей чертой АД является переход гравитационной энергии при аккреции (падении) вещества в тепло с последующим излучением. Весьма часто понятием ``аккреционный диск'' пользуются в более широком смысле, имея в виду газовый диск, для которого существен процесс падения вещества на аккумулирующий центр, например, когда речь идет о протозвездных, протопланетных, галактических дисках.

Таким образом, с точки зрения физических процессов любые результаты, относящиеся к газовым дискам, в полной мере относятся к АД. Поэтому различие между главами 4 и 5 весьма условно и связано в основном с видом внешнего потенциала. Для дисков, вращающихся вокруг компактных объектов массой $M_1$, он является ньютоновым для точечного тела ( $\Phi \propto - M_1/r$), и самогравитация газа несущественна, за исключением, пожалуй, дисков вокруг сверхмассивных черных дыр. В случае галактических газовых дисков внешний потенциал определяется распределением звездного вещества (см. (2.44)), и самогравитация часто оказывается весьма важным фактором.

Одной из ключевых проблем физики аккреционных дисков является вопрос о механизмах отвода углового момента, обеспечивающих падение вещества на гравитирующий центр. Достигнутый прогресс в немалой степени основан на вязких моделях аккреционных дисков (АД) [ 252, 253, 280, 304, 689, 790, 814, 842 и др.]. По-видимому, первым рассмотрел механизм переноса момента импульса из-за действия вязкости еще C.F.$\,$von Weizsacker в 1948 г. [892]. Другой возможностью является отвод углового момента посредством крупномасштабных спиральных ударных волн, которые могут формироваться при перетекании газа в тесных двойных системах (ТДС) с учетом гравитационного влияния звезды-донора [16, 302, 829, 830].

Вязкие модели АД, в основе которых лежит гипотеза о турбулентной вязкости ($\alpha$-модель Шакуры и Сюняева [221, 790] и ее многочисленные модификации), позволяют объяснить многие наблюдения в ТДС, активных галактических ядрах (АГЯ), протозвездных дисках, однако удовлетворительной нефеноменологической модели турбулентной вязкости в настоящее время нет. Вязкость может играть важную роль в газовых галактических дисках [34, 184]. Следует отметить, что, несмотря на достижения последних лет [271, 304, 653, 696, 845, 857], активное изучение турбулентности в аккрецирующих системах только начинается и основные результаты еще впереди.

Имеется много работ, посвященных линейному анализу устойчивости вязких АД [36, 149, 205, 209, 230, 245, 270, 271, 303, 335, 538, 580, 617, 643, 791, 884, 885, 842 и др.]. Однако вопрос о последствиях развития неустойчивостей на существенно нелинейных стадиях только начинает рассматриваться в рамках численного гидродинамического моделирования [212, 233, 237, 272, 510, 647, 857].

Проблема конвекции в газовых дисках восходит к работам Вейцзеккера [892] и Сафронова [183]. Развитие ADAF-моделей (адвективно доминирующие аккреционные течения) для объяснения низкой светимости рентгеновских двойных и активных галактических ядер с черными дырами [231, 508, 659] привело к построению CDAF-моделей (конвективно-доминирующие аккреционные течения), в которых при слабой вязкости формируются течения с сильной турбулентностью из-за развития конвективной неустойчивости [272, 506, 663, 837]. Другим примером вертикальной конвекции в газовом диске является неустойчивость в радиационно-доминирующей области [237, 303]. При учете ионизации водорода получаются неустойчивые вертикальные распределения температуры, приводящие к конвективному переносу тепла [514, 624, 814].

Характерной чертой аккрецирующих систем является широкий спектр нестационарных проявлений (карликовые новые различных типов: SS Лебедя, SU Большой Медведицы, Z Жирафа; рентгеновские пульсары; рентгеновские барстеры I и II типов; квазипериодические осцилляции рентгеновских источников; переменность излучения у кандидатов в черные дыры: Cyg X-1, GX 339-4, Cir X-1). Для объяснения некоторых переменностей излучения важным представляется изучение нестационарных течений, которые обусловлены развитием гидродинамических (МГД) неустойчивостей.

Наличие магнитного поля существенно усложняет гидродинамические течения в гравитационном поле компактного объекта, где разнообразие неустойчивостей плазмы, по-видимому, не меньше, чем для термоядерной плазмы [112]. Магнитные поля играют важную роль в динамике ряда аккреционно-струйных систем вблизи компактных звезд. В частности, влияют на структуру струй, характер и темп аккреции [17, 304, 653, 696].

Активное изучение турбулентности в аккрецирующих системах должно привести, в частности, к построению нефеноменологических моделей турбулентной вязкости [15, 149, 304, 354, 570]. Существенная доля газа в галактиках, протопланетных и протозвездных дисках и подавляющая его часть в аккрецирующих системах представляет собой полностью или частично ионизованную плазму. В плазме может существовать большое число неустойчивых мод, развитие которых эффективно турбулизует вещество. Без преувеличения можно сказать, что турбулентность является естественным состоянием для плазмы. Причем, одновременно могут работать много механизмов, приводящих к турбулентности через нелинейную стадию неустойчивостей.

В первом параграфе рассмотрены модели осесимметричных аккреционных дисков, в основе которых лежит предположение о переносе углового момента по диску за счет вязкости. В § 4.2 представлены в некотором смысле альтернативные результаты, связанные с возможностью образования спиральных волн (в том числе ударных) в газовых дисках. В § 4.3 рассмотрены некоторые магнитогазодинамические неустойчивости, которые при определенных условиях могут играть важную роль в аккреционных дисках. Интерес к такого рода исследованиям связан с тем, что неустойчивости могут приводить к эффективной турбулизации вещества в АД, определять механизм проникновения плазмы на замагниченные компактные объекты, приводить к различного рода нестационарным процессам.

4.1 Осесимметpичная дисковая аккpеция

4.1.1 Диффузионное приближение

Запишем уравнения нестационарной осесимметричной дисковой аккреции и попутно определим условия их применимости.

Сама идея механизма, приводящего к переносу углового момента в осесимметричном диске и, следовательно, аккреции вещества, достаточно проста. Она заключается в том, что в дифференциально вращающейся среде происходит взаимодействие между соседними слоями, связанное, например, с существованием магнитного поля, турбулентности, молекулярной или радиационной вязкости и т. п.^4.1. При типичных условиях молекулярная вязкость не может обеспечить величину темпа аккреции, вытекающую из наблюдений, а радиация сама является следствием аккреции. Магнитное поле попадает в диск, например, вместе с веществом, вытекающим из нормальной звезды. При определенных условиях величина магнитного поля может достигать $H^2/8\pi \sim \rho\,G\,M_1/r$, что является достаточным для аккреции.

В основе рассмотренной ниже осесимметричной модели лежит, как и в главе 4, предположение о том, что толщина диска $2h$ везде мала по сравнению с радиальной координатой $r$.

Запишем уравнения газодинамики с учетом диссипации для осесимметричного диска вокруг центрального тела массой $M_1$. Закон сохранения массы запишем с учетом источника вещества

$${\Vert\s \o \Vert t} + {\Vert \o r\Vert r}\Bigl(r\,\s\,u\B... ...= {1 \o 2\,\pi\,r} {\Vert\dot M_e \o \Vert r}. %\eqno(5.1.1)$$ (4.1)

Здесь $\s = 2\,h\,\overline\rho$ -- поверхностная плотность, $u$ -- радиальная скорость, $\dot M_e(r,t)$ -- заданная функция, описывающая темп притока вещества в АД. В законе сохранения момента импульса

$${\Vert \o \Vert t}\Bigl(\s\,l\Bigr) + {\Vert \o r\Vert r}\... ..._e\Vert\dot M_e \o 2\,\pi\,r\,\Vert r} - \s\,B %\eqno(5.1.2)$$ (4.2)

величина $l = rv$ есть удельный момент импульса вещества, находящегося на расстоянии $r$, $l_e(r,t)$ -- удельный момент импульса того вещества, которое определяется источником $M_e$, компонента тензора вязких напряжений, проинтегрированная по $z$-координате, имеет вид

$$W_{r\varphi} = \s\,\nu\left({\Vert v \o \Vert r} - {v \o r}\right),$$ (4.3)

где $\nu$ -- первая (сдвиговая) кинематическая вязкость. В (4.2) последнее слагаемое описывает эффективное приливное взаимодействие, обусловленное второй компонентой массы $M_2$ и $B$ -- удельный темп потери углового момента веществом.

Турбулентную среду можно при описании крупномасштабного движения рассматривать как жидкость, обладающую турбулентной вязкостью $\nu_t$, отличной от истинной (молекулярной) кинематической вязкости [92]. Развитую турбулентность можно рассматривать как иерархию турбулентных пульсаций, различающихся пространственными масштабами. Если через $v_t$ и $l_t$ обозначить соответственно скорость и масштаб основного (наиболее крупномасштабного) турбулентного движения, то по порядку величины можно записать

$$\nu_t \sim v_t\,l_t. %\eqno(5.1.4)$$ (4.4)

В предположении об изотропной турбулентности естественно предположить $l_t \,\lee\, h$. Поскольку сверхзвуковые турбулентные пульсации быстро диссипируют, то $v_t \,\lee\, c_s$. В результате имеем

$$\nu_t \,\lee\, c_s\,h. %\eqno(5.1.5)$$ (4.5)

В 1973 г. Шакура и Сюняев [790] отмечали: ``Сейчас, в отсутствие полной теории турбулентности, с одной стороны, и наблюдаемых данных существования турбулентности в дисках -- с другой, мы можем только предположить ее присутствие''. Это замечание не устарело по прошествии более трех десятков лет. Несмотря на это, модели аккреционных дисков с турбулентной вязкостью достаточно подробно разработаны и весьма популярны. Вопрос о вязкости в теории аккреционных дисков является, пожалуй, наиболее неясным. Как уже упоминалось, не вызывает сомнения, что молекулярная вязкость $\nu_{mol} \simeq c_s l_{ft}$ ($l_{ft}$ - длина свободного пробега) не может обеспечить необходимый темп аккреции. Поэтому, когда в теории дисковой аккреции идет речь о вязкости, обычно подразумевается турбулентная вязкость $\nu_t \gg \nu_{mol}$.

В отсутствие сколько-нибудь законченной теории турбулентности чрезвычайно важной стала работа Шакуры [221], который, по меткому выражению, ``свел все наше незнание турбулентности к одному безразмерному параметру $\alpha$'' [100]. Для тензора вязких напряжений принимается

$$W_{r\varphi} = -\,\alpha\,2\,h\,P, %\eqno(5.1.6)$$ (4.6)

где $P$ -- усредненное по $z$-координате давление; $\alpha$ -- свободный параметр. Модели, в основе которых лежит соотношение типа (4.6), принято называть $\alpha$-моделями или стандартными моделями АД.

В уравнении состояния вещества будем учитывать газовое, радиационное и магнитное давление

$$P = P_{gas} + P_{rad} + P_m, %\eqno(5.1.7)$$ (4.7)

где $P_{gas} = {\cal R}\,\s\,T/(2\,h\,\mu)$, $P_{rad} = a\,T^4/3$, $T$ -- температура, ${\cal R} = 8,\!31\,\cdot\,10^7$ эрг/(моль$\cdot$K), постоянная излучения эрг/(см$^3\cdot$K$^4$), $\mu$ -- молярная масса (в случае полностью ионизованной водородной плазмы $\mu = 0,\!5$). Для определения магнитного давления $P_m$ требуются дополнительные предположения о структуре магнитного поля.

Рассмотрим структуру диска в $z$-направлении. Будем полагать, что в $z$-направлении вещество находится в гидростатическом равновесии (см. ())

$${\Vert{\cal P}(z) \o \Vert z} = \rho(z){\Vert \o \Vert z}{G\,M_1 \o \sqrt{r^2 + z^2}} \,. %\qquad(a)$$ (4.8)

Уравнение баланса энергии в $z$-направлении имеет вид

$${1 \o \rho}{dF \o dz} = {F_\nu \o \s} \,, %\qquad(b)$$ (4.9)

где $F = - {\Oo{c \o 3\,k\,\rho}{d \o dz}(aT^4)}$ -- поток лучистой энергии; $F_\nu$ -- источник энергии, $k(\rho,T)$ -- непрозрачность; $c$ -- скорость света. Совместное решение уравнений (4.8, 4.9) определяет $z$-структуру тонкого диска. Учитывая выполнение $h \ll r$, условие гидростатического равновесия (4.8) можно записать в виде ^4.2

$$p = 2hP = C^2 \,\s\,h^2\,{G\,M_1 \o r^3}, %\eqno(5.1.9})$$ (4.10)

где $C^2\sim 1$ (см. п. ).

Обратимся к радиальной компоненте уравнения движения ()

$${\Vert(\s\,u) \o \Vert t} + {\Vert(r\,\s\,u^2) \o r\Vert r... ...{\Vert(r\,W_{rr}) \o r\Vert r} - {W_{\varphi\varphi} \o r} ,$$ (4.11)

где $W_{rr} = 2\,\s\,\nu\,(\Vert u/\Vert r) + \s\,(\mu - \nu)\,{\rm div}(\vec u)$; $\mu(\sim\nu)$ -- вторая кинематическая вязкость; $W_{\varphi\varphi} = 2\s\nu u/r + \s(\mu-\nu)\,{\rm div}(\vec u)$. Воспользуемся очевидной оценкой, вытекающей из (4.10):

$$c_s \simeq \sqrt{{p\o \s}} \simeq {h \o r}\,v \,, %\eqno(5.1.11)$$ (4.12)

тогда из (4.1), (4.2) с учетом (4.5) и (4.12) следует

$$\vert u \vert \sim \nu/r \sim c_s h/r \,. %\eqno(5.1.12)$$ (4.13)

Таким образом, с точностью до малой величины $(h/r)^2$ в радиальном направлении (см. (4.11)) имеется баланс только гравитационной и центробежной сил и можно записать

$$v = r\,\Omega = \sqrt{G\,M_1 \o r} \,. %\eqno(5.1.13)$$ (4.14)

Следует подчеркнуть приближенный характер соотношения (4.14), что особенно важнo при изучении динамики звуковых возмущений в плоскости диска. Диск не является строго кеплеровским, в противном случае в нем не могли бы распространяться крупномасштабные ($kh>1$) звуковые волны. Отличие скорости вращения от кеплеровской обусловлено практически только давлением, поскольку диссипативные члены в (4.11) малы $[\propto (h/r)^4]$.

Используя (4.14), исключим из (4.1) и (4.2) радиальную скорость, в результате получим эволюционное уравнение для поверхностной плотности

$${\Vert\s \o \Vert t}\! =\! {3\Vert \o r\Vert r}\biggl[\! {... ...o \Vert r} \!\biggr\} + {2\Vert \o r\Vert r}{B \s \o \Omega},$$ (4.15)

где $r_e = l^2_e / G\,M_1$ есть кеплеровский радиус, соответствующий удельному угловому моменту $l_e$, а последний определяется величиной $q = M_2/M_1$. Основным слагаемым, определяющим структуру диска, является в (4.15) первое слагаемое, обусловленное вязкостью. Остальные члены уравнения могут играть роль только во внешней области диска. Таким образом, имеем типичное уравнение диффузионного типа (что особенно наглядно видно из (4.37)). Для вязкости с учетом (4.3), (4.6), (4.10) можно записать

$$\nu = {2 \o 3}\,\alpha\,h^2\,\Omega. %\eqno(5.1.15)$$ (4.16)

Исходя из (4.1) и (4.2), нетрудно получить выражение для радиальной скорости

$$u = {1 \o \s}\left\{-{3 \o \sqrt{r} }{\Vert \o \Vert r}(\s... ...o \Vert r} - {2\,\s\,B \o r\,\Omega}\right\}. %\eqno(5.1.16)$$ (4.17)

В стационарном диске без учета двух последних слагаемых выражение для скорости принимает простой вид $u = -3\,\nu/(2\,r)$.

Рассмотрим уравнение баланса энергии для вещества диска в следующем виде:

$$\s\,T\,\left({\Vert S \o \Vert t} + u{\Vert S \o \Vert r}\... ...\nu - F_{rad} + F_{conv} + F_{\varkappa} + Y, %\eqno(5.1.17)$$ (4.18)

где $S$ -- удельная энтропия вещества. Справа стоят слагаемые, описывающие энергию, которая поступает на единицу площади (или уходит с нее) за единицу времени. Энергия выделяется вследствие вязкой диссипации:

$$F_\nu = {W^2_{r\varphi} \o \s\,\nu} = \nu\,\s\left(r{d\Omega \o dr}\right)^2 = {9 \o 4}\nu\,\s\,\Omega^2 \,.$$ (4.19)

Излучение с поверхности описывается членом $F_{rad}$. В рамках модели оптически толстого АД с учетом двух его сторон можно записать [222]:

$$F_{rad} = {4 \o 3}{a\,c\,T^4 \o \s\,k}.$$ (4.20)

Энергия может переноситься конвекцией в вертикальном направлении $F_{conv}$ (см., например, [814]). Величина $F_{\varkappa}$ обусловлена диффузионным переносом энергии вдоль радиальной координаты. Вещество, втекающее в диск, приносит не только момент импульса, но и энергию, что учитывает слагаемое $Y$. Эта энергия связана прежде всего с переходом в тепло кинетической энергии газовой струи, втекающей во внешнюю область АД с образованием ударной волны.

В качестве основных источников непрозрачности вещества в АД обычно принимают [222]: рассеяние на свободных электронах

$$k_{es} = 0,\!4\,f_e \ {\textrm{см}}^2/{\textrm{ г}},\qquad... ... = 0,\!4\,{\textrm{ см}}^2/{\textrm{ г}} \,; %\eqno(5.1.20a)$$ (4.21)

свободно-свободное поглощение (обратное тормозное излучение)

$$k_{ff} = 0,\!645\cdot 10^{23}{f_e f_i Z \o A}g_{ff}\,\rho\,T^{-7/2} {\textrm{ см}}^2/{\textrm{ г}}, %\eqno(5.1.20б)$$ (4.22)

где $g_{ff} \simeq 1$; $f_e$ -- число свободных электронов, приходящихся на барион; $f_i$ -- доля атомов, ионизованных до заданного состояния $Z$ -- заряд, $A$ -- атомный вес.

Таким образом, задача определения эволюции АД в рамках построенной модели сводится к совместному решению системы уравнений (4.15)$\div$(4.21 или 4.22), (4.7), (4.10). В результате для определенного начального состояния рассчитывается зависимость параметров диска ($\!\s$, $\!T$, $\!h$, $\!\nu$, $\!k$, $\!u$,...) от времени и радиальной координаты. При таком подходе функции $\dot M_e(r,t)$, $Y(r,t)$, $B(r,t)$ должны задаваться из каких-то дополнительных соображений. Процедура ввода в диск вещества и энергии описана в [280, 281, 282]. Учету приливного взаимодействия в рамках осесимметричных моделей посвящены работы [592, 691, 719, 814].

Характерными временами дисковой аккреции являются: динамическое время $t_{\varphi} \simeq 1/\Omega$, время установления гидростатического равновесия в $z$-направлении $t_z \simeq h/c_s$, тепловой временной масштаб $t_{th} \simeq \s\,c^2_s /F_{rad}$ и обусловленное вязкостью $t_{\nu} \simeq r/u \simeq (r/h)^2 t_{\varphi}/\alpha$. Соотношения между этими величинами зависят от радиальной координаты, но в целом можно считать

$$t_{\varphi} \simeq t_z \,\lee\, t_{th} \ll t_{\nu}.$$ (4.23)

В силу этих оценок, поскольку , часто пренебрегают теплоемкостью и ограничиваются следующим приближением (4.18): $\Oo F_{\nu} - F_{rad} + F_{conv} + F_{\varkappa} + Y = 0.$

Остановимся несколько подробнее на учете приливного взаимодействия. Под действием вязких сил диск в отсутствие приливного взаимодействия неограниченно ``расползается'' в радиальном направлении. Радиус диска стремится к бесконечности. Угловой момент благодаря трению отводится из внутренних областей диска во внешние и при наличии объекта-донора диск может простираться только до радиуса, на котором приливное взаимодействие способно отвести весь угловой момент из диска. В стационарном случае общий угловой момент вещества есть величина постоянная, следовательно, должен выполняться баланс

$$\dot M_e \sqrt{G\,M_1 r_e} = \int^R_{r_1}\s\,B\,2\pi r\,dr \,,$$ (4.24)

где $r_1$, $R$ -- внутренний и внешний радиусы диска соответственно. С учетом приливной силы АД оказывается конечным в радиальном направлении.

Большинство работ, в которых учитывалось в (4.15) приливное взаимодействие, основывается на результатах Папалоизу и Прингла [691]. Удобно величину $B$ представить в виде

$$B = {r^2\Omega^2 \o \pi\,q^2}\,f(r) \,, %\eqno(5.1.23)$$ (4.25)

где $q = M_1 / M_2$, $M_2$ -- масса объекта-донора, а $f(r)$ -- известная безразмерная функция [691]. Смак [814] использовал простую аппроксимацию

$$f = {\rm const}\cdot r^5 \,;$$ (4.26)

эта формула описывает результаты расчетов Папалоизу и Прингла [691] для случая $q \simeq 1$. Для отношения масс $q \gg 1$ можно записать

$$f \sim {1 \o [(a/r) - 1]^2 [1 - (\omega/\Omega)]^2} \,,$$ (4.27)

где $\omega^2 = G\,(M_1 + M_2)/a^3$, $a$ -- расстояние между объектами.

4.1.2 Стационарные модели

При изучении стационарных моделей достаточно во всех эволюционных уравнениях п. 4.1.1 положить . В рамках некоторых предположений о виде законов непрозрачности, вязкости и переноса энергии удается получить аналитические зависимости параметров АД от радиальной координаты. Обычно, следуя Шакуре и Сюняеву [790], диск разбивается на несколько областей, в которых преимущественную роль играют определенные процессы переноса излучения ($k_{es}$ или ) и давление ( или $P_{rad}$), приливное взаимодействие не учитывается, пpинимается $F_\nu = F_{rad}$.

1. В самой внутренней радиационно доминирующей ( $p_{rad} \gg p_{gas}$) зоне диска ``a'' с учетом выполнения $k_T \gg k_{ff}$ имеем:

$$h = 3\cdot 10^6 \dot m\,m\,{\cal L} \,\, \textrm{см},\qu... ...10^{10}\alpha\,\dot m^2 R^{-5/2}{\cal L} \,\, \textrm{см/с},$$ (4.28)

здесь $m = M_1/M_{\odot}$, $\dot m = \dot M /(m\cdot 3\cdot 10^{-8} {\rm M}_{\odot}/{\rm год})$, $R = r/(3R_g) = r/(9m\,{\rm км})$, ${\cal L} = R_1^{-1/2} - R^{-1/2}$. Появление функции ${\cal L}(r)$ обусловлено граничным условием $W_{r\varphi}(R = R_1) = 0$.

2. Для области ``b'', где $p_{gas} \gg p_{rad}$ и $k_T \gg k_{ff}$, можно записать:

$$\s = 0,\!81\cdot 10^5 \dot m^{3/5} m^{1/5}\alpha^{-4/5} R^{-3/5}\, {\cal L}^{3/5} \,\, \textrm{г/см}^2\,$$

$$T = 3,\!1\cdot 10^8 (\dot m\,{\cal L})^{2/5}(\alpha\,m)^{-1/5}R^{-9/10} \,\, \textrm{K},$$

$$h = 1,\!64\cdot 10^4 (\dot m\,{\cal L})^{1/5}m^{9/10}\alpha^{-1/10} R^{21/20} \,\, \textrm{ см},$$

$$u_r = 4,\!1\cdot 10^6\alpha^{4/5}\dot m^{2/5} m^{-1/5} R^{-2/5}{\cal L}^{-3/5} \,\, \textrm{см/с}.$$ (4.29)

3. Во внешней области диска ``c'' ( $p_{gas} \gg p_{rad}$, $k_T \ll k_{ff}$):

$$\s = 3,\!3\cdot 10^5 \alpha^{-4/5}\dot m^{7/10} m^{1/5} R^{-3/4} {\cal L}^{7/10} \,\, \textrm{г/см}^2,$$

$$T = 7,\!6\cdot 10^7\alpha^{-1/5}\dot m^{3/10} m^{-1/5}R^{-3/4}{\cal L}^{3/10} \,\, \textrm{K},$$

$$h = 8,\!1\cdot 10^3\alpha^{-1/10}\dot m^{3/20}m^{9/10}R^{9/8}{\cal L}^{3/20} \,\, \textrm{см},$$

$$u_r = 10^6\alpha^{4/5}\dot m^{3/10} m^{-1/5} R^{-1/4} {\cal L}^{-7/10} \,\, \textrm{см/с}. %\eqno(5.1.25{\it в})$$ (4.30)

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

В настоящее время имеется большое количество работ, посвященных определению радиальной структуры АД при тех или иных условиях. Описание некоторых моделей можно найти в работах [144, 209, 304].

4.1.3 Модели с вертикальной конвекцией

Наряду с лучистым переносом энергия из внутренних слоев диска может выноситься наружу конвективным движением (см. (4.18)). Однако радиационный перенос всегда имеет место при высокой температуре, а конвекция может и не возникать. Роль и условия возникновения конвекции в приложении к внутреннему строению звезд подробно исследованы (см., например, [21]). Имеется ряд факторов, не позволяющих формально перенести результаты теории звезд на аккреционные диски. Прежде всего отличается зависимость силы тяжести от координаты и имеется сильная вязкость. На важность учета конвекции в газовых дисках было указано Пачинским [687]; Вила [877, 878] рассматривал конвективные модели холодных дисков в катаклизмических двойных и массивных горячих дисках, а Лин и Папалоизу [582] исследовали протопланетный диск с конвекцией. Построению нестационарных моделей, основанных на детальном расчете $z$-структуры, посвящены работы [392, 616, 813, 814] и др. В приложении к галактическим ядрам конвективную неустойчивость рассматривали Минешига и Осаки [638].

Для конвективного потока тепла можно записать

$$F_{conv} = C_p\,\rho\,v\,l\,\left\{\left\vert{dT \o dz}\right\vert - \left\vert{dT \o dz}\right\vert_{ад}\right\},$$ (4.31)

здесь $v$, $l$ -- соответственно скорость восхождения и линейные размеры конвективных ячеек газа при перемешивании; -- теплоемкость. Как видно из (4.31), решающую роль играет превышение истинного градиента температуры над адиабатическим.

Мейер и Мейер-Хофмейстер [615] обнаружили неустойчивое распределение температуры, обусловленное ионизацией водорода. При ионизации водорода ( $T \simeq 10^4$ К) непрозрачность $k(\rho,T)$ сильно меняется по величине, так что радиационный механизм переноса энергии не может уже обеспечить необходимый темп. Образуется резкий перепад температуры между внутренними слоями, где в основном генерируется тепло, и внешними, из которых происходит высвечивание. В результате возникают условия для конвективного движения. При детальном расчете переноса излучения, естественно, нельзя пользоваться простыми соотношениями для непрозрачности типа (4.21), (4.22), тем более, что возможно нарушение приближения оптически толстого диска [616, 814]. Обычно используют таблицы непрозрачности (см., например, [21]).

Figure: Зависимость эффективной температуры поверхности АД $T_\textrm{эфф}$ от поверхностной плотности $\s$ при различных значениях $\alpha$ и $r$ (в единицах $10^{10}$ см) в случае центрального объекта массы $1$ M$_{\odot }$. Кривые получены: $S$ -- Смаком [814]; $MH$  -- Мейером и Мейером-Хофмейстером [615]; $MO$ --  Минешигой и Осаки [637]

=0.57

Типичная зависимость эффективной температуры поверхности диска $T_\textrm{эфф}(z = \pm h)$ от величины поверхностной плотности вещества показана на рис. 4.1. Каждая кривая имеет характерную $S\!$-форму -- состоит из трех участков: холодный $AB$, горячий $CD$ и переходный $BC$. Область $BC$ соответствует неустойчивости, тогда как горячие и холодные решения термически устойчивы. Поскольку вещество при $T \,\gee\, 10^4$ К практически полностью ионизовано, горячий участок определяется достаточно уверенно. Этого нельзя сказать об участке $AB$. В области низких температур закон непрозрачности известен хуже, холодный участок может быть оптически тонким, именно поэтому ветвь $AB$ менее определена и имеются отличия у разных авторов [876].

4.1.4 Модели карликовых новых

Характерной чертой карликовых новых (звезд типа U-Gem и дp.) является их нестационарное поведение (п. 1.5). В настоящее время предложено немало механизмов для объяснения феномена карликовых новых и в целом катаклизмических переменных [33]. Взрывной характер поведения многих систем обусловлен, по-видимому, нестационарным режимом дисковой аккреции. Весь вопрос заключается в определении местонахождения ``клапана'', который ``открывается'' на определенное время, что в конечном счете приводит к вспышке. Можно выделить три типа механизмов:

1. Причина нестационарности связана с нормальной звездой. Вследствие нестабильности вытекания вещества из красного карлика возникают квазипериодические колебания в накоплении вещества диском и тем самым в светимости газового диска [279]. Поскольку газ в АД поступает определенными порциями, то и нестационарная аккреция (и, как следствие, вспышка) есть просто отклик диска на меняющиеся внешние условия. В рамках изложенного в п. 4.1.1 подхода величина темпа поступления вещества в АД $\dot M_e(t)$ является функцией времени, которая должна быть задана.

2. Механизм квазипериодической активности может находиться в самом аккреционном диске [813, 814]. Если величина вязкости в диске достаточно мала, то во внешней области диска происходит накопление вещества. Когда плотность в диске достигает определенного критического значения, то в силу каких-то причин (развития неустойчивости, турбулизации среды) резко возрастает темп аккреции, что приводит к вспышке. Таким образом, клапан находится в самом АД.

3. Отсутствие аккреции между вспышками можно объяснить эффектом ``пропеллера'' [511]. Магнитное поле быстро вращающегося белого карлика препятствует падению вещества на его поверхность. Происходит накопление вещества вблизи границы магнитосферы, которое приводит к медленному приближению внутренней границы диска к белому карлику. Поскольку при этом из-за твердотельного характера вращения скорость движения силовых линий уменьшается, то в определенный момент эффект ``пропеллера'' исчезает, что приводит к мощной аккреции (вспышке). После этого система оказывается в исходном состоянии. Таким образом, в рамках описанного сценария клапаном является магнитное поле самого аккрецирующего объекта.

Неустойчивые АД. Обсудим результаты моделирования нестационарных АД, основанного на учете конвективной неустойчивости, рассмотренной в п. 4.1.3 [392, 616, 814, 876]. Весьма полное исследование было проведено Смаком [814]. Численно решались уравнения дисковой аккреции (п. 4.1.1) с учетом конвективного переноса в $z$-направлении (п. 4.1.3). В течение всего времени расчета темп поступления вещества в АД $\dot M_e$ и параметр $\alpha$ оставались постоянными. Основной интерес представляют временные зависимости светимости диска $L_d(t)$ и внешнего радиуса диска $R(t)$, поскольку эти величины являются наблюдаемыми (см. п. 1.5). Почти во всех случаях получены квазипериодические режимы аккреции. При этом наблюдаются два типа решений: тип А -- первоначально неустойчивость возникает во внешней зоне диска и, захватывая все более внутренние области, распространяется к центру АД (рис. 4.2$\!$а). Форма кривой светимости во время вспышки несимметрична, а повторяемость почти строго периодическая. Для типа B характерно возникновение неустойчивости вначале во внутренней области [ $r \simeq (0,\!1\div 0,\!2)\,R$] и последующее распространение внутрь и наружу АД ^4.3. При этом неустойчивость может в некоторых случаях не достигать внешней области диска (рис. 4.2$\!$б). В целом вспышки типа $B$ являются менее регулярными. При прочих равных условиях вспышки типа $A$ характерны для более высокого темпа притока вещества. Возможны комбинированные вспышки $AB$, когда одновременно возникают неустойчивости во внешней и внутренней областях АД.

Figure: =0.9 Временная эволюция светимости $L_d$ и внешнего радиуса $R$ АД [814]. Кружок указывает на местоположение зарождения неустойчивости, пунктирная линия показывает распространение неустойчивости по АД: а -- pешения типа A; б -- решения типа B

=0.34

Природа вспышечной активности обоих типов становится ясной при рассмотрении эволюции поверхностной плотности между активными фазами. В этот период величина поверхностной плотности меньше стационарного значения -- происходит накопление вещества. Когда поверхностная плотность достигает определенного значения, соответствующего критической величине -- точке поворота на кривой $T_\textrm{эфф}(\s)$ (см. рис. 4.1, точка $B$), развивается неустойчивость -- начинается вспышка. Высвечивающая энергия есть гравитационная энергия, то есть увеличивается радиальный поток вещества. Происходит распространение области неустойчивости в обе стороны от первоначального очага. Темп аккреции увеличивается, запасы вещества в диске уменьшаются и в конечном счете условия для конвекции исчезают, что приводит к прекращению вспышки.

Отметим отличительную особенность в поведении $R(t)$: в случае вспышки типа $A$ радиус диска увеличивается на $\,\lee\,$ 20 %, для вспышки типа $B$ характерно очень малое изменение величины $R$ ($\!\!\lee\,$ 7 %), а в некоторых случаях ``волна неустойчивости'' даже не доходит из внутренней области во внешнюю.

Конвективная неустойчивость поперек плоскости диска не является единственно возможной неустойчивостью, приводящей к квазипериодическому режиму аккреции. К аналогичным последствиям приводит рассмотренная авторами [145] градиентно-энтропийная неустойчивость в плоскости диска (см.  п.  и § 4.3). Для развития ГЭ-неустойчивости необходимо выполнение определенных соотношений между характерными масштабами неоднородностей поверхностной плотности $L_\s$ и температуры $L_T$. Пусть в начальный момент времени градиенты величин $\s(r)$ и $T(r)$ таковы, что АД является ГЭ-устойчивым. Однако стационарно поступающее на внешний край АД вещество ($\dot M_e =$ const) и нагрев за счет ``яркого пятна'' увеличивают градиенты поверхностной плотности и температуры, что в конечном счете создает необходимые для ГЭ-неустойчивости условия. Рост амплитуды возмущений в АД приводит к увеличению уровня турбулентной вязкости. Накопленное во внешней части АД вещество эффективно аккрецирует на компактный объект (вспышка)$\!$, и градиенты $\s(r)$ и $T(r)$ эволюционируют к значениям, при которых АД становится ГЭ-устойчивым. Не поддерживаемая неустойчивостью турбулентная вязкость и определяемые ею аккреционные процессы затухают, но стационарно поступающее на внешний край АД вещество подготавливает систему к новому аккреционному циклу (вспышке) (ход кривой $L_d(t)$ аналогичен изображенному на рис. 4.2).

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

Figure: Зависимость темпа притока вещества в АД $\dot M_e$ от времени и соответствующий отклик светимости $L_d$ [280]

=0.34

Модели с нестационарным притоком массы. Рассмотрим отклик аккреционного диска на увеличение темпа притока вещества $M_e(t)$ (обсуждение причин нестационарности величины $M_e$ выходит за рамки данной книги ^4.4). Басом и Принглом были проведены подробные расчеты нестационарного АД без учета приливного взаимодействия с функцией $M_e(t)$, типичный вид которой изображен на рис. 4.3$\!$а. На рис. 4.3$\!$б показано изменение светимости АД, представляющее собой отклик на внешнее воздействие [$\dot M_e(t)$] при постоянном значении параметра $\alpha$. Движение вещества в радиальном направлении обусловлено действием вязкости (из (4.17) следует оценка $u \sim -\nu /r$), то есть фактически величиной $\alpha$. Анализ динамики процессов и сравнение с наблюдениями позволяет оценить значение параметра $\alpha$. Удовлетворительное согласие достигается при $0,\!01 \,\lee\, \alpha \,\lee\, 0,\!5$.

Включение в расчет приливного взаимодействия в форме (4.26$\!$а) позволило Ливио и Вербунту [592] исследовать динамику радиуса диска $R(t)$, вызванную изменением темпа перетекания вещества. Пусть эволюция диска подчиняется уравнению (4.24), которое, как мы помним, является следствием закона сохранения момента движения. Перепишем (4.24) с учетом (4.25):

$$\dot M_e\sqrt{G\,M_1 r_e} = {2\,G\,M_1 \o q^2}\int^R_{r_1} \s\,f\,dr. %\eqno(5.1.27)$$ (4.32)

Левая часть соотношения (4.32) является линейной относительно $\dot M_e$. Под знаком интеграла только величина $\s$ зависит от $\dot M_e$; примем, что $\s \propto \dot M^n_e$. Из (4.29, 4.30) следует, что $n = 0,\!6;\, 0,\!7$. Эти значения являются достаточно типичными [679]. Если $n < 1$, то с ростом величины $\dot M_e$ будет увеличиваться радиус диска $R$. Основываясь на этом, можно качественно описать реакцию диска на увеличение $\dot M_e$. Пусть значению $\dot M_{e0}$ соответствует равновесный радиус $R_0$, а $\dot M_{e1}$ -- величина $R_1$. Если $\dot M_{e0} < \dot M_{e1}$, то $R_0 < R_1\!$. Втекающее в диск вещество обладает удельным угловым моментом $l_e$, которому соответствует радиус $r_e = l_e^2 / GM_1 < R_0$. Следовательно, при увеличении величины $\dot M_e$ от $\dot M_{e0}$ до $\dot M_{e1}$ в начальный момент диск сжимается до радиуса $R_2$ ( $r_e < R_2 < R_0$). В результате вместе с ростом плотности вещества во внешней области диска увеличивается и вязкость ( $\nu \sim \alpha\Omega h^2 \propto \s$). Тем самым внешняя граница диска отодвигается до нового равновесного радиуса $R_1$. В течение этого времени масса и светимость АД увеличиваются. Если затем темп притока вещества вернуть к прежнему значению $\dot M_{e0}\!$, то диск примет первоначальные размеры. Характерные времена определяются вязкими процессами. На рис. 4.4 показаны результаты численного решения уравнений нестационарной аккреции.

Figure: Эволюция АД вследствие изменения притока вещества $\dot M_e(t)$ [592]. Показаны функции: а -- $\dot M_e(t)$ и темп аккреции на центральный объект $\dot M(t)$; б -- $R(t)$

=0.35

Наиболее важным различием между рассмотренными двумя моделями (неустойчивый диск и переменный темп поступления вещества) является наличие во втором случае короткого промежутка времени, когда размер диска резко уменьшается и только потом возрастает. Имеются данные, свидетельствующие о такой особенности у Z Cha (см. п. 1.5.1).

Рассмотpенные выше модели являются пpедельно пpостыми, что связано в пеpвую очеpедь с феноменологическим подходом в постpоении вязких моделей АД, в основе котоpых лежат соотношения типа (4.16). Разумеется, пpи изменении состояния диска (плотности, темпеpатуpы и т. п.) вязкость может эволюциониpовать с существенной задеpжкой, наличие и наpастание мелкомасштабных магнитных полей также может игpать важную pоль [336].

4.1.5 Автомодельные нестационарные решения

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

$$F = - W_{r\varphi}r^2,\qquad \dot M = - 2\,\pi\,r\,\s\,u \,.$$ (4.33)

Величина $\dot M$ равна массе, проходящей через радиус $r$ за единицу времени. Если пренебречь влиянием второй компоненты в системе, то уравнение (4.1) при переходе к переменной $l = \sqrt{GM_1 r}$ примет вид

$${\Vert\s \o \Vert t} = {1 \o 2\,\pi\,r}{\Vert\dot M \o \Ve... ...(G\,M_1)^2 \o l^3}{\Vert\dot M \o \Vert l} \,. %\eqno(5.1.29)$$ (4.34)

Аналогично уравнение (4.2): $\Oo\dot M(r,t) = 2\,\pi\,{\Vert F \o \Vert l}$. Подставляя последнее в (4.34), получаем:

$${\Vert\s \o \Vert t} = {1 \o 2} {(G\,M_1)^2 \o l^3}{\Vert^2 F \o \Vert l^2} \,.$$ (4.35)

Для определения еще одной связи между $\s$ и $F$ можно воспользоваться уравнением баланса энергии (4.18) в приближении

$$F_\nu = F_{rad}. %\eqno(5.1.32)$$ (4.36)

В результате эволюционное уравнение примет вид

$${\Vert F \o \Vert t} = A {F^m \o l^n}{\Vert^2 F \o \Vert l^2},$$ (4.37)

где $A$ -- постоянная, значения показателей $m$ и $n$ зависят от выбора конкретной модели (законов вязкости и непрозрачности, уравнения состояния вещества). Так, например, если газовое давление преобладает над радиационным, то $m = 2/5$, $n = 6/5$ для случая $k_T \gg k_{ff}$ и $m = 0,\!3$, $n = 0,\!8$ в обратном пределе $k_{ff} \gg k_T$.

Пусть в начальный момент времени во внешней области диска на радиусе $r_0$ вещество находится в виде кольца, которое в последующем аккрецирует на компактный объект. В рамках автомодельного подхода получены три стадии аккреции. За время первой стадии вещество доходит до аккрецирующего объекта. На второй стадии вещество аккрецирует, темп аккреции $\dot M(t)$ и светимость $L(t)$ со временем растут:

$$\dot M \propto t^{2,47}\qquad (k_{ff} \gg k_T)\,; \quad \dot M \propto t^{1,66}\qquad (k_{ff} \ll k_T)\,. %\eqno(5.1.34)$$ (4.38)

Третья стадия характеризуется убыванием со временем величин $\dot M(t)$ и $L(t)$:

Во время стадии энергия в основном высвечивается из области $r \sim r_0$ и она мала в силу малости величины $r_1 /r_0$. В течение стадий II и III большая часть энергии уходит из области $r_1$ и оказывается пропорциональной темпу аккреции.

4.1.6 ``Толстые'' аккреционные диски

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

$$\s_T = {8\,\pi \o 3}\left({e^2 \o m_e c^2}\right)^2 = 0,\!665\cdot 10^{-24}\,\textrm{см}^2 \,. %\eqno(5.1.36)$$ (4.40)

В среднем фотон при столкновении передает весь свой импульс электрону, а затем и протону благодаря электростатической связи. Пусть светимость равна $L$ (эрг/с). Так как энергия фотона $\sim pc$, число фотонов, пересекающих в единицу времени единичную площадку, равно $N = L/(4\pi r^2 pc)$. За единицу времени электрон испытывает $N\s_T$ столкновений. Пусть в среднем за одно столкновение фотон передает электрону импульс $p$. Поскольку сила, действующая на электрон, есть скорость передачи импульса, то

$$f_{rad} = {L\,\s_T \o 4\,\pi\,r^2 c} \,. %\eqno(5.1.37)$$ (4.41)

Сила тяготения, действующая на протон, а следовательно, и на электрон, также пропорциональна $1/r^2$. Поэтому в приближении сферической симметрии существует критическая светимость (которую называют эддингтоновским пределом), определяемая балансом градиента радиационного давления и силы тяжести:

$$L_{crit} = {4\,\pi\,c\,G\,M_1 m_p \o \s_T} = 1,\!3\cdot 10^{38}(M_1/M_\odot) \quad \textrm{эрг/с} \,. %\eqno(5.1.38)$$ (4.42)

Критическому значению светимости соответствует критическое значение темпа аккреции $\dot M_{crit}$.

В случае сверхкритической дисковой аккреции при приближении к центральному телу диск перестает быть тонким ($r \sim 2h$). Радиус, на котором светимость достигает $L_{crit}$, называют радиусом сферизации:

$$r_\textrm{сф} = {G\,M_1\,\dot M_{crit} \o L_{crit}} = {\do... ... \o M_\odot/\textrm{год}}\, (\textrm{cм}) \,. %\eqno(5.1.39)$$ (4.43)

Аналогичная оценка получается из равенства $r_\textrm{сф} = 2h$ для величины $h$, определяемой (4.28). Из области $r < r_\textrm{сф}$ часть вещества начинает истекать в виде квазисферической оболочки. Как считают Шакура и Сюняев, в случае $\dot M > \dot M_{crit}$ реализуется самосогласованный режим аккреции, обеспечивающий светимость, лишь ненамного превышающую эддингтоновский предел. При этом в области $r > r_\textrm{сф}$ диск остается тонким и его структура аналогична докритическому режиму. Для черной дыры примем $r_\textrm{сф} = 10$ км, тогда из (4.43) получаем

$$\dot M_{crit} = {4\,\pi\,r_\textrm{сф} m_p\, c \o \s_T} = {3 \o 2} 10^{-8} M_\odot/\textrm{год} \,. %\eqno(5.1.40)$$ (4.44)

Обсудим структуру ``толстого диска'', близкого к гидростатическому равновесию, находящегося на балансе градиента давления, центробежной и гравитационной сил [227, 247, 285, 442, 516, 579, 688]. Прежде всего заметим, что при наличии существенного градиента давления скорость вращения вещества может весьма сильно отличаться от кеплеровской, и тем самым общий темп диссипирующей энергии может быть много меньше, чем значение $\sim G M_1\dot M /r_1$, характерное для тонких АД.

В теории аккреционных дисков вопрос о вязкости является наиболее невыясненным, и только благодаря предположению о малой толщине АД ($h/r \ll 1$) удается построить достаточно правдоподобные и простые модели. В случае же толстого диска вопрос о величине вязкости значительно усложняется, и в работах [227, 516, 688] построены модели толстых дисков без вязкости. Позже в ряде исследований [142, 143, 247, 285, 442, 579] рассматривалась область $r_1 \ll r \ll r_\textrm{сф}$ с учетом диссипативных процессов. Не проводя подробных выкладок, укажем только принципиальный подход в построении такого рода моделей в сферической системе координат $(r,\varphi,\Theta)$. Данный подход состоит в построении стационарных осесимметричных $(\Vert/\Vert\varphi \equiv 0)$ решений, характеризуемых взаимосогласованным распределением параметров течения газа в центральном гравитационном поле вдоль радиальной и меридиональной координаты $\Theta$^4.6. Это оказывается возможным, если соотношения между радиальными компонентами действующих на вращающееся вещество сил (гравитационной, центробежной и обусловленной градиентом давления) остаются постоянными вдоль радиальной координаты. В этом случае $\Theta$-распределения всех величин на различных расстояниях от центра будут подобны друг другу, а поток вещества чисто радиальным. Решение уравнений газодинамики ищут в форме $f(r,\Theta) = f_*(\Theta)\,r^{-a_f}$, где $f = \{\rho,\vec u,T,s,p,...\}$, $a_f=$ const. При этом параметры $a_f$ выбирают так, чтобы в каждом из уравнений все слагаемые были пропорциональны одной и той же степени радиальной координаты. Такие решения естественно называть автомодельными (ср. с п. 4.2.2). В результате проделанной процедуры уравнения в частных производных сводятся к обыкновенным дифференциальным уравнениям относительно переменной $\Theta$ (или $t \equiv z/r$, как у Лианга [579], в цилиндрической системе координат). Затем задаются определенными граничными условиями при $\Theta = \pi / 2$ и $\Theta = \Theta_0$ $(0 \le \Theta_0 < \pi/2)$; полученная таким образом краевая задача решается численно. В результате определяются распределения всех физических величин на плоскости $(r,\Theta)$.

Одним из главных достоинств автомодельных решений является простота их построения, что позволяет изучать влияние разнообразных факторов на структуру аккрецирующего течения. Однако, по-видимому, в рамках автомодельного приближения можно говорить только о получении качественной картины и трудно ожидать надежных количественных результатов. Кроме того, существуют проблемы, связанные с устойчивостью моделей толстых АД [671]. Заметим, что областью приложения моделей диссипативных толстых дисков может являться зона $r_1 \ll r \ll r_\textrm{сф}$ в ``тонких АД''.

Учитывая естественные трудности в построении аналитических (не автомодельных) двумерных моделей, весьма популярно численное моделирование осесимметричных течений в поле черной дыры, включая релятивистские эффекты. Уже первые расчеты показали, что в процессе падения холодного вещества с ненулевым моментом импульса происходит разогрев газа стоячей ударной волной, что приводит к формированию толстого ($\!h \sim r\!$) диска, поддерживаемого градиентом давления [900]. Фактически возникает тор, вращающийся с некеплеровской скоростью. В работе [467] моделировались области, наиболее близкие к черной дыре, в случае достаточно малого момента импульса аккрецирующего вещества. Газ довольно быстро затягивается по спиральной траектории в черную дыру, не испытывая, как считают авторы, вязкого взаимодействия.

ADAF- и CDAF-модели. В рамках стандартной $\alpha$-модели АД [790] имеется ряд проблем как теоретического плана, так и объяснения данных наблюдений. К последним относится объяснение низкой светимости рентгеновских двойных и активных галактических ядер с черными дырами ([662], см. обзор [664] и ссылки). Одним из наиболее изученных объектов такого типа является источник в центре Галактики Sgr $A^\ast$ [619]. Уменьшение светимости обеспечивается в так назывемых ADAF-моделях (адвективно доминирующие аккреционные течения) [231]. Они более горячие и эддингтоновская светимость достигается при меньшем темпе аккреции. Наблюдения за рядом объектов, например A0620-00, GRO J1655-40, центральный источник в активной галактике M 106, допускают интерпретацию спектральных кривых в рамках ADAF-моделей [457, 572, 661].

Начало исследованию ADAF-моделей положила работа [659], в которой были построены автомодельные решения стационарной аккреции, характерной особенностью которых является аккреция при слабом излучении или даже в отсутствие охлаждения за счет радиационных потерь. Основным механизмом переноса энергии выступает радиальный адвективный поток, а перенос углового момента обеспечивается $\alpha$-вязкостью. Структура течения оказывается близкой к сферической и вращение газа существенно отличается от кеплеровского закона. При наличии сильной вязкости ( $\alpha\gee 0,3$) численное моделирование дает результаты, сходные с ADAF-моделями [508]. Однако в расчетах с маленьким значением параметра $\alpha$ получаются течения с сильной турбулентностью из-за развития конвективной неустойчивости, которая существенно изменяет пространственную структуру течения [508, 837]. Такие модели получили название ``конвективно-доминирующие аккреционные течения'' (CDAF) [272].

Важнейшими свойствами данных моделей являются генерация конвективной турбулентности, перенос углового момента к центру, компенсирующий поток углового момента наружу вследствие вязкости (или, например, магнитовращающей неустойчивости [245, 270, 665 и ссылки там]), совместно с потоком тепловой энергии вдоль радиуса. Численное моделирование демонстрирует низкий темп аккреции и значительный поток энергии наружу за счет сильной радиальной конвекции. Фундаментальной проблемой моделей CDAF является наличие потока углового момента внутрь [272].

Для конвективно-доминирующего аккреционного потока удалось построить автомодельные решения в рамках приближения пути перемешивания [663], а в рамках численного моделирования изучена конвекция в газовом торе в плоскости $r$-$z$ [506]. Основные результаты по ADAF и CDAF получены в рамках осесимметричных моделей [506, 507, 508, 659, 660, 837]. Адвективное охлаждение начинает играть важную роль, когда диск перестает быть тонким $h\sim r$ [689]. Отношение скорости звука к кеплеровской скорости вращения в цитированных выше численных моделях достигает .

Переход к трехмерным ADAF-моделям, по-видимому, качественно не меняет выводы осесимметричных расчетов [509]. Отметим, что важный результат о конвективном переносе углового момента внутрь был обнаружен еще при построении моделей быстровращающихся звезд с конвективным ядром [19].

Модели конвективно-доминирующих аккреционных течений (CDAF) выглядят привлекательными для объяснения низкой светимости аккрецирующих черных дыр в рентгеновских двойных и галактических ядрах при наличии сверхмассивной черной дыры.

Наличие магнитного поля (MHD CDAF) может существенно изменять свойства течения, в частности, конвекция может приводить к потоку углового момента как внутрь, так и на периферию вдоль радиальной координаты, что, по-видимому, связано с влиянием магнитовращательной неустойчивости [272, 510]. В модели газового осесимметричного толстого диска при наличии магнитного поля показано, что возмущения с длиной волны, превышающей вертикальную шкалу диска, остаются конвективно неустойчивыми [665].

Отметим, что внутренние радиационно-доминирующие области АД могут быть неустойчивыми относительно вертикальной конвекции и в достаточно тонком диске [303]. Нелинейная стадия такой конвекции в тонком осесимметричном АД изучена для радиационно-доминирующей зоны стандартной модели АД для $r$-$z$-возмущений в работе [237].

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

4.2 Hеосесимметpичная дисковая аккpеция

4.2.1 Газодинамическое моделирование перетекания вещества в ТДС. Условия образования диска

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

Возникает ряд интересных вопросов: при каких условиях в ТДС возникает АД? Какая часть вещества теряется системой? Будет ли вещество аккрецировать без вязкости? Задачи такого рода являются для любителей аналитических решений практически неразрешимыми в силу нестационарности и неодномерности. И почти единственный выход -- численное моделирование.

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

Из самых общих соображений ясно, что при аккреции в ТДС возможно возникновение ударных волн. Бирман [297], по-видимому, был первым, кто в рамках гидродинамического подхода рассмотрел течение газа в близкой двойной системе в режиме звездного ветра. Методом характеристик было рассмотрено только сверхзвуковое течение. Заведомо такое решение не может содержать ударных волн. В работе [824] получена коническая ударная волна за аккрецирующим объектом. Однако используемый метод конечных разностей, имеющий первый порядок точности, приводит к слишком большой численной вязкости. Кроме того, декартова сетка не позволяет правильно задать граничные условия на поверхности обеих звезд.

В работах [609, 611, 766, 767, 768, 769, 830] применялись численные схемы второго порядка на криволинейной сетке, координатные линии которой близки к изолиниям эффективного потенциала системы, состоящей из двух тел ( $q = M_2 / M_1 = 1$), находящихся на расстоянии $a$ друг от друга и вращающихся с угловой скоростью $\Omega_0$. Одна из звезд заполняет свою критическую область Роша, а радиус другой не превышает $0,\!03 a$. Эффекты, связанные с охлаждением, нагревом, вязкостью ^4.7 и магнитными полями, не принимались во внимание. На поверхности нормальной звезды задавались значения плотности $\rho_0$ и скорости звука $c_0$.

Figure: Характерное расположение ударных волн для различных значений $c_0 /a\Omega _0$: a -- $c_0 \ll a\Omega _0$; б -- $c_0 \simeq a\Omega _0$; в -- $c_0 \,\gee\, a\Omega_0$

=0.35

Исследованию течений при различных $\beta = c_0 /a\Omega_0$ посвящены работы [609, 767] для $q = 1$. Если скорость звука мала ( $\beta \,\lee\, 0,\!75$), то вокруг компактного объекта возникает диск с двумя спиральными ударными волнами ^4.8 (рис. 4.5$\!$а). Максимальное число Маха $M=\vert \vec u \vert/c_s$ не превышает $3$. В случае $0,\!75\,\lee\,\beta\,\lee\,0,\!9$ происходит перестройка течения: диск становится менее выраженным, при этом остается только одна спиральная ударная волна (рис. 4.5$\!$б). При значениях $\beta$, лежащих в области $0,\!9\div 1,\!0$, возникает коническая ударная волна (рис. 4.5$\!$в), внутри конуса течение становится существенно дозвуковым. При дальнейшем увеличении скорости звука ( $\beta \,\gee\, 1$) образуется ярко выраженный режим звездного ветра. С ростом скорости звука на поверхности звезды-донора угол между ударными волнами становится меньше. Проходя через коническую ударную волну, скорость газа сильно уменьшается и часть его аккрецирует на компактный объект. Большая часть вещества из системы уходит. Похожие результаты получены в работе [795].

Таким образом, тип аккрецирующего течения (истечение с образованием диска или в форме звездного ветра с возникновением конической ударной волны) в системе с заполнившей свою полость Роша звездой-донором определяется значением параметра $c_0 /a\Omega _0$. Типичной для рассматриваемых систем является оценка $a\Omega_0 \sim 10^7$ см/с, что соответствует температуре $10^6$ К. В отсутствие звездной короны температура истекающего из звезды вещества много меньше $10^6$ K. Следовательно, наиболее вероятен режим истечения через внутреннюю точку Лагранжа с образованием аккреционного диска вокруг компактного объекта. При перетекании вещества через внутреннюю точку Лагранжа велик удельный угловой момент вещества, что приводит к образованию диска. В случае звездного ветра удельный угловой момент достаточно мал и диск не образуется [511].

В работах [16, 302, 648] изучалась структура перетекания газа в ТДС в рамках трехмерных моделей. Подробно исследованы эффекты взаимодействия газового потока через внутреннюю точку Лагранжа $L_1$ (см. рис. 4.5). Результаты расчетов демонстрируют более сложный характер течения по сравнению с двумерными. Формируются потоки газа с небольшим удельным угловым моментом выше плоскости системы, которые приводят к образованию горячей короны. Трехмерное численное моделирование позволило рассмотреть изгибные неустойчивости АД [647].

4.2.2 Автомодельные ударные волны

Предположение о том, что в газовых дисках, вращающихся вокруг компактных объектов, могут возникать спиральные ударные волны, высказывалось неоднократно [597, 805]. Притягательность их изучения связана с тем, что спиральные ударные волны могут переносить угловой момент из внутренних областей диска во внешние. В тонких аккреционных дисках ($h/r \,\ll\, 1$) течение является сверхзвуковым, что допускает возможность возникновения ударных волн. Причинами возникновения ударных волн могут являться вторая компонента в системе либо асимметричная магнитосфера вокруг компактного объекта [625]. Благодаря диссипативным процессам на фронте волны, вещество может по спирали падать на центр.

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

$${\Vert \o r\,\Vert r}(h\,\rho\,u\,r) + {\Vert \o r\,\Vert\varphi}(h\,\rho\,v) = 0\,, %(5.2.1)$$ (4.45)

$$\rho\left(u{\Vert u \o \Vert r} + {v\Vert u \o r\Vert\varp... ...ht) = -{\Vert p \o \Vert r} - \rho{G\,M_1 \o r^2}\,, %(5.2.2)$$ (4.46)

$$\rho\left(u{\Vert v \o \Vert r} + {v\Vert v \o r\Vert\varp... ...{u\,v \o r}\right) = -{\Vert p \o r\Vert\varphi}\,. %(5.2.3)$$ (4.47)

Будем полагать, что полутолщина $h$ может зависеть только от радиальной координаты. В частности, изучим два случая:

$$h = h_0 = {\rm const}, %\eqno(5.2.4)$$ (4.48)

Первый случай соответствует строго двумерному течению в плоскости диска. Во-втором принимается, что газ находится в гидростатическом равновесии в $z$-направлении [ср. с (4.12)]. В (4.49) величина $\hat{T}(r)$ -- хаpактеpная температура на расстоянии $r$, $\Omega(r) = \sqrt{ G\,М_1/r^3 }$. Уравнение для энтропии с учетом излучения с поверхности диска запишем в форме

$$h\,\rho\,T\,\left(u{\Vert S \o \Vert r} + {v\,\Vert S \o r\,\Vert\varphi}\right) = -\s_{SB} T^4_s \,, %\eqno(5.2.6)$$ (4.50)

где $T_s$ -- температура поверхности диска, $\s_{SB}$ -- постоянная Стефана-Больцмана. Для энтропии в случае идеального газа имеем:

$$S = c_V \ln\{p/\rho^\gamma\} \,. %\eqno(5.2.7)$$ (4.51)

В случае термодинамического равновесия для большой оптической толщины $\tau \gg 1$ можно принять [186]:

$$T^4 = {3 \o 4}\tau\,T^4_s,\qquad \tau = h\,\hat{\rho}\,k \,.$$ (4.52)

Считаем, что непрозрачность $k$ и хаpактеpная плотность $\hat{\rho}$ являются функциями только радиальной координаты. Азимутальную скорость представим в виде

$$v = r\,\Omega + v_1 \,. %\eqno(5.2.9)$$ (4.53)

Пользуемся безразмерными величинами, которые пометим сверху значком `` $\>\>\tilde {}\>\>$'':

$$(u,v_1) = \Omega_0 r_0(\tilde u,\tilde v);\quad h = r_0\t... ...ho_0\tilde p; \quad\rho = \rho_0\tilde\rho\,, %\eqno(5.2.10)$$ (4.54)

где $r_0$ -- некоторый характерный радиус; $\rho_0$ -- произвольный масштаб плотности. Введем новые ``спиральные'' координаты [829]:

$$x = r/r_0\,, \quad %(5.2.11) \Psi = \varphi + \beta(x)\,, %(5.2.12)$$ (4.55)

функция $\beta(x)$ будет определена ниже, но ясно, что случай $\Psi =$ const определяет спираль. Считаем, что непрозрачность является функцией только радиальной координаты

$$k = k_0 x^\nu. %\eqno(5.2.13)$$ (4.56)

С учетом (4.53)$\div$(4.55) уравнения (4.45)$\div$(4.47) примут вид:

$$\tilde\rho\biggl[\biggl(x^{-3/2} + {\tilde v \o x}\biggr)... ...x} - {d\beta \o dx}{\Vert\tilde p \o \Vert\Psi}\,, %(5.2.14)$$ (4.57)

$$\tilde\rho\biggl[\biggl(x^{-3/2} + {\tilde v \o x}\biggr)... ...\o x}\biggr] = -{\Vert\tilde p \o x\,\Vert\Psi}\,, %(5.2.15)$$ (4.58)

$${\Vert \o \Vert x}(\tilde h\,\tilde \rho\,\tilde u\,x) + {... ...\o dx}\tilde u\,x + \tilde v + x^{-1/2}\biggr)\biggr] = 0\,.$$ (4.59)

Решения ищем в автомодельном виде $q_i = x^n Q_i(\Psi)$, $n = n_i$. Величины $n_i$ определяются из условия независимости уравнений от переменной $x$. Для этого необходимо положить:

$$\tilde u = x^{1/2}U(\Psi),\quad \tilde v = x^{-1/2}V(\Psi)\,,$$

$$\tilde\rho = x^{-n}R(\Psi), \tilde p = x^{-n-1}P(\Psi),\... ...ta \o dx} = x^{-1}B,\quad \tilde h = x^m C\,, %\eqno(5.2.17)$$ (4.60)

Figure: К вопросу об автомодельной аккреции

=0.6

где $B$, $C$, $m$, $n$ -- постоянные. Если $\Theta$ есть угол между касательной к спирали ($\Psi =$ const) и радиальным направлением (рис. 4.6), то

$$B = {\rm tg}(\Theta). %\eqno(5.2.18)$$ (4.61)

Нетрудно заметить, что случай $B =$ const соответствует логарифмической спирали. Показатель $m$ равен нулю для диска постоянной толщины (4.48). В случае закона (4.49) $m=1$ и для температуры справедливо . Если выбрать в качестве характерного значения температуру при $\Psi = 0$, то

$$C = \sqrt{ P(0)/R(0) } \,. %\eqno(5.2.19)$$ (4.62)

Как видим, параметр $C$ есть изотермическая скорость звука в единицах кеплеровской скорости $\Omega_0 r_0$ или обратное число Маха. С учетом (4.60) уравнения (4.57)$\div$(4.59) принимают вид обыкновенных дифференциальных уравнений относительно $\Psi$:

$$W\,U' - {1 \o 2}U^2 - V^2 - 2V = (n+1){P \o R} - B{P' \o R},$$ (4.63)

$$(R\,W)' + \left({1 \o 2} + m - n\right)\,R\,U = 0,% \eqno(5.2.22)$$ (4.65)

здесь $W = 1 + V + B\,U$. Определим величину $n$. Проинтегрируем (4.65), в результате получим

$$(R\,W)\Big\vert^{2\pi}_0 + \left({1 \o 2} + m - n\right)\, \int^{2\pi}_0 R\,U\,d\Psi = 0. %\eqno(5.2.23)$$ (4.66)

Первое слагаемое в (4.66) равно нулю. Второй член определяет темп аккреции, который отличен от нуля в случае наличия ударных волн. Следовательно,

Итак, в случае $h =$ const $n = 1/2$, а для $h \propto x$ $n = 3/2$. С учетом (4.67) уравнение непрерывности (4.65) можно переписать в виде

$$(R\,W)' = 0. %\eqno(5.2.25)$$ (4.68)

Для сохранения автомодельности уравнения (4.50), с учетом (4.51), (4.52), (4.56), необходимо положить $\nu = - 1 / 2$, тогда закон изменения энергии принимает вид:

$$W\left({P' \o P} - \gamma{R' \o R}\right) + U\,(\gamma\,n - n - 1) + F{P^3 \o R^4} = 0, %\eqno(5.2.26)$$ (4.69)

где

$$F = {4 \o 3}{\s_{SB} T^3_0 \o \rho^2_0 k_0 c_V\Omega_0 r^2_0 C^2}.$$ (4.70)

Здесь $T_0{\cal R} /\mu = p_0 /\rho_0 = \Omega^2_0 r^2_0$.

Систему уравнений (4.63 $\div$ 4.70), (4.69) относительно неизвестных $U$, $V$, $P$, $R$ необходимо дополнить граничными условиями. Рассмотрим $N$ одинаковых ударных волн, разделенных фиксированным углом $\Delta \Psi$, тогда решения должны быть периодичны с периодом $\Delta \Psi = 2\pi/N$. Запишем выражения для нормальной и касательной к линии $\Psi = const$ компонент скорости (рис. 4.6):

$$V_\Vert = \sin\Theta\left\{{U \o B} - 1 - V\right\} = (1 + B^2)^{-1/2} \Bigl\{U - B\,(1 + V)\Bigr\}, %\eqno(5.2.28)$$ (4.71)

$$V_\perp = \cos\Theta\Bigl\{1 + V + B\,U\Bigr\} = (1 + B^2)^{-1/2}W. %\eqno(5.2.29)$$ (4.72)

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

$$\Bigl\{R\,W\Bigr\} = 0\,, %(5.2.30)$$ (4.73)

$$\Bigl\{V_\Vert\Bigr\} = 0 \,, %(5.2.31)$$ (4.74)

$$\Bigl\{P + R\,W^2/(1 + B^2)\Bigr\} = 0\,, %(5.2.32)$$ (4.75)

$$\Bigl\{{\gamma \o \gamma - 1}\,{P \o R} + {1 \o 2}\,{W^2 \o 1 + B^2}\Bigr\} =0\,. %(5.2.33)$$ (4.76)

Здесь через $\{A\}$ обозначена разность значений величины $A$ по разные стороны фронта ударной волны. Условие (4.73) в силу (4.68) удовлетворяется автоматически. Условия (4.74), (4.75) и (4.76) определяют только три постоянные интегрирования. Для получения четвертого условия сведем систему (4.63), (4.64), (4.68), (4.69) к уравнению

$$\left[(1 + B^2)\,c^2 - W^2\right]\,V' = {1 \o 2}UW(V + 1) ... ... - {n + 1\o \gamma}\right) + F\,{P^4\o R^5}\,, %\eqno(5.2.34)$$ (4.77)

где $Q = \frac{1}{W}{\Oo\biggl[{1 \o 2}\,U^2 + V^2 + 2V +{1 \o 2}\,B\,U\, (V + 1) + (n + 1)\,\frac{c^2}{\gamma}\biggr]}$,  $c^2= \gamma P/R$ -- безразмерная адиабатическая скорость звука. Уравнение (4.77) при выполнении $(1 + B^2)\,c^2 = W^2$ (или с учетом (4.72) $V^2_\perp = c^2$]) имеет сингулярность, которая соответствует наличию звуковой точки при $\Psi = \Psi_s$. Равенство нулю правой части (4.77) в звуковой точке $\Psi_s$ дает четвертое условие

$$\left[{1 \o 2}\,U\,W\,(V + 1) - B\,Q\,c^2 - U\,c^2\left(n ... ...\o R^5}\right]_{\Big\vert \Psi=\Psi_s} = 0\,. %\eqno(5.2.35)$$ (4.78)

Условия (4.74 $\div$ 4.76), (4.78) определяют постоянные интегрирования уравнений (4.63), (4.64), (4.68), (4.69).

В предельном случае большого числа ударных волн можно решить задачу аналитически [829]. Рассмотрим только адиабатическое течение [$F \equiv 0$ в (4.69)]. Если исходить из малости параметра $\epsilon = 1/N \ll 1$ и предположения о том, что функции $U,~V,~R,~P$ являются линейными между ударными волнами, то можно записать соотношение между $\gamma$ и углом спирали $\Theta$ [829]:

Figure: Зависимость угла $\Theta$ от $\gamma$ для разного числа ударных волн [символ ``$\infty$'' соответствует аналитическому решению (4.40)]. Тонкие линии относятся к модели диска с постоянной толщиной, жирные линии -- к модели диска с (4.49)

=0.6

Численный подход к решению сформулированной выше задачи позволяет рассматривать произвольное число ударных волн, в том числе с учетом радиационных потерь. Результаты такого рода расчетов приведены на рис. 4.7. Включение радиационных потерь позволяет оценить эффективный $\alpha$-параметр, фигурирующий в ``вязких осесимметричных моделях'' (см. § 4.1). Если определить средний радиальный поток

$$\langle u\rangle = \int u\rho\,d\varphi \bigg/ \int\rho\,d\varphi$$ (4.80)

и темп аккреции

$$\dot M = - 2\,\pi\,r\,\s\,\langle u\rangle \,, %\eqno(5.2.38)$$ (4.81)

то, сравнивая с результатом, вытекающим из стандартной модели АД

$$\dot M = 3\,\pi\,\s\,\nu_\textrm{эфф} = 3\,\pi\,\s\,\alpha_\textrm{эфф}\Omega\, h^2 \,, %\eqno(5.2.39)$$ (4.82)

получим

Figure: Зависимость $\alpha_\textrm{эфф}$ от угла $\Theta$ для двух ударных волн

=0.6

На рис. 4.8 показан коэффициент $\alpha_\textrm{эфф}$ как функция угла $\Theta$. В случае двух ударных волн ($N=2$) имеется максимум при $\Theta \simeq \pi/4$ и $\max\{\alpha_\textrm{эфф}\} = 1,\!02\cdot 10^{-2}$.

4.2.3 Спиральные ударные волны в ТДС.
Газодинамическое моделирование

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

Прежде всего, в работах [766, 767] было показано, что:

1) газ теряется нормальной звездой через окрестность точки $L_1$ в форме сверхзвуковой струи (см. рис. 1.3);

2) основная часть вещества вращается вокруг компактного объекта в форме аккреционного кольца/диска;

3) в результате приливного взаимодействия образуется две или три спиралевидные ударные волны (УВ);

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

5) система может терять значительную часть вещества через точку $L_2$;

Figure: =10000 Характерное поведение величин $\dot M_e(t)$ и $\dot M(t)$ по результатам экспериментов

=0.32

6) величина темпа потери вещества оптической звездой $\dot M_e(t)$ может достаточно сильно осциллировать, в то время как темп аккреции $\dot M(t)$ является более гладкой функцией (рис. 4.9);

7) отношение $\dot M /\dot M_e$ сильно зависит от параметров системы и составляет $\sim 20\div 90$ %.

Причиной возникновения ударных волн является вторая компонента, то есть генератор находится на периферии АД, тем самым возникает вопрос о том, как близко к аккрецирующему объекту могут простираться УВ. Для решения этой проблемы была проведена серия экспериментов [768], в которых размер компактного объекта равнялся $r_* = 0,\!01\,a$. Поскольку для тесных двойных с периодом от нескольких часов до дней величина $a$ составляет $\sim 10^{10} \div 10^{11}$ см, то $r_* \simeq 10^8 \div 10^9$ см, что соответствует радиусу белого карлика. Если компактным объектом является нейтронная звезда с магнитным полем $\sim 10^{12}$ Гс, то диск разрушается на расстоянии $\sim 5 \cdot 10^8$ см [401]. Расчеты убедительно продемонстрировали, что ударные волны простираются вплоть до $r = r_*$.

Обсудим влияние численной вязкости. Используемые численные схемы для решения уравнений газодинамики имеют II порядок точности и дают схемное число Рейнольдса $R_c \sim (r/\delta r)^2$ ($\delta r$ -- размер ячейки). Вблизи компактного объекта $R_c \sim 100$. Таким образом, угловой момент отводится наружу и газ падает на центр даже в случае осесимметричного потенциала (без ударных волн). Эффект численной вязкости можно снизить, уменьшая величину $\delta r$. Для ответа на вопрос: какая часть углового момента теряется в УВ, был поставлен эксперимент [766], в котором в момент времени $t = 22/\Omega_0$ (диск находится в состоянии квазистационара) каждая пространственная ячейка в радиальном направлении делилась пополам и расчет продолжался до $t = 27/\Omega_0$. В целом глобальная структура течения не изменялась, а усредненная величина $\dot M /\dot M_e$ уменьшалась от $0,\!9$ до $0,\!7$. Таким образом, по оценкам авторов, около 60$\div$ 70 % общих потерь углового момента связано с ударными волнами.

Figure: =10000 Зависимость $\tau$ от $q$ по результатам экспериментов [611]

=0.6

Процесс аккреции удобно характеризовать временем аккреции $\tau = {\cal M}/\dot M$ (${\cal M}$ -- масса диска). На рис. 4.10 показана экспериментальная зависимость величины $\tau$ от отношения масс компонент $q = M_2/M_1$ [611]. Горизонтальная линия соответствует осесимметричной модели ($q = 0$), в которой аккреция полностью обусловлена численной вязкостью. В рамках вязкой стандартной модели АД величина $\tau = r^2 /\nu$ (п. 4.1.1). Для вязкости $\nu = \alpha\,c^2_s/\Omega$ имеем

$$\tau = {1 \o \alpha\,\Omega}\left({\Omega\,r \o c_s}\right)^2 = {M^2 \o \alpha\,\Omega},% \eqno(5.2.41)$$ (4.84)

где $M$ -- число Маха. Принимая для внешних областей $M \sim 3$ и обращаясь к рис. 4.10, для $q = 10^{-3} \div 1$ получим $\alpha_\textrm{эфф} \simeq 5 \cdot 10^{-4} \div 4 \cdot 10^{-1}$.

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

Обсудим результаты, вытекающие из описанного выше газодинамического моделирования, в сравнении с автомодельными решениями (п. 4.2.2). Из рис. 4.7 видно, что стационарные автомодельные решения, содержащие две спиральные УВ в диске постоянной толщины, невозможны для $\gamma > 1,\!6$. Численное моделирование при $\gamma =5/3$ приводит к сильно осциллирующим течениям (см. рис. 4.9), то есть стационарные решения также не получаются. При $\gamma = 1,\!2$ осцилляции малы (рис. 4.9), и непосредственное сопоставление угла спирали $\Theta$ автомодельной волны с экспериментальными результатами дает удовлетворительное согласие. Сравнению результатов численного моделирования ударных волн, автомодельных решений и стандартной теории дисковой аккреции посвящена работа [611]. Зависимость угла закрутки УВ $\Theta$ от показателя адиабаты ^4.9 показана на рис. 4.11. Измерения $\Theta$ относятся к внутренней зоне АД, где влияние второй компоненты минимально. В области $\gamma \,\lee\, 1,\!4$ имеется хорошее согласие. В численных экспериментах при $\gamma \,\gee\, 1,\!55$ две стационарные ударные волны не появлялись, автомодельный подход также запрещает их существование при $\gamma > 1,\!595$ (см. рис. 4.7, 4.11). В области $1,\!595 < \gamma < \sqrt 3$ возможны стационарные решения с числом УВ больше двух. На рис. 4.12 показаны радиальные зависимости числа Маха ударной волны $M$. Наблюдается существенное различие по сравнению с автомодельными решениями во внешней области АД, которое уменьшается при приближении к центру. Такое поведение, по-видимому, вызвано тем, что приливное взаимодействие при построении автомодельных решений не учитывалось.

Figure: Зависимости угла $\Theta$ от $\gamma$. Сплошная линия -- автомодельное решение для двух УВ ($N=2$), пунктирная линия -- для $N =~\infty$ (см.  рис.5.7). Кружком показаны результаты по данным численных экспериментов, приходящих к стационарному состоянию, а звездочкой -- для моделей, далеких от стационарного состояния [611]

=0.6

Figure: Зависимость числа Маха ударной волны от радиальной координаты. Линиями показаны автомодельные решения

=0.6

Удивительным, на первый взгляд, аспектом вышеописанных результатов является возможность аккреции без радиационных потерь. В противоположность этому в рамках стандартной дисковой аккреции вся диссипирующая энергия высвечивается. В связи с этим рассмотрим автомодельные волны с радиацией (п. 4.2.2). Зависимость температуры от темпа аккреции в случае гидростатического равновесия в $z$-направлении показана на рис. 4.13. Здесь безразмерный темп аккреции $\dot m$ определен следующим образом:

Figure: Зависимость безразмерной температуры $y = T/T_\textrm{вир}$ $[T_\textrm{вир} = \mu\,G\,M_1 /({\cal R}\,r)]$ от темпа аккреции для дисков с автомодельными ударными волнами для различных показателей адиабаты. Сплошная линия -- с учетом радиационных потерь, пунктирная -- адиабатическая аккреция

=0.6

$$\dot m = \dot M/K,\quad K = 4\sqrt{\pi\,a\,c\,(\gamma - 1... ...}}\right)^2 (G\,M_1)^{7/4} k^{-1/2} r^{-1/4}, %\eqno(5.2.42)$$ (4.85)

где $a$ -- постоянная излучения; $\mu$ -- средняя молярная масса; $k$  -- непрозрачность. Существует критическое значение величины $\gamma = \gamma_{crit} \simeq \sqrt 6 - 1 \simeq 1,\!45$. При $\gamma < \gamma_{crit}$ температура достигает конечной величины, в то время как темп аккреции не ограничен. В случае высокого темпа аккреции течение становится настолько оптически непрозрачным, что радиационные потери не играют роли, и решения асимптотически стремятся к адиабатическим решениям (горизонтальная пунктирная линия). Следуя [611], покажем, что $\alpha$-моделям присуще аналогичное поведение.

Уравнение (4.18), выражающее баланс энергии, запишем для стационарного случая в виде

$$\s\,T\,u\,{dS \o dr} = {9 \o 4}\,\s\,\nu\,\Omega^2 - {4\,a\,c\,T^4 \o 3\,\s\,k}. %\eqno(5.2.43)$$ (4.86)

С помощью соотношений $\dot M = -2\pi r\s u$, $u = - 3\nu/(2r)$ исключим $\s$ из (4.86)

$$\dot M^2\left({2 \o 3}\,T\,r\,{dS \o dr} + \Omega^2 r^2\right) = {16\,\pi^2 a\,c\,T^4 \o 3\,k}\,\nu\,r^2. %\eqno(5.2.44)$$ (4.87)

Если предположить $k = k_0 r^{-1/2}$, то возможны автомодельные решения с $h \propto r$, $p \propto r^{-5/2}$, $\rho \propto r^{-3/2}$. Для идеального газа, принимая во внимание (4.16), запишем уравнение для величины $\dot M$:

$$\dot M^2\left[y\,(\gamma - 5/3) + \gamma - 1\right] = {32\... ... {\cal R}\right)^4 (G\,M_1)^{7/2}\alpha\,y^5,% \eqno(5.2.45)$$ (4.88)

где $y = {\cal R}\,T\,r/(\mu\,G\,M_1)$ -- отношение температуры диска к вириальной температуре. С учетом (4.85) соотношение (4.88) можно переписать в виде

$$\dot m = {\sqrt{ 2\,\pi\,\alpha} \o 3}{y^{5/2} \o \sqrt{(\gamma - 5/3)\,y + \gamma - 1}}. %\eqno(5.2.46)$$ (4.89)

Зависимости $y(\dot m)$ аналогичны случаю с ударными волнами -- ниже критической величины $\gamma$ температура диска стремится к своему асимптотическому значению. В предельном случае адиабатической аккреции ( $a \rightarrow 0$) из (4.88) имеем

$$y = {\gamma - 1 \o 5/3 - \gamma}. %\eqno(5.2.47)$$ (4.90)

Из (4.90) следует, что адиабатическая аккреция возможна только для $\gamma < 5/3$. Для значений $\gamma$, близких к единице, температура диска мала $y \ll 1$. Как видим, существует максимальное значение показателя адиабаты, ниже которого аккреция может идти с произвольным темпом. В случае $\gamma > \gamma_{crit}$ имеется максимально возможный темп аккреции. Значение $\gamma_{crit}$ определяется выбором модели АД. Заметим, что и для сферической аккреции существует критическое значение величины $\gamma$ и $\gamma_{crit} = 5/3$ [313]. Последнее связано с тем, что гравитационное поле в обоих случаях одинаковое.

Совершенно очевидно, что обсуждаемые здесь ударные волны весьма сходны с изучаемыми в теории спиральной структуры галактик. Рассмотрим (возможно слабый) источник неосесимметричных возмущений во внешней области диска. Им может являться не только вторая компонента, но и, например, какая-либо неустойчивость. Такое возмущение распространяется по диску, принимая спиральную форму благодаря дифференциальности вращения. Волны в такой ситуации переносят угловой момент, взаимодействуя с веществом диска. Эта проблема широко обсуждалась в приложении к галактикам [658, 802, 853].

4.3 Hеустойчивости в аккpеционных дисках

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

Вещество в АД представляет собой полностью или частично ионизованную плазму. В плазме может существовать большое число неустойчивых мод [112, 113], развитие которых эффективно турбулизует вещество. Можно сказать, что турбулентное состояние является естественным для плазмы. Как показывают оценки, такие моды в основном мелкомасштабны по сравнению с толщиной диска. Но для последовательного решения проблемы устойчивости наши представления о физических условиях на этих масштабах недостаточны. По этой же причине автоматический перенос достижений физики плазмы на изучаемые здесь системы не дает реалистичной картины.

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

4.3.1 Неустойчивость радиационно-доминирующей области

Рассмотрим стандартную модель аккреционного диска ($\alpha$-модель, см. § 4.1). Температура газа в АД растет с приближением к аккрецирующему объекту. Действительно, исходя из баланса энергии, выделяющейся вследствие действия вязкости (4.19) и уносимой излучением (4.20) в приближении чернотельного излучения, имеем

$${1 \o 2}\, W_{r\varphi} r {d\Omega \o dr} = {3 \o 8\pi}\do... ... \sqrt{ {r_1 \o r} }\,\,\right] = \s_{SB} T^4, %\eqno(5.3.1)$$ (4.91)

где $\s_{SB}$ -- постоянная Стефана-Больцмана; $r_1$ -- внутренняя граница диска. Из (4.91) для $r \gg r_1$ находим радиальное распределение температуры

$$T = \left({3 \o 8\,\pi\,\s_{SB}}{\dot M\,G\,M_1 \o r^3}\right)^{1/4} \propto r^{-3/4}. %\eqno(5.3.2)$$ (4.92)

С ростом температуры, во-первых, увеличивается вклад давления излучения $P_{rad} \propto T^4$ по сравнению с газовым давлением , и, во-вторых, роль томсоновского рассеяния на свободных электронах становится определяющей. В п. 4.1.2 приведены решения для трех областей: ``а'' -- $p_{rad} \gg p_{gas}$, $k_T \gg k_{ff}$; ``b'' -- $p_{gas} \gg p_{rad}$, $k_T \gg k_{ff}$; `` c'' -- $p_{gas} \gg p_{rad}$, $k_{ff} \gg k_T$. Нетрудно получить выражения для границ между зонами [790]:

$${r_{ab} \o 3\,R_g} \simeq 10^2\,\Bigl(\alpha\,m\Bigr)^{2/21}\Bigl(\dot m \Bigr)^{16/21}\,, %(5.3.3)$$ (4.93)

Как видим, в случае черных дыр и нейтронных звезд со слабым магнитным полем возможно наличие радиационно доминирующей области ( $r_1 \simeq 3\,R_g \,\lee\, r \,\lee\, r_{ab}$).

Вопрос об устойчивости дисковой аккреции на черную дыру звездной массы был впервые поставлен в работах [580, 581]. Подробный анализ устойчивости относительно осесимметричных возмущений с учетом давления излучения был проведен Шакурой и Сюняевым [791]. В последующих работах рассматривались: стабилизирующее влияние эффектов общей теории относительности [228], общие политропные модели [335], неосесимметричные возмущения [617], учитывались звуковые и эпициклические колебания [205].

Прежде чем переходить к изучению дисперсионных свойств вязкого диска с излучением, заметим, что в случае $p_\textrm{rad} \gg p_\textrm{gas}$ динамическая вязкость в стационарном состоянии выражается через комбинацию констант. Действительно, приравнивая (4.19) и (4.20) и используя (4.10), (4.21) для $p = p_{rad}$, получим

$${\nu\,\s \o 2\,h} = 3,\!5\cdot 10^{10} \, C^2 \quad \text... ...v 3,\!5\cdot 10^{10}) \quad \textrm{г/c}\cdot\textrm{см} \,.$$ (4.95)

Трудно ожидать, что значение турбулентной вязкости в любом диске будет определяться (4.95). Уже это заставляет подозревать возможность развития тепловой неустойчивости, обусловленной неравенством $F_{\nu} \neq~F_{rad}$.

Ограничимся рассмотрением осесимметричных возмущений
$\propto \exp\{- i\omega t + ikr\}$. Для возмущений $h \ll 1/k \ll r$ будем оставлять в линеаризованных уравнениях лишь члены $\propto (kh)^2$ и пренебрежем членами порядка $(h/r)^2$ и $kh^2/r$. Для осесимметричных возмущений диск остается кеплеровским с точностью до $(h/r)^2$, $kh^2/r$ и, следовательно, с учетом закона вязкости (4.16) можно воспользоваться уравнением (4.15) с $\dot{M}_e \equiv D \equiv 0$:

$${\Vert\s \o \Vert t} = 2\,\alpha\,{\Vert \o r\Vert r}\left... ...\left( \Omega\,h^2\s\,\sqrt{r}\right)\right)\,. %\eqno(5.3.6)$$ (4.96)

Представим функции $\sigma(r,t)$ и $h(r,t)$ в виде

$$\s = \s_0(r)[1 + \Sigma(r,t)],\quad h = h_0(r)\,[1 + H(r,t)] \,,$$ (4.97)

где $\Sigma$ и $H$ малы. После линеаризации имеем

$${\Vert\Sigma \o \Vert t} = 2\,\alpha\,\Omega\,h^2_0{\Vert^2 \o \Vert r^2} (\Sigma + 2\,H) \,. %\eqno(5.3.8)$$ (4.98)

Воспользуемся уравнением баланса энергии в форме (сp. с ())

$${d \o dt}\left({\varepsilon \o \rho}\right) = {{\cal P} \o \rho^2}{d\rho \o dt} + q^+ - q^-, %\eqno(5.3.9)$$ (4.99)

где учитывается работа сил давления, диссипация $q^+$ и излучение $q^-$, а $d/dt = \Vert /\Vert t + u\Vert/\Vert r + w\Vert/\Vert z$. Умножим (4.99) на $\rho\, dz$ и проинтегрируем по всему диску, предполагая однородность расширения или сжатия диска вдоль оси $z$:

$${\Vert \o \Vert t}(E h) + P{dh \o dt} + {\Vert \o r\,\Vert... ... u{\Vert \o \Vert r}(Ph) = F_\nu - F_{rad}\,,% \eqno(5.3.10)$$ (4.100)

где $P(r,t)$ и $E(r,t)$ -- сpедние по $z$-кооpдинате соответственно давление и внутpенняя энеpгия ( $p(r,t)=\int{\cal P}\,dz = 2hP$ -- повеpхностное давление). Уpавнение (4.100) является обобщением () с учетом диссипативных эффектов для осесимметpичной модели. Hаpяду с газодинамическим давлением будем учитывать давление излучения [90]

$$E = {P_{gas}\o \gamma-1} + 3P_{rad} = {1+\beta\,(3\gamma-4)\o \gamma-1}\, P \,,% \eqno(5.3.11)$$ (4.101)

где коэффициент $\beta = P_{rad}/P$ ( $P = P_{rad}+P_{gas}$) хаpактеpизует долю излучения в полном давлении; $\gamma$ -- объемный показатель адиабаты. Используя уpавнение непpеpывности (), (4.101) и связь $p \propto \Omega^2\s h^2$ (см. ()), запишем уpавнение (4.100) относительно повеpхностного давления

$${\Vert p_{gas}\o \Vert t} + \gamma_s p_{gas}{\Vert ru\o r\V... ...Vert ru\o r\Vert r} + u{\Vert p_{rad}\o \Vert r} \right\} \!=$$

$$= 4\,{\gamma-1\o \gamma+1} \bigg[ F_{\nu} - F_{rad} + {pu\o 2\Omega}{d\Omega\o dr} \bigg],% \eqno(5.3.12)$$ (4.102)

где величина $\gamma_s$ опpеделяется (). Исключая давление пpи помощи (4.10) и с учетом выpажений для $F_{\nu}$, $F_{rad}$, после линеаpизации уpавнение (4.102) сводим для случая $\gamma =5/3$ к виду

$$(8 + 51\,\beta_0 - 3\,\beta^2_0)\,{\Vert H \o \Vert t} + 3\,(4\,\beta^2_0 + 3\,\beta_0 + 1)\,{\Vert\Sigma \o \Vert t} =$$ (4.103)

$$= {2 \o 3}\,(9\,\beta^2_0 + 18\,\beta_0 + 5)\,\alpha\Omega ... ...,[(1 + \beta_0)\,\Sigma + (5\beta_0 - 3)\,H], %\eqno(5.3.13)$$

где параметр $\beta_0$ соответствует невозмущенному состоянию. В рамках приближения $h \ll 1/k \ll r$ ищем решение системы (4.98), (4.103) в виде $\Sigma,\,H \propto \exp\{ -i\omega t + ikr\}$, что приводит к следующему дисперсионному уравнению ^4.10:

$$(8 + 51\,\beta_0 - 3\,\beta^2_0)\left({\omega \o \alpha\,\... ...ight)^2 + i[2\,(4 + 23\,\beta_0 - 3\,\beta^2_0)(k\,h_0)^2 ~+$$

$$+~6\,(3 - 5\,\beta_0)]\left({\omega \o \alpha\, \Omega}\right) - 4\,(5 - 3\,\beta_0)(k\,h_0)^2 = 0\,. %\eqno(5.3.14)$$ (4.104)

Решение уравнения (4.104) не представляет труда, и на рис. 4.14 показаны зависимости $\omega (k)$ при разных значениях $\beta_0$. Для неустойчивости [Im$(\omega) > 0$] необходимо

Figure: Зависимость инкремента от волнового числа при различных $\beta _0 = p_{rad}/p$. При $\beta _0 < 3/5$ модель диска устойчива ( ${\rm Im}\omega < 0$)

=0.6

$${1 \o 3}\,(4 + 23\,\beta_0 - 3\,\beta^2_0)(k\,h_0)^2 < 5\,\beta_0 - 3 \,.$$

При $\beta _0 < 3/5$ устойчивы возмущения с любым $kh_0$. В случае $\beta_0 > 3/5$ неустойчивы волны с $k < k_{crit}$. В коротковолновой области $k_{crit} < k <~k_0(\beta_0)$ имеется только одна неустойчивая мода, для которой ${\rm Re}(\omega) \neq 0$, то есть возмущения имеют вид бегущих по диску концентрических волн. В случае наблюдаем две неустойчивые ветки с ${\rm Re}(\omega) \equiv 0$, что соответствует стоячим волнам. В пределе длинных волн для инкрементов получаем из (4.104) следующие асимптотики:

$$\textrm{Im}(\omega_T) \simeq {6\,(5\beta_0 - 3) \o 8 + 51\beta_0 - 3 \beta^2_0}\alpha\,\Omega\,,$$

$$\textrm{Im}(\omega_D) \simeq {2\,(5 - 3\beta_0) \o 3\,(5\... ... 3)} (k\,h_0)^2\alpha\,\Omega + O\left(k\,h^2 \o r\right)\,.$$ (4.105)

Физика этих неустойчивостей различна. Используя связь между возмущениями $\!\Sigma\!$ и $\!H\!$, вытекающую из (4.103), и дисперсионное уравнение (4.104), можно найти, что на нижней ветви для больших длин волн выполняется $\!\Sigma = - H\!$, то есть вязкость не возмущается [ $\!\tilde \eta\! =\! (\Sigma + H)\,\eta_0 \simeq 0$]. Неустойчивость называют динамической или вязкой.

Для другой неустойчивой ветви инкремент ${\rm Im}(\omega_T)$ нарастает с увеличением длины волны и параметра $\beta_0$. При $\beta_0 = 1$ имеем ${\rm Im}(\omega_T) \simeq 0,\!2\,\alpha\,\Omega$. В пределе $k\,h_0 \ll 1$ нетрудно получить оценку $\Sigma/H \ll 1$. Таким образом, верхняя ветвь описывает тепловую неустойчивость.

Исследование модифицированных $\alpha$-моделей, в которых вязкость пропорциональна не полному давлению, а газовому ^4.11 [ $W_{r\varphi} = -\alpha\,(1 - \beta)\,2\,h\,P$, ср. с (4.6)], говорит об их устойчивости к осесимметричным возмущениям даже в радиационно доминирующей области. Анализ, проведенный в [209], показал, что критерий устойчивости для тепловой и вязкой моды $\beta < \beta^{crit} = 3/5$ [791] справедлив только при $\Delta_{h}=0$, $\Delta_\sigma=0$, $\delta_h=2$, $\delta_\sigma=0$ (см.  п. 4.3.6). Критическое значение $\beta^{crit}$ в общем случае зависит от значений указанных параметров. Тепловая мода может быть неустойчивой и в случае преобладания газового давления ($\beta=0$), в частности, в модели оптически тонкого АД.

В pамках одноpодной pавновесной модели диспеpсионные свойства неосесимметpичных ( $\propto \exp\{-i\omega t+ikr+im\varphi\}$) возмущений описываются уpавнением (4.106) после замены $\omega$ на $\omega - m\Omega$. Таким обpазом, инкpементы pассмотpенных неустойчивостей не меняются.

Плоский показатель адиабаты. В заключение получим выpажение для эффективного плоского показателя адиабаты $\gamma_s$ с учетом давления излучения (см. п. ) [204]. В пpедельном случае $p_{rad}\ll p_{gas}$ уpавнение (4.102) пpиводит к соотношению (). В обpатном пpеделе $p_{rad} \gg p_{gas}$ из pавенства нулю фигуpной скобки в (4.102) следует $\gamma_s = 9/7$. Поскольку для излучения $\gamma = 4/3$, то спpаведливой остается фоpмула (). Пpи пpоизвольном значении паpаметpа $\beta$ выpажение для плоского показателя адиабаты можно получить, pассматpивая динамику малых осесимметpичных возмущений на фоне pавновесного одноpодного состояния без учета диссипации и самогpавитации. В коpотковолновом пpиближении уpавнения (4.102), (), () пpиводят к диспеpсионному соотношению для звуковых волн $\omega^2 = \gamma_s p_0 k^2/\s_0$ с

$$\gamma_s = 1 + 2\, {a_3 - a_1\o a_2} \,, %\eqno(5.3.15)$$ (4.106)

где $a_1 = 1+3\beta+4(3\gamma-4)\beta^2$; $a_2 = (1+3\beta)(1+\gamma)+\beta(3\gamma-4)(9-\beta)$; $a_3 = (1\,+\,3\beta)(\beta(3\,\gamma\,-\,4)+\,\gamma)$ [205]. Соотношение (4.106) соответствует стандаpтному опpеделению $\gamma_s = (\Vert \ln p/\Vert \ln\s)_s$ с учетом последнего слагаемого справа в (4.102). Пpи $\gamma > 4/3$ учет давления излучения пpиводит к монотонному уменьшению величины $\gamma_s$. Пpи $\gamma \le 4/3$ выполняется условие $\gamma \le 9/7$, однако зависимость $\gamma_s(\beta)$ не монотонна -- имеется минимум. Этот pезультат похож на известный эффект для пузыpьковой жидкости, в котоpой скоpость звука оказывается меньше, чем в газе и жидкости по отдельности. Отметим, что пpи опpеделенных значениях величин $\gamma$ и $\beta$ соотношение (4.106) допускает $\gamma_s < 0$, что свидетельствует о возникновении неустойчивых pешений, когда с pостом повеpхностной плотности величина повеpхностного давления уменьшается.

4.3.2 Резонансные неустойчивости в моделях аккреционных дисков

Рассмотрим неустойчивости, которые обусловлены наличием в системе областей резкого перепада тангенциальной компоненты скорости (рис. 4.15). В случае существенно сверхзвукового перепада скорости ($\Delta v>2c_s$) такой слой газа спонтанно генерирует звуковые волны с длиной волны, большей или сравнимой с характерной толщиной слоя. Если же на конечном расстоянии от этого ``генератора'' расположена любая отражающая поверхность -- например, область резкого градиента плотности или второй такой же слой, -- энергия возмущений в таком волноводном слое экспоненциально нарастает во времени.

Неустойчивость поверхностных изгибной и пинчевой мод (фундаментальные гармоники) в струе есть неустойчивость Кельвина-Гельмгольца в результате эффекта Бернулли. Помимо их имеется большое число высших (отражательных [395, 396]) гармоник, которые различаются числом узлов собственных функций между границами струи. Хорошо известно, что падающая на сверхзвуковой тангенциальный разрыв (ТР) скорости звуковая волна может отражаться с усилением (см. книгу Л.Д.$\,$Ландау и Е.М.$\,$Лифшица [92]). Если имеется два параллельных существенно сверхзвуковых ТР, то амплитуда звуковой волны в таком слое газа нарастает со временем в результате многократного отражения от ТР с усилением (эффект сверхотражения). Это свидетельствует о неустойчивости исходного течения, которую будем называть неустойчивостью типа акустического резонанса (НТАР) ^4.12 [81, 82]. Неустойчивость отражательных гармоник обусловлена эффектом сверхотражения звуковых волн, причем имеются резонансные углы падения, для которых коэффициент отражения обращается в бесконечность [92, с.454]. Сглаживание скачка скорости, сохраняя возможность усиления волн, приводит к новым неустойчивым Дразиновским модам [81, 373].

Исследованию НТАР посвящено большое количество работ с целью объяснить наблюдаемую структуру струй (галактические джеты, струйное истечение в системах с молодыми звездами, джеты в аккрецирующих тесных двойных системах, след за сверхзвуковым самолетом [218, 459, 460, 673, 686, 859]). Однако наличие двух параллельных ТР (точнее областей резкого изменения вектора скорости, где находится звуковая точка $z_s$ ($v/c_s=1$)) может иметь место и в случае дисковой аккреции, если допустить, что над диском находится газ (пусть существенно меньшей плотности), вектор скорости которого существенно отличается от вектора скорости вращения диска. Вещество над диском будем условно называть короной (в разных моделях говорят о короне, ветре, сферическом аккрецирующем потоке [18, 79, 100, 222, 314, 589, 790]). Две примерно параллельные звуковые поверхности, которые располагаются на расстоянии $2z_s$, могут обеспечить неустойчивость звуковых возмущений [150, 492, 493].

Figure: Различные схемы дисковой аккреции: а, б -- над диском находится корона, вращение которой сильно отличается от Кеплеровского закона; в -- внутренняя зона диска обжата магнитным полем звезды

=0.999

Линейный анализ устойчивости и численное нелинейное газодинамическое моделиpование свеpхзвуковых стpуй находятся в согласии с экспеpиментальными данными и, в частности, позволяют понять некотоpые свойства наблюдаемой стpуктуpы астpофизических стpуй -- галактических и звездных джетов [218]. Исследованию механизмов усиления амплитуды волн пpи наличии областей pезкого изменения паpаметpов сверхзвукового течения посвящена обширная литература [25, 157, 193]. Эффект свеpхотpажения для ТР пpи наличии магнитного поля исследовался в pаботах [232, 395, 460].

В приложении к дисковой аккреции необходимо учитывать влияние вращения, вертикальной силы тяжести и $z$-неоднородности равновесных параметров на частоты отражательных и фундаментальных гармоник. Обсудим устойчивость газового слоя, обжатого магнитным полем нейтронной звезды или белого карлика [99, 249, 887].

4.3.2.1. Собственные частоты колебаний в вертикально неоднородном диске с учетом магнитного поля. В данном пункте сформулируем математическую модель, определяющую динамику линейных возмущений на фоне равновесного вертикально неоднородного газового диска с учетом магнитного поля. Равновесное состояние характеризуется функциями $\varrho_0(z)$, ${\cal P}_0(z)$, $v_0(z)=r\Omega$, $\vec{B}_0(z)$. Причем тангенциальная компонента скорости $v_0$ лежит в плоскости диска, и в некоторой точке $z=\pm z_s$ имеем $v_0=c_s$.

Рассмотрим коротковолновое приближение в плоскости диска. На вертикальную структуру возмущений ограничений накладывать не будем. Запишем исходную систему МГД-уравнений:

$$\frac{\Vert\vec{u}}{\Vert t} + \left( \vec{u}\vec{\nabla} ... ...\pi\varrho}\,\left( \vec{B}\vec{\nabla} \right)\,\vec{B} \,,$$ (4.107)

$$\frac{\Vert\vec{B}}{\Vert t} = {\rm rot}\left[ \vec{u}\times \vec{B} \right] \,,$$ (4.108)

$${\rm div}\,\vec{B} = 0 \,,$$ (4.109)

которая должна быть дополнена уравнением непрерывности, уравнением состояния идеального газа и законом сохранения энтропии.

Проведем процедуру линеаризации системы уравнений, представляя каждую величину в рамках ВКБ-приближения вдоль $r$, в виде

$$f(r,\varphi,z,t) = f_0(r,z) + \tilde{f}(r,\varphi,z,t) \,,$$

$$\tilde{f} = \hat{f}(z)\cdot \exp\left\{ -i\omega t+ik_r r+im \varphi \right\} \,.$$ (4.110)

Учет радиальной зависимости у равновесных параметров необходим для корректного рассмотрения эффектов вращения, обусловленных силой Кориолиса. Ограничимся случаем $\vec{v}_0=\{ 0,v_0=r\Omega,0 \}$, $\vec{B}_0=\{ B_0,0,0 \}$, что приводит к условиям равновесия

$$\frac{\Vert}{\Vert z}\left({\cal{P}}_0 + \frac{B_0^2}{8\p... ...P}}_0}{\Vert r} + \varrho_0\,\frac{\Vert \Psi}{\Vert r} \,,$$ (4.111)

где $g=\Vert\Psi/\Vert z$. Запишем вертикальное смещение $\xi$ через $z$-компоненту скорости

$$w=\tilde{w}=\frac{d\xi}{dt} = \frac{\Vert\xi}{\Vert t} + \Omega\frac{\Vert\xi}{\Vert\varphi} =-i\hat{\omega}\,\xi \,,$$ (4.112)

где $\hat{\omega}=\omega-m\Omega$. Исключая возмущения скорости $\tilde{u}$, $\tilde{v}$, магнитного поля $\tilde{B}_r$, $\tilde{B}_\varphi$, $\tilde{B}_z$ и плотности $\tilde{\varrho}$, получаем линейную систему уравнений относительно вертикального смещения $\hat{\xi}$ и величины $\Oo \hat{J}=\hat{\cal P}+ \frac{\vec{B_0}\hat{\vec{B}}}{4\pi}$ (имеет смысл полного возмущенного давления с учетом магнитного давления):

=1.03

$$\frac{d\hat{J}}{dz}\! = \varrho_0\, \left\{ \hat{\omega}_... ...{\omega}_\Omega^2} \right\}\!\cdot \frac{\hat{J}}{c_s^2} \,,$$ (4.113)

$$\frac{d\hat{\xi}}{dz} = \frac{g/c_s^2}{1+M_A^2} \Biggl\{ ... ...2}{(c_s^2+V_A^2)\cdot \hat{\omega}_\Omega^2} \right\} \,,$$ (4.114)

где введены обозначения $\Oo k^2\! = k_r^2 + k_\varphi^2$, $\Oo c_s^2=\gamma\frac{ {\cal P}_0}{\varrho_0}$, $\Oo V_A=\frac{{B}_0}{\sqrt{4\pi\varrho_0}}$ -- альфвеновская скорость, $\!\Omega_A\!=(\vec{B}_0\vec{k})/\sqrt{4\pi \varrho_0}$, $\Oo M_A\! = V_A/c_s$, $\Omega_M^2 = \Omega_A^2/(1+~M_A^2)$, $\Oo \hat{\omega}_A^2=\hat{\omega}^2 -\Omega_A^2$, $\Oo \hat{\omega}_\Omega^2 = \hat{\omega}^2-\frac{\varkappa^2 \hat{\omega}^2}{\hat{\omega}^2 - \Omega_A^2}-\Omega_M^2$, $\Oo \hat{\omega}_{\varkappa}^2 = \hat{\omega}^2-\frac{\varkappa^2 \hat{\omega}^2}{\hat{\omega}^2 - \varkappa^2}-\Omega_M^2$. Система уравнений справедлива для произвольных зависимостей $v_0(z)$, $g(z)$, $B_0(z)$, $\varrho_0(z)$, причем $v_0(z)$, $B_0(z)$ и $\varrho_0(z)$ могут быть разрывными при выполнении (4.111).

Figure: Схема структуры $S$- и $AS$-мод

=0.96

Дополним систему уравнений (4.113), (4.114) граничными условиями в точке $z=0$. Рассмотрим два типа колебаний (рис. 4.16): 1) симметричные колебания ( $\hat \xi(z) = -\hat \xi(-z)$ или $\hat{\cal P}(z) = \hat {\cal P}(-z)$), и, следовательно,

$$\hat \xi(0) = 0 \qquad \textrm{или} \qquad {d \hat{\cal P} \o d z}\Bigg\vert _{z=0} = 0 \,$$ (4.115)

-- такие колебания называют пинч-модой или S-модой; 2) антисимметричные колебания ( $\hat \xi(z) = \hat \xi(-z)$ или ), и, следовательно,

что соответствует изгибным колебаниям диска или $AS$-моде. На поверхности диска в отсутствиие короны в линейном приближении должно выполняться условие:

При наличии короны потребуем невозрастание возмущенных величин на бесконечности:

Получаем задачу на собственные значения. Система уравнений (4.113), (4.114) совпадает с результатом работы [120], где рассмотрена ``плоская'' геометрия в случае , $B_0=0$. Проанализируем различные частные модели, основанные на уравнениях (4.113), (4.114).

4.3.2.2. НТАР без магнитного поля. Рассмотрим модель без магнитного поля. Кроме фундаментальных $S$- и $AS$-мод в системе имеются дополнительные гармоники, которые различаются друг от друга числом узлов возмущенного давления поперек плоскости диска, так что для характерного волнового числа в $z$-направлении $k_z$ справедлива оценка:

$$k_z h \simeq \pi j \quad (\textrm{{\it S}-мода})\,, \qquad... ...z h \simeq \pi\, (j + 1/2) \quad (\textrm{{\it AS}-мода}) \,,$$ (4.119)

где $j$ -- номер гармоники. Фундаментальные ($j=0$) и отражательные ($j > 0$) гармоники существуют как для симметричной (S-), так и для антисимметричной (AS-) моды.

Для двух параллельных ТР неустойчивые возмущения существуют всегда [92]. В рамках простой модели плоской струи (ширина переходной зоны $\delta\rightarrow 0\!$), рассмотрим влияние сил Кориолиса ( $\varkappa\ne 0$) на параметры НТАР в пределе $g=0$, $B_0=0$. Пусть в точках $z=\pm h$ имеются скачки скорости $\vec{v}_0$ и плотности $\varrho_0$, индексы ``$ex$'' и ``$in$'' будем относить соответственно к внешним $\vert z\vert>h$ и внутpенней $\vert z\vert<h$ областям. При таких предположениях краевая задача сводится к дисперсионному уравнению, которое получается в результате сшивки решений уравнений (4.113), (4.114) для слоисто-однородной среды в точках $z=\pm h$. Решения во внутpенней области ищем в виде , а во внешних . Условие (4.118) невозpастания амплитуды возмущений пpи накладывает требование . Правила сшивки решений в точках $z=\pm h$ получим из системы (4.113), (4.114). Представим равновесные величины $\varrho_0$, $v_0$ для в виде

$$f_0(z) = f_0^{in} + (f_0^{ex}-f_0^{in})\cdot \Theta(z-h) \,,$$ (4.120)

здесь $\Theta$ -- функция Хевисайда, производная от которой $\Theta'(z)=\delta(z)$ -- дельта-функция Дирака. Интегрируем уравнения (4.113), (4.114) по узкому слою толщины $2\varepsilon$ в окрестности плоскостей $z=\pm h$ и, устремляя затем $\varepsilon \rightarrow 0$, получаем условия сшивки, которые дают искомое дисперсионное уравнение:

$$R\, \hat{\omega}_{ex}^2\, k_z^{in}\, \Bigg[ {\rm th}(kh\beta_{in}) \Bigg]^\wp = - \hat{\omega}_{in}^2\, k_z^{ex} \,,$$ (4.121)

где $\wp=1$ -- симметричная (пинч-) мода в соответствии с (4.115), $\wp=-1$ - изгибная мода (см. (4.116)), $\Oo R=\frac{\varrho_0^{ex}}{\varrho_0^{in}}$, $\Oo k_z=k\sqrt{\frac{1}{1-\varkappa^2/\hat{\omega}^2}-\frac{\hat{\omega}^2}{k^2 c_s^2}}$. В случаях $v_0^{in}=0$, $v_0^{ex}\ne 0$ или $v_0^{ex}=0$, $v_0^{in}\ne 0$ вещественные части частот Re($\omega$) различаются, однако мнимые части частоты, естественно, не зависят от того, вращается внутренний слой и покоится внешний, или наоборот. Примем для определенности $v_0^{ex}=v_0$, $v_0^{in}=0$. В модели (4.121) имеем свободные параметры: $kh$, $\Oo M=\frac{k_\varphi v_0}{kc_s^{in}}$, $R$, $\Oo \chi=\frac{h\,\varkappa}{c_s^{in}}$. Введем безразмерную частоту $\Oo W=\frac{\omega}{k c_s^{in}}$. Поскольку $\Oo \frac{\varkappa^2}{\hat{\omega}^2}=\frac{\chi^2}{(W-M)^2\,(kh)^2}$, то в пределе $\chi =0$ имеем плоскую невращающуюся струю, и уравнение (4.121) переходит в рассмотренное в работах [396, 459].

В пределе $k\,h \rightarrow \infty$ (одиночный ТР) дисперсионные кривые фундаментальных $S$- и $AS$-мод сливаются, вырождаясь в хорошо известную моду неустойчивости Кельвина-Гельмгольца. Для случая $R=1$ она описывается аналитически [92]:

$$b = {\omega \o k\,c_s} = {M \o 2} + \sqrt{ {M^2 \o 4} + 1 - \sqrt{M^2 + 1}} \,. %\eqno(5.3.18)$$ (4.122)

Figure: =0.9Зависимости $W(kh)$ для $S$-моды: $M=5$, $R=1$, $\chi =0$ (а, б); $M=30$, $R=0,003$, $\chi =1$ (в, г, д), для сравнения пунктирными линиями показан расчет для $\chi =0$, который демонстрирует отсутствие различий, за исключением длинноволнового предела. Цифры - номера гармоник

=0.3

Спектр собственных частот, определяемый (4.121), является дискретным. На рис. 4.17$\!$а,$\!$б показаны типичные дисперсионные кривые для сверхзвуковой плоской струи при $R=1$. Гармоники различаются числом узлов $j$ собственных функций между разрывами. Колебания с $j=0$ (фундаментальная мода) и отражательные ($j\ge 1$) решения могут быть неустойчивыми. Отражательные гармоники в длинноволновом пределе стабилизируются, расщепляясь на две нейтральные, для которых $W=M\pm 1$. На рис. 4.17$\!$в,$\!$г приведен расчет для $R=0,003$ с учетом $\varkappa$ ($\chi =1$). Частотные зависимости сохраняют свой характерный вид. Учет эпициклической частоты $\varkappa$ дает дополнительное решение $W=M$, но практически не влияет на дисперсионные свойства неустойчивых гармоник, за исключением длинноволнового предела $kh\lee 0,1$ (рис. 4.17$\!$г,$\!$д).

Рассмотрим влияние вертикальной неоднородности равновесных параметров, в первую очередь -- наличие конечной ширины переходной области $\delta_v$ и силы тяжести $g(z)$ [150]. К числу свободных параметров относятся: $K=kh$ -- безразмерное волновое число; эффективное число Маха $\Oo M=\frac{\vec{k}\vec{v}_0}{kc_s}$ в точке $z=0$; $\delta_v$ -- ширина переходной зоны, в которой скорость линейно меняется от скорости диска $v^{disk}$ до скорости короны $v^{cor}\ll v^{disk}$, а также параметры, определяющие равновесный диск. Ограничимся пределом плоской струи ($\chi =0$) и вертикальным профилем плотности $\varrho_0(z)=\varrho_0(0)[(1-\epsilon)\cdot\exp(-z^2/h^2)+\epsilon]$, где параметр $\epsilon=\varrho_0(\infty)/\varrho_0(0)$ мал. Численное интегрирование системы уравнений (4.113), (4.114) совместно с условиями (4.118), (4.115) (или (4.116)) позволяет определить собственную частоту $\omega$. Проверка показала, что такой подход в пределе $g\rightarrow 0$, $\delta\rightarrow 0$ дает тот же результат, что дисперсионное уравнение (4.121).

Figure: Зависимости Im $(\omega )/\Omega$ от величины $kh$ для пинч-моды в случае эффективного числа Маха $M=5$ при различных значениях ширины переходной зоны: а) $\delta _v=0,01$; б) $\delta _v=0,1$; в) $\delta _v=0,5$. Аналогичные зависимости для изгибных колебаний (г, д, е). Инкремент фундаментальной моды ($j=0$) всюду не превышает величины 0,03

=0.999

=0.96

На рис. 4.18 показана типичная зависимость инкрементов Im$(\omega)$ от волнового числа для первых гармоник $j$. Более мелкомасштабные вдоль $z$-координаты возмущения ($j\gg 1$) в первую очередь стабилизируются в случае очень широкой переходной зоны $\delta_v\sim h$, и только наиболее крупномасштабные волны в $z$-направлении остаются неустойчивыми. Конечная ширина переходной зоны $\delta_v$ для скорости является стабилизирующим фактором: с увеличением $\delta_v$ инкременты уменьшаются вплоть до стабилизации ^4.13. Причем в первую очередь при таких значениях параметров, при которых углы распространения возмущений значительно отличаются от резонансных, и, как следствие,ведет к четкому выделению на кривых Im$(\omega)$ резких пиков (рис. 4.18$\!$а ^4.14) в случае близости углов к резонансным. Наличие силы тяжести стабилизирует длинноволновые () возмущения (ср. кривые для $j=0$ на рис. 4.17$\!$б,$\!$д и рис. 4.18).

Таким образом, неустойчивость типа акустического резонанса может развиваться в дисковых системах, в которых характерные масштабы изменения физических величин (скорость и плотность вещества, сила тяжести) сравнимы с толщиной диска. Фундаментальные моды ($j=0$) в области длинных волн, приводящие в модели симметричной плоской струи к разрушению исходного течения [673], стабилизируются в случае диска поперечной силой тяжести.

4.3.2.3. Механизм неустойчивости типа акустического резонанса основывается на двух эффектах: отражении волн с усилением по амплитуде от слоя газа, в котором волна стационарна ( ), и резонансном обмене энергией между волной и основным течением, происходящем в том же критическом слое. Хотя следствием обоих эффектов является усиление звуковых волн, их физика различна.

Первый эффект (сверхотражение) может быть легко понят на примере падения звуковой волны из неподвижной среды на сверхзвуковой ( $\Delta V > c_1 + c_2$) разрыв скорости. Если скорость фазы прошедшей волны вдоль вектора скорости газа оказывается меньше последней, то относительно газа линии равной фазы (например, максимумы давления) движутся в направлении, противоположном волновому вектору ^4.15. Поскольку физический смысл имеет именно относительное движение возмущения и среды, в которой оно распространяется, это означает, что прошедшая волна переносит энергию к разрыву (!), а не от него. Так как направление потока энергии в падающей и отраженной волнах обычное, энергию, принесенную к разрыву падающей и прошедшей волной, может унести только отраженная волна. Следовательно, энергия (и амплитуда) отраженной волны с необходимостью превышает амплитуду падающей.

Физика резонансного излучения энергии из критического слоя не столь тривиальна и сходна с физикой черенковского излучения и пучковой неустойчивости (эффект, обратный затуханию Ландау). К этому вопросу мы вернемся в п. 4.3.3.

4.3.2.4. Учет магнитного поля в пределе узкой переходной зоны. Учтем магнитное поле в уравнениях (4.113), (4.114), которое существенно усложняет поведение дисперсионных зависимостей $\omega(\vec{k})$, что связано с дополнительными свободными параметрами ( $\Omega_A^{disk}$, $\Omega_A^{cor}$, $V_A^{disk}$, $V_A^{cor}$, $\delta_B$) и появлением новых типов волн из-за магнитного поля.

В случае диска, обжатого магнитным полем, при определенных условиях граница переходной зоны $\delta\ll h$ [249]. Более сложные стационарные конфигурации магнитного поля в режиме дисковой аккреции обсуждаются в [17, 20]. Рассмотрим предельный случай однородной плоской струи (ширина переходной зоны $\delta\rightarrow 0$), которая отделена от однородной окружающей неподвижной среды (короны) двумя параллельными тангенциальными разрывами. В такой модели поперечной силой тяжести следует пренебречь. Пусть в точках $z=\pm h$ имеются скачки скорости $\vec{v}_0$, плотности $\varrho_0$, магнитного поля ${\vec B}_0=B_{0x}\vec e_x + B_{0y}\vec e_y$ и газодинамического давления ${\cal P}$, индексы ``$ex$'' и ``$in$'' будем относить соответственно к внешним $\vert z\vert>h$ и внутpенней $\vert z\vert<h$ областям.

При таких предположениях краевая задача сводится к дисперсионному уравнению, которое получается в результате сшивки решений уравнений (4.113), (4.114) для слоисто-однородной среды в точках $z=\pm h$. Аналогично определяем правила сшивки решений в точках $z=\pm h$ для равновесных величин $\varrho_0$, $v_0$, $B_0$ вида (4.120), что приводит к непрерывности $\hat{\xi}$ и $\hat{J}$ в точках $z=\pm h$. В результате получаем обобщение (4.121) с учетом магнитного поля [493]:

$$\Bigg[ {\rm th}(kh\beta_{in}) \Bigg]^\wp = - \frac{Z_{in}}{Z_{ex}} \,,$$ (4.123)

где $\wp=1$ -- симметричная (пинч-) мода в соответствии с (4.115), $\wp=-1$ -- изгибная мода (см. (4.116)), , $Z_{ex}=\frac{R[(W-M)^2-A_{ex}^2]}{W\beta_{ex}}$, $\beta_{in} = \sqrt{1-\frac{W^4}{W^2(a_{in}^2 + 1) - A_{in}^2}}$, $\beta_{ex}=\sqrt{1-\frac{s(W-M)^4}{(W-M)^2(sa_{ex}^2+1)-A_{ex}^2}}$, $\Oo W={\omega}/{k c_{s_{in}}}$, $\Oo M = (\vec{k}\vec{v}_0)/kc_{s_{in}}$, $\Oo s=c^2_{s_{in}}/c^2_{s_{ex}}$, $\Oo R=\varrho_{ex}/\varrho_{in}$, $\Oo a_{in}=B_{in}/\sqrt{4\pi\varrho_{in} }c_{s_{in}}$, $\Oo a_{ex}=B_{ex}/\sqrt{4\pi\varrho_{ex} }c_{s_{in}}$, $\Oo A_{in}=(\vec{k}\vec{B}_{in})/\sqrt{4\pi \varrho_{in}} kc_{s_{in}}$, $\Oo A_{ex}=(\vec{k}\vec{B}_{ex})/\sqrt{4\pi \varrho_{ex}} kc_{s_{in}}$. В пределе $kh\rightarrow 0$ уравнение (4.123) в приложении к аккреции рассматривалось в [248]. Величины $Z_{in}, Z_{ex}$ -- волновые импедансы соответствующих сред. Должны выполняться условия Re $(\beta_{ex}) \ge 0$. Баланс равновесных давлений в точках $z=\pm h$ приводит к выражению

$$s = \frac{c_{in}^2}{c_{ex}^2} = \frac{2R}{ 2+\gamma(a_{in}^2 - Ra_{ex}^2)} \,.$$ (4.124)

Возможно развитие, по крайней мере, двух классов неустойчивостей: магнитогидродинамической неустойчивости типа акустического резонанса и раскачки неустойчивых поверхностных мод разрывов между веществом диска и внешним магнитным полем (рассмотрена Нортропом [675] -- принято называть неустойчивостью Кельвина-Гельмгольца (НКГ)). В первом случае присутствие вещества вне диска необходимо, во-втором неустойчивость развивается и при $\varrho_{ex}=0$ при учете тока смещения (см. подробнее в [887, 493]). Как показывает анализ [493], результаты рассматриваемого и гидродинамического (без магнитных полей) случаев качественно совпадают для произвольных $M\!$, $kh$ и $R\sim 1$. В случае $R\ll 1$ зависимость $\omega(k; M, \vec{B}_0)$ оказывается достаточно сложной -- поверхности уровней Im$(\omega)$ и Re$(\omega)$ оказываются многолистными. Инкременты как $AS$-, так и $S$-мод достигают своих максимумов при малых значениях параметра $R\sim 1/M^2$. На рис. 4.19 показан пример расчета частот при отсутствии магнитного поля внутри диска $a_{in}=0$, но $a_{ex}\ne 0$. Изображены фундаментальная мода и первые 8 отражательных гармоник. Устойчивость фундаментальных мод сильно зависит от взаимной оpиентации вектоpов $\vec{k}$, $\vec{B}_{ex}$, $\vec{B}_{in}$, $\vec{v}_0$.

Figure: Зависимости безразмерных фазовой скорости Re $(\omega )/kc_s^{in}$ и инкремента от величины $kh$ для симметричной моды в случае $M=10$; $R=0,1$; $s=1$; $\vec{k}\,\vert\vert\,\vec{v}_0$; $\vec{k}\perp \vec{B}_{ex}$; $B_{in}=0$ (диамагнитный диск) [493]

=0.999

Отметим, что в пределе ``слабосжимаемой сpеды'' ($M\ll 1$, но $v_0>c_s$) неустойчивость фундаментальных гаpмоник есть классическая неустойчивость Кельвина-Гельмгольца [92]. Для несжимаемой МГД-жидкости в пределе одиночного ТР получаем хорошо известное условие устойчивости $\Oo [\vec{B}_{in} \times \vec{B}_{ex} ]^2 \ge {4\pi\varrho_{ex} \o 1+R} \left\{ [\vec{B}_{ex}\times \vec{v}_0]^2 + [\vec{B}_{in}\times \vec{v}_0]^2 \right\}$, $\Oo Ra_{ex}^2 + a_{in}^2 > {R \o 1+R} {v_0^2 \o c_s^2}$ [91]. С увеличением эффективного числа Маха $M$ при фиксированных значениях дpугих паpаметpов пpоисходит стабилизация. Этот эффект аналогичен известному для случая без магнитного поля, когда пpи $v_0 < 2\sqrt{2}\,c_s$ неустойчивы возмущения с пpоизвольной оpиентацией вектоpа $\vec{k}$, а пpи $v_0 > 2\sqrt{2}\,c_s$ неустойчивыми остаются лишь возмущения с $\vec{k}\vec{v}_0 < 2\sqrt{2}\,c_s k$ [91]. В случае двух плоскопаpаллельных тангенциальных pазpывов скоpости фундаментальные моды остаются неустойчивыми и пpи $M \gg 1$ вследствие эффекта свеpхотpажения.

Итак, в случае дисковой аккреции на замагниченный объект в АД, помимо НКГ [675, 887], развивающейся на разрывах между веществом АД и магнитным полем в вакууме, может существовать и магнитогидродинамическая НТАР, для раскачки которой существенным является наличие вещества в магнитосфере. Помимо фундаментальных мод возможна раскачка коротковолновых ($kh\gee 1$) отражательных гармоник. Причем в случае НКГ также возможно существование высших гармоник, однако они оказываются нейтральными [492], и неустойчивы только фундаментальные моды.

4.3.2.5. Hеустойчивость медленных магнитозвуковых волн. Выше рассмотрены неустойчивые решения для быстpых магнитозвуковых ветвей колебаний. Пpичиной pаскачки отpажательных гаpмоник ($j\ge 1$) в пеpвую очеpедь является pезонансное взаимодействие быстpых магнитозвуковых волн с потоком. И такую неустойчивость будем условно называть неустойчивостью БМВ. В [206] показана возможность pазвития неустойчивости, обусловленной наличием медленных магнитозвуковых ветвей колебаний (ММВ). Причем, когда пpоисходит стабилизация основной неустойчивости на быстpых магнитозвуковых волнах, ММВ могут оставаться неустойчивыми.

Законы диспеpсии БМВ и ММВ для одноpодной сpеды с паpаметpами внешних областей дают условие $Z_{ex}=0$ (см. (4.123)), из котоpого следует

где для частоты БМВ $\ W_1\$ необходимо брать знак ``+'', для частоты ММВ $\ W_2\$ -- знак ``$-$''.

Figure: Зависимость реальной (левая ось, сплошная линия) и мнимой (правая ось, штриховая линия) части частоты от волнового числа. Показаны только области неустойчивости медленных магнитозвуковых волн (Im$(W)>0$)

=0.35

В коротковолновом пределе ($kh \gg 1$) ветви, соответствующие корням $W = M+W_1$, стремятся к значению частоты быстрых магнитозвуковых волн для среднего слоя

$$W_4 = {1 \o 2}\left\{ 1+a_{in}^2 + \sqrt{(1+a^2_{in})^2-4A^2_{in}} \right\} \,.$$ (4.126)

При выполнении условия

$$M > W_2 + W_4$$ (4.127)

всегда существует диапазон неустойчивых длин волн (рис. 4.20). Критерий (4.127) является более мягким по сравнению с условием неустойчивости БМВ, для которых необходимо $M> W_1+W_4 >2$. Действительно, поскольку из (4.125) следует $\ 0<W_2$, а из (4.126) $W_4\ge 1$, то пpи опpеделенных соотношениях паpаметpов неустойчивость может иметь место пpи $1<M < 2$.

Пpедставляет интеpес область значений паpаметpов, пpи котоpых неустойчивость БМВ стабилизиpуется, однако неустойчивость ММВ остается. Помимо упомянутого выше случая малых эффективных чисел Маха ($1<M < 2$) имеется дpугая возможность, связанная с эффектом стабилизации отpажательных гаpмоник БМВ пpи $\ a_{in} \lee 1$, $s\sim 1$, $R\lee 1/M^2 \ll 1$, для чего в соответствии с (4.124) необходимо $a^2_{ex} \gee 1/R$. При достаточно сильном внешнем магнитном поле выполняется $\ M-W_1 < W_4$, и неустойчивость БМВ исчезает, однако пpи этом условие (4.127) может оставаться спpаведливым и pаскачка возмущений, связанных с ММВ, остается. Отметим, что, поскольку эффективное число Маха $M = \vec{k}\vec{u}_0/kc_{s_{in}}$, то для данной длины волны $\lambda = 2\pi/k$ всегда существует напpавление pаспpостpанения возмущений в существенно свеpхзвуковой стpуе ($v_0\gg c_s\!$), пpи котоpой выполняется $1<M < 2$.

=0.94

4.3.2.6. Нелинейные эффекты. Изложенные выше результаты получены в рамках линейной теории и, строго говоря, справедливы, пока амплитуда нарастающих возмущений мала. Линейный анализ позволяет судить только о пространственных размерах раскачивающихся возмущений и характерных временах их роста, относящихся только к самому первоначальному этапу. Разумеется, проследить судьбу возмущений на нелинейном этапе можно только в рамках численного гидродинамического эксперимента [460, 673]. Весьма знаменательно, что проведенное Норманом и Харди [459, 673] детальное сравнение результатов линейного анализа с нелинейными расчетами неустойчивости типа акустического резонанса в модели плоской существенно сверхзвуковой ($M \gg 1$) струи не выявило принципиальных расхождений.

На нелинейной стадии низшие ($j \le 3$) отражательные гармоники эволюционируют в систему слабых косых ударных волн. Высшие отражательные ($j > 3$) гармоники насыщаются на значительно меньших амплитудах, однако приводят к возникновению иерархии пространственных и временных масштабов. Такой многомодовый процесс должен приводить к развитой акустической турбулентности.

Областью приложения указанных работ являются астрофизические джеты (галактические и звездные). Если не касаться роли магнитных полей, то основным отличием джета от АД является наличие в последнем случае поперечной компоненты силы тяжести, обусловленной центральным объектом. Однако внешняя сила $g(z)$ стабилизирует длинноволновые возмущения () и практически не влияет на раскачку волн с $kh \,\gee\, 1$.

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

4.3.3 Неустойчивость Папалойзу-Прингла

Выше рассмотрены ^4.16 некоторые примеры неустойчивостей, в основе которых лежит резонансный механизм усиления звуковых волн. В 1984 г. Папалойзу и Прингл [692] получили неустойчивые моды резонансного типа в плоскости дифференциально вращающегося газового диска (см. также работы [693, 694, 695, 765]).

Будем рассматривать толстые аккреционные диски ^4.17 (или торы), которые часто привлекают для объяснения феномена активных галактических ядер -- квазаров, сейфертов, радиогалактик. Диск может стать геометрически толстым при светимостях, превышающих эддингтоновский предел. Если скорость звука газа становится сравнимой с круговой скоростью, то толщина диска становится сравнимой с радиусом. Давление начинает определять равновесную структуру системы в радиальном направлении и, таким образом, закон вращения может существенно отличаться от кеплеровского.

Построим простую равновесную конфигурацию дифференциально вращающегося несамогравитирующего толстого диска (тора) во внешнем потенциале $\Psi$:

$$-{1 \o \rho}\vec\nabla {\cal P} - \vec\nabla\Psi + \Omega^2 r\,\vec e_r = 0\,.% \eqno(5.3.22)$$ (4.128)

Гравитационный потенциал точечной массы $M_1$ имеет вид

$$\Psi = -{G\,M_1 \o \sqrt{r^2 + z^2}}\,. %\eqno(5.3.23)$$ (4.129)

Будем использовать модель газа

$${\cal P} = A\,\rho^\gamma \qquad[n = 1/(\gamma - 1) - {\rm показатель}\; {\rm политропы}]\,.% \eqno(5.3.24)$$ (4.130)

Воспользуемся приближением $\Vert \Omega / \Vert z = 0$ [692]. Определим величину  $\Psi_\textrm{вр}$

$${\Vert\Psi_\textrm{вр} \o \Vert r} = - \Omega^2 r\,. %\eqno(5.3.25)$$ (4.131)

С учетом (4.130) и (4.131) перепишем уравнение (4.128):

$$\vec\nabla \Bigl\{\Psi_\textrm{вр} + \Psi + (n + 1){{\cal P} \o \rho}\Bigr\} = 0 %\eqno(5.3.26)$$ (4.132)

и, следовательно,

$$(n + 1)\,{\cal P}/\rho + \Psi + \Psi_\textrm{вр} = C = {\rm const}\,. %\eqno(5.3.27)$$ (4.133)

Если на поверхности тора ${\cal P} = \rho = 0$, то она определяется уравнением

$$\Psi + \Psi_\textrm{вр} = C\,. %\eqno(5.3.28)$$ (4.134)

Наиболее простой для анализа является модель с постоянным угловым моментом $l \equiv r^2\Omega = {\rm const}$. Из (4.131) следует, что $\Psi_\textrm{вр} = l^2 / (2 r^2)$ и, принимая во внимание (4.129), уравнение (4.133) сводим к виду

$${{\cal P} \o \rho} = {G\,M_1 \o (n + 1)\,r_0}\left\{{r_0 \... ... + z^2}} - {1 \o 2}\left(r_0 \o r\right)^2 - C\,'\right\}\,,$$ (4.135)

где $r_0 = l^2 / GM_1$, $C\,' = - Cr_0 / GM_1$. На поверхности тора давление обращается в нуль и уравнение (4.135) дает радиальные координаты $r_{\pm}$, которые ограничивают тор в плоскости $z=0$:

$$r_{\pm} = r_0/\Bigl(1 \mp \sqrt{1 - 2\,C\,'}\Bigr). %\eqno(5.3.30)$$ (4.136)

На рис. 4.21 показаны равновесные конфигурации для различных $C\,'\!\!$.

Figure: Поверхность нулевой плотности тора в плоскости $(r,z)$ для различных $C\,'$: 1 -- $C\,'=~0,\!495$; 2 -- $C\,' = 0,\!444$; 3 -- $C\,' = 0,\!25$

=0.6

Диск является толстым и это вынуждает нас использовать трехмерные уравнения газодинамики в цилиндрической системе координат. Проделаем стандартную процедуру линеаризации этих уравнений -- проводя выкладки в духе главы 4, получим для возмущенных величин $\propto \exp\{im\varphi~-~i\omega t\}$ (выделенных значком `` $\tilde {}$ ''):

$$\tilde u_r = {i \o \hat\omega^2 - \varkappa^2}\left(\hat\... ...t r} - {\varkappa^2\hat\omega\,m\,W \o 2\,r\,\Omega}\right),$$ (4.137)

$$\tilde u_\varphi = {1 \o \hat\omega^2 - \varkappa^2}\left... ...\o \hat\omega}{d\Omega \o dr}\right)\right\}, %\eqno(5.3.32)$$ (4.138)

$$\tilde u_z = i{\Vert W \o \Vert z}, %\eqno(5.3.33)$$ (4.139)

$$i\,\hat\omega\,\tilde\rho = {\rm div}(\rho_0 \tilde{\vec u}),$$ (4.140)

$$\tilde {\cal P}/{\cal P}_0 = \gamma\,\tilde\rho/\rho_0,$$ (4.141)

где $W = -\tilde {\cal P}/(\rho_0\hat\omega)$; $\hat\omega = \omega - m\Omega$. Используя (4.137)$\!$-$\!$ (4.139) и (4.141), исключим из (4.140) возмущенные плотность и скорость, в результате получим уравнение для величины $W$ [692]:

$${D^2 \o r}{\Vert \o \Vert r}\biggl\{{\rho_0 r \o D}\biggl(\... ...arkappa^2{\hat\omega\,m\,W \o 2\,r\,\Omega}\biggr)\biggr\}~+$$

$$+~\rho_0 D\biggl\{- {m^2 W \o r^2}\biggl({\varkappa^2 r \o ... ... \hat\omega\,m \o 2\,r\,\Omega}{\Vert W \o \Vert r}\biggr\}~+$$

$$+~D^2 {\Vert \o \Vert z}\biggl(\rho_0{\Vert W \o \Vert z}\... ...t\omega^2 \rho^2_0 W/\gamma\,{\cal P}^2_0 = 0, %\eqno(5.3.36)$$ (4.142)

где $D \equiv \hat\omega^2 - \varkappa^2$.

Уравнение (4.142) описывает динамику малых возмущений в торе с произвольными равновесными зависимостями $\rho_0(r,z)$, $\Omega(r)$. Анализ уравнения (4.142) был проведен в работах [229, 692, 693, 695, 765] в рамках как аналитического подхода, так и численного.

В пределе коротких волн в $z$- и $r$-направлениях можно искать решения в виде

$$W = W_0 \exp\{i\,(k_r r + k_z z)\}. %\eqno(5.3.37)$$ (4.143)

В этом случае уравнение (4.142) дает дисперсионное уравнение

$$\hat\omega^2 = \varkappa^2 k^2_z/(k^2_z + k^2_r), %\eqno(5.3.38)$$ (4.144)

и для устойчивости достаточно $\varkappa^2 = 4\Omega^2\,[1 + rd\Omega / (2\Omega\,dr)] \ge 0$.

Естественно, основной интерес представляет изучение возможности нарастания возмущений, то есть получение условий, при которых $Im(\omega) > 0$. Прежде всего обсудим роль сжимаемости газа. Ограничимся рамками модели с постоянным удельным угловым моментом $l = r^2\Omega =$ const. Поскольку в этом случае , то (4.142) принимает следующий вид:

$${\Vert \o r\Vert r}\left(\rho_0 r {\Vert W \o \Vert r}\rig... ...right) = -{\hat\omega^2 \rho_0 \o c^2_s} W\,. %\eqno(5.3.39)$$ (4.145)

Определим величину $\omega_c \equiv \omega - m\Omega_c$, где $\Omega_c$ -- некоторое фиксированное значение угловой скорости. Умножим уравнение (4.145) на $rW^*$ ($W^*$ -- комплексно сопряженная $W$) и проинтегрируем по плоскости ($r,z$). В результате получим

$$A\,\omega^2_c - B\,\omega_c - C = 0\,, %\eqno(5.3.40)$$ (4.146)

где

$$A = \int\vert W\vert^2 \rho_0 c^{-2}_s r\,dr\,dz\,, %\eqno(5.3.41)$$ (4.147)

$$B = \int\vert W\vert^2 2m\,(\Omega - \Omega_c)\,\rho_0 c^{-2}_s r\,dr \,dz\,, %\eqno(5.3.42)$$ (4.148)

$$C = \int\rho_0\left[\biggl\vert {\Vert W \o \Vert r}\biggr... ...a_c)^2 \o c^2_s}\right)\right]\,r\, dr\,dz\,. %\eqno(5.3.43)$$ (4.149)

Запишем решение уравнения (4.146)

$$\omega_c = \left(B \pm \sqrt{ B^2 + 4\,A\,C }\right)\big/(2\,A)\,.$$ (4.150)

Рост возмущений [Im$(\omega) > 0$] возможен только в случае больших отрицательных значений величины $C$. Из (4.149) следует, что всегда существует такое критическое значение скорости звука, при котором $B^2 + 4\,A\,C = 0$. В несжимаемом пределе ( $c_s \rightarrow \infty$) возмущения заведомо являются устойчивыми. Сжимаемый газ может быть неустойчивым, причем чем меньше величина $c_s$, тем более благоприятны условия для развития неустойчивости. Второй вывод заключается в том, что осесимметричные возмущения ($m=0$) являются устойчивыми.

Для объяснения природы неустойчивых глобальных мод отметим, что в газовом диске (или торе) имеются три характерных радиуса, определяемых соотношением

$$\omega - m\,\Omega(r) = l\,\varkappa(r)\,, %\eqno(5.3.45)$$ (4.151)

где $l=0$ соответствует радиусу коротации; $l=-1$ -- внутреннему, а $l=1$ -- внешнему линдбладовскому резонансу. Между линдбладовскими резонансами лежит запрещенная для нейтральных волн область, где решения (4.145) с Im$(\omega)=0$ не могут осциллировать по радиальной координате, а снаружи, наоборот, расположены так называемые колебательные полости.

В основе неустойчивости лежат два механизма:

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

2). Излучение углового момента и энергии из коротационного резонанса. Физика данного эффекта основана на резонансном обмене между волной и течением и сходна с затуханием (усилением) Ландау. Отметим, что для звездного диска обмен энергией на радиусе коротации подробно разобран в монографии Поляченко и Фридмана [163], однако он существенно отличается в силу бесстолкновительности плазмы звезд и pазличия в знаках плотности энеpгии акустических волн в газе и волн, pаскачивающихся на гpавитационной ветви колебаний звездного диска.

В работе [694] было показано, что для политропного газа из системы (4.137)-(4.141) может быть получено уравнение баланса углового момента с источником:

$${\Vert h \o \Vert t} + {\Vert (rf) \o r\,\Vert r} = -\,{{\... ... \o dr}\left( {1 \o \s_0 r}{dH \o dr}\right), %\eqno(5.3.46)$$ (4.152)

где -- плотность углового момента и $\ f=r\,\s_0\, \cdot\langle\tilde u_r\tilde u_\varphi \rangle_\varphi$ -- плотность радиального потока углового момента возмущений, осредненные по азимутальному углу, $H = \Omega (r)\,r^2$ -- распределение углового момента вещества в диске. Если правая часть (4.152) положительна, она описывает излучение, а если отрицательна -- поглощение углового момента возмущений в единицу времени в окрестности коротационного радиуса (где $\vert \omega - m\Omega \vert$ минимален).

Внутри радиуса коротации волновой узор вращается медленнее газа, а снаружи -- быстрее; поэтому угловой момент волны отрицателен во внутренней области и положителен во внешней. Этим обусловлена возможность развития неустойчивости как в случае излучения, так и в случае поглощения углового момента возмущений на коротации, и различная локализация неустойчивых возмущений в диске. В первом случае угловой момент излучается на коротации и передается наружной ``положительной'' моде, а во втором -угловой момент отбирается у моды внутри радиуса коротации, которая и без того обладает отрицательным угловым моментом, из-за чего последний растет по абсолютной величине ^4.18. Из приведенного описания механизмов неустойчивости понятно, что скорость роста возмущений определяется скоростью переноса энергии и углового момента звуковой волной, то есть рост возмущений происходит в динамической шкале времени:

$${\rm Im}(\omega) \sim c_s / L\,, %\eqno(5.3.47)$$ (4.153)

где $L$ -- характерный радиальный масштаб в диске.

В заключение отметим, что резонансные неустойчивости газовых дисков и торов являются многомодовыми -- для каждой из азимутальных мод (с фиксированным номером $m$) удается обнаружить до 30 неустойчивых гармоник, различающихся числом и расположением узлов функции возмущенного давления по радиальной координате и скоростью вращения волнового узора (радиусом коротации). Указанная многомодовость, с одной стороны, и относительная малость инкремента -- с другой (см. (4.153)), приводят к возможности одновременного существования в диске большого числа таких гармоник ^4.19. Поэтому, несмотря на то, что каждая спиральная волна имеет малую амплитуду, суммарный перенос этими волнами углового момента на периферию диска может оказаться значительным.

4.3.4 Неустойчивости, обусловленные радиальной неоднородностью
давления и энтропии

В предыдущих двух пунктах были рассмотрены неустойчивости, вызванные неоднородностью скорости движения вещества. В п.  изучена градиентно-энтропийная неустойчивость (ГЭН) в плоскости неоднородного тонкого газового диска. Прямыми аналогами ГЭН являются конвективная и Рэлея-Тейлора неустойчивости. Они могут развиваться в неоднородной среде, находящейся в поле тяжести ^4.20. В связи с этим необходимо сказать о встречающемся в литературе утверждении о том, что, поскольку вещество, вращающееся по кеплеровским орбитам, невесомо, то неустойчивость конвективного типа в случае дисковой аккреции не играет роли. Разумеется, это справедливо для случая постоянного давления, но такая модель представляется достаточно искусственной. Рассмотрим стационарный равновесный неоднородный диск без радиального движения и вязкости, для которого баланс сил в радиальном направлении определяется (). Хотя два последних слагаемых дают малый вклад в равновесие, но их учет является необходимым для развития неустойчивостей в плоскости диска. Удобно ввести удельную равновесную силу, для которой с учетом уравнения () запишем

$$g \equiv \frac{d\Phi}{dr} - \frac{v_0^2}{r} + D\frac{p_0}... ...a_0}\frac{d\ln\Omega_z}{dr}= - \frac{dp_0}{\sigma_0\,dr} \,,$$ (4.154)

она мала в случае тонкого диска ( $g\ll d\Psi/dr \simeq r\Omega^2$), но этого оказывается достаточно для развития конвективных неустойчивостей с характерным временем роста возмущений $\tau = 1/{\rm Im}(\omega) \sim (h^2 k/r)\Omega$ (эту оценку нетрудно получить из () или (4.157)). Для равновесных параметров $f=\{ \sigma_0, p_0, v_0=r\Omega, \ldots \}$ определим масштабы радиальной неоднородности $L_f=(d\,\ln(f)/d\,r)^{-1}$.

Возможность раскачки ГЭН определяется, как мы видели в п. , соотношением характерных масштабов неоднородности поверхностной плотности и температуры ($L_\s$ и $L_T$). К настоящему времени построено уже довольно много осесимметричных стационарных моделей АД, отличающихся в конечном счете различным радиальным распределением поверхностной плотности и температуры и, следовательно, значениями $L_\s$ и $L_T$ [ $L_f = \{d(\ln f)/dr\}^{-1}$] [144].

Представим все функции в виде $f(r,\varphi,t)=f_0(r)+\t{f}(r,\f,t)$ и линеаризуем уравнения относительно возмущенных величин $\t{f}$ с учетом (4.154). В рамках ВКБ-приближения считаем для возмущений $\t{f}=f_1\cdot \exp\{ -i\,\omega t+i\,k_r r + i\,m\f \}$. Условие существования нетривиальных решений для линеаризованной системы приводит к дисперсионному уравнению 4-й степени относительно частоты $\omega$. Если рассматривать отдельно радиационно-доминирующий диск $p_{rad} \gg p_{gas}$ и обратный случай $p_{rad}\ll p_{gas}$, то в результате дисперсионное уравнение принимает вид:

$$\hat{\omega}^4 - \hat{\omega}^2\cdot \left[ \varkappa^2 +... ...\right] - 2\,s\,\Omega\hat{\omega}kc_T^2 \,\Gamma \, \times$$ (4.155)

$$\times {d \o dr} \ln \left( {2\Omega\,p_0^{2/\Gamma}\o \v... ... dr} \ln {p_0\o \sigma_0^\Gamma\Omega_z^{\Gamma-1}} = 0 \,,$$

где $s=\sin(\Theta)=k_\varphi/k=k_\f/\sqrt{k_r^2+k_\f^2}$ определяет степень неосесимметричности возмущений, $k_\f=m/r$. В случае $p_{rad} \gg p_{gas}$ необходимо считать $\gamma = 4/3$ и $\Gamma=9/7$ (см. ()). Строго говоря, (4.155) справедливо для $k_r r\gg 1$, $k_r L_f \gg 1$ (неоднородность учитывалась в первом порядке), $\vert\sin(\Theta)\vert\ll 1$.

В пределе изэнтропической модели $\Oo {p_0\o \sigma_0^\Gamma\Omega_z^{\Gamma-1}} ={\rm const}$ (см. подробнее § ) порядок уравнения (4.155) понижается, поскольку энтропийная ветвь колебаний вырождается $\hat{\omega}=0$. В общем случае уравнение (4.155) описывает две высокочастотные ветви колебаний (акустические), для которых приближенно можно записать $\hat{\omega}^2\sim\varkappa^2+k^2c_s^2$, и две низкочастотные (энтропийную и вихревую).

Формальный переход в (4.155) к невращающейся среде $\Omega=0$ и $\Omega_z={\rm const}$ дает для низкочастотных волн $\Oo \omega^2 = {\sin^2(\Theta)\o \gamma^2}{1\o L_p}\left({\gamma\o L_\sigma} - {1\o L_p}\right)\cdot c_s^2$, и для устойчивости необходимо:

$${1\o L_p}\left({\Gamma\o L_\sigma} - {1\o L_p}\right) > 0\,,$$ (4.156)

что в точности совпадает с условием конвективной устойчивости в неоднородной среде [92]. Отметим, что в осесимметричной модели с учетом вертикальных движений устойчивость определяется критерием Hoiland для вращающегося потока газа [506] и рассматривалась также в [183]. Возможность развития конвективной неустойчивости во вращающемся газовом диске в приложении к проблеме спиральной структуры обсуждалась в [155].

Границы конвективной устойчивости. Уравнение (4.155) позволяет определить границы конвективной неустойчивости с учетом дифференциальности вращения ( $1/L_\Omega\ne 0$) и конечной толщины диска в главном приближении ( $1/L_z\equiv d\,\ln\Omega_z/dr \ne 0$). Перепишем (4.155) в пределе низкочастотных волн ( $\vert\hat{\omega}^2\vert\ll \Omega^2$, малость $\vert\omega\vert/\Omega$ не требуется):

$$\Lambda\cdot \nu^2 + \delta\cdot\sin(\Theta)\cdot\left[ {2... ...{r\o L_\Omega}- 2\,{r\o L_{\varkappa}} \right]\cdot \nu\, +$$

$$+ \sin^2(\Theta)\, {\delta^2\o\Gamma^2}\cdot \left[ {... ...\right] \cdot\left[ {r\o L_p} + {r\o L_z} \right] = 0 \,,$$ (4.157)

где $\nu=\hat{\omega}/\Omega$, , $\delta=2kc_s^2/(\Omega^2\,r)$. Пpинимая во внимание оценку для полутолщины диска $h\sim c_s/\Omega$, параметр $\Lambda$ может принимать значения от $\Lambda \sim 1$ в длинноволновом пределе до $\Lambda\lee 100$ для возмущений, у которых длина волны сравнима с толщиной диска. Для многих стационарных моделей АД можно принять степенную зависимость равновесных параметров диска от радиальной координаты: $r/L_f={\rm const}$. При выполнении условия ${\rm Im}(\nu)> 0$ имеем неустойчивые решения, рассмотренные в [144] в случае $r/L_z=0$.

Figure: Границы конвективной устойчивости для моделей: линия 1 -- $r/L_\Omega =r/L_z=-3/2$, $\gamma =5/3$, $\Lambda =10$. В областях A,$\!$C возмущения неустойчивы, B -- зона устойчивости; 2 -- для невращающейся атмосферы в соответствии с (4.156) (в областях I, III возмущения неустойчивы; II, IV -- зоны устойчивости); 3 -- для случая $r/L_z=0$

=0.44

На рис. 4.22 на плоскости параметров $r/L_\sigma$ и $r/L_p$ изображены границы конвективной неустойчивости, определенные из условия ${\rm Im}(\nu)=0$, для базовой модели $r/L_\Omega =r/L_z=-3/2$, $\gamma =5/3$, $\Lambda =10$. Здесь же для сравнения приведены границы устойчивости для невращающейся среды, для которой равновесие обеспечивается только внешней силой и градиентом давления, а также для модели $r/L_z=0$. Как видим, вращение и конечная толщина диска заметно изменяют условия существования конвективной неустойчивости. Причем, в зависимости от знаков $L_\sigma$ и $L_p$ области на плоскости ( $r/L_\sigma,r/L_p$), благоприятные для возникновения градиентной неустойчивости, могут как увеличиваться, так и уменьшаться. Для равновесных распределений с $r/L_\sigma<0$ и $r/L_p<0$ дифференциальность вращения и конечная толщина дисков являются стабилизирующими факторами. При любых значениях $r/L_\sigma$ с уменьшением $\gamma$ критическое значение $\vert r/L_p\vert$ становится меньше. Этот эффект согласуется с критерием (4.156).

Параметр $\Lambda$ характеризует пространственную структуру возмущений. Большие значения параметра $\Lambda\gg 1$ достигаются для коротковолновых в радиальном направлении волн $k\sim \Omega/c_s$. Величина $\Lambda$ выше для возмущений с большим азимутальным номером $m$. Следует подчеркнуть, что наиболее неустойчивыми с точки зрения уравнений (4.155), (4.157) являются предельно неосесимметричные возмущения $\sin(\Theta) = 1$ (поскольку ${\rm Im}\,\omega\propto \sin\{\Theta\}$), для которых заведомо нарушаются приближения, лежащие в основе дисперсионных уравнений. Границы конвективной неустойчивости в плоскости диска зависят от степени неосесимметричности возмущений только через параметр $\Lambda(\Theta)$. Наиболее сильно параметр $\Lambda$ влияет на модели с $r/L_p>0$. И в случае $r/L_p>0$ неустойчивыми могут быть только мелкомасштабные возмущения.

Наличие предельных переходов свидетельствует, что физический механизм, обусловливающий рост возмущений со временем, аналогичен классической конвективной неустойчивости при наличии градиента энтропии, сонаправленного с внешней силой, либо неустойчивости Рэлея-Тейлора, если градиенты равновесных давления и плотности имеют разные знаки. Для неустойчивости, котоpая в основном имеет $r/L_p<0$, следует pазличать две области. В одной повеpхностная плотность убывает с увеличением pадиальной кооpдинаты ($r/L_\sigma<0$), пpямым аналогом такой неустойчивости является конвективная в одноpодном поле тяжести. В дpугом случае повеpхностная плотность pастет с удалением от центpа ($r/L_\sigma>0$), что соответствует неустойчивости Рэлея-Тейлоpа.

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

Нелинейная стадия радиальной конвекции. Проведенный выше анализ подразумевает выполнение, как минимум, трех условий: 1) линейность; 2) малость длины возмущения в радиальном направлении по сравнению с масштабом неоднородности; 3) волны должны быть сильнозакрученными (ограничение связано с дифференциальностью вращения). Инкремент пропорционален $\sin(\Theta)$, поэтому необходимо изучить влияние сильной дифференциальности вращения на формирующиеся на линейной стадии конвективные ячейки. Обсудим результаты нелинейной эволюции конвективно неустойчивого диска методом TVD-E (см. [754]) без учета вязкости, ограничившись $\Omega_z={\rm const}$, используя безразмерные координаты и время: $t=1$ -- период обращения на радиусе $r=1$ [541].

В начальный момент времени задаем степенные зависимости плотности и давления $(r/L_p={\rm const}, r/L_\sigma={\rm const})$. Если начальное равновесное состояние, определяемое функциями $p_0(r)$ и $\sigma_0(r)\!$, обеспечивает устойчивость диска в соответствии с (4.155) ( ${\rm Im}(\omega)=0$), то роста возмущений со временем в численных моделях не наблюдается.

Figure: Изолинии поверхностной плотности ${\sigma }/\sigma _0$ в два разных момента времени $t_1=15000$ (a), $t_2=41000$ (б)

=0.999

Рассмотрим модель с $r/L_p=-3/2$ и $r/L_\sigma=-1/2$, $\gamma =5/3$, которая попадает в неустойчивую область в соответствии с (4.155). Независимо от амплитуды начального возмущения происходит формирование растущих со временем волн, имеющих характерную спирально-ячеистую структуру. На рис. 4.23 показаны изолинии отношения плотности ${\sigma}(r,\f)$ к равновесному значению $\sigma_0(r)$ в два момента времени $t_1=15\,000$ и $t_2=41\,000$. На начальной стадии формируются типичные конвективные ячейки (см. рис. 4.23$\!$a) с небольшой относительной амплитудой поверхностной плотности $\vert{\sigma}-\sigma_0\vert/\sigma_0\lee 5~\%$. Со временем происходит рост амплитуды возмущений и усложнение пространственной структуры из-за дифференциальности вращения диска (см. рис. 4.23$\!$b).

Закон вращения слабо отличается от кеплеровского $\Omega_K$. Поле возмущений скорости наглядно демонстрирует вихревой характер течения. Параметры $k$ и $m$, характеризующие пространственную структуру возмущений, независимы в рамках линейного анализа. И формально, при прочих равных, наиболее неустойчивы волны с углом закрутки $\Theta\rightarrow 90^\circ$. В численной модели нарастают волны с $\Theta\simeq 15^\circ\div 25^\circ$ в результате компромисса между неустойчивостью и дифференциальностью вращения. Структура неустойчивых возмущений на нелинейной стадии определяется параметрами модели (${\cal M}$, $m$, $\delta_{p}$, $\delta_\sigma$, $\gamma$), причем азимутальное число $m$ мы можем изменять, варьируя начальные возмущения вдоль угла $\varphi$. С ростом азимутального числа возмущения становятся более мелкомасштабными и в радиальном направлении. Аналогичный эффект имеем при увеличении числа Маха.

Из линейного анализа следует, что инкремент неустойчивости пропорционален ${\rm Im}\,\omega\propto m$, что в целом подтверждается на начальной стадии эволюции возмущений. На рис. 4.24 показаны зависимости от времени максимальных значений относительного возмущения для выбранной конвективной ячейки для различных $m$ и ${\cal M}$. За характерные времена $t^{(sat)}\simeq 350{\cal M}$ происходит нарастание возмущений до существенно нелинейной стадии, близкое к насыщению. Амплитуда возмущений с малым азимутальным числом нарастает медленнее.

Figure: Зависимость максимального значения возмущенной плотности $\t{\sigma}/\sigma_0$ для конвективной ячейки от времени для различных моделей: ${\cal M}=10$, 1 -- $m=2$, 2 -- $m=4$, 3 -- $m=6$, 4 -- $m=8$, 5 -- $m=12$. Чем больше номер гармоники, тем выше темп роста возмущений, который следует экспоненциальному закону вплоть до насыщения

=0.51

Сверхзвуковое натекание газа на область повышенной плотности в зоне спиральной волны приводит к формированию слабых ударных волн на задней кромке. Радиальные профили течения аналогичны полученным в п. 4.3.6 при рассмотрении диссипативной неустойчивости, поскольку волны являются сильно закрученными (рис. 4.25).

Figure: Радиальные зависимости $\sigma /\sigma _0$, $p/p_0$, $c_s/c_{s0}$, $u$ (левая ось ординат), $v/v_0$ (правая ось) для азимутального угла $\varphi =0$ в момент времени $t=2\cdot 10^4$. Указаны положения ударных волн

=0.5

Конвекция в плоскости диска приводит в среднем к падению вещества на гравитационный центр. Поток вещества $\Oo\dot{M}=r\,\int_0^{2\pi} \sigma u\,d\varphi$ в целом отрицателен. Поток удельного углового момента направлен наружу и среднее радиальное движение к центру связано с отводом углового момента спиральными волнами, что обеспечивает аккрецию в результате перераспределения вещества.

Полученные результаты, сделанные в рамках квазикеплеровского диска, в основном сохраняются в случае законов вращения с $r/L_\Omega > -3/2$. В приложении к газовым подсистемам дисковых галактик значения параметра $r/L_\Omega$ лежат от 0 (близкая к твердотельной центральная зона) до $-1$ -- кривая типа ``плато'' характерна для большинства галактик в существенной части диска.

4.3.5 Приливная неустойчивость

В газовом диске, вращающемся вокруг компактного объекта массой $M_1$, может развиваться неустойчивость, связанная с приливным влиянием со стороны второго компонента массой $M_2$ тесной двойной системы. Физику этой неустойчивости легко понять, рассматривая движение пробной частицы в гравитационном поле двойной системы.

Пусть в случае $M_2 = 0$ пробная частица движется по периодической орбите. Как известно, малые колебания вблизи траектории происходят с эпициклической частотой $\varkappa = \varkappa_0\!$, и отклонения описываются обычным уравнением гармонического осциллятора (см. п. 1.1.3)

$${d^2\xi \o dt^2} + \varkappa^2 \xi = 0\,. %\eqno(5.3.49)$$ (4.158)

С учетом нормальной звезды ($M_2 \neq 0$), вращающейся с периодом $T_0 = 2\pi /\Omega_0$, частота колебаний в (4.158) будет периодической функцией времени и можно считать

$$\varkappa^2 = \varkappa^2_0\, \Bigg[ 1 + \sum\limits^{\infty}_{n=1}\varepsilon_n\, \cos(n\,\Omega\, t)\Bigg]\,,$$ (4.159)

где $\Omega$ -- средняя угловая скорость при движении частицы по периодической орбите, коэффициенты $\varepsilon_n$ определяются потенциалом второй звезды [484]. В этом случае уравнение (4.158) описывает параметрический резонанс [92]. При определенных соотношениях между собственной частотой и частотой вынуждающей силы отклонение $\xi$ начинает быстро нарастать со временем -- развивается приливная неустойчивость. Поскольку период двойной системы $T_0$ зависит от относительной массы возмущающего тела $\mu = M_2 / (M_1 + M_2$), то и параметры $\varepsilon_n$ определяются величиной $\mu$, причем отличны от нуля только нечетные члены (из-за симметрии потенциала). Наиболее интенсивным является резонанс с $\varepsilon_3 = 1,\!85\mu$, и, следовательно, приливная неустойчивость возникает, когда период пробной частицы составляет $\simeq 1/3$ от периода второй звезды в инерциальной системе ^4.21 [484, 485, 894].

В случае кеплеровского диска период обращения частиц растет с радиусом и существует значение $r_0$, при котором наступает резонанс. Как мы видели в п. 4.1.1, 4.1.4, из-за приливного взаимодействия в ТДС диск имеет конечный размер. Внешний радиус диска $R$ определяется отношением $q = M_2/M_1$ -- чем больше $q$, тем меньше величина $R$. Для значений $q \ll 1$ диск простирается за радиус $r_0$. В результате во внешней области диска развивается приливная неустойчивость, приводящая к возникновению медленно прецессирующего эллиптического диска. Период обращения такого эллиптического диска на 3-6 % превышает орбитальный период двойной. При $q \simeq 1/4$ выполняется равенство $R \simeq r_0$. Поэтому в системах с $q \ge 1/4$ приливная неустойчивость не развивается.

Возникновение приливной неустойчивости, обусловленной параметрическим резонансом $3:1$ между орбитальным движением газа в диске и орбитальным вращением двойной, подробно исследовано в рамках численного моделирования газового диска в ТДС [484, 485, 894]. Возможно, что так называемые ``супергорбы'' в системах типа SU UMa (периодическое увеличение блеска кривой на 20-30 % с периодом, несколько превышающим орбитальный период во время супервспышек (см. п. 1.5.1)) вызваны эллиптическим диском.

4.3.6 Диссипативно-акустическая неустойчивость

Динамика длинноволновых возмущений. Хаpактеpной особенностью стандаpтной модели АД является наличие зависимости туpбулентных диссипативных коэффициентов от паpаметpов диска. В §  pассматpивалась динамика малых возмущений без учета возмущения величины вязкости. Обсудим в pамках пpедельно пpостой модели влияние зависимости динамической вязкости от плотности ( $\eta = \s\nu \propto \s^{1+\delta}$) на хаpактеp звуковых колебаний в плоскости диска.

С учетом возмущения динамической вязкости $\eta_1/\eta_0 = \break (1+\delta)\,\s_1/\s_0$ для кеплеpовского диска система уpавнений для относительно длинноволновых возмущений ()$\!$-$\!$() пpинимает вид:

$$-i\omega\s_1 + i\s_0 ku_1 = 0\,, %(5.3.51)$$ (4.160)

$$-i\omega u_1 - \ \ 2\Omega v_1 = -ikc_T^2p_1/p_0 -\lambda_0k^2u_1\,, %(5.3.52)$$ (4.161)

$$-i\omega v_1 + \ {1\o 2}\,\Omega u_1 = -k^2\nu_0 v_1 - ik\,{3\o 2}\,\Omega\nu_0\eta_1/\eta_0 \,. %(5.3.53)$$ (4.162)

В пеpвом пpиближении будем считать $p_1/p_0 =\gamma_s \s_1/\s_0$ (в pаботах [205, 209] в уpавнении теплового баланса учитывались диссипативные фактоpы). В pезультате получаем диспеpсионное уpавнение

$$\omega^3 + ik^2(\lambda_0+\nu_0)\omega^2 - (\Omega^2+k^2c... ... - ik^2\nu_0\Omega^2[3(1+\delta) - c_s^2k^2/\Omega^2] = 0\,.$$ (4.163)

Будем искать выpажение для частоты длинноволновых звуковых колебаний в виде $\omega/\Omega = 1+\varepsilon k^2$ ( $\vert \varepsilon k^2\vert \ll 1$). Таким обpазом, для мнимой части частоты имеем

$${\rm Im}(\omega) \simeq {k^2 \nu_0\o 2}\, \big( 2+3\delta-\lambda_0/\nu_0 \big) \,. %\eqno(5.3.55)$$ (4.164)

Пpи выполнении условия $(2+3\delta)\nu_0 > \lambda_0$ имеем неустойчивость акустических колебаний. Отметим, что для $\alpha$-модели в области доминиpования газового давления и непpозpачности $k_T$ имеем $\delta=2/3$ (см. § 4.1). Hаpастание возмущений целиком обусловлено зависимостью динамической вязкости $\eta$ от паpаметpов диска. Учет тепловых пpоцессов, излучения, давления излучения, зависимости вязкости от темпеpатуpы, вертикальной структуры диска и вертикальных колебаний изменяет значения инкрементов, сохраняя возможность неустойчивости $Im(\omega) > 0$ [205, 209, 538].

Зависимости вязкости и непрозрачности от поверхностной плотности и толщины диска ( $\nu \propto \sigma^{\delta_\sigma} h^{\delta_h}$, $\bar\kappa \propto \sigma^{\Delta_h} h^{\Delta_h}$) оказывают решающее влияние на диссипативные неустойчивости акустических, а также вязкой и тепловой (см. п. 4.3.1) ветвей колебаний. Построенные к настоящему времени многочисленные модели АД оказываются как в устойчивой, так и в неустойчивой областях по параметрам $\Delta_{h}$, $\Delta_\sigma$, $\delta_h$, $\delta_\sigma$ для различных ветвей колебаний [209].

Неустойчивость коротковолновых акустических волн. Модель тонкого диска накладывает ограничение на длину волны $\lambda \gg h$ и частоты возмущений [36]. Поэтому для корректного исследования динамики возмущений с $\lambda \lee h$ необходимо рассматривать вертикальные движения и $z$-структуру АД. Модель тонкого диска позволяет рассматривать только пинч-колебания ($S$-мода). Для них возмущенное давление является симметричной функцией $z$-координаты, а вертикальные смещения газа не сдвигают центр массы в диске относительно плоскости симметрии $z=0$. Таким образом, изгибные колебания ($AS$-мода) исключаются из рассмотрения. Кроме того, двумерные модели не позволяют изучать высокочастотные (отражательные) гармоники с $\omega\gg \Omega$. Изучим динамику акустических возмущений с учетом вертикальных движений. Кроме того, попутно рассмотрим вопрос о пределах применимости модели тонкого диска.

Задача определения собственных частот. Аналогично п. 4.3.2.1 рассмотрим математическую модель, позволяющую определять собственные частоты для различных неустойчивых мод в $z$-неоднородном вязком осесимметричном диске [538]. В отличие от (4.113), (4.114), пренебрежем магнитным полем, но учтем вязкость.

Будем считать, что равновесная скорость в диске имеет только $r$- и $f$- компоненты: ${\bf\vec u}_0(r,z) = (u_0,\, v_0, \, 0 )$. Для компонент тензора вязких напряжений примем:

$$\Pi_{ij} = - \alpha_{ij}\,{\cal P} \,,$$ (4.165)

где $i,j = (r,\f, z)$. Параметры $\alpha_{rz}$, $\alpha_{\f z}$ и $\alpha_{r \f}$ определяют уровень турбулентности в диске. Причем $\alpha_{rz}$ и $\alpha_{\f z}$ обусловлены сдвиговым характером течения в $z$-направлении, а величина $\alpha_{r \f}$ связана с дифференциальностью вращения в плоскости диска и совпадает с $\alpha$-параметром стандартной теории дисковой аккреции [221]. С учетом оценок $\alpha_{r\f}{\,\bf :\,}\alpha_{rz}{\,\bf :\,}\alpha_{\f z}{\,\bf :\,}\alpha_{r... ...ert{\Vert u_0\o\Vert r}\right\vert {\,\bf :\,} \left\vert{u_0\o r}\right\vert$, $h \ll r$, ограничимся в уравнениях движения только компонентой $\Pi_{r\f}$. Считаем, что равновесный баланс сил определяется системой уравнений:

$${v_0^2 \o r} = {1\o \varrho_0}{\Vert {\cal P}_0 \o \Vert ... ... P}_0 \o \Vert z} = - \varrho_0 \,{\Vert \Psi \o \Vert z} \,.$$ (4.166)

Равновесные функции представим в виде: $f_0(r,z) = f_{01}(r)\cdot f_{02}(z)$. Воспользуемся решениями п. . Соотношение

$$c_s^2(r,z) = {\gamma {\cal P}_0\o \varrho_0} = c_s^2(r,0)\... ...\textrm{где}}\ \ c_s^2(r,0) = {\gamma \o a}\Omega_K^2 h^2 \,$$ (4.167)

определяет адиабатическую скорость звука ($\gamma$ -- показатель адиабаты, в общем случае $\gamma \ne \ell$). Для равновесных скоростей $v_0$ и $u_0$ из уравнений (4.166) с учетом (4.167) получим:

$$v_0(r,z) = r\,\Omega(r,z) = \sqrt{r^2 \Omega_K^2 + \ell_{\cal{P}}\, c_s^2/\gamma} \,,$$ (4.168)

$$u_0(r,z) = -{\alpha (2+\ell_{\cal P})\o \gamma (1+\ell_{v})}\,{c_s^2 \o v_0} \,,$$ (4.169)

где параметры $\ell_{\cal P} = \Vert\, (\ln\, {\cal P}_0) / \Vert\, (\ln\, r)$, $\ell_\varrho = \Vert\, (\ln\, \varrho_0)/ \Vert\, (\ln\, r)$, $\ell_v = \Vert (\ln v_0) / \Vert (\ln r)=1-n$ характеризуют радиальные неоднородности равновесных величин.

После стандартной процедуры линеаризации получим систему уравнений относительно возмущенных величин:

$${\Vert \t \varrho \o \Vert t} + {\Vert \o r \Vert r}\bigg[... ..._0 \t u)\bigg] + {\Vert ( \varrho_0 \t w )\o \Vert z} = 0 \,,$$ (4.170)

$${\Vert \t u\o \Vert t} + u_0{\Vert \t u\o \Vert r} + \t u{... ...rho \o \varrho_0}\,{\Vert {\cal P}_0 \o \Vert r} \right) \,,$$ (4.171)

$${\Vert \t v \o \Vert t} + u_0{\Vert (r \t v) \o r \Vert r}... ...arrho_0}{\Vert (r^2 {\cal P}_0 ) \o r^2 \Vert r} \right) \, ,$$ (4.172)

здесь $g \equiv \Vert \Psi /\Vert z$.

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

$${\Vert \t s \o \Vert t} + u_0{\Vert \t s \o \Vert r} + \t u{\Vert s_0 \o \Vert r} + \t w{\Vert s_0 \o \Vert z} = 0 \,,$$ (4.174)

здесь $s_0$- удельная энтропия равновесного газа. Исключим возмущение энтропии $\t{s}$ из (4.174) с помощью уравнения состояния $s = s({\cal P},{\varrho})$, которое в линейном приближении запишем:

$$\t s = \left({\Vert s\o \Vert {\cal P} }\right)_{\!\!\!\Oo... ...P}\o {\cal P}_0 } - c_{\cal P}{\t \varrho \o \varrho_0 } \,,$$ (4.175)

где $c_V = T\,(\Vert s / \Vert T)_\varrho$, $c_{\cal P} = T\,(\Vert s / \Vert T)_{\cal P}$ -- удельные теплоемкости при постоянных плотности и давлении. Введем функцию смещения вещества в вертикальном направлении $\xi$:

$$w = {d\xi\o dt}= {\Vert\xi\o\Vert t}+u\,{\Vert\xi\o\Vert r}+v\,{\Vert\xi\o r\Vert\f} \,.$$ (4.176)

С учетом ВКБ-приближения в радиальном направлении ($kr \gg 1$, $k$ -- радиальное волновое число) представим решения в виде:

$$\t{f}(r,z,t) = \hat{f}(z)\cdot \exp\{- i\, \omega t + i\, k r \} \,.$$ (4.177)

С учетом (4.175) и (4.177), система (4.170)-(4.174) сводится к двум обыкновенным дифференциальным уравнениям относительно амплитуд возмущенных давления $\hat{{\cal P}}(z)$ и вертикального смещения $\hat \xi(z)$ от плоскости $z=0$:

$${d \hat \xi \o d z} = {D_1(\hat\omega,k,z) \o \hat \omega ... ... \varrho_0} + {g\o c_s^2}\cdot\hat\xi \, , \ \quad\quad\quad$$ (4.178)

$${d \hat{\cal P} \o d z} = \varrho_0 \left\{ \hat\omega^2-{... ...z \right\}\cdot \hat\xi - {g\o c_s^2}\cdot \hat{\cal P} \, ,$$ (4.179)

здесь , $\varkappa = \Omega \sqrt{2(2-n)}$ -- эпициклическая частота, $\Oo D_1 = - (\hat\omega^3 - \hat\omega(\varkappa^2+k^2 c_s^2) - 2i\alpha\Omega\,k^2 c_s^2 )$, $\Oo S_z = {1\o c_V}{\Vert s_0 \o \Vert z} = (\ln {\cal P}_0)' - \ell\, (\ln \varrho_0)'$, штрих означает дифференцирование по $z$-координате ( ), $\hat \xi$ -- комплексная амплитуда материального $z$-смещения от равновесного состояния, так что с учетом (4.176) выполняется $\t w = d\t \xi / dt = -i\hat \omega \,\t \xi$. В более общем виде с учетом радиальных неоднородностей система уравнений на $\hat{\xi}$, $\hat{\cal P}$ приведена в [538]. Отметим, что ВКБ-приближение вдоль $r$ требует, чтобы в бездиссипативном пределе коэффициенты в уравнениях типа (4.178), (4.179) были вещественными. В невязком приближении без магнитного поля имеем согласие с (4.113), (4.114), а без учета вращения имеем предельный переход к уравнениям, рассмотренным в [120].

Численно решая краевую задачу для системы уравнений (4.178), (4.179) с учетом граничных условий (4.117) и (4.115) для $S$-моды (или (4.116) для $AS$-моды), мы находим собственные значения комплексной частоты $\omega$ для заданных распределений равновесных параметров диска вдоль $z$-координаты. Одновременно получаем собственные функции $\hat{\xi}(z)$ и $\hat{\cal P}(z)$, которые с учетом (4.170)-(4.173) дают вертикальные зависимости для всех остальных возмущенных величин. Вышеизложенную схему будем называть ``3D-модель'', а модель тонкого диска -- ``2D-модель''. Будем пользоваться безразмерной частотой $W = \omega / \Omega_K$ и безразмерным волновым числом $K=kh$. Если специально не оговаривается, то параметры принимают значения: $\alpha = 0,\!2$; $\gamma=\ell = 5/3$; $h/r= 0,\!05$; $n = 3/2$.

Фундаментальные $S$- и $AS$-моды в 2D- и 3D-моделях. Предварительно, для коротковолновых возмущений в вертикальном направлении формально подставим $\hat{\cal P}, \hat{\xi}\propto\exp[\,i\,k_z z\,]$ в (4.178), (4.179) и, пренебрегая неоднородностью равновесных величин, запишем дисперсионное уравнение:

$$\omega^4 - \omega^2\, [\varkappa^2 + c_s^2\, (k^2 + k_z^2)... ...Omega c_s^2 k^2 \,\omega + \varkappa^2 c_s^2\, k_z^2 = 0 \,,$$ (4.180)

здесь $k_z$ есть волновое число в $z$-направлении. Естественно, что уравнение (4.180) при $k_z=0$ совпадает с результатом [209] для 2D-модели в соответствующем приближении в случае замены $\gamma$ на плоский показатель адиабаты $\Gamma$ (4.106), (), и квадрат адиабатической скорости звука в плоскости диска должен определяться следующим выражением:

$$c_s^2 = \Gamma \, {\int_0^h {\cal P}(z)\,dz \o \int_0^h \varrho(z)\,dz} \,.$$ (4.181)

Очевидно, что эти колебания соответствуют $S$-моде (4.115).

При $\alpha > 0$ уравнение (4.180) определяет две неустойчивые акустические ветви колебаний и две затухающие моды (вязкую и тепловую). Затухание последних связано с отсутствием радиационного давления, диссипации и охлаждения в тепловом уравнении. При учете этих факторов вязкая и тепловая низкочастотные колебательные ветви могут быть неустойчивыми (см. п. 4.3.1).

Figure: Зависимости величины $\Vert W/\Vert K$ от $K$ в случае: 1 -- 2D-модель с $\Gamma$, вычисляемой по (); 2 -- 2D-модель с $\Gamma =\Gamma _2$; 3 -- 3D-модель

=0.27

На рис. 4.26 показаны зависимости $\Oo{\Vert W\o\Vert K}$ от $K$ в 2D- и 3D-модели, демонстрирующие слабое различие между собой. Там же приведен результат вычисления с плоским показателем адиабаты $\Oo\Gamma_2={3\gamma-1-\gamma\,\delta\o 1+\gamma - \delta}$ ( $\delta=\varkappa^2/\Omega_K^2=1$), полученным при более аккуратном учете вертикальных движений [78]. Линейная зависимость $\Oo{\Vert W\o\Vert K}$ указывает на то, что $\omega\propto k^2$ в широких пределах в соответствии с асимптотикой (4.164).

Таким образом, прямое сравнение дает удовлетворительное согласие между моделью тонкого диска и результатами решения 3D-задачи. Собственные частоты колебаний, полученные из уравнения (4.180) при $k_z=0$ с учетом (4.181) и посредством решения краевой задачи (4.178)-(4.117), слабо различаются в области $kh \lee 3$ (рис. 4.26). Заметные отличия появляются только при $kh \gee 3$. Отметим, что в этой области заведомо нарушается формальное условие применимости 2D-модели (). В случае малых значений параметра $\gamma$ расхождения возникают при больших $kh$, поскольку характерный масштаб неоднородности равновесных параметров в вертикальном направлении увеличивается с уменьшением $\gamma$ [538]. В коротковолновой области инкремент и фазовая скорость возмущений в рамках 3D-модели меньше, чем в 2D-модели, что связано с неоднородным $z$-распределением равновесных величин и поперечной силы тяжести. Инкремент неустойчивости линейно пропорционален $\alpha$-параметру, а значение волнового числа $k$, при котором появляются различия между моделями, слабо зависит от значения $\alpha$-параметра.

Обсуждаемая здесь низкочастотная мода ( $\omega\simeq\sqrt{\Omega^2+k^2c_s^2}$) слабо зависит от $z$-структуры, она имеет место и в 2D-, и в 3D-моделях, поэтому ее естественно называть фундаментальной модой. В случае $\alpha \ll 1$ для фундаментальной моды из (4.180) получаем

$$\omega^2 = \varkappa^2 + k^2 c_s^2 + i \alpha {\Omega c_s^2 k^2 \o \varkappa^2 + c_s^2 k^2}\,,$$ (4.182)

что согласуется с результатами расчета по 3D-модели.

Таким образом, прямая проверка показывает, что 2D-модель пригодна для изучения диссипативно-акустической неустойчивости с характерными пространственными масштабами возмущений в плоскости диска в пределах $\sim(2h \div r)$.

=1.2

Дисперсионные свойства высокочастотных неустойчивых мод. Газовый диск является своеобразным волноводом (см. п. 4.3.2), где может распространяться дискретный набор волн, различающихся вертикальным волновым числом, так что между $z=-h$ и укладывается целое или полуцелое число длин волн. Вертикальная неоднородность равновесных параметров приводит к более сложной картине, поскольку гармонические волны не являются собственными решениями краевой задачи (4.178), (4.179).

Figure: Зависимость $W$ от $kh$ для $S$- и $AS$-мод при различных значениях $\alpha$. Сплошная линия соответствует $\alpha = 0,\!2$, пунктирная линия -- $\alpha = 0,\!1$. Числа на рисунке указывают номера гармоник $j$. На рисунках б) и д) сплошная и пунктирная линии совпадают

=0.3

Кроме фундаментальных $S$- и $AS$-мод в вязких аккреционных дисках может генерироваться любое число неустойчивых гармоник, обусловленных диссипативным механизмом [538]. На рис. 4.27 изображены зависимости собственной частоты $\omega$ от радиального волнового числа $k$ для фундаментальной моды $j=0$ и первых четырех отражательных гармоник $j=1,2,3,4$. Неустойчивыми являются как пинч-колебания, так и изгибные моды. Инкременты рассматриваемых неустойчивостей возрастают с увеличением $k$ и достигают максимума при некоторых значениях волнового числа. Максимум ${\rm Im}(\omega)$ сдвигается в более коротковолновую область с ростом номера гармоники $j$, причем значение максимума увеличивается с уменьшением характерного масштаба возмущений в $z$-направлении. Инкремент у всех гармоник увеличивается с ростом $\alpha$-параметра, а ${\rm Re}(\omega)$ практически не зависит от $\alpha$.

Для каждой гаpмоники $j$ pеальная часть частоты изгибных колебаний в области длинных волн больше соответствующей пинч-моды, что следует уже из пpиближенного диспеpсионного уpавнения $\omega^2 \simeq \varkappa^2 + (k^2 + k_z^2)\, c_s^2$ с учетом $k_z h \simeq \pi j$ для $S$-моды и $k_z h \simeq \pi (j+1/2)$ для $AS$-колебаний. Пинч-возмущения в целом более неустойчивы по сpавнению с изгибными. При больших значениях $kh$ pазличие между $S$- и $AS$-колебаниями исчезает. Причем с pостом номеpа гармоники $j$ это пpоисходит пpи больших значениях pадиального волнового числа $k$. Однако необходимо отметить, что очень коротковолновые возмущения ( $kh \ge 2\div 10$) должны быть устойчивыми из-за наличия радиального градиента возмущенной скорости в вязком тензоре напряжений, который не учитывается в законе (4.165). С учетом этого фактора максимальные значения инкрементов будут уменьшаться с ростом номера гармоники $j$, и очень мелкомасштабные волны (как по $r$-, так и по $z$-координате) не должны проявляться на нелинейной стадии.

В данном параграфе продемонстрирована принципиальная возможность неустойчивости высокочастотных акустических волн в дифференциально вращающемся газовом диске для бегущих возмущений. Они пребывают в системе ограниченное время, и наличие положительного инкремента не означает еще, что возмущения достигнут нелинейной стадии развития. Возмущения распространяются в диске со скоростью $\Oo{\Vert \omega / \Vert k} \sim c_s kh / \sqrt{1+k^2h^2} < c_s$, и характерное время пребывания усиливающихся волн в диске равно $\sim r/h \sim 10^2$ периодам вращения диска. Если принять для оценки ${\rm Im}(\omega) \sim 0,\!1 \Omega$, то для $kh>1$ возможен заметный рост амплитуды волны ( $\exp \{0,\!1\cdot r/h \}$). Неосесимметричные колебания могут играть более важную роль, потому что они находятся в системе дольше, чем осесимметричные волны, и, следовательно, нелинейная стадия развития может заведомо достигаться.

Нелинейное моделирование вязких АД. В работе [212] изучена нелинейная динамика рассмотренной выше диссипативно-акустической неустойчивости с использованием различных численных схем решения уравнений гидродинамики.

Figure: Радиальные зависимости относительных возмущений поверхностной плотности в различные моменты времени для базовой модели. В момент времени $t=2100$ (а) относительная амплитуда поверхностной плотности не превышает 2 % и ударные волны отсутствуют, но амплитуда возмущений увеличивается с удалением от генератора ($r_g=208$). Начиная с момента времени $t\gee 3500$ возникают условия для формирования ударных волн в наиболее удаленных от источника возмущенных областях $r\gee 270$ и $r\lee 130$ (б). К моменту времени $t\simeq 6000$ возмущения достигают внутренней ($r_{in}$) и внешней ($r_{out}$) границ (в)

=0.35

Оказалось, что при наличии начальных возмущений с амплитудой $\gee 0,\!1$ % от равновесных значений акустические колебания могут нарастать до слабых ударных волн за время прохождения звуковой волны по диску (рис. 4.28). В результате в диске формируется нестационарная система мелкомасштабных ударных волн, относительная амплитуда которых может превышать 50 % в случае $\alpha\gee 0,\!2$ (рис. 4.29). Следствием указанных волновых движений в АД является нестационарная компонента светимости диска с амплитудой в несколько процентов (рис. 4.30).

Figure: Радиальные распределения поверхностной плотности $\sigma$ в разные моменты времени аккрецирующего кольца в случае $\alpha=0,\!3$, ${\cal M}=200$. Пунктирная линия показывает начальный профиль

=0.44

Эффективность рассмотренных процессов (формирование ударных волн и квазипериодические осцилляции светимости) сильно зависит от уровня турбулентной вязкости, и образование системы волн значительной амплитуды без учета эффектов отражения от границ в данных моделях оказывается возможным только при $\alpha\gee 0,\!05$.

Figure: Долговременная зависимость светимости неустойчивого кольца, которое расположено в области $70 \lee r\lee 350$ в модели с ${\cal {M}}=200$

=0.999

<< 3. Численное моделиpование звездных 5. Замечания о пpиpоде >>