Electromagnetic emission from circumbinary disk of merging black holes

Abstract

In the paper a scenario of an electromagnetic response formation from the merging of two black holes is considered. In this scenario it’s assumed that the binary black hole is surrounded by an accretion disk. As a result of the black holes merging and mass loss, the accretion disk experiences a disturbance, which is accompanied by shock waves propagation of sufficiently high intensity. Heating of matter by shock waves leads to a sharp increase in the flux of electromagnetic radiation from the disk. The paper includes a calculated light curve, radiation spectrum, and the estimation of characteristic duration of the flare. This method can be used to discover of electromagnetic responses from gravitational-wave events, which registered by the LIGO (Laser Interferometer Gravitational-Wave Observatory) and the Virgo detectors. Supporting the registration of gravitational-wave events by observations in the electromagnetic channel has far-reaching prospects, since it corresponds to the multi-messenger approach to the study of astrophysical objects.

Key words
gravitational waves; black holes; binary stars; accretion disks

Introduction

In classical observational astronomy, two important “channels“ or “windows“ were used to study objects. The first channel is associated with detection of electromagnetic radiation from space objects in the entire spectral range, from radio to gamma emission. The second channel is determined by various kinds of particles coming to Earth from outer space. They include, in particular, cosmic rays and neutrinos. In 2015 gravitational waves (Abbott et al. 2016a) were experimentally detected. As a result of this discovery, gravitational wave astronomy has become another crucial channel for exploring the Universe. Observational data, obtained from objects in different channels, allow us to get a much more of useful information about the source. After all, in this case, we simultaneously advance research from the point of view of fundamentally different physical processes! This is not simply a multi-frequency method of investigation, when the object is simultaneously observed in different ranges of the spectrum. This is the so-called “multi-messenger“ approach. Using such complex research methods, we can achieve a much deeper understanding of the nature of astrophysical objects. Obviously, this type of research has great and far-reaching prospects.

According to the theoretical concepts (Landau & Lifshitz 1975, Misner et al. 1973), binary black holes (Clark & Eardley 1977) have the highest intensity of gravitational-wave radiation, and binary black holes of stellar mass (Lipunov et al. 1997, Tutukov & Yungelson 1993) have the highest probability of detection. Since 2015 LIGO and Virgo detectors have discovered several gravitational-wave events accompanying the merging of pairs of black holes (Abbott et al. 2016a, b, 2017a, b, c, 2019). All these objects were located at distances from several hundred to several thousand Mps.

From the point of view of multi-messenger astronomy, the following important problem arises, having both a theoretical and an observational aspect. At the disposal of researchers, there is a means of detecting gravitational waves, the source of which are merging black holes. Is it possible to see the electromagnetic response accompanying this event? What energy and what spectral characteristics should it have? What is its characteristic duration in time?

It follows from the General theory of relativity that as a result of energy loss of a binary black hole, when gravitational waves are emitted, the mass of the merged product turns out to be several percent less than the original mass (Misner et al. 1973). In addition, the loss of momentum due to a deviation from symmetry can cause the object to receive a recoil pulse — “kick“ (Peres 1962) and as a result it gains speed of up to 1000 km/s (Bekenstein 1973). When a binary system is surrounded by an envelope or accretion disk, this surrounding matter experiences a perturbation (Bode & Phinney 2007) as a result of the merging, which can lead to increased luminosity and quasi-periodic flares. The black hole recoil effect and the subsequent gravitational perturbation can excite shock waves in the accretion disk (Megevand et al. 2009), and the metric perturbations themselves cause mechanical stresses in the disk, which can dissipate on the viscous time scale (Kocsis & Loeb 2008) or turn into shock waves (Cherepashchuk 2016).

Attempts to study these questions for the case of merging supermassive black holes in the galactic nuclei have been made repeatedly by many authors (see, e.g., Megevand et al., 2009, O’Neill et al. 2009, Corrales et al. 2010, Rosotti et al. 2012). In particular, the works (Lippai et al. 2008, MacFadyen & Milosavljević 2008, Milosavljević & Phinney 2005,Shields & Bonning 2008, Schnittman & Krolik 2008) studied perturbations of both geodesic orbits of particles and accretion disk in the framework of simple analytical models. The application to disks in Active Galactic Nuclei is considered in (McKernan et al. 2019). If a binary black hole with total mass of ${10}^6{M}_{\odot }$ after the merging loses 5% of its mass, the luminosity of the accretion disk increases by an order of magnitude, reaching the value of ${10}^{43}$ erg/s (Corrales et al. 2010). The effect of the recoil pulse can lead to the same increase in luminosity, but its value depends significantly on poorly defined parameters of the merging binary black hole (Fitchett 1983, Pietilä et al. 1995). For stellar mass black holes, estimates of accretion disk luminosity growth were obtained up to ${10}^{42}$–${10}^{43}$ erg/s with the maximum radiation energy in the X-ray range (de Mink & King 2017), however, it was assumed that the resulting waves are not transformed into shock waves.

Let’s consider the problem of the accretion disk perturbation around merging black holes. The disturbance occurs as a result of a decrease in the mass of the central object and it can cause shock waves in the disk, the propagation of which heats disk matter. This should lead to a sharp increase in brightness of the disk, which can be regarded as the electromagnetic response of the system to the effect of merging a binary black hole. An important advantage of this approach is that all the necessary parameters of the model, and, consequently, the parameters of the electromagnetic response, can be found directly from observations.

The electromagnetic response of the accretion disk around a stellar-mass binary black hole resulting from decrease in the mass of the central object after the merging was investigated in our papers (Bisikalo et al. 2019, Bisikalo et al. 2019). In (Bisikalo et al. 2019), a 1D numerical model was developed and series of calculations were performed for disks with dominant gas pressure. In (Bisikalo et al. 2019) we developed a more detailed 1.5D numerical model that takes into account the radiation pressure, as well as the adiabatic heating and cooling of the disk due to the vertical motion of matter in the disturbed disk. For models, with parameters corresponded to the event GW170814 (Abbott et al. 2017c), bolometric light curves, electromagnetic radiation spectra were calculated, and flare duration estimates were obtained. As it turned out, the main part of the electromagnetic radiation energy is displayed in the X-ray and gamma-ray ranges, where the luminosity reaches ${10}^{43}$–${10}^{45}$ erg/s. We have shown that such a flare, whose duration reaches hundreds of seconds, can be registered by tools that currently exist.

The analysis of the results of numerical calculations carried out in papers (Bisikalo et al. 2019, Bisikalo et al. 2019) allowed us to conclude that for a certain type of distribution density and temperature in the original undisturbed accretion disk, the propagation of shock waves in the disturbed disk becomes self-similar over time. In particular, the distance travelled by the shock wave front in time $t$ is described by the power law $t^{2/3}$, and its speed falls in time as $t^{-1/3}$. In (Bisikalo & Zhilkin 2020), we obtained a complete self-similar solution to this problem using the necessary information about shock waves taken from numerical calculations. From the self-similar solution, we can obtain algebraic integrals of angular momentum and adiabaticity, as well as find asymptotic relations for self-similar functions. This analytical model based on a self-similar solution allows to estimate approximately the value of electromagnetic response of a gravitational-wave event without performing time-consuming numerical calculations.

We emphasize that our numerical and analytical (self-similar) solutions assume spherical symmetry of the central source. In case of accounting recoil momentum it is violated and therefore our approach becomes inapplicable. It seems reasonable to solve such problems using three-dimensional numerical models. In this paper, we discuss the results of our work that are of interest to researchers who makes observations of astrophysical events (including black hole mergings) in the electromagnetic channel in various ranges of the spectrum.

The paper is organized as follows. The next section describes the problem statement and discusses approaches to solve it. The following section describes the results obtained using numerical calculations and a self-similar solution. The last section contains the main conclusions of the paper.

PROBLEM STATEMENT

Basic equations

In the completed series of observations of O1 (from September 12, 2015 to January 19, 2016) and O2 (from November 30, 2016 to August 25, 2017), LIGO and Virgo gravitational wave detectors discovered 10 events of merging binary black holes (Abbott et al. 2019), of which only 5 do not have artefacts. These events are shown in the Table I along with estimates of the component masses, the fraction of mass loss, and the distance to the corresponding object. We can see from the table, that all these objects correspond to black holes of stellar mass. At the same time, the share of mass loss when black holes merges is about (2–6)%. In this paper we will consider the event GW170814 as the object under study, corresponding to the merging of a binary black hole with component masses of 30.5 $M_{\odot }$ and 25.3 $M_{\odot }$. For certainty, we will assume that the mass loss in this event was 5%.

Table I
Parameters of merging binary black holes according to data gravitational-wave observations. The columns indicate (from left to right): the name of the object, the mass of components, the relative change in mass, the distance to object, and reference.

Object	$M_1/M_{\odot }$, $M_2/M_{\odot }$	$\Delta M/(M_1+M_2)$, %	$D$, Mps	Reference
GW150914	${36}_{-4}^{+5}$, ${29}_{-4}^{+4}$	3.3–5.2	${410}_{-180}^{+160}$	(Abbott et al. 2016a)
GW151226	${14.2}_{-3.7}^{+8.3}$, ${7.5}_{-2.3}^{+2.3}$	2.9–5.4	${440}_{-190}^{+180}$	(Abbott et al. 2016b)
GW170104	${31.2}_{-6.0}^{+8.4}$, ${19.4}_{-5.9}^{+8.3}$	2.3–5.7	${880}_{-390}^{+450}$	(Abbott et al. 2017a)
GW170608	${12}_{-2}^{+7}$, $7_{-2}^{+2}$	2.8–5.1	${340}_{-140}^{+140}$	(Abbott et al. 2017b)
GW170814	${30.5}_{-3.0}^{+5.7}$, ${25.3}_{-4.2}^{+2.8}$	4.0–5.8	${540}_{-210}^{+130}$	(Abbott et al. 2017c)

