АСТРОНОМИЧЕСКИЙ ЖУРНАЛ, 2019, том 96, № 1, с. 3-16

УДК 524.3-59

О ВОЗМОЖНЫХ ЭЛЕКТРОМАГНИТНЫХ ПРОЯВЛЕНИЯХ

СЛИВАЮЩИХСЯ ЧЕРНЫХ ДЫР

Институт астрономии РАН, Москва, Россия

Поступила в редакцию 01.08.2018 г.; принята в печать 25.08.2018 г.

Рассматривается сценарий слияния двух черных дыр, окруженных аккреционным диском. В резуль-

тате излучения гравитационных волн масса центрального объекта уменьшается, и аккреционный

диск испытывает возмущение. По результатам расчетов показано, что основным следствием этого

возмущения является формирование ударной волны, распространяющейся от центра к периферии

диска. Рассчитана кривая блеска и оценена продолжительность вспышки, в предположении, что

при снижении светимости до исходного значения вспышка заканчивается. Показано, что если масса

сливающейся двойной составляет 55 M_⊙ (как в случае события GW170814), то вспышка от ударной

волны вызовет рост болометрической светимости диска на 4-6 порядков, до 10⁴⁵ эрг/с (абсолютная

звездная величина -23.8^m). При расстоянии до источника 540 Мпк и разумных предположени-

ях о параметрах аккреционного диска оказывается, что видимая яркость вспышки в максимуме

спектральной плотности потока должна составить 12.8^m-14.2^m, длительность вспышки — несколько

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

В спектральной полосе инструмента EPIC обсерватории XMM-Newton или телескопа eROSITA

обсерватории Спектр-РГ (0.3-10 кэВ) рост светимости составит 3-4 порядка (7.5^m-10^m), вплоть

до 10⁴⁴ эрг/с, что соответствует видимой звездной величине примерно 17^m. В полосе наблюдения

инструмента IBIS обсерватории INTEGRAL (20 кэВ - 10 МэВ) светимость максимальна и составит

10⁴⁴-10⁴⁵ эрг/с, что соответствует видимому потоку 10-4 фотонов/см²/с/кэВ на длине волны

∼100 кэВ. Начиная с дальнего ультрафиолетового диапазона, в сторону б ´ольших длин волн, поярчание

практически отсутствует — на длине волны 10 эВ светимость вырастает примерно в два раза и

соответствует абсолютной величине -6^m и видимой 32^m.

DOI: 10.1134/S000462991901002X

1. ВВЕДЕНИЕ

ной ЧД. Очевидно, что такой диск не может об-

разоваться на стадии звезд главной последова-

К моменту написания этой работы (июль 2018 г.)

тельности или звезд-гигантов — предшественни-

гравитационный детектор LIGO зафиксировал

ков ЧД — поскольку сильный звездный ветер и

5 слияний черных дыр звездных масс. Слияния

взрывы сверхновых разрушат диск. Образование

сопровождались излучением

2.3%-5.8% энер-

аккреционного диска на стадии двойной ЧД пред-

гии покоя исходной двойной системы в виде

ставляется возможным, если время жизни систе-

гравитационных волн. Параметры сливающихся

мы до слияния достаточно велико, причем в этом

объектов приведены в табл. 1 (массы объектов,

случае могут быть эффективны те же механизмы,

доля теряемой массы, расстояние) [1-5].

которые приводят к формированию диска вокруг

В исследованиях черных дыр (ЧД) часто пред-

одиночной ЧД. Время жизни двойной ЧД до сли-

полагается, что они окружены аккреционными дис-

яния можно оценить как характерное время потери

ками. Источником вещества диска может быть

углового момента тесной двойной системой с кру-

молекулярное облако, через которое пролетают

говыми орбитами [6]:

ЧД, ветер от соседних астрофизических объек-

(

)₄

M3⊙

тов, остатки разрушенной приливным воздействи-

τ_GWR = (10⁸ лет)

(1)

ем звезды и т.п. Вполне реальна возможность су-

R_⊙ M₁M₂(M₁ + M₂)

ществования аккреционного диска и вокруг двой-

где a — межкомпонентное расстояние; M₁ и M₂ —

^*E-mail: bisikalo@inasan.ru

массы компонентов. Чтобы двойная ЧД с пара-

^**E-mail: zhilkin@inasan.ru

метрами объекта GW170814 [5] (последнее на се-

^***E-mail: kurbatov@inasan.ru

годняшний момент зарегистрированное событие)

БИСИКАЛО и др.

Таблица 1. Параметры двойных сливающихся ЧД по данным гравитационно-волновых наблюдений. Столбцы,

слева направо: название объекта, массы компонентов, относительное изменение массы, расстояние до объекта

Объект

M₁, M₂ [M_⊙]

ΔM/(M₁ + M₂) [%]

D [Мпк]

GW150914 [1]

36+5-4, 29+4-4

3.3-5.2

410+160-180

GW151226 [2]

14.2+8.3-3.7, 7.5+2.3-2.3

2.9-5.5

440+180-190

GW170104 [3]

31.2+8.4-6.0, 19.4+8.3-5.9

2.3-5.7

880+450-390

GW170608 [4]

12+7-2, 7+2-2

2.8-5.1

340+140-140

GW170814 [5]

30.5+5.7-3.0, 25.3+2.8-4.2

4.0-5.8

540+130-210

слилась за хаббловское время, 13.7 × 109 лет,

который может привести к возникновению сверх-

необходимо, чтобы ее начальное межкомпонентное

звуковых возмущений [11]. В-третьих, изменения

расстояние не превышало 48 R_⊙. Таким образом,

пространственно-временной метрики в гравитаци-

при любом разумном предположении об исходном

онных волнах может непосредственно вызывать

расстоянии между компонентами двойной ЧД, вре-

в диске механические напряжения, которые дис-

мя жизни системы будет достаточным для форми-

сипируют на вязкой шкале времени [12]. Все эти

рования циркумбинарного аккреционного диска.

явления приводят к электромагнитному отклику

аккреционного диска: поярчанию, квазипериоди-

Отличия аккреционного диска вокруг двойной

ческим вариациям светимости.

звезды от диска вокруг одиночной звезды неве-

лики [7, 8] и заключаются, главным образом, в

Большое число работ, посвященных этому яв-

довольно большом внутреннем радиусе диска — на

лению, касались сверхмассивных черных дыр в

уровне 1.5-2 расстояний между компонентами си-

ядрах сливающихся галактик (см., например, [13-

стемы. Процесс стремительного слияния двойной

16]). Так, если двойная ЧД суммарной массой

черной дыры начинается, когда расстояние между

10⁶ M_⊙ после слияния теряет 5% своей массы, то

компонентами составляет примерно a ≈ 6r_g, где

светимость аккреционного диска должна вырасти

r_g — гравитационный радиус аккретора, т.е. внут-

на порядок, до 10⁴³ эрг/с [15]. Эффект импульса

ренний радиус аккреционного диска r_in составляет

отдачи приводит к сравнимому росту светимости,

порядка 10r_g. Это позволяет предположить, что

однако его величина существенным образом зави-

диск вокруг двойной ЧД звездной массы является

сит от плохо определенных параметров сливаю-

щейся двойной [17, 18].

