Black Hole–Disk Interactions in Magnetically Arrested Active Galactic Nuclei: General Relativistic Magnetohydrodynamic Simulations Using a Time-dependent, Binary Metric

Perturber objects interacting with supermassive black hole accretion disks are often invoked to explain observed quasiperiodic behavior in active galactic nuclei (AGN). We present global, 3D general relativistic magnetohydrodynamic (GRMHD) simulations of black holes on inclined orbits colliding with magnetically arrested thick AGN disks using a binary black hole spacetime with mass ratio 0.1. We do this by implementing an approximate time-dependent binary black hole metric into the GRMHD Athena++ code. The secondary enhances the unbound mass outflow rate 2–4 times above that provided by the disk in quasiperiodic outbursts, eventually merging into a more continuous outflow at larger distances. We present a simple analytic model that qualitatively agrees well with this result and can be used to extrapolate to unexplored regions of parameter space. We show self-consistently for the first time that spin–orbit coupling between the primary black hole spin and the binary orbital angular momentum causes the accretion disk and jet directions to precess significantly (by 60°–80°) on long timescales (e.g., ∼20 times the binary orbital period). Because this effect may be the only way for thick AGN disks to consistently precess, it could provide strong evidence of a secondary black hole companion if observed in such a system. Besides this new phenomenology, the time-average properties of the disk and accretion rates onto the primary are only marginally altered by the presence of the secondary, consistent with our estimate for a perturbed thick disk. This situation might drastically change in cooled thin disks.