When studying single black holes, it is often assumed that they are surrounded by accretion disks. In our model, it is assumed that the system consists of two black holes surrounded by a circumbinary disk. Obviously, such a disk cannot be formed at the stage of main-sequence stars or giant stars (precursors of black holes), because a strong stellar wind and supernova explosions will destroy it (Cherepashchuk 2016). If the lifetime of a binary system before merging is long, several scenarios for forming a disk can be implemented. For example, the source of the disk matter may be an interstellar molecular cloud, wind from a third star in a hierarchical system, or the remains of a star destroyed by tidal action. As can be seen from the analysis of the considered scenarios, the probability of disk formation is quite high and directly depends on the lifetime of the system. The lifetime of a binary black hole before merging can be estimated as the characteristic time of angular momentum loss by a close binary system with circular orbits. To merge a binary black hole with the parameters of object GW170814 in Hubble time, it is necessary that initial inter-component distance does not exceed 48 $R_{\odot }$. Thus, under any reasonable assumption about the initial distance between the components of a binary black hole of stellar mass, the lifetime of the system will be sufficient for the formation of a circumbinary accretion disk.

For single black holes, disks can exist starting from several gravitational radii. It is known (see e.g. (Misner et al. 1973)) that there are no stable circular orbits inside the zone of \(3R_ {\rm g}\), where \(R_{\rm g} = 2GM/c^2\) — gravitational radius, $G$ — the gravitational constant, $c$ — speed of light. It means that the internal radius of the accretion disk around a single black hole should be larger than \(3R_ {\rm g}\). In the model under consideration it is assumed that a binary black hole is surrounded by a common accretion disk. The inner radius of the accretion disk around a binary black hole is about twice the distance between the centers of black holes. This conclusion follows from estimates of the stability of orbits relative to tidal forces (Bisikalo et al. 2019, Bisikalo & Zhilkin 2020) and is in good agreement with the results of numerical simulation of binary systems with circular orbits of components whose mass ratio is close to unity (Artymowicz & Lubow 1994, Bowen et al. 2019, Kaigorodov et al. 2010, Sytov et al. 2011).

The characteristic lifetime of a binary black hole is long, since the reduction of the inter-component distance due to the emission of gravitational waves is very slow. The reduction time of separation $A$ from 3 $R_{\odot }$ to 3 gravitational radii is approximately 10${}^6$ years, so the circumbinary disk has time to adjust to these slow changes and can be considered stationary at any time. The final stage of merging black holes happen very quickly — in a fraction of a second. The process of rapid merging of black holes begins when the distance between them is \(A = 3R_{{\rm g}, 1} + 3R_{{\rm g},2} = 3R_{\rm g}\), where \(R_{{\rm g}, 1}\) and \(R_{{\rm g}, 2}\) — gravitational radii of components, \(R_{\rm g}\) — gravitational radius calculated from their combined mass $M=M_1+M_2$ (gravitational radius of accretor) (Pretorius 2005, Campanelli et al. 2006, Baker et al. 2006, Cohen et al. 2012). From the point of view of the disk, the final stage of merging can be considered almost instantaneous. Therefore, we can assume that directly before the merger of black holes, the internal radius of the disk \(r_ {\rm d}\) is about 6 times greater than the gravitational radius of the accretor \(R_{\rm g}\). This allows us to ignore relativistic effects and use the theory of Shakura & Sunyaev 1973 when describing the properties of such an accretion disk.

In this model, the accretion rate $\dot{M}$ is determined by the rate of dissipation of the angular momentum that occurs due to the turbulent viscosity. The intensity of turbulence is described by the dimensionless parameter $\alpha$, whose value does not exceed unit. The critical accretion rate is determined by the value of the Eddington luminosity

$${\dot{M}}_{\mathrm{c}\mathrm{r}}=1.8\cdot {10}^{-9}\frac{m}{\chi }\frac{M_{\odot }}{\mathrm{year}},$$(1)

where $m=M/M_{\odot }$ — the dimensionless mass of the central object. The parameter $\chi$ determines the accretion efficiency and in the case of a single non-rotating black hole is approximately equal to 0.06–0.08 (Lipunov 1992). Rotation increases the efficiency of stationary accretion to 0.32 (Bardeen 1970, Thorne 1974). In episodic accretion scenarios, the efficiency can reach 0.43 (Li & Paczyński 2000). However, in our model, the specific value of this parameter is not critical, since the dimensionless accretion rate \(\dot{m} = \dot{M}/\dot{M}_{\rm cr}\) plays a determining role.

In the model of Shakura & Sunyaev 1973, three characteristic zones (A, B, and C) are distinguished in the accretion disk, in which the distribution of hydrodynamic quantities is determined by the relative contribution of gas pressure and radiation pressure, as well as the opacity mechanism. In zone A, the radiation pressure prevails, and in zones B and C, it can be neglected compared to gas pressure. In zones A and B, opacity is determined by Thomson scattering, and in zone C, opacity is provided by free-free transitions. The vertical structure of the disk in these zones differs quite significantly, but the radial dependences of density and temperature are close to power laws with almost the same degree indices. Therefore, in all our calculations we use the following power distributions of density and temperature in an undisturbed disk:

$${\rho }_0\left(r\right)={\rho }_{*}{\left(\frac{r}{r_{*}}\right)}^{-k_d},T_0\left(r\right)=T_{*}{\left(\frac{r}{r_{*}}\right)}^{-k_t},$$(2)

where ${\rho }_{*}$ and $T_{*}$ — the corresponding values at the point $r=r_{*}$. In particular, for zone B, the indices $k_d=33/20$, $k_t=9/10$. The characteristic radius $r_{*}$ can be defined as the internal radius of the disk \(r_{\rm d}\) or as the inner radius of the zone under consideration. For example, if zone B is considered, the radius of the inner boundary $r_{AB}$ of this zone is determined from the condition equality of radiation pressure and gas pressure (Shakura & Sunyaev 1973),

$$r_{AB}=159.3{\alpha }^{2/21}{m}^{2/21}{\dot{m}}^{16/21}{R}_{\mathrm{g}}.$$(3)

The model of Shakura & Sunyaev 1973 is only used to describe the undisturbed state of an accretion disk. Since the process of relaxation of the disk after merging black holes and mass loss of the central object is relatively fast, the further evolution of the disk can be considered in the approximation of non-dissipative (without taking into account the viscosity) gravitational gas dynamics (Bisikalo et al. 2019, Bisikalo et al. 2019).

The equations of gravitational gas dynamics describing the evolution of a geometrically thin accretion disk after black holes merging can be written in the axisymmetric approximation in cylindrical coordinates ($r$, $\phi$, $z$) as follows:

$$\frac{\partial \rho }{\partial t}+v_r\frac{\partial \rho }{\partial r}+v_z\frac{\partial \rho }{\partial z}+\rho [\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_r\right)+\frac{\partial v_z}{\partial z}]=0,$$(4)

$$\frac{\partial v_r}{\partial t}+v_r\frac{\partial v_r}{\partial r}+v_z\frac{\partial v_r}{\partial z}-\frac{v_{\phi }^2}{r}+\frac{1}{\rho }\frac{\partial P}{\partial r}+(1-\xi )\frac{GM}{r^2}=0,$$(5)

$$\frac{\partial v_z}{\partial t}+v_r\frac{\partial v_z}{\partial r}+v_z\frac{\partial v_z}{\partial z}+\frac{1}{\rho }\frac{\partial P}{\partial z}+(1-\xi )\frac{GM}{r^3}z=0,$$(6)

$$\frac{\partial v_{\phi }}{\partial t}+v_r\frac{\partial v_{\phi }}{\partial r}+v_z\frac{\partial v_{\phi }}{\partial z}+\frac{v_r{v}_{\phi }}{r}=0,$$(7)

$$\frac{\partial \epsilon }{\partial t}+v_r\frac{\partial \epsilon }{\partial r}+v_z\frac{\partial \epsilon }{\partial z}+\frac{P}{\rho }[\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_r\right)+\frac{\partial v_z}{\partial z}]=0,$$(8)

where $\rho$ — density, $v_r$ — radial velocity, $v_z$ — vertical velocity, $v_{\phi }$ — azimuthal velocity, $P$ — pressure, $\epsilon$ — specific internal energy, $M$ — total mass of a binary black hole before merging, $\xi$ — fraction of the mass of black holes radiated as gravitational waves.

In the case of geometrically thin disks, we can neglect the detailed description of the vertical structure to study the effects of shock heating and adiabatic expansion. Instead, it can be used a simple approximation when the vertical velocity $v_z$ at any given time remains proportional to the height $z$ (analogous to Hubble’s law in cosmology):

$$v_z(r,z,t)=\psi (r,t)z.$$(9)

In the plane of symmetry of the disk, the vertical velocity component $v_z=0$. In our model we will write all the equations (4)–(8) for the plane of symmetry of the disk. In this case the term $v_z\partial /\partial z$ in convective derivatives disappears in all equations, and the equations of continuity (4) and energy (8) will take the form:

$$\frac{\partial \rho }{\partial t}+v_r\frac{\partial \rho }{\partial r}+\rho \frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_r\right)+\rho \psi =0,$$(10)

$$\frac{\partial \epsilon }{\partial t}+v_r\frac{\partial \epsilon }{\partial r}+\frac{P}{\rho }\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_r\right)+\frac{P\psi }{\rho }=0.$$(11)