стандартным (нерелятивистским) и его свойства

могут быть описаны при помощи общепринятой

В качестве примера работы, относящейся к ЧД

модели α-диска Шакуры-Сюняева [9]. Использо-

звездных масс, можно отметить статью [19], где

вание этой модели позволяет сразу оценить свети-

были получены приближенные оценки светимо-

мость исходного диска: при довольно консерватив-

сти возмущенного диска как функции параметров

ных предположениях о его параметрах (темпера-

диска и двойной системы. Характерная величина

тура внутренней части диска порядка 10⁷ К, темп

светимости оказалась ≲10⁴³ эрг/с с максимумом

аккреции не превышает нескольких сотых от кри-

в рентгеновской области электромагнитного спек-

тического значения) светимость в мягком и среднем

тра.

рентгеновском диапазоне составит 1040±1 эрг/с.

Основной целью статьи является исследование

Во многих теоретических работах было по-

отклика аккреционного диска вокруг двойной ЧД

казано, что слияние двойной ЧД и последую-

звездной массы на уменьшение массы ЧД после

щее гравитационно-волновое излучение непо-

слияния. Для исследования этого эффекта была

средственно ведет к нескольким эффектам для

разработана численная модель и проведена се-

аккреционного диска. Во-первых, потеря массы

рия расчетов для дисков с доминирующим газо-

центральным компактным объектом после слияния

вым давлением. Был выполнен анализ примени-

приводит к нарушению (квази-)кеплеровского

мости модели, рассчитаны спектры электромагнит-

равновесия в диске и последующему возбуждению

ного излучения и получены оценки длительности

волн плотности большой амплитуды

[10]. Во-

вспышки. Как оказалось, основная часть энер-

вторых, в зависимости от отношения масс слива-

гии электромагнитного излучения высвечивается в

рентгеновском и гамма-диапазонах, где светимость

ющихся объектов и ориентации их собственных

угловых моментов, гравитационно-волновое из-

достигает 10⁴⁵ эрг/с. Были определены условия

лучение может быть асимметричным. В резуль-

регистрации электромагнитного сигнала обсерва-

тате вещество диска испытает импульс отдачи,

ториями XMM-Newton и INTEGRAL.

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

О ВОЗМОЖНЫХ ЭЛЕКТРОМАГНИТНЫХ ПРОЯВЛЕНИЯХ

Статья организована следующим образом. В

верна и для любой двойной ЧД. Действительно,

следующем разделе дано описание параметров ста-

гравитационный радиус аккретора может быть за-

ционарных аккреционных дисков, используемых в

писан в виде

данной работе. В третьем разделе описана модель

2GM

r_g

для численного решения нестационарной задачи.

r_g =

или

= 4.2 × 10-6m.

(3)

c²

R_⊙

В четвертом разделе приведены расчеты спектра

электромагнитного излучения и рассмотрены воз-

Двойная ЧД начинает сливаться, когда расстояние

можности наблюдения моделируемых систем. Вы-

между ее компонентами составляет примерно a ≈

воды содержатся в пятом разделе.

≈ 6r_g. Вблизи внутренней границы аккреционного

диска течение газа возмущено приливной силой

со стороны двойной ЧД. Положение r_in можно

2. ПОСТАНОВКА ЗАДАЧИ

оценить из условия равенства приливной силы и

Одной из универсальных моделей стационарной

гравитационной силы от точечного источника той

же массы:

дисковой аккреции является стандартная модель

α-диска Шакуры-Сюняева [9]. Эта модель описы-

2GM

a/2

(4)

вает геометрически тонкий диск, в котором дисси-

r²

r_in

r²

пация и перенос углового момента обеспечиваются

турбулентной вязкостью. Турбулентность харак-

Отсюда следует, что r_in = 6r_g, что дает оценку

теризуется единственным параметром α, значение

нижней границы для внутреннего радиуса диска,

которого не превышает единицы. Темп диссипации

так как газодинамические эффекты приводят к уве-

углового момента влияет как на скорость аккре-

личению r_in. Ниже мы будем полагать r_in = 10r_g.

ции, так и на темп энерговыделения в диске, и, в

Эта оценка согласуется с результатами [25] для

конечном итоге, задает распределение плотности и

двойной системы с круговыми орбитами компонен-

температуры по диску.

тов и отношением их масс, близким к единице. Для

Структура аккреционного диска в стандарт-

m = 55 будем иметь r_in = 1.6 × 10⁸ см.

ной модели Шакуры-Сюняева зависит от следу-

В зависимости от величины темпа аккреции m

ющих параметров: безразмерная масса аккрето-

(а так же от α и m, но в меньшей степени), в

ра m ≡ M/M_⊙, безразмерный темп аккреции m˙ ≡

диске можно выделить несколько зон, в которых

≡M˙ / Mcr, эффективность отвода углового момен-

преобладает либо газовое, либо радиативное дав-

та (параметр турбулентности) α; эффективность

ление и которые различаются механизмами погло-

высвечивания гравитационной энергии (эффектив-

щения излучения. В данной работе мы ограничимся

ность аккреции) η; внутренний радиус аккрецион-

рассмотрением дисков, в которых в равновесном

ного диска r_in.

состоянии преобладает газовое давление. Соглас-

В рассматриваемом случае двойной ЧД перед

но [9], это так называемые зоны “B” и “C”. Условие,

при котором в диске отсутствует область, где су-

слиянием масса аккретора равна суммарной мас-

се компонентов двойной системы, M ≡ M₁ + M₂.

щественно давление излучения (зона “A”), выглядит

Критический темп аккреции определяется величи-

так¹:

ной эддингтоновской светимости [9]:

m< m_low ≡

(5)

M_cr

170 (αm)1/8

0.06

= 3 × 10-8

(2)

M_⊙/год

Зависимость от произведения αm довольно сла-

бая. Так, для m = 55 и

0.001 ≤ α ≤ 1 имеем

Параметр η в случае одиночной невращающейся

0.025 ≲ m˙_low ≲ 0.059.

черной дыры примерно равен 0.06-0.08 [20]. Учет

Характерная шкала температур в зонах “B” и

вращения повышает эффективность стационарной

“C” порядка 10⁶-10⁷ К, а плотности — порядка

аккреции до ∼0.32 [21, 22]. В сценарии эпизоди-

ческой аккреции эффективность может достигать

1 г/см³. В работе Шакуры и Сюняева были получе-

0.43 [23]. Значение этого параметра, однако, не

ны радиальные профили объемной плотности газа,

критично в рассматриваемой модели, т.к. опреде-

усредненной по толщине диска, а также поверх-

ляющую роль играет величина безразмерного тем-

ностной температуры (на единичной оптической

па аккреции m.

глубине). Температура во внутренних слоях дис-

ка зависит от эффективности переноса тепловой

Предположение о размере внутреннего радиуса

диска было получено из численного моделирования

¹Множитель в правой части выражения (5) отличается

аккреции в конкретной системе двойной ЧД [24].

соответствующего множителя в формуле (2.18) из рабо-

Из расчета следовало, что внутренний радиус дис-

ты [9]: во-первых, в оригинальной статье полагалось rin =

ка приблизительно равен двум межкомпонентным

= 3rg; во-вторых, мы уточнили этот множитель числен-

