GRMHD study of accreting massive black hole binaries in astrophysical environment: a review

We present an overview of recent numerical advances in the theoretical characterization of massive binary black hole (MBBH) mergers in astrophysical environments. These systems are among the loudest sources of gravitational waves (GWs) in the universe and particularly promising candidates for multimessenger astronomy. Coincident detection of GWs and electromagnetic (EM) signals from merging MBBHs is at the frontier of contemporary astrophysics. One major challenge in observational efforts searching for these systems is the scarcity of strong predictions for EM signals arising before, during, and after merger. Therefore, a great effort in theoretical work to-date has been to characterize EM counterparts emerging from MBBHs concurrently to the GW signal, aiming to determine distinctive observational features that will guide and assist EM observations. To produce sharp EM predictions of MBBH mergers it is key to model the binary inspiral down to coalescence in a full general relativistic fashion by solving Einstein's field equations coupled with the magnetohydrodynamics equations that govern the evolution of the accreting plasma in strong-gravity. We review the general relativistic numerical investigations that have explored the astrophysical manifestations of MBBH mergers in different environments and focused on predicting potentially observable smoking-gun EM signatures that accompany the gravitational signal.


INTRODUCTION
The study of massive binary black hole (MBBH) evolution dates back as far as the early eighties.In their pioneering work, Begelman, Blandford, & Rees (1980) explored the possibility that following a galaxy coalescence event, the active nucleus of the post-merger galaxy may harbor two massive black holes (MBHs) in orbit about each other.
Once two galaxies merge, the MBHs in their nuclei sink to the center of the merger remnant via a number of mechanisms such as dynamical friction via interaction with surrounding stars, individual interactions with singular stars, and torques exerted by a circumbinary gaseous disk (see, e.g., Dotti et al. 2012, for a review).Following these interactions, the binary will enter the so-called gravitational wave (GW) regime where it will evolve until merger.Gravitational waves are perturbations of the geometry of spacetime produced by the acceleration of massive objects (Dirkes 2018), which propagate away from the sources at the speed of light.
Since the first detection of GWs by the Laser Interferometer Gravitational-Wave Observatory (LIGO)/Virgo collaboration in 2015 (Abbott et al. 2016), GW observations have opened a new way to observe and characterize binary systems of compact objects, such as (stellar-mass) black hole binaries (Abbott et al. 2021a,b), neutron-star binaries (Abbott et al. 2017), and black hole-neutron star binaries (Abbott et al. 2021c).
The binary evolution in the GW regime (Fig. 1) is commonly divided into three phases, i.e. inspiral, merger, and ringdown (for further knowledge on the features of propagating gravitational radiation and GW sources we defer the reader to Maggiore (2007)).The classification into three (approximately defined) evolutionary stages is motivated by the different techniques used to describe each phase: (i) The inspiral is the most protracted stage of the GWdriven orbital evolution, in which the binary is wellseparated and the black holes' inward radial migration timescale is much longer than their orbital timescale t orb .This stage can be treated analytically and is modelled using Post-Newtonian (PN) techniques.
(ii) The second stage of the GW-driven binary evolution includes the late stages of the inspiral and the merger Figure 1.The three phases in the temporal evolution of a binary system.In the inspiral phase (i), the black holes' inward radial migration timescale is much longer than their orbital timescale t orb .In the merger phase (ii), the two objects merge into one and the GW signal reaches its peak.In the ringdown phase (iii), the resulting object relaxes to a stationary Kerr configuration by emitting GWs at characteristic complex frequencies.
(or coalescence), in which the black holes plunge together and merge, while the GW emission reaches its peak.This phase (in which strong field, fully nonlinear gravity plays a paramount role) must be treated numerically.
(iii) The third phase is ringdown, in which the merger remnant settles down into a Kerr configuration by emitting GWs at characteristic complex frequencies.It can be modelled analytically with perturbation methods.
Exploring the second stage of the GW-driven binary evolution requires solving on a computer the full Einstein's equations of general relativity in which we assumed the geometric unit system for which c = G = 1.This is a formidable task and took numerical relativists several decades and a great deal of trials and errors to overcome.In the past two decades, a number of remarkable breakthroughs has yielded robust and accurate full GR simulations of inspiralling and merging binary black holes, contributing to lay the foundations of gravitational wave astrophysics (Pretorius 2005;Baker et al. 2006;Campanelli et al. 2006a,b;Berti et al. 2007).More recently, an increasing number of general relativistic magnetohydrodynamic (GRMHD) numerical studies have begun to explore the strong-field dynamics of the gas in the vicinity of massive binary black holes, investigating the mechanisms that might potentially drive electromagnetic (EM) emission concurrent to the gravitational radiation (Palenzuela et al. 2010b;Noble et al. 2012;Farris et al. 2010Farris et al. , 2012;;Gold et al. 2014a,b;Giacomazzo et al. 2012;Kelly et al. 2017;Cattorini et al. 2021).
We direct the reader to Gold (2019) for a review of theoretical works on accreting black holes in the strong gravitational regime (see also Schnittman 2013).Gold (2019) presents an overview of the theoretical predictions on these systems, focusing on results from simulations of black hole binaries surrounded by a circumbinary accretion disk.Also, the recent review by Bogdanović et al. (2022) provides a comprehensive introduction to the current body of knowledge about the nature of EM counterparts to MBBH mergers, both on the observational and theoretical sides.
The present work offers a complementary revision of the subject, presenting the current state-of-the-art of GRMHD simulations and discussing the diverse environments in which a massive binary may be embedded in the late stages of the GW-driven inspiral.