INTRODUCTION
Active galactic nuclei (AGN) produce copious amounts of electromagnetic radiation with strong variabilities occurring on timescales of minutes to years at all observed frequencies (Ulrich, Maraschi & Urry 1997;Czerny 2004;Ricci & Trakhtenbrot 2023;Sartori et al. 2018).Although this variation is mainly stochastic (dominated by red noise) some of these sources seem to exhibit (quasi-)periodic emission (Graham et al. 2015;Charisi et al. 2016;D'Orazio & Charisi 2023).A natural mechanism to explain possible periodic behavior in AGN is the presence of a binary system in the central engine, e.g., a supermassive binary black hole (Begelman, Blandford & Rees 1980).In this scenario, (magneto-)hydrodynamical accretion processes (Noble et al. 2021) and kinematic effects such as Doppler boosting (D'Orazio, Haiman & Schiminovich 2015) might link the AGN periodicity with the binary properties.In particular, repeated collisions between a smaller binary companion and the accretion disk surrounding the central supermassive black hole has long been discussed as a possible mechanism for quasiperiodic behavior (Zentsova 1983;Lehto & Valtonen 1996;Ivanov, Igumenshchev & Novikov 1998;Komossa 2006;Dai, Fuerst & Blandford 2010;Franchini et al. 2023;Linial & Metzger 2023a, 2024;Pasham et al. 2024).Typically it is proposed that the secondary object is either a star or a smaller mass black hole, which impacts the disk, ejects matter, and/or heats up the disk material via bow shocks.For instance, the source OJ287 shows consistent quasi-periodic outbursts at optical wavelengths and is supposed to host a supermassive binary black hole system (though the masses of the black holes are uncertain, Sillanpaa et al. 1988;Dey et al. 2018;Laine et al. 2020;Komossa et al. 2023).Moreover, one possible explanation for recent observations of quasi-periodic Xray eruptions in AGN (Miniutti et al. 2019) is an association with flaring emission due to disk collisions (Linial & Metzger 2023a;Franchini et al. 2023;Tagawa & Haiman 2023;Xian et al. 2021;Zhou et al. 2024).
The case where the secondary object is a black hole is particularly interesting.For a given effective size, black holes are much more massive than stars and so they can have a stronger influence on the primary black hole at close separations through effects like spin-orbit coupling, where the orbital angular momentum and black hole spin directions will oscillate as a function of time.In fact, this has been proposed as a possible mechanism for periodic radio variability in AGNs (von Fellenberg et al. 2023) and optical variability in blazars (Abraham & Romero 1999;Romero et al. 2000;Britzen et al. 2023) If this model is correct, then it would mean that these quasi-periodic observational features are also the electromagnetic counterparts to low frequency gravitational wave sources expected to be detected by the Laser Interferometer Space Antenna (LISA, Flanagan & Hughes 1998;Amaro-Seoane et al. 2007;Berry & Gair 2013) at merger, and pulsar timing arrays during the inspiral regimes.Modelling these systems in detail could then have important implications for multi-messenger astronomy.
Because of the strong gravitational field and nonlinear magnetohydrodynamics involved, self-consistent models of black hole-disk interaction can only be accurately built using 3D General-Relativistic Magnetohydrodynamical (GRMHD) simulations.There are essentially two different approaches to simulate binary black hole accretion systems in General Relativity (for detailed reviews, see Gold 2019;Cattorini & Giacomazzo 2024).The most accurate approach is to couple MHD to numerical relativity (e.g., Farris et al. 2012;Giacomazzo et al. 2012;Gold et al. 2014;Paschalidis et al. 2021); this method, however, is computationally expensive and numerically challenging.In practice these simulations can only be evolved for a short amount of time (typically ∼ tens of orbits) close to merger and are usually far from steady-state.A computationally cheaper and numerically simpler approach is to approximate the dynamical spacetime by some semi-analytic expression for use in existing GRMHD codes (e.g., Noble et al. 2012;Bowen et al. 2018;Gold et al. 2014;Mundim et al. 2014;Lopez Armengol et al. 2021;Avara et al. 2023).The approximate metric must be constructed with care to be globally applicable and an accurate representation of the evolving spacetime.The accuracy can also become poor at times close to merger.Lopez Armengol et al. (2021) and Combi et al. (2021) recently proposed an approximate metric constructed by superimposing two boosted Kerr black holes, demonstrating it to be accurate during the inspiral regime leading up to merger and welldefined for the entire domain.The robustness of the metric in GRMHD was exhibited in a follow-up work, where it was applied to an equal mass binary system with spinning black holes (Combi et al. 2022).
Using semi-analytical approximations for the binary black hole spacetime as in Combi et al. (2021) allows for a much larger exploration of the vast parameter space in binary black hole accretion.This is not only because avoiding numerical relativity makes the computation faster, but it will allow the community to take full advantage of the numerous advances made in single black hole GRMHD codes if these are adapted for time-dependent metrics.For instance, in the past couple of decades, the GRMHD accretion community has been able to explore effects such as varying black hole spin (Event Horizon Telescope Collaboration et al. 2019;Akiyama et al. 2022), varying magnetic flux supply (Tchekhovskoy, Narayan & McKinney 2011;Narayan et al. 2012), varying disk tilt (McKinney, Tchekhovskoy & Blandford 2013;Liska et al. 2018;White, Quataert & Blaes 2019;Chatterjee et al. 2020), including radiative effects (Ryan, Dolence & Gammie 2015;White et al. 2023), including nonideal physics (Ressler et al. 2015;Foucart et al. 2017;Chandra, Foucart & Gammie 2017;Ripperda et al. 2019), and studying a variety of initial conditions informed by larger scales (Ressler et al. 2020(Ressler et al. , 2021;;Kaaz et al. 2023;Cho et al. 2023;Lalakos et al. 2023).Furthermore, since the user base for GRMHD is currently much larger than that for numerical relativity (e.g., Porth et al. 2019), it could encourage more researchers to study binary black hole systems.
Recently, Suková et al. (2021) investigated the collisions of spherical objects with AGN disks using GRMHD simulations.The simulations were essentially agnostic to the type of secondary object (e.g., a star or black hole), and simply enforced that a region with a particular "influence radius" be comoving with the object.The authors found that the presence of the secondary could significantly modify the structure of the accretion flow and produce strong outbursts of relativistic outflow 1-2 orders of magnitude larger than the "background" outflow rate.Such strong features in the accretion and outflow properties would seem more than enough to explain some of the observed quasi-periodic behavior in AGN.However, the simulations of Suková et al. (2021) were primarily performed in 2D axisymmetry, with only one 3D perturber simulation with a limited runtime (5000 primary black hole light crossing times) and initialized with 2D data.Since accretion flows are known to behave much differently in 3D than in 2D and the geometry of the perturber scenario is fundamentally three dimensional, it is therefore important to revisit this problem using full 3D simulations (and a fully general relativistic treatment of secondary black holes).This is especially true for magnetically arrested accretion disks (MADs), where the axisymmetry causes the flow to be fully halted in 2D, while in 3D the gas can penetrate the magnetic barrier (McKinney, Tchekhovskoy & Blandford 2012;White, Stone & Quataert 2019;Ripperda, Bacchini & Philippov 2020;Ripperda et al. 2022).
Here we present our implementation of a time-dependent binary black hole metric into Athena++ (White, Stone & Gammie 2016;Stone et al. 2020) and present a series of appropriate test problems for time-dependent metrics.We use a superposed metric as in Combi et al. (2021), generalized for arbitrary spins and eccentricities that result from solving the Post-Newtonian equations for the evolution of the black holes.Athena++ is particularly suited for the binary black hole accretion problem with small mass-ratios due to its ability to use adaptive mesh refinement and excellent scalability.As a conservative code using constrained transport for magnetic fields it conserves mass, energy, momentum, and magnetic flux to machine precision.It also now has full support for radiation, being the first GRMHD code to directly solve the GR Boltzmann transport equation (White et al. 2023).Finally, it is widely used, public, and a ported version for use on GPUs will soon be publicly available (AthenaK).
Although our implementation is general to any binary black hole configuration, in this work we focus on the study of small mass ratio (0.1) inspirals in galactic centers with thick disks.In doing so we seek to improve on past work studying black hole-disk interactions with simulations.We do this by evolving a 3D torus around a single supermassive black hole long enough to reach equilibrium out to ∼ 100 gravitational radii and then introducing a secondary black hole at some initial distance and inclination to the disk.The secondary then moves through the disk based on the solution to the post-Newtonian orbital equations and alters the accretion and outflow of the disk due to gravitational, electromagnetic, and plasma interactions.We also study how the effects of spin-orbit coupling alter the accretion flow as the primary black hole spin axis changes direction as a function of time and study the accretion properties of the secondary.We defer the study of radiatively cooled thin disks (more common for AGN) to later work.
This paper is organized as follows.§2 details the analytic and numerical methods we use to simulate binary black hole systems, §3 presents some useful analytic estimates, §4 describes our simulation results for small mass ratio inspirals, §5 compares our simulations to past work, §6 discusses the limitations of our work, and §7 concludes.

Time-Dependent Metrics in Athena++
We start with the publicly available multi-purpose fluid dynamics code Athena++ (Stone et al. 2020), particularly the extension for GRMHD (White, Stone & Gammie 2016), which solves the conservative equations in an arbitrary spacetime set by the choice of coordinates and metric through the "GRUser" module.In principle, this module can accom-modate any choice of time-independent metric as long as the spatial derivatives (used to compute the connection coefficients) are supplied.We extend this framework to support time-dependent metrics (in a similar way to Noble et al. 2012).In the 3 + 1 conservative formulation of GRMHD (e.g., Gammie, McKinney & Tóth 2003), the mass, energy, momentum, and magnetic field equations in a coordinate basis (and Lorentz-Heaviside units) are where is the stress-energy tensor, t is the coordinate time, x i are the spatial coordinates (x, y, and z), g is the determinant of the metric, ρ is the rest-frame mass density, u µ is the gas fourvelocity, Γ α βκ is the connection coefficient, γ is the adiabatic index, P is the rest-frame gas pressure, B i is the magnetic field, b µ is the four-magnetic field defined via and g µν is the metric.Here the indexes i (not to be confused with the orbital inclination used later in this work) are limited to the spatial directions (i ∈ {1, 2, 3}), while Greek indexes include the time direction (µ ∈ {0, 1, 2, 3}).The connection coefficients can be written in lowered form as As formulated, equations (1-4) are valid for both timeindependent and time-dependent metrics.Athena++, however, like all GRMHD codes designed for applications to stationary metrics, assumes that the terms proportional to ∂g µν /∂ t in Equation (4) are 0. We are therefore required to generalize the calculation of Γ αβκ , or, more precisely, T α β Γ β να .This is simplified by using the symmetry of the stress-energy tensor, T µν = T νµ , the symmetry of the metric, g µν = g νµ , and the fact that the last two terms on the right hand side of Equation ( 4) are antisymmetric in the first and last index of Γ αβκ , so that This means that the nonzero ∂g βα /∂t only explicitly appears in the GRMHD equations as a source term in the conserved energy, T t t .The other equations remain unchanged, although √ −g can now depend on time.
The set of equations (1) are solved by Athena++ using a finite volume method, that is, integrated over space and time for each three-dimensional cell and time step.For instance, integrated over volume, and assuming uniform spacing in discretization, the left-hand side of the hydrodynamic conservation equations in (1) become (in one spatial dimension without loss of generality): where U is a conserved variable (not including B i ), F x1 is the associated flux in the x 1 direction, S is the source term, n, n + 1/2, and n + 1 denote the time at the initial, half, and final stages of a given time step, i and i ± 1/2 denote the cell center and cell faces of the ith grid cell, ∆t is the time step, denotes an average over the area of the face with normal in the x 1 direction at a fixed time, and denotes an average over the cell volume at fixed time.Similarly, the magnetic field equation in the x 1 direction becomes (in 3D): where are the electric fields (i.e., the covariant version of the standard E = −v × B), and denotes an average over the length of an x 2 edge at a fixed time.
Note that for static metrics, the factors of √ −g on the lefthand side of Equations ( 6) and ( 9) cancel, but for general time-dependent metrics they need to be accounted for.
In Appendix A.1, we present a 1D test of a particular time-dependent metric and demonstrate the code's secondorder convergence when the metric is updated every subtimestep.When updated only every timestep or any multiple of a timestep the code converges at first order.

Approximate Binary Metric
To establish notation and for explicitness, in this subsection, we review the approximate binary black hole metric presented in Combi et al. (2021) in which the authors superimpose two boosted individual Kerr metrics 1 .
As in Lopez Armengol et al. (2021), we will use Kerr-Schild gauge of the superposition but following the construction in Combi et al. (2021).In the rest frame of an isolated Kerr black hole, the metric in Cartesian Kerr-Schild coordinates (CKS, Kerr 1963) is where η µν is the Minkowski metric, M BH is the black hole mass, a i = (a X , a Y , a Z ) is the dimensionless black hole spin in the X, Y, and Z directions, and a = a 2 X + a 2 Y + a 2 Z .This useful form of the KS metric for arbitrary spin direction was presented in the Appendix of Ma et al. (2021) for the first time (as far as we know).Here, X i = (X, Y, Z) are the black hole rest frame spatial coordinates.Now consider a Kerr black hole moving on a trajectory given by s µ (τ) = [t τ (τ), s x (t τ ), s y (t τ ), s z (t τ )], where τ is the proper time of the black hole.This is related to the time in the inertial lab frame, where , and β i (t) = ds i (t)/dt.Following Combi et al. (2021; see also Mashhoon & Muench 2002), for every space time event x µ = (t, x, y, z) in the lab frame we construct a coordinate system centered on the black hole with a proper time τ such that the given event and the black hole trajectory are simultaneous.The lab frame time corresponding to this τ, t τ (τ), is not the same as the lab frame event time t unless the black hole is not moving.Mathematically, this corresponds to the the relation: where Λ µ ν is the instantaneous Lorentz transformation: where n i = n i = β i /β and the time dependence of β and Γ is assumed implicitly.Equation ( 16) defines the relationship between X µ = (τ, X, Y, Z) and x µ = (t, x, y, z) at every point in time and space given the black hole trajectory.Note that we have chosen coordinates where the lab frame x, y, and z axes are aligned with the black hole frame X, Y, and Z axes (since the spin axis of the black hole can be in an arbitrary direction this choice comes without loss of generality).
Equation ( 16) results in a nonlinear set of equations (see Equations B6 and B7 in Appendix B), which, given a black hole trajectory can be solved numerically at each lab frame location and time to obtain t τ and thus determine X µ in terms of x µ .Doing so, however, would introduce a significant extra computational cost for computing the metric and potentially reducing the overall speed of the simulations (though it is not obvious by how much).Not only that, but the coordinates become ill-defined for larger distances from the black hole if the trajectory is accelerating.Primarily motivated by this coordinate issue, instead of solving the nonlinear coordinate transformation we make the approximation that s i (t τ ) ≈ s i (t) + β i (t)(t τ − t) and β i (t τ ) ≈ β i (t).In general, this is a good approximation close to the black hole where |t τ − t| is small but it has the potential to be highly inaccurate at large distances from the black hole when |t τ −t| is large.For the specific case of an orbiting black hole, however, there is an upper limit to the error incurred due to this approximation on the rest frame distance from the black hole.This error is roughly equal to the typical β 2 of the orbit (e.g., ≲ 5% for a black hole on a circular orbit around a much larger black hole at a separation of ≳ 20 light crossing times of the larger black hole).This is because the relative difference in black hole position for any two given points along the orbit become negligible at larger distances from the system.The potential ≲ 5% error in distance could cause small differences in the gas distribution at large distances from the black hole, however, these differences would be smaller than the uncertainties in, e.g., the initial conditions of black hole accretion flows.
Using this approximation, Equations (B6-B7) become: and where β i and Γ are now evaluated at t instead of t τ .Taking the derivatives of , one can show that dx µ = Λ µ α (t)dX α , where Λ µ ν is the standard Lorentz transformation given by Equation ( 17).
Using this boost and the coordinate transformation we obtain our approximate binary metric: where the subscripts 1 and 2 indicate the primary and secondary black holes.We obtain the inverse metric, g µν , as well as the spatial and temporal derivatives numerically.
In Appendices A.2 and A.3 we present two different tests involving moving black holes using this metric, demonstrating our code's convergence properties and the accuracy of our implementation.

Post-Newtonian Orbits of the Black Holes
The above metric has as an input the black hole trajectories, which have to be solved for independently.To do this, we use the public code CBwaves (Csizmadia et al. 2012) to evolve the trajectories of the two black holes starting from a given eccentricity, separation distance, and initial spin directions and magnitude with respect to the orbital plane.CBwaves is a fast C code that solves the post-Newtonian equations of motion (Blanchet 2014) up to the 4th PN order and includes all the relevant acceleration terms, radiation reaction, spin-spin coupling, and spin-orbit coupling.It has been tested in Csizmadia et al. (2012) and is one of the few public codes to solve the PN equations (a coupled set of ordinary differential equations with lengthy terms) in full generality.The outputs of the code are the two black hole 3D positions, 3D velocities, and 3D spins as a function of time up to merger.For simplicity, in this work, we always use zero initial eccentricity but allow for the black hole spin and angular momentum vectors to be misaligned.This can lead to significant precession due to spin-orbit coupling, as we shall see in some of our simulations.We define the initial inclination angle, i 0 , such that for i 0 = 0, the orbit is initially prograde with the accretion disk in the midplane, and for i 0 = 90 • the orbit is initially perpendicular to the midplane and moves clockwise in the y-z plane.We distinguish between i as the current inclination of the orbit (always defined with respect to the fixed x-y plane) and i 0 , the initial inclination of the orbit.

Further Approximations Specific to This Work
In this work, we apply the boosted binary metric to a system in which the mass ratio, M 2 /M 1 ≡ q ≪ 1, where M 1,2 are the masses of the primary and secondary black hole, respectively.Therefore, we approximate the primary black hole as stationary and located at the origin of the lab frame.This approximation would break down for higher mass ratios, q ≲ 1, at which point the qualitative picture of a perturber black hole impacting an established AGN disk becomes inaccurate.We also neglect self-gravity of the accretion flow and radiative effects.This formally limits the applicability of our study to lower mass accretion rates where these effects are negligible.At higher accretion rates self-gravity could introduce enough dynamical friction and drag to alter the black hole orbits, while radiation would significantly cool the disk and reduce its scale height, altering the flow dynamics substantially.The latter regime we plan on studying in future work.

Algorithmic Details
For stability purposes, the time step used by Athena++ is set by the light crossing time across the shortest x, y, or z length of a cell.In particular, we choose a Courant-Friedrichs-Lewy (CFL) number of 0.3.Because of this, the timestep can be significantly shorter than the characteristic time for the metric to change (e.g., an orbital time for a black hole binary metric, which for the orbits we consider in this work is proportional to v ≲ 0.1c).Taking advantage of this fact, we only update the metric every 10 timesteps.This means that the MHD equations are first solved including the time-dependent source terms described at the beginning of §2 but otherwise as if the metric were time-independent.That is, the conversion between conserved and primitive variables (and vice-versa) is done using the metric at the time of the last update, not the metric at the current time of a given step or substep of the algorithm.
As a result, the code is formally only first-order accurate in time.However, when the metric changes at a rate much slower than that of the fluid (i.e., when the black holes move at velocities much slower than the characteristic GRMHD velocities), the errors incurred by the first-order metric update can still be much less than the errors incurred by the secondorder GRMHD evolution.Quantifying this more precisely is difficult and is likely different for every simulation and choice of parameters.That said, for two 3D GRMHD accretion problems with moving black holes that are similar in several ways to our target problem (see Appendix A.2,Appendix A.3, we have found that updating the metric once every 10 timesteps still results in satisfactory agreement with the expected solutions.This is true even though the black hole velocities in those test problems are much higher (0.9c) than those we expect in our simulations (≲ 0.25c for orbits with separations ≳ 20r g ).Moreover, when discontinuities or shocks are present in the flow (as they are in turbulent accretion simulations in general but especially for the Bondi-Hoyle-Lyttleton-type flows we expect near each black hole) all methods reduce to first-order anyway.

Floors and Patches to the Metric
To avoid coordinate singularities within the event horizon of the black holes we modify each black hole's coordinates when r ≤ 0.8r H (defined in each black hole's rest frame), where r H = M BH + M 2 BH − a 2 is the event horizon radius for an isolated black hole.For r ≤ 0.8r H we first calculate where l µ is defined in Equation ( 13).Then we set r = 0.8r H and recalculate from the new r and old θ and φ.Since this coordinate modification is only applied well within the event horizon it should have no effect on the simulation outside the horizon and it prevents occasional NaNs from crashing the simulation.Within the horizon of each black hole, we also set the gas to be moving along with the black hole by setting u µ = Λ µ ν (u ′ rest ) ν , where (u ′ rest ) µ is the four velocity of a stationary observer in the instantaneous black hole rest frame: where α ≡ 1/ −g tt and Λ µ ν is the Lorentz transformation defined in Equation ( 17).This helps prevent gas and magnetic fields from within the event horizon from 'leaking' out into the rest of the computational domain as the black hole moves across the grid.In particular, without enforcing this velocity condition we have found that 'magnetic explosions' caused by unphysically large magnetic fields leaking out of the horizon can ruin the simulation.
For the MHD quantities, the density floor is 10 −6 (r/r g ) −3/2 and the pressure floor is 3.33 × 10 −9 (r/r g ) −5/2 , with σ ≡ b 2 /ρ ≤ 100 and β ≥ 0.001 enforced via additional density and pressure floors, respectively.Here β is the ratio between the thermal and magnetic pressures while b 2 is twice the magnetic pressure in Lorentz-Heaviside units.Additionally, the velocity of the gas is limited such that the maximum Lorentz factor is 50.The radial power law indices of the pressure and density floors are chosen to be consistent with spherical Bondi-type accretion flows appropriate for a nonrotating, low density atmosphere surrounding the accretion flow.The precise magnitudes of these floors have negligible effects on the accretion flow (Porth et al. 2019) and are chosen to be several orders of magnitude less than the initial density maximum (=1 in code units).The σ and plasma β limits help prevent primitive inversion failures in strongly magnetized regions while the limit on Lorentz factor helps Figure 1.2D slice of density overplotted with magnetic field lines at t = 100,000M for the single a = 0.9375 black hole torus simulation used as the 'initial' state of our small mass ratio binary simulations.The black circle at the center represents the event horizon.The flow is thick, turbulent, magnetically arrested, and is associated with a powerful jet in the direction of the spin axis (+z).To this simulation, we introduce perturber black holes on various orbits to study the effect they have on outflows and accretion.The initial projected locations of the secondaries are marked by thick 'X's.
localize failures; the values we use are based on those found to be fairly robust in GRMHD torus simulations (Porth et al. 2019).The resulting Lorentz factor in the very low density/highly magnetized jet of the simulations can directly depend on these limits and thus should not be over-interpreted.

Initial Conditions
Before adding the secondary black hole into the system, we run a single black hole simulation with a stationary Kerr metric initialized with a Fishbone & Moncrief (1976) torus with inner radius 20r g and pressure maximum at 41r g .Note that we define r g = GM 1 /c 2 in terms of the mass of the primary, which also sets the timescale r g /c = M 1 .A large, single loop of magnetic field is seeded in the torus with max P/ max P B = 100, where P B = b µ b µ /2 is the magnetic pressure.The black hole has dimensionless spin a = 0.937 in the +z direction.
The grid encompasses (1000r g ) 3 and includes a 128 3 base resolution with 8 additional levels of static mesh refinement (SMR), increasing the resolution by a factor of 2 every factor of ≈ 2 in radius.The finest resolution is concentrated within −6.25r g ≤ x, y, z ≤ 6.25r g with cell size ≈ (0.1r g ) 3 .This resolution is comparable or better than the highest resolution (192 3 ) spherical modified Kerr-Schild simulations in the Event Horizon Code Comparison Project (Porth et al. 2019) that were found to be converged for most fluid quantities.Specifically, at (r = 12r g , θ = π/2), our simulations are better resolved in most of the domain by a factor of ∼ 1.5 in the radial direction and ∼ 1.875 in the azimuthal direction but less resolved by a factor of ∼ 0.75 in the θ direction at the midplane where the modified Kerr-Schild coordinates focus the highest resolution.Within −3.125r g ≤ x, y, z ≤ 3.125r g our grid becomes comparably less resolved by a factor of ∼ 2 than the rest of our simulations because we do not place a 9th level of mesh refinement in this region (as would be required for effectively logarithmic radial spacing).We do this to save computational cost since our focus in this paper is predominantly on the flow at larger radii and there are still many cells contained within the event horizon.
We use piecewise-linear reconstruction and the HLLE Riemann solver.
The simulation is run for 100,000 M 1 to obtain the initial conditions for our binary simulations and then an additional 50,000 M 1 for comparison purposes.
Figure 1 plots a 2D poloidal density slice at this time, showing a thick, turbulent accretion flow with a narrow jet in the z direction.This flow is magnetically arrested and in equilibrium out to 70-100 r g , as shown in the dimensionless black hole flux, ϕ BH vs. time and net mass accretion rate vs. radius in Figure 2.These are defined as where u r and B r are the radial component of the fourvelocity and magnetic three-vector (converted from Cartesian to spherical CKS coordinates), dΩ = √ −gdθdφ, and the expressions are evaluated at r = 5r g 2 .

Introducing The Secondary Black Hole
At t = 100,000M 1 , we instantaneously change the metric in the simulation from that of a stationary, single black hole to the binary metric described in §2.2, starting with the initial conditions of the Post-Newtonian orbit of the secondary ( §2.3).These initial conditions are chosen such that the perturber initially is located at x = y = 0 and z = r BH,0 .This location is chosen to be within the jet so that any artifact of the sudden addition of the secondary black hole has a negligible effect on the accretion flow.
To conserve fluid quantities, after the instantaneous change in the metric we rescale each conserved variable U ∈ [ρu t , T t t , T t i ] and the magnetic field B i via 2 Ṁ and ϕ BH at r = 5r g are very similar to Ṁ and ϕ BH at the event horizon but less noisy.This ensures that the conserved energy, momentum, mass, and the divergence of the magnetic field are the same before and after the introduction of the secondary.Though we see no obvious artifact of the instantaneous introduction of the secondary in simulation quantities, we argue that even if such artifacts are present they would have a negligible effect on our results.Firstly, since we study small mass ratios, q ≪ 1, the metric only significantly changes very close to the initial location of the secondary black hole (within a few r g ), where the matter related quantities are predominantly set by the floors anyway.Secondly, the flow in this region is rapidly outflowing and so any potential artificial feature would be quickly swept away to larger radii and out of the domain of interest.
In addition to the 8 levels of SMR centered on the primary, we add another 8 levels of adaptive mesh refinement (AMR) centered on the secondary.The highest level of refinement is contained within where X 2 , Y 2 , Z 2 are the secondary's rest frame coordinates (with the secondary as the origin), the second highest level of refinement is contained within |X 2 |, |Y 2 |, |Z 2 | ≤ 6.25qM 1 , and so on.More precisely, the nth level of AMR is contained within , where n max is the maximum AMR level.This results in a cell size of (0.1r g ) 3 at the maximum refinement level (which means that the secondary horizon radius is ∼ 2 cells in length for q = 0.1 and a nonspinning secondary).As an example of how this works in our simulations, Figure 3 shows the grid structure at a representative time in our r BH,0 = 80r g , i 0 = 90 • simulation, plotted over a 2D contour of density in the y-z plane.Each 16 3 block of cells is outlined by a yellow square.This demonstrates how our grid effectively focuses resolution on the two black holes and resolves multiple scales.
Once the secondary is introduced, we run the simulations an additional 50,000 M 1 for r BH,0 = 80r g and 40,000 M 1 for r BH,0 = 20r g .This time is sufficient for ≳ 10 orbits for secondary black holes located at r BH,0 ≲ 80r g and long enough to see spin-orbit effects for orbits around r BH,0 ∼ 20r g .

Suite of Runs
The primary goals of this work are to 1) demonstrate the basic properties of a thick accretion disk around a supermassive black hole perturbed by a smaller black hole on an inclined orbit and 2) to describe, test, and demonstrate capabilities of the new time-dependent metric version of Athena++.We do not therefore seek to either simulate an exhaustive sweep of parameter space nor do we specifically focus on a target astrophysical system.Instead, we choose a select few simulations to run that we expect can represent some of the more general possibilities in such a system.Namely, we fix q = 0.1, and use two different initial black hole sepa-rations, r BH,0 = 20r g and r BH,0 = 80r g .We also use two different initial orbital inclinations, i 0 = 90 • (i.e., an orbit initially passing perpendicularly through the accretion disk), and i 0 = 45 • .These orbits are initiated as quasi-circular by using an eccentricity reduction procedure; in particular, we change the initial velocities iteratively until the quantity e = (r max −r min )/(r max +r min ) is below 0.01, where e is the eccentricity, and r min/max is the minimum/maximum distances between the two black holes.All simulations use a 2 = 0, that is, the secondary black hole is non-spinning.
We choose to focus on q = 0.1 because it is the highest mass ratio at which we feel our approximation of a stationary primary black hole is justified, while for much smaller mass ratios we have found the effects of the secondary on the primary accretion flow to be almost undetectable.We choose r BH,0 = 20r g and r BH,0 = 80r g because these represent the minimum and maximum initial separation distances at which we can reasonably trust our results.For r BH,0 < 20r g our metric approximation of a linear superposition of two boosted Kerr metrics becomes poor, while for r BH,0 > 80r g the secondary would be traveling through regions of the primary accretion flow that have not yet reached a steady state.We choose i 0 = 90 • and i 0 = 45 • to bracket the two extremes of orbits still in the "collision regime" for our thick primary accretion disk.i 0 = 90 • orbits are completely perpendicular to the disk and thus the impact velocity of the secondary is  and i 0 = 45 • (black, dashed).Bottom: r BH,0 = 80r g for i 0 = 90 • (red, solid) and i 0 = 45 • (black, dashed), where the inset zooms in on the time covered by our simulations.For orbits with initial separations of r BH,0 = 20r g , the timescale for the orbital radius to significantly change is ∼ 10 4 M 1 , much shorter than the ∼ 10 6 M 1 timescale for the r BH,0 = 80r g orbits.Since 10 6 M 1 is much longer than the length of our simulations, for r BH,0 = 80r g we safely assume Keplerian circular orbits that do not change in time.For r BH,0 = 20r g simulations, however, we include the full orbital evolution and can expect to see significant changes on the timescales we can simulate.The differences between i 0 = 90 • and i 0 = 45 • binary separation distances are relatively small until close to merger (∼ 4 × 10 4 for r BH,0 = 20r g and ∼ 9 × 10 6 for r BH,0 = 80r g ).maximized.i 0 = 45 • orbits on the other hand only graze the edge of the disk with much smaller tangential velocity (as mentioned in §3.2).
The resulting binary separation distances as a function of time given by solving the PN equations using CBwaves for these 4 different orbits are plotted in Figure 4.For r BH,0 = 20r g , merger would happen after ≈ 4.5-5.5 ×10 4 M 1 (earlier for i 0 = 90 • , later for i 0 = 45 • ), with significant changes to Figure 5. Angles and distances describing the orbits in our simulations.r BH is the separation between the primary and secondary black holes, a 1 is the primary black hole spin, l 2 is the secondary's orbital angular momentum, i is the angle between the x-y plane and the orbital plane (where i = 0 corresponds to clockwise motion in the x-y plane), and θ a is the angle between the primary black hole spin direction and the z-axis (the initial primary black hole spin direction).
the binary separation happening on timescales of ∼ 10 4 M 1 (note that these are comparable to the 4 × 10 4 M 1 runtime of our simulations).For r BH,0 = 80r g , merger would happen after ≈ 10.3-10.7 ×10 6 M 1 (earlier for i 0 = 90 • , later for i 0 = 45 • ), with significant changes to the binary separation happening on timescales of ∼ 10 6 M 1 (note that these are much longer than the 4 × 10 4 M 1 runtime of our simulations).The same timescales are seen in the evolution of the orbital angular momentum of the secondary and primary black hole spin, which we quantify using the inclination of the orbit, i, as well as the angle between the primary black hole spin and its initial direction along the z-axis, θ a .These are defined as and where a 1,i is the spin vector of the primary, a 1 = (a 1,i a i 1 ) 1/2 , l 2,i is the specific angular momentum of the secondary (e.g., 2 ) 1/2 , x 2 , y 2 , and z 2 are the x, y, and z positions of the secondary, and v 2,i are the velocities of the secondary.These angles are diagrammed schematically in Figure 5.
The angles i and θ a are plotted in Figure 6 for all four orbital configurations.For r BH,0 = 80r g , the direction of the primary spin of the black hole changes by ≈ 45 • for i 0 = 45 • in the first ∼ 10 6 M 1 and ≈ 90 • for i 0 = 90 • in the first ∼ 1.5 × 10 6 M 1 , while the direction of the orbital angular momentum changes similar amounts during the same times.For r BH,0 = 20r g , the direction of the primary spin of the black hole changes by ≈ 30 • for i 0 = 45 • in the first ∼ 10 4 M 1 and ≈ 55 • for i 0 = 90 • in the first ∼ 1.5 × 10 4 M 1 , while the direction of the orbital angular momentum changes by ≈ 60 • for i 0 = 45 • and ≈ 125 • for i 0 = 90 • during the same times.The orbital and spin directions then continue to oscillate back and forth from the initial values on shorter and shorter timescales until merger (which happens after 3-6 oscillations for r BH,0 = 20r g and ∼ 20 oscillations for r BH,0 = 80r g ).
Since we can only reasonably simulate timescales ≲ 10 5 M 1 , for the r BH,0 = 80r g simulations we neglect orbital changes and assume circular Keplerian orbits, evolving for ≈ 9 orbits.For r BH,0 = 20r g , however, we could in principle simulate all the way to merger, although at that point the approximation used in superimposing the two black hole metrics without any interaction terms would break down.Instead, for r BH,0 = 20r g we simulate up to separations of ∼ 14r g (∼ 85 orbits) and ∼ 12r g (∼ 79 orbits) for i 0 = 45 • and i 0 = 90 • , respectively, so we see significant changes in both the primary black hole spin and orbital directions throughout our simulations.
We emphasize that the orbital and spin evolutions used in our simulations that we have just described do not include any fluid effects like drag or dynamical friction.Instead, they represent the solution to the Post-Newtonian orbital equations for two black holes in a vacuum.This is a good approximation for moderate and low accretion rates, but fluid effects could become important for the highest accretion rates (close to and exceeding Eddington for either the primary or the secondary), at least for larger separations (e.g., our 80r g case) where there is significant time for accumulated drag and friction to affect the secondary before merger.Incorporating such effects in the orbital evolution of the binary would require coupling the PN orbital equations to the GRMHD simulation via additional source terms that account for the accretion of linear/angular momentum and the gravitational effects of the stress-energy tensor of the surrounding plasma.Naively one might expect drag and friction to reduce the orbital velocity of the secondary (and thus increasing the accretion rate onto the secondary) and reduce the time it takes for the binary to merge.On the other hand, previous studies have found that flows around compact objects with significant outflows can have negative dynamical friction (Gruzinov, Levin & Matzner 2020; Li et al. 2020;Kaaz et al. 2023).In that case the orbital velocity of the secondary and the time it takes for the binary to merge could increase.The cumulative long term effect that drag and friction have on binary evolution is still controversial and requires further numerical study.We also note the additional challenge that this fluid back-reaction is likely only significant on timescales much longer than the orbital time during which the secondary can sample several different regions and realizations of the flow.

