Астрономический журнал, 2021, T. 98, № 11, стр. 901-921

1. ВВЕДЕНИЕ

С 2015 г., когда были экспериментально обнаружены гравитационные волны [1], начинается новая страница в астрономии. Возможность использовать “мультиканальный” (multi-messenger) подход позволяет достичь более глубокого понимания природы астрофизических объектов. Согласно теоретическим представлениям [2, 3], наибольшей интенсивностью гравитационно-волнового излучения обладают двойные черные дыры [4], а наибольшую вероятность обнаружения имеют двойные черные дыры звездных масс [5, 6]. И действительно, начиная с 2015 г. на детекторах LIGO и Virgo было обнаружено несколько гравитационно-волновых событий, сопровождающих слияние пар черных дыр [1, 7–11]. Возникает вопрос: а как нам добиться реализации “мультиканального” подхода для этих событий, или, другими словами, возникает ли электромагнитный отклик при слиянии черных дыр и можно ли его наблюдать?

В исследованиях черных дыр часто делается предположение о том, что эти объекты окружены аккреционными дисками. Сценарии формирования диска вокруг двойной сверхмассивной черной дыры, а также вокруг двойной черной дыры звездной массы указывались нами в работах [12, 13]. Рассмотрим задачу про возмущение аккреционного диска или оболочки вокруг сливающихся черных дыр. Согласно ОТО потеря энергии двойной черной дыры при излучении гравитационных волн приводит к формированию конечного объекта с массой на несколько процентов меньше исходной [3]. Кроме того, потеря импульса из-за отклонения от симметрии может привести к тому, что объект получает импульс отдачи (kick) [14] и в результате этого он приобретает скорость до 1000 км/с [15]. Вещество, окружающее сливающиеся черные дыры, испытывает возмущение [16], которое может привести к увеличению светимости и квазипериодическим вспышкам. Эффект отдачи черной дыры и последующее за этим гравитационное возмущение могут возбуждать в аккреционном диске ударные волны [17], а сами возмущения метрики вызывают в диске механические напряжения, которые могут диссипировать на вязкой шкале времени [18] или превращаться в ударные волны [19]. Нагрев вещества диска ударными волнами должен приводить к резкому увеличению яркости диска, что можно расценивать как электромагнитный отклик системы на эффект слияния двойной черной дыры.

Для случая слияния сверхмассивных черных дыр в ядрах галактик с полной массой ${{10}^{6}}\,{{M}_{ \odot }}$ потеря массы в 5% приводит к росту светимости аккреционного диска до ${{10}^{{43}}}$ эрг/с [20]. Эффект импульса отдачи может приводить к такому же росту светимости, однако его величина существенным образом зависит от плохо определенных параметров сливающейся двойной черной дыры [21, 22]. Отклик примерно такой же величины можно ожидать и для черных дыр звездных масс. В работе [23] было показано, что рост светимости аккреционного диска может достигать ${{10}^{{42}}}{\kern 1pt} - {\kern 1pt} {{10}^{{43}}}$ эрг/с, при этом максимум энергии излучения будет расположен в рентгеновской части спектра даже для случая волн плотности вместо ударных. В наших работах [12, 24] для моделей, параметры которых соответствовали событию GW170814 [10], были рассчитаны болометрические кривые блеска, спектры электромагнитного излучения и получены оценки длительности вспышки. Как оказалось, основная часть энергии электромагнитного излучения высвечивается в рентгеновском и гамма-диапазонах, и такого рода вспышки могут быть зарегистрированы инструментами, существующими в настоящее время.

Анализ результатов численных расчетов, проведенных в работах [12, 24], позволил заключить, что при определенном характере распределения плотности и температуры в исходном невозмущенном аккреционном диске режим распространения ударных волн в возмущенном диске с течением времени становится автомодельным. В  частности, расстояние, пройденное фронтом ударной волны за время $t$, описывается степенным законом ${{t}^{{2/3}}}$, а ее скорость затухает со временем как ${{t}^{{ - 1/3}}}$. В работе [13] мы получили полное автомодельное решение этой задачи с использованием необходимой информации об ударных волнах, взятой из численных расчетов. Из автомодельного решения можно получить алгебраические интегралы углового момента и адиабатичности, а также найти асимптотические соотношения для автомодельных функций. Такая аналитическая модель на основе автомодельного решения позволяет приближенно оценивать величину электромагнитного отклика гравитационно-волнового события без проведения трудоемких численных расчетов. В данной работе мы уточнили разработанную ранее модель путем учета вертикального расширения диска, приводящие к адиабатическому охлаждению вещества за ударными волнами. Использованная методика позволяет приблизить модель к реальной физической картине и получить достоверные оценки электромагнитного отклика от сливающихся черных дыр звездных масс.