NUMERICAL RELATIVITY SIMULATIONS OF BINARY BLACK HOLES
Most numerical relativity simulations start with decomposing the four-dimensional spacetime into a set of non-intersecting three-dimensional spatial hypersurfaces parametrized by a time function t.In the 1960s, Arnowitt, Deser, & Misner (1962) pioneered this "3+1" approach with their Hamiltonian formulation of GR.During the following years, a number of attempts to numerically evolve the head-on collision of two equal-mass black holes using the ADM formalism were carried out (Eppley 1975;Smarr 1975;Smarr et al. 1976;Smarr 1977).These simulations were twodimensional and employed an axisymmetric approach, managing to evolve the systems to the collision and to extract information about the emitted gravitational signal.To this aim, they employed the most powerful computers in their day (as an example, the code presented in Smarr (1975) ran on the UT-CDC system, a mainframe calculator manufactured by Control Data Corporation with a 60 bit central processor performing at 10 MHz).Unfortunately, the next steps -threedimensional simulations of orbiting binaries -were not feasible at the time, due to insufficient computational resources and numerical instabilities, e.g., the "grid-stretching" problem, which leads to divergent metric coefficients (Brügmann 1999).Consequently, the numerical exploration of merging binary black holes laid dormant for years.
The tides turned during the late 1990s and into the early 2000s, when a number of important results were accomplished.These key developments include: new grid-based methods for representing black holes such as punctures (Brandt & Brügmann 1997) and excision (Seidel & Suen 1992); recognition of the importance of the hyperbolicity in formulating the Einstein's equations (Bona & Massó 1992;Abrahams et al. 1997); the improved BSSNOK formulation of the Einstein's equations (Nakamura et al. 1987;Shibata & Nakamura 1995;Baumgarte & Shapiro 1998); singularityavoiding coordinate conditions (Bona et al. 1995;Alcubierre et al. 2003); the Cactus computational framework (Goodale et al. 2003), including several codes apt to the numerical evolution of general relativistic systems.
Eventually, in 2005, the pioneering work by Pretorius (2005) ushered in the first stable numerical evolution of merging BBHs in vacuum.Later in 2005, the groups at University of Texas at Brownville (UTB) and NASA's Goddard Space Flight Center independently developed a new technique, called "moving punctures", that also produced successful simulations of binary black hole mergers (Baker et al. 2006;Campanelli et al. 2006a).Since then, outstanding progress has been made to explore numerically the late stages of binary evolution and it is now possible to perform stable and accurate general relativistic simulations of BBHs in vacuum exploring a wide range of parameter space (Healy & Lousto 2022), including mass ratios1 as small as q = 1/128 (Lousto & Healy 2020), and nearly extremal spin parameters χ = S i /m2 i = 0.99 (Lovelace et al. 2010;Zlochower et al. 2017).
After the initial breakthrough with equal-mass, nonspinning black holes, the remarkably robust "moving puncture" method was soon applied to a wide variety of systems, including unequal masses (Berti et al. 2007), eccentric orbits (Hinder et al. 2008), and spinning BHs (Campanelli et al. 2006b), and generic black hole binaries, i.e., binaries containing unequal-mass black holes with misaligned spins (Campanelli et al. 2007).

MAGNETOHYDRODYNAMIC SIMULATIONS OF MBBHS
While numerical relativity had initially been devoted to applications to the direct detection of GWs, much of the recent theoretical work has focused on predicting potentially observable electromagnetic (EM) signals emerging from MBBHs.In fact, since galaxy mergers are usually followed by gas inflows toward the center of the merger remnant, the gravitational signal emitted by MBBH mergers is expected to be accompanied by an EM counterpart generated by the gas accreting onto the merging binary (Bogdanović et al. 2022;Amaro-Seoane et al. 2022).
One step toward the general relativistic modelling of the accretion flows onto merging MBBHs was taken by van Meter et al. (2010).This work produced three-dimensional general relativistic simulations of merging binary black holes and modelled the environment around the binary as a flow of pressureless matter made up by non-interacting point particles (i.e., particles of non-zero rest mass in geodesic motion) evolved with particle tracing techniques.While pressurelessmatter simulations can be computationally efficient and offer certain insights, they lack the physical depth required to explore realistic accretion flows, which involve the interaction of gas, radiation, and magnetic fields.The exploratory simulations by van Meter et al. (2010) lacked the required sophistication needed to model the gaseous environment in the vicinity of merging binaries.In actual accretion flows onto black holes, gas often experiences compression and heating as it spirals inward, leading to the formation of shocks.Hydrodynamic (HD) simulations excel in capturing shock physics, including the dissipation of kinetic energy into thermal energy.Conversely, pressureless-matter simulations are unable to accurately model these shock formations and the associated energy dissipation processes.Furthermore, the realistic portrayal of accretion flows necessitates the consideration of the interaction between gas and magnetic fields.This intricate coupling can be effectively addressed through magnetohydrodynamic (MHD) calculations.
At present, our understanding of the features of these accretion flows is in its infancy.The structure of the accretion flows around coalescing MBBHs largely depends on the angular momentum content of the accreting gas conveyed in the galactic merger and on its thermodynamical state.Two limiting physically motivated scenarios bracket the range of properties of accreting fluids around MBBHs, i.e. the circumbinary disk and the gas cloud models (Fig. 2).In the following we outline the properties of both scenarios and describe the results from recent state-of-the-art, general relativistic numerical simulations.
In the following, we will often express distances in units of the total mass of the system M.In GRMHD simulations (as long as that the total mass of the fluid is negligible with respect to the mass of the binary), one can neglect source terms in the spacetime evolution equations.If this is the case, the total mass of the system M sets2 the length unit l = GM/c 2 and the time unit t = GM/c 3 .In geometric units, l = t = M.