BLACK-HOLE DISK INTERACTION: ANALYTICAL CONSIDERATIONS
The passage of a smaller, secondary black hole through a thick accretion disk surrounding a more massive primary black hole can be compared to a Bondi-Hoyle-Lyttleton-type flow where a uniform wind impacts a massive object (Hoyle & Lyttleton 1939;Bondi & Hoyle 1944;Edgar 2004; quite different from a BH-disk collision in thin cooled disks cf.Ivanov, Igumenshchev & Novikov 1998).This is true specifically for inclined orbits on small spatial scales close to the secondary and short time-scales where the orbital motion of the wind can be approximated as linear.
It is thus instructive to consider that solution in the context of our simulations.Doing so allows us to get rough estimates of what we expect to find in the simulations and gives us a conceptual framework to interpret our results.As the simulations confirm the basic paradigm described by the model, it also allows us to make predictions about regions of parameter space that we have not simulated.
For the purposes of this section, we use the variable r for the radial distance away from the primary and R for the radial distance from the secondary.

90 Degree Inclined Orbits
Assuming that the secondary is on a fixed circular orbit with inclination i = 90 • around the primary at r BH , then the asymptotic impact velocity with respect to the secondary black hole is , where α kep is the rotational velocity of the accretion disk divided by the Keplerian velocity (α kep ≈ 0.5 for radiatively inefficient MADs, see Narayan et al. 2012;Ressler, Quataert & Stone 2020), and the asymptotic impact density is ρ(r = r BH ) ≈ ρ H (r H /r) 0.8 , where r H is the event horizon radius of the primary, ρ H ≡ ρ(r = r H ), and we have assumed that the density scales as r −0.8 in the radial range of interest (Xu 2023; we will show later that this is a good assumption in our simulations).For simplicity we also have taken the density to be independent of angle.The asymptotic sound speed is expected to be some fixed fraction of the Keplerian velocity, which we measure to be ≈ 0.3 for r ≳ 10r g in our simulations (c 2 s,∞ ≈ 0.3GM 1 /r BH ).The accretion radius is then This can be compared with the Hill radius (inside of which the gravity of the secondary dominates over the gravity of the primary): R Hill ≈ {q/[3(1 + q)]} 1/3 r BH .For these parameters, when q ≲ 0.34, R BHL < R Hill and thus the effective influence radius of the secondary is determined by R BHL .Gas outside this radius will be relatively unaffected by the secondary black hole while gas inside this radius will be accreted.The Figure 6.Angle between the primary black hole spin vector and its initial direction, θ a (green solid), and the angle between the orbital angular momentum vector and the x-y plane, i (black dashed), vs. time for the orbital choices in our simulations (all with q = 0.1).Top left: r BH,0 = 20r g , i 0 = 90 • .Top right: r BH,0 = 80r g , i 0 = 90 • .Bottom left: r BH,0 = 20r g , i 0 = 45 • .Bottom right: r BH,0 = 80r g , i 0 = 45 • .In all simulations, both the primary black hole spin and orbital angular momentum change directions significantly over time.i changes by 45 • -130 • while θ a peaks at around 30 • -90 • .Generally, the r BH,0 = 20r g orbits have stronger amplitude variation in i, while the r BH,0 = 80r g orbits have stronger amplitude variation in θ a .The timescales for these changes are the same as for the binary separation distances shown in Figure 4 so that we expect to see significant change in primary spin and orbital angular momentum in r BH,0 = 20r g simulations but not in r BH,0 = 80r g simulations due to the run times.For comparison, the insets in the r BH = 80r g plots show θ a and i zoomed in on the timescale of our simulations.approximate accretion rate onto the secondary is (30) so we expect the secondary's size of influence and accretion rate to increase with orbital radius.
The accretion disk is not spherically symmetric, however, but has magnetically dominated, matter-deficient polar regions.For a rapidly spinning black hole in a MAD state as we study here, there will also be a powerful, electromagnetically dominated jet pushing outwards.As the secondary passes through the disk into the pole, we might expect it to bring with it the amount of mass contained within R ≤ R BHL .This will be ∝ q 3 r 2.2 BH .
(31) Now this mass will be deposited into the jet region in a time so that the passage of the secondary from the disk into the polar region should increase the unbound outflow rate (as-suming it is outside of the stagnation surface) by .
(33) Now, we can compare this with the expected scaling of the unbound outflow rate for the accretion disk itself (i.e., material blown off the disk in the process of accretion, not the highly relativistic jet material).In radiatively inefficient flows with significant outflows, the mass inflow ( Ṁin ) and outflow rates will be roughly equal in magnitude with the inflow speeds being some fraction of the Keplerian speed v kep (Begelman 2012).Then Ṁunbound ≈ | Ṁin | ∝ ρv kep r 2 ∝ r 0.7 , which scales the same way with r as ∆ Ṁunbound scales with r BH .We have confirmed that this scaling holds in the single black hole simulation described in §2.7 for r ≳ 10r g (not plotted).Thus we expect that the ratio between ∆ Ṁunbound and Ṁunbound will be similar if measured at r = r BH for all r BH .
To estimate the impact the secondary might have on the primary accretion flow, as an upper limit we can think of the secondary as effectively screening a fraction of the inflowing material, determined by the area that it sweeps out in the disk over the course of its orbit on a spherical shell located at r = r BH .This area is where H is the scale height of the disk and we have used R BHL ≪ r BH .The effective area of the inflowing accretion disk is similarly A rough estimate of the amount by which the net accretion rate could change is then Therefore, for i 0 = 90 • orbits we expect the secondary to have a minimal effect on the net accretion flow of the primary (as we show later).
In this brief analysis we have neglected many considerations that might be important in the simulations, including magnetic fields (which can change the structure of the Bondi-Hoyle-Lyttleton accretion flow, Kaaz et al. 2023; Gracia-Linares & Guzmán 2023), the velocity gradient in the wind provided by the accretion disk (which can induce turbulence and also change the structure of the flow, Xu & Stone 2019), turbulence (which could introduce stochastic variability to the predicted quantities), the time-dependent nature of accretion (which could introduce secular variability to the predicted quantities), and the variation of density with angle (which could lead to smaller-than-predicted mass outbursts since the density on the surface of the disk is smaller than the midplane).These approximate values, however, give us a good set of comparisons for our numerical results.