Статья организована следующим образом. В разделе 2 описаны постановка задачи и основные уравнения. В разделе 3 представлены автомодельные уравнения и проведен их анализ. В разделе 4 описаны результаты расчетов. В разделе 5 обсуждаются характеристики электромагнитного отклика. В разделе 6 содержатся основные выводы работы. В Приложении приведены некоторые детали численного метода и аналитического решения задачи.

2. ОСНОВНЫЕ УРАВНЕНИЯ

С помощью гравитационно-волновых детекторов LIGO и Virgo были проведены две завершенные серии наблюдений O1 (с 12 сентября 2015 г. по 19 января 2016 г.) и O2 (с 30 ноября 2016 г. по 26 августа 2017 г.). Третья серия наблюдений O3 началась 1 апреля 2019 г. и должна была завершиться 30 апреля 2020 г. Однако в конце марта 2020 г. наблюдения были прекращены из-за пандемии коронавируса. В завершенных сериях наблюдений O1 и O2 было обнаружено 10 событий слияния двойных черных дыр [11], из которых только 5 не имеют артефактов. Эти события представлены в табл. 1 вместе с оценками масс компонентов, доли теряемой массы и расстояния до соответствующего объекта. Как видно из таблицы, все эти объекты соответствуют черным дырам звездной массы. При этом доля теряемой массы при слиянии черных дыр составляет величину порядка 2–6%. В качестве исследуемого объекта в данной работе будем рассматривать событие GW170814, соответствующее слиянию двойной черной дыры с массами компонентов 30.5 и 25.3 ${{M}_{ \odot }}$. Для определенности будем считать, что потеря массы в этом событии составляла 5%.

При исследовании черных дыр часто предполагается, что они окружены аккреционными дисками. Внутренний радиус аккреционного диска вокруг одиночной черной дыры равен $3{{R}_{{\text{g}}}}$, где ${{R}_{{\text{g}}}} = 2GM{\text{/}}{{c}^{2}}$ – гравитацонный радиус, $G$ – гравитационная постоянная, $c$ – скорость света. Внутри этой зоны устойчивых круговых орбит не существует [3]. Предположим, что двойная черная дыра также окружена общим аккреционным диском. Конечная стадия слияния черных дыр происходит очень быстро, за доли секунды. Процесс стремительного слияния черных дыр начинается, когда расстояние между ними составляет $A = 3{{R}_{{{\text{g}},1}}} + 3{{R}_{{{\text{g}},2}}} = 3{{R}_{{\text{g}}}}$, где ${{R}_{{{\text{g}},1}}}$ и ${{R}_{{{\text{g}},2}}}$ – гравитационные радиусы компонентов, ${{R}_{{\text{g}}}}$ – гравитационный радиус, вычисленный по их суммарной массе $M = {{M}_{1}} + {{M}_{2}}$ (гравитационный радиус аккретора) [25–28]. Внутренний радиус аккреционного диска вокруг двойной черной дыры примерно в два раза превышает расстояние между центрами черных дыр. Такой вывод следует из оценок устойчивости орбит относительно приливных сил [12, 13] и хорошо согласуется с результатами численного моделирования двойных систем с круговыми орбитами компонентов, отношение масс которых близко к единице [29–32].

До момента слияния сближение черных дыр происходит медленно и поэтому, благодаря действию вязких сил, аккреционный диск успевает подстроиться к изменяющемуся межкомпонентному расстоянию. Следовательно, можно считать, что непосредственно перед слиянием черных дыр внутренний радиус диска ${{r}_{{\text{d}}}}$ примерно в 6 раз превосходит гравитационный радиус аккретора ${{R}_{{\text{g}}}}$. Это позволяет при описании свойств такого аккреционного диска пренебрегать релятивистскими эффектами и использовать теорию Шакуры и Сюняева [33]. В этой модели скорость аккреции $\dot {M}$ определяется темпом диссипации углового момента, возникающего вследствие турбулентной вязкости. Интенсивность турбулентности описывается безразмерным параметром $\alpha $, значение которого не превышает единицы. Критический темп аккреции определяется величиной эддингтоновской светимости