Circumbinary disk scenario
Most theoretical studies in the past decade have focused on the circumbinary disk (CBD) model.If, at the point of formation of the gravitationally bound binary, the two MBHs have a mass ratio larger than a few percents, theoretical studies indicate that the binary torques can clear a central low density cavity with a radius corresponding to approximately twice the binary semimajor axis (see, e.g., D 'Orazio et al. 2013'Orazio et al. , 2016)).Simulations also indicate that despite strong binary torques, accretion into the central cavity continues unhindered relative to the single MBH case (e.g., D 'Orazio et al. 2013;Farris et al. 2014;Shi & Krolik 2015).
Circumbinary black hole accretion simulations have been conducted adopting diverse methods (see Gold 2019, for a recent review of circumbinary accretion simulations incorporating relativistic effects).In the following, we will focus on numerical relativity simulations that solve the GRMHD equations over dynamical spacetimes that are evolved either by solving the full set of Einstein's equations or with approximate techniques such as the so called Post-Newtonian (PN) gravity, which describes the spacetime by a highorder PN approximation.One advantage of evolving threedimensional GRMHD simulations with PN gravity is that they do not carry the burden of evolving the spacetime by solving the full set of Einstein's equations.Therefore, a larger amount of computational resources may be dedicated to the MHD study of the CBD regions (Noble et al. 2012;Bowen et al. 2018;Noble et al. 2021), including more realistic physical assumptions and covering sufficiently wide time scales in order to account for accurate variability analysis.
Conversely, a full general relativistic treatment of the spacetime evolution (Farris et al. 2012;Gold et al. 2014a,b;Giacomazzo et al. 2012;Paschalidis et al. 2021) is crucial to resolve the physics in the proximity of the BH horizons (e.g., the launching of jets from the interaction of magnetic fields with the horizons).

POST-NEWTONIAN MHD SIMULATIONS
In the last decade, a variety of numerical groups have developed a number of high-end MHD simulations of circumbinary accretion onto MBBHs employing diverse successful approximations for the spacetime metric evolution.The pio-neering investigation by Noble et al. (2012) presented a simulation of circumbinary accretion onto an equal-mass binary, evolving the GRMHD equations with the Harm3D code (Noble et al. 2009) over a changing spacetime described by a high order Post-Newtonian (PN) approximation.This simulation provided a firm step toward the general relativistic investigation of a CBD surrounding a MBBH and established reasonable prior conditions for the gas that feeds these systems.To produce a quasisteady disk structure, Noble et al. (2012) first ran a "secularly-evolving" configuration (RunSE) that kept the binary semimajor axis artificially fixed at a separation R = 20M, followed by an "inspiral" run (RunIn) that evolved the binary down to a separation of 8M.
The RunSE evolution exploited the helical Killing symmetry (Klein 2004), i.e. a particular symmetry of the spacetime that describes a binary held at fixed separation.This symmetry ensures that, while the separation is fixed, the spacetime is invariant in a frame corotating with the binary.Within such an asymmetric spacetime, the circumbinary disk is initialized by an axisymmetric gas distribution supported by pressure and rotation, and in hydrostatic equilibrium with a specific angular momentum profile (see also Chakrabarti 1985a,b;Villiers et al. 2003).In their work, Noble et al. (2012) followed the evolution of the gas surface density and found that torques driven by the binary can drive matter toward the inner binary domain (this central region was excised from the computations to avoid evolving the gas in the vicinity of the black holes).Also, they found that a "lump" of matter (Shi et al. 2012) forms in the disk, whose density contrast with respect to the adjacent regions grows steadily in time (see Fig. 3).The subsequent investigation by Noble et al. (2021) presented a suite of 3D GRMHD simulations with the same PN spacetime as model "RunSE" by Noble et al. (2012) and carried out a parameter exploration surveying different mass ratios and mass/magnetic fluxes in the accretion disks, aiming to linking the properties of the CBDs to the GW signal.
The circumbinary accretion and disk-feeding mechanisms were extensively explored by Bowen et al. (2018Bowen et al. ( , 2019)), which presented MHD simulation in which a CBD around a binary feeds mass to individual accretion disks ("minidisks") around each MBH, and found that the streams of matter feeding the mini-disks comes into phase with the lump at the inner edge of the circumbinary disk.The simulation presented in these works evolve an equal-mass binary for approximately 12 orbital periods starting at an initial separation of 20M, thus allowing for a detailed examination of the asymmetric accretion flow and streams which deposit matter onto the mini-disks.The authors found that a number of aspects of the binary late inspiral can be characterized by quasiperiodic oscillations associated with the interaction of the lump with the individual MBHs.For instance, they identify quasiperiodic oscillations in the mass-flux of the accre- tion stream, passing (or "sloshing") from near one MBH to near the other and back again.Bowen et al. (2018) suggested that the quasiperiodic modulations in the mini-disks structure may be related to characteristic time-dependent signatures in the EM emission of the binary.
To make predictions on the spectrum and the time dependence of the EM signal, D'Ascoli et al. ( 2018) post-processed the data reported in Bowen et al. (2018) with the ray-tracing code BOTHROS (Noble et al. 2007) adopting the fast-light approximation, in which the spacetime is frozen and light propagates on the fixed background metric.Their predictions are based on the assumption that (i) where the gas is optically thick (in the system's photosphere), it radiates a local thermal spectrum and (ii) where the gas is optically thin (mostly on the top and bottom surfaces of the disks), it radiates hard X-rays (in a manner similar to AGN).The spectra produced using BOTHROS can then be described in terms of two components, i.e. a thermal UV/soft X-ray portion and a coronal hard X-ray spectrum.The CBD, the accretion streams, and the mini-disks contribute to both spectral components (see Fig. 4); still, the emitted power is dominated by the ther-mal UV (with only ∼ 1% radiated in the hard X-rays).Also, the majority of the flux was found to be originating from the CBD rather than from the mini-disks (∼65% and ∼25%, respectively).Notably, the spectral energy distribution in Fig. 4 displays features consistent with predictions by the model of Tanaka et al. (2012), which addresses the dynamical state and thermal emission features of accretion flow onto a MBBH in a CBD.D'Ascoli et al. ( 2018) observe that, depending on the observer's viewing angle, modulations in the hard X-rays fluxwith frequencies comparable to the orbital frequency -might arise due to a number of causes (i.e., Doppler beaming, gravitational lensing, periodic interaction of individual BHs with the lump).
There are also other approaches that may be adopted to evolve the relativistic MHD equations over an approximate background spacetime.Mundim et al. (2014) constructed a global, approximate metric that describes a nonspinning, equal-mass binary in quasicircular orbit by asymptotically matching different approximate solutions to Einstein's equations (see also Bowen et al. 2018).A complementary treatment was developed by Ireland et al. (2016), who developed an approximate spacetime which describes the inspiral of nonprecessing, spinning BBHs by matching BH perturbation theory to PN formalism.More recently, Lopez Armengol et al. ( 2021) introduced a new approach for evolving binaries of spinning BHs, building the approximate metric as a linear superposition of two boosted Kerr-Schild BHs.This approximate spacetime (called superimposed Kerr-Schild, or SKS) was employed to produce longterm GRMHD circumbinary simulations of aligned-spinning binary black holes orbiting at fixed separation, in which the CBD is initialized by a Fishbone-Moncrief solution (Fishbone & Moncrief 1976) and evolved solving the GRMHD equations over the background SKS metric.This system is then evolved for 266 binary periods (corresponding to approximately t = 150 × 10 3 M), allowing for the investigation of the role of spin on the circumbinary accretion and luminosity.Specifically, spins anti-aligned (aligned) with the orbital angular momentum of the binary result in an enhanced (reduced) mass accretion rate with respect to the nonspinning case.(called superposed harmonic PN, or SHPN) is valid at every point of the computational grid, hence dispensing with the need for artificial sink terms or excision.The SHPN metric was adopted by Combi et al. (2022) in a GRMHD evolution of a mini-disk circumbinary accretion onto alignedspinning BHs.The results of this simulation provided insight on the role that spin plays in the accretion rate, inflow time, and periodicities, as well as in the structure of the minidisks.Data from the spinning MBBHs simulation performed by Combi et al. (2022) and from the nonspinning configurations by Bowen et al. (2018Bowen et al. ( , 2019) ) have been recently postprocessed by Gutiérrez et al. (2022), who extended the work done by D' Ascoli et al. (2018) to much longer simulations.
The longer duration allowed the system to relax and, consequently, to identify periodic features in the EM emission after the system has reached a quasi-steady state.These periodicities were found to be associated with the lump's dynamics, and were predicted to be present at various wavelengths.
These recent investigations pose major progress in the study of more realistic accretion flows onto massive black hole binaries, as MBHs are expected to have nonzero spin.In fact, there is observational evidence suggesting that MBHs have grown primarily by efficient accretion (Marconi et al. 2004) and are expected to acquire a non-vanishing spin depending on whether accretion is prograde or retrograde, coherent or chaotic, as discussed extensively in the literature (see, e.g., Refs.Gammie et al. 2004;King et al. 2005;Berti & Volonteri 2008;Dotti et al. 2013;Sesana et al. 2014;Izquierdo-Villalba et al. 2020).As MBBHs approach merger, spins can have a major impact on the evolution of spacetime: they can alter the orbital motion, induce precession and nutation, tilt the orbital orientation and flip their sign (see, e.g., Campanelli et al. 2006b;Campanelli et al. 2007;Healy & Lousto 2018;Gangardt et al. 2021).They have a key role in the mass accretion onto BHs and BBHs (Krolik et al. 2005).Also, they are an essential ingredient for the production of relativistic outflows (Blandford & Znajek 1977).
Investigating the role of spin in determining the gas dynamics and EM emission during the late-inspiral and merger stages of massive black hole binary evolution requires a fully general relativistic treatment of the spacetime evolution.In the following section, we give an overview of GRMHD simulations of circumbinary accretion onto MBBHs that employ fully general relativistic calculations for the evolution of the spacetime metric.