More General Expressions
The above analysis can also be done for orbital planes closer to the midplane of the disk.This will have the effect of either increasing or decreasing v ∞ in the frame of the secondary depending on whether the orbit is prograde or retrograde to the accretion disk.It will also increase the time it takes to deposit matter outside the disk (for orbits sufficiently inclined that the secondary still passes out of the disk) because the component of the velocity perpendicular to the disk, v perp , will be reduced.Both v ∞ and v perp depend on the particular location of the secondary in its orbit when it crosses the surface of the disk.However, we can approximately evaluate them when the secondary crosses the midplane as } and v perp = sin(i)v orbit .These expressions are approximately valid if the accretion disk is not too thick (≲ 30 • above and below the midplane).We can also parameterize the sound speed of the disk as c 2 s,∞ = α 2 s GM 1 /r BH and the rotational velocity of the disk as v 2 rot = α 2 kep GM 1 /r BH .Repeating the same calculation as in the previous subsection, this results in and Substituting in i 0 = 45 • , α kep = 0.5, and α 2 s = 0.3, we find that R BHL , ṀBHL , and ∆ Ṁunbound are larger than the i 0 = 90 • expressions by factors of ≈ 1.8, 2.5, and 2.4, respectively.Note, however, that Equation (39) for ∆ Ṁunbound crucially depends on the assumption that the orbit of the secondary brings it out of the disk into the polar region.If the disk is too thick or the orbit not inclined enough the actual value of ∆ Ṁunbound will be much less.This is true in our simulations for i 0 = 45 • (note the thickness of the disk in Figure 1), where the orbit only grazes the edge of the disk instead of plunging out into the polar region.Thus we might expect the relative ∆ Ṁunbound to be smaller, though it is not obvious by how much.Note additionally that for low inclinations R Hill can become less than R BHL .For these parameters, at q = 0.1 this happens for −24 • ≲ i ≲ 24 • (meaning that for all the simulations in this work, the Bondi-Hoyle-Lyttleton radius determines the influence radius).
The general expression for the area swept out in the disk by the orbit of the secondary on a spherical shell located at r BH is where the top expression is used when the orbit of the secondary passes out of the disk at some point in its orbit and the second expression is used when the orbit is entirely contained within the disk.For the thick disks we study in this work, the term related to the scale height is relatively close to 1, so for all orbits the maximum |∆ Ṁ/ Ṁ| predicted by Equation ( 41) is ≈ 0.6q ≪ 1 .
It is interesting to note that Equation (41) predicts that for a thinner disk with H/r BH ≪ 1, secondaries with significantly inclined orbits would still have a relatively small effect on the primary accretion disk.On the other hand, for orbits with low inclination, 2π| sin(i)| ≲ R BHL /r BH , Equation (41) predicts |∆ Ṁ/ Ṁ| ∼ 1, at which point the disk structure would likely significantly change and this approximation would break down.If we assume α kep ≈ 1 and α s ≈ 0 for a thin disk, then this would happen for −18 • ≲ i ≲ 18 • when q = 0.1.However, we again emphasize that this estimate is simplistic and neglects thermal effects that could be significant, especially for thin disks.The problem of a secondary black hole impacting a thin accretion disk has also been studied analytically with significantly more detail in other works (e.g., Ivanov, Igumenshchev & Novikov 1998;Ivanov, Papaloizou & Polnarev 1999;Pihajoki 2016).Our numerical method will be able to study such systems when combined with optically thin radiative cooling or a full treatment of radiation (e.g., White et al. 2023).

GRMHD RESULTS
To highlight the basic features of our binary simulations, Figure 7 shows the time evolution of mass density in our r BH,0 = 80r g , i 0 = 90 • simulation.The secondary black hole travels supersonically through the accretion disk, forming a bow shock that propagates through the flow.As the secondary continues along its orbit and crosses out of the disk into the jet region, it carries with it a small amount of matter that gets deposited into the funnel region and subsequently blown away by the jet.Evidence for the fact that the gas is being accelerated by the jet is in the fact that the time-averaged electromagnetic outflow energy (e.g., the Poynting flux, not shown) is reduced in the binary simulations when compared with the single black hole simulation3 .This process continues in a periodic or quasi-periodic way for the duration of the simulation.
Figure 8 zooms in on the secondary black holes in the lab frame as they pass through the midplane for our four simulations, displaying 2D contours of the square sound speed normalized to the square Keplerian velocity, c 2 s /v 2 kep (∝ T g r, where T g is the gas temperature), in the x = x BH plane (parallel to y-z) at representative times for our four binary simulations (except for the r BH,0 = 20r g , i 0 = 45 • simulation which is plotted in the y = y BH plane parallel to x-z to better display the motion of the secondary through the gas).Around the secondary, the flow resembles a turbulent Bondi-Hoyle-Lyttleton-type flow with a shock front and wake.The shocks are stronger for orbits with i 0 = 90 • than i 0 = 45 • due to the reduced relative motion of the gas in the latter, but all simulations have moderate shocks with temperature increasing by factors of ∼ a few from the average temperature at the orbital radius.These temperatures are consistent with Bondi-Hoyle-Lyttleton simulations for flows with Mach numbers moderately greater than 1 and ≲ 2. Another consequence of the reduced relative velocities of the gas is that the i 0 = 45 • simulations have larger secondary influence radii than their i 0 = 90 • counterparts.The sizes of the influence radii of the secondaries in the r BH,0 = 80r g simulations are also about 4 times larger than those in the r BH,0 = 20r g simulations due to the slower orbital velocity.

Effects on the Primary Accretion Flow
In Figure 7 neither the shock from the secondary nor the ejection of matter from the disk seem to dramatically affect the accretion flow dynamics onto the primary.To quantify this more directly, we plot the accretion rate onto the primary black hole and dimensionless magnetic flux threading the primary black hole vs. time in Figure 9  the binary simulations compared with the single black hole simulation.These quantities are calculated in the same way as in Equations ( 24) and ( 25), that is, using a single black hole Kerr metric (which is a good approximation near the primary event horizon).The quantities for the binary simulation show statistically almost identical behavior to those in the single black hole simulation.[13,13,14,14,13] for the [(r BH,0 = 80r g i 0 = 90 • ), (r BH,0 = 80r g i 0 = 45 • ), (r BH,0 = 20r g i 0 = 90 • ), (r BH,0 = 20r g i 0 = 45 • ), single black hole] simulations, which display only slight differences.In Figure 9 there are also no clear signatures of the periodicity of the secondary orbit.
The time and angle-averaged radial profiles in the binary simulations are also quite similar to those in the single black hole simulation as seen in Figure 10 for the density and square sound speed, where we perform time averages over the range 120,000-140,000 M 1 .For radii near and within the r in these units.The contours are centered on the secondary black holes as they pass through the midplane of the accretion disk during representative times in our simulations.The white circles represent the influence radii predicted from Equation (37).The moving black hole creates a bow shock in the gas of varying shapes and sizes.All simulations generally have moderate shocks with temperatures (note that c 2 s is directly proportional to temperature) increasing by factors of ∼ a few from the average gas temperature at the orbital radius.The i 0 = 45 • (right) simulations generally display a larger secondary influence radii but slightly weaker shocks than the i 0 = 90 • simulations (left) due to the reduced relative gas motion.Similarly, the r BH,0 = 80r g simulations (top) display a larger secondary influence radius than the r BH,0 = 20r g simulations (bottom).All simulations show good agreement with the predicted influence radii from Equation (37).
orbital radius of the secondary, the density is decreased relative to the single black hole simulation by ∼ 20%.This is caused by a combination of matter being expelled from the disk by the secondary and by matter being accreted onto the secondary.The density in all simulations agrees reasonably well with an r ∝ r −0.8 dependence between 3r g ≲ r ≲ 80r g as used for our analytic estimates in §3 (Xu 2023).The temperatures (directly proportional to c 2 s ) of the binary simulations are slightly hotter (by ∼ 20%) in the bulk of the disk for r ≳ 10r g caused by the bow shock propagating through the flow.This is particularly evident near the orbital radii where there are small peaks in temperature.The r BH,0 = 80r g , i 0 = 45 • simulation has an especially prominent peak at the orbital radius because it spends a large fraction of time within the disk and so the shocked temperature contributes more to the time-averaged temperature at that radius.
These findings agree with our analytical estimates in §3.2, particularly Equations ( 36) and ( 41), where we argued that the effect of the secondary on the primary accretion flow would be quite small for all orbital inclinations.

Quasi-Periodic Outbursts/Eruptions
Even though the primary accretion flow dynamics are not significantly affected by the secondary black hole, there are, .Accretion rate, Ṁ, and dimensionless magnetic flux threading the primary black hole, ϕ BH , as a function of time for the single black hole simulation (black dotted line) compared with the q = 0.1, r BH,0 = 80r g , i 0 = 90 • (solid blue line) and q = 0.1, r BH,0 = 20r g , i 0 = 90 • (dashed red line) binary simulations.Neither the accretion rate nor the magnetic flux show any clear signature of the secondary orbit's periodicity.In fact, both quantities are very similar in the binary and single black hole simulations, with mainly stochastic differences.This shows that the effect of the secondary is small on the near-horizon dynamics of the flow.We argue that this is because the area swept out in the disk by the orbit of the secondary on a spherical shell is small compared with the area of the disk itself at the orbital radius.Note that Ṁ and ϕ BH for the i 0 = 45 • simulations are similar and are not plotted to avoid clutter.
however, clear signatures of the secondary in the outflow.Zooming in on times when the secondary passes out of the disk into the polar regions, we plot 2D contours of mass density in the x = x BH plane (parallel to y-z) centered on the secondary black hole as it crosses into the jet in Figure 11 for two particular representative time series in the r BH,0 = 80r g , i 0 = 90 • and r BH,0 = 20r g , i 0 = 90 • simulations.In both simulations, as the secondary passes from the disk/jet boundary region to the jet itself, it brings with it a blob that subse- .Time and angle-averaged radial profiles of density multiplied by r 0.8 (top set of lines), ⟨ρ⟩r 0.8 , and mass-weighted square sound speed divided by the square Keplerian speed (bottom set of lines), ⟨c 2 s /v 2 kep ⟩ ρ for our four different binary simulations compared with the single black hole simulation (black solid line).The secondary black hole has only a small effect on the average density and temperature, providing a small increase in temperature around the orbital radius (caused by the bow shock following the secondary) and a small decrease in density near and at smaller radii than orbital radius (caused by the expulsion of mass from the disk).
quently expands, shears out, and gets blown away by the jet.The sizes of the initial blobs are consistent with the Bondi-Hoyle-Lyttleton radius (Equation 29), roughly 10r g and 2.5r g in radius for r BH,0 = 80r g and r BH,0 = 20r g when mapped to a sphere.The temperatures of the blobs are approximately virial, with c s ≈ v kep .
To measure the effect of these ejected blobs on the outflow, we particularly focus on the unbound mass outflow rate, defined as where u r is the radial component of the four velocity, converted from CKS coordinates in the rest frame of the primary, is the relativistic Bernoulli parameter (where Be > 0 implies unbound material, Penna, Kulkarni & Narayan 2013), dΩ = √ −g sin(θ)dθdφ, and θ and φ are the polar and azimuthal angles converted from CKS coordinates in the rest frame of the primary.Ṁunbound is shown as a function of time in Figure 12 for our binary simulations measured at r = 2r BH,0 , compared to the same quantities in the single black hole simulation measured at the same radii.The unbound outflow rates display clear quasi-periodic peaks on timescales comparable to an orbital period (≈ 4500 M 1 for r BH,0 = 80r g and ≈ 560 M 1 for r BH,0 = 20r g ), with durations ∼ 2000 M 1 for r BH,0 = 80r g and ∼ 200 M 1 (or ∼ 1/2 an orbital period for each).Each "outburst" is of varying intensity when compared to the single black hole simulation.For some of the peaks, Ṁunbound is only a few percent higher than the corresponding Ṁunbound in the single black hole simulation, while at others it reaches ∼ 2-4 times Ṁunbound in the single black hole simulation.The fact that the peaks for both r BH,0 values are about the same magnitude relative to the single black hole simulation is expected based on our analysis in §3.Also consistent with our analysis is that the absolute magnitude of the peaks in Ṁunbound are larger for r BH,0 = 80r g than for r BH,0 = 20r g (when measured at the same multiple of r BH,0 ).The former peaks are larger by a factor of ∼ 2-3, which is in good agreement with our estimate (Equation 39) of ∼ 4 0.7 ≈ 2.6.The peaks in the i 0 = 45 .Unbound mass outflow rate, Ṁunbound , plotted vs. time elapsed since the secondary was introduced, t − t 0 , in our r BH,0 = 80r g (top) and r BH,0 = 20r g (bottom) simulations for i 0 = 90 • (blue solid line) and i 0 = 45 • (green dashed line), compared to the single black hole simulation (black dotted line).Ṁunbound is measured at r = 160r g for r BH,0 = 80r g and r = 40r g for r BH,0 = 20r g .The curves have been normalized by the time-averaged accretion rate over the interval (0 ≤ t − t 0 ≤ 40,000M 1 ).The secondary increases the unbound mass outflow rate at quasi-period intervals of order the orbital period (≈ 4500 M 1 for r BH,0 = 80r g and ≈ 560 M 1 for r BH,0 = 20r g ) by factors of ∼ 2-4 compared to the same quantity in the single black hole simulation.Note that as in all GRMHD torus simulations there is a secular decline in all outflow rates due to the mass in the torus depleting over time.
generally about the same magnitude or slightly lower on average than those in the i 0 = 90 • simulations.This is because the orbits for i 0 = 45 • do not bring the secondary fully out of the disk into the polar regions and are thus less efficient at depositing matter there (even though their influence radii are larger due to the reduced relative gas speed, see §3.2).
To more quantitatively measure the frequencies of the outbursts in unbound outflow rates we compute power spec- .Power spectra of unbound mass outflow rates, Ṁunbound , for orbits with i 0 = 90 • (blue solid) and i 0 = 45 • (green dashed) compared to the single black hole simulation (black dotted).Each spectrum is normalized such that the highest peak is 1.Top: r BH,0 = 80r g , with Ṁunbound measured at r = 160r g and the x-axis measuring frequency in units of the orbital frequency at r = 80r g .Bottom: r BH,0 = 20r g , with Ṁunbound measured at r = 40r g and the x-axis measuring frequency in units of the orbital frequency at r = 20r g .For r BH,0 = 80r g , the i 0 = 90 • PSD shows clear peaks at the orbital frequency and twice the orbital frequency, with other smaller peaks.The i = 45 • does have a peak at the orbital frequency but it is subdominant when compared with the lower frequency peak also seen in the single black hole simulation.For r BH,0 = 20r g , there are significant peaks at the orbital frequency for both i 0 = 90 • and i 0 = 45 • , but these are lower than some peaks at lower frequencies.This is likely due to propagation effects that spread out mass ejected by the secondaries in radius (and correspondingly, time in Ṁunbound ).tra of the Ṁunbound curves by first linearly de-trending the data and then using Welch's method, which averages periodograms for overlapping windows of the data.We plot the resulting power spectra in Figure 13 using a window size of 20, 000M 1 and scale the frequency to the orbital frequency, f kep (r = r BH,0 ).Due to the limited sample size of the data and the complex variability of the primary MAD accretion disk, the power spectra for the five simulations display a number of peaks at different frequencies, many of which are sensitive to the precise method/averaging used to compute the periodogram.Because of this, specific features in Figure 13 should not be over-interpreted, rather, our focus is on the general behavior.Ṁunbound for the single black hole simulation generally shows the most power at lower frequencies (e.g., periods ∼ 17,000 M 1 at r = 160r g and ∼ 8000 M 1 at r = 40r g ).At 2r BH,0 , The binary simulations generally show a peak at the orbital frequency of the binary and sometimes twice that frequency (i.e., every half orbit), with the cleanest example being r BH,0 = 80r g , i 0 = 90 • , where the two highest peaks are located at f kep (r = 80r g ) and 2 f kep (r = 80r g ).The r BH,0 = 80r g , i 0 = 45 • power spectrum also shows a peak at f kep (r = 80r g ), but it is subdominant compared to the lower frequency peak also seen in the single black hole simulation, likely due to the fact that the orbit does not bring the secondary fully out of the disk to create as distinctive outbursts as i 0 = 90 • .The r BH,0 = 20r g simulations show a diverse set of frequencies that stand out in the Ṁunbound power spectrum.For i 0 = 90 • the orbital frequency is the third highest peak, with the highest located at a little over half the orbital frequency (or a period of ∼ 1000 M 1 ), while for i 0 = 45 • the orbital frequency is the second highest peak, with the highest located at ≈ 0.4 times the orbital frequency (or a period of ∼ 1400 M 1 ).Both simulations have several other prominent peaks located at lower frequencies.This is at least partly due to propagation effects.As the secondary brings matter into the polar regions, the jet accelerates and unbinds this matter and it is propelled to larger radii.The gas then can spread out in radius (resulting in peaks in Ṁunbound with a broader spread in time) and even catch up with previously ejected matter (resulting in merged peaks in Ṁunbound showing up in lower frequencies in the power spectrum).
The unbound outflow rates we have studied thus far have been measured at r = 2r BH,0 .Perhaps a more observationally meaningful measurement of outflow is to measure Ṁunbound at infinity or very large distances from the primary black hole.Practically, we have only a limited range of radial distances included in our simulation (r ≲ 1600r g ), and so use r = 5r BH,0 as a proxy for larger radii.Ṁunbound at r = 5r BH,0 and the associated power spectra are shown for the i 0 = 90 • simulations in Figure 14 (again compared to the single black hole simulation values at each radius).When averaged over time, the unbound outflow rates in the binary simulations are now 20-30% larger than the single black hole unbound outflow rates at the same radii.For r BH,0 = 80r g the peaks seen in Figure 12  r BH,0 = 80r g , i 0 = 90°r BH,0 = 20r g , i 0 = 90°F igure 14.Top: Unbound mass outflow rate, Ṁunbound , plotted vs. time elapsed since the secondary was introduced, t − t 0 , in our r BH,0 = 80r g , i 0 = 90 • (red solid line) and r BH,0 = 20r g , i 0 = 90 • (blue dashed line) simulations measured at r = 5r BH,0 compared to the single black hole simulation measured at the same radii (red and blue dotted lines).Ṁunbound normalized by the time-averaged accretion rate over the interval (0 ≤ t − t 0 ≤ 40,000M 1 ).Bottom: Power spectrum for Ṁunbound in the r BH,0 = 80r g , i 0 = 90 • and r BH,0 = 20r g , i 0 = 90 • simulations, each normalized such that the highest peak is 1.Note that the frequency resolution compared with the orbital frequency at r BH,0 , f kep (r = r BH,0 ), is much higher for r BH,0 = 20r g than r BH,0 = 80r g , resulting in smoother and broader peaks for the r BH,0 = 80r g power spectrum.Compared with Figures 12 and 13, power has shifted to lower frequencies than the orbital frequencies at r BH,0 .Peaks in Ṁunbound have become more broad in time.We propose that this is due to propagation effects as the mass outflow provided by the secondary spreads out in radius as it is accelerated outwards by the jet.For r BH,0 = 80r g this results in a more continuous increase in unbound outflow rates compared with the single black hole simulation, while for r BH,0 = 20r g there are still distinct, relatively high-contrast peaks.As in all GRMHD torus simulations, there is a secular decline in all outflow rates due to the mass in the torus depleting over time.all times).For r BH,0 = 20r g , on the other hand, there are still peaks with relatively high contrast (factors of 2-3), with much of the Ṁunbound curve lying close to the single black hole Ṁunbound curve.In terms of the power spectra, this behavior corresponds to a shift in power to lower frequencies.The r BH,0 = 80r g power spectrum still shows noticeable peaks at f kep (r = r BH,0 ) and 2 f kep (r = r BH,0 ), but they are shorter than the peak at low frequencies.The r BH,0 = 20r g power spectrum has almost no peak at f kep (r = r BH,0 ) but instead has several peaks ≲ 0.5 f kep (r = r BH,0 ) (or periods ≳ 1100M 1 ).Again, we interpret this behavior as the dispersion of the unbound matter provided by the secondary as it travels outwards in radius.Indeed, starting from r BH,0 and plotting Ṁunbound at progressively larger radii (not shown) we see that the peaks spread out in time and even merge together.Eventually the quasi-periodicity in the unbound outflow rate disappears entirely for r ≳ 10r BH,0 (also not shown), with the power spectra shifting to the lowest frequencies.The resulting Ṁuunbound is then a more continuous 20-30% increase above the "background" unbound outflow rate.This means that we expect quasi-periodic signatures in the outflow to be only significant for a limited range of radii that depend on the orbital period: r BH,0 ≲ r ≲ 10r BH,0 .
We note that the fact that the peaks in the power spectra for r BH,0 = 80r g are broader and smoother in Figure 14 than those for r BH,0 = 20r g is due to the fact that f kep (r = 80r g ) is only ∼ 4 times the frequency resolution (i.e., the orbital period, ∼ 4500M 1 , is ∼ 1/4 the window size of 20,000M 1 ).