где $m = M{\text{/}}{{M}_{ \odot }}$ – безразмерная масса центрального объекта. Параметр $\chi $ определяет эффективность аккреции и в случае одиночной невращающейся черной дыры равен примерно 0.06–0.08 [34]. Учет вращения повышает эффективность стационарной аккреции до 0.32 [35, 36]. В сценарии эпизодической аккреции эффективность может достигать 0.43 [37]. Однако в нашей модели конкретное значение этого параметра не является критичным, поскольку определяющую роль играет величина безразмерного темпа аккреции $\dot {m} = \dot {M}{\text{/}}{{\dot {M}}_{{{\text{cr}}}}}$.

В модели [33] в аккреционном диске выделяются три характерные зоны (A, B и C), в которых распределение гидродинамических величин определяется относительным вкладом газового давления и давления излучения, а также механизмом непрозрачности. В зоне A преобладает давление излучения, а в зонах B и C им можно пренебречь по сравнению с газовым давлением. В зонах A и B непрозрачность определяется томсоновским рассеянием, а в зоне C непрозрачность обеспечивается свободно-свободными переходами. Вертикальная структура диска в этих зонах различается довольно сильно, но радиальные зависимости плотности и температуры оказываются близкими к степенным с почти одинаковыми показателями степени.

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

где параметр $\alpha $ определяет эффективность переноса углового момента в диске (параметр Шакуры–Сюняева).

Поэтому в нашей модели для задания начальных условий мы используем следующие степенные распределения плотности и температуры в невозмущенном диске:

где ${{\rho }_{*}}$ и ${{T}_{*}}$ – соответствующие значения в точке $r = {{r}_{*}}$. Как следует из выражений (2) и (3), для зоны B показатели ${{k}_{d}} = 33{\text{/}}20$, ${{k}_{t}} = 9{\text{/}}10$. Начальные профили остальных величин будут представлены ниже. Модель Шакуры–Сюняева [33] используется только для описания невозмущенного состояния аккреционного диска. Поскольку процесс релаксации диска после слияния черных дыр и потери массы центрального объекта происходит относительно быстро, то дальнейшую эволюцию диска можно моделировать в приближении бездиссипативной (без учета вязкости) гравитационной газовой динамики [12, 24].

Уравнения гравитационной газовой динамики, описывающие эволюцию геометрически тонкого аккреционного диска после слияния черных дыр, в осесимметричном случае в цилиндрических координатах ($r$, $\varphi $, $z$) можно записать в следующем виде:

Здесь $\rho $ – плотность, ${{{v}}_{r}}$ – радиальная скорость, ${{{v}}_{z}}$ – вертикальная скорость, ${{{v}}_{\varphi }}$ – азимутальная скорость, $P$ – давление, $\varepsilon $ – удельная внутренняя энергия, $M$ – полная масса двойной черной дыры до слияния, $\xi $ – доля массы черных дыр, излученная в виде гравитационных волн.

В случае геометрически тонких дисков мы можем пренебречь подробным описанием вертикальной структуры для исследования эффектов ударного нагрева и адиабатического расширения. Вместо этого можно использовать простое приближение, когда вертикальная скорость ${{{v}}_{z}}$ в любой момент времени остается пропорциональной высоте $z$ (аналог закона Хаббла в космологии):

В плоскости симметрии диска вертикальный компонент скорости ${{{v}}_{z}} = 0$. В нашей модели все уравнения (6)–(10) мы записываем для плоскости симметрии диска. В этом случае члены вида ${{{v}}_{z}}\partial {\text{/}}\partial z$ в конвективных производных исчезают во всех уравнениях. С учетом этого замечания и приближения (11) уравнения (6)–(10) (кроме уравнения (8)) примут вид: Последнее слагаемое в уравнении энергии (15) описывает адиабатический нагрев диска, обусловленный его вертикальными движениями.

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

где $H$ – локальная полутолщина диска. Тогда из уравнения для вертикальной скорости (8) с помощью соотношения (11) можно получить уравнения для коэффициента $\psi $, Подставляя $z = H(r,t)$ в соотношение (11), находим С другой стороны, мы можем написать Поэтому локальная полутолщина диска $H(r,t)$ удовлетворяет уравнению