расстояниям. Нетрудно показать, что эта оценка

ным расчетом.

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

БИСИКАЛО и др.

энергии. В условиях зон “B” и “C” поглощение

Зададим равновесные профили температуры и

излучения определяется как томсоновским рассе-

плотности как степенные функции радиальной ко-

янием, так и свободно-свободными переходами.

ординаты:

Их сечения равны соответственно σ_T = 0.40 см²/г

)-kt

⁽r

и σ_ff = 0.11N/T7/2 см²/г, где N — концентрация.

T =T_∗

(12)

r_∗

Оптическую толщину диска в вертикальном на-

)-kd

правлении можно оценить как τ = (σ_T + σ_ff)Σ ≳

⁽r

ρ=ρ_∗

(13)

≳ 10⁵, причем σ_T ≫ σ_ff. В условиях большой оп-

r_∗

тической толщины поверхностная температура T_s

(на оптической глубине порядка единицы) и темпе-

где r_∗ = r_in. Характерные масштабы температуры

и плотности в этих распределениях зависят от

ратура в середине диска, T , связаны соотношени-

параметров α-модели. Для зоны “B” имеем [9]

ем [26]

3τ

T_∗ = 8.3 × 10⁷(αm)-1/5m˙2/5 =

(14)

T⁴ ≈

T^4s.

(6)

= 2.7 × 10⁷ К, k_t = 9/10,

Таким образом, если перенос тепловой энергии

ρ_∗ = 0.76(αm)-7/10 m2/5 =

(15)

осуществляется путем лучистой теплопроводности,

то T/T_s ≳ 10.

= 0.33 г/см³, k_d = 33/20.

Если же за перенос энергии в вертикальном на-

В этом выражении подразумевается температура

правлении отвечает конвекция, перепад темпера-

внутри диска, а не на его поверхности. Множи-

туры можно оценить, полагая распределение газа

тель (11) уже учтен. Для зоны “C”

изэнтропическим:

T_∗ = 2.9 × 10⁷(αm)-1/5 m3/10 =

(16)

⁽ρ)γ-1,

(7)

= 1.3 × 10⁷ К, k_t = 3/4,

T_s

ρ_s

ρ_∗ = 3.8(αm)-7/10 m11/20 =

(17)

где ρ и ρ_s — объемная плотность газа соответ-

ственно в середине и на поверхности диска; γ —

= 1.0 г/см³, k_d = 15/8.

показатель адиабаты. Согласно оценкам [9], ρ/ρ_s ≈

≈ 2, откуда для одноатомного газа (γ = 5/3) полу-

Как видно, профили в зонах “B” и “C”, хотя

чаем T/T_s ≈ 1.6.

и имеют различный наклон, но нормированы на

схожие значения температуры. Чтобы не услож-

Характерное время лучистой теплопроводности

нять модель, далее будет подразумеваться, что диск

есть

состоит только из зоны “B”.

τ c_s

t_rad ∼ τ

∼

(8)

Диапазон допустимых значений величин T_∗ и ρ_∗

Ω c

довольно широк и ограничен только условием (5).

где τ — оптическая толщина диска (см. выше); H ∼

Если положить, что α = 0.001 есть наименьшее ре-

∼ c_s/Ω — полутолщина диска; c_s — скорость зву-

алистичное значение параметра турбулентности, то

ка; Ω — кеплеровская частота. Время конвекции

подставляя m = 55 и m˙ = m˙_low в выражения (14),

оценим как

(15), получим T_∗ = 4.8 × 10⁷ К, ρ_∗ = 1.87 г/см³ —

эти значения температуры и плотности следует

t_conv ∼

∼

(9)

рассматривать как верхнюю оценку соответствую-

αc_s

αΩ

щих величин. Нижнюю оценку получим, положив

Их отношение показывает, какой из процессов

α = 1, m= 0.01. В этом случае T∗ = 5.9 × 106 К,

переноса энергии более эффективен,

ρ_∗ = 7.29 × 10-3 г/см³.

(

)1/2

Распределение углового момента по диску близ-

t_rad

c_s

∼ ατ

∼ 400α

(10)

ко к кеплеровскому, при этом вклад давления газа

t_conv

10⁷ К

и радиальной составляющей кинетической энергии

Видно, что в зонах “B” и “C” конвекция более

в радиальный баланс сил — мал (порядка квадрата

эффективно переносит энергию. Исходя из оценок,

отношения толщины диска к радиальному масшта-

сделанных выше, далее будем полагать, что значе-

бу [27, 28]). Тепловой структурой диска, однако,

ния температуры в середине и на поверхности диска

нельзя пренебречь, она определяет интенсивность

относятся друг к другу как

ударной волны и важна для оценки реакции диска

на потерю массы аккретором. Напротив, радиаль-

≡ 2.

(11)

ная составляющая скорости представляется несу-

T_s

щественной в этом процессе.

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

О ВОЗМОЖНЫХ ЭЛЕКТРОМАГНИТНЫХ ПРОЯВЛЕНИЯХ

Суммируем постановку задачи. Дан осесиммет-

В области диска, где давлением излучения можно

ричныый аккреционный диск в гравитационном по-

пренебречь по сравнению с газовым давлением,

ле точечной массы. Распределение термодинами-

уравнение состояния определяется соотношения-

ческих величин по диску соответствует зоне “B”

ми:

модели Шакуры и Сюняева [9]. Центробежная сила

P = P(ρ,T) = AρT,

(25)

в диске полностью уравновешивает силу грави-

тации и газового давления (радиальная скорость

ε = ε(ρ,T) =

(26)

всюду равна нулю). В начальный момент вре-

γ-1

мени масса аккретора уменьшается на 5%. Это

состояние диска неравновесно. Дальнейшая эво-

где A = k_B/m, k_B — постоянная Больцмана, m =

люция состояния рассчитывается в приближении

= 0.5m_p — средняя молекулярная масса, m_p —

бездиссипативной газовой динамики. Нас интере-

масса протона, γ = 5/3 — показатель адиабаты,

сует решение в средней плоскости диска, при этом

а средний молекулярный вес для полностью

вертикальной структурой диска мы пренебрегаем

ионизованной водородной плазмы принят равным

(основания для этого будут даны в п. 5). В каче-

0.5.

стве начальных условий примем распределения для

Положим начальные распределения температу-

температуры и плотности (12)-(15).

ры и плотности, T(r,0) и ρ(r,0), степенными функ-

В следующем разделе будут даны описание

циями радиуса, (12) и (13). Из этих определений и

уравнений и процедура их численного решения.

из уравнений состояния (25), (26) можно получить

соответствующие распределения для давления и

удельной внутренней энергии:

3. ЧИСЛЕННАЯ МОДЕЛЬ

)-kd-k

⁽r

В осесимметричном случае в цилиндрических

P (r, 0) = ρ_∗c_∗

(27)

r_∗

координатах вместо радиальной координаты r

удобно ввести массовую лагранжеву координату q,

c∗2

которая удовлетворяет соотношению dq = r²ρdr.

ε(r, 0) =

⁽r)-kt ,

(28)

В результате уравнения газовой динамики, опи-

γ-1

r_∗

сывающие эволюцию аккреционного диска после