FULL GENERAL RELATIVISTIC MHD SIMULATIONS
The studies discussed in the previous section have produced GRMHD simulations of circumbinary accretions flows onto binaries that are evolved by Post-Newtonian gravitation, often excluding the inner region in the proximity of the BHs and stopping the evolution at a binary separation where the PN approximation becomes inaccurate (e.g., the RunIn configuration by Noble et al. 2012, evolves the binary down to r = 8M).To evolve a black hole binary beyond this regime, one needs to carry out fully non-linear general relativistic simulations, using numerical relativity techniques to evolve the spacetime metric by solving Einstein's equations with no approximation.This is the case of the innovative work carried out by Farris et al. (2011Farris et al. ( , 2012)), which presented full general relativistic hydrodynamic (GRHD) and GRMHD simulations of equal-mass, nonspinning binary black hole mergers in a circumbinary disk.These simulations evolved both the pre-decoupling and post-decoupling phases of a binary-disk system.
Throughout the pre-decoupling stage (or "early-inspiral" epoch), the inspiral time scale is much longer than the orbital time scale.This fact can be exploited by neglecting the shrinking of the binary radius; to this extent, the binary is initialized adopting conformal thin-sandwich (CTS) initial conditions for a binary in quasicircular orbit and evolving the MHD equations in a metric that is quasistationary in the rotating frame of the binary (analogously to the RunSE simulation by Noble et al. 2012).Such metric is evolved employing CTS lapse and shift functions via a simple coordinate rotation, allowing for an accurate solution of the evolving spacetime without the computational burden of a full evolution of the Einstein's equations.This first stage of the binary evolution allows for the relaxation of the CBD and provides realistic initial data, employed to start the evolution of the late-inspiral.The disk is initialized by an equi-librium solution for a stationary disk around a single Kerr BH (Chakrabarti 1985a).Having allowed for the relaxation of the disk to a quasistationary state, the simulation of lateinspiral and merger stages is then produced by solving the Einstein's field equations with the BSSNOK formalism with moving puncture gauge conditions (Campanelli et al. 2006a;Baker et al. 2006).
The work by Farris et al. (2012) presented GRMHD simulations of circumbinary accretion onto BBHs that account both for the dynamical spacetime and the BH horizons and laid the groundwork for subsequent studies.In 2014, the numerical investigation presented in Gold et al. (2014a,b) extended this work by producing full GRMHD simulations of circumbinary accretion onto unequal-mass BBHs (with mass-ratio q ≡ m 2 /m 1 ∈ [0.1, 1], where m 2 < m 1 .).This suite of simulations considered both the pre-decoupling phase (Gold et al. 2014a) and the late-inspiral and merger phase (Gold et al. 2014b).Solving the full set of Einstein's field equations along with the GRMHD equations allowed Gold et al. (2014b) to explore purely relativistic features of MBBHs, such as the formation and evolution of jets resulting from twisting the magnetic field lines anchored in the plasma, thus highlighting the importance of full GR in discovering phenomena that may be entirely missed in Newtonian and Post-Newtonian studies.The subsequent investigation by Khan et al. (2018) tested the sensitivity of the GRMHD features of circumbinary accretion flows to the initial disk model, carrying out simulations of unequal-mass (q = 29/36) binaries in magnetized accretion disk and surveying disk configurations characterized by different scale heights, physical extent, and magnetic-to-gas-pressure ratios.This work has shown that the general relativistic features of magnetized accretion discovered in Gold et al. (2014a,b) are robust and little sensitive to the choice of the initial disk model.
Another feature of circumbinary accretion in which GR plays a key role is the formation and stability of mini-disks.The simulations by Bowen et al. (2018Bowen et al. ( , 2019) ) are initialized with mini-disks (set up by quasi-hydrostationary torii around each BH) in addition to a CBD, but do not report steadystate mini-disks (see previous section).Gold et al. (2014b) suggested that mini-disks could form whenever there are stable circular orbits within the Hill sphere 3 around each BH.This hypothesis was tested by Paschalidis et al. (2021), who produced full GRMHD simulations of accreting binaries of equal-mass, spinning black holes and focused on the role of 3 The Hill spheres are regions of the spacetime where gravity is dominated by each of the binary components and are defined by the Hill radius r ± Hill = 0.5 M ± /3M ∓ 1/3 a, where M ± is the mass of the primary (secondary) BH, and a is the orbital separation.spin in the formation and stability of mini-disks.This work evolved four binary configurations with black holes either non-spinning or with spin aligned (antialigned) with the orbital angular momentum.To study the formation and evolution of mini-disks, the authors track the rest-mass accretion rates, the density profiles, and the mass contained within the Hill spheres of each BH.They observe that black holes with spin aligned (antialigned) to the orbital angular momentum contain more (less) mass within the Hill sphere and accrete less (more) matter.This finding is consistent with the expectation that the size of the innermost-stable-circular-orbit (ISCO) radius has a significant impact on the ability of a black hole to maintain mass within the Hill radius and form mini-disks4 .
The accretion streams from the CBD form mini-disks whenever there are stable circular orbits around each black hole's Hill sphere; thus, persistent mini-disks are expected to be present for large radial separations between the BHs, corresponding to Hill radii larger than the ISCO radius.As the binary shrinks due to the emission of GWs, the separation between the BHs will reach a threshold at r Hill ∼ r ISCO beyond which no stable mini-disk is able to exist (as the Hill radius decreases linearly with the binary separation).Since the spin of a BH influences the radius of the ISCO, spin is found to play a significant role in whether mini-disks can form at rela-tivistic orbital separations.This work also has important implications for future GRMHD simulations of circumbinary systems, as it establishes that the initial binary separation should respect the criterion r Hill > r ISCO to allow for the formation of mini-disks.In addition to the mini-disk dynamics, Paschalidis et al. (2021) also found that the Poynting luminosity associated with the collimated jet outflows launched by the binary systems is significantly larger when spinning BHs are involved.
The simulations presented in Paschalidis et al. (2021) were extended in Bright & Paschalidis (2023), in which the structure of the mini-disks is examined in greater detail.We display in Fig. 5 the rest mass density in the equatorial plane in the aligned-spin model χ ++ from the simulations by Bright & Paschalidis (2023), which exhibits clear persistent minidisks around each of the black holes.This work also explored the periodic nature of the accretion rate and mass contained within the Hill spheres around each black hole, and found that the presence of persistent mini-disks is not only associated with lower accretion rates, but also with a suppression in the strength of the quasiperiodic modulation in the accretion rate.This result suggests that binaries at large separations may have dampened modulations, and it is not clear whether they will be able to exhibit any observable periodicity.Thus, quasiperiodic behavior in the observed EM signal arising from modulations in the mass accretion rate may not be a smoking gun for the existence of a binary.The recent numerical investigation by Ruiz et al. (2023) has provided a comprehensive exploration of the parameter space characterizing circumbinary accretion presenting full GRMHD simulations of magnetized accretion disks onto unequal-mass binaries of aligned-and misaligned-spinning MBBHs.This work investigates the impact of the binary mass-ratio and the black hole spins on the EM signatures that may emerge from the system.Considering different spin configurations (which may lie either along the initial orbital plane or an angle θ = π/4 above it) allows to explore the impact of spin on the direction of the jet-like emission that is launched from the individual inspiraling black holes and the remnant.The results of Ruiz et al. (2023) show that the reorientation of the remnant's spin axis drives a sudden change of the magnetic-field-lines orientation, which, however, is dampened in ∆t ≤ 15 M.This result indicates that spin-realignment at merger is not sufficient to induce helical distortions of the EM jet produced by spinning black holes.The production and orientation of jets has also been studied in the recent work by Cattorini et al. (2023), which presents GRMHD calculations of merging unequalmass binaries of misaligned-spinning black holes embedded in clouds of magnetized matter.This study shows that jets are oriented along the direction of the black hole spin axis (both across the inspiral and after merger) for distances ≤ 5 M; at larger separations, the jets lose memory of the spin di- rections and are orientated along the z-axis, corresponding to the initial orientation of the magnetic field.In Fig. 6, we show a comparison between the general relativistic studies by Ruiz et al. (2023) and Cattorini et al. (2023), displaying two-dimensional (meridional) slices of the magneticto-fluid energy density (or pressure) parameter for unequalmass, misaligned-spinning binary configurations.
Despite considering different environments, the results presented in these studies are consistent.