Spin-Orbit Coupling
In this subsection, we highlight the features of the r BH,0 = 20r g simulations where the spin-orbit coupling between the binary orbit and the primary black hole spin is evident on simulated timescales.For instance, Figure 15 shows a volume rendering of the accretion flow at five different times and two different spatial scales for the r BH,0 = 20r g , i 0 = 90 • simulation.Regions with high σ = b 2 /ρ and regions with high ρ are highlighted with green/blue and red, respectively.Initially, when the perturber is first introduced, the accretion flow, electromagnetically dominated jet, and the primary black hole spin axis are all aligned.As the primary black hole spin tilts at later times due to relativistic spin-orbit coupling with the perturber, both the accretion flow and the jet adjust so that near the peak tilt (see Figure 6) the jet and disk angular momentum are mostly aligned with the new spin axis.As the spin returns to its original orientation and this process repeats, the disk and jet orientation on these scales for the most part trace the black hole spin axis (especially at small radii).This is because the primary black hole spin axis evolves on relatively long timescales (see Figure 6) so that at these radii the flow can adjust approximately adiabatically.To see this more quantitatively, Figure 16 compares the tilt of the primary black hole spin axis, θ a , to the tilt angle of the accretion flow as a function of time for the two r BH,0 = 20r g simulations at r = 5r g , r = 50r g , and r = 120r g .We define this tilt angle with respect to the z-direction (i.e., the original spin axis before the perturber is introduced) as (e.g., White, Quataert & Blaes 2019) where and ⟨⟩ denotes an angle average.In Figure 16, θ tilt at r = 5r g for both simulations displays significant stochastic temporal variability but on average follows the θ a curve, i.e., the disk angular momentum axis quickly aligns with the primary black hole spin axis.This corresponds to a peak average tilt of ∼ 60 • for i 0 = 90 • and ∼ 30 • for i 0 = 45 • , though the i 0 = 45 • simulation has times with larger θ tilt (e.g., ∼ 40 • ).At this radius the angular momentum of the gas varies by as much as ≈ 20 • on short (≲ 1000M 1 ) timescales.At larger radii the gas angular momentum is not only less variable but slower to respond to the change in the primary black hole axis.In both simulations the angular momentum direction at r = 50r g and r = 120r g first tilts along with the primary black hole spin axis but at roughly half the rate of the disk at smaller radii.When the primary black hole spin axis starts returning to its original value, however, the gas does not similarly return to its original orientation.Instead, it remains tilted for the rest of the simulations, effectively saturating at θ tilt ∼ 30 • for i 0 = 90 • and θ tilt ∼ 15 • for i 0 = 45 • .This is likely because the timescale for the primary black hole spin to change is smaller than the timescale for the larger radii flow to align and so the larger radii flow effectively sees a time-averaged black hole spin direction (∼ 30 • from the z-axis for i 0 = 90 • and ∼ 15 • from the z-axis for i 0 = 45 • ).
We perform a similar analysis with the jet in Figure 17, which shows how the tilt angle of the jet, θ jet , changes as a function of time for r = 10r g , r = 100r g , and r = 800r g compared to θ a .Here we define θ jet by measuring the angle between the z-axis and the σ-weighted position vector of either the upper or lower jet.For simplicity, we only plot θ jet for the upper jet in Figure 17, but the quantity looks similar for the lower jet.Generally, for both simulations the behavior of θ jet is similar to θ tilt at r = 5r g , that is, the jet direction mostly follows the black hole spin axis with added stochastic variability.The jet at larger radii (r = 100r g and r = 800r g ) like the disk also similarly lags behind the black hole spin axis initially, with a slower change in direction than at smaller radii (though to a lesser extent than the disk, with peak tilt angles around ∼ 40-50 • and ∼ 15-20 • for i 0 = 90 • and i 0 = 45 • , respectively).Unlike the disk, however, the jet at large radii does return to the initial orientation and fol- r BH,0 = 20r g i 0 = 45°F igure 16.Angle between the angle-averaged angular momentum of the gas and its initial direction, θ tilt , vs. time measured at r = 5r g (solid gold), r = 50r g (dashed purple), and r = 120r g (dotted teal), compared to the tilt angle of the primary black hole spin axis, θ a (solid black), for the r BH,0 = 20r g simulations with i 0 = 90 • (top) and i 0 = 45 • (bottom).As the black hole spin axis changes due to spin-orbit coupling (see Figure 6), frame-dragging and torques from magnetic fields cause the flow quickly to align with the gas at small radii r ≈ 5r g so that θ tilt (r = 5r g ) mostly traces θ a with added stochasticity.At larger radii (r ≈ 50r g and above), however, the gas tilt slowly rises to a saturation value of ∼ 30 • and ∼ 15 • for i 0 = 90 • and i = 45 • , respectively, and never returns back to θ tilt = 0.This is likely because the flow at these radii effectively sees the time-averaged central black hole spin.low the black hole spin axis, at least in part.This is likely because changes in the jet propagate at roughly the speed of light which is ≫ than the average radial velocity in the bulk of the disk The fact that the secondary causes such a dramatic change in the orientation of the accretion flow and jet in these two r BH,0 = 20r g simulations makes it all the more remarkable that their accretion rates and dimensionless black hole flux values were so similar to the single black hole simulation in r BH,0 = 20r g i 0 = 45°F igure 17.Angle between the electromagnetically dominated upper jet (that is, the jet with polar angle θ < π/2) and its initial direction, θ tilt , vs. time measured at r = 10r g (solid gold), r = 100r g (dashed purple), and r = 800r g (dotted teal), compared with the primary black hole spin tilt angle, θ a (solid black), for the r BH,0 = 20r g simulations with i 0 = 90 • (top) and i 0 = 45 • (bottom).The jets at all radii tend to align with the black hole spin but at a rate that decreases with increasing distance from the central black hole.Unlike the disk (Figure 16), the jet does tend to return to θ jet = 0 with the black hole spin due to the shorter timescales associated with the relativistic outflow speeds.
Figure 9.In fact, the time and angle-averaged gas quantities are also almost unchanged as we showed in Figure 10.This is likely because the timescales for the central black hole to tilt are so much longer than the dynamical times of the accretion flow at near horizon scales that the flow can adiabatically adjust as it tilts.
Observationally speaking, the jet precession we see in our binary simulations could generally be detected (or at least inferred) from radio observations of jet morphologies.The most direct signature of precession would be quasi-periodic variations in the radio position angles and fluxes of observed jets (e.g., Britzen et al. 2018;Britzen et al. 2023;Cui et al. 2023;von Fellenberg et al. 2023).The latter effect is caused by variations in the relativistic Doppler beaming as the jet points more or less toward the Earth.Other evidence for precession could be more subtle, such as the appearance of what looks like multiple jet components at different locations and directions (caused by the jet changing direction over time; Lister et al. 2013;Nandi et al. 2021) or variations in uncollimated outflow features (Falceta-Gonc ¸alves et al. 2010;Britzen et al. 2019;Krause et al. 2019).

Accretion and Magnetic Flux Accumulation on the Secondary
In this section, we describe the accretion properties onto the secondary black hole.To do this, we boost the simulation data into the rest frame of the secondary using the coordinate transformations described in §2.2 and then convert and interpolate the data onto a spherical grid in local Kerr-Schild coordinates.The transformation into spherical Kerr-Schild coordinates utilizes the single-black hole expressions (neglecting the effects of the primary on the metric), which is only appropriate for small distances from the secondary.We then calculate the accretion rate (relative to the time-averaged single black hole accretion rate) and dimensionless flux threading the secondary's horizon in the usual way (Equations 24 and 25).Figures 18 and 19 show these two quantities for our r BH,0 = 80r g and r BH,0 = 20r g simulations, respectively.For both i 0 = 90 • simulations, the accretion rate and black hole flux vary from 0 (when the secondary passes through the jet) to some peak value (when the secondary passes through the midplane of the disk) and then back to 0 every half orbit.The peak values of | Ṁs |/|⟨ Ṁp ⟩| are generally larger for r BH,0 = 80r g , i 0 = 90 • (∼ 0.3-0.6)than r BH,0 = 20r g , i 0 = 90 • (∼ 0.1-0.6).This is because the r BH,0 = 80r g black hole is moving slower than the r BH,0 = 20r g black hole and so it can accrete more gas even though it is surrounded by lower densities.We showed in §3 that these two effects compete to result in a scaling of Ṁs ∝ r 0.7 BH,0 , which is (80/20) 0.7 ≈ 2.6 for r BH,0 = 20r g and r BH,0 = 80r g , consistent with our findings.The peak values of | Ṁs |/|⟨ Ṁp ⟩| are significantly higher for i 0 = 45 • than i 0 = 90 • for both sets of simulations by factors of 2-3.This is because the black holes with i 0 = 45 • orbits are travelling with prograde motion through the accretion disk and so the net velocity of the gas in the frame of the secondary is reduced, leading to a larger accretion radius and accretion rate.The magnitude of this increase is consistent with our analysis in §3.2, where we predicted that the accretion rates for i 0 = 45 • would be ≈ 2.5 times the i 0 = 90 • accretion rates.
Note that since the secondary black hole is 10 times less massive than the primary, for all simulations the Eddington ratio for the secondary is larger than the primary (at times by as much as a factor ≳ 10).Accretion properties of the secondary black hole for r BH,0 = 80r g for i 0 = 90 • (solid blue ) and i 0 = 45 • (dashed green).Top: Accretion rate normalized to the time-averaged accretion rate onto the primary, Ṁs /|⟨ Ṁp ⟩|.Bottom: Dimensionless black hole flux, ϕ BH , calculated using the time-averaged accretion rate onto the secondary.The secondary accretes a significant amount of gas while passing through the disk (∼ 0.2 Ṁp for i 0 = 90 • and ∼ 0.6 Ṁp for i 0 = 45 • ).Relative flux accumulation is also significant (ϕ bH ∼ 6-10 for both inclinations), but does not attain to the MAD level.Note that the total magnetic flux, Φ BH , is larger in the i 0 = 45 • simulation than the i 0 = 90 • simulation by roughly the same factor as | Ṁs |, resulting in similar values for ϕ BH .
In terms of magnetic flux, all simulations show a significant amount of flux accumulation but not enough to reach the MAD state.All simulations show similar time-averaged values of ϕ BH ∼ 5-8, with quasi-periodic variation.There is no indication that the secondaries in any simulation will eventually accumulate enough magnetic flux to become MAD.This is partially because the accretion rate is not steady; as the secondary passes away from the midplane, the reduced accretion rate allows magnetic flux to be expelled.Top: Accretion rate normalized to the time-averaged accretion rate onto the primary, Ṁs /|⟨ Ṁp ⟩|.Bottom: Dimensionless black hole flux, ϕ BH , calculated using the time-averaged accretion rate onto the secondary.The secondary accretes less than the r BH,0 = 80r g simulations, but still a significant amount (∼ 0.1-0.2Ṁp for i 0 = 90 • and ∼ 0.3-0.5 Ṁp for i 0 = 45 • ).The accumulated flux is roughly comparable to the r BH,0 = 80r g simulations, ϕ bH ∼ 5-8.Note that the total magnetic flux, Φ BH , is larger in the i 0 = 45 • simulation than the i 0 = 90 • simulation by roughly the same factor as | Ṁs |, resulting in similar values for ϕ BH .