где обозначено c∗2 = AT_∗. Считая, что начальная

слияния черных дыр, можно записать в виде:

радиальная скорость равна нулю, из условия ради-

⁽1⁾

∂

ального равновесия в диске

(rv_r) ,

(18)

dt ρ

∂q

∂P

l²

(29)

dv_r

∂P

v2ϕ

ρ ∂r

r²

r³

= -r

- (1 - ξ)

(19)

∂q

r²

можно получить

[

]

dv_ϕ

v_rv_ϕ

)

(20)

⁽r

⁽r)2-kt

l² = r∗2c_∗

μ²

- (k_d + k_t)

r_∗

dε

∂

= -P

(rv_r) ,

(21)

(30)

∂q

или

=v_r.

(22)

[

(

)-1

]1/2

⁽r)-kt

v_ϕ = c_∗ μ²

- (k_d + k_t)

Здесь ρ — плотность, v_r — радиальная скорость,

r_∗

v_ϕ — азимутальная скорость, P — давление, ε —

(31)

удельная внутренняя энергия, M — полная масса

двойной черной дыры до их слияния, ξ — доля мас-

где обозначено

сы черных дыр, излученная в виде гравитационных

(

)-1/2

⁽ GM )1/2

T_∗

волн. Отметим, что вместо уравнения (18) можно

μ=

≡ 165

(32)

использовать эквивалентное уравнение

r_∗c∗2

10⁷ К

∂r

По своему смыслу параметр μ равен числу Маха на

(23)

радиусе r = r_∗.

∂q

Для последующих оценок нам также понадобит-

Кроме того, из уравнений (20) и (22) следует закон

ся характерная полутолщина диска. Определим ее

сохранения удельного углового момента l = rv_ϕ,

как вертикальную шкалу равновесного диска:

= 0.

(24)

H =

(AT )1/2 =

(33)

v_φ

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

БИСИКАЛО и др.

[

)

kt-3

Таблица 2. Некоторые характеристики дисков. Столб-

⁽r

⁽r)-2]-1/2

=r_∗ μ²

− (k_d + k_t)

цы: число Маха, температура в центральной части

r_∗

невозмущенного диска, скачок температуры за ударной

волной

При построении вычислительного алгоритма

удобно перейти к безразмерным переменным.

T_∗ [К]

T/T(t = 0)

Это достигается путем выделения характерных

100

2.7 × 10⁷

7.5

масштабов по схеме: f → f₀f, где f — какая-либо

величина, а f₀ — ее размерный масштаб. Чтобы не

150

1.2 × 10⁷

15.5

загромождать запись, безразмерные переменные

200

0.43 × 10⁷

будем далее обозначать теми же символами, что и

сами исходные размерные переменные. В качестве

размерных масштабов выберем следующие вели-

статьи, зарегистрированному событию слияния

чины:

двойной ЧД. Графики безразмерной плотности,

r₀ = r_∗, ρ₀ = ρ_∗, v₀ = c_∗,

(34)

температуры, радиальной и тангенциальной ско-

ростей приведены на рис. 1-4. Все три расчета

t₀ = r_∗/c_∗, P₀ = ρ_∗c∗2, ε₀ = c∗2, q₀ = r∗2ρ_∗.

обнаружили качественно схожие свойства течения:

В безразмерных переменных уравнение для ра-

в центральной области диска возникает одна,

либо две ударные волны, в зависимости от μ.

диальной скорости (19), а также уравнения состо-

Волна распространяется по диску с замедлением.

яния (25), (26) примут вид:

Плотность на ударной волне испытывает скачок в

dv_r

∂P

v2ϕ

(1 - ξ)μ

три и более раз, температура — примерно от 7.5

= -r

(35)

∂q

r²

до 27 раз, по сравнению с исходным значением

(см. табл. 2). Радиальная скорость испытывает

осцилляции с амплитудой порядка скорости звука

P = ρT, ε =

(36)

и более. Тангенциальная скорость меняется слабо.

γ-1

Характеристики течения перед ударной волной

Остальные уравнения после соответствующих пе-

остаются практически постоянными на протяже-

реобозначений своего вида не изменят.

нии времени расчета.

Для численного решения этих уравнений ис-

Течение газа в возмущенной части диска пред-

пользовалась неявная полностью консервативная

ставляет собой совокупность волн плотности, тем-

разностная схема Самарского-Попова [29], неко-

пературы и радиальной скорости большой ампли-

торые детали которой описаны в Приложении.

туды. Для μ = 200 и 150 на внешней границе воз-

Следует подчеркнуть, что в этой схеме выполня-

мущенной области возникает одна ударная вол-

ются не только разностные аналоги законов со-

на с предвестником (каустика на распределении

хранения массы, импульса и энергии, но и допол-

плотности и температуры), при μ = 100 каустика

нительные соотношения, описывающие баланс по

превращается в ударную волну. Мы также провели

определенным видам энергии. Кроме того, в нашем

расчеты для чисел Маха μ = 50 и μ = 10 (не по-

случае в схеме выполняется разностный аналог за-

казаны на рисунках). Интересно то, что для малых

кона сохранения удельного углового момента (24).

значений μ в диске имеется лишь одна ударная

В безразмерной формулировке задача имеет

волна с предвестником.

только два параметра: ξ и μ. Долю теряемой

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

массы зададим, пользуясь оценкой, полученной

функция времени не зависит от числа Маха в диске,

по гравитационно-волновым наблюдениям [5], ξ =

μ, и допускает простую аппроксимацию² (рис. 5):

= 0.05. В предыдущем разделе мы показали, что в

реалистичном диапазоне значений параметра тур-

r_D = 1.76r_∗t2/3.

(37)

булентности α и безразмерного темпа аккреции

Отсюда можно получить скорость ударной волны,

m характерные значения температуры T_∗ лежат в

)-1/2

интервале 5.9 × 10⁶-4.8 × 10⁷ К. Этому интервалу

dr_D

⁽r_D

= 1.2r_∗t-1/3

= 1.59r_∗

(38)

соответствует интервал чисел Маха μ = 75-215. В

r_∗

численном расчете мы задавали определенные зна-

чения μ и по ним, в соответствии с (32), получали

²Несовпадение графиков для μ = 100 и остальных связано

T_∗, c_∗ и t_∗.

с тем, что в программном коде осуществляется поиск УВ

со стороны правой (внешней) части расчетной области. В

Мы рассчитали три численные модели: для

то же время для этого значения числа Маха в течении воз-

μ = 100,150,200. Эти значения μ, а также до-

никают две ударных волны (см. рис. 2), причем положение

ля теряемой массы ξ, соответствовали объекту

внутренней совпадает с положением УВ для μ = 150 и

GW170814 [5] — последнему, на момент написания

μ = 200.

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

О ВОЗМОЖНЫХ ЭЛЕКТРОМАГНИТНЫХ ПРОЯВЛЕНИЯХ

ρ/ρ_*

ρ/ρ(t = 0)

10⁰

μ = 100, t = 49.7 [s]

μ = 100

μ = 150, t = 56.7 [s]

μ = 200, t = 57.2 [s]

10^-1

10^-2

μ = 150

10^-3

μ = 200

10^-4

10^-5

10⁰

10¹

10²

10³

