Fully General Relativistic Magnetohydrodynamic Simulations of Accretion Flows onto Spinning Massive Black Hole Binary Mergers
Federico Cattorini, Bruno Giacomazzo, Francesco Haardt, Monica Colpi
FFully General Relativistic Magnetohydrodynamic Simulations of Accretion Flows onto SpinningMassive Black Hole Binary Mergers
Federico Cattorini , , ∗ Bruno Giacomazzo , , , Francesco Haardt , , , and Monica Colpi , Dipartimento di Scienza e Alta Tecnologia, Universit´a degli Studi dell’Insubria, Via Valleggio 11, I-22100, Como, Italy INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy Dipartimento di Fisica G. Occhialini, Universit`a di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy and INAF, Osservatorio Astronomico di Brera, Via E. Bianchi 46, I-23807 Merate, Italy (Dated: March 1, 2021)We perform the first suite of fully general relativistic magnetohydrodynamic simulations of spinning massiveblack hole binary mergers. We consider binary black holes with spins of different magnitudes aligned to theorbital angular momentum, which are immersed in a hot, magnetized gas cloud. We investigate the effect of thespin and degree of magnetization (defined through the fluid parameter β − ≡ p mag /p fluid ) on the properties ofthe accretion flow. We find that magnetized accretion flows are characterized by more turbulent dynamics, as themagnetic field lines are twisted and compressed during the late inspiral. Post-merger, the polar regions aroundthe spin axis of the remnant Kerr black hole are magnetically dominated, and the magnetic field strength isincreased by a factor ∼ (independently from the initial value of β − ). The magnetized gas in the equatorialplane acquires higher angular momentum, and settles in a thin circular structure around the black hole. We findthat mass accretion rates of magnetized configurations are generally smaller than in the unmagnetized cases byup to a factor ∼
3. Black hole spins have also a suppressing effect on the accretion rate, as large as ∼ ∼ with increasing spin, regardless of the initial level of magnetization ofthe fluid. Our results stress the importance of taking into account both spins and magnetic fields when studyingaccretion processes onto merging massive black holes. I. INTRODUCTION
Massive black hole binary (MBHB) mergers are naturaloutcome of galaxy collisions [1, 2], and are among the mostpowerful sources of gravitational waves (GWs) which will bedetected by future space-based interferometers such as LISA[3]. These mergers may occur in gas rich environment [4–10], leading to the intriguing possibility of concurrent electro-magnetic (EM) emission observable by traditional astronom-ical facilities. Observing these powerful events both in theEM and the GW windows will provide unique opportunitiesfor multimessenger astronomy. A major goal of forthcom-ing multi-band EM observations (e.g.,
Athena [11]) is to ob-serve and study the EM counterparts to LISA MBHB coales-cences: detecting the EM signal emitted alongside an ongoingmerger will let us probe the existence of multiple disk struc-tures around the massive black holes (MBHs) [12–15] and,possibly, the launch of relativistic jets during the inspiral andwhen the new MBH has formed [16, 17]. Concurrent observa-tion of EM counterparts to GW events will help illuminatingthe physical processes that power quasars, and offer new op-portunities for testing the propagation of GWs on cosmologi-cal scales, e.g., measuring the differences in the arrival timesof light and GWs, or inferring the redshift z versus luminos-ity distance d L relation without resorting to EM distance scalecalibrators [18–20].Our knowledge of the properties of the EM signals emerg-ing during a MBHB merger is still incomplete, despite recentadvances [13, 21–24]. Predictions on this EM emission de-pend on the fueling rate, on the hydrodynamical, geometrical ∗ [email protected] and radiative properties of the accreting magnetized gas, andon the MBHs masses and spins.The development of numerical relativity (NR) simulationsof these powerful events is required to advance our theoreticalunderstanding of the physical mechanisms which drive EMsignals associated to GW detections. A jump in the predictivepower of NR simulations will allow better predictions for theEM spectrum rising during the late inspiral and coalescence ofMBHBs, providing guidance to future observations and max-imizing the scientific return of LISA.The structure of the accretion flows around coalescing MB-HBs largely depends on the angular momentum content of theaccreting gas conveyed in the galactic merger, and on its ther-modynamical state. Two limiting scenarios bracket the rangeof physical properties of accreting fluids around MBHBs:(i) The circumbinary disk (CBD) model, where a rotation-ally supported disk surrounds the binary, and in whichviscous and gravitational torques balance to clear a cen-tral cavity at twice the MBHB separation [25]. Nu-merical simulations show that the system evolves into anon-axisymmetric configuration with the cavity becom-ing highly lopsided and filled with a tenuous, shockedplasma, in part ejected against the disk wall where itloses angular momentum to feed the MBHs. This leadsto the formation of two narrow streams, which peri-odically convey mass onto the MBHs in the form oftransient “mini-disks” that persist down to coalescence[e.g., 13, 15, 26–31]. The first simulations of equal-mass, non-spinning binaries in magnetized CBDs wereperformed by [26] (adopting high-order PN approxima-tions) and [32] (in full general relativity).(ii) If the surrounding gas is hot, tenuous, and not rotation- a r X i v : . [ a s t r o - ph . H E ] F e b ally supported, the MBHs may find themselves embed-ded in a turbulent and radiatively inefficient accretionflow [33, 34]. In this scenario, the gas is unable to coolefficiently, and thus the energy is stored in the accre-tion flow as thermal energy instead of being radiated.We refer to this scenario as the gas cloud model [35–37]. The first general relativistic hydrodynamical simu-lations of merging equal-mass binaries in unmagnetizedgas clouds were carried out by [35] and [36] (the lat-ter considered both non-spinning and parallel-spin bi-naries). These works established that the phases of late-inspiral and merger are accompanied by a gradual risein the emitted bremsstrahlung luminosity, followed bya sudden drop-off corresponding to the post-merger ac-cretion of the shock-heated gas. In a subsequent work[37], the impact of misaligned spins and unequal massratios on the physics of hot accretion flows was inves-tigated, and it was found that less symmetric systemsresult in lower luminosity and delayed emission fromthe regions near the BHs.In the present work we consider the hot gas cloud model.We perform the first general-relativistic magnetohydrody-namic (GRMHD) simulations of merging spinning BHs im-mersed in an initially homogeneous fluid, and examine howmagnetic fields and spins affect the dynamics of the gas andthe Poynting luminosity emission. Our simulations revise thescenario analyzed in Giacomazzo et al. 2012 [38] ( Gi12 here-after) and Kelly et al. 2018 [22] ( Ke17 hereafter), and ex-plore the behavior of moderately magnetized accretion flows(MMAFs) onto binaries of MBHs within the ideal MHD limit.The analysis of moderately magnetized plasma bridges thestudy of unmagnetized gaseous environments [e.g., 35, 36]and results obtained in the force-free (FF) regime [e.g., 16],which approximate magnetically dominated plasma (i.e., flu-ids for which β − ≡ p mag /p fluid (cid:29) ).The simulations of Gi12 were the first to study the nature ofMMAFs around equal-mass, non-spinning black hole binaries(set at an initial separation of 8.48 M , where M is the totalmass of the binary), solving the GRMHD equations with the WhiskyMHD code [39, 40]. Gi12 considered two modelsfor the gas cloud surrounding the binary, both with an initiallyuniform rest-mass density ρ : a unmagnetized plasma, and aplasma threaded by an initially uniform magnetic field with aninitial ratio of magnetic-to-fluid pressure β − equal to 0.025.Their results showed that MMAFs exhibit different dynamicscompared to unmagnetized accretion flows, and can lead tostrong, collimated EM emission.The results of Gi12 were farther extended by Ke17, whocovered a broader collection of physical scenarios adoptingthe IllinoisGRMHD code [41, 42] to solve the GRMHDequations. The simulations of Ke17 consider equal-massbinaries of non-spinning black holes with initial separationscovering values between . M and . M . Evolving higher-separation binaries allowed to better resolve the timing fea-tures of the EM (Poynting) emission. Several configurationsdiffering only in the initial magnetic field value b were alsoevolved, showing that the level of the Poynting luminosityreached during the inspiral is little sensitive to the initial mag- netic field strength.In this work, we progress studying the scenario examinedin Gi12- Ke17, and carry out the first 3-dimensional GRMHDsimulations of merging BHs including spins. Extending thestudy to BHs with non-zero spin is key when considering bi-naries of massive black holes. The motivation is astrophysi-cal, as there is observational evidence that MBHs have grownprimarily by efficient accretion [43], and are expected to ac-quire a non-vanishing spin depending on whether accretion isprograde or retrograde, coherent or chaotic, as discussed ex-tensively in the literature [see, e.g., 44–49].Two recent works presented preliminary results describ-ing GRMHD pre-merger simulations of spinning binary blackholes. Lopez Armengol et al. [50] construct an approximatespace-time metric, which they name “Superimposed Kerr-Schild” (SKS) metric, and carry out simulations of circumbi-nary accretion onto binary systems with separation fixed at M , and spin parameters a = (0 , ± . . They find that spincan significantly affect the circumbinary accretion via frame-dragging effects, enhancing or reducing it according to thesign of the spin-orbit coupling. Paschalidis et al. [51] performfully general relativistic MHD simulations of BHBs, and con-sider binary configurations of spinning BHs set at an initialdistance d = 20 M , with spin parameters a = (0 , ± . .Their work addresses the formation and dynamics of mini-disks. In particular, they demonstrate the impact of spin inallowing the formation of mini-disks. Both forementioned in-vestigations focus on the late stages of BHB inspiral, i.e. thepre-merger phase. By contrast, our work considers differentenvironment (i.e., the gas cloud model) and examines merging binary systems, treating both pre- and post-merger phases.Our simulations consider binary equal-mass BHs withequal spins, both aligned with the orbital angular momen-tum, and with spin dimensionless parameters of magnitude a , = (0 , . , . , immersed in a uniform plasma with dif-ferent initial degrees of magnetization. The binary evolu-tions are carried out with the Einstein Toolkit [52]on adaptive-mesh refinement (AMR) grids provided by theCarpet driver [53]. The space-time metric evolution is ob-tained using the Kranc-based McLachlan [54, 55] thorn,adopting the BSSN [56–58] formalism. We adopt the “mov-ing puncture” method [59], and our initial metric data areof the Bowen-York type [60], conditioned to satisfy the con-straint equations using the
TwoPunctures thorn [61]. TheGRMHD equations were solved with the
IllinoisGRMHD code [41, 42].The main structure of the paper is as follows. In Sec. II, wegive a brief description of the numerical methods adopted inour simulations. The initial configuration of our binary evolu-tions are described in Sec. III. In Sec. IV, we present resultsfrom all models considered: the dynamics of the plasma sur-rounding the BHs across evolution (IV A), the magnetic fieldenhancement and the formation of magnetically-dominatedregions (IV B), the mass accretion rate both during the orbital http://einsteintoolkit.org evolution and in the post-merger (IV C), the development ofstrong Poynting flux emission (IV D). II. NUMERICAL METHODS
We consider three families of simulations, each defined byBHs spin parameter a , = (0 , . , . ; for each family, werun three simulations characterized by different degrees of ini-tial “magnetization”, i.e., different values of initial magnetic-to-gas pressure ratio β − . All runs consider black holes im-mersed in an adiabatic gas with initial uniform density andpressure. We take the gas to be either unmagnetized (B0 mod-els) or moderately magnetized (B1 and B2 models), set withinitial uniform magnetic field aligned with the total angularmomentum of the system.In this section we give a brief overview of the mathematicalsetup used for producing the simulations discussed in the fol-lowing. For more detailed discussion on the numerical frame-work adopted for evolving BH binaries in general relativitysee, e.g., [62, 63]. A. Evolution of Gravitational Fields
All the equations presented below are in geometrized units( G = c = 1 ). In these units, Einstein’s field equations ofgeneral relativity are G µν = 8 πT µν (1)where G µν is the Einstein tensor and T µν the total stress-energy tensor. For a magnetized fluid, the stress-energy tensoris the sum of matter and EM components: T µν = T µν matter + T µν EM , (2) T µν matter = ρhu µ u ν + p fluid g µν , (3) T µν EM = b (cid:18) u µ u ν + 12 g µν (cid:19) − b µ b ν , (4)where g µν is the metric tensor, ρ the rest-mass density, u µ the four-velocity of the fluid, h the specific enthalpy, p fluid the fluid pressure and b µ the magnetic four-vector. The space-time metric in standard 3+1 form is ds = − α dt + γ ij (cid:0) dx i + β i dt (cid:1) (cid:0) dx j + β j dt (cid:1) (5)where α is the lapse function, β i the i -th component of theshift vector and γ ij the spatial metric. The extrinsic curvature K ij is given by ( ∂ t − L β ) γ ij = − αK ij (6)with L β denoting the Lie derivative with respect to β . Weevolve the metric variables ( γ ij , K ij ) using the BSSN formu-lation (BSSN evolution and constraint equations are summa-rized in [57, 58]). TABLE I. BBH initial data parameters and derived quantities in codeunits of the GRMHD runs: initial puncture separation d and linearmomentum components p x & p y , dimensionless spin parameter a ofeach BH, merger time t Merger and initial ratio of magnetic-to-fluidpressure β − .Run d p x p y a , t Merger β − B0S0 0B1S0 12.038 5.26 · − · − · − Our metric evolution equations do not include matter sourceterms, since for all the simulations considered in this workwe assume that the total mass of the fluid is negligible withrespect to the mass of the two BHs, M fluid (cid:28) M BHs (i.e.,we evolve the Einstein equations in vacuum). We adopt the“1+log” slicing condition for the lapse and a “hyperbolicgamma-driving” condition for the shift [59].In general, non-radiative GRMHD simulations are scale-free. Thus, we will use length and time units that scale withthe total mass of the system M . All our simulations areevolved setting c = G = M = 1 . The unitary values ofmass, length and time in code units correspond to ˆ m = M = 2 · M g (7) ˆ l = Gc M = 1 . · M cm (8) ˆ t = Gc M = 4 . M s (9)where M ≡ M/ M (cid:12) . Since we assume M fluid (cid:28) M BHs ,the fluid contribution to Eq. (1) can be ignored, and we canset T µν ≈ . This means that we are free to rescale an ap-propriate set of the fluid field variables independently of thegeometric scaling that arises by the condition M = 1 . B. Evolution of magnetohydrodynamic fields
The GRMHD equations and constraint equations are de-rived from1. the conservation of baryon number ∇ µ ( ρu µ ) = 0 , (10)2. the conservation of energy-momentum ∇ µ T µν = 0 , (11)where T µν = T µν matter + T µν EM , and3. the homogenous Maxwell’s equations ∇ ν F ∗ µν = 1 √− g ∂ ν (cid:0) √− gF ∗ µν (cid:1) = 0 , (12)where F µν is the Faraday tensor, F ∗ µν its dual and g the determinant of g µν .The IllinoisGRMHD code evolves a set of conservativeMHD fields C ≡ (cid:110) ρ ∗ , ˜ τ , (cid:126)S, (cid:126)B (cid:111) solving the coupled Einstein-Maxwell equations. It assumes a perfect fluid stress-energytensor for the matter, and infinite conductivity (ideal MHDlimit). The vectors of the “conservative variables” C dependdirectly on the “primitive variables” P ≡ (cid:8) ρ, p fluid , v i , B i (cid:9) where ρ is the rest-mass density, p fluid is the fluid pressure, v i ≡ u i /u are the components of the fluid three-velocity and B i are the spatial components of magnetic field B µ measuredby Eulerian observers. C. Magnetized Accretion Flows
At present, our understanding of accretion flow propertiesaround merging MBHBs is uncertain. On very small scales(such as the ones we consider), it is not possible to uniquelydefine initial conditions for the gas in the vicinity of mergingbinaries. Following Gi12- Ke17, we choose to evolve ourmodels in a simple environment consisting in a homogenous,ideal gas with initial uniform rest-mass density ρ , which hasan initially uniform magnetic field (aligned with the orbital an-gular momentum) and fills the entire computational domain.An ideal gas with adiabatic index γ has a pressure p fluid = ρ(cid:15) ( γ − (13)(where ρ is the rest-mass density, (cid:15) is the specific internal en-ergy, and γ is the adiabatic index), and a specific enthalpy h = (1 + (cid:15) ) + p fluid ρ = 1 + γ(cid:15) (14)We also choose our gas to obey the polytropic equation ofstate p fluid = κρ Γ , (15)with a polytropic index Γ = 4 / and a polytropic constant κ to be assigned. We assume that the adiabatic index of the fluid γ is coincident with the polytropic index Γ ; hence, the specificinternal energy of the gas can be expressed as (cid:15) = κρ Γ − Γ − (16)The speed of sound in the gas is c s = (cid:115) h (cid:18) ∂p fluid ∂ρ + p fluid ρ ∂p fluid ∂(cid:15) (cid:19) = (cid:114) (Γ − (cid:15) + p fluid /ρ ) h = (cid:115) Γ p fluid hρ (17) We express the magnetic field with the magnetic four-vector b µ [see, e.g., 64]: b µ = 1 √ πα (cid:18) u m B m , B i + ( u m B m ) u i u (cid:19) (18)where repeated Latin indices indicate sums over spatial com-ponents only.The relativistic Alfv´en velocity of a magnetized plasma[65] is defined as v Alf = (cid:115) b ρ (1 + (cid:15) ) + p fluid + b = (cid:115) b ρ (1 + Γ (cid:15) ) + b = (cid:115) b ρ + 4 p fluid + b (19)where the second line holds for polytropic fluids with Γ =4 / . D. Diagnostics
To explore the effects of spin and magnetic field strengthon the dynamics of the accreting gas we track the evolution ofthe following quantities:• rest-mass density ρ (normalized to its initial value ρ );• Newtonian Mach number M ≡ v /c s , where v is de-fined as the velocity magnitude of the fluid on the or-bital plane v = (v x + v y ) / (20)and c s is the speed of sound in the medium (Eq. (17));• Newtonian angular velocity Ω fluid of the fluid about theorbital axis, defined as Ω fluid = xv y − yv x ( x + y ) . (21)The quantity Ω fluid has the dimension of ˆ t − .Given a test particle orbiting a Kerr BH of mass M and spinparameter a , the coordinate angular frequency of a circularorbit (for those values of r for which circular orbits exist) is[see, e.g., 66] Ω ± circ = ± M / r / ∓ aM / (22)where r is the areal radius in KBL coordinates, and the sign+(-) refers to corotating (counter-rotating) orbits. We define a circularity parameter ω as ω ≡ Ω fluid / Ω +circ (23)Following the evolution of these diagnostics allows us to bet-ter interpret the results of each simulation (e.g., monitoring M and ω will help us studying the degree of rotation inducedon the gas by the inspiralling BHs).In Sec. IV D we explore how the evolution of the magneticfields affects the possible emission of EM signals. As waspointed out in [16, 67] (for electrovacuum) and [68, 69] (forforce-free plasma), the later inspiral and merger of massiveBH binaries immersed in a magnetized gas may be connectedwith an EM counterpart in the form of a jet, which could bepotentially visible at large distances. In this work we studythis strong and collimated electromagnetic emission lookingprimarily at the Poynting vector [22, 38]. It is calculated as S i ≡ αT i EM , = α (cid:18) b u i u + 12 b g i − b i b (cid:19) (24) III. INITIAL DATA
We reexamine the setup of Ke17 performing simulationsof MMAFs onto binaries of equal-mass BHs. The individ-ual mass of each BH in code units is M BH = M/ . ,where M = 1 is the total mass of the system. Our binaries areimmersed in an initially uniform, radiation-dominated poly-tropic fluid ( p = κρ Γ0 , with ρ = 1 , κ = 0 . , Γ = 4 / ). Tocapture the effect of the individual BHs spins on the accretionflows, we evolve binaries of spinning BHs with parallel spinsaligned with the orbital axis, and adimensional spin parame-ters a = a = (0 , . , . .We adopt a cubical domain given by [ − M, M ] and employ AMR with N = 11 levels of refinement. Thecoarsest resolution is ∆ x c = 64 M/ , and the finest one is ∆ x f = ∆ x c · − N = M/ . All our simulations could beeasily rescaled to consider systems of binary black holes witha total mass M = 2 × M (cid:12) , and immersed in a gas withuniform initial rest-mass density ρ = 10 − g cm − . Thesevalues are consistent with the approximation T µν ≈ , sincethey yield M fluid /M BH ∼ − (see Sec. II A).The BHs rotate around each other starting on quasi-circularorbits at an initial separation d (cid:39) M . Our quasi-circularinitial data are obtained from larger scale PN evolutions (thePN equations are evolved from a larger separation, ∼ M ,to the distance we begin our full GR runs with). In table I wegive the initial data for the 9 configurations presented in thispaper. A. Initial plasma configuration
We evolve MBHBs immersed in a hot plasma, which isthreaded by an initially uniform magnetic field parallel to thebinary angular momentum, i.e. B i = (0 , , B z ) . The mag-netic field is assumed to be anchored to a distant circumbinarydisk located outside the computational domain. This initialconfiguration of the magnetic field is analogous to that imple-mented in previous works (e.g., [16, 68, 69], Gi12, Ke17).Whilst simplistic, our choice of the initial plasma configura-tion is sufficiently clear to aid in pinpointing the effects of sub-tle physical processes (e.g., the spins) on the accretion flows. We set three different initial plasma configurations, which arechosen so that β − ≡ p mag p fluid = (B0 runs) . (B1 runs) . (B2 runs) (25)or, equivalently, ζ ≡ u mag u fluid = (B0 runs) . (B1 runs) . (B2 runs) (26)where β − and ζ are the adimensional magnetic-to-fluid pres-sure ratio and magnetic-to-fluid energy density ratio, respec-tively, and u mag = p mag = B π = b ,u fluid = ρc , p fluid = κρ Γ . In Table II we list the initial uniform GRMHD field values forthe 3 sets of configurations B0, B1 and B2. The value of theinitial fluid rest-mass density ρ = 10 − g cm − , along witha specific choice of ζ ( β − ), uniquely fixes the correspondingphysical values of the initial magnetic field strenght B andof the initial fluid temperature T . In physical units, theadimensional ratios (25) and (26) are β − ≡ p mag p fluid = B πκρ Γ0 (27) ζ ≡ u mag u fluid = B πc ρ (28)Therefore, given a specific value of the adimensional parame-ter ζ , one has B = (cid:112) πc ρ ζ = (B0 runs) . × G (B1 runs) . × G (B2 runs) (29)To calculate the initial physical temperature of the accretionflow we use the equation T = µm p k B p ρ = µm p k B κρ Γ − (30) Run ρ κ ζ β − v alf B0 1 0.2 0.0 0.0 -B1 1 0.2 5e-3 2.5e-2 7.4e-2B2 1 0.2 6.3e-2 0.31 0.26TABLE II. Initial uniform GRMHD field values for the 3 sets of con-figurations B0, B1 and B2: rest-mass density ρ , polytropic constant κ , magnetic-to-gas energy density (pressure) ratio ζ ( β − ), Alfv´envelocity v alf . The values are in code units. − − − x [M] − − − y [ M ] t = 306.67 [M] . . . . . . ρ / ρ − − − x [M] − − − y [ M ] t = 1863.33 [M] . . . . . . ρ / ρ − − − x [M] − − − y [ M ] t = 2500.00 [M] . . . . . . ρ / ρ − − − x [M] − z [ M ] t = 306.67 [M] . . . . . . ρ / ρ − − − x [M] − z [ M ] t = 1863.33 [M] . . . . . . ρ / ρ − − − x [M] − z [ M ] t = 2500.00 [M] . . . . . . ρ / ρ FIG. 1. Rest-mass density ρ (normalized to its initial value ρ ) on the xy (top panels) and xz planes (bottom panels) for the B2S3 configuration( a = 0 . , β − = 0 . ). The snapshots were taken, respectively, after ∼ ∼ ∼ M after themerger. The regions inside the BH horizons have been masked out. where µ is the mean molecular weight, m p is the proton mass, κ is the polytropic constant and k B is the Boltzmann constant.In code units, we set κ = 0 . (conforming with Gi12, Ke17).To assign the physical value of κ in cgs units (which, for a Γ = 4 / polytrope, has the dimensions of g − / cm s − ) weproceed as follows: combining Eqs. (27) and (28), we find κ = ζ β − c ρ − Γ0 (31)Inserting the values in code units for κ, c and ρ in Eq. (31)yields ζ /β − = 0 . , which is adimensional, and thus inde-pendent of the units of measure. Therefore, entering the cgsvalues of c and ρ in (31), we get κ ∼ . × ρ − / − g − / cm s − (32)where ρ − ≡ ρ / − g cm − . Employing Eq. (30)with µ = 1 / , ρ = 10 − g cm − and κ ∼ . × g − / cm s − , we find that the initial temperature of theaccretion flow is T ∼ . × K (33)Our physical values of the initial magnetic field magni-tudes B and initial temperature T are consistent with thoseadopted in other general relativistic simulations of hot accre-tion flows onto MBHBs, e.g. [36, 37], Gi12, Ke17. IV. RESULTS
With our work we probe the physics of MMAFs onto bina-ries of spinning BHs, evolving a number of simulations which cover a range of black hole spins and gas magnetization. Fol-lowing the work of Ke17, we aim at exploring the astrophys-ical processes which may give rise to electromagnetic coun-terparts to GWs, by studying the near-zone mechanisms thatcould drive EM emission. We investigate the role of the BHspins and of magnetic fields on the gas dynamics, exploringhow those parameters affect the rest-mass density evolution,as well as the velocity of the fluid in the vicinity of the binary.As a channel of EM emission we consider the Poyntingflux, which may provide a powerful supply of energy that canbe converted to strong EM emission farther from the BHs [70–72].To make contact with the results of Ke17, we evolve similarconfigurations (their canonical configuration is an equal-massbinary with d = 14 . M , ρ = 1 and β − = 0 . , in apolytropic gas with Γ = 4 / and κ = 0 . ). A. Gas Dynamics
Figures 1-2 show the evolution on the orbital plane xy andon the polar plane xz of the rest-mass density ρ (normalizedto its initial value ρ ) for the B2S3 ( β − = 0 . , a , = 0 . )and B0S3 ( β − = 0 , a , = 0 . ) configurations. We do notshow snapshots for S0, S6 cases since they qualitatively lookvery similar to S3 models.The evolution of the unmagnetized model B0S3 (Fig. 2) issimilar to B0 configuration (no magnetic fields, non-spinningBHs) in Gi12, with the production of two denser gas wakesduring the inspiral and the formation of a central spinning BH − − − x [M] − − − y [ M ] t = 306.67 [M] . . . . . . ρ / ρ − − − x [M] − − − y [ M ] t = 1863.33 [M] . . . . . . ρ / ρ − − − x [M] − − − y [ M ] t = 2500.00 [M] . . . . . ρ / ρ . − − − x [M] − z [ M ] t = 306.67 [M] . . . . . . ρ / ρ − − − x [M] − z [ M ] t = 1863.33 [M] . . . . . . ρ / ρ − − − x [M] − z [ M ] t = 2500.00 [M] . . . . . . ρ / ρ FIG. 2. Rest-mass density ρ (normalized to its initial value ρ ) on the xy (top panels) and xz planes (bottom panels) for the B0S3 configuration( a = 0 . , β − = 0 ). The snapshots were taken, respectively, after ∼ ∼ ∼ M after themerger. The regions inside the BH horizons have been masked out. after merger. Throughout the evolution, the two inspirallingBHs are surrounded by spherical overdensities of matter ac-creting onto the horizons; after merger, the final BH is ringedby an almost-isotropycal, high-density, spherical distributionof accreting matter (Fig. 2, right column).The magnetized models exhibit different features (for acomparison, see e.g. the B2 configuration in Gi12, Fig.1, andthe b1e-1 configuration in Ke17, Fig. 3 and 4). In all ourmagnetized simulations the density close to each BH and inthe regions connecting them is larger compared to the unmag-netized cases. We found that the rest-mass overdensities nearthe BHs in the B1 (B2) models are ∼
50% larger than those inthe B0 (B1) models; conversely, the individual BH spins showno effect on the enhancement of ρ . In the magnetized models,the regions close to the BHs reveal the presence of turbulencein the fluid which is absent in the unmagnetized configurations(see, e.g., the top panels in Fig. 1, which display snapshots ofthe rest-mass density on the orbital plane xy for the magne-tized model B2S3). Figures 3-4 highlight the differences inthe dynamical evolution of the accretion flows between theunmagnetized model B0S6 ( β − = 0 , a , = 0 . and mag-netized model B2S6 ( β − = 0 . , a , = 0 . . For eachconfiguration, Fig. 3 displays a 2-dimensional snapshot ofthe Mach number field M on the equatorial plane xy , takenapproximately 1 orbit prior to coalescence. In the unmag-netized case B0S6 (top panel), the fluid is mostly subsonic.The two spiral fronts of the shock waves travel at transonicspeed through the inspiral, and are present all the way downto merger. In the magnetized case B2S6 (bottom panel), theshock fronts are hardly visible. The motion of the fluid is more chaotic, and the gas speed in the regions close to the BHs issupersonic.Differences between the unmagnetized and magnetizedconfigurations are noticeable also post-merger. In Fig. 4 weshow 2-dimensional snapshots of the circularity parameter ω (Eq. (23)) for the B0S6 and B2S6 models. Both snapshotswere taken ∼ M after coalescence, and display the mag-nitudes of ω and the fluid velocity fields around remnant KerrBHs with spin parameter a (cid:39) . [73]. The shaded ar-eas mark the regions within the innermost stable circular orbit[66]. In the unmagnetized model, the accretion flow on theequatorial plane is nearly radial at distances ≥ M , and thecircularity parameter at r ISCO is < . . Conversely, in themagnetized case the fluid exhibits a higher degree of rotationin the xy plane, and the φ -averaged circularity at r ISCO is ∼ . . B. Evolution of Magnetic Fields
During the evolution of the B1 and B2 models, the initiallyweak magnetic fields are dragged along each BH, and soonbecome dynamically important in the polar regions close tothe horizons (see Fig. 5, left panel). The magnetic field linesare twisted and compressed, producing a magnification of themagnetic field strength. After the coalescence, the magneticfield strength in the polar regions surrounding of the remnantBH is amplified by a factor ∼ ; this amplification is ob-served in all magnetized configurations (in agreement withGi12), and is little-to-no sensitive to the initial β − parameter − x [M] − y [ M ] − M a c h N u m b e r ℳ − − x [M] − − y [ M ] − M a c h N u m b e r ℳ FIG. 3. Mach number field M ≡ v /c s on the equatorial plane forB0S6 (top) and B2S6 (bottom) models. Arrows denote velocity vec-tors. The snapshots were taken ∼ orbit before merger. The BHsinteriors have been masked out. and to the individual BH spins.In Fig. 5 we show the evolution of the magnetic-to-gaspressure ratio β − on the xz plane for the B2S3 model. Afteras short as one orbit (left panel), we see that the polar regionsclose to the individual horizons are magnetically dominated(i.e., they have larger values of β − than the initial condi-tions). After 8 orbits (central panel) these regions becomemore pronounced, and outline two vertical areas which aredepleted of gas. After merger, a magnetically dominated fun-nel is created around the spin axis of the remnant BH (rightpanel). Along this region, the magnetic field strength is in-creased by a factor ∼ , contributing considerably to thetotal pressure in the gas.The effect of the magnetically dominated regions on theplasma distribution is noticeable in Fig. 6, where we compare − − x [M] − − y [ M ] − − − − − − c i r c u l a r i t y ω − − x [M] − − y [ M ] − − − − − − c i r c u l a r i t y ω FIG. 4. Circularity parameter ω on the equatorial plane for B0S6(top) and B2S6 (bottom) models. Arrows denote velocity vectors.The snapshots were taken ∼ M after merger. The BHs interiorshave been masked out. The shaded areas denote the regions withinthe innermost stable circular orbit ( r < R ISCO ) for a Kerr BH withspin parameter a (cid:39) . (see, e.g., Eq. (2.21) in [66]). the rest-mass distributions on the xz plane after merger (atthe same time t = 2500 M ) for the B0S3 (unmagnetized) andB2S3 (magnetized) models: the evolution of B0S3 results ina spherical distribution of matter accreting onto the final BH,whereas the endpoint of B2S3 evolution is the formation of athin, “disk-like” structure around the BH (see also Fig. 1-2,right panels). C. Mass Accretion Rate
An important diagnostic of our simulations is the flux ofrest-mass across the horizons of each BH. To study the massaccretion rate onto the BH horizons, we use the
Outflow − −
10 0 10 20 x [M] z [ M ] t = 306.67 [M] 1010 −− − β − − −
10 0 10 20 x [M] z [ M ] t = 1863.33 [M] 1010 −− − β − − −
10 0 10 20 x [M] z [ M ] t = 2500.00 [M] 10 − − − β − FIG. 5. Evolution of the magnetic to gas pressure β − = b / p fluid on the xz plane for the B2S3 configuration. The snapshots were taken,respectively, after ∼ ∼ ∼ M after the merger. − z [ M ] − − − x [M] − z [ M ] t = 2500.00 [M] . . . . . . FIG. 6. Rest-mass density ρ (normalized to its initial value ρ ) onthe xz plane for the B0S3 (top panel) and B2S3 (bottom panel) con-figurations, at t ∼ M after the merger. The regions inside theBH horizons have been masked out. thorn [74], which computes the flow of rest-mass densityacross a given spherical surface (e.g., in our case, across eachapparent horizon). This quantity is calculated via ˙ M = − (cid:73) S √ γDv i dσ i , (34)where D ≡ ραu is the fluid density measured in the observerframe (i.e. ρW , where W is the Lorentz factor), and σ i is theordinary (flat) space directed area element of the surface en-closing the horizon. Figures 7-9 show the time evolution ofthe mass accretion rates ˙ M onto the BH horizons for each bi-nary system. We plot the evolution of ˙ M for the unmagnetized t [ M · hours] . . . . . . . ˙ M [ ρ − M · M (cid:12) / y r ] a = 0 . a = 0 . a = 0 . B0 FIG. 7. Accretion rate ˙ M in solar masses per year onto the blackhole horizons for the unmagnetized (B0) models. The dotted linesmark the merger times for the non-spinning (magenta), a , = 0 . (blue), a , = 0 . (green) configurations. (B0) runs, the β − = 0 . (B1) runs, and the β − = 0 . (B2) runs. In each plot (i.e., for each level of magnetization),we compare the values of ˙ M for the three different spin con-figurations. The vertical, dotted lines mark the time of coa-lescence for each spin configuration (as expected, the mergerof spinning BHs is delayed as a result of the hang-up mecha-nism [75], which delays or prompts the coalescence accordingto the sign of the spin-orbit coupling).The quantities in Fig. 7-9 are scaled from code to physicalunits as follows: since ˙ M generally scales as ρM (g cm − ),we multiply the rate in code units ˙ M c . u . by a factor G c − (g − cm s − ) to obtain the rate in cgs units as ˙ M cgs = 6 . × ˙ M c . u . ρ − M g s − (35)0 t [ M · hours] . . . . . . . ˙ M [ ρ − M · M (cid:12) / y r ] a = 0 . a = 0 . a = 0 . B1 FIG. 8. Accretion rate ˙ M in solar masses per year onto the blackhole horizons for the β − = 0 . (B1) models. The dotted linesmark the merger times for the non-spinning (magenta), a , = 0 . (blue), a , = 0 . (green) configurations. t [ M · hours] . . . . . . . ˙ M [ ρ − M · M (cid:12) / y r ] a = 0 . a = 0 . a = 0 . B2 FIG. 9. Accretion rate ˙ M in solar masses per year onto the blackhole horizons for the β − = 0 . (B2) models. The dotted linesmark the merger times for the non-spinning (magenta), a , = 0 . (blue), a , = 0 . (green) configurations. and the rate in solar masses per year as ˙ M M (cid:12) yr − = 4 . × − ˙ M c . u . ρ − M M (cid:12) yr − (36)For each level of magnetization, the estimates of ˙ M share anumber of common features:• for the B0 configurations, ˙ M shows ( i ) an early max-imum as the gas surrounding the binaries establishesa quasi-equilibrium flow with the orbital motion, fol- t − t merger [ M · hours] . . . . . . . ˙ M [ ρ − M · M (cid:2) / y r ] a = 0 . a = 0 . a = 0 . (cid:4567) (cid:16)(cid:20) (cid:3) (cid:4567) (cid:16)(cid:20) (cid:3) (cid:4567) (cid:16)(cid:20) (cid:3) (cid:29)(cid:3)(cid:37)(cid:20)(cid:3)(cid:80)(cid:82)(cid:71)(cid:72)(cid:79)(cid:86)(cid:29)(cid:3)(cid:37)(cid:21)(cid:3)(cid:80)(cid:82)(cid:71)(cid:72)(cid:79)(cid:86)(cid:29)(cid:3)(cid:37)(cid:20)(cid:3)(cid:80)(cid:82)(cid:71)(cid:72)(cid:79)(cid:86)(cid:29)(cid:3)(cid:37)(cid:21)(cid:3)(cid:80)(cid:82)(cid:71)(cid:72)(cid:79)(cid:86) FIG. 10. Post-merger accretion rate ˙ M in solar masses per year ontothe remnant black hole horizons for the β − = 0 . (B1, straightlines) and β − = 0 . (B2, dotted lines) models. The colored arrowsdenote transitions from lower to higher levels of magnetization β − for same-spin configurations. lowed by ( ii ) a steady growth, that reaches its peak atmerger ( iii ) . After the coalescence, the accretion ratessettle to constant values ( iv ) .• For the B1/B2 configurations, ˙ M shows the same initialtransient as B0 ( i ) , followed by a steep decrease ( ii ) ,after which it settles to quasi-constant values whichslowly decline prior to merger ( iii ) . Just before the co-alescence, the flows drop ( iv ) , and jump upon merger ( v ) as the apparent horizons join discontinuously.The accretion rates of the magnetized configurations are gen-erally smaller than in the unmagnetized cases by a factor ∼ ˙ M for themagnetized models B1 and B2. We focus on the post-mergeraccretion onto the remnant Kerr BHs. Same-colored lines de-note same-spin models, whereas straight (dotted) lines standfor B1 (B2) models. We see that a higher initial magnetization(B2 to B1) has a suppressing effect on ˙ M , which is reducedby ∼
27% for spin parameter a , = 0 , by ∼
20% for spinparameter a , = 0 . , and by ∼
13% for a , = 0 . .Conversely, for a given value of β − , we find that ˙ M isreduced (compared to the non-spinning case) by ∼
23% for a , = 0 . , and by ∼
48% for a , = 0 . . D. Poynting Luminosity
Several works [16, 67–69] have shown that the interactionof orbiting MBHs with ambient magnetic fields results in theconversion of some of the BH energy into EM energy in theform of collimated regions of Poynting flux. Such regionsmay generate relativistic outflows [70–72], and through a cas-cade of matter interaction yield strong EM emission. All our1 − −
10 0 10 20 x [M] z [ M ] t = 306.67 [M] 10 − − − − − L P o y n t − −
10 0 10 20 x [M] z [ M ] t = 1863.33 [M] 10 − − − − − L P o y n t − −
10 0 10 20 x [M] z [ M ] t = 2500.00 [M] 10 − − − − − L P o y n t FIG. 11. Evolution of the z component of the Poynting vector (code units) on the xz plane for the B2S3 configuration. The snapshots weretaken, respectively, after ∼ ∼ ∼ M after the merger. t [ M · hours] L P o y n / L ( @ R = M ) a = 0 . a = 0 . a = 0 . B1 t [ M · hours] L P o y n / L ( @ R = M ) a = 0 . a = 0 . a = 0 . B2 FIG. 12. L Poyn ( z components) in units of L ≡ . × ρ − M erg s − for the B1 ( β − = 0 . , left) and B2 ( β − = 0 . , right)models, extracted on a coordinate sphere of radius R = 30 M . The dotted lines mark the merger times for the non-spinning (magenta), a = 0 . (blue), a = 0 . (green) configurations. magnetized simulations develop strong flows of electromag-netic energy in the form of Poynting flux; the Poynting fluxluminosity can be computed as a surface integral across a two-sphere at a large distance (see Appendix A): L Poynt ≈ lim R →∞ R (cid:114) π S z (1 , (37)where S z (1 , is the dominant ( l, m ) = (1 , spherical modeof the Poynting vector (Eq. (24)). Following the evolutionof L Poynt helps us measuring the amount of potential emis-sion on time scales comparable to the merger time. To this ex-tent, we extract the luminosity on a coordinate sphere of radius R ext ; we set the extraction radius at R ext = 30 M as was donein Ke17 (in Gi12, extraction was carried out at R ext = 10 M ,but the initial binary separation was ∼
30% smaller than in oursimulations). This choice allows us to avoid spurious effects due to the orbital motion of the BHs.In Fig. 11 we show the evolution of the z component of thePoynting vector on the polar plane xz for the B2S3 configura-tion. As in the simulations of Gi12- Ke17, the Poynting fluxemission in our simulations is largely collimated and parallelto the orbital angular momentum and to the spin of the post-merger BH. In Fig. 12 we display the Poynting flux luminositycomputed for each of the 6 magnetized models. On the left,we show the B1 configurations, i.e. those with β − = 0 . ;on the right the B2 configurations, with β − = 0 . . The val-ues of L Poynt are in units of L ≡ . × ρ − M ergs − (see Appendix B). The values of L Poynt which we observeare consistent with the EM power generated by the Blandford-Znajek [70] mechanism [see, e.g., Eq. (4.50) in 76]: L BZ ∼ erg s − ( a ) (cid:18) M M (cid:12) (cid:19) (cid:18) B G (cid:19) (38)2The main difference between the B1 and B2 configurationsis the time at which the modes reach the extraction sphere at30M. This is what we expected: both configurations evolvea magnetic field which is initially “dynamically weak”, i.e.the inertia of the plasma is larger than the magnetic field en-ergy. The lower the value of β − (B1 configuration), thestronger the initial magnetic field B must become in orderto surmount the fluid pressure. The development of a strongermagnetic field requires more time, thus a lower β − impliesa longer time for launching a jet [71], which is in agreementwith our simulations.We find that the peak luminosity which is reached shortlyafter merger is sensitive to the BH spins, and is enhanced bya factor of ∼ ∼ a =0 . ( a = 0 . ) with respect to the non-spinning case. Thisintensification is a general feature, and does not depend onthe initial level of magnetization of the gas. Nevertheless, thequalitative behaviour of all models is very similar. V. SUMMARY AND CONCLUSIONS
To expand our understanding of the physical processeswhich arise in the vicinity of merging massive black hole bina-ries, we carry out GRMHD simulations of equal-mass, spin-ning MBHB mergers in hot, magnetized environments. Wehave for the first time investigated the role of individual BHspins in the evolution of magnetic fields and gas dynamics.We evolve a set of 9 simulations covering a range of initiallyuniform, moderately magnetized fluids with different initialmagnetic-to-gas pressure ratios. For each magnetization level,we study distinct spin configurations defined by adimensionalspin parameters a = (0 , . , . .Our results offer some insight on the role of spin and mag-netization in the magnetohydrodynamical properties of hot ac-cretion flows around merging MBHBs, and on the physicalmechanisms which may provide electromagnetic counterpartsto future LISA observations. We have shown that across theorbital evolution, the magnetic field can be distorted by themotion of the BHs and significantly increase its strength, de-veloping magnetically dominated structures in the polar re-gions above each BH, and ultimately producing a magneti-cally dominated funnel around the spin axis of the remnantBH. In general, the dynamics of a magnetized fluid is differ-ent than in the unmagnetized case, even if the fluid is initiallynot magnetically dominated. The accretion flow in magne-tized environments yields turbulent motion in the gas near theinspiralling BHs, eventually leading to the formation of a thin,disk-like structure rotating on the equatorial plane of the rem-nant Kerr BH. These results are consistent with previous sim-ulations of non-spinning binaries.We find that mass accretion rates onto BH horizons inmagnetized fluids are generally smaller than in unmagnetizedcases by a factor ∼ ∼ ∼ β − of the gas.Technical limitations in our analysis prevent more detailedpredictions. Our simulations do not account for the emissionof radiation by the plasma, which would determine the magni-tude of the accretion luminosity and the shape of the spectra.Therefore, the accretion flows that we describe lack any radia-tive mechanism, including cooling and feedback.In this paper we made an attempt to extract physicallyrelevant information by evolving our simulations in simplegaseous environments, which help us highlighting the effectsof spins and magnetization on the accretion flows and emittedPoynting luminosity. While this choice may be useful to iden-tify the subtle effects of different spins and degrees of mag-netization, it is not clear how well this simplistic environmentcan stand in for real accretion flows, which realistically pos-sess angular momentum support and carry dynamical effectsfrom radiation flows. Assessing these limitations motivate ourfuture work. We aim at extending our exploration of the pa-rameter space of merging MBHBs investigating less symmet-rical systems. We will consider binaries with spins which areanti-aligned with the orbital angular momentum, as well asgeneric misaligned spin configurations. Also, we will con-sider binary systems with unequal mass ratios, which are thenatural outcome of galaxy mergers in cosmological simula-tions (as shown, e.g., in [77]).These improvements will let us question the effects of sur-rounding (pre-merger) material in powering electromagneticcounterparts to gravitational wave events, and the effect of thegravitational recoil imparted to the newly formed MBH on theshock-heated gas along the MBH trajectory. ACKNOWLEDGMENTS
We thank Bernard Kelly for useful comments and sugges-tions. All simulations were performed on GALILEO andMARCONI machines at CINECA (Bologna, Italy). Some ofthe numerical calculations have been made possible through aCINECA-INFN agreement, providing access to resources onMARCONI (allocation INF20 teongrav). F.C. acknowledgesCINECA award under the ISCRA initiative, for the availabil-ity of HPC resources on GALILEO (ISCRA-C project num-ber HP10CP7PQ1, allocation IsC83 HIGHSPIN). M.C. andF.H. acknowledge funding from MIUR under the grant PRIN2017-MB8AEZ.
Appendix A: Relation between Poynting Vector and EM flux
In the main text, we calculate the Poynting emissionthrough the (1, 0) spherical harmonic S z (1 , of the z -3component of the Poynting vector. The quantity S z (1 , isclosely related to the EM luminosity computed in the pioneer-ing work by Palenzuela et al. (2010) [16], where the emittedluminosity is determined in terms of the outgoing Newman-Penrose radiative scalar Φ = F µν n µ n ∗ ν . The square of Φ is connected to the electromagnetic energy flux: the EM lumi-nosity is given by the integral L EM = dE EM dt = lim r →∞ (cid:73) r π | Φ | d Ω . (A1)The quantity | Φ | is proportional to the radial component ofthe Poynting vector. Assuming that Φ is calculated on a Kerrbackground using the Kinnersley tetrad, we have (see, e.g.,[78] and Eq. (5.13) in [79]) dE EM dt = lim r →∞ (cid:73) r T r d Ω = lim r →∞ (cid:73) r π | Φ | d Ω (A2)In the 3+1 formulation of space-time, the quantity T r may beexpressed as [80, 81] (see also Eq. (24) in the main text) T r = − α e rjk E j H k = 1 α S r , (A3)where e ijk = √ γ(cid:15) ijk is the Levi-Civita pseudo-tensor associ-ated to the spatial 3-metric γ . As r converges to the numericalradial coordinate at large distances, we have α → , and theemitted Poynting luminosity can thus be expressed as L Poynt ≡ lim r →∞ (cid:73) r S r d Ω = lim r →∞ √ πr S r (0 , (A4)where S r (0 , is the ( l, m ) = (0 , spherical mode of S r . Torelate this quantity to the dominant (1, 0) spherical harmonicof S z , we assume that the Poynting flux is dominated by emis-sion along the z -axis. Then, we can write S r ∼ S z cos θ (A5) The (0 , and (1 , spherical harmonics modes of S r arerelated by S r (0 , = S r (1 , √ θ (A6)Therefore, combining Eqs. (A4) and (A6), we find lim R →∞ R (cid:114) π S z (1 , (A7)which is formula (37) used in our study.The assumption (A5) is supported by previous simulationsof BHBs in gaseous and magnetized environments. For exam-ple, Ke17 compared the evolution of both S z (1 , and √ S r (0 , and found that, except that for some difference in the initialgauge relaxation, the two signals closely agree. Appendix B: Converting L Poynt from code to cgs units
In the GRMHD simulations presented here, the Poyntingluminosity scales as L Poynt = ρ M F ( t/M ; (cid:15) , ζ ) (B1)where (cid:15) is the initial specific internal energy, ζ ≡ u mag /u fluid the initial magnetic-to-fluid energy density ratioand F is a dimensionless function of time (for more details,see Sec. 3 in Ke17).Equation (B1) is in code units, where c = G = 1 . Toconvert this relation to cgs units, we need to multiply by afactor G /c ≈ . × − g − cm s − , and we obtain L Poynt ( t ) =1 . × − (cid:18) ρ − (cid:19) (cid:18) M (cid:19) × F ( t ; (cid:15) , ζ ) erg s − (B2)If we want to scale with our canonical density ρ = 10 − gcm − , and for a system of two BHs of M = M = 10 M (cid:12) (i.e., M (cid:39) . × g), we find L Poynt ( t ) = 2 . × ρ − M F ( t ; (cid:15) , ζ ) erg s − = L ρ − M F ( t ; (cid:15) , ζ ) erg s − (B3)where ρ − ≡ ρ / (10 − g cm − ) and M ≡ M/ (10 M (cid:12) ) .The quantity L is the normalization factor used in Sec. IV D,Fig. 12. [1] M. C. Begelman, R. D. Blandford, and M. J. Rees, Nature ,307 (1980).[2] J. Kormendy and L. C. Ho, Ann. Rev. Astron. Astrophys. ,511 (2013).[3] P. Amaro-Seoane et al. , arXiv e-prints , arXiv:1702.00786(2017), arXiv:1702.00786 [astro-ph.IM]. [4] J. E. Barnes and L. Hernquist, Annu. Rev. Astron. Astrophys. , 705 (1992).[5] J. E. Barnes and L. Hernquist, Astroph. J. , 115 (1996).[6] L. Mayer, S. Kazantzidis, P. Madau, M. Colpi, T. Quinn, andJ. Wadsley, Science , 1874 (2007), arXiv:0706.1562 [astro-ph]. [7] M. Dotti, M. Ruszkowski, L. Paredi, M. Colpi, M. Volonteri,and F. Haardt, Mon. Not. R. Astron. Soc. , 1640 (2009).[8] M. Dotti, A. Sesana, and R. Decarli, Adv. Astron. , 1 (2012).[9] D. Chapon, L. Mayer, and R. Teyssier, Mon. Not. R. Astron.Soc. , 3114 (2013), arXiv:1110.6086 [astro-ph.GA].[10] M. Colpi, Space Science Reviews , 189 (2014).[11] S. McGee, A. Sesana, and A. Vecchio, Nature Astronomy10.1038/s41550-019-0969-7 (2020).[12] Z. Haiman, B. Kocsis, K. Menou, Z. Lippai, and Z. Frei, Clas-sical Quantum Gravity , 094032 (2009).[13] Y. Tang, Z. Haiman, and A. MacFadyen, Mon. Not. R. Astron.Soc. , 2249 (2018).[14] D. B. Bowen, M. Campanelli, J. H. Krolik, V. Mewes, and S. C.Noble, Astrophys. J. , 42 (2017).[15] D. B. Bowen, V. Mewes, M. Campanelli, S. C. Noble, J. H.Krolik, and M. Zilh˜ao, Astrophys. J. , L17 (2018).[16] C. Palenzuela, L. Lehner, and S. L. Liebling, Science , 927(2010).[17] A. Khan, V. Paschalidis, M. Ruiz, and S. L. Shapiro, Phys. Rev.D , 044036 (2018), arXiv:1801.02624 [astro-ph.HE].[18] B. F. Schutz, Nature (London) , 310 (1986).[19] B. Kocsis, Z. Haiman, and K. Menou, Astrophys. J. , 870(2008).[20] N. Tamanini, C. Caprini, E. Barausse, A. Sesana, A. Klein, andA. Petiteau, J. Cosm. Astrop. Phys. , 002 (2016).[21] C. Roedig, J. H. Krolik, and M. C. Miller, Astrophys. J. ,115 (2014).[22] B. J. Kelly, J. G. Baker, Z. B. Etienne, B. Giacomazzo, andJ. Schnittman, Phys. Rev. D , 123003 (2017).[23] S. d’Ascoli, S. C. Noble, D. B. Bowen, M. Campanelli,J. H. Krolik, and V. Mewes, Astrophys. J. , 140 (2018),arXiv:1806.05697 [astro-ph.HE].[24] C. Yuan, K. Murase, B. T. Zhang, S. S. Kimura, andP. M´esz´aros, arXiv e-prints , arXiv:2101.05788 (2021),arXiv:2101.05788 [astro-ph.HE].[25] M. Milosavljevi´c and E. S. Phinney, Astrophys. J. , L93(2005).[26] S. C. Noble, B. C. Mundim, H. Nakano, J. H. Krolik, M. Cam-panelli, Y. Zlochower, and N. Yunes, Astrophys. J. , 51(2012).[27] D. J. D’Orazio, Z. Haiman, and A. MacFadyen, Mon. Not. R.Astron. Soc. , 2997 (2013).[28] B. D. Farris, P. Duffell, A. I. MacFadyen, and Z. Haiman, As-trophys. J. , 134 (2014).[29] B. D. Farris, P. Duffell, A. I. MacFadyen, and Z. Haiman, Mon.Not. R. Astron. Soc.: Letters , L80 (2015).[30] Y. Tang, A. MacFadyen, and Z. Haiman, Mon. Not. R. Astron.Soc. , 4258 (2017).[31] D. B. Bowen, V. Mewes, S. C. Noble, M. Avara, M. Campanelli,and J. H. Krolik, Astrophys. J. , 76 (2019).[32] B. D. Farris, R. Gold, V. Paschalidis, Z. B. Etienne, and S. L.Shapiro, Phys. Rev. Lett. , 221102 (2012).[33] S. Ichimaru, Astrophys. J. , 840 (1977).[34] R. Narayan and I. Yi, Astrophys. J. Lett. , L13 (1994),arXiv:astro-ph/9403052 [astro-ph].[35] B. D. Farris, Y. T. Liu, and S. L. Shapiro, Phys. Rev. D ,084008 (2010).[36] T. Bode, R. Haas, T. Bogdanovi´c, P. Laguna, and D. Shoemaker,Astrophys. J. , 1117 (2010).[37] T. Bode, T. Bogdanovi´c, R. Haas, J. Healy, P. Laguna, andD. Shoemaker, Astrophys. J. , 45 (2012).[38] B. Giacomazzo, J. G. Baker, M. C. Miller, C. S. Reynolds, andJ. R. van Meter, Astrophys. J. , L15 (2012).[39] B. Giacomazzo and L. Rezzolla, Classical Quantum Gravity , S235 (2007).[40] B. Giacomazzo, L. Rezzolla, and L. Baiotti, Phys. Rev. D ,044014 (2011).[41] S. C. Noble, C. F. Gammie, J. C. McKinney, and L. D. Zanna,Astrophys. J. , 626 (2006).[42] Z. B. Etienne, V. Paschalidis, R. Haas, P. M¨osta, and S. L.Shapiro, Classical Quantum Gravity , 175009 (2015).[43] A. Marconi, G. Risaliti, R. Gilli, L. K. Hunt, R. Maiolino,and M. Salvati, Mon. Not. R. Astron. Soc. , 169 (2004),arXiv:astro-ph/0311619 [astro-ph].[44] C. F. Gammie, S. L. Shapiro, and J. C. McKinney, Astrophys.J. , 312 (2004).[45] A. R. King, S. H. Lubow, G. I. Ogilvie, and J. E. Pringle, Mon.Not. R. Astron. Soc. , 49 (2005), arXiv:astro-ph/0507098[astro-ph].[46] E. Berti and M. Volonteri, Astroph. J. , 822 (2008),arXiv:0802.0025 [astro-ph].[47] M. Dotti, M. Colpi, S. Pallini, A. Perego, and M. Volonteri,Astrophys. J. , 68 (2013), arXiv:1211.4871 [astro-ph.CO].[48] A. Sesana, E. Barausse, M. Dotti, and E. M. Rossi, Astrophys.J. , 104 (2014), arXiv:1402.7088 [astro-ph.CO].[49] D. Izquierdo-Villalba, S. Bonoli, M. Dotti, A. Sesana, Y. Rosas-Guevara, and D. Spinoso, Mon. Not. R. Astron. Soc. , 4681(2020).[50] F. G. Lopez Armengol, L. Combi, M. Campanelli, S. C. No-ble, J. H. Krolik, D. B. Bowen, M. J. Avara, V. Mewes,and H. Nakano, arXiv e-prints , arXiv:2102.00243 (2021),arXiv:2102.00243 [astro-ph.HE].[51] V. Paschalidis, M. Ruiz, M. Ruiz, and R. Gold, arXiv e-prints ,arXiv:2102.06712 (2021), arXiv:2102.06712 [astro-ph.HE].[52] F. L¨offler et al. , Classical Quantum Gravity , 115001 (2012).[53] E. Schnetter, S. H. Hawley, and I. Hawke, Classical QuantumGravity , 1465 (2004).[54] S. Husa, I. Hinder, and C. Lechner, Comput. Phys. Commun. , 983 (2006).[55] D. Brown, P. Diener, O. Sarbach, E. Schnetter, and M. Tiglio,Phys. Rev. D , 044023 (2009).[56] T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor. Phys.Supp. , 1 (1987).[57] M. Shibata and T. Nakamura, Phys. Rev. D , 5428 (1995).[58] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D , 024007(1998).[59] J. R. van Meter, J. G. Baker, M. Koppitz, and D.-I. Choi, Phys.Rev. D , 124011 (2006).[60] J. M. Bowen and J. W. York, Phys. Rev. D , 2047 (1980).[61] M. Ansorg, B. Br¨ugmann, and W. Tichy, Phys. Rev. D ,064011 (2004).[62] T. W. Baumgarte and S. L. Shapiro, Numerical Relativity: Solv-ing Einstein’s Equations on the Computer (Cambridge Univer-sity Press, 2010).[63] L. Rezzolla and O. Zanotti,
Relativistic Hydrodynamics (Ox-ford University Press, 2013).[64] M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens, Phys.Rev. D , 024028 (2005).[65] M. Gedalin, Phys. Rev. E , 4354 (1993).[66] J. M. Bardeen, W. H. Press, and S. A. Teukolsky, Astrophys. J. , 347 (1972).[67] P. M¨osta, C. Palenzuela, L. Rezzolla, L. Lehner, S. Yoshida, andD. Pollney, Phys. Rev. D , 064017 (2010).[68] C. Palenzuela, T. Garrett, L. Lehner, and S. L. Liebling, Phys.Rev. D , 044045 (2010).[69] P. M¨osta, D. Alic, L. Rezzolla, O. Zanotti, and C. Palenzuela,Astrophys. J. , L32 (2012).[70] R. D. Blandford and R. L. Znajek, Mon. Not. R. Astron. Soc. , 433 (1977).[71] V. Paschalidis, M. Ruiz, and S. L. Shapiro, Astrophys. J. ,L14 (2015).[72] M. Ruiz, R. N. Lang, V. Paschalidis, and S. L. Shapiro, Astro-phys. J. , L6 (2016).[73] M. Colpi and A. Sesana, Gravitational Wave Sources in the Eraof Multi-Band Gravitational Wave Astronomy, in An Overviewof Gravitational Waves: Theory (2017) pp. 43–140.[74] R. Haas, Cactus code thorn outflow, https://bitbucket.org/einsteintoolkit/einsteinanalysis/src/master/Outflow/ (2009).[75] M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys. Rev. D , 084023 (2006).[76] K. S. Thorne, R. H. Price, and D. A. MacDonald, Black holes:The membrane paradigm (1986).[77] M. Volonteri, H. Pfister, R. S. Beckmann, Y. Dubois, M. Colpi,C. J. Conselice, M. Dotti, G. Martin, R. Jackson, K. Kraljic,C. Pichon, M. Trebitsch, S. K. Yi, J. Devriendt, and S. Peirani,Mon. Not. R. Astron. Soc. , 2219 (2020).[78] S. A. Teukolsky, Phys. Rev. Lett. , 1114 (1972).[79] S. A. Teukolsky, Astrophys. J. , 635 (1973).[80] S. S. Komissarov, Mon. Not. R. Astron. Soc. , 427 (2004).[81] D. Alic, P. M¨osta, L. Rezzolla, O. Zanotti, and J. L. Jaramillo,Astrophys. J.754