<figure><img src="image_1.png" /></figure>

HAL
open science

# The Structural and Orbital Effects of Active Galactic Nuclei Feedback on SMBH Binaries Embedded in Gaseous Circumbinary Disks

Juan D'etigny, Andrés Escala, Joakim Rosdhal

► To cite this version:

Juan D'etigny, Andrés Escala, Joakim Rosdhal. The Structural and Orbital Effects of Active Galactic Nuclei Feedback on SMBH Binaries Embedded in Gaseous Circumbinary Disks. The Astrophysical Journal, 2024, 971 (1), pp.38. 10.3847/1538-4357/ad5352. hal-04679353

HAL Id: hal-04679353

https://hal.science/hal-04679353v1

Submitted on 30 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# The Structural and Orbital Effects of Active Galactic Nuclei Feedback on SMBH Binaries Embedded in Gaseous Circumbinary Disks

Juan d'Etigny¹, Andrés Escala¹, and Joakim Rosdhal²

¹Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile

²Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France
Received 2023 March 20; revised 2024 May 14; accepted 2024 May 30; published 2024 August 5

## Abstract

Using subparsec-scale-resolution radiation+hydrodynamical adaptive mesh refinement simulations deployed with the RAMSES code, we study the dynamics of supermassive black hole (SMBH) binaries embedded in gaseous nuclear circumbinary disks, where we investigate the effects of active galactic nucleus feedback on the SMBH binaries' migration behavior and disk structure. The radiative feedback effects are modeled by injecting photons that interact with the gas, through the adoption of a grid of BH emission spectra. We run simulations with initial conditions that lead by pure gravity plus hydrodynamics both to the formation of a low-density tidal cavity and to systems where gas-viscous diffusion is efficient enough to maintain a sizable gas reservoir surrounding the binary. For gap-forming binaries we find that orbital evolution is unchanged with the inclusion of feedback, but ionizing radiation photoevaporates gas that is at the outer edge of the low-density region. For non-gap-forming systems we find that when feedback is included a strong initial disruption of the circumbinary disk is followed by an eventual stabilization of the medium that can usher a return to a fast binary migration regime. All of this is possible as a result of how our simulations capture the ionization states of the nuclear disk region and how this affects the coupling efficiency decrease with respect to the radiative feedback.

*Unified Astronomy Thesaurus concepts: Black hole physics (159); Hydrodynamics (1963); Quasars (1319); Radiative transfer simulations (1967); Accretion (14)*

## 1. Introduction