Accretion Disk Perturber Simulations
There have not been many studies of impacts of smaller objects with supermassive black hole accretion disks in GRMHD.The most relevant to the current work is Suková et al. (2021;hereafter S21), in which the authors simulated the passage of objects (including stars and black holes) of various sizes through a MAD accretion disk (see also Pasham et al. 2024 where similar simulations are used to make a case for the existence of a secondary black hole in the source ASASSN-20qc).This was done by enforcing that all the gas within a specified radial distance from the center of the objects has the same velocity as the objects, which are on geodesic orbits calculated alongside the simulations.The authors investigated a range of orbital distances (10-50 r g ) and influence radii (0.1-10 r g ).Note that the latter quantity does not necessarily correspond to the radius of the object itself but the radial range where the secondary has a significant effect on the accretion flow.Even for an influence radius of 1r g , in 2D S21 found that the secondary significantly altered the accretion flow, resulting in quasi-periodic oscillations of the accretion rate onto the primary, effectively shutting off accretion with every passage of the object through the midplane.This same motion produced quasi-periodic relativistic mass outflow rates with peaks that were 1-2 orders of magnitude higher than the "background" mass outflow rates.In the one 3D simulation4 with a secondary, the authors note that these effects are greatly diminished because the object now has a more realistic size in the azimuthal direction; in 2D the perturber was essentially a ring extending across the full 2π in azimuth.
The black holes in our simulations have influence radii self-consistently set by the dynamical interaction between gravity and the MHD fluid, but we have roughly estimated them as ∼ 2-3 r g , 4-5 r g , 10 r g , and 18-19 r g for (r BH,0 = 20r g , i 0 = 90 • ), (r BH,0 = 20r g , i 0 = 45 • ), (r BH,0 = 80r g , i 0 = 90 • ), and (r BH,0 = 80r g , i 0 = 45 • ), respectively.These are all larger than the fiducial 1r g used in S21, yet we see almost no effect on the resulting primary accretion rates and the magnetic flux threading the black hole, and only marginal changes in the time and angle-averaged radial profiles of gas quantities.We do however, see quasi-periodic outflows caused by the secondaries similar to those of S21, though the peak outflows are only 2-4 times the "background" outflow rates provided by the disk.This is consistent with their 3D result.Our peaks are also much more variable in magnitude and shape likely due to the increased turbulence and variability in the 3D disk compared with 2D.
The biggest differences between our simulations and S21 are 1) all of our simulations are fully 3D and 2) we focus specifically on black hole perturber and utilize a full treatment of the resulting binary metric.1) is particularly important for a number of reasons.First, there is no realistic way to treat a ballistic spherical object moving through an azimuthally symmetric accretion flow.The 2D perturber will act as ring with a substantially enhanced effect on the flow, as noted by S21.The 3rd dimension also allows us to study inclined orbits in a more straight-forward way.Furthermore, 2D, MRI-driven accretion flows are unrealistic when run for any significant length of time (more than a few thousand M) because the MRI is not sustainable in axisymmetry (Cowling 1933).2) allows us to self-consistently study the effects of the black holes on the accretion flow and to investigate the effects of spin-orbit coupling on the primary black hole.For example, the influence radii of the black holes are set by a combination of the orbital parameters, the accretion flow dynamics, and the secondary-to-primary mass ratio.The precession of the primary accretion disk caused by spin-orbit coupling may also be the most significant observable effect of the secondary for certain parameters.

Tilted Disk Simulations
The interaction of the accretion disk in our r BH,0 = 20r g simulations with the precessing primary black hole spin axis (due to spin-orbit coupling with the orbit of the secondary black hole) is in some ways similar to the interaction of an incoming accretion disk tilted with respect to a fixed black hole spin axis.For thick disks, such a situation has been studied by several authors in GRMHD (e.g., Fragile et al. 2007;McKinney, Tchekhovskoy & Blandford 2013;Liska et al. 2018;White, Quataert & Blaes 2019;Ressler, White & Quataert 2023;Chatterjee et al. 2023).These simulations include both SANE and MAD disks, which were found to have different alignment properties.In particular, the strong magnetic fields rotating with the black hole in MAD disks are very efficient at aligning accretion flows and jets (called magneto-spin alignment in McKinney, Tchekhovskoy & Blandford 2013), at least for a ≳ 0.9 and misalignments ≲ 60 • (larger misalignments tend to inhibit the MAD state, Ressler, White & Quataert 2023;Chatterjee et al. 2023).Alignment in MAD disks is seen in the simulations out to at least ≳ 100r g (Ressler, White & Quataert 2023;Chatterjee et al. 2023) and in reality could reach even larger distances on longer timescales.The inner accretion flow (r ≲ 10-20r g ) tends to align on timescales of ∼ 10 4 M for misalignments of ≲ 60 • (see Figure 14 in Ressler, White & Quataert 2023 and Figure 3 of Chatterjee et al. 2023).The accretion flow at larger and larger radii takes progessively longer times to align (e.g., ∼ 2 × 10 4 M at r = 50r g ).Jets in tilted MAD disks tend to align on even shorter timescales and out to several hundred r g (e.g., Figure 15 in Ressler, White & Quataert 2023).
Thick SANE disks, on the contrary, do not align efficiently but instead tend to form standing shocks as the gas accretes from the disk onto the black hole (Fragile et al. 2007;White, Quataert & Blaes 2019) and perhaps precess about the spin axis due to the Lense & Thirring (1918) effect in the azimuthal direction (defined with respect to the black hole spin axis, Liska et al. 2018), though it is argued in Chatterjee et al. (2023) that this precession is short lived and the true steady state of tilted SANE flows is instead a warped disk without precession.
Our simulations all contain MAD disks.The biggest difference between our study and the aforementioned works on tilted disks (apart from the presence of the secondary black hole) is that the misalignment between the angular momentum of the disk and the primary black hole spin axis is introduced gradually via the slowly changing spin axis instead of suddenly.Despite this, the properties of the disk alignment agree very well with previously studied tilted MAD simulations.There is, however, an observable difference between the single black hole case and the binary case.Because thick MAD disks align so well, precession is not seen in single black hole simulations.For binary simulations with spin-orbit coupling, however, precession would be observed even with strong alignment because the primary spin axis is changing with time.Moreover, these disks would also be persistently warped (i.e., the angular momentum vector changes with radius) because the outer part of the accretion flow aligns with the time-averaged primary spin axis while the inner part of the accretion flow aligns with the instantaneous spin axis.

LIMITATIONS OF OUR STUDY
The simulations we have presented are formally applicable only to low luminosity supermassive accretion flows where radiative cooling is inefficient and the disk is thick due to thermal pressure support.However, most observed AGN are in the radiative efficient regime where the disk is either thin from rapid cooling or thick from radiative pressure support.We expect that the effect of a secondary on a thin disk to be more dramatic than a thick disc because it will impact a larger fraction of its volume (see our analytic argument in §3.2).At very small disc thicknesses, the process of accretion and ejection will no longer be well approximated by our Bondi-Hoyle-Lyttleton framework described in §3 because a spatially extended wind will no longer be a good approximation in the frame of the secondary.Orbits of low inclination may also be able to have more significant unbound mass outflows because they will fully enter and exit the polar regions (unlike our i 0 = 45 • simulations).On the other hand, in the radiative pressure dominated, thick disk regime, it is reasonable to expect that many of our qualitative conclusions may still hold due to the similar geometry of the system compared to thick, non-radiative disks.The biggest difference would likely be the significantly lower gas temperatures.These are particularly important for determining the emission associated with the ejection of the gas into the polar region, which is determined by a combination of geometry and the photospheric temperature of the "blobs" (Franchini et al. 2023).For the highest accretion rates, self-gravity of the accreting gas and dynamical friction on the secondary may also become important, which would affect the orbit of the secondary.
In this work, we have only considered MAD accretion flows.While recent work has shown that MAD disks may be relatively common for supermassive black holes (e.g., Event Horizon Telescope Collaboration et al. 2019;Akiyama et al. 2022;Ressler et al. 2020;Liska, Tchekhovskoy & Quataert 2020), there are still likely many AGN either in less magnetized states or with more toroidally dominated magnetic fields given that most do not display obvious jet signatures (whereas MAD disks generally have strong jets).The powerful magnetic fields in MAD disks are known to be more efficient at aligning the accretion disk to tilted black hole spins than the weaker magnetic fields in less magnetized disks (McKinney, Tchekhovskoy & Blandford 2013;White, Quataert & Blaes 2019;Ressler, White & Quataert 2023;Chatterjee et al. 2023).Thus the disk and jet in a SANE flow may not as closely mirror the black hole spin axis as our simulations.The outflow and jet from MAD disks are also much stronger than SANE flows by factors of 10-100 (Tchekhovskoy, Narayan & McKinney 2011;Ressler, White & Quataert 2023), which could alter the behavior of the blobs ejected by the secondary black hole.
We have also only considered rapidly rotating primary black holes.In less rapidly rotating black hole systems, the accretion flow would be less affected by spin-orbit coupling because the alignment would be less efficient.The jets would also be weaker than those in our a = 0.9375 simulations (jet power is a strong function of a, ∝ a 2 -a 4 ; Blandford & Znajek 1977;McKinney 2005;Tanabe & Nagataki 2008;Tchekhovskoy, Narayan & McKinney 2010) and perhaps be less efficient at blowing away gas brought into the polar regions by the secondary.
Similarly, we have focused only on secondary black holes with no spin.More rapidly rotating secondaries would likely produce their own jets which could have a stronger effect on the primary accretion flow.The strength of these jets would depend on the amount of magnetic flux accreted by the secondary; for secondaries accreting flux at the MAD level, their jets would inject ∼ | Ṁs c 2 | worth of power into the primary accretion flow, which can be a sizable fraction of ∼ | Ṁp c 2 | (i.e., comparable to the primary jet power, see Figures 18 and 19).The effect that such a deposition of energy would have on the primary flow is unclear, but it could be dramatic.Depending on the spin orientation with respect to the orbit and the primary black hole spin, the secondary spin would change direction due to spin orbit coupling in a similar way to the primary black hole spin (though on a much shorter timescale due to the lower relative mass).The jets from the primary and secondary could also in principle interact with each other (Molnar et al. 2017;Volonteri et al. 2022;Gutiérrez et al. 2023; though this may only be significant for mass ratios closer to one where the jets are of comparable strength and size).In addition, a rapidly spinning secondary could torque the surrounding accretion flow and alter the dynamics in that region more prominently than a non-spinning secondary.
Here we have presented only one choice of mass ratio, q = 0.1.For mass ratios higher than this our assumption of a stationary primary black hole and steady state accretion disk would start to break down and the system would no longer be in the perturber regime.We have performed simulations with significantly smaller mass ratios and found the effect of the secondary on outflows and accretion to be almost neglible and thus not particularly interesting to present.Smaller mass ratios also make it more computationally demanding to resolve the secondary's influence region.
Finally, we have only considered primary accretion flows that are aligned with the primary black hole spin axis.In reality, supermassive black hole accretion disks could be significantly tilted, which would reduce the amount of material in the polar regions, reduce the strength of the jet, and change the geometry of the accretion flow (Fragile et al. 2007;White, Quataert & Blaes 2019;Chatterjee et al. 2020;Ressler, White & Quataert 2023;Chatterjee et al. 2023).The consequences of these combined effects could have on our results are not clear, though if the jet is weaker and more entrained with matter the outflows produced by the secondary black hole may be relatively weaker.
Future work should explore a larger amount of parameter space by investigating disks of different thicknesses, tilts, magnetic field strengths (i.e., SANE vs. MAD flows), and black hole spins.