10¹

10²

r/r_*

Рис. 1. Радиальное распределение плотности в моделях с различными значениями числа Маха μ. Моменты времени

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

T/T_*

T/T(t = 0)

10²

μ = 100, t = 49.7 [s]

μ = 100

μ = 150, t = 56.7 [s]

10¹

μ = 200, t = 57.2 [s]

10¹

10⁰

μ = 150

10¹

10^-1

10⁰

μ = 200

10^-2

10¹

10⁰

10^-3

10⁰

10¹

10²

10³

10¹

10²

r/r_*

Рис. 2. Радиальное распределение температуры. Обозначения те же, что и на рис. 1.

4. СПЕКТР ЭЛЕКТРОМАГНИТНОГО

температуры примерно 2 × 10⁸ К, что соответствует

ИЗЛУЧЕНИЯ И КРИВЫЕ БЛЕСКА

тепловой энергии, в расчете на одну частицу, по-

рядка 17.2 кэВ. В условиях немагнитного аккреци-

Ударная волна, возникающая в диске, приводит

онного диска спектр электромагнитного излучения

к росту температуры в 7.5-27 раз, в зависимости

в этой области энергий определяется двумя про-

от числа Маха (рис. 2, табл. 2). В результате этого

цессами [30]: тормозным излучением электронов

газ во внутренней части диска разогревается до

и комптоновским рассеянием. Последний меха-

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

БИСИКАЛО и др.

v_r/c_*

μ = 100, t = 49.7 [s]

μ = 100

μ = 150, t = 56.7 [s]

μ = 200, t = 57.2 [s]

-1

-2

μ = 150

-1

-2

-1

μ = 200

-1

-2

10⁰

10¹

10²

10³

10¹

10²

r/r_*

Рис. 3. Радиальное распределение радиальной скорости. Обозначения те же, что и на рис. 1.

v_ϕ/c_*

/v_ϕ(t = 0)

103

μ = 100, t = 49.7 [s]

μ = 100

μ = 150, t = 56.7 [s]

1.00

μ = 200, t = 57.2 [s]

0.95

102

0.90

μ = 150

1.00

0.95

101

0.90

μ = 200

1.00

0.95

100

0.90

10⁰

10¹

10²

10³

10¹

10²

r/r_*

Рис. 4. Радиальное распределение тангенциальной скорости. Обозначения те же, что и на рис. 1.

низм может приводить к росту энергии фотонов

зрачен [30, 32]. Оценим оптическую толщину диска:

и искажению спектра — снижению интенсивности

τ_ν = (σ_T + σff,ν)Σ.

(39)

на частотах планковского максимума и появлению

“хвоста” в области высоких энергий [31, 32]. Од-

Подставляя выражение для сечений поглощения,

нако комптонизация спектра эффективна, только

например из работы [31], можно показать, что в

если слой, через который проходит излучение, про-

центральной части диска (ρ ∼ 1 г/см³, T ∼ 10⁸ К,

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

О ВОЗМОЖНЫХ ЭЛЕКТРОМАГНИТНЫХ ПРОЯВЛЕНИЯХ

r_D [cm]

L [erg/s]

10⁴⁶

μ = 100

μ = 150

10¹⁰

10⁴⁵

μ = 150

μ = 200

1.76r_*t^2/3

10⁴⁴

10⁴³

10⁹

10⁴²

10⁴¹

10⁸

10⁴⁰

10³⁹

10⁷

10³⁸

0.1

100

0.01

0.1

t [s]

Рис. 5. Положение ударной волны как функция време-

Рис. 6. Болометрические кривые блеска в моделях с

ни. Серым цветом показана аппроксимация.

различными значениями числа Маха.

H ∼ 10⁶ см) оптическая толщина во всем спектре

табл. 2). Из оценки максимальной светимости и с

много больше единицы и на длинах волн короче

помощью формулы (32) получаем

10 нм определяется комптоновским рассеянием.

)₂