The last term in the energy equation (11) describes the adiabatic heating of the disk due to its vertical motion.

Suppose that the pressure in the disk at any time is determined by the expression

$$P(r,z,t)=P(r,t)e^{-z^2/H^2(r,t)},$$(12)

where $H$ is a local half-thickness of the disk. Then from the equation for the vertical velocity (6) with the help of relation (9) we can obtain the equation for the coefficient $\psi$,

$$\frac{\partial \psi }{\partial t}+v_r\frac{\partial \psi }{\partial r}+{\psi }^2-\frac{2P}{\rho H^2}+(1-\xi )\frac{GM}{r^3}=0.$$(13)

Substituting $z=H(r,t)$ into the relation (9), we find

$$v_z(r,H(r,t),t)=\psi (r,t)H(r,t).$$(14)

On the other hand, we can write

$$v_z(r,H(r,t),t)=\frac{\partial H}{\partial t}+v_r\frac{\partial H}{\partial r}.$$(15)

Therefore the local half-thickness of the disk $H(r,t)$ satisfies the equation

$$\frac{\partial H}{\partial t}+v_r\frac{\partial H}{\partial r}-\psi H=0.$$(16)

To complete the resulting system of equations, it is necessary to take into account the equations of state that determine the dependence of the pressure $P$ and the internal energy $\epsilon$ on density $\rho$ and temperature $T$. Since in this paper we consider only the zone B of the accretion disk, we do not take into account the radiation pressure. However, it should be noted that in strongly disturbed regions of this zone, where shock waves of high intensity occur, the temperature can increase so much that the effects of radiation pressure can play a significant role. The equations of state can be written as:

$$P=A_{\mathrm{g}\mathrm{a}\mathrm{s}}\rho T,\epsilon =\frac{A_{\mathrm{g}\mathrm{a}\mathrm{s}}T}{\gamma -1},$$(17)

where \(A_{\rm gas} = 2 k_{\rm B} / m_ {\rm p}\) — the gas constant, \(k_{\rm B}\) — the Boltzmann constant, \(m_{\rm p}\) — proton mass, $\gamma$ — adiabatic index. For a fully ionized hydrogen plasma, the mean molecular weight should be assumed to be $1/2$, and the adiabatic index $\gamma =5/3$.

Numerical model

For the numerical solution of the obtained equations, it is convenient to rewrite them in dimensionless form. To do this, we will select the following values as dimensional scales: $r_0=r_{*}$, $v_0=c_{*}$, $t_0=r_{*}/c_{*}$, ${\rho }_0={\rho }_{*}$, $P_0={\rho }_{*}{c}_{*}^2$, ${\epsilon }_0=c_{*}^2$, $T_0=T_{*}$. Here we use the notation for the speed of sound $c_{*}$, which is determined by the ratio $P_{*}={\rho }_{*}{c}_{*}^2$, where $P_{*}$ is the pressure at the point $r_{*}$.

To simplify the writing, when describing a numerical model, dimensionless quantities will be denoted by the same symbols as the corresponding dimensional ones. Then the system of equations obtained in the previous section will take the following form:

$$\frac{d\rho }{dt}+\rho \frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_r\right)+\rho \psi =0,$$(18)

$$\frac{d{v}_r}{dt}-\frac{v_{\phi }^2}{r}+\frac{1}{\rho }\frac{\partial P}{\partial r}+(1-\xi )\frac{{\mu }^2}{r^2}=0,$$(19)

$$\frac{d{v}_{\phi }}{dt}+\frac{v_r{v}_{\phi }}{r}=0,$$(20)

$$\frac{d\epsilon }{dt}+\frac{P}{\rho }\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_r\right)+\frac{P\psi }{\rho }=0,$$(21)

$$\frac{dr}{dt}=v_r,$$(22)

$$\frac{d\psi }{dt}+{\psi }^2-\frac{2P}{\rho H^2}+(1-\xi )\frac{{\mu }^2}{r^3}=0.$$(23)

$$\frac{dH}{dt}=\psi H,$$(24)

$$P=\rho T,$$(25)

$$\epsilon =\frac{T}{\gamma -1}.$$(26)

In these equations we use the notation for the full derivative is used,

$$\frac{d}{dt}=\frac{\partial }{\partial t}+v_r\frac{\partial }{\partial r}.$$(27)

The parameter

$$\mu =\sqrt{\frac{GM}{r_{*}{c}_{*}^2}}$$(28)

corresponds to the Mach number at the point $r=r_{*}$.

To implement a numerical method for solving this system of equations, we need to pass to the mass Lagrangian variable. However, the presence of an additional term $\rho \psi$ in the continuity equation (18) implies that the mass of the fluid column of unit height is not conserved:

$$\frac{d}{dt}\underset{0}{\overset{r}{\int }}\rho rdr\ne 0.$$(29)

Therefore, this value cannot be used as a mass Lagrangian variable. To solve this problem, we use the method of virtual fluid. Namely, along with the real fluid with density $\rho$, which satisfies the equation (18), we will consider a virtual fluid with a density $\tilde{\rho }$ that satisfies the equation

$$\frac{d\tilde{\rho }}{dt}+\tilde{\rho }\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_r\right)=0,$$(30)

that does not include the corresponded term $\tilde{\rho }\psi$. It is obvious that the virtual fluid describes the same disk, but it is in strict vertical hydrostatic equilibrium.

In this case, we can define a mass Lagrangian variable of the following view:

$$q=\underset{0}{\overset{r}{\int }}\tilde{\rho }rdr.$$(31)

It is easy to see that the relations

$$\frac{\partial }{\partial r}=r\tilde{\rho }\frac{\partial }{\partial q},\frac{1}{\tilde{\rho }}=r\frac{\partial r}{\partial q}$$(32)

are fulfilled. In Lagrangian variables, the equation (30) takes the form:

$$\frac{d}{dt}\left(\frac{1}{\tilde{\rho }}\right)=\frac{\partial }{\partial q}\left(r{v}_r\right).$$(33)

It is easy to check that taking into account the equation for the radial coordinate (22), the equation (33) and the second equation in (32) are equivalent.

Denote by $\eta$ the ratio of densities of real and virtual fluids,

$$\eta =\frac{\rho }{\tilde{\rho }}.$$(34)

Then the second relation in (32) can be rewritten as

$$\eta =\rho r\frac{\partial r}{\partial q}.$$(35)

The equations (18), (19) and (21) in mass Lagrangian variables will take the form:

$$\frac{d}{dt}\left(\frac{1}{\rho }\right)=\frac{1}{\eta }\frac{\partial }{\partial q}\left(r{v}_r\right)+\frac{\psi }{\rho },$$(36)

$$\frac{d{v}_r}{dt}=-\frac{r}{\eta }\frac{\partial P}{\partial q}+\frac{v_{\phi }^2}{r}-(1-\xi )\frac{{\mu }^2}{r^2},$$(37)

$$\frac{d\epsilon }{dt}=-\frac{P}{\eta }\frac{\partial }{\partial q}\left(r{v}_r\right)-\frac{P}{\rho }\psi .$$(38)

The other equations in the system (18)–(26) remain the same.

Let’s introduce a difference mesh in the computational domain in the coordinate system associated with the mass Lagrangian coordinate $q$, the structure of which is determined by the distribution of coordinates of nodes $q_i$, where the index $i=0,1,\dots ,N$. The mesh step $\Delta q_{i+1/2}=q_{i+1}-q_i$ is chosen so that the logarithmic scale of the distance between neighboring nodes is the same, $ln\left(q_{i+1}/q_i\right)=\mathrm{\text{const}}$. We will refer the radial coordinate $r$ and the velocities $v_r$, $v_{\phi }$ at time $t^n$ to the mesh nodes: $r_i^n$, $v_{r,i}^n$, $v_{\phi ,i}^n$, and the thermodynamic values $\rho$, $P$, $\epsilon$, $T$ as well as the values $\eta$, $\psi$, $H$ to the centers of the cells which are numbered with half-integer indices: ${\rho }_{i+1/2}^n$, $p_{i+1/2}^n$, ${\epsilon }_{i+1/2}^n$, $t_{i+1/2}^n$, ${\eta }_{i+1/2}^n$, ${\psi }_{i+1/2}^n$, $h_{i+1/2}^n$ (Samarskii & Popov 1980).

To simplify further expressions we will use the following notation for finite difference operators:

$${\left({\Delta }_{t}f\right)}^{n+1/2}=\frac{f^{n+1}-f^n}{\Delta t},{\left({\Delta }_{q}f\right)}_{i+1/2}=\frac{f_{i+1}-f_i}{\Delta q_{i+1/2}},{\left({\Delta }_{q}f\right)}_i=\frac{f_{i+1/2}-f_{i-1/2}}{\Delta q_i},$$(39)

where $\Delta q_i=(\Delta q_{i+1/2}+\Delta q_{i-1/2})/2$. At the same time, wherever this does not cause misunderstandings, indices from operators can be omitted.

The equations (18)–(26) can be approximated using the following finite-difference scheme:

$${\Delta }_t\frac{1}{{\rho }_{i+1/2}}=\frac{1}{{\eta }_{i+1/2}^{n+1/2}}{\left[{\Delta }_q\left(r^{n+1/2}{v}_r^{n+1/2}\right)\right]}_{i+1/2}+\frac{{\psi }_{i+1/2}^{n+1/2}}{{\rho }_{i+1/2}^{n+1/2}}.$$(40)

$${\Delta }_t{v}_{r,i}=-\frac{r_i^{n+1/2}}{{\eta }_i^{n+1/2}}{\left({\Delta }_q{\Pi }^{n+\sigma }\right)}_i+\frac{{\left(v_{\phi ,i}^{n+1/2}\right)}^2}{r_i^{n+1/2}}-(1-\xi )\frac{{\mu }^2}{r_i^{n+1}{r}_i^n}.$$(41)