Для замыкания полученной системы уравнений (12)–(15), (17), (20) следует учесть уравнения состояния, определяющие зависимости давления $P$ и внутренней энергии $\varepsilon $ от плотности $\rho $ и температуры $T$. Поскольку в данной работе мы рассматриваем только зону B аккреционного диска, то давление излучения мы не учитываем. Однако следует отметить, что в сильно возмущенных областях этой зоны, где возникают ударные волны большой интенсивности, температура может возрастать настолько, что эффекты давления излучения могут играть существенную роль [12]. Уравнения состояния можно записать в виде:

где ${{A}_{{{\text{gas}}}}} = 2{{k}_{{\text{B}}}}{\text{/}}{{m}_{{\text{p}}}}$ – газовая постоянная, ${{k}_{{\text{B}}}}$ – постоянная Больцмана, ${{m}_{{\text{p}}}}$ – масса протона, $\gamma $ – показатель адиабаты. Для полностью ионизованной водородной плазмы средний молекулярный вес следует принять равным $1{\text{/}}2$, а показатель адиабаты $\gamma = 5{\text{/}}3$.

В качестве начальных условий для плотности и температуры используются выражения (4) и (5). Начальные распределения для давления и внутренней энергии можно получить из уравнений состояния (21) и (22),

где квадрат скорости звука определяется как $c_{*}^{2} = {{A}_{{{\text{gas}}}}}{{T}_{*}}$. Кроме того, считаем, что в начальный момент ${{v}_{r}}(r,0) = 0$, $\psi (r,0) = 0$. Начальные условия для азимутальной скорости и полутолщины диска можно получить из стационарных решений уравнений (13) и (17), соответствующих невозмущенному диску, когда параметр $\xi = 0$: Безразмерный параметр определяет число Маха на радиусе ${{r}_{*}}$, где ${{t}_{*}} = {{r}_{*}}{\text{/}}{{c}_{*}}$ – характерный масштаб времени.

Для численного решения описанных выше уравнений использовалась неявная полностью консервативная разностная схема Самарского–Попова [38], некоторые детали которой описаны в [39] (см. Приложение A). Следует подчеркнуть, что в этой схеме выполняются не только разностные аналоги законов сохранения массы, импульса и энергии, но и дополнительные соотношения, описывающие баланс по определенным видам энергии. Кроме того, в нашем случае в схеме выполняется разностный аналог закона сохранения удельного углового момента. Численная модель определяется следующими параметрами: показателями ${{k}_{d}}$ и ${{k}_{t}}$ начальных профилей плотности и температуры, долей теряемой массы $\xi $ и числом Маха $\mu $. Отметим, что в частном случае $\psi = 0$ эта численная модель переходит в более простую модель, рассмотренную нами в работе [24], в которой не учитываются эффекты, обусловленные адиабатическим расширением диска.

Будем считать, что безразмерная масса $m = 55$ соответствует гравитационно-волновому событию GW170814 (см. табл. 1). Выберем в качестве характерного радиуса ${{r}_{*}}$, который определяет пространственный масштаб исследуемой задачи, внутренний радиус невозмущенного диска, ${{r}_{{\text{d}}}} = 6{{R}_{{\text{g}}}}$. Тогда с учетом распределения температуры в невозмущенном диске (5) из (27) находим

Если зафиксировать здесь значение параметра $\dot {m}$, то в пределе при $\alpha \to 0$ мы можем получить сколь угодно малые значения параметра $\mu $. Наоборот, если мы зафиксируем значения параметра $\alpha $, то в пределе при $\dot {m} \to 0$ можно получить сколь угодно большие значения параметра $\mu $. Например, в случае $\alpha = 0.01$, $\dot {m} = 1$ получаем значение $\mu = 49.23$. В случае $\alpha = 0.01$, $\dot {m} = 0.01$ получаем значение $\mu = 123.65$. Отсюда видно, что возможные значения параметра $\mu $ лежат в довольно широких пределах.

Поскольку мы рассматриваем только зону B аккреционного диска, то характерный радиус ${{r}_{*}}$ можно определить и как внутренний радиус этой зоны ${{r}_{{{\text{AB}}}}}$. Величину этого радиуса можно найти из условия равенства давления излучения и газового давления [33],