( μ

Это свойство выполняется и во внешних частях

f = 27

(42)

200

диска. Искажение спектра в оптически тонком слое

вблизи поверхности диска также несуществен-

Напомним, что в данной работе мы основыва-

но, поскольку отношение геометрической толщины

лись на гидродинамической модели, без учета ра-

этого слоя к толщине диска (его можно оценить как

диативных процессов. На рис. 2 и 6 видно, что по-

сле ударного скачка температура газа практически

τ-1ν) пренебрежимо мало на всех частотах. Таким

образом, можно полагать, что излучение с поверх-

не меняется. Это приближение заведомо неспра-

ведливо на больших временах, т.к. процесс ради-

ности диска имеет всюду планковский спектр³:

ативного охлаждения должен приводить к осты-

2hν

ванию диска после прохождения ударной волны и

B_ν =

(40)

постепенному установлению нового равновесного

c²

e^hν/kBTs - 1

состояния, где темп высвечивания энергии будет

где T_s — поверхностная температура (см. п. 2 и

равен скорости вязкого нагрева. Оценим время

соотношение (11)).

радиативного охлаждения как время диффузии из-

На рис. 6 показаны кривые изменения боломет-

лучения:

рической светимости диска во времени, рассчитан-

ной как

t_rad ∼ (σ_T + σ_ff)Σ

(43)

∫∞

где σ_T = 0.40 см²/г и σ_ff = 0.11N/T7/2 см²/г; N —

L(t) = 2π drrσ_SBT^4s(t, r).

(41)

концентрация; Σ = 2ρH. Например, в окрестности

r∗

внутреннего радиуса диска имеем ρ ∼ 1 г/см³, T ∼

Примечательно то, что разные значения параметра

∼ 10⁸ К и H ∼ 10⁶ см, откуда получаем t_rad ≈

μ привели в итоге к одной и той же светимости дис-

≈ 31.6 с. Более общее выражение имеет вид (при

ка, порядка 10⁴⁵ эрг/с. Аналитически светимость

таких температурах сечение σ_ff не вносит вклада в

поглощение)

можно оценить как L ∼ πr∗2σ_SBf⁴T∗4, где f — ска-

чок температуры на ударной волне (третья колонка

2σ_Tρ_∗H∗2

⁽r)3-(kd+kt)

t_rad =

(44)

r_∗

³Комптонизация электромагнитного излучения эффектив-

на в зоне “A” стандартного α-диска [9]. В используемой

Характерную ширину горячей области за ударной

нами модели диска зона “A” отсутствует, а в зоне “B”

волной можно оценить как

поверхностная плотность на 4-5 порядков больше, чем

можно оценить для зоны “A”.

Δr_D = Dt_rad,

(45)

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

БИСИКАЛО и др.

λ [nm]

10⁴

10³

10²

10¹

10⁰

10^-1

10^-2

10^-3

10^-27

10⁵⁴

μ = 100

10²⁷

10^-1

μ = 150

10⁵³

10^-28

μ = 200

10²⁶

10^-2

10^-29

10⁵²

10²⁵

10^-3

10^-30

10⁵¹

10²⁴

10^-4

10^-31

10⁵⁰

10²³

10^-5

10^-32

10⁴⁹

10²²

10^-6

10^-33

10⁴⁸

10²¹

10^-7

10^-34

10⁴⁷

10²⁰

10^-8

10⁴⁶

10¹⁴

10¹⁵

10¹⁶

10¹⁷

10¹⁸

10¹⁹

10²⁰

10²¹

ν [Hz]

10⁰

10¹

10²

10³

10⁴

10⁵

10⁶

hν [eV]

Рис. 7. Спектр электромагнитного излучения. Внутренние вертикальные шкалы: спектральная плотность потока (слева)

и оценка потока фотонов в расчете на единичный логарифмический интервал частот (справа). Внешние вертикальные

шкалы соответствуют источнику, удаленному на расстояние 540 Мпк: спектральная плотность потока через единичную

площадку (слева) и оценка потока фотонов через единичную площадку (справа). Серым цветом обозначены спектры

невозмущенных дисков для соответствующих значений μ.

где D — скорость ударной волны (38). Подставив

На рис. 7 видно, основная часть энергии излучается

r_D вместо r в выражения (44) и (45), можно

в рентгеновском и гамма-диапазоне.

увидеть, что ширина горячей области Δr_D слабо

Для гамма-диапазона удобно использовать

зависит от времени (как t-0.075) и приближенно

оценку числа фотонов, излучаемых за единицу

равна

времени в единичном логарифмическом интервале

Δr_D ≈ 42.4r_∗.

(46)

частот (правая шкала на рис. 7):

∫

Это означает, что эффекты радиативного охлажде-

Fν′

ния не должны сильно изменить распределение

dν^′

(49)

hν^′

температуры на рис. 2. Время уменьшения темпе-

ратуры до исходного значения оценим следующим

∫

образом:

Fν′

= ln 10

d(lg ν^′)

∼ ln 10

t_cool = t_rad

(47)

T (t = 0)

На внешних вертикальных шкалах отложены ве-

Если температура за ударной волной возрастает

личины потока через площадку 1 см² при условии,

в 7.5 - 27 раз, то время спадания кривой блеска

что наблюдаемый объект удален на расстояние

составит от 4 до 14 минут.

540 Мпк (система GW170814 [5]). Для этого поток

Спектральная плотность потока излучения с

корректировался на множитель

одной стороны диска, в перпендикулярном направ-

лении:

Ω=

= 3.6 × 10-55 см-2.

(50)

∞

(540 Мпк)²

∫

Fν = 2π² drrBν[Ts(r)].

(48)

Аналогичным образом можем приближенно со-

r∗

считать энергию, излучаемую за единицу времени

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

О ВОЗМОЖНЫХ ЭЛЕКТРОМАГНИТНЫХ ПРОЯВЛЕНИЯХ

λ [nm]

10⁴

10³

10²

10¹

10⁰

10^-1

10^-2

10^-3

10^-7

-30

10⁴⁷

μ = 100

10^-8

μ = 150

10⁴⁶

μ = 200

-25

10^-9

10⁴⁵

10^-10

10⁴⁴

-20

10^-11

10⁴³

10^-12

10⁴²

-15

10^-13

10⁴¹

10^-14

10⁴⁰

-10

10^-15

10³⁹

10^-16

10³⁸

-17

-5

10³⁷

10¹⁴

10¹⁵

10¹⁶

10¹⁷

10¹⁸

10¹⁹

10²⁰

10²¹

ν [Hz]

10⁰

10¹

10²

10³

10⁴

10⁵

10⁶

hν [eV]

Рис. 8. Оценка светимости в расчете на единичный логарифмический интервал частот. Внутренние вертикальные

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

соответствуют источнику, удаленному на расстояние 540 Мпк: поток через единичную площадку детектора (слева) и

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

значений μ.

в единичном логарифмическом интервале частот

инструмента EPIC обсерватории XMM-Newton

(рис. 8):

(уровень 5σ с экспозицией 10⁵ с) [35].

∫

Поток излучения 10⁴⁵ эрг/с и выше, на частотах

dν^′F_ν′ ∼ ln 10 · νF_ν .

(51)

10-300 кэВ, могут обеспечить “горячие” α-диски

(∼10⁹ К) вокруг ЧД звездных масс [32, 36], либо

аккрецирующая сверхмассивная ЧД. Наблюдае-

Эту величину можно рассматривать как оценку

мая плотность потока в жестком рентгеновском

светимости в полосе частот шириной один порядок

диапазоне будет лишь на один-два порядка мень-

величины.

ше светимости крабовидной туманности (1 краб =

Рис. 7 и 8 показывают, что основная часть

= 2.4 × 10-8 эрг/см²/с).

энергии излучения лежит в диапазоне среднего и

В гамма-диапазоне поток излучения на много

жесткого рентгеновского излучения (1- 100 кэВ),

порядков превышает поток от невозмущенного

а также гамма-излучения (≳100 кэВ).

диска и сравним с потоком от мягких гамма-

В среднем рентгеновском диапазоне, до 10 кэВ,

источников, ≳10²⁶ эрг/с/Гц [37]. На длине вол-

светимость газа, разогретого ударной волной, рав-

ны

100

кэВ наблюдаемый поток составляет

на 10⁴³-10⁴⁵ эрг/с, что на 2-4 порядка (в за-

10-2

фотонов/см²/с, в расчете на единичный

висимости от μ) превышает светимость невозму-

логарифмический интервал длин волн (dex), или

щенного диска, которая, в свою очередь, сопо-

10-4 фотонов/см²/с/кэВ. Отметим, что предель-

ставима со светимостью рентгеновских пульса-

ная чувствительность инструмента IBIS обсерва-

ров [20] и ультраярких рентгеновских источни-

тории INTEGRAL на этой длине волны составляет

ков [33, 34]. После коррекции на расстояние до

2.85 × 10-6 фотонов/см²/с/кэВ (уровень 3σ с

источника GW170814 [5], 540 Мпк, наблюдаемый

экспозицией 10⁵ с) [38, 39].

поток от возмущенного диска в этой области спек-

тра составит ≳10-9 эрг/см²/с, что на несколь-

В экстремальном ультрафиолетовом диапазоне,

ко порядков превышает порог чувствительности

10-10² эВ, поярчание диска может превышать по-

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

БИСИКАЛО и др.

рядок величины или 2.5^m и достигать 3 × 10⁴⁰-3 ×

аккретора. Это различие сильнее во внутренней

× 10⁴¹ эрг/с. В дальнем ультрафиолете поярчание

части диска, что должно приводить к распростра-

нению по диску сильного возмущения.

не превышает двух раз. На б ´ольших длинах волн

поярчание практически отсутствует. В ближнем

Мы применили эти соображения к объекту

ультрафиолетовом диапазоне светимость диска на

GW170814 — двойной ЧД суммарной массой око-

4-5 порядков ниже светимости типичной сверхно-

ло 55 M_⊙, удаленной на расстояние 540 Мпк. В ка-

вой (10⁴¹ эрг/с).

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

Оценки светимости, приведенные выше, явля-

стандартная модель Шакуры и Сюняева. В диске

преобладало газовое давление, для этого темпера-

ются весьма приближенными также и вследствие

тура исходного диска полагалась достаточно низ-

неточно определенного расстояния до двойной си-

стемы. Так, минимальная и максимальная оценки

кой, ≲3 × 10⁷ К. Мы нашли, что из внутренней об-

расстояния до источника GW170814 [5] различа-

ласти аккреционного диска начинает распростра-

ются в два раза (см. табл. 1). При этом изменению

няться возмущение, которое быстро становится

расстояния в два раза соответствует изменение

ударной волной. Температура за ударной волной

видимой яркости в четыре раза или изменению

испытывает скачок в 7.5-27 раз, в зависимости от

видимой величины на 1.5^m.

начальной температуры диска.

Расчет болометрических кривых блеска пока-

зал, что светимость диска возрастает на 4-6 по-

5. ОБСУЖДЕНИЕ И ВЫВОДЫ

рядков, до 10⁴⁵ эрг/с, что соответствует абсо-

В численной модели, описанной в п. 3, мы пре-

лютной звездной величине -23.8^m. Максимум по-

небрегли вертикальной структурой диска. Сейчас,

тока излучения лежит в рентгеновском и гамма-

по результатам расчетов, мы можем убедиться в

диапазонах. В спектральной полосе инструмента

справедливости этого приближения. После того,

EPIC обсерватории XMM-Newton или инструмен-

как черная дыра потеряла долю ξ своей массы, ве-

та eROSITA обсерватории Спектр-РГ рост све-

щество диска приобрело ускорение в вертикальном

тимости, по сравнению с невозмущенным дис-

направлении величиной vz ∼ ξHΩ2. Характерное

ком, составит 3-4 порядка (7.5^m-10^m в абсо-

время отклика диска есть

лютных величинах); это соответствует видимой

)3/2