$${\Delta }_t{v}_{\phi ,i}=-\frac{v_{r,i}^{n+1/2}{v}_{\phi ,i}^{n+1/2}}{r_i^{n+1/2}},$$(42)

$${\Delta }_t{\epsilon }_{i+1/2}=-\frac{{\Pi }_{i+1/2}^{n+\sigma }}{{\eta }_{i+1/2}^{n+1/2}}{\left[{\Delta }_q\left(r^{n+1/2}{v}_r^{n+1/2}\right)\right]}_{i+1/2}-\frac{P_{i+1/2}^{n+1/2}{\psi }_{i+1/2}^{n+1/2}}{{\rho }_{i+1/2}^{n+1/2}}.$$(43)

$${\Delta }_t{r}_i=v_{r,i}^{n+1/2},$$(44)

$${\Delta }_t{\psi }_{i+1/2}=-{\left({\psi }_{i+1/2}^{n+1/2}\right)}^2+\frac{2P_{i+1/2}^{n+1/2}}{{\rho }_{i+1/2}^{n+1/2}{\left(H_{i+1/2}^{n+1/2}\right)}^2}-(1-\xi )\frac{{\mu }^2}{{\left(r_{i+1/2}^{n+1/2}\right)}^3},$$(45)

$${\Delta }_t{H}_{i+1/2}={\psi }_{i+1/2}^{n+1/2}{H}_{i+1/2}^{n+1/2},$$(46)

$$P_{i+1/2}^{n+1}={\rho }_{i+1/2}^{n+1}{T}_{i+1/2}^{n+1},$$(47)

$${\epsilon }_{i+1/2}^{n+1}=\frac{T_{i+1/2}^{n+1}}{\gamma -1},$$(48)

$${\eta }_{i+1/2}^{n+1}=\frac{1}{2}{\rho }_{i+1/2}^{n+1}{\left[{\Delta }_q{\left(r^{n+1}\right)}^2\right]}_{i+1/2},$$(49)

$${\Pi }_{i+1/2}^{n+1}=P_{i+1/2}^{n+1}+{\omega }_{i+1/2}^{n+1},$$(50)

$${\omega }_{i+1/2}^{n+1}=\Omega ({\rho }_{i+1/2}^{n+1},r_i^{n+1},r_{i+1}^{n+1},v_i^{n+1},v_{i+1}^{n+1}),$$(51)

where

$${\Pi }_{i+1/2}^{n+\sigma }=\sigma {\Pi }_{i+1/2}^{n+1}+(1-\sigma ){\Pi }_{i+1/2}^n,$$(52)

$$r_{i+1/2}=\frac{1}{2}(r_i+r_{i+1}).$$(53)

The parameter σ characterizes the degree of implicitness of the scheme and varies from 0 to 1. This parameter determines the order of approximation in time. In the case of σ = 1/2, the scheme has a second order of approximation in $\Delta t$, and in other cases — the first order. The value $\omega$ describes the artificial viscosity, that is necessary for a more correct description of solutions with shock waves. The explicit form of the function $\Omega$ is defined by a specific artificial viscosity model. In our calculations we used a linear viscosity of the following form

$$\omega =-\frac{\nu \rho }{2}(\frac{\partial v_r}{\partial r}-\left|\frac{\partial v_r}{\partial r}\right|),$$(54)

where the coefficient ${\nu }_{i+1/2}={\nu }_0\Delta q_{i+1/2}$. The value ${\nu }_0$ was set in the range from 3 to 5.

It can be shown that this difference scheme is fully conservative (Samarskii & Popov 1980) in the sense that it ensures not only a balance of the total energy, but also of the individual types of energy (thermal, kinetic, rotational, and gravitational). In addition, it can easily be shown that the difference relation satisfies in the scheme

$${\Delta }_t\left(r_i{v}_{\phi ,i}\right)=0,$$(55)

describing the law of conservation of angular momentum $l=r{v}_{\phi }$. This means that the following equality is implied

$$r_i^n{v}_{\phi ,i}^n=r_i^0{v}_{\phi ,i}^0=l_i,$$(56)

where $l_i$ is the mesh function for specific angular momentum. It follows

$$v_{\phi ,i}^n=\frac{l_i}{r_i^n}.$$(57)

The system of algebraic equations (40)–(51) is non-linear. To solve this system, we use a combination of Newton’s method and the sweep method. The expressions (2) are used as initial conditions for density and temperature. The initial distribution of pressure and the internal energy can be obtained from the equations of state (25) and (26). In addition, we assume that at the initial moment $v_r(r,0)=0$, $\psi (r,0)=0$, $\eta (r,0)=1$. Initial conditions for the azimuthal velocity and half-thickness of the disk can be obtained from stationary solutions of the equations (19) and (23):

$$v_{\phi }(r,0)={[\frac{{\mu }^2}{r}-(k_d+k_t)r^{-k_t}]}^{1/2},$$(58)

$$H(r,0)={\left[\frac{2P(r,0)}{\rho (r,0)}\frac{r^3}{{\mu }^2}\right]}^{1/2}=\frac{\sqrt{2}}{\mu }{r}^{\frac{3-k_t}{2}}.$$(59)

The numerical model is defined by the following parameters: the parameters $k_d$ and $k_t$ of the initial density and temperature profiles, the fraction of mass loss $\xi$, and the Mach number $\mu$. Note that in the special case $\psi =0$, this numerical model turns into the simpler one we considered in paper (Bisikalo et al. 2019), which does not take into account the effects caused by adiabatic disk vertical expansion.

Analytical approach

The analysis of numerical solutions we obtained in (Bisikalo et al. 2019) allows us to conclude that the shock waves arising in a perturbed disk have self-similar nature. In particular, the distance travelled by the shock wave front during the time $t$ is described by the power law $t^{2/3}$, and its velocity decreases over time as $t^{-1/3}$. In (Bisikalo et al. 2019), we explored the possibility of constructing an appropriate self-similar solution. However, it was not possible to obtain it completely due to the difficulties associated with analytical definition for positions of shock waves. A complete self-similar solution is obtained in (Bisikalo & Zhilkin 2020) based on a comparison of self-similar and numerical solutions.

A self-similar solution describes the evolution of a perturbed disk when shock waves travel far enough away from its internal boundary. We can assume that this region corresponds to zone B of the accretion disk, and the radiation pressure already does not significantly affect the dynamics of matter behind the shock wave front. In (Bisikalo & Zhilkin 2020), we also did not take into account the effects associated with adiabatic disk expansion when constructing a self-similar solution. This simplified numerical model we considered in (Bisikalo et al. 2019).

The initial distributions of all values with the exception of azimuthal velocity are described by power functions. The initial distribution for angular momentum (in dimensional values)

$$l_0\left(r\right)=r_{*}{c}_{*}{[{\mu }^2-(k_d+k_t){\left(\frac{r}{r_{*}}\right)}^{1-k_t}]}^{1/2}{\left(\frac{r}{r_{*}}\right)}^{1/2}.$$(60)

This expression shows that in the special case of $k_t=1$ (which is close to $9/10$ in the model of Shakura & Sunyaev 1973 for zone B), this distribution also turns out to be power-law,

$$l_0\left(r\right)=r_{*}{c}_{*}{({\mu }^2-k_d-1)}^{1/2}{\left(\frac{r}{r_{*}}\right)}^{1/2}.$$(61)

Note that this expression is real if ${\mu }^2>k_d+1$. For example, in the case of $k_d=33/20$, it should be $\mu >1.63$.

We will look for a solution in the following form (Sedov 1982):

$$\rho (r,t)={\rho }_0\left(r\right)\sigma \left(\lambda \right),$$(62)

$$v_r(r,t)=\frac{r}{t}u\left(\lambda \right),$$(63)

$$l(r,t)=l_0\left(r\right)\Lambda \left(\lambda \right),$$(64)

$$\epsilon (r,t)={\epsilon }_0\left(r\right)w\left(\lambda \right),$$(65)

where the functions ${\rho }_0\left(r\right)$ and $l_0\left(r\right)$ are defined by the expressions (2) and (61), and the function

$${\epsilon }_0\left(r\right)=\frac{c_{*}^2}{\gamma -1}{\left(\frac{r}{r_{*}}\right)}^{-1}$$(66)

corresponds to the initial temperature distribution (2) with the index $k_t=1$. The self-similar functions σ, $u$, $l$, and $w$ depend on the self-similar variable $\lambda$, which we select as

$$\lambda =\frac{r}{r_{*}}{\left(\frac{t_{*}}{t}\right)}^{\delta },$$(67)

where $\delta$ is the self-similarity index, and $t_{*}=r_{*}/c_{*}$. This choice of a self-similar variable implies that the limit at $\lambda \to \infty$ corresponds to the undisturbed state of the disk. Therefore, in the self-similar solution in the limit at $\lambda \to \infty$ functions σ $\to 1$, $u\to 0$, $\Lambda \to 1$, $w\to 1$. The $\delta$ parameter should be defined from the self-similarity condition of the original equations (18)–(21) (in dimensional form for the case of $\psi =0$).

In terms of self-similar variables (62)–(65) and (67), the equations (18)–(21) can be written as:

$$(u-\delta )\frac{dln\sigma }{dln\lambda }+\frac{du}{dln\lambda }+(2-k_d)u=0,$$(68)