Gas cloud scenario
The balance of heating and cooling in the accretion flow in the vicinity of a MBBH can significantly alter towards the merger, when the gas is expected to be shock-heated and permeated by energetic radiation.If radiative cooling of the accretion flow around the MBBH is inefficient, it would resemble a hot and tenuous radiatively inefficient accretion flow (RIAF Ichimaru 1977;Narayan & Yi 1994).A RIAF tends to be hot, optically thin and geometrically thick (h/r ≳ 1), and less luminous than its radiatively efficient counterpart.The electron and ion plasma in a RIAF can form a twotemperature flow; most of the energy generated by accretion is stored within the ion plasma as thermal energy, while the electron plasma cools more efficiently and is responsible for the features of the emitted radiation.
MBBH mergers in RIAFs have been explored in simulations by different groups.General relativistic hydrodynamic (GRHD) simulations of RIAF-like accretion flows were carried out independently by Bode et al. (2010) and Farris et al. (2010).These works presented a suite of full GR simulations of equal-mass, nonspinning binary black hole mergers immersed in a gaseous environment.They found that because of high thermal velocities, the radial inflow speeds of gas are comparable to the orbital speed at a given radius.Thus, in a hot accretion flow, binary torques are not capable of clearing a central cavity because the gas ejected by the binary is refilled on a dynamical timescale.
Aiming at identifying possible characteristic variability of the emitted EM signals, these works track the evolution of the rest-mass accretion rate, as well as the EM luminosity due to bremsstrahlung emission.Farris et al. (2010) also considers synchrotron emission assuming the presence of a small-scale, turbulent magnetic field with a magnetic-to-gaspressure ratio β −1 = 0.1.Both works established that the merger of equal-mass, nonspinning binaries in hot accretion flows features an EM flare followed by a sudden drop-off in luminosity.Additionally, Bode et al. (2010) argued that, depending on the relative orientation between the observer and the binary, the variation in the beaming pattern of the binary could potentially give rise to modulations in the detected luminosity resulting in an "EM chirp" analogous to the gravitational one.The subsequent follow-up by Bode et al. (2012) extended those early results investigating the effects of varying the magnitude and orientation of the individual BH spins and the binary mass-ratio.Binaries with q = 0.5 were found to exhibit lower and narrower luminosity peaks relative to the equal-mass case; also, q = 0.5 binary black holes with spins misaligned with respect to the orbital angular momentum were found to produce peaks at even lower values.These results are interpreted as the consequence of the lower effi-ciency by unequal-mass/precessing binary torques at driving shocks in the gas.
The early GRHD studies of MBBHs in gas clouds neglected the important role that magnetic fields are likely to play in driving the plasma dynamics, in the radiation emission mechanisms, and in forming relativistic jets.EM fields were first included in general relativistic electrovacuum simulations by Palenzuela et al. (2009); Mösta et al. (2010) and in general relativistic force-free electrodynamics (GRFFE) simulations by Palenzuela et al. (2010b,a); Mösta et al. (2012); Alic et al. (2012).These ground-breaking studies explored binary black hole mergers in a magnetically dominated plasma.Such a scenario is illustrative of conditions that could arise if a massive binary decouples from the CBD, but continues to interact with the magnetic field anchored to it.Noticeably, these simulations present evidence for the formation of separate, strongly collimated jets around each BH during the inspiral; at merger, these two jets would coalesce into a single jet around the remnant Kerr BH.The post-merger behavior of this EM outflow is well represented by the Blandford-Znajek process (Blandford & Znajek 1977).
More recent studies explored the behavior of moderately magnetized accretion flows (MMAF) around binary black holes in an ideal-GRMHD context.Giacomazzo et al. (2012) presented ideal-GRMHD simulations of merging equal-mass, nonspinning binary black holes immersed in a MMAF using the WhiskyMHD code (Giacomazzo & Rezzolla 2007).Giacomazzo et al. (2012) considered two models for the gas cloud surrounding the binary, both with an initially uniform rest-mass density ρ 0 : an unmagnetized plasma, and a plasma threaded by an initially uniform magnetic field with an initial ratio of magnetic-to-fluid pressure β −1 equal to 0.025.Though this study was limited to evolve a limited number of orbits (∼3) before merger, it showed a rapid amplification of the magnetic field of approximately two orders of magnitude.This results showed that MMAFs exhibit different dynamics compared to unmagnetized accretion flows.The accretion flow in magnetized environments yields turbulent motion in the gas near the inspiralling BHs, eventually leading to the formation of a thin, disk-like structure rotating on the equatorial plane of the remnant Kerr BH.Also, MMAFs can lead to strong, collimated EM emission: in fact, the Poynting luminosity resulting from their magnetized binary configuration was estimated to be approximately four orders of magnitude larger than the luminosity obtained in the GRFFE simulations by Palenzuela et al. (2010b); Mösta et al. (2012).This indicates that the dynamics of the lateinspiral and merger of BBHs has a strong role in driving the magnetic fields in its surroundings.
The results of Giacomazzo et al. (2012) were extended by Kelly et al. (2017), who covered a broader collection of physical scenarios adopting the IllinoisGRMHD code (Noble et al. 2006;Etienne et al. 2015) to solve the GRMHD equations.The simulations of Kelly et al. (2017) consider equal-mass binaries of nonspinning black holes with initial separations larger than those considered in Giacomazzo et al. (2012).Evolving higher-separation binaries allowed the plasma in the strong-field regions to establish a quasiequilibrium flow with the binary motion; in addition, it allows for better resolution of the timing features of the EM emission.Several configurations differing only in the initial magnetic field value were evolved, showing that the level of the Poynting luminosity reached during the inspiral and merger is little sensitive to the initial magnetic field strength.Also, Kelly et al. (2017) carried out a simplified calculation of the direct EM emission generated from the plasma by post-processing the GRMHD data with the Monte Carlo radiation transport code PANDURATA (Schnittman & Krolik 2013).Kelly et al. (2017) assumed a simplified emission model of thermal synchrotron and found that the locally generated EM luminosity remain steady throughout the inspiral, exhibits a small burst before merger, and decreases by ∼50% afterwards.Based on this simple emission model, the direct emission luminosity is orders of magnitude lower than the Poynting luminosity.
The astrophysical scenario examined by Giacomazzo et al. (2012); Kelly et al. (2017) was further investigated in Cattorini et al. (2021Cattorini et al. ( , 2022)), which have produced three suites of GRMHD simulations of MBBHs immersed in a gas cloud and expanded the parameter space by considering binaries of spinning BHs with general spin orientation and unequal-mass ratios.Cattorini et al. (2021) evolve a suite of nine simulations covering a range of initially uniform, moderately magnetized fluids with different initial magnetic-to-gas-pressure ratio β −1 0 .For each value of β −1 0 , distinct spin configurations defined by adimensional spin parameters a = (0, 0.3, 0.6) are considered.For a given initial magnetization, it is shown that (aligned) spins of the individual BHs have a suppressing effect on the accretion rate as large as ∼48%.In addition, spin is found to affect the peak luminosity reached shortly after merger, which is enhanced by up to a factor of ∼2.5 for binaries of spinning BHs compared to the nonspinning models.This intensification does not depend on the initial level of magnetization β −1 0 of the gas.The follow-up study by Cattorini et al. (2022) evolved five equal-mass binary configurations characterized by differing spin-alignment, and investigate how the spin inclination modifies the magnetohydrodynamical behavior of the plasma during the binary late-inspiral and merger.Cattorini et al. (2022) identify quasiperiodic modulations in the mass accretion rate, which appear to be related to the spin-alignment.To  explore the variability of these timing features, frequency analysis is carried out calculating the power spectral density of the accretion rate for the five configurations.Furthermore, they perform wavelet power spectral density (WPSD) analysis on one misaligned-spinning simulation and identify a "chirp" in the accretion rate that strongly resembles the gravitational one (Fig. 7).
The link between spin-alignment and quasiperiodicity in the accretion rate has been further investigated in the recent work by Cattorini et al. (2023), which produced several unequal-mass models and considered a broader family of spin misalignment and found a relation between the orientation of the individual spins relative to the orbital angular momentum and the amplitude of the quasiperiodic oscillations of the accretion rate.Also, they calculated the spinprecession parameter χ p (Schmidt et al. 2015;Gerosa et al. 2021) for their misaligned-spinning models and found linear trend connecting binary systems with larger χ p (and hence, higher precession) to larger amplitude of the modulations in the accretion rate.This finding could provide a step toward posing an additional means to characterize spin precession and determine it observationally.