DISCUSSION AND CONCLUSIONS
We have presented the results of 3D simulations of small mass ratio (q = 0.1) binary black holes in the inspiral regime at the centers of galaxies, considering four different scenarios: secondaries with orbits at initial separations of 20r g and 80r g with initial inclinations of 90 • and 45 • from the midplane.We did this by implementing the time-dependent approximate metric from Combi et al. (2021) describing a superimposition of boosted Kerr black holes into Athena++.This metric, though approximate, fully captures general relativistic effects such as spin-orbit coupling on both the primary and secondary black holes, the inspiral of the secondary, the event horizon contractions when the black holes are moving close to the speed of light, and the Blandford & Znajek (1977) process that can drive electromagnetically dominated jets.The secondaries are introduced into an accretion flow around a supermassive black hole that has already evolved for 10 5 M (Figure 1), long enough to come into inflow equilibrium out to ∼ 100r g (Figure 2).
We find that even for a relatively high secondary mass of 1/10 the mass of the primary, the overall effects on the pri-mary accretion flow are small as measured by the accretion rate, dimensionless magnetic flux threading the black hole (Figure 9), and the angle averaged radial quantities (Figure 10).This is because the area on a spherical shell swept up in the disk by the secondary during an orbit at the orbital radius is much smaller than the overall area of the disk at that radius.In other words, the fraction of the disk affected by the secondary as determined by the orbit and influence radius is very small (Figure 7).
The secondary black holes in our simulations do, however, produce quasi-periodic outbursts with periods ≳ the orbital period as seen in the unbound mass outflow rates (Figures 12).These outbursts are caused by the secondary bringing gas from the disk into the polar regions, where the electromagnetically dominated jet then blows this gas away (Figure 11).When measured at distances from the primary of twice the secondary orbital radius, the peaks of the mass unbound outflow rates are ∼ 2-4 times the "background" unbound outflow rate provided by the accretion disk.The orbital frequency is clearly seen in the power spectra of the unbound outflow rates at these distances (Figure 13), especially for secondaries with i 0 = 90 • orbits.However, there are also several other prominent frequencies, particularly lower frequencies and especially for r BH,0 = 20r g .At distances farther away from the primary (e.g., 5r BH,0 ), peaks in unbound outflow rates become less frequent and more spread out in time (Figure 14), with power shifting to frequencies lower than half the orbital frequency.At even larger radii (e.g., ≳ 10r BH,0 ), the variability in the unbound outflow rates shifts to even lower frequencies and instead of having welldefined peaks is characterized by a more continuous distribution in time consistently ∼ 20-30% larger than the single black hole unbound outflow rates.This means that quasiperiodic features are only present in the unbound outflow rates for r BH,0 ≲ r ≲ 10r BH,0 .We interpret this as a result of the mass ejected from the primary accretion disk by the secondary not propagating/accelerating as a coherent structure as it is propelled to larger radii by the jet.Instead, it spreads out in space and can even catch up and merge with previously ejected matter.Eventually, by about r ≈ 10r BH,0 the initially discrete chunks of outflow have diffused into a more continuous wind.
Because the accretion disk used in this work is so thick, orbits with i 0 = 45 • do not fully escape the disk into the polar regions and so their ability to deposit matter into the jet is reduced.At the same time, the secondaries on these orbits have larger influence radii because the relative velocity of the gas in the frame of the secondary is reduced.These combined effects results in in similar outburst magnitudes as orbits with i 0 = 90 • (Figure 12).Our simple analytic model predicts that for thinner disks in which the secondary fully passes into the polar regions, secondaries with i 0 = 45 • or-bital inclinations will have larger outbursts than secondaries with i 0 = 90 • orbital inclinations ( §3.2).
For smaller r BH,0 (e.g., r BH,0 = 20r g ), spin-orbit coupling between the primary black hole spin and the orbit of the secondary can cause the spin direction of the primary to change by up to 60 • on timescales we can simulate (∼ 10 4 M 1 , Figure 6).This results in both the disk and jet of the primary tilting along with the spin axis (Figures 15,16,and 17) in an approximately adiabatic way.As the spin axis returns to its original orientation, the jet and inner region of the accretion disk (r ≈ 5r g ) also return to their original orientation.The larger radii accretion disk (r ≳ 50r g ), however, remains tilted by ∼ half of the peak spin tilt, even as the primary spin axis continues to change.This is because the dynamical timescales in the disk at these radii are longer than the timescale for the spin to change, resulting in the larger radii flow seeing an effective time-average of the primary spin.Previous work in tilted accretion flows around single black holes has found that thick disks tend to either align or warp but not consistently precess (see §5.2).Thus observations of precession in thick AGN disks would be strong evidence for the presence of a secondary black hole companion.
The secondary black holes themselves accrete a significant amount of mass from the primary accretion disk, from 0.1 times the primary accretion rate up to greater than the primary accretion rate (Figures 18 and 19), meaning that the Eddington ratio for the secondaries is larger than the Eddington ratio for the primary.Accretion rates vary in a quasiperiodic manner, with peak accretion occuring as the secondary passes through the midplane of the primary accretion disk.Accretion rates generally increase with increasing r BH,0 (roughly as r 0.7 BH,0 ), though there is significant overall timevariability in the magnitude of the individual peaks.Accretion rates onto secondaries with i 0 = 45 • orbits are 2-3 times higher than those onto secondaries with i 0 = 90 • orbits because the relative velocity of the gas in the accretion disk in the frame of the secondary is reduced.The black hole flux threading the secondaries is, on average, ∼ 15% of the MAD limit for all orbital configurations.We see no indication that the secondaries will ever become MAD because of the frequent dips in accretion rate during which the magnetic flux can leak out of the black hole.
Directly predicting the electromagnetic emission from our simulations is difficult without a detailed treatment of radiation, which we reserve for future work. in particular, there are several features of our simulations could imprint themselves on observations.For instance, the quasi-periodic ejection of matter from the disk to the jet could directly result in quasi-periodic emission if the hot gas radiates significantly between 1-10 r BH,0 .This could happen via thermal emission from the gas soon after it is expelled, or by interactions between jet/outflow material and the disk or a surrounding medium.Alternatively, emission could be generated via the shock front of the secondary and the associated accretion onto the secondary.Note, however, that red noise contamination can also result in apparent periodic features in AGN light curves and thus make distinguishing true periodic signals more challenging.The slow precession of the disk and jet from spin-orbit coupling between the orbit of the secondary and the primary black hole spin axis could be detectable through direct jet imaging or by variations in Doppler effects.Finally, the matter brought into the polar regions by the secondary induces an interaction between the magnetic fields in the highly magnetized jet and in the mildly magnetized matter from the secondary.The turbulence resulting from this interaction may cause dissipation of strong magnetic field, e.g., through magnetic reconnection, that could potentially power bright and fast flares.More accurately studying this possibility would require much higher resolution simulations.
Our results are valid for thick accretion disks associated with low luminosity AGN.They may also be qualitatively applicable to thick accretion disks associated with ∼ Eddington and super-Eddington accretion where the gas becomes radiation pressure dominated due to the similar geometry.Future studies of the effect of secondary black holes on thin accretion disks will require a proper treatment of radiative cooling.Future work should also explore a larger range of parameter space in this binary black hole perturber scenario by investigating SANE disks, different black hole spins, tilted disks, and different mass ratios.
This work has implications for observed quasi-periodic outflows/outbursts/eruptions and the detection of companion supermassive or intermediate mass black holes.We have shown that the signatures of the latter may be subtle even when the mass ratio is relatively high (q = 0.1) if the disk is thick.
We thank the referee for their careful reading of the manuscript and useful comments.We acknowledge the support of the Natural Sciences and  The computational resources and services used in this work were partially provided by facilities supported by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government -department EWI.This research was supported in part by the NSF through XSEDE computational time allocation TG-AST200005 on Stampede2 and Bridges-2.This work was made possible by computing time granted by UCB on the Savio cluster.Figure 20 shows the resulting density profile at four different times for N = 2048.The initial overdensity feature moves to the left and then to the right (asymmetrically) before returning to the center, essentially overlapping with the initial condition.
Plotted in Figure 21 is the L1 norm of the density error, vs. resolution for three different frequencies of updating the metric.The most accurate updates the metric every subtimestep of the algorithm (i.e., every half time step), the next most accurate updates the metric every full timestep, and the least accurate updates the metric ever 10 full timesteps.Figure 21 demonstrates that all three of these methods converge, though at different orders.As expected, updating every substep of the algorithm results in asymptotically 2nd order convergence with resolution, reproducing the convergence properties of a static metric flow in Athena++ (White, Stone & Gammie 2016).Updating the metric only every full timestep formally reduces the order of the algorithm to 1st order, however, we find that the for this test problem the convergence is noticeably less than 2nd order but beats 1st order at least up until N = 2048, the highest resolution we simulate.This is because spatial errors in the gradient of the gas can still dominate the temporal errors at low resolution.Finally, updating  when the metric is updated ever half time-step (green circles), every full time step (magenta triangles), and every 10 timesteps (blue squares).Updating the metric every half timestep results in asymptotically 2nd order convergence while updating every 10 timesteps asymptotes to 1st order convergence.The curve for updating ever timestep is very close to the updating every half time step curve until relatively high resolution (N ≳ 512) at which point the convergence flattens out to slower than 2nd order but faster than 1st order.At higher resolution we expect this curve to continue to flatten until it reaches ∝ r −1 .The L1 norm of the error in pressure and magnetic field look similar.
the metric only once every 10 timesteps results in asymptotically 1st order convergence.
The absolute values of the errors in Figure 21 are specific to this simple test problem and should not be generalized.Instead, our goal of this test is to demonstrate that the metric update described in §2 converges properly (at 2nd or 1st order depending on whether the metric is updated every subtimestep or just every full timestep or less).In particular, we demonstrate that even updating the metric every 10 timesteps (the method employed in the main body of this work) results in 1st order convergence.

A.2. Boosted Bondi Accretion
Another test of our implementation is a single black hole accreting spherically in its rest frame moving through a grid at some velocity v BH .To do this, we take the spherically symmetric accretion solution for a nonspinning Schwarzschild black hole (Hawley, Smarr & Wilson 1984, see also Section 4.4 of White, Stone & Gammie 2016) and then boost into the frame described by the trajectory s x = s y = 0, s z = v BH t + z 0 .This is a special case of the metric described in §2.2, where M 2 = 0 and the transformation is exact since v BH is constant in time.The resulting solution is length-contracted in the direction of motion (z) and thus no longer spherically sym- Left: initial conditions.Right: final solution at t = 100M.For plotting purposes the z-axis is shifted so that the black hole is located at the origin.Yellow lines demarcate the regions that are "active" in the simulation.Regions outside the outermost ellipse and inside the innermost ellipse are fixed to the analytic solution.The inner black surface covers the event horizon of the black hole.Note that the boosted Bondi solution is no longer spherical, and regions of constant r ′ (radial distance from the black hole in its rest frame) are ellipses in the lab frame.The differences between the final and initial states are indistinguishable by eye.metric in the lab frame but only trivially time-dependent in that it moves uniformly with a constant velocity.
We simulate this by initializing a −100r g < x, y, z < 100r g box with the boosted Bondi solution for z 0 = 80r g and v BH = 0.9.We use 5 levels of nested adaptive mesh refinement centered on the moving black hole.The cells within r ′ ≤ 3 and r ′ ≥ 10, where r ′ is the radial distance from the black hole in its rest frame, are fixed to the analytic solution and the simulation is run for 100 M. The resulting mass density is shown in Figure 22 plotted alongside the initial condition, showing that the simulation accurately reproduces the analytic solution.
Quantitatively, we show the convergence properties of the L1 norm of the density error at t = 100 M in Figure 23, comparing a simulation that updates the metric every half timestep to one that updates the metric every 10 timesteps.Both simulations do very well at maintaining the Bondi solution, even at relatively low resolution (e.g., L1 norms of density errors ≲ 3 × 10 −3 even at a base resolution of 32 3 ).The errors in both simulations converge to zero and are within a factor of ∼ 2 of each other.The simulation that updates the metric every 10 timesteps converges roughly at the expected rate of N −1 .The simulation that updates the metric every half timestep, however, converges at a slower rate than ex-  for a simulation that updates the metric every half timestep (blue triangles) compared to a simulation that updates the metric ever 10 timesteps (green circles).Note that the resolution N on the x-axis is the base resolution in one direction (which is the same in every direction).On top of this, there are 5 additional levels of adaptive mesh refinement for each simulation.At the highest base resolution tested (N = 256), the simulation that updates the metric every half timestep converges at a rate steeper than N −1 but shallower than N −2 .The lower-than-expected order of convergence is likely due to the nature of the moving boundary condition.The simulation that updates the metric every 10 timestep does quite well in comparison, converging at approximately r −1 and only a factor of ∼ 2 higher error than the simulation that updates the metric every half timestep.Both errors are quite low at all resolution (≲ 3 × 10 −4 ).
pected, just moderately steeper than N −1 but shallower than N −2 .There are at least three possible causes of this. 1) The inner/outer boundaries at r ′ = 3r g and r ′ = 10r g , which move along with the black hole so that cells are continuously entering and exiting the computational domain.2) This same boundary is ellipsoidal on a Cartesian grid, so it is misaligned with the coordinate faces and is located at slightly different locations at different resolution.3) AMR, which can focus resolution at slightly different locations at different base resolution.The errors introduced by these three factors, of course, converge to zero with increasing resolution, but they converge at a rate ∝ N −1 (at least for the latter two, since the errors are independent of time and ∝ ∆x) .If these errors are dominant that would explain why the simulations that update every 10 timesteps and those that update every half timestep converge at a similar rate.Given the relative complexity of this test problem compared to a stationary, face-aligned boundary and uniform grid test problem, however, we do not consider it worrisome that the convergence rate does not reach 2nd order, especially since 1) the errors in both simulations are small and still converge to zero and 2) we update the metric every 10 timesteps in our target problem 3 for a black hole/gas moving at 0.9c and base resolution of 128 3 (including 5 levels of adaptive mesh refinement for the stationary black hole simulations and 7 levels fo adaptive mesh refinement for the moving black hole simulation).Left: simulation for a moving black hole through an initially stationary gas with the metric updated every 10 timesteps at t = 300M.Right: simulation for a stationary black hole with gas moving at an initially constant velocity at t ′ = 300M.For plotting purposes the z-axis is shifted so that the black hole is located at the origin.The stationary black hole simulation is also boosted into the frame of the moving black hole for ease of comparison.The inner black surface covers the event horizon of the black hole (which is ellipsoidal in this frame).At these late times the flow has settled into a steady state solution of an extended shock front trailing the black hole.The differences between the solutions in the two simulations are almost indistinguishable by eye.
(see §2) and thus would not expect to see N −2 convergence anyway.

A.3. Boosted Bondi-Hoyle-Lyttleton Accretion
For this test, we compare a simulation of a stationary black hole embedded in uniform gas moving with constant velocity to a simulation of a moving black hole through a stationary, uniform gas.The latter is simply a boosted version of the former, and so we can use the metric of §2.2 with a constant velocity and M 2 = 0.This problem, in general, doesn't have a known analytic solution, and is in fact still being actively researched (e.g., Blakely & Nikiforakis 2015;Gracia-Linares & Guzmán 2023;Kaaz et al. 2023).It is, however, a good test to see whether the time-dependency of the metric is accurately captured by our simulations in a complex problem.
For the stationary black hole simulation, we use a grid of size (100r g ) 3 with resolution of 128 3 and 5 additional levels of static mesh refinement centered on the black hole.The gas is initialized with ρ = 1, P = 5, u x = u y = 0, and u z /u t = −v BH everywhere except for regions with r ′ ≤ 5r g , where ρ and P are set to the floors and the three-velocities are set to zero.These parameters describe a mildly supersonic flow, with Mach number M ≈ 1.15 (Mignone & McKinney 2007).
For the moving black hole simulation, we use a grid of size (400r g ) 3 with resolution of 128 3 and 7 additional levels of adaptive mesh refinement centered on the black hole (achieving the same effective resolution as the stationary black hole case).The gas is initialized with ρ = 1, P = 5 and threevelocities set to zero everywhere except for regions with r ≤ 5r g , where ρ and P are set to the floors.The black hole is initially located at x = 0, y = 0, and z = −380r g moving at v BH = 0.9 in the +z direction.The metric is updated every 10 timesteps.
Both simulations are run for 300 M in their own reference frame, at which point the solutions reach a steady state.2D slices of the density in these solutions are shown in Figure 24.When boosted for comparison, these two solutions agree very well, both showing a wide shock front extending behind the black hole as it moves through the uniform medium (or as the uniform wind impacts the black hole in the stationary black hole reference frame).The agreement between the solutions is strong evidence that our treatment of the time-dependent, boosted metric is accurate for complex problems even when the metric is only updated every 10 timesteps.
B. NONLINEAR COORDINATE EQUATION Equation ( 16) in §2.2 results in the following nonlinear set of equations: and