$$\begin{array}{c}\hfill (u-\delta )\frac{du}{dln\lambda }+u(u-1)+\frac{w}{{\lambda }^3}(\frac{dln\sigma }{dln\lambda }+\frac{dlnw}{dln\lambda }-k_d-1)-\hfill \\ \hfill -({\mu }^2-k_d-1)\frac{{\Lambda }^2}{{\lambda }^3}+(1-\xi )\frac{{\mu }^2}{{\lambda }^3}=0,\hfill \end{array}$$(69)

$$(u-\delta )\frac{dln\Lambda }{dln\lambda }+\frac{u}{2}=0,$$(70)

$$(u-\delta )\frac{dlnw}{dln\lambda }+(\gamma -1)\frac{du}{dln\lambda }+(2\gamma -3)u=0.$$(71)

All dimensional coefficients are completely reduced in the case of $\delta =2/3$.

The shock wave front in the self-similar solution corresponds to a constant value of the self-similar variable \(\lambda_s = {\rm const}\). It follows that the position of the shock wave

$$r_s={\lambda }_s{r}_{*}{\left(\frac{t}{t_{*}}\right)}^{\delta }\propto t^{\delta },$$(72)

and its velocity

$$D=\frac{d{r}_s}{dt}=\delta \frac{r_s}{t}\propto t^{\delta -1}.$$(73)

The Hugoniot conditions must be performed on the shock wave, which connecting the value before the discontinuity (we mark it with index 1) and after the discontinuity (we mark it with index 2). A simple analysis leads to the following relations for self-similar functions:

$${\sigma }_2{V}_2={\sigma }_1{V}_1,$$(74)

$${\Lambda }_2={\Lambda }_1,$$(75)

$$V_2=\frac{\gamma -1}{\gamma +1}\frac{1}{V_1}(V_1^2+\frac{2}{\gamma -1}{W}_1),$$(76)

$$W_2={\left(\frac{\gamma -1}{\gamma +1}\right)}^2\frac{1}{V_1^2}(V_1^2+\frac{2}{\gamma -1}{W}_1)(\frac{2\gamma }{\gamma -1}{V}_1^2-W_1),$$(77)

where the notation for auxiliary values is used

$$V=u-\delta ,W=\frac{\gamma w}{{\lambda }^3}.$$(78)

From these conditions, in particular, it follows that the mass flow density $j=$ σ $V$ and the angular momentum $\Lambda$ remain continuous when passing through the front shock wave. Since $j\ne 0$ on shock waves, $u$ cannot equal $\delta$ either to the left or to the right of the discontinuity, $u\ne \delta$. Moreover, all shock waves in the disk propagate from the center to the periphery. Hence, $j<0$ and therefore should be $u<\delta$.

From the self-similar equations (68)–(71), two algebraic integrals can be obtained that express the laws of conservation of angular momentum and entropy. These integrals have the following form:

$$\frac{\sigma |u-\delta |}{{\mathrm{\Lambda }}^{2(2-k_d)}}=\mathrm{const},$$(79)

$$\frac{{\sigma }^{2\gamma -3}}{w^{2-k_d}}|u-\delta {|}^{k_d(\gamma -1)-1}=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$$(80)

Since the values $\sigma V$ and $\Lambda$ remain continuous on shock waves, the integration constant in the first algebraic integral (79) has the same value in the entire computational domain. Using the values of self-similar functions σ $=1$, $u=0$, $\Lambda =1$ in the limit at $\lambda \to 1$, we find that this constant is equal to the self-similarity index $\delta$. So we can write

$${\Lambda }^{2(2-k_d)}=\frac{\sigma }{\delta }|u-\delta |.$$(81)

This algebraic integral expresses the law of conservation of angular momentum in a self-similar form. It allows to express explicitly the $\Lambda$ function via the σ and $u$ functions. In contrast to the law of conservation of angular momentum, the law of conservation of entropy holds only in the smoothness region of the solution, since when passing through the shock wave front, as follows from the Hugoniot conditions (74)–(77), the entropy of the gas is not conserved. Therefore, the integration constant in (80) will be different in each smoothness region. In the outer region of smoothness (before the first shock wave), the integration constant can again be determined using the limit values at $\lambda \to \infty$ of self-similar functions. It turns out to be equal ${\delta }^{k_d(\gamma -1)-1}$. So we can write:

$$w^{2-k_d}{\delta }^{k_d(\gamma -1)-1}={\sigma }^{2\gamma -3}{|u-\delta |}^{k_d(\gamma -1)-1}.$$(82)

This relation in self-similar form is the algebraic integral of adiabaticity. It allows (in the outer region of the perturbed disk) to express explicitly the function $\epsilon$ via the functions σ and $u$.

The self-similar solution allows to find the asymptotic behaviour of self-similar functions in the limits at $\lambda \to \infty$ and at $\lambda \to 0$. The first asymptote describes the growth of perturbations in the accretion disk at small times. On the other hand, it corresponds to the state of the disk at any time at large distances $r\gg r_{*}$ from the central object. The analysis shows (Bisikalo & Zhilkin 2020) that in the limit at $\lambda \to \infty$ functions

$$u=\frac{\xi {\mu }^2}{{\lambda }^3},\sigma =1+\frac{1}{2}(k_d+1)\frac{\xi {\mu }^2}{{\lambda }^3}.$$(83)

In the dimensional values, we find:

$$v_r=\xi GM\frac{t}{r^2},$$(84)

$$\frac{\rho -{\rho }_0}{{\rho }_0}=\frac{1}{2}(k_d+1)\xi GM\frac{t^2}{r^3}.$$(85)

In other words, velocity perturbations grows linearly (proportional to $t$), and density perturbations does according to the quadratic law (proportional to $t^2$). At this moment, the rate of growth of perturbations is directly proportional to the value $\xi$. The corresponding expressions for the functions $\Lambda$ and $w$ can be obtained from here using the algebraic integrals of angular momentum (81) and adiabaticity (82).

The second asymptote of self-similar functions in the opposite limit at $\lambda \to 0$ describes the transition of a perturbed accretion disk to a new stationary state at large times $t\gg t_{*}$. On the other hand, due to the self-similar nature of the solution, it corresponds to the state of the disk close to the central object at any time. Let’s assume that in the limit at $\lambda \to 0$, the self-similar functions $u\to 0$, σ$\to$σ_$0$, $w\to w_0$, $\Lambda \to {\Lambda }_0$. Then from the equation (69), we can obtain the relation

$$(k_d+1)w_0+({\mu }^2-k_d-1){\Lambda }_0^2=(1-\xi ){\mu }^2.$$(86)

The algebraic integral of angular momentum (81) allows to write an expression

$${\Lambda }_0^2={\sigma }_0^{1/(2-k_d)}.$$(87)

From these two equations, $w_0$ can be expressed in term of σ_$0$. The remaining parameter σ_$0$ cannot be found analytically directly from the self-similar solution due to the presence of shock waves in it. But it is not difficult to determine it from the numerical solution of self-similar equations.

We built a self-similar solution using the following methodology. First, for a sufficiently large value of the self-similar variable $\lambda ={\lambda }_{max}$, the values $\sigma \left({\lambda }_{max}\right)$ and $u\left({\lambda }_{max}\right)$ were calculated by the formula (83), describing the asymptotes of these quantities in the limit at $\lambda \to \infty$. Values of $w\left({\lambda }_{max}\right)$ and $\Lambda \left({\lambda }_{max}\right)$ were calculated using the algebraic integrals (82) and (81). These values were used as initial conditions. Next, by setting a certain (in general case, it can be a variable) step of the self-similar variable $\Delta \lambda <0$, using the Runge-Kutta numerical scheme of the 4th order, we can find the values for all quantities for $\lambda <{\lambda }_{max}$.

The number of shock waves and the positions of their fronts ${\lambda }_s$, where the index $s$ determines the shock wave number, cannot be found from the self-similar solution itself. However, this necessary information about shock waves can be taken from the corresponding numerical solution. As soon as in the process of numerical integration of self-similar equations, the self-similar variable $\lambda$ becomes equal to ${\lambda }_1$, the value of the functions at the point $\lambda +\Delta \lambda$ is calculated using the Hugoniot conditions (74)–(77). After this, the integration procedure continues until the next shock wave.

RESULTS OF CALCULATIONS

Numerical simulation

As an example, we present the results of numerical calculation for a model with the parameter $\mu =100$. This model we was considering in (Bisikalo & Zhilkin 2020), as one of the physical models corresponding to the GW170814 event. In this case, the dimensionless accretion rate $\dot{m}=0.05$, and the characteristic radius $r_{*}$ was considered as the internal radius of the disk \(r_ {\rm d} = 6 R_{\rm g}\). If the turbulence parameter $\alpha =0.01$, then from the model of Shakura & Sunyaev (1973) we can find the following parameters at a radius of $r_{*}$: density \(\rho_* = 1.03~{\rm g}/{\rm cm}^3\), temperature \(T_* = 5.65 \cdot 10^7~{\rm K}\), speed of sound \(c_* = 9.66 \cdot 10^7~{\rm cm}/{\rm s}\), characteristic time \(t_* = 1.01~{\rm s}\). The initial state of the accretion disk was determined by power distributions of density and temperature (2) with indices $k_d=33/20$ and $k_t=9/10$, which correspond to zone B. The relative change in the mass of the central object as a result of merging black holes was set to equal $\xi =0.05$. The calculation was performed for a mesh with the number of cells $N=5000$. As mentioned above, to increase the accuracy of calculation in the central regions of the disk, the mesh step in mass Lagrangian coordinates was set to that the distances between neighboring nodes were the same on a logarithmic scale.