Supermassive black holes (SMBHs), through the conversion of accretion into energy, are considered the source of power behind the observed luminosity in active galactic nuclei (AGNs; Salpeter 1964; Zel'dovich & Novikov 1967). Not only are these the engines behind active nuclei, but SMBHs are also found in the overwhelming majority of galaxies that show the presence of a bulge (Magorrian et al. 1998; Gültekin et al. 2009).

The study of BH binaries and coalescence becomes apparent at first approximation, with the combination of two factors. First, in the current paradigm of the evolution of our Universe, galactic mergers are a common occurrence (White & Frenk 1991). Second, in balanced galaxy mergers the SMBHs will sink by dynamical friction into the nuclear region of the merger remnant (see, e.g., Chandrasekhar 1943; or see Colpi 2014 for a more updated view).

The “paired” BHs in a common bulge mark the first of three stages set by the paradigm first established by Begelman et al. (1980), who first entertained the possibility of BH mergers in galactic nuclei. This is followed by what is sometimes referred to as the “hardening” phase, where the BHs become bound as a binary and have to shrink their orbital separation by other means than dynamical friction (e.g., three-body interactions with the stellar background). The third and final phase is defined by the point at which the BHs are close enough for gravitational radiation to become an efficient mechanism for angular momentum loss (this is expected to happen at

separations comparable to $a_{GW} \sim 10^{-3}(M_{BH}/10^6 M_{\odot})$ pc), which will thereafter drive the binary to coalescence in short order (Peters 1964). BHs in the SMBH range of mass that are merging will be the primary source of study in the resolution range of the Laser Interferometer Space Antenna (Amaro-Seoane et al. 2017).

The work in this paper involves the intermediate hardening phase. Classically, without the presence of gas, with dynamical friction becoming less efficient as a driving mechanism, and with three-body interactions becoming unsustainable at a high enough rate, the orbital shrinkage is stalled, and without additional proposed physical considerations, it enters what is called “the last parsec problem” (Quinlan 1996; Milosavljević & Merritt 2003). The refilling of the loss cone by two-body relaxation may sometimes solve this stalling on sensible timescales, and as of today, mechanisms such as nonaxisymmetric potentials (Berczik et al. 2006) or three-body perturbations (e.g., from a third BH) have been proposed as ways to cope with the last parsec problem, but no definite solution has been found to surmount this major bottleneck in the theory of SMBH binaries embedded in purely stellar environments. It is in this context that the inclusion of the hydrodynamical and gravitational interactions between the binary and a nuclear circumbinary gas disk appears as an additional model ingredient to foster orbital decay.

The possible presence of gas in the merger remnant's nuclear regions is supported by both observational (Sanders & Mirabel 1996) and numerical work (Barnes & Hernquist 1992; Mihos & Hernquist 1995), and in some cases it may be the main driver of a binary merger. The interaction of binaries and gaseous systems has been studied in a variety of contexts and setups (e.g., Goldreich & Tremaine 1980; Lin & Papaloizou 1986; Armitage & Natarajan 2002), and more specifically for BH binaries (Escala et al. 2004, 2005, 2006; Dotti et al. 2007;

<figure><img src="image_4.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Cuadra et al. 2009). To summarize some important takeaways from these works, we have gas coupling with the binary, extracting angular momentum. This is contingent on the gas budget remaining high, which depends on the efficiency of viscous diffusion and its capacity to redistribute this extracted angular momentum. If diffusion is not able to keep up, strong tidal and resonant forces may excavate a low-density gaseous cavity that hampers orbital decay (as the gas budget depletes in the vicinity of the binary).

When this tidal cavity is formed, the systems enters a slow shrinking regime where orbital separation evolves on time-scales of $\frac{1}{a_{bin}}\frac{da_{bin}}{dt} \gtrsim \frac{1}{10^3} t_{orb}$ (del Valle & Escala 2012, 2013; where $a_{bin}$ is the binary's orbital separation and $t_{orb}$ is the orbital timescale). On the other hand, if diffusion is high in comparison to the gap excavation rate, the system will be on a fast-shrinking regime with $\frac{1}{a_{bin}}\frac{da_{bin}}{dt} \sim \frac{1}{t_{orb}}$ (del Valle & Escala 2012, 2013). In this paper we work with setups that are predicted to naturally evolve toward both types of regimes based on the initial conditions selection criterion from del Valle & Escala (2012, 2013), which predicts gap formation (or the lack of such) for balanced binary systems ($\frac{M_{BH1}}{M_{BH2}} = q \sim 1$).

As mentioned, the effects that tidal cavity formation has on orbital shrinkage have been studied extensively, but not with an important consideration given to radiation and feedback. Therefore, with our binary+disk systems, the main objective of this work is to measure in depth the effects that feedback spurred from SMBH accretion has on both the inspiral of the binary and also its impact on disk structure. AGN feedback in the context of SMBH binaries (Lima et al. 2017; del Valle & Volonteri 2018, hereafter VV18) has been found to affect orbital evolution timescales in differing degrees, where feedback processes have been modeled by direct-heating and jet production recipes. As the effects of direct energy injection recipes from AGNs have been found to be very model dependent and highly sensible to coupling calibrations (Thacker et al. 2014; Prieto et al. 2021), we try to improve on this problem by putting specific emphasis on the modeling of radiative feedback processes by coupling radiative transfer physics to the gravity and hydrodynamics, as photons are released through an appropriate emission spectrum.

The paper is organized as follows: In Section 2, we describe the numerical details of the experiment. In Section 3, we show our results without radiation hydrodynamics. In Section 4, we show our results with the coupled radiation hydrodynamics. Finally, in Section 5, we present our discussion and conclusions.

## 2. Methodology and Numerical Simulation

In this paper we model binary+disk systems with a detailed AGN feedback description; for this we perform a suite of simulations using the hydrodynamical adaptive mesh refinement code RAMSES (Teyssier 2002), and more specifically, we employ the radiation-hydrodynamic (RHD) extension RAMSES-RT (Rosdahl et al. 2013; Rosdahl & Teyssier 2015).

RAMSES employs a second-order Godunov scheme for the hydrodynamics, where we used the HLLC solver (Toro 2019) with a MinMod variation diminishing scheme in order to recover the interpolated variables from the cell values that are assigned in the integration scheme. The gravitational field is calculated with finite-difference methods, after a coarse-level

fast Fourier transform (Hockney & Eastwood 2021) for the coarse-level estimation of the Poisson equation.

The additional coupling of radiative transfer is done by approximating the propagation of photons at different frequency bands and computing their interaction with hydrogen and helium through heating, photoionization, and momentum transfer. At each hydrodynamical time step, the pure advection of photons is solved by employing a first-order moment method, using the M1 closure for the Eddington pressure tensor (Rosdahl et al. 2013), which is then solved with a Godunov scheme. This has additional steps in, including the diffusion and multiscattering of IR radiation (Rosdahl & Teyssier 2015). For numerical convergence issues in the RT module, radiation is usually set to travel at a fraction of the speed of light (Rosdahl et al. 2013), for which we set a fraction of $c_r = f_c \cdot c = 5 \times 10^{-3}c$.

Our implementation of SMBHs is done through the use of sink particles and their implementation in RAMSES (as seen in Bleuler & Teyssier 2014). Although in the sink particle paradigm we will be modeling SMBHs as collisionless particles, we do not use the Particle Mesh solver designed for the dark matter component. Instead, a spherical uniform distribution of test particles (“cloud particles”) is placed around each sink. These cloud particles are evenly spaced within the sphere (with roughly eight cloud particles per grid cell) and follow the sink particle as a rigid body. These cloud particles are used to probe and interact with the gas distribution around the sink in order to distribute the accretion and the ejection of mass, momentum, and energy.

### 2.1. Mesh Resolution and Refinement

In the oct-based data structure of RAMSES, the cell width is defined by its refinement level $l$, with $\Delta x = L_{box}^{l-1}$, where our simulation box size is $L_{box} = 260$ pc. We refine from a coarse level of $l_{min} = 7$, which corresponds to a resolution of $\Delta x_{coarse} \sim 4$ pc, up to a fine level of $l_{max} = 16$, corresponding to a maximum resolution $\Delta x_{fine} \sim 0.008$ pc.

For refinement of the mesh we employ three different criteria: (i) We use a quasi-Lagrangian scheme where we set a mass resolution $m_*$ such that the initial x-y plane of the disk is expected to be refined up to level $l = 10$ given the density distribution of the initial disk model (see Section 2.2); cells are then refined when they have mass values larger than $m_*$; (ii) Cell width is also forced at every level to satisfy that the Jeans length is always adequately sampled as per the requirements in Truelove et al. (1997) set to address artificial fragmentation. (iii) Lastly, we enforce that the cells that contain our SMBHs are always at the maximal refinement levels in order to correctly model BH dynamics and their gravitational influence (Lupi et al. 2015).

### 2.2. Initial Conditions

Our simulation suite follows the initial conditions from VV18, who used a criterion that predicts whether hydrodynamical gap formation will occur in binary+disk initial conditions (del Valle & Escala 2012, 2013), this criterion may be parameterized in $((v_{bin}/v), (c_s/v)(h/a_{bin}))$ space, as it reads

$$ \frac{\Delta t_{open}}{\Delta t_{close}} = \frac{1}{f_i} \left(\frac{v}{v_{bin}}\right)^2 \left(\frac{c_s}{v}\right) \left(\frac{h}{a_{bin}}\right). \quad (1) $$
Here $\Delta t_{\text{open}}$ is the gap-opening timescale and $\Delta t_{\text{close}}$ is the gap-filling timescale. In addition, $v$ and $v_{\text{bin}}$ are the binary–disk gas tangential velocity and the binary’s Keplerian velocity, respectively, $c_s$ is the sound speed of the gas, $a_{\text{bin}}$ is the binary separation, and $h$ is the scale height of the disk. Parameter $f_1$ is the relative strength between the gravitational and viscous torques, which itself depends on the geometry of the density perturbation and was fitted in del Valle & Escala (2013; $f_1 = 0.35$). While also expressing the Toomre parameter $Q$ in this phase space, VV18 fixed the setup configurations in a manner such that it also guarantees disk stability (Figure 1 in said work).

With this partition of phase space, we will have four setups, all gravitationally stable against rotation, of which two are expected to form a tidal cavity and two will not. The two setups that will show gap formation are expected to exhibit slow binary migration when no AGN feedback is included, and the two setups that do not form a low-density cavity are expected to maintain a gas budget that will be big enough to promote fast binary migration.

Our disks are embedded in a potential field that imitates the presence of a central stellar bulge, which serves both to stabilize the gaseous disk and to better model the gravitational impact the bulge may have in gas outflows. The bulge gravitational potential is modeled with a Plummer profile (Plummer 1911).

Even though we use a strictly analytic recipe for disk initialization, the potential field of the system including the bulge–disk–binary is nonaxisymmetric, which means that there is no good a priori definition for the circular velocity of the binary–disk system. To cope with this, the binary is first approximated with a spherical mass distribution (with a radius equal to the initial binary separation $a_{\text{bin}}$), and the system is relaxed for ~30 orbits, after which the binary is initialized in the disk. Additional information about the homogeneous spherical potential can be found in del Valle & Escala (2012).

## 2.3. Radiative and Mechanical AGN Feedback Models

For AGN feedback, we first specify how accretion onto the BHs is described and then show how said accretion rates convert to feedback. The most commonly used BH accretion model is that of Bondi accretion (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952), which operates at a characteristic length scale given by $r_{\text{acc}} = \frac{2GM}{c_s(\infty)^2} \approx 1 \times 10^{-4} \left(\frac{M}{M_{\odot}}\right) \left(\frac{10^4 \text{ K}}{T(\infty)}\right)$ pc, which in practice is too coarse for our resolutions. Furthermore, if we also add turbulent velocity to the calculation of the influence radius, the spatial scales will remain orders of magnitude above the length scales of our problem. As we have a resolution that allows us to use a more complex accretion recipe, we then employ a version of threshold accretion similar to the one originally implemented in Bleuler & Teyssier (2014), but with additional checks that are consistent with the relatively resolved thin-disk accretion regime we are dealing with. The gas inside the spatial region defined by the sink particle cloud is accreted if its specific angular momentum falls below the value that characterizes a circular orbit for the biggest Shakura–Sunyaev disk that is stable against self-gravity. This specific angular momentum threshold is set by the value $l_{\text{thr}} = \sqrt{GMR_d}$, where $R_d$ is the characteristic size of the disk

we use for said threshold (Kolykhalov & Sunyaev 1980):

$$R_d = 3R_g r_d \approx 2.64 R_g (\alpha \dot{m}^2 r_0^{1/2} ct / R_g)^{1/4}, \quad (2)$$

with $R_g = \frac{2GM}{c_s^2}$ being the Schwarzschild radius. For $\dot{m}$ we initially assume a nominal set by half of the Eddington accretion rate $\dot{m} = \frac{1}{2}\dot{m}_{\text{Edd}} = \frac{L_{\text{Edd}}}{\epsilon_{\gamma}c^2} = \frac{4\pi GMm_p}{\epsilon_{\gamma}c\sigma_{\tau}}$, and the subsequent rate that is used inside the calculation for $R_d$ is set to be the actual accretion previously outputted by the simulation. If we quote the values of $\alpha = 1$, $r_0 = 9$, and $t = 10^7$ yr from Kolykhalov & Sunyaev (1980), we can readily calculate $R_d$.

The net energy that is released from the BHs in feedback form is proportional to this accretion rate through $L_{\text{bol}} = \eta \dot{m} c^2$, with $\eta = 0.1$ the radiative efficiency parameter (see, e.g., Shakura & Sunyaev 1973). The output of this energy onto the simulation grid is done through two main mechanisms that operate in differing amounts depending on the accretion rates; these two mechanisms are the radiative mode (sometimes known as quasar mode) feedback, which is powered by photons, and the mechanical mode (sometimes known as radio mode) feedback, which is powered by the production of jets perpendicular to the accretion flow.

Mechanical feedback is thought to be relevant when the accretion flow is optically thin and at low accretion values in Eddington units of $\dot{M}/\dot{M}_{\text{Edd}} \lesssim 10^{-3}$ (Narayan et al. 1995; Blandford & Begelman 1999; Yuan & Narayan 2014). Jets are generated perpendicular to the rotation axis of the BH as hydrodynamical source terms that accelerate the immediate environment of the sink particles, and they are implemented as in Dubois et al. (2012), where half of the AGN luminosity is converted to this form (and the other half is maintained at radiative mode).

Radiative feedback mode has already been run in RAMSES-RT (Bieri et al. 2017; Cielo et al. 2018; Costa et al. 2022), but only for galactic-scale simulations, and with a scale-appropriate fixed AGN spectrum from Sazonov et al. (2004). Instead of this, we generate spectra that follow the high-energy AGN distribution models from disk reflection recipes (Pounds et al. 1990; Ross & Fabian 1993), which we deploy from spectrum generation code XILLVER (García et al. 2013). In these models the illuminating X-ray continuum is generated by Compton upscattering of thermal photons by electrons in a hot corona or jet base, which itself is modeled by a power-law spectrum of the form $F(E) = AE^{\Gamma+1} \exp(-E/E_c)$. The power-law index $\Gamma$ and the ionization parameter $\xi = \frac{4\pi F_e}{n_e}$ (Tarter et al. 1969) are two of the biggest factors for our spectrum geometry, and as such we choose to input variable spectra that are updated depending on the values of these two parameters reported by the simulation outputs (we extract the power-law index using the empiric formula from Yang et al. 2015). Further details of our implementation of the radiative spectrum model are found in Appendix A.

This model implementation is a major contrast to the way feedback is modeled in VV18, where luminosity is converted to heating through an analytic approximation from Sazonov et al. (2004), with an additional X-ray pressure term. Vertical winds are also included on the condition of ejecting half of the material that is eligible for accretion in the same way as done in Choi et al. (2012).
<table><caption>**Table 1**
Summary of Main Setup Parameters</caption><thead><tr><th>Name</th><th><i>M</i><sub>bin</sub></th><th><i>M</i><sub>d</sub>/<i>M</i><sub>bin</sub></th><th><i>a</i><sub>bin</sub></th><th>Gap</th><th>AGN</th></tr></thead><tbody><tr><td>GAP-a3</td><td>1</td><td>0.01</td><td>3</td><td>yes</td><td>no</td></tr><tr><td>GAP-a3-AGN</td><td>1</td><td>0.01</td><td>3</td><td>yes</td><td>yes</td></tr><tr><td>GAP-a2</td><td>1</td><td>1</td><td>2</td><td>yes</td><td>no</td></tr><tr><td>GAP-a2-AGN</td><td>1</td><td>1</td><td>2</td><td>yes</td><td>yes</td></tr><tr><td>a2</td><td>0.1</td><td>10</td><td>2</td><td>no</td><td>no</td></tr><tr><td>a2-AGN</td><td>0.1</td><td>10</td><td>2</td><td>no</td><td>yes</td></tr><tr><td>a1</td><td>0.01</td><td>100</td><td>1</td><td>no</td><td>no</td></tr><tr><td>a1-AGN</td><td>0.01</td><td>100</td><td>1</td><td>no</td><td>yes</td></tr></tbody></table>

**Note.** Masses ($M_{\text{bin}}$, $M_{\text{disk}}$) (binary and disk masses) are in $10^6 M_{\odot}$, and $a_{\text{bin}}$ (binary separation) is in parsecs. The Gap column specifies whether initial conditions are predicted to result in gap formation, and AGN specifies whether the simulation employs radiation hydrodynamics and AGN feedback. Simulations a2 and a2-AGN correspond, as in VV18, to runs that start at 7 pc but feedback is turned on at a separation of 2 pc, when we start the comparative analysis.

## 2.4. Simulation Parameters

All of our circumnuclear disks initially follow a Mestel profile radially (truncated to be constant at $R < a_{\text{bin}}$) and an isothermal profile vertically having an $R_{\text{disk}} = 45$ pc radius, $H = 5$ pc, and a homogeneous distribution temperature of $2 \times 10^4$ K. The external potential contribution by the stellar bulge is modeled by a potential that arises from a Plummer sphere profile (Plummer 1911),

$$ \rho_{\text{bulge}}(r) = \frac{3M_b}{4a_b^3} \left(1 + \frac{r^2}{a_b^2}\right)^{-5/2}, \quad (3) $$

with a Plummer length scale of $a_p = 65$ pc and a bulge mass scale $M_b$ consistent with the SMBH–bulge mass relation $M_{\text{BH}} \sim 10^{-3} M_b$ (Häring & Rix 2004).

The simulation suite includes eight different physical setups, with four different initial conditions that are run with or without radiative feedback, as well as different feedback settings. These physical setups separate into two systems that are inside the dynamical cavity-forming regimes in the del Valle & Escala (2012, 2013) paradigm (namely GAP-a2 and GAP-a3) and two systems that are not (systems a2 and a1). The mass of the circumbinary disk is $10^6 M_{\odot}$ for all setups except for one of the cavity-forming conditions, which, along with initial binary separation $a_{\text{bin}}$ and binary mass value, is summarized in Table 1, where we also tag the naming of the respective versions of these setups that will include feedback.

We show in Table 2 the physical quantities that are derived from these setups, which will set the inner physical scales of each system. Here we include the threshold radius $R_d^0$, the Schwarzschild radius $R_g$, the orbital time $t_{\text{orb}}$, and the pure orbital speed $v_{\text{bin}}$ (which is used for the calculation of the real initial circular velocity of the binary or the Toomre parameter). We see in these overall values that the ratio between the unresolved estimated accretion disk radius and the Schwarzschild radius is generally around $R_d^0/R_g \sim 10^3$, which is relevant in illustrating how we will not be able to capture the scales at which the gravitational wave regime starts operating.

Table 3 includes the six photon groups used for the modeling of radiation in the simulation. The photon groups consist of one IR group (0.01–1 eV); one optical group (1–13.5 eV); two groups of ionizing UV photons (UV$_1$ from 24.6 to 54.4 eV and

<table><caption>**Table 2**
Tabulated Values Necessary for the Accretion Rate and Physical Scales</caption><thead><tr><th>Name</th><th><i>R</i><sub>d</sub><sup>0</sup></th><th><i>R</i><sub>g</sub></th><th><i>v</i><sub>bin</sub></th><th><i>t</i><sub>orb</sub></th></tr></thead><tbody><tr><td>GAP-a3</td><td>4.96 × 10<sup>-5</sup></td><td>4.787 × 10<sup>-8</sup></td><td>5.356 × 10<sup>6</sup></td><td>54.78</td></tr><tr><td>GAP-a2</td><td>4.96 × 10<sup>-5</sup></td><td>4.787 × 10<sup>-8</sup></td><td>6.559 × 10<sup>6</sup></td><td>29.82</td></tr><tr><td>a2</td><td>2.79 × 10<sup>-3</sup></td><td>4.787 × 10<sup>-9</sup></td><td>1.109 × 10<sup>6</sup></td><td>617.4</td></tr><tr><td>a1</td><td>1.57 × 10<sup>-7</sup></td><td>4.787 × 10<sup>-10</sup></td><td>9.276 × 10<sup>5</sup></td><td>105.4</td></tr></tbody></table>

**Note.** The physical units are parsecs for $R_d^0$ and $R_g$, cm s$^{-1}$ for $v_{\text{bin}}$, and kiloyears for $t_{\text{orb}}$.

<table><caption>**Table 3**
Properties of the Photon Groups Used in the Simulations</caption><thead><tr><th>Photon</th><th><i>E</i><sub>min</sub></th><th><i>E</i><sub>max</sub></th><th><i>r</i><sub>c</sub></th><th><i>κ</i><sub>sc</sub></th></tr></thead><tbody><tr><td>IR</td><td>0.01</td><td>1</td><td>0</td><td>1</td></tr><tr><td>Optical</td><td>1</td><td>24.6</td><td>10<sup>3</sup></td><td>0</td></tr><tr><td>UV<sub>1</sub></td><td>24.6</td><td>54.4</td><td>10<sup>3</sup></td><td>0</td></tr><tr><td>UV<sub>2</sub></td><td>54.4</td><td>130</td><td>10<sup>3</sup></td><td>0</td></tr><tr><td>Soft X-ray</td><td>130</td><td>10<sup>3</sup></td><td>10<sup>3</sup></td><td>0</td></tr><tr><td>Hard X-ray</td><td>10<sup>3</sup></td><td>10<sup>6</sup></td><td>10<sup>3</sup></td><td>0</td></tr></tbody></table>

**Note.** Columns show the name; minimum and maximum energies, in eV; absorption opacity (Planck), in cm$^2$ g$^{-1}$; and dust-scattering opacity (Rosse-land), in cm$^2$ g$^{-1}$.

UV$_2$ in the band 54.4–130 eV) that have boundaries that coincide with the ionization energy peaks of H I, He I, and He II; and finally two X-ray groups that represent soft X-rays (130–10$^3$ eV) and hard X-rays (1–10$^3$ keV). We specify these energy boundaries and accompany them with the corresponding group absorption opacities and scattering opacities.

## 3. Results without AGN Feedback

We first begin analyzing our runs done without radiation, where we explore binary evolution in a purely hydrodynamical +gravitational paradigm. We also check for consistency with VV18 in order to validate the subsequent comparison of the AGN radiative feedback models. For this posterior analysis we need certain information from the binary+disk dynamics/ interaction in order to ascertain and quantify the impact that our feedback model and the inclusion of photons in the dynamics have.

### 3.1. Non-gap-forming Setups

We start studying the results in setups where no tidal gap is expected to form, namely the systems a1 and a2. The setup parameters have corresponding values of $M_{\text{bin}} = 10^5 M_{\odot}$, $M_{\text{disk}} = 10^6 M_{\odot}$, and an initial SMBH separation of $a_0 = 7$ pc for a2 and $M_{\text{bin}} = 10^4 M_{\odot}$, $M_{\text{disk}} = 10^6 M_{\odot}$, and an initial SMBH separation of $a_0 = 1$ pc for a1.

In the long run we predict no important tidal cavity formation on these setups, but there are some differences: For a2, its place on the parameter space regarding gap formation indicates that the expected perturbation of the disk is nontrivial, as it is close to the regime-partitioning line. In contrast, a1 is expected to show fewer signs of disk perturbation, as viscous diffusion will be very efficient with respect to angular momentum distribution in the circumbinary disk. For both of these setups we expect to see fast orbital migration. This is especially true for a1, as it starts out with a
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Orbital separation of the binary for all our no-feedback setups. Systems a1 and a2 show coalescence in a few units of $t_{\text{orb}}$. Meanwhile, systems GAP-a3 and GAP-a2 show very little deviation from the binary orbit.</figcaption></figure>

much lower initial separation for the binary and a smaller fraction between the BHs and the disk.

In Figure 1 we show the BH separation evolution in time and see how BH coalescence occurs at $1-1.5t_{\text{orb}}$, exhibiting an extremely fast evolution when compared to the system's equilibrium orbital timescales. This fast orbital migration was originally studied in Escala et al. (2004, 2005), where it was shown that it is due to a tidally generated ellipsoidal-shape deformation (for equal-mass massive BHs; Escala et al. 2005) or to a pear-shape deformation (for unequal ones, see Escala et al. 2006), which lags behind the binary with an offset angle, producing the torques responsible for this fast and self-similar migration. Escala et al. (2005, 2006) found that this transition happens when the two spheres of influence ($R_{\text{inf}}$) of gas bound to each BH merge (shown in Figure 2 for run a1), approximately at $1.5 R_{\text{inf}}$ of the secondary (Escala et al. 2005, 2006).

The fast migration mechanism seen in runs a1 and a2 is not restricted to simulations of binary MBHs. This is ubiquitous for binaries in general, when the condition of binary separation $<1.5R_{\text{inf}}$ is fulfilled. For example, in common-envelope simulations similar tidal distortions of ellipsoidal/pear shapes are found (Passy et al. 2011; Ricker & Taam 2012), associated with extremely fast (i.e., dynamical) migration timescales, which are generally confused with an “out-of-equilibrium” configuration (i.e., the “plunge-in” phase; Ivanova & Nandez 2016; Ivanova et al. 2020). This fast and self-similar migration can also be seen as a tidal instability since it can produce a merger in less than one (initial) orbital time, and it is the most likely candidate mechanism to be responsible for a large proportion of binary mergers in the Universe, in a mass range that goes from stellar to supermassive.

## 3.2. Gap-forming Setups

We begin the results analysis for our setups where the initial conditions are primed for tidal gap formation. As expected, these setups do form tidal gaps in a few orbital times and show slow migration times, both for GAP-a3 (which has $M_{\text{bin}} = 10^6 M_{\odot}$, $M_{\text{disk}} = 10^4 M_{\odot}$, and an initial SMBH separation of $a_0 = 3$ pc; see Table 1) and for GAP-a2 ($M_{\text{bin}} = 10^6 M_{\odot}$, $M_{\text{disk}} = 10^6 M_{\odot}$, and an initial SMBH separation of $a_0 = 2$ pc). From how the parameter space is partitioned, we expect a slower gap clearing in GAP-a2 than in GAP-a3 (this is intuitive: the disk's structure is more affected for setups in which the mass of the BHs is higher relative to the medium), as

it is closer to the phase-space partitioning line that separates gap-forming capacity (Figure 1 in del Valle & Volonteri 2018), but it may be more noticeable, as the density contrasts are bigger with the more massive disk present here.

In Figure 1 we confirm how systems without the presence of a low-density cavity exhibit fast migration up to BH coalescence in just a few orbital times, and setups that form a tidal gap show next to no orbital migration, as the gas that will be available for angular momentum removal of the binary is too low.

Both gap-forming and non-gap-forming systems were in complete accordance with the runs done in VV18, which then validates our further analysis and comparisons regarding the inclusion and refining of AGN models.

# 4. Results with AGN Feedback

We now proceed to analyze the same setups done in the previous section, with the addition of AGN feedback and radiation-coupled hydrodynamics. We again separate our study into gap-forming and non-gap-forming initial condition runs, to isolate and measure the effects of our model applied to both contexts separately.

## 4.1. Non-gap-forming Setups

Our non-gap-forming systems are in the regime that is predicted to be affected the most by feedback, by virtue of adding a disruptive mechanism in the disk that could potentially make it so that the system clears out the gas reservoir close to the binary that is needed to remove its angular momentum and promote the binary's merger, in what we saw is a fast migration regime.

In the beginning the perturbation from feedback becomes apparent, specifically through the radiative mode, as a “feed-back bubble” is carved out (see Figure 3, where we show density slices that illustrate the evolution of said bubble). This is driven by two factors: accretion is initially quite rapid, and therefore the feedback’s net luminosity is high; second, as we begin with gas that is fully nonionized, the radiation couples easily to the gas and thus transfers its energy very efficiently. The impact and scale of this radiation bubble differ depending on how massive the disk is in relation to the SMBH mass, where the ability of the system a1-AGN binary’s feedback to perturb the gaseous medium on a disk-wide scale is much weaker than in system a2-AGN. This bubble eventually suppresses luminosity and accretion, until an upturn is observed (in Figure 4 we see the evolution of AGN luminosity in time, where these trends may be observed).

Other clear feedback bubbles are created on smaller scales during this simulation, which usually corresponds to accretion feedback coupling to the medium, as it tends both to mix with the outer nonionized gas from the disk and to smooth out density gradients owing to gas diffusion (see, e.g., the smaller cavities at the nuclear region in Figure 3 at $t = 1.8t_{\text{orb}}$). We see eventually that the gas forms a more or less stable configuration, for instance, we see in Figure 5 (which shows the evolution of mass distribution in the vertical axis) how the z-axis matter distribution stabilizes back into a geometrically thinner distribution after the initial shock from the feedback bubble.

The structure repairing as gas falls back into the disk is timed with how radiation will progressively couple less efficiently
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Final stages leading up to BH coalescence for system a1, where we see the overlapping spheres of influence for the BHs and the characteristic ellipsoidal density enhancement associated with this fast migration.</figcaption></figure>

with the gas as it becomes ionized and the mean free path (MFP) becomes large in the nuclear region where photons are denser. This is relevant when considering that this reopens the possibility that the system will enter a new stable configuration that is non-gap-forming and where the tidal instability mentioned in Section 3 will operate, eventually causing the binary to fall into the characteristic fast coalescence regime.

A decreasing feedback coupling efficiency is evident when looking at how the increasing luminosity rates in Figure 4 at around ~2.5$t_{orb}$ are not matched with a corresponding gas ejection in Figure 3; this could, in principle, be due to radiation leaking through low-density channels in the density distribution, but it is easy to discard this idea, as we see that the density

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Density slices in face-on and edge-on view for run a2-AGN at times $t = 0.2t_{orb}$, $0.8t_{orb}$, $1.8t_{orb}$, and $2.4t_{orb}$. The density scales for these setups are in the density range of $\rho \in (10^{-23}, 10^{-18})$ g cm$^{-3}$. We see the initially blown feedback bubble and its eventual trend back to a stable disk.</figcaption></figure>

distribution becomes progressively smoother and more compact on the X-Y plane; thus, it must be mainly the fact that ionized material is confined to a central disk region up to where recombination/ionization balances out with the radiation field in a way akin to a Strömgren sphere.

We can see the interplay of ionization states and gas density of the disk of simulation a2-AGN at time $t = 3t_{orb}$ in Figure 6, where we see that ionized hydrogen and helium are confined into a central region showing a much sharper density contrast with respect to the actual net gas density. This Strömgren-like region where ionization is contained is important, since it suggests that the disk is stabilizing toward an equilibrium in which there is an inner highly ionized region where the
<figure><img src="image_8.png" /><figcaption>**Figure 4.** AGN luminosity by photon emission bands and total luminosity for run a2-AGN. Absolute emitted luminosity shows values close to $L \sim 10^{42} \text{ erg s}^{-1}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Evolution of the vertical concentration of mass throughout time. We see that, because of feedback, mass is initially vertically diffused on average, after which the disk vertical structure is gradually regained to approximate that of the original setup.</figcaption></figure>

effective cross section is much lower (as such, the corresponding MFP will be much higher) and the corresponding radiation force will be small.

It is useful to look at the MFP map that corresponds to Figure 6, as it gives an actual measure of the expected flow of radiation and how much it is impeded at this ionization threshold region. MPF codifies the distance a photon would be

able to traverse given the information of the cross section $\sigma$ and number density $n$ at a point in space, i.e., $l = \frac{1}{\sigma n}$. As seen in Figure 7, inner regions of the disk have MFPs that are >100 pc, which means that radiation will almost definitely penetrate onto the optically thick regions, which are marked by sharp decreases to low values on the order of $l \sim 10^{-2}$ pc. These low-opacity regions are clearly correlated initially to the low-density excavated region and then with the repairing ionized central region.

In Figure 8, which presents the number density of UV$_2$ photons (which are consistently, alongside X$_2$ photons, the biggest energy contributors), we see more confirmation on the density contrast of photons in the inner ionized region, but we can also appreciate how photons have a very strong tendency to stream freely from the disk in the vertical axis. Another thing we appreciate is how even though it is known that radiation usually prefers to propagate through lower-density regions, this is not primarily evidenced in the face-on slices, as radiation maps tend to have a smoother gradient than density. We will only map this photon group, but it is worth mentioning that were we to map out IR radiation, we would see a much larger encompassing and diffuse photon presence that softly follows the same behavior of ionizing radiation, and if we were to map out hard X-rays, we would basically see no absorption from the gas at all.

At this point we see that feedback by itself may not maintain a “feedback cavity” without the system carving a low-density region by dynamical-only effects already, as it becomes unfeasible when accounting for the diminishing coupling efficiency seen when accounting for ionization. This decreasing coupling efficiency is not captured in models that directly translate luminosity into momentum or heating. As already mentioned, such AGN feedback implementations show a large range of energy coupling efficiencies that are usually calibrated by context-dependent scaling relationships.

In summary, our simulations shows feedback bubbles that become progressively smaller to the point where they become irrelevant after gas in the central region of the disk becomes mostly ionized. After a couple orbital times, we do not see relevant feedback cavities formed and gas stabilization is free to occur. This is explained by the fact that, after a certain level of ionization, coming across big amounts of “coupling-capable” material becomes harder, as mixing between ionized and nonionized gas will not occur violently. This would be challenged, for instance, in the case of having big (neutral) gaseous streams falling into the central regions of the circumbinary disk, which is a possibility in the overall paradigm of nuclear galactic regions (sometimes quoted as a possible facilitating mechanism to promote binary migration), but this is beyond the scope of this paper.

Finally, as we have mentioned, the repair of the disk structure as time passes and coupling efficiency decreases should be accompanied with the binary showing a tendency to coalesce in a fashion that slowly returns to approximating what is seen in Figure 1. In Figure 9, we see the relative decrease (although smaller than what is seen in the corresponding non-radiation run) in the peak binary separation, indicating how the environment is decreasing the BH's angular momentum.

## 4.2. Gap-forming Setups

Gap-forming setups already tend to push out gas by virtue of pure hydrodynamical and gravitational interactions. In general,
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Top row: H II, He II, and He III number density at *t* = 3*t<sub>orb</sub>*. Bottom row: density slices at *t* = 3*t<sub>orb</sub>*. All for run a2-AGN.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** MFP of photons for a2-AGN at *t* = 3*t<sub>orb</sub>*.</figcaption></figure>

feedback mechanisms should aid in the removal of gas in the inner disk’s regions, which will mean that the gas budget should remain low, as for the case without jets and radiation. The fact that this happens indicates that a slow migration regime should be maintained, which we see in the binary separation graph in Figure 14, where there is practically no binary migration.

The original gap formation occurs owing to the tidally induced spiral-wave-patterned wakes the BHs carve onto the disk gas (del Valle & Escala 2012, 2013), which is affected on shorter timescales by the isotropically propagated radiation feedback. This interplay and the isotropic perturbation of the

<figure><img src="image_12.png" /><figcaption>**Figure 8.** UV₂ photons for a2-AGN at *t* = 0.2*t<sub>orb</sub>* and *t* = 3*t<sub>orb</sub>*.</figcaption></figure>

gas mean, at least on first approximation, that gap formation is not as smooth as in the purely hydrodynamical+gravity runs, as the clear spiral-patterned formation of tidal cavity will be disturbed. The fact that accretion is inhibited may be seen by looking at luminosity rates in Figure 10.

We consider in Figure 11 the density slices for system GAP-a2, where we not only see the central low-density region but also specifically see, on the edge-on view, how there is now also a conical-shaped outward ejection flow. This swept-up
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Orbital separation of the binary for a2-AGN.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** AGN luminosity by photon emission bands and total (each panel represents one of the two BHs).</figcaption></figure>

flow is primarily expelled from the low-density cavity's inner limit mostly as singly ionized gas, with a large fraction at speeds that exceed the escape velocity of the system. This steady expulsion of material steadily occurs, as gas, which is not strongly bound to the disk, is not maintained as highly ionized at the cavity's edge, both because there is mixing between ion species in this region and mostly because radiation is too spread out to ionize too deeply into the denser regions of the disk (see Figure 12) and, as such, couples efficiently to the highly energetic radiation. In Figure 12 we show density maps of the ion distributions, where we see how hydrogen ionization traces most expelled material, singly ionized helium traces material that is ejected from the cavity's edge, and finally we get to see some perpendicular double-ionization gas expulsion

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Density slices in X-Y and X-Z view for the GAP-a2 setup at time $t = 2t_{orb}$. The color gradient in the figure goes for $\rho \in (10^{-23}, 10^{-18})$ g cm$^{-3}$.</figcaption></figure>

closer to the binary/disk's center, which comes from the jet component of feedback. We argue that feedback models that only heat up or inject energy directly to the medium should not show the same ejection of boundary material, as they will not generate the same decoupling highly energetic ionizing radiation.

This process in the context of protoplanetary disk theory is commonly known as “disk photoevaporation” (see, e.g., Alexander et al. 2006) and has been widely studied in disk systems with highly energetic and ionizing radiation sources at their centers. If we were to quantify the amount of ejected material at the point of Figures 11 and 12, we would have to measure the dynamically unbound material in the disk (this, at the same time, allows us to correctly see the actual bound
<figure><img src="image_16.png" /><figcaption>Figure 12. Ionization slices for GAP-a2 setup at time $t = 2t_{orb}$.</figcaption></figure>

disk’s structure). First, we need to tag unbound gas that fulfills $v(\vec{r}) \cdot \hat{r} > v_{esc}(r)$, for which we have to approximate the escape velocity (as we do not have a closed analytic expression); we do this by taking

$$v_{esc}(R) = \sqrt{2GM(R)/R}; \quad M(R) = M_{\text{bin}} + M_b(R) + M_d(R),$$

where $M_b(R)$ is the mass inside radius $R$ for our bulge, which comes from a Plummer profile $M_b(R) = M_{\text{Plummer}}(<R) = M_0 \frac{R^3}{(R^2 + a^2)^{3/2}}$, and $M_d(R)$ is the accumulated mass from the disk at radius $R$, which is summed up from the simulation. If we plot the radial density profile with these considerations, we get Figure 13.

The amount of ejected mass is not high enough to pose a structural risk to the system, as disk evaporation in our setups is a self-regulated process, where, as the disk thins out and the gap is markedly formed, accretion will slow down and evaporation will naturally slow down as well. This may also mean that the eventual pseudoequilibrium state of this self-regulated disk will have a bigger gap and thinner disk that we would see without the presence of feedback.

We now finally may see the binary separation graph in Figure 14 for all simulations with feedback, as we did for nonfeedback setups. Here we appreciate how gap-forming setups are indeed not affected in the lack of orbital decay by the inclusion of radiation feedback. Non-gap-forming setups are expected to eventually display sharp orbital decay. This last hypothesis takes a lot of simulation time, as our radiation-coupled setups are quite computationally intensive, and the modeling up to the orbital times needed to capture the full orbital decay is left for future work.

## 5. Discussion/Conclusions

In this work we have studied the effects of AGN feedback in the context of gas-rich SMBH binaries in their transition from middle coalescence stages to their final states, where gravitational waves become important. For this we ran 3D radiation-hydrodynamical simulations with the RAMSES-RT code. There are few simulations of SMBH binaries embedded in gaseous disks at this resolution with the addition of SMBH feedback, and none with coupled radiation, which renders this work a big step in understanding the impact these phenomena have in the evolution of BH binary coalescence.

Our simulation suite took the initial conditions from del Valle & Volonteri (2018), which were selected to be two disk +binary systems in which a tidal gap is expected to happen by gravito-hydrodynamical effects (let us call this group 1) and two systems in which no gap is expected to form (group 2).

1. We begin by going over the results involving group 1 (slow migration setups). The evolution of binary orbit shrinkage without AGN feedback consistently corresponds to their predicted behavior by construction, and little binary separation evolution is seen in the analyzed time window. With the inclusion of feedback, we did not see any relevant changes in their orbital evolution, as although the nature of the low-density region is affected, this region still effectively emerges and lowers the available gas reservoir that will be able to couple to the binary in a way that extracts angular momentum.

Feedback has different effects on the overall disk structure and dynamics for systems in group 1. As mentioned, there is still formation of a low-density central region around the binary, since the outward-exerted pressure should not complicate gap formation, and although the cavity now forms more quickly, it is now formed with less of a clear propagating spiral-wave pattern induced by the BH wakes. Disks in this regime also now show photoevaporation, where material from the inner edge of the tidal cavity is expelled at speeds above escape velocity owing to ionizing radiation coupling very efficiently to gas that is not strongly bound. This mechanism is self-regulated. Some mechanical feedback is seen in these setups, as accretion goes below the threshold, which enables vertical jet formation,
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Density profile when considering unbound material for GAP-a2-AGN at $t = 2t_{orb}$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Orbital separation of the binary for all our feedback-including setups. Systems a1-AGN and a2-AGN show oscillating separations that should eventually start exhibiting orbital decay (which starts for a2-AGN at around three periods). Systems GAP-a3-AGN and GAP-a2-AGN, similarly to their nonfeedback counterparts, show very little deviation from the binary orbit.</figcaption></figure>

but is not a major feature, as the volume of sustained ejected gas is not too high and is less than or comparable to the amounts of gas that are already being unbounded by radiation.

2. For group 2, with setups where no gap is predicted to be excavated and a fast migration is expected, the inclusion of feedback is in the short term potentially disruptive for orbital evolution, but in the long term, after the ionization area is established around the binary, the initial feedback blown bubbles are erased, and the density gradient is “repaired” back toward equilibrium, whereby the gas may be maintained in order to promote BH coalescence. This means that although coalescence is indeed delayed, with the timescales in which the disks stabilize back (i.e., a couple of orbital timescales), a fast migration regime could still be reached. These observations for group 2 are a direct result of how radiative feedback interacts in the long term with gas, which will mean that, just as in our following points regarding the feedback implementation, it shows different behavior than the orbital evolution that is found when implementing a “direct-heating” feedback recipe. One main implication of fast migration not being the end result in this branch of setups is that we now do

not need additional mechanisms for BH coalescence (like the introduction of a third BH or big gas inflows from the outside of the CMD), making coalescence a more natural outcome in the hardening stages of binaries.

The eventual stabilization of the disk due to lowering energy coupling efficiency does not correspond to the findings shown in simulations with AGN recipes that do not include photons (VV18), where group 2 setups behave similarly to tidal gap-forming systems owing to the AGN carving out a “feedback cavity.” We attribute this differing result (and the lack of the photoevaporation mentioned for group 1) to how direct-heating implementations of AGN feedback do not take into account the coupling efficiency of radiation to gas in a consistent manner (e.g., as mentioned in Prieto et al. 2021, order-of-magnitude differences may be seen in how luminosity translates into heating of gas depending on the context and resolution of different simulations) and how ionization states may effectively swing this coupling efficiency depending on the context.

Due to the unique context and the methods employed in this paper, some additional avenues of work and improvements may be highlighted. A main point of interest that is left for work in preparation is how as the introduction of radiation changes accretion flow geometry overall, it impacts the evolution of BH spin magnitude and alignment. This issue comes up in the findings of LIGO merger detections, being a critical constraint on how BH spins align (or potentially counteralign) with respect to the binary orbit angular momentum.

Going beyond our setups, the deployment of RAMSES-RT with a dynamic BH/sink particle emission spectra could be useful in a broader variety of simulations. For instance, for common-envelope setups with binaries having massive objects, radiation effects have been investigated (e.g., Reichardt et al. 2020), but, as per the author's knowledge, not with simulations that employ detailed coupled radiation interactions. Other avenues of research include different contexts in which one would like to see a detailed BH emission spectrum interact with its environment, like in high-resolution simulations of AGN evolution or tidal disruption events.

## Acknowledgments

We are grateful to Luciano del Valle for making his initialization available to us. We thank the anonymous referee for the comments, which helped to improve the manuscript quality. Both J.D. and A.E. would like to acknowledge partial support from the Center for Astrophysics and Associated Technologies CATA (FB210003), Millennium Nucleus NCN19_058 (TITANs), and Proyecto Regular Fondecyt (grant 1181663). Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02).

## Data Availability

The data underlying this article will be shared on reasonable request to the corresponding author.
## Appendix A
AGN Spectrum Generation

In this work the role of AGN feedback takes a central role, especially radiative mode feedback. We now expand further on how this is implemented for our framework.

As we mentioned, AGN radiative feedback has already been run on RAMSES-RT, which has been done by emitting net energy in proportion to the accretion rate to the sink particle through $L_{bol} = \eta \dot{M}c^2$, and then this is distributed in photons by using a fixed AGN emission spectrum (the usual choice of being the generic spectrum from Sazonov et al. 2004). There are two problems with this approach, where it pertains to our work: The first is that emission spectra such as the ones used in the literature are only justifiable for galactic-scale simulation resolutions, for they come as empirically derived averages of both obscured and unobscured spectra that already have reprocessed radiation embedded in their shapes. The other reason is that as we have both a high resolution and an interest in seeing how the system's dynamics interact with feedback, a fixed spectrum is a model oversimplification that we may correct by using a dynamic spectrum that evolves given the state of its respective BH.

The idea of dynamic emission spectra depending on the given sink particle conditions is not novel and has been deployed since RAMSES-RT's creation. This is specifically done for stellar particles, for which the spectral energy distributions (SEDs) are read from a grid of model tables to do on-the-fly evaluation of the respective particle luminosities, using the mass, age, and metallicity of a given star (see Rosdahl et al. 2013). As the code infrastructure is already there, we ported the on-the-fly spectrum calculation method to our BH particles with our own tailored SED model.

The generation model itself comes mainly from a combination of the spectrum generation model XILLVER (García et al. 2013) and, in a smaller measure, the relativistic correction code RELLINE (Dauser et al. 2010). These two may be deployed from the newer code RELXILL (García et al. 2014). This code generates spectra from a disk reflection model that reproduces coronal radiation that is reprocessed in an optically thick disk together with the thermal emission of such a disk. We will not mention the detailed calculations done on such a reflection model (for the classical paradigms see Pounds et al. 1990; Ross & Fabian 1993; for an updated model see Garcia & Kallman 2010), but we do recall the fixed disk, coronal, and BH model features we used to calibrate such a subgrid recipe:

1. The continuum plus line emissivity index from the coronal model is set at the fixed value of $j=3$, which is the default value for a classical $\alpha$-disk paradigm.

2. Second, we specify the inner and outer radii of the disk in the reflection model. The inner radius is taken as the innermost stable circular orbit, and the outer radius is taken as the length $R_d^0$ we employ for our accretion method (Equation (2)).

3. Next, we adopt a BH spin parameter of $\chi=0.25$ (this spin value is consistent with the average parameter values established in King et al. 2008). In practice, varying the spin parameter only relevantly affects the SEDs for extremely high spin parameter values ($\xi \to 1$), where the emission of low-energy photons becomes much higher. The actual modeling of spin magnitude evolution of our BHs will be will be explored in future publications.

4. The energy cutoff for the power-law component of the spectrum is set to be $E_{cut} = 300$ [keV], a value consistent with models such as the ones from Ross & Fabian (2005) for quasar spectrum fittings.

5. Some secondary parameters are set to the default code values, such as the metal abundance $A_{Fe}$, limb darkening, and reflection fraction.

After accounting for all these factors, we are only left with the primary spectrum affecting parameters, namely the ionization parameter $\log(\xi)$ and the power-law index $\Gamma$. The ionization parameter $\xi = \frac{4\pi F_x}{n_e}$ is fixed such that it may be employed for closing and solving the ionization balance equations and therefore determining gas structure. The power-law index is used to define the spectrum of the upscattering thermal photons characterized by a power law of the form $F(E) \propto E^{\Gamma+1} \exp(-E/E_c)$.

The values for $\Gamma$ in quasars are usually found to be $1.8 \lesssim \Gamma \lesssim 2.2$ but may be found to be closer to 1 or 3 depending on conditions (our grid will take into account values in $\Gamma \in (1.3, 2.5)$), and ionization parameters are realistically bound between $\log(\xi) \in (0, 4)$.

Let us consider the variation of the $\Gamma$ parameter for the spectrum generation. We advocate following the evolution of the $\Gamma$ parameter through computing the total X-ray luminosity $L_X$ from the BH (scaled to Eddington units). For this we use the empirical relationship from Yang et al. (2015):

$$ \Gamma = (0.31 \pm 0.01)\log_{10}(L_x/L_{Edd}) + (2.48 \pm 0.02), \\ \forall \left(\frac{L_x}{L_{Edd}}\right) \ge 0.01 $$

$$ \Gamma = (-0.1 \pm 0.02)\log_{10}(L_x/L_{Edd}) + (1.27 \pm 0.03). $$

With this formula, we can graph the energy-scaled spectrum for 10 different luminosity fractions, and therefore different indices, at $\log(\xi) = 2$.

Figure 15 shows how for lower accretion, after the Compton bump (at around and above 2–10 keV), there is a bigger fraction of energy in the hard X-ray band compared to the higher fraction of energy in the soft X-rays for the higher accretion rates. We do not show it here, but this behavior repeats for all our binned $\xi$ values. Let us see now how, for a fixed $\Gamma$ power-law index, the spectral distributions change with a varying $\log(\xi) \in (0, 4)$, shown in Figure 16.

The variation of this parameter has a high overall effect in the shape of the spectrum. Line features are more prominent with lower values of $\log(\xi)$, getting up to the point of getting a very smooth distribution for values nearing the upper model limit of 4. To determine the ionization index from the simulation state, we use the fact that RT-enabled simulations track the electron number density, and as such, we only need the X-ray flux from the AGN to recover the needed value by modifying the formula to

$$ \xi = \frac{4\pi F_x}{n_e} = \frac{4\pi}{n_e} \frac{L_x}{4\pi r_{\text{sink}}^2} = \frac{L_x}{n_e r_{\text{sink}}^2}. $$

Since the correlation of $\Gamma$ and $\log(\xi)$ is big, we create the varying SED grid from binning at the $\Gamma \in (1.3, 2.5)$ and $\xi \in (0, 4)$ intervals, using a $10 \times 10$ grid. Note that for both of these we get the energy fraction of the X-ray photons $L_X$ directly from the spectrum emission we are assuming implicitly, meaning that when going from $t_{n \to n+1}$, we need a spectrum
<figure><img src="image_19.png" /><figcaption>Figure 15. Grid of 10 spectra in $E \cdot F(E)$ form with $\Gamma \in [1.76, 2.4]$. Vertical lines represent the photon energy bins.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 16. Grid of 10 spectra in $E \cdot F(E)$ form with $\log(\xi) \in [0, 4]$.</figcaption></figure>

shape depending on ($\Gamma_n, \xi_n$) to estimate ($\Gamma_{n+1}, \xi_{n+1}$). In practice for quick calculations, assuming that the X-ray fraction is 50% of the total emitted feedback flux is a good approximation. The list of our binned template spectra, energy fractions for the different photon groups, and the energy-weighted cross sections to ionization by H I, He I, and He II will not be included in this paper but may be provided upon request.

## Appendix B
### Stability and Thermal States of the Disk Setups

As a first approximation for stability, Toomre's classical stability parameter $Q = \frac{\zeta_n \kappa}{\pi \Sigma G}$ (Toomre 1964) does a fair job at evaluating the local linear stability of self-gravitating disks by only looking at the local dynamics, playing also a role in setting the mass scale in the formation of stellar clusters in

galactic disks (Escala & Larson 2008) and the star formation triggering at galactic scale (Chavarría et al. 2010; Escala 2011). We plot the profiles for this parameter in setups GAP-a2 and a2 in Figure 17. We see in these graphs (which maintain the same scales for our other setups) that, beyond any radial tendencies we observe, values stay exceedingly above the stability threshold $Q > 1$. This possibly stems from the fact that in general the nonlinear overall behavior of stability in disks is more complex than this, for instance, it is usual to treat the timescale on which the disk can radiate its thermal energy as inversely proportional to stress (Gammie 2001). This implication of the local description of thermal equilibrium gives a new stability parameter $\beta = \frac{\Omega_{\text{cool}}}{3}$. Now, this criterion is still limited, as it assumes that the dominant heating/cooling processes of the disk are spurred by turbulence, and just as with the Toomre parameter, it has been shown to be inadequate at predicting
<figure><img src="image_21.png" /><figcaption>**Figure 17.** The Toomre's $Q$ stability parameter as a function of the radial distance $R$ to the center of the disk. The left panel shows $Q$ for the GAP-a2 simulation as a function $R$ at two different times, with (RT-run) and without (Hydro-run) radiative transfer. The right panel shows $Q$ for the a2 simulation at two different times again, also with and without radiative transfer. In all eight cases, we observe that values stay exceedingly above the stability threshold $Q > 1$ at all radii (i.e. stable beyond any radial variation).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 18.** Cooling rate radial profiles with the separated contribution of different processes for simulation a2-AGN at times $0.2t_{orb}$ and $1.5t_{orb}$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 19.** Cooling rate radial profiles with the separated contribution of different processes for simulation GAP-2-AGN at times $1t_{orb}$ and $8t_{orb}$.</figcaption></figure>

fragmentation in disks with radiative transfer (Tsukamoto et al. 2015; Suazo et al. 2019). We retain the main idea of comparing the local heating rates to the cooling rates for the definition of a consistent stability criterion, where we will have the cooling time due to radiative losses and heating due to photons and viscous diffusion. The viscous dissipation rate is given by $Q\nu ds \approx rH\tau_{r\phi} \cdot d\Omega/dr \approx (9/4)\Omega^2\nu\Sigma$. With this, we will finally look at the ratio

$$ \Pi = \frac{t_{\text{cool}}}{t_{\text{heat}}} = \frac{\epsilon}{\dot{\epsilon}_{-}} \cdot \frac{\dot{\epsilon}_{+}}{\epsilon} = \frac{\dot{\epsilon}_{+}}{\dot{\epsilon}_{-}} $$

Before checking how the cooling/heating fraction behaves, let us first look at how cooling evolves on isolated terms. If we

were to graph each cooling rate $\Lambda$ [erg s⁻¹] associated with the cooling rates from the processes found in Katz et al. (1996) (we separate the rates depending on the ion species they are associated with, excluding bremsstrahlung and Compton cooling), for simulation a2-AGN we get the result shown in Figure 18.

At the same time, if we plot the same cooling contributions for GAP-a2-AGN, we get Figure 19.

With these graphs in hand, we may observe first how there is a clear indication of how, in terms of strength, processes that are associated with hydrogen (especially H1) are the main contributors to radiative cooling, whereas it is clear that Compton cooling is consistently the least important mechanism
<figure><img src="image_24.png" /><figcaption>**Figure 20.** Different radial profiles for our Π parameter for simulation a2-AGN. The image on the right is the same as that on the left, but zoomed in to 0–10 pc.</figcaption></figure>

(excepting instances where cooling ionized helium mechanisms on the outer regions of the disk is extremely low). We also see in the non-gap-forming system how, in the later stages, ionization fraction plays a big role. Between 5 and 10 pc, where the “Strömgren” area of ionization is located, overall changes not only to the rates but also to the separation between mechanisms are introduced (nonionized associated processes increase much more sharply than ionized associated processes). The final and probably most important takeaway is how cooling as an overall net physical process slows down as time progresses and, at the same time, becomes less centrally concentrated.

If we now move toward checking how cooling and heating balance out, let us see first how radial profiles for simulation a2-AGN look in Figure 20.

We first notice how our heating/cooling fraction (Π) is above 1 at every disk radius, in some ways following the line set by our exploration of the Toomre parameter, which was comfortably over the stability threshold. We see that our ionized gas region, where photons flow almost free, shows a sharp contrast to the nonionized region in terms of Π, where the parameter reaches values in the range of Π ∈ [1, 10]. The fact here that this fraction nears unity at the inner regions is consistent with how radiation is not an efficient heating vehicle when it tries to couple with ionized material. It is still above the threshold, which makes sense, as no fragmentation occurs, and we see an overall lowering of the region's values as photoheating becomes less efficient and gas begins falling back toward the disk (at $t = 0.75t_{orb}$, gas at the nuclear disk regions is still in its AGN-blown state, and as such quite hot). As threshold values are kept above unity and fragmentation does not occur, we can only ascertain for the moment that our criterion is consistent with this reality, but this still opens up an avenue for future work, where we could simulate setups that live near the fragmentation phase-space interface and then test the consistency of using thermal energy balances to quantify the stability tendencies of a system and potentially calibrate an actual stability threshold through Π.

References

Alexander, R. Clarke, C., & Pringle, J. 2006, MNRAS, 369, 229

Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv:1702.00786

Armitage, P. J., & Natarajan, P. 2002, ApJ, 567, L9

Barnes, J. E., & Hernquist, L. 1992, ARA&A, 30, 705

Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980, Natur, 287, 307

Berczik, P., Merritt, D., Spurzem, R., & Bischof, H.-P. 2006, ApJ, 642, L21

Bieri, R., Dubois, Y., Rosdahl, J., et al. 2017, MNRAS, 464, 1854

Blandford, R. D., & Begelman, M. C. 1999, MNRAS, 303, L1

Bleuler, A., & Teyssier, R. 2014, MNRAS, 445, 4015

Bondi, H. 1952, MNRAS, 112, 195

Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273

Chandrasekhar, S. 1943, ApJ, 97, 255

Chavarría, L., Mardones, D., Garay, G., et al. 2010, ApJ, 710, 583

Choi, E., Ostriker, J. P., Naab, T., & Johansson, P. H. 2012, ApJ, 754, 125

Cielo, S., Bieri, R., Volonteri, M., Wagner, A. Y., & Dubois, Y. 2018, MNRAS, 477, 1336

Colpi, M. 2014, SSRv, 183, 189

Costa, T., Arrigoni Battaia, F., Farina, E. P., et al. 2022, MNRAS, 517, 1767

Cuadra, J., Armitage, P., Alexander, R., & Begelman, M. 2009, MNRAS, 393, 1423

Dauser, T., Wilms, J., Reynolds, C., & Brenneman, L. 2010, MNRAS, 409, 1534

del Valle, L., & Escala, A. 2012, ApJ, 761, 31

del Valle, L., & Escala, A. 2013, ApJ, 780, 84

del Valle, L., & Volonteri, M. 2018, MNRAS, 480, 439

Dotti, M., Colpi, M., Haardt, F., & Mayer, L. 2007, MNRAS, 379, 956

Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662

Escala, A. 2011, ApJ, 735, 56

Escala, A., Larson, R., Coppi, P., & Mardones, D. 2006, RMxAC, 26, 141

Escala, A., & Larson, R. B. 2008, ApJ, 685, L31

Escala, A., Larson, R. B., Coppi, P. S., & Mardones, D. 2004, ApJ, 607, 765

Escala, A., Larson, R. B., Coppi, P. S., & Mardones, D. 2005, ApJ, 630, 152

Gammie, C. F. 2001, ApJ, 553, 174

García, J., Dauser, T., Lohfink, A., et al. 2014, ApJ, 782, 76

García, J., Dauser, T., Reynolds, C., et al. 2013, ApJ, 768, 146

Garcia, J., & Kallman, T. R. 2010, ApJ, 718, 695

Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425

Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 698, 198

Häring, N., & Rix, H.-W. 2004, ApJ, 604, L89

Hockney, R. W., & Eastwood, J. W. 2021, Computer Simulation Using Particles (Boca Roton, FL: CRC Press)

Hoyle, F., & Lyttleton, R. A. 1939, PCPS, 35, 405

Ivanova, N., Justham, S., Ricker, P., et al. 2020, Common Envelope Evolution (Bristol: IOP Publishing)

Ivanova, N., & Nandez, J. L. A. 2016, MNRAS, 462, 362

Katz, N., Weinberg, D. H., & Hernquist, L. 1996, ApJS, 105, 19

King, A., Pringle, J., & Hofmann, J. 2008, MNRAS, 385, 1621

Kolykhalov, P., & Syunyaev, R. 1980, SvAL, 6, 357

Lima, R. S., Mayer, L., Capelo, P. R., & Bellovary, J. M. 2017, ApJ, 838, 13

Lin, D. N., & Papaloizou, J. 1986, ApJ, 309, 846

Lupi, A., Haardt, F., & Dotti, M. 2015, MNRAS, 446, 1765

Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285

Mihos, C., & Hernquist, L. 1995, arXiv:astro-ph/9512099

Milosavljević, M., & Merritt, D. 2003, ApJ, 596, 860

Narayan, R., Yi, I., & Mahadevan, R. 1995, Natur, 374, 623

Passy, J.-C., Marco, O. D., Fryer, C. L., et al. 2011, ApJ, 744, 52

Peters, P. C. 1964, PhRv, 136, B1224

Plummer, H. C. 1911, MNRAS, 71, 460

Pounds, K., Nandra, K., Stewart, G., George, I., & Fabian, A. 1990, Natur, 344, 132

Prieto, J., Escala, A., Privon, G. C., & d'Etigny, J. 2021, MNRAS, 508, 3672

Quinlan, G. D. 1996, NewA, 1, 35

Reichardt, T. A., De Marco, O., Iaconi, R., Chamandy, L., & Price, D. J. 2020, MNRAS, 494, 5333
Ricker, P. M., & Taam, R. E. 2012, ApJ, 746, 74

Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188

Rosdahl, J., & Teyssier, R. 2015, MNRAS, 449, 4380

Ross, R., & Fabian, A. 1993, MNRAS, 261, 74

Ross, R., & Fabian, A. 2005, MNRAS, 358, 211

Salpeter, E. E. 1964, ApJ, 140, 796

Sanders, D., & Mirabel, I. 1996, ARA&A, 34, 749

Sazonov, S., Ostriker, J., & Sunyaev, R. 2004, MNRAS, 347, 144

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Suazo, M., Prieto, J., Escala, A., & Schleicher, D. R. G. 2019, ApJ, 885, 127

Tarter, C. B., Tucker, W. H., & Salpeter, E. E. 1969, ApJ, 156, 943

Teyssier, R. 2002, A&A, 385, 337

Thacker, R. J., MacMackin, C., Wurster, J., & Hobbs, A. 2014, MNRAS, 443, 1125

Toomre, A. 1964, ApJ, 139, 1217

Toro, E. F. 2019, ShWav, 29, 1065

Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179

Tsukamoto, Y., Takahashi, S. Z., Machida, M. N., & Inutsuka, S.-i. 2015, MNRAS, 446, 1175

White, S. D., & Frenk, C. S. 1991, ApJ, 379, 52

Yang, Q.-X., Xie, F.-G., Yuan, F., et al. 2015, MNRAS, 447, 1692

Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529

Zel'dovich, Y. B., & Novikov, I. D. 1967, SvA, 10, 602