звездной величине 17^m. В полосе наблюдения ин-

1/2

f1/2c_s

⁽r

струмента IBIS обсерватории INTEGRAL све-

t_z ∼

∼t_∗

(52)

v_z

ξμ r_∗

тимость максимальна и соответствует видимому

потоку 10-4 фотонов/см²/с/кэВ на длине волны

где f — фактор повышения температуры (см.

∼100 кэВ. В дальнем ультрафиолетовом диапа-

предыдущий раздел); здесь мы использовали

зоне, и в сторону б ´ольших длин волн, возмущение

определения из п. 3. Это время нужно сравнить

диска практически не приводит к его поярчанию:

со временем распространения ударной волны, t_D ∼

на длине волны 10 эВ светимость вырастает при-

∼ r/D. С учетом оценок для f и D из п. 4 получаем

мерно в два раза и соответствует видимой величине

t_z

32^m. Оценки видимых потоков в рентгеновском

∼

(53)

t_D

ξμ

диапазоне превышают порог чувствительности ин-

струмента EPIC на несколько порядков величины

Для ξ = 0.05 и μ ∼ 100 видно, что релаксация диска

(до четырех, в зависимости от экспозиции). Поток

в вертикальном направлении не оказывает влия-

в гамма-диапазоне примерно на полтора порядка

ния на распространение ударной волны. Можно

превышает порог чувствительности инструмента

показать, однако, что время t_z оказывается одного

IBIS.

порядка со временем t_cool, поэтому расчет излуче-

ния на больших временах потребует более полного

Использование реалистичных параметров ак-

учета вертикальной структуры диска.

креционного диска приводит к оценке времени спа-

дания кривой блеска порядка нескольких минут.

В данной работе мы рассчитали реакцию цир-

Интересно то, что болометрическая светимость

кумбинарного аккреционного диска на уменьше-

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

ние массы двойной черной дыры. Идея состоя-

определяется, по-видимому, массой аккретора и

ла в том, что сливающаяся черная дыра теря-

долей теряемой массы. При этом форма спектра

ет около 5% суммарной массы-энергии посред-

чувствительна к температуре: рост начальной тем-

ством излучения гравитационных волн. В результа-

пературы диска приводит к смещению спектра в

те этого вещество изначально равновесного квази-

коротковолновую область.

кеплеровского диска приобретает избыточный им-

пульс в направлении от аккретора. Импульс вы-

В следующих работах мы планируем продол-

зван дисбалансом между центробежной силой и

жить исследование эффекта потери массы слива-

уменьшившейся гравитационной силой со стороны

ющейся черной дырой для более полной модели

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

О ВОЗМОЖНЫХ ЭЛЕКТРОМАГНИТНЫХ ПРОЯВЛЕНИЯХ

аккреционных дисков, с учетом радиативного дав-

Tn+1i+1/2

ления и выхолаживания газа.

εn+1i+1/2 =

(61)

γ-1

Авторы благодарят А.В. Тутукова, Я.Н. Пав-

Πn+1i+1/2 = Pn+1i+1/2 + ωn+1i+1/2,

(62)

люченкова и О.Ю. Малкова за ценные замечания.

(

)

Работа Курбатова Е.П. поддержана Программой

ωn+1i+1/2 = Ω ρn+1i+1/2,rn+1i,rn+1i+1,vn+1i,vn+1

(63)

Президиума РАН №28 “Космос: исследования

i+1

фундаментальных процессов и их взаимосвязей”.

где обозначено

Приложение

Πn+σi+1/2 = σΠn+1i+1/2 + (1 - σ)Πni+1/2.

Параметр σ характеризует степень неявности схе-

РАЗНОСТНАЯ СХЕМА

мы и изменяется в пределах от 0 до 1. Этим

Зададим в расчетной области разностную сетку

параметром определяется порядок аппроксимации

q_i, i = 0,1,... ,N, координаты узлов которой свя-

по времени. В случае σ = 1/2 схема имеет вто-

заны с массовой лагранжевой координатой q. Шаг

рой порядок аппроксимации по Δt, а в остальных

сетки Δqi+1/2 = qi+1 - q_i выберем таким образом,

случаях — первый. Величина ω описывает искус-

чтобы в логарифмическом масштабе расстояния

ственную вязкость, необходимую для более кор-

между соседними узлами были одинаковыми, т.е.

ректного описания решений с ударными волнами

ln(qi+1/q_i) = const. Радиальную координату и ско-

и контактными разрывами. Явный вид функции Ω

рость в моменты времени tⁿ будем относить к узлам

определяется конкретной моделью искусственной

сетки: rⁿⁱ, v^nr,i, v^nϕ,i, а термодинамические величины

вязкости. В наших расчетах использовалась ли-

будем относить к центрам ячеек, которые нумеру-

нейная вязкость вида



)

ются полуцелыми индексами: ρni+1/2, Pni+1/2, εni+1/2,

∂vr

ηρ

⁽∂v

ω=-



(64)

Tni+1/2 [29]. Для упрощения дальнейших записей



∂r

 ∂r

будем использовать следующие обозначения для

где коэффициент ηi+1/2 = η₀Δqi+1/2. Значение η₀

операторов конечных разностей:

задавалась в диапазоне от 3 до 5.

fn+1 - f

(Δ_tf)n+1/2 =

(54)

Эта схема является полностью консервативной

Δt

в том смысле, что она обеспечивает не только

fi+1 - f_i

баланс полной энергии, но и отдельных ее видов

(Δ_qf)i+1/2 =

Δqi+1/2

(тепловой, кинетической, вращательной и грави-

тационной). Кроме того, нетрудно показать, что в

fi+1/2 - fi-1/2

(Δ_qf)_i =

данной схеме выполняется разностное соотноше-

Δq_i

ние