In Fig. 1–4 numerical distributions of density $\rho$, temperature $T$, radial $v_r$, and azimuthal $v_{\phi }$ velocities, coefficient $\psi$, half-thickness of the disk $H$ and coefficient $\eta$ at two moments of time $t=7.84t_{*}$ and $t=31.4t_{*}$ are presented. Dotted lines show the corresponding initial distributions of values. As we can see from the figures, merging black holes and mass loss of the central object leads to a perturbation of the accretion disk. This disturbance first occurs at the inner edge of the disk, and then propagates radially to its outer parts. The spatial perturbation profile has a complex structure and consists of a number of strong non-linear waves, some of which are shock waves. The profile analysis allows to identify at least two shock waves of sufficiently high intensity for this calculation run. It is interesting to note that the azimuthal velocity distribution $v_{\phi }$ (see the right panel in Fig. 2) and, consequently, the local specific angular momentum $l=r{v}_{\phi }$, practically does not change.

Figure 1
Radial density distribution $𝛒$ (left) and temperature $𝐓$ (right) at two time moments $𝐭=7.84{𝐭}_{*}$ and $𝐭=31.4{𝐭}_{*}$. Initial distributions are shown as a dotted line.

Figure 2
Radial distribution of the radial ${𝐯}_{𝐫}$ (left) and azimuthal ${𝐯}_{𝛗}$ (right) velocities at two time moments $𝐭=7.84{𝐭}_{*}$ and $𝐭=31.4{𝐭}_{*}$. Initial distributions are shown as a dotted line.

Figure 3
Radial distribution of the coefficient $𝛙$ (left) and the half-thickness of the disk $H$ (right) at two time moments $𝐭=7.84{𝐭}_{*}$ and $𝐭=31.4{𝐭}_{*}$. Initial distributions are shown as a dotted line.

Figure 4
Radial distribution of the coefficient $𝛈$ at two time moments $𝐭=7.84{𝐭}_{*}$ and $𝐭=31.4{𝐭}_{*}$. The initial distribution is shown as a dotted line.

In the region of perturbation, the disk structure becomes significantly non-stationary. The distribution of the radial velocity $v_r$ has an oscillatory nature and it changes sign many times (see the left panel in Fig. 2). This means that some parts of the disk move outward ($v_r>0$), while others move inward ($v_r<0$). The vertical structure of the disk experiences strong variations. The coefficient $\psi$, which determines the vertical velocity (see the equation (9)), changes rapidly over time with a pronounced oscillatory nature (see the left panel in Fig. 3). This leads to quasi-periodic changes in the half-thickness of the disk $H$ (see the right panel in Fig. 3), the amplitude of which gradually decays. After the passage of non-linear and shock waves, a new stationary state is formed in the disk, in which the radial and vertical velocities are zero. This new state is characterized by a lower density, higher temperature, and increased disk thickness. Analysis of the radial distribution of the coefficient $\eta$ (see Fig. 4) allows to conclude that as a result of the vertical expansion of the disk, the density in the plane of symmetry drops by more than 2 times compared to the situation when this effect is not taken into account ($\psi =0$).

The passage of shock waves in the disk leads to heating of matter. A new stationary state of the disk that begins to form in its internal part behind the shock waves is characterized by a higher temperature. As a result, the luminosity of the accretion disk increases dramatically. This phenomenon can be consider as an electromagnetic response of the object to the effect of merging black holes and loss of mass due to the radiation of gravitational waves.

The bolometric luminosity of an accretion disk (when viewed from one side) can be calculated using the expression

$$L(t)=2{\pi }^2\underset{r_{\mathrm{d}}}{\overset{\mathrm{\infty }}{\int }}drr{\sigma }_{\mathrm{S}\mathrm{B}}{T}_{\mathrm{e}\mathrm{f}\mathrm{f}}^4(r,t),$$(88)

where σ_SB — the Stefan–Boltzmann constant, \(T_{\rm eff}\) — local (at a given radius $r$) effective temperature. Effective temperature \(T_ {\rm eff}\) is not equal to the temperature in the plane of symmetry $T$, because the accretion disk in zone B is optically thick and therefore radiation comes from the surface layer. In this case, we can assume that the radiation of the disk in zone B has the Planck’s spectrum. The temperature of matter in this surface layer is determined by the processes of heat energy transfer in the vertical direction from the plane of symmetry of the disk to its surface.

In a geometrically thin accretion disk, a simple approximation can be used to estimate the effective temperature, in which it is assumed that \(T_{\rm eff} = f T\), where the coefficient $f$ is determined by a specific process of heat transfer in the vertical direction (Bisikalo et al. 2019, Bisikalo et al. 2019). If the heat energy transfer is carried out by radiant heat conduction, then it can be assumed that the coefficient $f=0.1$. If the energy transfer in vertical direction is determined by convection, then the temperature difference can be estimated by $f=0.5$.

The light curve $L\left(t\right)$, calculated numerically under these assumptions, is shown on the left panel of Fig. 5. Scale factor

$$L_{\ast }=2{\pi }^2{f}^4{r}_{\ast }^2{\sigma }_{\mathrm{S}\mathrm{B}}{T}_{\ast }^4]$$(89)

is determined by the type of heat transfer process. In this version of the calculation (parameter $\mu =100$), the bolometric luminosity of the disturbed disk increases by about three orders of magnitude compared to the luminosity in the undisturbed one. The luminosity grows very quickly, in a time of the order of $0.1t_{*}$, which for this model corresponds to the time of \(0.1~{\rm s}\). Therefore, this process has a flare nature.

Figure 5
Dependence of bolometric luminosity on time $𝐋\left(𝐭\right)$ (left). Electromagnetic radiation spectra ${𝐅}_{𝛎}$ for radiative (thin lines) and convective (bold lines) disks. Solid lines describe the spectra of disturbed disks (at time $𝐭=31.4{𝐭}_{*}$), and dotted lines show the corresponding spectra of undisturbed disks.

Due to the processes of radiative cooling, the heated matter of the disturbed disk should gradually cool down, which will lead to a slow decrease in luminosity. We did not take this process into account in the hydrodynamic model. The characteristic cooling time of the accretion disk can be estimated as follows. The volumetric cooling coefficient in the diffusion approximation of radiation transfer can be written as

$${\mathrm{\Lambda }}_{\mathrm{cool}}=\frac{4{\sigma }_{\mathrm{SB}}{T}^4}{3\kappa \rho H^2},$$(90)

where $\kappa$ is the opacity coefficient, which in zone B is determined by the Thomson scattering. Therefore, the typical cooling time

$$t_{\mathrm{cool}}=\frac{\rho \epsilon }{{\mathrm{\Lambda }}_{\mathrm{cool}}}.$$(91)

From the obtained numerical results, it follows that at the point $r=r_{*}$ the characteristic cooling time \(t_{\rm cool}\) is several units of $t_{*}$ and in this area cooling occurs quickly enough. However, as the distance increases, the cooling time \(t_{\rm cool}\) becomes large. For example, at a distance of $r=20r_{*}$, the cooling time is several hundred times greater than the typical time scale of the problem $t_{*}$. In future works, we plan to take these processes into account more accurately in our numerical model.

The spectral density of the radiation flux from one side of the disk in a perpendicular direction

$$F_{\nu }(t)=2{\pi }^2\underset{r_{\mathrm{d}}}{\overset{\mathrm{\infty }}{\int }}drr{B}_{\nu }(T_{\mathrm{e}\mathrm{f}\mathrm{f}}(r,t)),$$(92)

defined by the Planck function

$$B_{\nu }(T)=\frac{2h{\nu }^3}{c^2}{(e^{\frac{h\nu }{k_{\mathrm{B}}T}}-1)}^{-1},$$(93)

where $\nu$ — frequency, $h$ — the Planck constant. The obtained electromagnetic radiation spectra at the time $t=31.4t_{*}$ are presented on the right panel in Fig. 5. The spectra for radiative disks are shown by thin lines, for convective disks — by bold lines. At the same time, solid lines describe the spectra of disturbed disks, and the dotted ones correspond to the spectra of undisturbed disks. The values of the $F_{\nu }$ threads are given in the absolute units.

The main radiation energy is concentrated in the range from \(1~{\rm keV}\) to \(100~{\rm keV}\), which corresponds to soft and hard X-rays. In addition, a significant portion of the radiation energy comes in the gamma range (\(h\nu \ge 100~{\rm keV}\)). In both cases (radiative and convective disks), the maximum radiation flux during the flare increases by two orders of magnitude. At the same time, the frequency at which the maximum is reached practically does not shift. In the case of convective disks, the radiation flux is two orders of magnitude greater than that from radiative disks. In addition, the maximum value of the flow from convective disks in comparison with radiative disks is shifted to the region of hard X-ray radiation.

Self-similar solution

Self-similar solutions describing the accretion disk perturbation that occurs as a result of black hole merging and mass loss are presented in our paper (Bisikalo & Zhilkin 2020). To complete the picture, we calculated the models with different values of the parameter $\mu$, lying in a fairly wide range of $2\le \mu \le 200$. In this paper, for example, we describe only one model that corresponds to the value of the parameter $\mu =100$.

The self-similar functions $u\left(\lambda \right)$, $\sigma \left(\lambda \right)$, $w\left(\lambda \right)$ and $\Lambda \left(\lambda \right)$, obtained for this solution, shown in Fig. 6 and 7. As we can see from these figures, for this version of solution, three shock waves are formed in the disk, following each other. The first (external) wave has the highest intensity, and the third (internal) has the lowest one. The positions of shock waves are determined by the values of the self-similar variable ${\lambda }_1=5.905$, ${\lambda }_2=4.249$, ${\lambda }_3=3.628$. Note that the self-similar function $\Lambda \left(\lambda \right)$, describing the profile of the local specific angular momentum, changes continuously on shock waves in accordance with the Hugoniot conditions (75). A non-linear precursor propagates before the external shock wave, the position of which can be determined by a local maximum at the point $\lambda \approx 14$. In the limit at $\lambda \to 0$, the function $w\left(\lambda \right)$ tends to a constant value of $w_0=11.256$. This value determines the temperature of the disk in the region behind the shock waves. The corresponding limit value of the $\sigma \left(\lambda \right)$ function equal to σ₀ = 0.982.