Используя это выражение, снова из (5) и (27) находим Эта зависимость $\mu $ от $\alpha $ и $\dot {m}$ имеет такой же качественный характер, как и рассмотренная выше зависимость (28). Параметр $\mu $ становится сколь угодно малым в пределе при $\alpha \to 0$ и сколь угодно большим в пределе при $\dot {m} \to 0$. В частности, в случае $\alpha = 0.01$, $\dot {m} = 1$ получаем значение $\mu = 41.90$, а в случае $\alpha = 0.01$, $\dot {m} = 0.01$ параметр $\mu = 194.49$. Эти значения очень близки к тем, что мы получили в предыдущем случае. На основе полученных оценок в приведенных ниже расчетах мы варьировали значения безразмерного параметра $\mu $ в пределах от 50 до 200.

Как показывают наши расчеты [12, 13, 24], степень возмущения аккреционного диска в результате слияния черных дыр и потери массы за счет излучения гравитационных волн определяется величиной параметра $\mu $. Электромагнитный отклик от гравитационно-волнового события оказывается тем выше, чем больше значение $\mu $. В  связи с этим, следует заметить, что в случае сверхмассивных черных дыр значения параметра $\mu $ могут быть намного больше по сравнению с приведенными выше оценками для черных дыр звездной массы. В самом деле, выберем в качестве характерного радиуса ${{r}_{ * }}$ внутренний радиус невозмущенного диска. Тогда, задавая безразмерную массу, например, равной $m = {{10}^{8}}$, получаем

Отсюда в случае $\alpha = 0.01$, $\dot {m} = 1$ получаем значение $\mu = 208.05$, а в случае $\alpha = 0.01$, $\dot {m} = 0.01$ параметр $\mu = 522.59$. Для более массивных черных дыр можно получить еще большие значения параметра $\mu $.

3. АВТОМОДЕЛЬНОЕ РЕШЕНИЕ

4. РЕЗУЛЬТАТЫ РАСЧЕТОВ

В предыдущем разделе 3 мы показали, что автомодельное решение исследуемой задачи существует. При этом мы подробно исследовали некоторые его свойства. Результаты численных расчетов полной задачи позволяют утверждать, то с течением времени решение действительно сходится к автомодельному. Сходимость численного решения уравнений (33)–(36) к автомодельному решению демонстрирует рис. 1. На рисунке показаны численные зависимости $\sigma (\lambda )$ (слева) и $u(\lambda )$ (справа) для различных моментов времени $t{\text{/}}{{t}_{*}}$. Решения соответствуют значению параметра $\mu = 5$. В результате слияния черных дыр и потери массы аккреционный диск возмущается и в нем формируется ударная волна слабой интенсивности, распространение которой происходит по автомодельному закону $\lambda = 1.428$ с показателем автомодельности $\delta = 2{\text{/}}3$. При этом в области перед ударной волной решение становится автомодельным практически сразу. Расчеты для больших значений $\mu $ показывают, что выход на автомодельный режим происходит гораздо быстрее. Таким образом, наряду с численным решением имеет смысл рассматривать и автомодельное решение.

Мы провели серию численных расчетов возмущения аккреционного диска, возникающего в результате слияния черных дыр и потери массы. С учетом оценок, приведенных в конце раздела 2, мы задавали значения параметра $\mu $ в диапазоне $50 \leqslant \mu \leqslant 200$ (см. табл. 2).

В работе [13] мы рассматривали чисто одномерную модель возмущенного аккреционного диска и строили автомодельное решение с помощью следующей методики. Сначала для достаточно большого значения автомодельной переменной $\lambda = {{\lambda }_{{max}}}$ вычислялись значения автомодельных функций по формулам, которые описывают асимптотику решения в пределе при $\lambda \to \infty $. Эти значения использовались в качестве начальных условий. Далее, задавая некоторый шаг автомодельной переменной $\Delta \lambda < 0$, с помощью численного метода находились значения величин для $\lambda < {{\lambda }_{{max}}}$. Количество ударных волн и положения их фронтов ${{\lambda }_{s}}$, где индекс $s$ определяет номер ударной волны, из самого автомодельного решения простым способом найти не удается. Эту необходимую информацию об ударных волнах мы брали из соответствующего численного решения. Как только в процессе численного интегрирования автомодельных уравнений автомодельная переменная $\lambda $ становилась равной ${{\lambda }_{1}}$, значение функций в точке $\lambda + \Delta \lambda $ вычислялись с помощью условий Гюгонио (48)–(53). После этого процедура интегрирования продолжалась дальше до следующей ударной волны.