DISCUSSION
As one of the loudest sources of GWs, merging MBBH systems provide a unique dynamical probe of strong-field general relativity and a fertile ground for the observation of fundamental physics.In the next decade, these systems will be a major target of the Laser Interferometer Space Antenna (LISA, Amaro-Seoane et al. 2012, 2017;Amaro-Seoane et al. 2022), which will detect MBBHs in the (restframe) mass range ∼ [10 4 − 10 7 ] M ⊙ at rates between a few and about a hundred per year (Klein et al. 2016).Their gravitational signal is expected to be accompanied by an EM counterpart, which can be triggered by the accretion of gas onto the black holes across the inspiral, merger, and ringdown.To maximize the scientific return of LISA, several groundand space-based missions will complement GW detection in the EM domain, looking at various frequencies for spectral and timing features characteristic of a MBBH system.The detection of an EM counterpart would allow for an accurate identification of the GW source sky position (i.e., of the host galaxy), enabling the estimation of the source redshift either spectroscopically or photometrically.However, the multimessenger detection of the gravitational and EM radiation emitted by a LISA source will be no easy task, since LISA sources will be localized within a sky-area ranging from hundreds of deg 2 to fractions of a deg 2 (Mangiagli et al. 2020), depending on intrinsic features of the binary and its luminosity distance.To be univocally identified as the counterpart to a GW signal, the EM counterpart must (i) be bright enough to emit radiation above the detection limit of the EM observatory and (ii) have distinctive characteristics that would allow to discriminate among the other sources present in the observatory's field of view (e.g., the periodicities that may arise due to the dynamical coupling of the binary with its surroundings).
The features of such an emission largely depend on the amount of gas in the vicinity of the binary and on the geometry of the accretion flow.To date, the specification of the counterpart signal to a LISA MBBH source is still unclear, due both to the scarcity of observational evidence (De Rosa et al. 2019) and the complexity of producing accurate numerical simulations that take into account all the intricate physical mechanism characterizing these systems.To this aim, over the last years several numerical investigations of MBBHs in astrophysical environment have contributed to the theoretical characterization of these sources and advanced our understanding of the properties of the material accreting onto massive binaries, providing a fundamental step towards the accurate modeling of EM counterparts to strong LISA sources.As a massive binary separation decreases to scales comparable to the gravitational radius, gravity becomes the dominant force in driving the fate of the system.To produce accurate numerical evolution of the late-stage of MBBH inspiral and merger it is essential to solving the equations of general relativity and magnetohydrodynamics.
As the evolution of a MBBH system is driven by emission of gravitational radiation, the morphology of the accretion flows will likely depend on the balance of heating and cooling mechanisms in the gas.Two scenarios have efficiently been investigated by GRMHD simulations, i.e., the circumbinary disk, and the gas cloud models.Exploring different possibilities related to the binary environment allows for a deeper understanding of the physical mechanisms that can be associated to detectable EM signals.At present, theoretical efforts involving MBBHs in circumbinary disks have incorporated relativistic magnetohydrodynamics in the numerical evolution of the accreting gas and evolved the binary spacetime with a hybrid (PN-evolved) metric or by evolving the full set of general relativistic field equations.These studies adopted the ideal GRMHD limit to probe EM signals generated from magnetized, circumbinary accretion flows onto binaries of equal/unequal-mass, nonspinning/spinning BBHs.
The timing features of the accretion rate are found to display diverse quasiperiodic behavior depending on the binary environment.Full general relativistic explorations of circumbinary accretion observed that the accretion rate displays modulations with f ∼ 0.75 f orb .Conversely, circumbinary GRMHD simulations with PN gravity have shown that the accretion of mass onto the individual MBHs strongly depends on the beat frequency set by the orbital motion of the lump.
In contrast with numerical simulations of circumbinary accretion, GRMHD studies of MBBHs in magnetized gas clouds have focused on structureless uniform distributions of plasma threaded by an asymptotically uniform magnetic field, which are distinctive of hot and radiation dominated accretion flows.These studied have displayed that magnetic fields can be distorted and significantly increase their strength, developing EM jets in the polar regions above inspiraling MBHs and ultimately producing a jet around the spin axis of the remnant BH.The geometry of these structures in the vicinity of the black hole horizons depends on the system's mass ratio and the individual spins of the BHs.Quasiperiodic oscillations of the mass accretion rate with a dominant frequency f ∼ 2 f orb are observed, suggesting that such modulations (which may be associated with observable EM signals) are not exclusive of environments in which a binary is embedded in a circumbinary accretion disk.Noticeably, these quasiperiodicities are linked to the orientation of the individual spins; a higher spin-misalignment corresponds to sharper quasiperiodicities in the accretion rate and to larger modulation amplitudes.
Exploring general relativistic features of diverse geometries is key to develop a coherent picture of astrophysical accretion flows onto MBBHs.The ultimate goal of their numerical exploration is to provide accurate predictions of their electromagnetic signatures.The strongest and most valuable predictions would be a little sensitive to different details that may characterize a binary environment.In this regard, a robust validation of our results would arise from observing similar features emerging when considering different gaseous configurations (see, e.g., Fedrigo et al. 2023).A greater understanding of the fueling rate, geometry, and magnetic-field properties of the accreting gas will be essential in advancing the future of numerical relativity simulations and enhancing the accuracy of EM predictions.Additionally, it will be crucial to effectively incorporate various spatial scales to adequately study the evolution of the accreting gas from its early inspiral phase to the stage of merger.GRMHD simulations should also consider radiation processes in order to accurately predict EM light curves and spectra.Furthermore, they must account for different modes of accretion and incorporate radiation feedback (see, e.g., Sadowski & Gaspari 2017).Thus, the development of reliable radiation-transport methods in dynamical spacetime (general-relativistic radiative magnetohydrodynamics, GR-RMHD) is of utmost importance.