Figure 6
Self-similar functions $𝐮\left(𝛌\right)$ (left) and $𝛔\left(𝛌\right)$ (right) for model with parameter $𝛍=100$.

Figure 7
Self-similar functions $𝐰\left(𝛌\right)$ (left) and $𝚲\left(𝛌\right)$ (right) for model with parameter $𝛍=100$.

Using the constructed self-similar solution, it is possible to estimate the luminosity of a perturbed accretion disk. To do this, again assume that the effective temperature \(T_{\rm eff}(r, t) = f T(r, t)\), where the coefficient $f$ is determined by the type of heat transfer process in the vertical direction. The temperature in the plane of symmetry of the disk in a self-similar solution is defined by the expression

$$T(r,t)=T_{*}\frac{r_{*}}{r}w\left(\lambda \right).$$(94)

Hence, the bolometric luminosity of (88) can be represented as

\[\label{eq3-lc1-ss} L(t) = 2 \pi f^4 \sigma_{\rm SB} T_*^4 r_*^2 \left( \frac{t_*}{t} \right)^{4/3} \int\limits_{(t_*/t)^{2/3}}^{\infty} \frac{w^4(\lambda) d\lambda}{\lambda^3},\]

where we set \(r_{\rm d} = r_*\). The integral over the self-similar variable $\lambda$ can be split into the sum of two integrals over the intervals ${\left(t_{*}/t\right)}^{2/3}\le \lambda \le 1$ and $1\le \lambda \le \infty$. In the first integral, the function $w\left(\lambda \right)$ with good accuracy, can be replaced with its limit value $w_0$ at $\lambda \to 0$:

$$\underset{{\left(t_{*}/t\right)}^{2/3}}{\overset{1}{\int }}\frac{w^4\left(\lambda \right)d\lambda }{{\lambda }^3}=\frac{w_0^4}{2}[{\left(\frac{t}{t_{*}}\right)}^{4/3}-1].$$(96)

The second integral $I_{\mu }$ does not depend on time and can be calculated based on the obtained self-similar solution. Its specific value is determined by the parameter $\mu$. In the case of $\mu =100$, it turns out to be equal to $I_{\mu }=7789.75$. As a result, we find the following expression

$$L(t)=2\pi f^4{\sigma }_{\mathrm{S}\mathrm{B}}{T}_{\ast }^4{r}_{\ast }^2{\left(\frac{t_{\ast }}{t}\right)}^{4/3}\{\frac{w_0^4}{2}[{\left(\frac{t}{t_{\ast }}\right)}^{4/3}-1]+I_{\mu }\}.$$(97)

This function increases monotonically over time and tends to a constant value in the limit $t\gg t_{*}$. In this limit, the bolometric luminosity of the perturbed accretion disk is described by the expression

$$L_{max}=\pi w_0^4{f}^4{\sigma }_{\mathrm{S}\mathrm{B}}{T}_{\ast }^4{r}_{\ast }^2.$$(98)

The dependence $L\left(t\right)$, obtained in the self-similar solution for the parameter $\mu =100$, normalized to the limit value $L_{max}$, is shown on the left panel of Fig. 8. Analysis of this light curve allows to conclude that the luminosity of the disturbed disk reaches a constant value of $L_{max}$ in the order of several characteristic times $t_{*}$ (for the model $\mu =100$ this is a few seconds). The increase in luminosity is determined by the factor $w_0^4$, which for the model $\mu =100$ is equal to approximately 16000. In other words, the luminosity of the disk increases by more than 4 orders of magnitude. This is about 20 times more than in the numerical solution discussed in the previous section. However, we recall that in the self-similar solution, we did not take into account the effect of vertical expansion of the disk, which leads to adiabatic cooling of the perturbed region. With further evolution the heated matter will gradually cool down due to radiative cooling and the luminosity of the disturbed disk will slowly decrease. Characteristic cooling time \(t_{\rm cool}\) can be evaluated in the way described in the previous section (see the expressions (90) and (91)).

Figure 8
Self-similar dependence of bolometric luminosity on time $𝐋\left(𝐭\right)$ (left) for model with parameter $𝛍=100$. Spectra of electromagnetic radiation ${𝐅}_{𝛎}$ obtained from self-similar solution for radiative (thin lines) and convective (thick lines) disks. Solid lines describe the spectra of perturbed disks, and dotted lines describe the corresponding spectra of undisturbed disks.

The spectral density of the electromagnetic radiation flux $F_{\nu }\left(t\right)$ can be calculated using the formula (92). Corresponding spectra are shown on the right panel Fig. 8. Thin lines correspond to spectra of radiative disks, and bold lines correspond to convective disks. At the same time, solid lines show the spectra of perturbed disks, and dotted lines the corresponding spectra of undisturbed disks are shown. The nature of the spectra fully corresponds to what we obtained in the numerical solution (see right panel of Fig. 5). As in the numerical solution, the main energy of radiation turns out to be concentrated in the range of soft and hard X-rays, in which quantum energy \(1~{\rm keV} \le h\nu \le 100~{\rm keV}\). However, unlike numerical solution, the maximum radiation flux in both cases (radiative and convective disks) increases by four orders of magnitude during the flare. In this case, the frequency at which the maximum is reached is slightly shifted to the side of higher energy. The radiation flux from convective disks exceeds the flux from radiative ones more than two orders of magnitude, and the maximum the value of the flow from the convective disks is shifted to the hard X-ray radiation.

Thus, in the self-similar solution, we got a stronger (approximately by one-two orders of magnitude) electromagnetic response compared to numerical solution. As noted above, this is due to the fact that in the self-similar solution we did not take into account the effect of vertical expansion of the disk. From here it can be conclude that adiabatic cooling is an important effect, and if it is not taken into account it leads to underestimation of bolometric luminosity and the spectral radiation flux density by approximately one to two orders of magnitude. This should be kept in mind when evaluating the electromagnetic response using the self-similar solution.

SUMMARY

We have considered one possibility of an electromagnetic response from the merging of two black holes. The scenario is based on the assumption that a binary black hole is surrounded by an accretion disk. When black holes merge, the mass of the product decreases due to the radiation of gravitational waves. As a result, the accretion disk receives some additional energy, enough to bring it out of equilibrium. Perturbations of the disk leads to the formation of non-linear and high-intensity shock waves propagating from the central region to the periphery. Shock waves heat the matter, which leads to a sharp increase in the flux of electromagnetic radiation from the disk. Within this model, it is possible to calculate the light curve, the spectral density of the radiation flux, as well as to estimate the characteristic duration of the flare. Calculation results shown that in the case of object GW170814, the burst of luminosity from the accretion disk may be quite sufficient to register this response by modern observation tools.

Accompanying the registration of gravitational-wave events with observations in the electromagnetic channel corresponds to “multi-messenger“ approach to the study of astrophysical objects. Obviously, this allows to get a much larger amount of useful information, and with the help of such methods, we can achieve a deeper understanding of the nature of the phenomena under study. Such approaches have far reaching prospects and should be fully developed.

In the paper we presented the numerical and analytical solutions describing the perturbation of the accretion disk caused by black hole merging and mass loss. We considered the structure of the perturbed disk in the region (the so-called zone B), which is dominated by gas pressure, and the radiation pressure was neglected. In this case, the opacity of matter in the studied region of the disk is determined by the processes of Thomson scattering of electrons. In contrast to the analytical solution, in the numerical model we took into account the effect of vertical expansion of the disk, resulting in adiabatic cooling of matter heated by shock waves. The analytical model is based on a self-similar solution described in our recent paper (Bisikalo & Zhilkin 2020). The analysis of self-similar equations allows us to obtain algebraic integrals of angular momentum and adiabaticity that express the corresponding laws of conservation. In addition, we can find asymptotic relations for self-similar functions that describe the development of the initial perturbation and the transition of the disk to a new stationary state.

The use of the analytical model allows us to approximate the magnitude of the electromagnetic response of a gravitational-wave event without performing time-consuming numerical calculations. However, as shown by comparison with the results of a more complete numerical model, these estimates overstate the intensity of the electromagnetic response by about one or two orders of magnitude. This is due to the fact that the self-similar solution used does not take into account the effects of vertical expansion of the disk, which leads to adiabatic cooling of matter behind the shock waves.

The results of the research presented in this paper show that electromagnetic radiation resulting from the merging of black holes in the framework of the considered scenario is possible and can be observed by the existing tools. However, for the correct interpretation of the received electromagnetic responses it is necessary to get a self-consistent solution that takes into account the entire complex of physical processes. Unfortunately, self-similar solutions can be obtained only in a simpler statement of the problem. Therefore, it is necessary to develop appropriate numerical models.

A comparison of the results of numerical simulation with the self-similar solution presented in the paper shows that the numerical model used adequately reflects the physical processes that occur when a binary black hole merges. Therefore, it seems appropriate to develop the presented numerical model. This development can be related to an increase in the dimension of the problem (two-dimensional and three-dimensional approximation), and taking into account additional physical processes.

ACKNOWLEDGMENTS

This work was supported by the Russian Foundation for Basic Research (grant 19-29-11013).