В данной работе мы рассматриваем модель, соответствующую 1.5-мерному приближению, в котором учитываются вертикальные движения вещества в аккреционном диске. Это приводит к значительно более сложному характеру решения. В результате с помощью описанной выше методики получить полное автомодельное решение не удается. Проблема возникает вследствие того, что в данном случае использованная нами численная схема Рунге–Кутты 4-го порядка становится неустойчивой в окрестности особых точек (см. Приложение C), которые располагаются сразу за ударными волнами. Для корректного решения этой проблемы необходимо использовать какой-либо более специализированный численный метод. Поскольку у нас и так уже есть автомодельное решение, полученное с помощью численного моделирования полной задачи, то мы не стали этого делать. Автомодельное решение, полученное с помощью численного расчета полной задачи, с учетом найденных асимптотик, алгебраических интегралов углового момента и адиабатичности, а также информации об ударных волнах позволяет в полной мере оценивать величину и спектральные характеристики электромагнитного отклика от соответствующего гравитационно-волнового события.

Полученные решения для всех моделей представлены на рис. 2–7. Эти рисунки показывают, что во всех моделях в возмущенном аккреционном диске формируются три ударные волны, следующие друг за другом. Автомодельные координаты ${{\lambda }_{s}}$ этих ударных волн указаны в табл. 2. В модели с максимальным параметром $\mu = 200$ на самом деле формируется еще одна четвертая ударная волна, которой отвечает автомодельная координата ${{\lambda }_{s}} = 3.622$, но она имеет очень слабую интенсивность. Область сильного возмущения, где расположены ударные волны, сосредоточена в сравнительно узком интервале автомодельной переменной $2 < \lambda < 10$. С ростом значения параметра $\mu $ энергия возмущения перераспределяется между ударными волнами, которые осуществляют разгрузку этого возмущения. Угловой момент $\Lambda (\lambda )$ (см. рис. 4), а также автомодельные функции $\Psi (\lambda )$ (см. рис. 6) и $h(\lambda )$ (см. рис. 7) остаются непрерывными на ударных волнах, как и следует из полученных нами условий Гюгонио (49) и (53). Остальные автомодельные функции на ударных волнах испытывают скачки.

Анализ рисунков позволяет заключить, что интенсивность первой (самой внешней) ударной волны постепенно падает, в то время как интенсивности внутренних ударных волн постепенно нарастают. Особенно это касается второй (промежуточной) ударной волны. Так, в случае μ = 200 основная энергия возмущения сосредоточена во второй ударной волне. Интересно отметить, что перед каждой ударной волной располагаются нелинейные предвестники, распространяющиеся впереди них. Особенно их хорошо видно на распределении $w(\lambda )$ (см. рис. 5). При этом перед самой внешней ударной волной имеются сразу два таких предвестника, а предвестник перед самой внутренней ударной волной имеет слабую интенсивность.

Во внешней области перед ударными волнами автомодельное решение является гладким и при больших $\lambda $ описывается асимптотическими выражениями, полученными в разделе 3.4. Во внутренней области за ударными волнами при малых $\lambda $ автомодельное решение также переходит в соответствующую асимптотику при $\lambda \to 0$. По полученным результатам расчета непосредственной проверкой можно убедиться, что алгебраический интеграл углового момента (57) удовлетворяется во всей области изменения автомодельной переменной. В то же время алгебраический интеграл адиабатичности (61) строго совпадает с полученным решением только во внешней области до первой ударной волны.

В области за ударными волнами структура аккреционного диска становится существенно нестационарной. Распределение радиальной скорости $u(\lambda )$ имеет колебательный характер и она много раз меняет знак (см. рис. 3). Это означает, что одни части диска движутся наружу (${{v}_{r}} > 0$), а другие, наоборот, внутрь (${{v}_{r}} < 0$). Сильные вариации испытывает и вертикальная структура диска. Автомодельная функция $\Psi (\lambda )$, которая определяет вертикальную скорость (см. уравнения (11) и (37)), быстро изменяется во времени с явно выраженным колебательным характером (см. рис. 6). Это приводит к квазипериодическим изменениям полутолщины диска $h(\lambda )$ (см. рис. 7), амплитуда которых постепенно затухает.