Figure 2 .
Figure2.Properties of the single black hole accretion flow used as initial conditions for our binary simulations.Top: time-averaged accretion rate, ⟨ Ṁ⟩ vs. radius (averaged over 50,000-100,000 M 1 ).Bottom: dimensionless black hole flux threading the black hole, ϕ BH , vs. time.⟨ Ṁ⟩ is roughly constant out to ≈ 70-100 r g , indicating that the flow has reached inflow equilibrium out to this distance.The accretion flow is firmly in the MAD state, with saturated magnetic flux around ϕ BH ≈ 60(Tchekhovskoy, Narayan & McKinney 2011) and going through cycles of build-up followed by dissipation every ∼ 3,000 M.

Figure 3 .
Figure 3. Representative example of our static and adaptive mesh-refined grid.A mass density contour is plotted in the y-z plane for our r BH = 80 • , i 0 = 90 • simulation, with yellow lines demarcating blocks of 16 3 cells.The red and blue outlined subplots zoom in on the primary and secondary black holes, respectively.The mesh refinement effectively allows us to resolve multiple scales, particularly surrounding the two black holes.

Figure 4 .
Figure 4. Binary separation distance vs. time for the 4 different orbits used in this work.Top: r BH,0 = 20r g for i 0 = 90 • (red, solid)and i 0 = 45 • (black, dashed).Bottom: r BH,0 = 80r g for i 0 = 90 • (red, solid) and i 0 = 45 • (black, dashed), where the inset zooms in on the time covered by our simulations.For orbits with initial separations of r BH,0 = 20r g , the timescale for the orbital radius to significantly change is ∼ 10 4 M 1 , much shorter than the ∼ 10 6 M 1 timescale for the r BH,0 = 80r g orbits.Since 10 6 M 1 is much longer than the length of our simulations, for r BH,0 = 80r g we safely assume Keplerian circular orbits that do not change in time.For r BH,0 = 20r g simulations, however, we include the full orbital evolution and can expect to see significant changes on the timescales we can simulate.The differences between i 0 = 90 • and i 0 = 45 • binary separation distances are relatively small until close to merger (∼ 4 × 10 4 for r BH,0 = 20r g and ∼ 9 × 10 6 for r BH,0 = 80r g ).

Figure 7 .
Figure7.2D slice of density for our q = 0.1, r BH,0 = 80r g , and i 0 = 90 • binary black hole perturber simulation at four different times.The yellow circles represent the influence radii predicted from Equation (37).Starting from the upper left panel, and proceeding clockwise, the panels represent 980M 1 , 2080M 1 , 3250M 1 , and 4260M 1 since the secondary was introduced.As it passes through the accretion flow, the secondary creates a Bondi-Hoyle-Lyttleton-esque shock front (upper left and lower right panels; see also Appendix A.3) caused by its supersonic trajectory.Each time black hole passes through the corona of the disk into the jet, it carries with it some a small amount of matter that is then quickly swept away by the jet (e.g.upper right and lower left panel).For an animated version of this figure, see https://youtu.be/lU82eT18HLU.

Figure 8 .
Figure 8. 2D slices of the square sound speed divided by the square Keplerian velocity, c 2 s /v 2 kep , where c 2 s = γP/ρ and v 2 kep = 1/ √r in these units.The contours are centered on the secondary black holes as they pass through the midplane of the accretion disk during representative times in our simulations.The white circles represent the influence radii predicted from Equation (37).The moving black hole creates a bow shock in the gas of varying shapes and sizes.All simulations generally have moderate shocks with temperatures (note that c 2 s is directly proportional to temperature) increasing by factors of ∼ a few from the average gas temperature at the orbital radius.The i 0 = 45 • (right) simulations generally display a larger secondary influence radii but slightly weaker shocks than the i 0 = 90 • simulations (left) due to the reduced relative gas motion.Similarly, the r BH,0 = 80r g simulations (top) display a larger secondary influence radius than the r BH,0 = 20r g simulations (bottom).All simulations show good agreement with the predicted influence radii from Equation (37).
Figure9.Accretion rate, Ṁ, and dimensionless magnetic flux threading the primary black hole, ϕ BH , as a function of time for the single black hole simulation (black dotted line) compared with the q = 0.1, r BH,0 = 80r g , i 0 = 90 • (solid blue line) and q = 0.1, r BH,0 = 20r g , i 0 = 90 • (dashed red line) binary simulations.Neither the accretion rate nor the magnetic flux show any clear signature of the secondary orbit's periodicity.In fact, both quantities are very similar in the binary and single black hole simulations, with mainly stochastic differences.This shows that the effect of the secondary is small on the near-horizon dynamics of the flow.We argue that this is because the area swept out in the disk by the orbit of the secondary on a spherical shell is small compared with the area of the disk itself at the orbital radius.Note that Ṁ and ϕ BH for the i 0 = 45 • simulations are similar and are not plotted to avoid clutter.
Figure10.Time and angle-averaged radial profiles of density multiplied by r 0.8 (top set of lines), ⟨ρ⟩r 0.8 , and mass-weighted square sound speed divided by the square Keplerian speed (bottom set of lines), ⟨c 2 s /v 2 kep ⟩ ρ for our four different binary simulations compared with the single black hole simulation (black solid line).The secondary black hole has only a small effect on the average density and temperature, providing a small increase in temperature around the orbital radius (caused by the bow shock following the secondary) and a small decrease in density near and at smaller radii than orbital radius (caused by the expulsion of mass from the disk).

Figure 11 .
Figure 11.Evolution of an ejected blob in our r BH,0 = 80r g , i 0 = 90 • (top) and r BH,0 = 20r g , i 0 = 90 • (bottom) simulations.Plotted is a 2D slice of mass density in the x = x BH plane (parallel to y-z) at 8 different times for each simulation (time proceeds from top left to top right and then bottom left to bottom right).The white circles represent the influence radii predicted from Equation (37).Note that the two sets of contours are on different spatial scales.The secondary black hole pulls gas from the accretion disk into the polar region and then the electromagnetically dominated jet blows this gas away.The sizes of these blobs are consistent with the Bondi-Hoyle-Lyttleton radius predicted by Equation (37), ∼ 10r g in radius for r BH,0 = 80 r g and ∼ 2.5r g in radius for r BH,0 = 20r g .
Figure12.Unbound mass outflow rate, Ṁunbound , plotted vs. time elapsed since the secondary was introduced, t − t 0 , in our r BH,0 = 80r g (top) and r BH,0 = 20r g (bottom) simulations for i 0 = 90 • (blue solid line) and i 0 = 45 • (green dashed line), compared to the single black hole simulation (black dotted line).Ṁunbound is measured at r = 160r g for r BH,0 = 80r g and r = 40r g for r BH,0 = 20r g .The curves have been normalized by the time-averaged accretion rate over the interval (0 ≤ t − t 0 ≤ 40,000M 1 ).The secondary increases the unbound mass outflow rate at quasi-period intervals of order the orbital period (≈ 4500 M 1 for r BH,0 = 80r g and ≈ 560 M 1 for r BH,0 = 20r g ) by factors of ∼ 2-4 compared to the same quantity in the single black hole simulation.Note that as in all GRMHD torus simulations there is a secular decline in all outflow rates due to the mass in the torus depleting over time.
Figure13.Power spectra of unbound mass outflow rates, Ṁunbound , for orbits with i 0 = 90 • (blue solid) and i 0 = 45 • (green dashed) compared to the single black hole simulation (black dotted).Each spectrum is normalized such that the highest peak is 1.Top: r BH,0 = 80r g , with Ṁunbound measured at r = 160r g and the x-axis measuring frequency in units of the orbital frequency at r = 80r g .Bottom: r BH,0 = 20r g , with Ṁunbound measured at r = 40r g and the x-axis measuring frequency in units of the orbital frequency at r = 20r g .For r BH,0 = 80r g , the i 0 = 90 • PSD shows clear peaks at the orbital frequency and twice the orbital frequency, with other smaller peaks.The i = 45 • does have a peak at the orbital frequency but it is subdominant when compared with the lower frequency peak also seen in the single black hole simulation.For r BH,0 = 20r g , there are significant peaks at the orbital frequency for both i 0 = 90 • and i 0 = 45 • , but these are lower than some peaks at lower frequencies.This is likely due to propagation effects that spread out mass ejected by the secondaries in radius (and correspondingly, time in Ṁunbound ).
have become less distinct and spread out of over time (with the Ṁunbound curve now clearly rising above the single black hole Ṁunbound curve at almost

Figure 15 .
Figure 15.Volume renderings of our r BH,0 = 20r g , i 0 = 90 • simulation at 5 different times on a (100r g ) 3 (left) and (400r g ) 3 (right) scale, with the primary black hole spin axis direction indicated by a yellow arrow.In this figure, x is the horizontal direction (positive to the left) and z is the vertical direction (positive up), while y is into and out of the page (positive out of the page).The volume rendering highlights the highest σ = b 2 /ρ (green/blue) and density (red), showing how the disk and jet tilt and precess as a function of time due to binary spin-orbit coupling.Both disk and jet align tend to align with the changing black hole spin axis, particularly at smaller scales.For an animated version of this figure, see https://youtu.be/GdgAINSolOY.
Figure18.Accretion properties of the secondary black hole for r BH,0 = 80r g for i 0 = 90 • (solid blue ) and i 0 = 45 • (dashed green).Top: Accretion rate normalized to the time-averaged accretion rate onto the primary, Ṁs /|⟨ Ṁp ⟩|.Bottom: Dimensionless black hole flux, ϕ BH , calculated using the time-averaged accretion rate onto the secondary.The secondary accretes a significant amount of gas while passing through the disk (∼ 0.2 Ṁp for i 0 = 90 • and ∼ 0.6 Ṁp for i 0 = 45 • ).Relative flux accumulation is also significant (ϕ bH ∼ 6-10 for both inclinations), but does not attain to the MAD level.Note that the total magnetic flux, Φ BH , is larger in the i 0 = 45 • simulation than the i 0 = 90 • simulation by roughly the same factor as | Ṁs |, resulting in similar values for ϕ BH .

Figure 21 .
Figure21.L1 norm of the error in density, ρ, vs. resolution, Nm in our 1D accelerating reference frame test described in Appendix A.1 when the metric is updated ever half time-step (green circles), every full time step (magenta triangles), and every 10 timesteps (blue squares).Updating the metric every half timestep results in asymptotically 2nd order convergence while updating every 10 timesteps asymptotes to 1st order convergence.The curve for updating ever timestep is very close to the updating every half time step curve until relatively high resolution (N ≳ 512) at which point the convergence flattens out to slower than 2nd order but faster than 1st order.At higher resolution we expect this curve to continue to flatten until it reaches ∝ r −1 .The L1 norm of the error in pressure and magnetic field look similar.

Figure 22 .
Figure22.2D slice of density of the 3D boosted Bondi test described in Appendix A.2 for a black hole moving at 0.9c and base resolution of 128 3 (including 5 levels of adaptive mesh refinement).Left: initial conditions.Right: final solution at t = 100M.For plotting purposes the z-axis is shifted so that the black hole is located at the origin.Yellow lines demarcate the regions that are "active" in the simulation.Regions outside the outermost ellipse and inside the innermost ellipse are fixed to the analytic solution.The inner black surface covers the event horizon of the black hole.Note that the boosted Bondi solution is no longer spherical, and regions of constant r ′ (radial distance from the black hole in its rest frame) are ellipses in the lab frame.The differences between the final and initial states are indistinguishable by eye.

Figure 23 .
Figure23.L1 norm of density error in the boosted Bondi accretion test described in Appendix A.2 vs. 3D resolution at t = 100M for a simulation that updates the metric every half timestep (blue triangles) compared to a simulation that updates the metric ever 10 timesteps (green circles).Note that the resolution N on the x-axis is the base resolution in one direction (which is the same in every direction).On top of this, there are 5 additional levels of adaptive mesh refinement for each simulation.At the highest base resolution tested (N = 256), the simulation that updates the metric every half timestep converges at a rate steeper than N −1 but shallower than N −2 .The lower-than-expected order of convergence is likely due to the nature of the moving boundary condition.The simulation that updates the metric every 10 timestep does quite well in comparison, converging at approximately r −1 and only a factor of ∼ 2 higher error than the simulation that updates the metric every half timestep.Both errors are quite low at all resolution (≲ 3 × 10 −4 ).

Figure 24 .
Figure24.2D slice of density in the 3D boosted Bondi-Hoyle-Lyttleton test described in Appendix A.3 for a black hole/gas moving at 0.9c and base resolution of 128 3 (including 5 levels of adaptive mesh refinement for the stationary black hole simulations and 7 levels fo adaptive mesh refinement for the moving black hole simulation).Left: simulation for a moving black hole through an initially stationary gas with the metric updated every 10 timesteps at t = 300M.Right: simulation for a stationary black hole with gas moving at an initially constant velocity at t ′ = 300M.For plotting purposes the z-axis is shifted so that the black hole is located at the origin.The stationary black hole simulation is also boosted into the frame of the moving black hole for ease of comparison.The inner black surface covers the event horizon of the black hole (which is ellipsoidal in this frame).At these late times the flow has settled into a steady state solution of an extended shock front trailing the black hole.The differences between the solutions in the two simulations are almost indistinguishable by eye.
[0.85, 0.93, 0.85, 0.87, 0.86], and standard deviation of φ BH The five simulations have averageṀ = [2.2,2.1, 2.1, 2.1, 2.3], average φ BH =[56, 59, 56, 56, 58], standard deviation of Ṁ = Engineering Research Council of Canada (NSERC), [funding reference number 568580] Cette recherche a été financée par le Conseil de recherches en sciences naturelles et en génie du Canada (CRSNG), [numéro de référence 568580].B.R. is supported by the Natural Sciences & Engineering Research Council of Canada (NSERC) and by a grant from the Simons Foundation (MP-SCMPS-00001470).Research at the Flatiron Institute is supported by the Simons Foundation.L.C. and H.Y. are supported in part by Perimeter Institute for Theoretical Physics.Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Colleges and Universities.