• ABBOTT BP ET AL. 2016a. Astrophysical Implications of the Binary Black-hole Merger GW150914. ApJ 818(2): L22. doi:10.3847/2041-8205/818/2/L22
• ABBOTT BP ET AL. 2016b. GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence. Phys Rev Lett 116(24): 241103. doi:10.1103/PhysRevLett.116.241103
• ABBOTT BP ET AL. 2017a. GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2. Phys Rev Lett 118(22): 221101. doi:10.1103/PhysRevLett.118.221101
• ABBOTT BP ET AL. 2017b. GW170608: Observation of a 19 Solar-mass Binary Black Hole Coalescence. ApJ 851(2): L35. doi:10.3847/2041-8213/aa9f0c
• ABBOTT BP ET AL. 2017c. GW170814: A Three-Detector Observation of Gravitational Waves from a Binary Black Hole Coalescence. Phys Rev Lett 119(14): 141101. doi:10.1103/PhysRevLett.119.141101
• ABBOTT BP ET AL. 2019. GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs. Physical Review X 9(3): 031040. doi:10.1103/PhysRevX.9.031040
• ARTYMOWICZ P & LUBOW SH. 1994. Dynamics of Binary-Disk Interaction. I. Resonances and Disk Gap Sizes. ApJ 421: 651. doi:10.1086/173679
• BAKER JG, CENTRELLA J, CHOI DI, KOPPITZ M & VAN METER J. 2006. Gravitational-Wave Extraction from an Inspiraling Configuration of Merging Black Holes. Phys Rev Lett 96(11): 111102. doi:10.1103/PhysRevLett.96.111102
• BARDEEN JM. 1970. Kerr Metric Black Holes. Nature 226(5240): 64-65. doi:10.1038/226064a0
• BEKENSTEIN JD. 1973. Gravitational-Radiation Recoil and Runaway Black Holes. ApJ 183: 657-664. doi:10.1086/152255
• BISIKALO DV & ZHILKIN AG. 2020. Shock propagation in accretion discs around merging black holes: self-similar solution. MNRAS 494(4): 5520-5533. doi:10.1093/mnras/staa1088
• BISIKALO DV, ZHILKIN AG & KURBATOV EP. 2019. Possible Electromagnetic Manifestations of Merging Black Holes. Astronomy Reports 63(1): 1-14. doi:10.1134/S1063772919010025
• BISIKALO DV, ZHILKIN AG & KURBATOV EP. 2019. Possible electromagnetic manifestations of merging black holes. Phys Usp 62(11): 1136-1152. doi:10.3367/UFNe.2019.04.038591
• BODE N & PHINNEY S. 2007. Variability in Circumbinary Disks Following Massive Black Hole Mergers. In: APS April Meeting Abstracts. p. S1.010.
• BOWEN DB, MEWES V, NOBLE SC, AVARA M, CAMPANELLI M & KROLIK JH. 2019. Quasi-periodicity of Supermassive Binary Black Hole Accretion Approaching Merger. ApJ 879(2): 76. doi:10.3847/1538-4357/ab2453
• CAMPANELLI M, LOUSTO CO, MARRONETTI P & ZLOCHOWER Y. 2006. Accurate Evolutions of Orbiting Black-Hole Binaries without Excision. Phys Rev Lett 96(11): 111101. doi:10.1103/PhysRevLett.96.111101
• CHEREPASHCHUK AM. 2016. Discovery of gravitational waves: a new chapter in black hole studies. Physics Uspekhi 59(9): 910. doi:10.3367/UFNe.2016.03.037819
• CLARK JPA & EARDLEY DM. 1977. Evolution of close neutron star binaries. ApJ 215: 311-322. doi:10.1086/155360
• COHEN MI, KAPLAN JD & SCHEEL MA. 2012. Toroidal horizons in binary black hole inspirals. Phys Rev D 85(2): 024031. doi:10.1103/PhysRevD.85.024031
• CORRALES LR, HAIMAN Z & MACFADYEN A. 2010. Hydrodynamical response of a circumbinary gas disc to black hole recoil and mass loss. MNRAS 404(2): 947-962. doi:10.1111/j.1365-2966.2010.16324.x
• DE MINK SE & KING A. 2017. Electromagnetic Signals Following Stellar-mass Black Hole Mergers. ApJ 839(1): L7. doi:10.3847/2041-8213/aa67f3
• FITCHETT MJ. 1983. The influence of gravitational wave momentum losses on the centre of mass motion of a Newtonian binay system. MNRAS 203: 1049-1062. doi:10.1093/mnras/203.4.1049
• KAIGORODOV PV, BISIKALO DV, FATEEVA AM & SYTOV AY. 2010. Structure of the circumbinary envelope around a young binary system. Astronomy Reports 54(12): 1078-1083. doi:10.1134/S1063772910120024
• KOCSIS B & LOEB A. 2008. Brightening of an Accretion Disk due to Viscous Dissipation of Gravitational Waves during the Coalescence of Supermassive Black Holes. Phys Rev Lett 101(4): 041101. doi:10.1103/PhysRevLett.101.041101
• LANDAU LD & LIFSHITZ EM. 1975. The classical theory of fields. Oxford: Pergamon Press.
• LI LX & PACZYŃSKI B. 2000. Extracting Energy from Accretion into a Kerr Black Hole. ApJ 534(2): L197-L198. doi:10.1086/312678
• LIPPAI Z, FREI Z & HAIMAN Z. 2008. Prompt Shocks in the Gas Disk around a Recoiling Supermassive Black Hole Binary. ApJ 676(1): L5. doi:10.1086/587034
• LIPUNOV VM. 1992. Astrophysics of Neutron Stars. Springer.
• LIPUNOV VM, POSTNOV KA & PROKHOROV ME. 1997. Black holes and gravitational waves: Possibilities for simultaneous detection using first-generation laser interferometers. Astronomy Letters 23(4): 492-497.
• MACFADYEN AI & MILOSAVLJEVIĆ M. 2008. An Eccentric Circumbinary Accretion Disk and the Detection of Binary Massive Black Holes. ApJ 672(1): 83-93. doi:10.1086/523869
• MCKERNAN B, FORD KES, BARTOS I, GRAHAM MJ, LYRA W, MARKA S, MARKA Z, ROSS NP, STERN D & YANG Y. 2019. Ram-pressure Stripping of a Kicked Hill Sphere: Prompt Electromagnetic Emission from the Merger of Stellar Mass Black Holes in an AGN Accretion Disk. ApJ 884(2): L50. doi:10.3847/2041-8213/ab4886
• MEGEVAND M, ANDERSON M, FRANK J, HIRSCHMANN EW, LEHNER L, LIEBLING SL, MOTL PM & NEILSEN D. 2009. Perturbed disks get shocked: Binary black hole merger effects on accretion disks. Phys Rev D 80(2): 024012. doi:10.1103/PhysRevD.80.024012
• MILOSAVLJEVIĆ M & PHINNEY ES. 2005. The Afterglow of Massive Black Hole Coalescence. ApJ 622(2): L93-L96. doi:10.1086/429618
• MISNER CW, THORNE KS & WHEELER JA. 1973. Gravitation. San Francisco: W. H. Freeman.
• O’NEILL SM, MILLER MC, BOGDANOVIĆ T, REYNOLDS CS & SCHNITTMAN JD. 2009. Reaction of Accretion Disks to Abrupt Mass Loss During Binary Black Hole Merger. ApJ 700(1): 859-871. doi:10.1088/0004-637X/700/1/859
• PERES A. 1962. Classical Radiation Recoil. Physical Review 128(5): 2471-2475. doi:10.1103/PhysRev.128.2471. doi:10.1103/PhysRev.128.2471
• PIETILÄ H, HEINÄMÄKI P, MIKKOLA S & VALTONEN MJ. 1995. Anisotropic Gravitational Radiation in the Problems of Three and Four Black Holes. Celest Mech Dyn Astron 62(4): 377-394. doi:10.1007/BF00692287
• PRETORIUS F. 2005. Evolution of Binary Black-Hole Spacetimes. Phys Rev Lett 95(12): 121101. doi:10.1103/PhysRevLett.95.121101
• ROSOTTI GP, LODATO G & PRICE DJ. 2012. Response of a circumbinary accretion disc to black hole mass loss. MNRAS 425(3): 1958-1966. doi:10.1111/j.1365-2966.2012.21488.x
• SAMARSKII AA & POPOV IP. 1980. Difference methods for solving problems of gas dynamics 2nd revised and enlarged edition. In: Moscow Izdatel Nauka.
• SCHNITTMAN JD & KROLIK JH. 2008. The Infrared Afterglow of Supermassive Black Hole Mergers. ApJ 684(2): 835-844. doi:10.1086/590363
• SEDOV LI. 1982. Similarity and Dimensional Methods in Mechanics. Moscow: Mir Publ.
• SHAKURA NI & SUNYAEV RA. 1973. Black holes in binary systems. Observational appearance. A&A 24: 337.
• SHIELDS GA & BONNING EW. 2008. Powerful Flares from Recoiling Black Holes in Quasars. ApJ 682(2): 758-766. doi:10.1086/589427
• SYTOV AY, KAIGORODOV PV, FATEEVA AM & BISIKALO DV. 2011. Structure of the circumbinary envelopes of young binary stars with elliptical orbits. Astronomy Reports 55(9): 793-800. doi:10.1134/S1063772911090071
• THORNE KS. 1974. Disk-Accretion onto a Black Hole. II. Evolution of the Hole. ApJ 191: 507-520. doi:10.1086/152991
• TUTUKOV AV & YUNGELSON LR. 1993. The merger rate of neutron star and black hole binaries. MNRAS 260: 675-678. doi:10.1093/mnras/260.3.675

Publication Dates

• Publication in this collection
  17 May 2021
• Date of issue
  2021

History

• Received
  24 May 2020
• Accepted
  12 Oct 2020