После прохождения нелинейных и ударных волн в диске формируется новое стационарное состояние, в котором радиальная и вертикальная скорости равны нулю. Это новое состояние характеризуется более низкой плотностью, более высокой температурой и увеличенной толщиной диска. Анализ полученного решения (см. табл. 2) позволяет заключить, что в результате вертикального расширения диска плотность в плоскости симметрии во всех моделях (в модели с параметром $\mu = 50$ в $1.53$ раза) падает более, чем в 2 раза по сравнению с невозмущенным состоянием.

5. ЭЛЕКТРОМАГНИТНЫЙ ОТКЛИК

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

Болометрическую светимость аккреционного диска (если смотреть на него с одной стороны) можно рассчитать с помощью выражения

где ${{\sigma }_{{{\text{SB}}}}}$ – постоянная Стефана-Больцмана, ${{T}_{{{\text{eff}}}}}$ – локальная (на данном радиусе $r$) эффективная температура. Эффективная температура ${{T}_{{{\text{eff}}}}}$ не равна температуре в плоскости симметрии $T$, поскольку аккреционный диск в зоне B является оптически толстым и поэтому излучение идет от поверхностного слоя. При этом можно считать, что излучение диска в зоне B имеет планковский спектр. Температура вещества в этом поверхностном слое определяется процессами переноса тепловой энергии в вертикальном направлении от плоскости симметрии диска к его поверхности.

В геометрически тонком аккреционном диске для оценки эффективной температуры можно использовать простое приближение, в котором считается, что ${{T}_{{{\text{eff}}}}} = fT$, где коэффициент $f$ определяется конкретным процессом переноса тепловой энергии в вертикальном направлении [12, 24]. Если перенос тепловой энергии осуществляется путем лучистой теплопроводности, то можно принять, что коэффициент $f = 0.1$. Если же перенос энергии в вертикальном направлении определяется конвекцией, то перепад температуры можно оценить величиной $f = 0.5$.

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

Следовательно, болометрическая светимость (85) может быть представлена в виде где мы положили ${{r}_{{\text{d}}}} = {{r}_{*}}$. Интеграл по автомодельной переменной $\lambda $ можно разбить на сумму двух интегралов по интервалам ${{({{t}_{*}}{\text{/}}t)}^{{2/3}}} \leqslant \lambda \leqslant 1$ и $1 \leqslant \lambda \leqslant \infty $. В первом интеграле функцию $w(\lambda )$ с хорошей точностью можно заменить ее предельным при $\lambda \to 0$ значением ${{w}_{0}}$: Второй интеграл не зависит от времени и может быть рассчитан на основе полученного автомодельного решения. Его конкретное значение определяется параметром $\mu $ (см. последний столбец табл. 2). В результате получаем следующее выражение Эта функция монотонно возрастает со временем и в пределе $t \gg {{t}_{*}}$ стремится к постоянному значению

Полученные зависимости $L(t)$, нормированные на предельное значение ${{L}_{{max}}}$, для моделей с различными параметрами $\mu $ показаны на рис. 8. Анализ этих кривых блеска позволяет заключить, что светимость возмущенного диска за время порядка нескольких характерных времен ${{t}_{ * }}$ (это составляет несколько секунд) достигает постоянного значения ${{L}_{{max}}}$. Рост светимости определяется фактором $w_{0}^{4}$, значения которого для различных моделей приведены в табл. 2. С увеличением параметра $\mu $ интенсивность электромагнитного отклика существенно возрастает. Если в модели светимость возрастает всего на порядок, то в модели $\mu = 200$ рост светимости составляет уже более 4 порядков величины.

В дальнейшем нагретое вещество будет постепенно остывать за счет радиативного охлаждения и светимость возмущенного диска будет медленно уменьшаться. В гидродинамической модели и в автомодельном решении этот эффект мы не учитывали. Однако характерное время охлаждения ${{t}_{{{\text{cool}}}}}$ нетрудно оценить [13]. Объемный коэффициент охлаждения в диффузионном приближении переноса излучения может быть записан в виде