Figure 2 .
Figure 2. Artistic illustration of the two accretion scenarios discussed in this section.Left: circumbinary disk model in which binary torques create a central low-density cavity.The accretion onto the individual black holes is mediated by individual mini-disks around each BH.Right: gas cloud scenario in which MBBHs remain engulfed in the hot thin gas until coalescence.At larger radii, the geometrically thick accretion flow transitions into a cold, geometrically thin disk.The illustrations are not to scale.Figure adapted from Bogdanović et al. (2011).

Figure 3 .
Figure 3. Color contours of surface density in units of Σ 0 (i.e., the maximum surface density in the initial condition) as a function of radius and azimuthal angle in RunSE at four different times.The linear color scale emphasizes the growth of the "lump" in the inner disk.The times shown are t = 40000M (upper-left), t = 51963M (upper-right), t = 63926M (lower-left), and t = 75890M (lowerright), where M is the total gravitational mass of the system.Figure taken from Noble et al. (2012).

Figure 4 .
Figure 4. Time-averaged luminosity (νL ν ) spectrum obtained at d = 1000M and θ = 0 using data from the second orbit of the binary simulation produced by Bowen et al. (2018).Separate contributions from the mini-disk regions, the accretion streams, and the circumbinary region are displayed.Figure taken from D'Ascoli et al. (2018).

Figure 5 .
Figure 5. Rest mass density profile on the equatorial (xy) plane of the aligned-spin model χ ++ of Bright & Paschalidis (2023) illustrating the mini-disk structures around each black hole after ∼13 binary orbits.Figure taken from Bright & Paschalidis (2023).

Figure 7 .
Figure 7. (Top row, left): GW strain of the misaligned-spinning model UUMIS by Cattorini et al. (2022) normalized to its maximum value at merger.(Top row, right): time-dependent premerger mass accretion rate over both BH horizons for UUMIS run, also normalized to its maximum value.(Bottom row, left): time-frequency representation of the GW strain via wavelet PSD, showing the frequency increase of the signal over time (the so-called "chirp").(Bottom row, right): time-frequency representation of the accretion rate via wavelet PSD, showing similarity with the GW strain in the frequency increase over time.Time and frequency units are normalized to a binary system with total mass M = 10 6 M ⊙ , and M 6 ≡ M/10 6 M ⊙ .Figure taken from Cattorini et al. (2022).