где Δq_i = (Δqi+1/2 + Δqi-1/2)/2. При этом всюду,

Δ_t (r_iv_ϕ,i) = 0,

(65)

где это не вызывает недоразумений, индексы у

описывающее закон сохранения углового момента.

операторов будем опускать.

Это означает, что имеет место равенство

Уравнения (23), (35), (22), (21), (36) можно

аппроксимировать следующей разностной схемой:

rⁿⁱv^nϕ,i = r0iv0ϕ,i = l_i,

(66)

Δ_q(rn+1)²,

(55)

где l_i — сеточная функция для удельного углового

ρn+1i+1/2

момента. Отсюда следует, что

Δ_tv_r,i = -rn+1/2iΔ_qΠn+σ +

(56)

l_i

v^nϕ,i =

(67)

n+1/2

rⁿ

(v_ϕ

)²

(1 - ξ)μ²

rn+1/2i

rn+1irⁿⁱ

Система алгебраических уравнений (55)-(63)

является нелинейной. Для ее решения использует-

Δ_tr_i = vn+1/2r,i,

(57)

ся комбинация метода Ньютона и метода прогонки.

(

)

Δ_tεi+1/2 = -Πn+σi+1/2Δ_q rn+1/2vn+1/2

(58)

n+1/2

СПИСОК ЛИТЕРАТУРЫ

v_r

vn+1/2ϕ,i

Δ_tv_ϕ,i = -

(59)

1. B. P. Abbott, R. Abbott, T. D. Abbott,

M. R. Abernathy, F. Acernese, K. Ackley, C. Adams,

rn+1/2

T. Adams, P. Addesso, R. X. Adhikari, et al.,

Pn+1

=ρn+1

Tn+1

(60)

Astrophys. J. Lett. 818, L22 (2016).

i+1/2

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

БИСИКАЛО и др.

B. P. Abbott, R. Abbott, T. D. Abbott,

20.

V. M. Lipunov, G. B ¨orner, and R. S. Wadhwa,

M. R. Abernathy, F. Acernese, K. Ackley, C. Adams,

Astrophysics of Neutron Stars (Springer, 1992).

T. Adams, P. Addesso, R. X. Adhikari, et al., Physical

21.

J. M. Bardeen, Nature 226, 64 (1970).

Review Lett. 116(24), 241103 (2016).

22.

K. S. Thorne, Astrophys. J. 191, 507 (1974).

B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese,

23.

L.-X. Li and B. Paczy ´nski, Astrophys. J. Lett. 534,

K. Ackley, C. Adams, T. Adams, P. Addesso,

L197 (2000).

R. X. Adhikari, V. B. Adya, et al., Physical Review

24.

D. B. Bowen, V. Mewes, M. Campanelli, S. C. Noble,

Lett. 118(22), 221101 (2017).

J. H. Krolik, and M. Zilhгo, Astrophys. J. Lett. 853,

B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese,

L17 (2018).

K. Ackley, C. Adams, T. Adams, P. Addesso,

25.

P. Artymowicz and S. H. Lubow, Astrophys. J. 421,

R. X. Adhikari, V. B. Adya, et al., Astrophys. J. Lett.

651 (1994).

851, L35 (2017).

26.

R. Dong, E. Vorobyov, Y. Pavlyuchenkov, E. Chiang,

B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese,

and H. B. Liu, Astrophys. J. 823, 141 (2016).

K. Ackley, C. Adams, T. Adams, P. Addesso,

R. X. Adhikari, V. B. Adya, et al., Physical Review

27.

P. J. Armitage, e-Print arXiv:astro-ph/0701485

Lett. 119(14), 141101 (2017).

(2007).

A. V. Tutukov and A. M. Cherepashchuk, Astronomy

28.

E. P. Kurbatov and D. V. Bisikalo, Astronomy Reports

Reports 61, 833 (2017).

61, 475 (2017).

P. V. Kaigorodov, D. V. Bisikalo, A. M. Fateeva, and

29.

A. A. Samarskii and I. P. Popov, Difference methods

A. Y. Sytov, Astronomy Reports 54, 1078 (2010).

for solving problems of gas dynamics (2nd revised

A. Y. Sytov, P. V. Kaigorodov, A. M. Fateeva, and

and enlarged edition) (Moscow, Izdatel Nauka,

D. V. Bisikalo, Astronomy Reports 55, 793 (2011).

1980).

N. I. Shakura and R. A. Sunyaev, Astron. and

30.

V. L. Ginzburg, Theoretical physics and

Astrophys. 24, 337 (1973).

astrophysics (Pergamon, Oxford, 1979).

10.

N. Bode and S. Phinney, APS April Meeting

31.

A. F. Illarionov and R. A. Syunyaev, Soviet Ast. 16, 45

Abstracts, S1.010 (2007).

(1972).

11.

J. D. Bekenstein, Astrophys. J. 183, 657 (1973).

32.

G. G. Pavlov, Y. A. Shibanov, and P. M ´esz ´aros, Phys.

12.

B. Kocsis and A. Loeb, Physical Review Lett. 101(4),

Rep. 182, 187 (1989).

041101 (2008).

33.

H. Feng and P. Kaaret, Astrophys. J. Lett. 650, L75

13.

M. Megevand, M. Anderson, J. Frank,

(2006).

E. W. Hirschmann, L. Lehner, S. L. Liebling,

34.

S. Fabrika, ASP Conf. Ser. 510, 395 (2017).

P. M. Motl, and D. Neilsen, Physical Review D

35.

XMM-Newton Community Support Team,

80(2), 024012 (2009).

XMM-Newton Users Handbook. Issue 2.16 URL:

14.

S. M. O’Neill, M. C. Miller, T. Bogdanovi ´c,

https://xmm-tools.cosmos.esa.int/external/

C. S. Reynolds, and J. D. Schnittman, Astrophys. J.

xmm_user_support/documentation/uhb/XMM_

700, 859 (2009).

UHB.html.

15.

L. R. Corrales, Z. Haiman, and A. MacFadyen,

36.

L. Maraschi and S. Molendi, Astrophys. J. 353, 452

Monthly Not. Roy. Astron. Soc. 404, 947 (2010).

(1990).

16.

G. P. Rosotti, G. Lodato, and D. J. Price, Monthly

37.

R. L. C. Starling, K. Wiersema, A. J. Levan,

Not. Roy. Astron. Soc. 425, 1958 (2012).

T. Sakamoto, et al., Monthly Not. Roy. Astron. Soc.

17.

M. J. Fitchett, Monthly Not. Roy. Astron. Soc. 203,

411, 2792 (2011).

1049 (1983).

38.

P. Ubertini, IBIS: Imager on Board the INTEGRAL

18.

H. Pietil ¨a, P. Hein ¨am ¨aki, S. Mikkola, and

Satellite URL: https://www.cosmos.esa.int/

M. J. Valtonen, Celestial Mechanics and Dynamical

Astronomy 62, 377 (1995).

web/integral/instruments-ibis.

19.

S. E. de Mink and A. King, Astrophys. J. 839, 7

39.

P. Ubertini, F. Lebrun, G. Di Cocco, A. Bazzano, et

(2017).

al., Astron. and Astrophys. 411, L131 (2003).

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019