где $\kappa $ – коэффициент непрозрачности, который в зоне B определяется томсоновским рассеянием. Поэтому характерное время охлаждения Из полученных результатов следует, что в точке $r = {{r}_{*}}$ характерное время охлаждения ${{t}_{{{\text{cool}}}}}$ составляет несколько единиц ${{t}_{*}}$ и в этой области охлаждение происходит достаточно быстро. Однако с увеличением расстояния время охлаждения ${{t}_{{{\text{cool}}}}}$ становится бóльшим. Например, на расстоянии $r = 20{{r}_{*}}$ время охлаждения в несколько сотен раз превышает характерный временнóй масштаб задачи ${{t}_{*}}$.

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

определяется функцией Планка где $\nu $ – частота, $h$ – постоянная Планка. Для расчета спектра электромагнитного излучения необходимо задавать параметры аккреционного диска, а также значения величин на его внутренней границе. Поэтому мы положили ${{r}_{{\text{d}}}} = {{r}_{*}}$ и для определенности задали значение параметра Шакуры–Сюняева $\alpha = 0.01$. Тогда температуру ${{T}_{*}}$ на внутренней границе диска можно вычислять по формуле (3), где безразмерный темп аккреции $\dot {m}$ выражается через параметр $\mu $ из соотношения (28).

Полученные спектры электромагнитного излучения, соответствующие автомодельному решению, представлены на рис. 9. Спектры для радиативных дисков ($f = 0.1$) для различных моделей показаны на левой панели рис. 9, а для конвективных дисков ($f = 0.5$) на правой панели. При этом жирные линии описывают спектры возмущенных дисков, а тонкие линии соответствуют спектрам невозмущенных дисков. Значения потоков ${{F}_{\nu }}$ приведены в абсолютных величинах. Закрашенная вертикальная полоса указывает оптический диапазон спектра.

Как видно из рисунка, основная энергия излучения сосредоточена в диапазоне от 1 до 100 кэВ, что соответствует мягкому и жесткому рентгену. Кроме того, значительная доля энергии излучения приходит и в гамма-диапазоне ($ \geqslant {\kern 1pt} 100$ кэВ). В обоих случаях (радиативные и конвективные диски) максимальный поток излучения в процессе вспышки возрастает, как правило, на два порядка величины. При этом частота, на которой достигается максимум, смещается в область высоких частот. В случае конвективных дисков поток излучения на два порядка величины превышает поток от радиативных дисков. Кроме того, максимальное значение потока от конвективных дисков по сравнению с радиативными дисками смещено в область жесткого рентгеновского излучения. В видимой области спектра поток излучения возрастает примерно на порядок. При этом наибольший относительный эффект вспышки проявляется для наибольшего значения параметра $\mu = 200$. Поэтому следует ожидать, что в случае слияния сверхмассивных черных дыр поток излучения от возмущенного аккреционного диска должен возрастать в существенно большей степени, в том числе и в оптическом диапазоне спектра.

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

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

В работе исследовано автомодельное решение, описывающее эволюцию возмущенного аккреционного диска после слияния черных дыр и потери массы. При этом показатель автомодельности оказывается равным 2/3. Решение описывает структуру возмущенного диска в области (так называемая зона B), в которой доминирует газовое давление, а непрозрачность вещества определяется процессами томсоновского рассеяния электронов. В отличие от нашей предыдущей работы [13] в данной статье мы учли эффект, обусловленный адиабатическим охлаждением диска. Для построения автомодельного решения мы использовали результаты численного моделирования. Из автомодельных уравнений получены алгебраические интегралы углового момента и адиабатичности, выражающие соответствующие законы сохранения. Алгебраический интеграл углового момента удовлетворяется во всей области изменения автомодельной переменной. Алгебраический интеграл адиабатичности удовлетворяется только в областях гладкости решения, поскольку на ударных волнах энтропия испытывает скачок. Мы также нашли асимптотические соотношения для автомодельных функций, описывающие как развитие начального возмущения, так и переход диска к новому стационарному состоянию.

С помощью полученного автомодельного решения задачи можно приближенно оценивать величину электромагнитного отклика гравитационно-волнового события без проведения трудоемких численных расчетов. Анализ полученных автомодельных решений показывает, что для параметров, соответствующих событию GW170814, рост светимости диска может составить от одного до четырех порядков величины. В видимой области спектра поток излучения возрастает примерно на порядок. Этого может оказаться вполне достаточно для регистрации этого отклика современными наблюдательными инструментами.

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