A unified accretion disc model for supermassive black holes in galaxy formation simulations: method and implementation

ABSTRACT

It is well established that supermassive black hole (SMBH) feedback is crucial for regulating the evolution of massive, if not all, galaxies. However, modelling the interplay between SMBHs and their host galaxies is challenging due to the vast dynamic range. Previous simulations have utilized simple subgrid models for SMBH accretion, while recent advancements track the properties of the unresolved accretion disc, usually based on the thin α-disc model. However, this neglects accretion in the radiatively inefficient regime, expected to occur through a thick disc for a significant portion of an SMBH’s lifetime. To address this, we present a novel ‘unified’ accretion disc model for SMBHs, harnessing results from the analytical advection-dominated inflow–outflow solution (ADIOS) model and state-of-the-art general relativistic (radiation-)magnetohydrodynamics (GR(R)MHD) simulations. Going from low to high Eddington ratios, our model transitions from an ADIOS flow to a thin α-disc via a truncated disc, incorporating self-consistently SMBH spin evolution due to Lense–Thirring precession. Utilizing the moving mesh code arepo, we perform simulations of single and binary SMBHs within gaseous discs to validate our model and assess its impact. The disc state significantly affects observable luminosities, and we predict markedly different electromagnetic counterparts in SMBH binaries. Crucially, the assumed disc model shapes SMBH spin magnitudes and orientations, parameters that gravitational wave observatories like LISA and IPTA are poised to constrain. Our simulations emphasize the importance of accurately modelling SMBH accretion discs and spin evolution, as they modulate the available accretion power, profoundly shaping the interaction between SMBHs and their host galaxies.

1 INTRODUCTION

Galaxy formation in the ΛCDM universe is a hierarchical process whereby primordial density fluctuations collapse into dark matter haloes with cosmic time and grow via inflows and mergers.The gravitational pull of the dark matter leads the gas to collapse within these haloes, where it can radiatively cool and form stars. To prevent excessive star formation and reproduce the observed galaxy stellar mass function, theoretical models need to include feedback processes (see e.g. Somerville & Davé 2015; Vogelsberger et al. 2020, for recent reviews). Supernova feedback is commonly invoked for the regulation of the baryon cycle in low-mass galaxies, whilst the more energetic active galactic nuclei (AGN) feedback from supermassive black holes (SMBHs) is required to suppress star formation in massive galaxies.

Recent tantalizing observations suggest that AGN feedback could also play a role at the low-mass end of the galaxy population, with observational hints at AGN-driven outflows in nearby dwarfs (e.g. Penny et al. 2018; Manzano-King, Canalizo & Sales 2019; Liu et al. 2020; Davis et al. 2022; Aravindan et al. 2023) as well as indications of significant AGN activity from overmassive BHs in high-redshift low-mass galaxies (e.g. Harikane et al. 2023; Maiolino et al. 2023, 2024; Übler et al. 2024).

SMBHs are likely to reside in the centre of the majority of massive galaxies (see e.g. review by Kormendy & Ho 2013) and intermediate-mass black holes (IMBHs) may reside in the centres of the majority of massive dwarfs (Greene, Strader & Ho 2020). As these galaxies are expected to undergo several major mergers with similarly sized galaxies over their cosmic lifetime, IMBHs and SMBHs also grow via mergers in addition to gas accretion. The merging process of stellar-mass BHs has been recently directly observed with the Advanced LiGO and VIRGO detectors which have opened a gravitational wave window on our Universe, even likely pushing into the IMBH regime (Abbott et al. 2020). While these ground-based detectors are not sensitive to the range of gravitational waves emitted by merging SMBHs, IPTA and LISA in future will be able to probe them. This potential has been spectacularly demonstrated by the detection of a signal consistent with a gravitational wave background from SMBH mergers by members of IPTA (Verbiest et al. 2016; Agazie et al. 2024), NANOGrav (Agazie et al. 2023a), Parkes PTA (Reardon et al. 2023) and joint results from European and Indian PTAs (Antoniadis et al. 2023) marking the onset of the multimessenger era for galaxy formation.

SMBH feedback in galaxy formation is an inherently multiscale process spanning 14 orders of magnitude from the event horizon (∼10⁻⁶ pc for Sgr A*) to the cosmic web (∼10⁸ pc). This renders an ab-initio treatment computationally infeasible and hence cosmological simulations have to resort to including SMBH physics via so-called subgrid models (also see discussion in Curtis & Sijacki 2015). The modelling of the gas accretion onto the AGN is essential for obtaining realistic AGN feedback as the luminosity is directly tied to the accretion rate. What is more, accretion disc physics is crucial for accurately modelling the SMBH spin evolution which dictates the radiative efficiency and jet powers as well as influences gravitational recoils.

Most cosmological simulations (e.g. Dubois et al. 2014; Schaye et al. 2015; Sijacki et al. 2015; Pillepich et al. 2018) have adopted the ‘Bondi–Hoyle–Lyttleton’ model (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952) for black hole (BH) accretion, which allows us to infer SMBH growth rates based on the gas properties at the Bondi radius (∼5–5000 pc for SMBHs), significantly reducing the resolution requirements. Though note that many cosmological simulations do not even resolve the Bondi radius, especially for low-mass SMBHs, though additional refinement in the central region can circumvent these issues (e.g. Curtis & Sijacki 2016a). However, the Bondi model is very simplistic, assuming radial symmetry, neglecting angular momentum transfer, and cannot be self-consistently used to predict luminosities of quasars. Moreover, due to the strong dependence of the Bondi accretion rate on SMBH mass, the Bondi model makes it very difficult to grow light seeds or reproduce observations of bright AGN in dwarfs (e.g. Koudmani, Henden & Sijacki 2021; Haidar et al. 2022; Koudmani, Sijacki & Smith 2022), whilst likely overestimating the gas accretion rate for elliptical galaxies and galaxy clusters (e.g. Russell et al. 2015; Bambic et al. 2023; Guo et al. 2023). To address these limitations more advanced approaches have emerged to model AGN accretion processes following two main strategies.

In the first approach, the range of spatial or temporal scales is reduced to make the problem computationally tractable. Recent examples include the hydrodynamic simulations of an elliptical galaxy by Guo et al. (2023), resolving SMBH accretion from galaxy scales all the way down to scales similar to the black hole horizon, and recent general relativistic magnetohydrodynamics (GRMHD) simulations bridging the scales from the event horizon to the Bondi radius and beyond (Lalakos et al. 2022; Cho et al. 2023). Bridging from large scales inwards, the cosmological-zoom in simulation of a high-redshift quasar by Hopkins et al. (2024a, b) follow the gas flows from the cosmological environment to the scale of the SMBH’s accretion disc (∼10⁻⁴ pc) for ∼10⁴ yr. While these attempts can reach an impressive dynamical range, they are computationally prohibitively expensive to allow for a study of a representative SMBH population.

The second approach is to improve the physical realism of the subgrid models. For example, extensions to the Bondi model have been developed to take the gas angular momentum into account in hydrodynamical simulations, e.g. Krumholz, McKee & Klein (2005), Rosas-Guevara et al. (2015), or Curtis & Sijacki (2016b). Furthermore, the torque-based accretion model infers the gas accretion rates motivated by analytical calculations and numerical simulations of angular momentum transport and gas inflow in galaxies, from scales of ∼ kpc to deep inside the potential of the central SMBH (Hopkins & Quataert 2011; Anglés-Alcázar, Özel & Davé 2013; Anglés-Alcázar et al. 2017), recently extended by Rennehan et al. (2023) to include AGN feedback in different physical accretion regimes. However, the torque-based model does not track small-scale angular momentum flows which are crucial to model the SMBH spin evolution, which modulates the radiative efficiency and jet power.

As a compromise between these two main approaches, the accretion disc particle method has been introduced by Power, Nayakshin & King (2011) and further developed by a range of groups (e.g. Dubois et al. 2014; Fiacconi, Sijacki & Pringle 2018; Yuan et al. 2018; Beckmann et al. 2019; Bustamante & Springel 2019; Cenci et al. 2021; Huško et al. 2022; Tartėnas & Zubovas 2022; Bollati et al. 2023b; Massonneau et al. 2023; Huško et al. 2024), including self-consistent evolution of the SMBH spin. The common approach here is to significantly increase the resolution around the SMBH to directly measure the mass and angular momentum flows onto the accretion disc, typically requiring us to resolve scales of ∼10⁻² pc corresponding to the outer radius of the accretion disc. The accretion disc and SMBH properties are then evolved in a subgrid fashion. These models build on the rich body of work which has focused on developing AGN accretion disc prescriptions for semi-analytical models of galaxy formation (e.g. Volonteri et al. 2005; Berti & Volonteri 2008; Lagos, Padilla & Cora 2009; Fanidakis et al. 2011; Barausse 2012; Dotti et al. 2013; Volonteri et al. 2013; Sesana et al. 2014; Gaspari & Sa̧dowski 2017; Lagos et al. 2024).

However, these accretion disc models mostly focus on the thin disc case (the high-Eddington-ratio regime) for the steady state as this regime can be described analytically with the Shakura-Sunyaev disc model (Shakura & Sunyaev 1973), whilst there is no equivalent global analytical disc model for the low-Eddington-ratio regime (though see Fanidakis et al. 2011; Yuan et al. 2018). Recent cosmological simulations projects have started incorporating thick discs into their frameworks (e.g. Dubois et al. 2021; Huško et al. 2022, 2024), however, these generally have to make some approximations based on the Bondi rate as the mass and angular momentum flows onto the disc are not directly resolved.

Observations of X-ray binaries and AGN indicate that at low Eddington ratios (f_Edd ≲ 0.01) accretion occurs via a different mode than the standard thin disc, e.g. observations show hard X-ray spectra instead of soft black body-like spectra and low radiative efficiencies. The ADAF (advection dominated accretion flow) solution has been shown to have the properties required to provide a physical description of the observations (e.g. Narayan & McClintock 2008). Moreover, most X-ray binaries in the hard and quiescent state are deduced to accrete via the so-called truncated discs with an outer thin disc and an inner ADAF component (e.g. Yuan & Narayan 2014). There is also some observational evidence that this phenomenon could transfer over to AGN (e.g. Trump et al. 2011; Yu, Yuan & Ho 2011; Nemmen, Storchi-Bergmann & Eracleous 2014).

Theoretical arguments suggest that ADAFs should produce strong winds and this outflow behaviour can be modelled with the general advection-dominated inflow–outflow solution (ADIOS) for radiatively inefficient accretion flows (Blandford & Begelman 1999).

Building on these observations and theoretical frameworks, as well as results from general relativistic (radiation-)magnetohydrodynamics (GR(R)MHD) simulations, we have developed a unified accretion disc subgrid model that proceeds via the thin disc model from Fiacconi et al. (2018) at high Eddington ratios and via the ADIOS flow model at low Eddington ratios. For intermediate Eddington ratios, we smoothly model the transition between these two modes via a truncated disc. A unified model for accretion through a truncated thin disc and inner ADIOS flow does not currently exist for cosmological simulations, so our new method is an important advancement for modelling accretion onto AGN in general.

In this paper, we present the numerical method and implementation as well as a suite of validation simulations, including single SMBHs and SMBH binaries in gaseous discs. The latter crucially demonstrates the scientific relevance of our model as electromagnetic counterpart and spin evolution predictions are hugely sensitive to the disc state. Due to significant observational challenges there are currently only a handful of confirmed binary detections (Rodriguez et al. 2006; Bansal et al. 2017; Kharb, Lal & Merritt 2017). High-resolution simulations of SMBH binaries (Dittmann & Ryan 2022; Franchini, Lupi & Sesana 2022; Bollati et al. 2023a; Bourne et al. 2023; Siwek et al. 2023a; Siwek, Weinberger & Hernquist 2023b; Dittmann & Ryan 2024) are crucial as next-generation facilities are set to significantly extend these samples with dedicated surveys across the electromagnetic spectrum (D’Orazio & Charisi 2023).

The remainder of this paper is structured as follows. Section 2 presents the theoretical background as well as the main equations governing the unified accretion disc model. The simulation suite is presented in Section 3 and the results, with a special focus on a comparison between the thin disc model and unified disc model predictions, are presented in Section 4. We discuss the observational implications and possible extensions of our model in Section 5 and conclude in Section 6.

2 METHODOLOGY

In this Section, we provide the theoretical background as well as derivations of the relevant equations for our unified accretion disc model. Readers who would just like to have an overview of the methods without delving into the fine details are advised to just read Sections 2.1 and 2.6, in conjunction with Fig. 1 and Table 1.

2.1 Numerical implementation overview

Our unified accretion disc model evolves the SMBH mass, M_•, and angular momentum, J_•, as well as the accretion disc mass, M_d, and angular momentum, J_d, in a subgrid fashion, where the SMBH and its accretion disc are modelled as a composite SMBH–accretion disc particle (see Fig. 1, upper panel). This builds on previous accretion disc based SMBH growth models in galaxy formation (e.g. Power et al. 2011; Dubois et al. 2014; Fiacconi et al. 2018; Bustamante & Springel 2019; Cenci et al. 2021; Tartėnas & Zubovas 2022). In the majority of these cases, the thin α-disc model is used to evolve the SMBH and accretion disc properties where the disc is assumed to be in local thermal equilibrium, and can radiate its heat efficiently (Shakura & Sunyaev 1973). The kinematic viscosity ν is parametrized as ν = αc_sH, where α is a dimension-less parameter representing the efficiency of viscous transport, c_s is the sound speed, and H is the scale height. As the thin α-disc offers a global, steady-state analytical solution and allows for the self-consistent calculation of rest-mass-energy and angular momentum transfer at the innermost stable circular orbit (ISCO), it is the commonly chosen for subgrid accretion disc prescriptions.

However, the thin α-disc model makes several simplifying assumptions about the nature of the accretion flow – in particular, it is only valid in the radiatively efficient regime. At low Eddington ratios, in the radiatively inefficient regime, the accretion flow is expected to transition to a geometrically thick disc following the ADAF solution. Narayan & Yi (1995) noted that ADAFs will likely be associated with strong winds and Blandford & Begelman (1999) derived a family of self-similar solutions, advection-dominated inflow–outflow solutions ‘ADIOS’, with a wide range of outflow efficiencies. Analytical models on their own cannot predict the wind loss in the advection-dominated regime and hence the wind parameters have to be informed by numerical simulations. In this work, we aim to extend the thin α-disc model implementation presented in Fiacconi et al. (2018) to the radiatively inefficient regime based on the analytical ADIOS model and results from high-resolution GR(R)MHD simulations, with the transition between the two disc states modelled as a truncated accretion disc.

With our unified accretion disc model, the gas mass and angular momentum inflows at the scale of the outer accretion disc are directly measured from the hydro solver and added appropriately to the subgrid accretion disc. The properties of the SMBH and its surrounding accretion disc are then evolved by the subgrid model (see Fig. 1).

First, we summarize the relevant equations for the mass and spin evolution of the accretion disc particle system. The SMBH mass, M_•, evolves as:

where |$\dot{M}_\mathrm{BH,0}$| is the rest mass flux onto the SMBH and η is the radiative efficiency. The disc mass, M_d, evolution is described by:

where |$\dot{M}$| is the steady state rest mass flux through the disc and |$\dot{M}_{\rm in}$| is the gas mass inflow onto the disc (see Section 2.2.3 for details on how this quantity is calculated by the hydro solver). Note that if there is no mass-loss via a jet or a wind then |$\dot{M} = \dot{M}_\mathrm{BH,0}$|⁠. Here we distinguish between these two quantities because we model mass-loss through a wind in the radiatively inefficient regime, following the ADIOS model.

The angular momentum of the SMBH, J_•, evolves both due to accretion and due to Lense–Thirring precession (if misaligned from the disc angular momentum vector J_d):

where L_inner is the (spin-dependent) specific angular momentum at the inner boundary of the accretion flow and j_• and j_d are the angular momentum unit vectors (versors). The angular momentum of the disc then evolves as:

where |$\dot{\boldsymbol {J}}_{\rm in}$| is the gas angular momentum flow onto the disc, see Section 2.2.3 for a more detailed definition.

Here we assume that the outflows are occurring in the innermost region so that the winds carry away the specific angular momentum associated with the inner boundary of the accretion flow, which represents a lower limit for the angular momentum loss. Furthermore, we note that in the truncated disc state, there would be an additional torque between the inner ADIOS flow and outer thin disc. Fully accounting for this effect is beyond the scope of this work, however, we include the impact of the outer thin disc on the precession frequency of the inner flow following Bollimpalli, Fragile & Klużniak (2023) (see Section 2.5.2 for details).

With the equations above, the evolution of the SMBH–disc particle system can be fully specified given η, |$\dot{M}_\mathrm{BH,0}$|⁠, |$\dot{M}$|⁠, L_inner and |$\left.\frac{\mathrm{d}\boldsymbol {J}_\mathrm{\bullet }}{\mathrm{d}t}\right|_{\substack{\mathrm{LT}}}$|⁠.

For the thin disc model, all of these quantities can be obtained from global analytical models (see Fiacconi et al. 2018, for details), however, for the ADAF model only self-similar analytical solutions exist which are not valid at the inner or outer flow boundaries due to their scale-free nature. In particular η and L_inner cannot be calculated analytically in the ADAF regime. Instead the disc equations have to be evaluated numerically or, alternatively, these characteristic quantities can be extracted from GR(R)MHD simulations.

The remainder of the methodology section is structured as follows. First, we provide details on the numerical implementation, including parameter selection and the connection between the subgrid model and the hydrodynamical simulation, in Section 2.2. Subsequently, we outline the process by which the disc state is determined within our subgrid model in Section 2.3. We then elaborate on our analytical modelling and use of the latest GR(R)MHD simulations to track mass and angular momentum transfer for the ADAF and truncated disc states, respectively, which are detailed in Sections 2.4 and 2.5. Finally, we present a summary of our unified accretion disc model in Section 2.6.

2.2 Numerical implementation details

2.2.1 Conventions

In this paper, we adopt the convention whereby lower case r is a radius expressed in units of the Schwarzschild radius |$R_\mathrm{s}=\frac{2\mathrm{G}M_{\bullet }}{\mathrm{c}^{2}}$|⁠, where G is Newton’s gravitational constant, M_• is the SMBH mass, and c is the speed of light.

2.2.2 Accretion disc parameters and calibrations

It is now widely accepted that angular momentum transport in ionized accretion flows occurs via the magnetorotational instability (Balbus 1991; Balbus & Hawley 1998). MHD simulations (e.g. McKinney, Tchekhovskoy & Blandford 2012; Tchekhovskoy et al. 2012; Sa̧dowski et al. 2013; Ryan et al. 2017; Liska et al. 2018) model this self-consistently but for our hydro-only simulations we need to employ an α-like prescription for the viscous stress. For a fully ionized accretion disc, we would expect relatively high viscosities (e.g. Martin et al. 2019). Here we need to work with an effective viscosity prescription and we set α to the canonical value of α = 0.1 for the thin disc state.

In the thick disc regime, the disc scale height obeys H/R > α. As the disc is unresolved in our simulations, we fix this to a value which matches the high-resolution simulations that we use as inputs for our subgrid model, setting H/R = 0.3 (see Tchekhovskoy et al. 2012). We then calibrate α as a constrained parameter to ensure a smooth transition between the truncated and pure ADAF/ADIOS regime (see Section 2.4.1).

In the thin disc regime, we have H/R < α and we use the truncated scale height as a constrained parameter to calibrate the viscous time-scale of the disc to the time-scales obtained from the accretion model by Fiacconi et al. (2018), see also Section 2.4.2.

2.2.3 Inflow quantities

The resolved inflow quantities are calculated in analogy to Fiacconi et al. (2018). These fluxes are obtained from the hydro solver and connect our subgrid model with the hydrodynamical simulation. In brief, the inflow radius, r_in, is defined as the kernel-weighted average size of the gas cell neighbours of the SMBH within the smoothing length. The mass inflow rate, |$\dot{M}_\mathrm{in}$|⁠, is defined as the mass flux on to the SMBH particle, which is calculated via a kernel-weighted average of the mass flux provided by the neighbours within the smoothing length. Similarly L_in is the kernel-weighted average specific angular momentum. Following Bourne et al. (2023), we modify the inflow calculation for the angular momentum as this approach tends to overestimate the magnitude of L_in. Instead, we only take the direction of L_in and set the magnitude of the inflow specific angular momentum to the specific angular momentum of the subgrid accretion disc. In practice, we typically find reasonable agreement between the specific angular momentum of the subgrid disc and the specific angular momentum of gas on the scale of the subgrid disc.

For a given time-step,¹ Δt, the disc mass M_d, and disc angular momentum J_d then get updated as:

2.3 Disc mode

There are three disc modes: the thin disc (ThD), truncated disc (TrD) consisting of an inner ADIOS flow and outer thin disc, and a pure ADIOS flow.

Initially, one has to make an educated guess with regards to which state the system should be in. To make this guess effectively, we calculate an auxiliary Eddington fraction, f_{Edd, aux}, at the beginning of every time-step, based on the expected Eddington fraction if the system were in the thin disc state² which is derived in Fiacconi et al. (2018) following the thin α-disc model from Shakura & Sunyaev (1973).

with

and the Eddington rate defined as

where m_p is the proton mass, η_ref = 0.1 is the reference radiative efficiency, σ_T is the Thomson cross-section, and c is the speed of light. We then use f_{Edd, aux} to calculate the characteristic radii of the disc system and set the disc state accordingly. To determine the accretion mode of the disc system, we compare the radius of the thin disc innermost stable circular orbit (ISCO), r_ISCO, the truncation radius, r_tr, and the outer thin disc radius, r_ThD.

At intermediate Eddington ratios, the outer thin disc is expected to be truncated at a characteristic radius r_tr, transitioning into a hot accretion flow. The physical reason for this transition is not fully understood and there are several possible mechanisms such as the disc evaporation model, the turbulent diffusion model, and the large viscosity model – all of which predict that the truncation radius r_tr should decrease as f_Edd increases in agreement with observations of black hole binaries and low-luminosity AGN (see Yuan & Narayan 2014, and references therein). Also see Liu & Taam (2009) for a discussion on the application of disc evaporation models to AGN and Cho & Narayan (2022) for a generalized disc evaporation and transition model applicable to both X-ray binaries and AGN.

First, we estimate the radius that the thin disc would have in the absence of disc truncation using the expression from Fiacconi et al. (2018):

and following Yuan et al. (2018), we estimate the truncation radius as:

This expression is based on disc truncation by diffusion where the energy from the interior of the disc heats the gas above the virial temperature forming an inner hot flow and truncating the outer thin disc (see Honma 1996).

However, note that alternative truncation mechanisms give the same qualitative behaviour. In the case of the disc evaporation model, the specific evaporation rate generally decreases as a function of radius. Hence, as the Eddington ratio decreases, the truncation radius (where evaporation and gas inflow are balanced) increases. For magnetically truncated discs, we would expect similarly for the truncation radius to increase as the Eddington fraction decreases (whilst higher magnetic flux at a fixed Eddington fraction may increase the truncation radius, however, this is beyond the scope of this work, see Liska et al. 2022, for details). More holistically, global radiation MHD simulations around SMBHs have demonstrated the link between optical depth trends and disc truncation, where the optical depth increases with larger radii and higher Eddington ratios. Hence the truncation radius, the radius where the optical depth falls below one, decreases as the Eddington ratio increases (see Jiang et al. 2019).

Finally, r_ISCO is given by

where Λ(a) is a function of the SMBH spin a (see e.g. equation B2 in Fiacconi et al. 2018, for an expression for Λ(a)). Λ(a) is largest for retrograde spinning BHs (orbiting against SMBH spin), whilst highly spinning prograde BHs have the smallest ISCO.

These characteristic disc radii are shown as a function of Eddington fraction in Fig. 2. For reference, we also show the warp radius (crucial for determining the Lense–Thirring precession regime, see Section 2.5.2), the self-gravity radius beyond which the disc could fragment³ (see Fiacconi et al. 2018, for details), and the Bondi radius, as well as the typical resolution scale for the simulations presented in this work.

The ratio between the ISCO, truncation radius, and thin disc radius, then determines our disc state at each time-step, with the thin disc at high Eddington ratios, the truncated disc emerging at intermediate Eddington ratios and the pure ADIOS flow at low Eddington ratios (see Fig. 1 for an illustration of the three main disc states):

• Thin disc: If r_tr ≤ r_ISCO, the ADIOS flow component is within the ISCO of the thin disc, so we proceed with the standard thin disc model from Fiacconi et al. (2018) if the Eddington fraction is sufficiently high.

• Truncated disc: If r_ISCO < r_tr < r_ThD and |$\sqrt{\mathrm{G}M_{\bullet }R_\mathrm{tr}} \le J_\mathrm{d}/M_\mathrm{d}$|⁠, we evolve the disc system as a truncated disc.

• ADIOS flow: If r_tr ≥ r_ThD (i.e. the size of the ADIOS flow is larger than the size of the thin disc) or the gas does not have sufficient angular momentum to circularize at the truncation radius (⁠|$\sqrt{\mathrm{G}M_{\bullet }R_\mathrm{tr}} \gt J_\mathrm{d}/M_\mathrm{d}$|⁠), then we proceed with the pure ADIOS flow model. In this case, the size of the flow, r_ADIOS, may increase to larger sizes than the typical inflow scale, the Bondi radius. To avoid this complication, we limit the size of the ADIOS flow to r_Bondi, which also ensures that we do not expand the flow beyond the so-called ‘strong ADAF radius’, which provides a theoretical upper limit r_max beyond which the ADAF solution becomes invalid (see Narayan & Yi 1995).

Having developed a mechanism for determining the state of the accretion disc system, we are now ready to calculate the mass transfer, luminosity, and angular momentum transfer for each of the accretion disc configurations. For the thin disc case, we simply follow the equations from Fiacconi et al. (2018) to calculate the mass and spin evolution as well as SMBH luminosities. In Sections 2.4 and 2.5, we outline the equations governing the mass transfer and angular momentum transfer if we have a pure ADIOS flow or a truncated disc with an inner ADIOS flow.

2.4 Mass transfer

In the ADIOS model (Blandford & Begelman 1999) the accretion rate depends on the radius due to the wind loss. Therefore the definition of the Eddington fraction becomes more subtle as this will also change with radius. Here we follow Xie & Yuan (2012) and define the Eddington fraction of the ADIOS flows with respect to the mass accretion rate at the (outer) event horizon, |$r_\mathrm{H}(a)= \frac{1 + \sqrt{1-a^{2}}}{2}$|⁠, which represents the rest mass flux onto the SMBH, |$\dot{M}_\mathrm{BH,0}$|⁠. The Eddington fraction is then given by:

Once the state of the disc has been decided, we can then drop f_{Edd, aux} and use the actual Eddington fraction of the ADIOS flow to estimate the radiative efficiency, η_ADIOS. We model the radiative efficiency η_ADIOS, which is dependent on the ADIOS Eddington fraction f_Edd, using the analytical fitting function obtained by numerically solving the dynamical equations of a two-temperature plasma from Xie & Yuan (2012) corrected by the GR-R-MHD simulation results from Ryan et al. (2017) as compiled in Inayoshi et al. (2019):

The fitted values are a₀ = −0.807, a₁ = 0.27, a_n = 0 (n ≥ 2), b₀ = −1.749, b₁ = −0.267, b₂ = −0.07492, and b_n = 0 (n ≥ 3), see Inayoshi et al. (2019) for details.

Note that the radiative efficiency will likely also depend on the magnetic field threading the disc (see Liska et al. 2024), however, capturing these effects is beyond the scope of this work. In the future, we also plan to explore different radiative efficiency prescriptions based on the weak magnetic field ‘SANE’ state and saturated magnetic field ‘MAD’ state.

Using the radiative efficiency, we can then calculate the luminosity of the ADIOS flow as

and the SMBH growth rate as:

The rest mass flux onto the SMBH, |$\dot{M}_\mathrm{BH,0}$|⁠, will depend on the disc model. In particular, we obtain different rest mass fluxes for the inner ADIOS flow of the truncated disc model and the ‘pure’ ADIOS state. This is mainly due to different accretion rates at the outer edge of the ADIOS flow (in the truncated disc case this is set by the feeding rate at the inner edge of the thin disc component) and different wind mass-losses (in the truncated disc case these tend to be lower due to the smaller size of the ADIOS flow). In Sections 2.4.1 and 2.4.2, we outline how |$\dot{M}_\mathrm{BH,0}$| is calculated for the pure ADIOS flow model and the truncated disc model, respectively. Note that in this paper, we include the possibility of wind loss according to the ADIOS model, however, we do not inject this wind mass.⁴ This makes the set-up well-defined and allows us to better isolate the effect of the different disc states as live wind injection will complicate the thermodynamics. Coupling our unified accretion model with wind and jet feedback is beyond the scope of this paper and will be addressed in future work.

2.4.1 Pure ADIOS flow model

To estimate the mass accretion rate through the pure ADIOS flow, we need to estimate the viscous time-scale in the advection-dominated regime, integrating the inverse of the radial velocity from the inner to the outer edge of the ADIOS flow. As we limit the size of the accretion flow to the Bondi radius r_Bondi, we obtain r_ADIOS as:

We then estimate the mass flow rate at r_ADIOS as:

where τ_{visc, ADIOS} is the viscous time-scale of the ADIOS flow. Since the disc mass M_d is updated as |$M_\mathrm{d} \rightarrow M_\mathrm{d} + \dot{M}_\mathrm{in} \Delta t$|⁠, this mass accretion rate at r_ADIOS is simply the inflow rate rescaled by the viscous time-scale and corrected for disc draining.

The viscous time-scale may be estimated from the radial velocity of the ADIOS flow v_{r, ADIOS} as:

We use the radial velocity formula from Narayan & Yi (1995) in the set-ups with and without wind loss, since v_r is not (strongly) dependent on density and therefore not very sensitive to the wind (Yuan, Bu & Wu 2012), to obtain:

where α_ADIOS is the viscosity parameter for the ADIOS flow. Note that there are significant uncertainties in the value of this parameter (see e.g. McKinney et al. 2012), which depends on the accumulated magnetic flux, disc thickness, and SMBH spin. Therefore, we calibrate this value to achieve a smooth transition in accretion rates between the truncated disc regime and the pure ADIOS regime, choosing α_ADIOS such that the mass flow rate through the pure ADIOS flow for |$r_\mathrm{tr} \xrightarrow {} r_\mathrm{ThD}$| tends to the mass flow rate through the outer thin disc for the same limit (see Section 2.4.2). Following Blandford & Begelman (1999), we can then use the accretion rate at r_ADIOS to calculate the mass accretion rate at the event horizon (r_H):

where s characterizes the importance of the outflow. This parameter cannot exceed s = 1 for energetic reasons and s = 0 corresponds to no outflow. For our ADIOS-model-based simulations, we set s = 0.4 for consistency with Xie & Yuan (2012) who assume this value in their numerical determination of the radiative efficiency of the outflow. We also run a reference simulation without wind loss setting s = 0.0.

2.4.2 Truncated disc model

In the truncated disc state, the inner hot flow is being fed by the outer thin disc, which limits the SMBH mass growth rate. The variability of the inner hot flow could, in principle, differ from the outer thin disc. Thick discs tend to exhibit higher variability (e.g. Lalakos et al. 2024; Liska et al. 2024), which may be explained by a variable α parameter (e.g. Turner & Reynolds 2023), and have shorter viscous time-scales. In our model, we set the mass flow rate through the inner hot flow (modulo wind mass-loss) equal to the feeding rate from the outer thin disc, for simplicity. This will not affect the overall SMBH mass growth rate but may impact electromagnetic counterpart predictions and we discuss this in more detail in Section 4.2.

The derivation of the thin disc mass accretion rate in Fiacconi et al. (2018) relies on the inner disc radius being much smaller than the outer disc radius. This is generally true when the inner radius is at the ISCO, however, this assumption may no longer be valid for the truncated disc model when the truncation radius can reach a significant fraction of the thin disc radius. For this reason, we only use the formula for |$\dot{M}_\mathrm{ThD}$| from Fiacconi et al. (2018) if r_tr/r_ThD < 10⁻². In this case, the thin disc accretion rate is given by:

If r_tr/r_ThD ≥ 10⁻², we estimate the thin disc accretion rate via the viscous time-scale τ_visc:

Here we evaluate the viscous time-scale based on the radial velocity for the thin α-disc (also see Pringle 1981)

where ν = αc_sH is the kinematic viscosity. Integrating this expression between the truncation radius and the outer thin disc radius yields:

The thin disc mass M_ThD, adjusted for the truncation radius, is given by (Fiacconi et al. 2018):

The derivation of the (unresolved) scale height for the thin disc model is significantly modified in the presence of an inner hot flow (e.g. Różańska et al. 1999) and will likely not apply for the heavily truncated discs where we employ the viscous time-scale formula – therefore we treat H/R as a free parameter in this regime (though we enforce H/R < α). We then calibrate the scale height such that |$\frac{M_\mathrm{ThD}}{\tau _\mathrm{visc}}$| matches |$\dot{M}_\mathrm{ThD}$| for r_tr = 0.01r_ThD.

To calculate the feeding rate at the truncation radius, the radiative losses in the thin disc need to be estimated. Note that for the ‘pure’ thin disc model the radiative efficiency is given by:

as r_ISCO defines an effective surface for the maximum feasible energy extraction from infalling particles. In the thin disc case, we can only extract energy from the thin disc component down to r_tr, so the expression for η_TrD then becomes:

The ADIOS flow feeding rate at the truncation radius is therefore set by the thin disc mass accretion rate |$\dot{M}_\mathrm{ThD}$| and the radiative efficiency η_TrD. In analogy to the previous section, we can then use this feeding rate to estimate the accretion rate at the event horizon (ADIOS model):

where we assume s = 0.4 for the simulations with wind loss and set s = 0.0 for the comparison simulations without wind loss.

Fig. 3 shows the Eddington fraction of the unified accretion disc model f_Edd as a function of the Eddington fraction in the thin disc model f_{Edd, ThD} for |$M_{\bullet } = 10^{6} \ \mathrm{{\rm M}_{\odot }}$|⁠. The dark purple and pink lines correspond to the pure ADIOS flow model without and with wind loss, respectively. The steeper slope in this regime is introduced by the dependence of the truncation radius on the thin disc Eddington fraction to the second power, which leads to a rapid increase in the viscous time-scale as f_{Edd, ThD} decreases. The gradient shifts at f_{Edd, ThD} ≲ 10⁻⁵, where the ADIOS radius becomes limited by the Bondi radius. The red and orange lines show the truncated disc Eddington fractions without and with wind loss, respectively. Without wind loss, this closely corresponds to the thin disc rate (indicated by the grey dashed line) as the inner hot flow rate is set by the outer feeding rate. Note the slight gradient change at f_{Edd, ThD} ≲ 10⁻³ as the thin disc accretion rate prescription switches to being estimated by the viscous time-scale. Finally, for f_{Edd, ThD} ≳ 10⁻², the system switches to the pure thin disc state (blue line).

2.5 Angular momentum transfer

In this Section, we describe how to model the angular momentum transfer for each disc state both due to accretion (see Section 2.5.1) and due to Lense–Thirring precession (see Section 2.5.2). Note that, as for the mass transfer, if we are in the thin disc regime, we evolve the angular momenta of the disc and the SMBH according to the equations from Fiacconi et al. (2018). Below we outline our scheme for the pure ADIOS flow and truncated disc cases.

2.5.1 Accretion

For both the pure ADIOS flow and the truncated disc model with an inner ADIOS flow, we model the evolution of the SMBH angular momentum, J_•, due to accretion as:

where j_• and j_d are the angular momentum unit vectors (‘versors’) of the SMBH and the disc, respectively, and L_ADIOS is the specific angular momentum of the ADIOS flow. The latter is evaluated at the event horizon because for thick discs the stress is non-zero down to the horizon (see e.g. Sa̧dowski et al. 2016). To evaluate this expression, we need to know the value of L_ADIOS, which we would expect to be below the thin disc value due to additional (gas and magnetic) pressure support. We obtain L_ADIOS from the analysis by Lowell et al. (2024) where they estimate the contributions to the spin-up function from both electromagnetic (jet) and hydrodynamic (disc) components using the GRMHD simulations of an ADIOS flow in the MAD state by Tchekhovskoy et al. (2012). They find that there is no significant dependence of the specific angular momentum L_ADIOS on the spin parameter and obtain an approximately constant value L_ADIOS ∼ 0.86 in geometrical units for the radiatively inefficient regime.

2.5.2 Lense–Thirring precession

Next, we focus on the angular momentum transfer between the SMBH and the disc due to the Lense–Thirring effect (Lense & Thirring 1918). Lense–Thirring precession is a frame-dragging effect which manifests itself as nodal precession of orbits whose angular momentum is misaligned with the SMBH spin. In the weak-field limit, the precession frequency, ω_LT, is given by (Lense & Thirring 1918):

Accretion discs are warped by the differential nature of the Lense–Thirring precession frequency. How these warps propagate depends on the type of accretion disc.

In the diffusive regime (α > H/R, thin disc case), the warps are communicated by viscosity resulting into the so-called Bardeen-Petterson configuration (Bardeen & Petterson 1975), where the disc is aligned with the SMBH equatorial plane at small radii (r < r_warp) and retains its initial tilt at large radii (r > r_warp) with a smooth warp in between at r_warp which is given by (see Fiacconi et al. 2018, for derivation):

where ξ is a parameter |$\mathcal {O}(1)$| which is related to the ratio of the radial viscosity ν₁ and vertical viscosity ν₂ as ν₁/ν₂ = ξ/(2α²) and can be determined numerically (Papaloizou & Pringle 1983; Lodato & Pringle 2007).

In the wave-like regime (α < H/R, ADIOS case), warps are communicated by pressure waves. These so-called bending waves result in a smooth warp far from the SMBH, whilst radial tilt oscillations are obtained close to the BH, as the wavelength of these pressure waves is a strong function of radius (see Ingram & Motta 2019, for details).

For our pure ADIOS flow model, the system is in the wave-like regime. The whole flow, which is strongly coupled by the pressure waves (the local sound crossing time-scale is much shorter than the precession time-scale), then remains misaligned and precesses as a solid body (see Ingram, Done & Fragile 2009; Ingram & Done 2011; Ingram & Motta 2019).

For our truncated disc model, the extent of the truncation radius determines whether we are in the diffusive or the wave-like regime.

If r_tr < r_warp, the thin disc feeds the ADIOS flow component aligned material, thereby twisting the ADIOS flow into alignment. Also see Liska et al. (2018) who find in their numerical simulations that the angular momentum of the inner hot flow is negligible so that the thin disc can easily affect its alignment.

If r_tr ≥ r_warp, the thin disc feeds misaligned material to the inner flow, so there is no alignment in the inner region and the inner ADIOS flow precesses as a solid body just as in the pure ADIOS flow case.

Following Ingram & Motta (2019), we can calculate the angular frequency of solid body precession, ω_prec, as:

where |$\mathcal {L}(R) = \Sigma R^{2} \Omega _{\phi }$| is the angular momentum per unit area. Note that since the |$\mathcal {L}(R)$| term appears in both the numerator and the denominator, we only need to know how this term scales with R, but not the absolute value. Using the scaling relations for the ADIOS model from Yuan & Narayan (2014), we have that the angular frequency scales as Ω_ϕ ∼ R^3/2 and the volume density scales as ρ ∼ R^{−3/2 + s}. Under the assumption of a constant H/R ratio, the surface density then scales as Σ ∼ R^{−1/2 + s}. Using the weak-field formula for ω_LT (equation 31), the integrals in equation (33) can then be evaluated analytically (Fragile et al. 2007; Ingram et al. 2009).

Note that with s = 0.4, the surface density Σ is nearly constant. However, GRMHD simulations predict that for misaligned discs, the density should drop off sharply in the region dominated by radial tilt oscillations. We therefore set the inner radius to be the bending wave radius r_bw = 3.0(H/R)^−4/5a^2/5 and the outer radius is set to be equal to the extent of the ADIOS flow.

The torque due to Lense–Thirring precession is then given by:

where in the pure ADIOS flow model, we simply have J_{d, ADIOS} = J_d. In the truncated disc case, we assume that the Σ ∼ R^−3/4 scaling still holds for the thin disc component, so that we can scale the thin disc angular momentum as:

and consequently the ADIOS angular momentum is given by

In addition to disc precession, we also expect a disc alignment process – much slower than the Bardeen-Petterson alignment in the thin disc case, with the thick disc aligning on the accretion time-scale |$t_\mathrm{acc} = M_\mathrm{BH}/\dot{M}_\mathrm{BH,0}$| (see e.g. King et al. 2005; Volonteri et al. 2005, for analytical models of disc alignment) as well as the GRMHD simulations by Liska et al. (2023a), which also find a global alignment mode for the disc-corona-jet system on the accretion time-scale. The full Lense–Thirring torque due to disc precession and alignment is then given by:

In Fig. 4 we show the solid body precession time-scales (blue-shaded regions) compared to the sound crossing time-scale (orange line, see Motta et al. 2018, for derivation). From equation (34), we can see that the precession time-scales depends on ω_prec, the angle between the versors and the ratio of the ADIOS angular momentum to BH angular momentum. Indeed from this equation, we can derive that the rate of change of the azimuthal angle goes as |$\frac{\mathrm{d}\phi }{\mathrm{d}t} = \omega _\mathrm{prec} \frac{J_\mathrm{d,ADIOS}}{J_\mathrm{d}}$| when the angular momentum of the BH dominates (recovering |$\frac{\mathrm{d}\phi }{\mathrm{d}t} = \omega _\mathrm{prec}$| in the pure ADIOS regime). The |$\frac{\mathrm{d}\phi }{\mathrm{d}t} = \omega _\mathrm{prec}$| case provides the lower limit on the precession time-scales in Fig. 4 and the upper limit (used in our implementation) is given by the truncated-disc-corrected precession time-scale, in all cases assuming J_d/J_BH = 0.9, matching the typical angular momentum ratios in our binary SMBH simulations. The light-blue shaded region corresponds to the precession time-scales without considering additional torques between the inner hot flow and outer truncated thin disc, whilst the dark-blue shaded region shows the time-scales corrected by 95 per cent, following Bollimpalli et al. (2023) who found that solid body precession may be slowed down by up to 95 per cent for truncated discs due to additional torques. In our model, we include this slowed down precession as an optional module to test the impact of reduced precession rates in this regime.

Finally, we note that in the MAD regime, jets can force the inner part of the ADIOS flow to align rapidly with the SMBH spin in GRMHD simulations (e.g. McKinney et al. 2013). However, Liska et al. (2018) investigated a range of initial conditions for GRMHD simulations of tilted black hole discs and found that the importance of this effect crucially depends on the magnetic field configuration in the initial conditions. As we do not model the jet component, we do not include electromagnetic alignment due to jets here. However, we will consider this effect in future work when coupling our accretion disc model with an AGN jet subgrid prescription.

2.6 Summary of the model

As discussed, the evolution of the SMBH and its accretion disc within our unified model can be fully specified given the radiative efficiency, η, rest mass flow rate onto the SMBH, |$\dot{M}_\mathrm{BH,0}$|⁠, rest mass flow rate through the disc, |$\dot{M}$|⁠, specific angular momentum at the inner boundary, L_inner, and Lense–Thirring torque, |$\left.\frac{\mathrm{d}\boldsymbol {J}_\mathrm{\bullet }}{\mathrm{d}t}\right|_{\substack{\mathrm{LT}}}$|⁠. These quantities are listed for each state of the disc system in Table 1. We also list the relevant references for all of the equations.

Note that the derivations of the thin disc equations can be found in Fiacconi et al. (2018), in particular the derivation of the Bardeen-Petterson torque for Lense–Thirring precession in the thin disc regime as listed in Table 1 (here J_• is the magnitude of the SMBH angular momentum and τ_align is the time-scale for the torque to modify the SMBH angular momentum). Note that if the warp radius exceeds the thin disc radius, instantaneous alignment of SMBH and disc angular momenta is assumed for the thin disc case, as listed in Table 1.

3 SIMULATIONS

To validate our new accretion disc implementation, we perform idealized simulations of a single or binary SMBH embedded in a gaseous disc, representing the inner region of a galaxy (see Fig. 5).

These test cases allow us to assess the sensitivity of our accretion disc subgrid model to different parameter choices and modelling assumptions, such as disc wind loss or reduced solid body precession, without being limited by the computational expense of full galaxy formation simulations. Ultimately, our accretion model has been developed with galaxy simulations in mind, however, exploring the accretion disc subgrid model in these set-ups is beyond the scope of this methodology-focused paper.

The single BH embedded in a circumnuclear disc (CND) represents a baseline case so that we can assess the impact of different modelling choices in a relatively steady environment. The binary SMBHs embedded in the circumbinary disc (CBD) allow us to study the behaviour of our new accretion disc subgrid model for strongly fluctuating gas inflow patterns, especially for the misaligned binary set-up. Furthermore, the binary simulations demonstrate a crucial science case for these types of SMBH accretion models, in particular with regards to electromagnetic counterpart and gravitational recoil predictions.

In this Section, we outline the different physical configurations and model variations explored in this paper. The set-up of the gaseous disc is based on the initial conditions (ICs) presented in Bourne et al. (2023) and for completeness we provide a brief summary in Section 3.1. We then outline how these ICs are evolved for the single SMBH set-up (see Section 3.2) and the binary SMBH set-up (see Section 3.3), including the subgrid model variations explored in both cases.

3.1 Gaseous disc set-up

The gaseous disc for both the single and binary set-ups extends from R_in = 4 pc to R_out = 14 pc and follows a surface density profile of

where α = 2 and the normalization Σ₀ is set so that the total mass of the gaseous disc constitutes 10 per cent of the single/binary SMBH mass. As discussed in Bourne et al. (2023), this set-up represents an idealized version of the inner region of a post-merger galaxy, which offers ideal conditions for disc formation via the circularization of infalling gas clouds.

The vertical structure is set up so that the aspect ratio of the gaseous disc, H/R, is approximately constant and the disc is in a stable configuration (initial Toomre parameter set to Q = 1.5).

We then evolve this system as an ideal gas with adiabatic index γ = 5/3 and allow the gas to cool following a β-cooling prescription, where the local cooling time-scale T_cool is proportional to the orbital time-scale by a factor β = 10. This ensures that the disc can settle into a marginally stable configuration (see Bourne et al. 2023, for details).

The target gas mass resolution is set to |$m_\mathrm{target}=0.2~\mathrm{{\rm M}_{\odot }}$|⁠. To increase the spatial resolution in the gaseous disc cavity, and allow for accurate measurements of the mass and angular momentum inflow rates onto the subgrid accretion disc, we apply the super-Lagrangian refinement technique (Curtis & Sijacki 2015) to the gas cells in the cavity (within R_ref = 3 pc). To this end, we define an additional (de-)refinement criterion so that the maximum allowed cell size |$R^\mathrm{cell}_\mathrm{max}$| decreases linearly within the cavity from |$R^\mathrm{cell}_\mathrm{max}=0.2$| pc at R_ref (typical size of the mesh cell just inside the cavity) to |$R^\mathrm{cell}_\mathrm{min}=0.01$| pc at the centre (chosen to match the gravitational softening of the BH). The cell sizes are allowed to be within a factor of two of these target radii and, to avoid excessive refinement, we set a minimum cell mass of m_min = 10⁻⁵m_target. Note that the super-Lagrangian refinement technique is not merely advantageous for the subgrid model, but also allows us to resolve important features in the low-density cavity (which is intrinsically difficult for Lagrangian codes); in particular we can resolve the formation of streams and mini discs around the binary SMBHs, and follow accurately how they torque the binary.

3.2 Single SMBH simulations

3.2.1 Initial conditions and relaxation

For all single SMBH simulation set-ups, we place a SMBH at the centre of the CND, with mass |$M_\mathrm{\bullet } = 2 \times 10^6 \ \mathrm{{\rm M}_{\odot }}$|⁠. Following a relaxation time of ∼20 Myr, we then activate the super-Lagrangian refinement and further relax the system for ∼6 Myr. Note that we keep the overall relaxation time deliberately as short as possible. This allows for the disappearance of transient features, whilst testing our subgrid model before the cavity has been refilled with gas due to the absence of binary torques, as we are mainly interested in the low-accretion-rate regime [see also the counter-rotating binaries in Bourne et al. (2023)]. Appendix  A summarizes the main features of the single SMBH relaxation procedure.

3.2.2 Subgrid model variations

For the subgrid model, we need to specify the initial masses and angular momenta of the accretion disc and SMBHs.

The SMBH spin depends on the system’s previous evolution with different accretion flows and feedback physics leading to either spin-up or spin-down (see Narayan et al. 2022; Lowell et al. 2024, for recent compilations). Following Bourne et al. (2023), we assume that both SMBHs are initially highly spinning with a spin value of a = 0.9, consistent with current observational constraints on |$\sim 10^{6}~\mathrm{{\rm M}_{\odot }}$| SMBHs (Reynolds 2019). The orientation is chosen randomly, however, to be able to make direct comparisons between different physics runs, the initially randomly chosen orientation for the first set-up, i.e. θ = 62°, is kept the same for all accretion disc physics variations. The initial angular momentum versor of the subgrid accretion disc is aligned with the angular momentum of the surrounding gas inflows, i.e. with the z-axis.

For the subgrid accretion disc evolution, we test the original thin disc model (single_thindisc) as well as the unified model with and without wind loss (single_uniadios_winds and single_uniadios). For these tests, we choose the initial disc mass as |$M_\mathrm{d} = 10^{-3} M_\mathrm{\bullet }=2000 \ \mathrm{{\rm M}_{\odot }}$| and f_{Edd, initial, TD} = 0.003, which leads to good agreement between the external inflow rate and the Eddington fraction for the ADIOS model, see Fig. 7.

Furthermore, we also test two variations of the initial disc mass set-up, rescaling the initial Eddington ratio according to equation (8). In the first case, we employ a significantly higher disc mass, |$M_\mathrm{d}=7500 \ \mathrm{{\rm M}_{\odot }}$|⁠, and the unified model without wind loss (single_uniadios_massive), and in the second case, we employ a significantly lower disc mass, |$M_\mathrm{d}=500 \ \mathrm{{\rm M}_{\odot }}$|⁠, and the unified model both without and with wind loss (single_uniadios_light and single_uniadios_winds_light). All of these runs and their parameters are listed in Table 2.

The different disc mass choices here represent different previous evolution histories of the system and allow us to test the behaviour of our model in a low-inflow environment.

3.3 Binary set-up

3.3.1 Initial conditions and relaxation

The total binary mass is set to |$M_{\rm bin} = 2 \times 10^{6}~\mathrm{{\rm M}_{\odot }}$| and the binary is initially placed on a Keplerian orbit with semimajor axis a_orbit = 2 pc, resulting in an orbital period of T_bin = 0.187 Myr.

For the binary orbits, we consider two configurations. First, we simulate a simple symmetric set-up with an equal mass ratio, q = 1, zero eccentricity and with the binary orbit aligned with the plane of the gaseous disc (labelled as bin_q1i0). Secondly, we simulate binaries with a mass ratio q = 3, zero eccentricity, and with the binary orbit inclined by i = 45° with respect to the CBD (labelled as bin_q3i45). These set-ups are equivalent to q01e00i00 and q03e00i45mod in Bourne et al. (2023), respectively.

The initial conditions are evolved for 50 binary orbits corresponding to about 1.5T_cool at the outer edge of the CBD without any additional refinement and without any accretion onto the SMBH–accretion disc particle. Following the cooling relaxation period, the set-up is further relaxed with the super-Lagrangian refinement for another 50 binary orbits, so that the total relaxation period amounts to approximately three cooling time-scales, allowing the CBD to reach a nearly steady state (see Bourne et al. 2023). The resulting snapshot is used as initial conditions to be further evolved for 500 binary orbits, testing the unified accretion disc model (varying the assumptions about wind loss and LT precession) as well as reference runs with the thin disc model from Fiacconi et al. (2018), as described in Section 3.3.2.

3.3.2 Subgrid model variations

For the subgrid modelling, we again need to specify the initial masses and angular momenta of the accretion disc and SMBHs. Following the single BH set-up and Bourne et al. (2023), we assume that both SMBHs are initially highly spinning with a spin value of a = 0.9. The orientation is chosen randomly, however, to be able to make direct comparisons between different physics runs, the initially randomly chosen orientation for the first set-up, i.e. θ_primary = 31° and θ_secondary = 62°, is kept the same for all accretion disc physics variations. Note, however, that this initial random orientation is different from the initial SMBH versor orientations in Bourne et al. (2023), so that our thin disc model runs display some small differences compared to the Bourne et al. (2023) set-ups.

For the subgrid accretion disc, the initial mass is set to M_d = 10⁻³M_•. We align the initial angular momentum versor with the angular momentum of the surrounding mini disc and we choose the ratio between the magnitudes of the disc and SMBH angular momenta such that the initial Eddington ratio, as set by equation (8), approximately matches the initial inflow rate from the mini disc (with f_{Edd, initial1, 2} = 10⁻³ for bin_q1i0 and f_{Edd, initial1} = 2 × 10⁻³ and f_{Edd, initial2} = 6 × 10⁻³ for bin_q3i45).

For both binary configurations, we run five different simulations. First, we repeat the thin-disc only set-up from Bourne et al. (2023), though as noted above, we have a different initial orientation of the SMBH versors. Secondly, we rerun both binary configurations using the ADIOS model without wind loss (uniadios) with efficient solid body precession for truncated discs following Ingram et al. (2009) as well as reduced solid body precession (slowprec) following Bollimpalli et al. (2023). Finally, we repeat the same two precession cases, where we also account for wind loss in the hot flow (uniadios_winds). All of these simulations are listed in Table 3.

4 RESULTS

The unified accretion disc model determines the evolution of the SMBH accretion rates, luminosities, and spins across different radiative regimes. In this section, we analyse all of these aspects in detail for both the single SMBH set-ups, which allow us to test the impact of our subgrid model choices in a simpler setting without the binary torques, as well as the binary SMBHs where we obtain more complex dependencies between the variations in the resolved inflows and the response of the subgrid model, in particular for the spin alignment.

Fig. 5 shows gas surface density maps of our three main simulation set-ups: the single SMBH simulation (single, left panel), the equal-mass, aligned binaries (bin_q1i0, middle panel), and the unequal-mass, misaligned binaries (bin_q3i45, right panel). We show the single BH simulation at t = 50 Myr to demonstrate how the cavity is filled with relatively low-density gas in the single SMBH simulation. For the binary simulations, on the other hand, streams and mini discs form in the cavity that feed gas and angular momentum to the SMBHs, here shown at the beginning of the simulation before the CBD is torqued by the misaligned binary (also see Bourne et al. 2023).

First, we consider the disc state transitions throughout the simulations in Section 4.1 and then go on to analyse the gas accretion and luminosities in Section 4.2. Finally, we investigate the SMBH spin and disc angular momentum evolution due to gas inflows and Lense–Thirring precession in Section 4.3.

4.1 Disc state transitions

We begin our analysis of the simulation suite by considering the different characteristic disc radii which determine the disc state transitions. Fig. 6 shows the thin disc radius r_ThD (dashed lines), the truncation radius r_tr (solid lines), and the warp radius r_warp (dotted lines) as a function of time. The grey-shaded region indicates the area inside the ISCO. The panel titles indicate the initial conditions and the colour-coding of the lines indicates the subgrid model employed as listed in the legend.

First, we consider the single SMBH set-up. With the standard initial disc mass (⁠|$M_\mathrm{d} = 10^{-3} M_{\bullet } = 2000 \ \mathrm{{\rm M}_{\odot }}$|⁠), the disc is in the truncated disc state for the whole simulation (r_ISCO < r_tr < r_ThD), with the truncation radius larger than the warp radius. Hence the outer truncated thin disc is fully misaligned and feeds misaligned material to the inner hot flow leading to solid body precession (see Section 4.3.1 for details on the precession behaviour in the single SMBH case). Note that the models with and without wind loss have the same mass flow rate through the disc for a given SMBH mass and spin (as the wind loss is at the expense of the SMBH growth and is assumed to occur at the ISCO). Given the insignificant changes in SMBH properties for the simulations considered here, these models hence have the same disc mass flow rates and therefore also virtually identical disc radii.

With the massive disc initialization (uniadios_massive, |$M_\mathrm{d} = 7500 \ \mathrm{{\rm M}_{\odot }}$|⁠), however, the disc spends a short initial period of approximately 1 Myr in the thin disc state with the truncation radius inside the ISCO. Afterwards, the disc is rapidly depleted and changes to a truncated disc with steadily increasing truncation radius. At t ∼ 10 Myr, the truncation radius starts exceeding the warp radius, so that the torque regime changes from Bardeen–Petterson alignment to solid body precession.⁵

Next, we consider the disc radii for the equal-mass, aligned binary (bin_q1i0, second panel). As the primary and secondary SMBH findings closely follow one another, we just show the disc radii evolution of the secondary here and comment on any (small) differences in the text. As with the standard single SMBH configuration, this simulation is in the truncated disc state for the whole duration of the simulation with r_ISCO < r_tr < r_ThD. Initially the truncation radius is smaller than the warp radius, so that the system is in the Bardeen–Petterson configuration. As the truncation radius increases the uniadios set-up (purple line) enters the solid body precession regime for the remainder of the simulation. The uniadios_winds set-up (pink line), however, briefly dips back into the Bardeen–Petterson regime as the truncation radius decreases due to the slower SMBH mass growth rate [which increases f_{Edd, ThD}, see equation (8)]. Note that this temporary decrease does not occur for the primary SMBH due to small stochastic differences in their evolution.

For the unequal-mass ratio, misaligned binaries, we show the primary and secondary SMBH in separate panels. Due to the asymmetry of the set-up, there are significant quantitative differences in the disc radii evolution, although the qualitative behaviour is very similar. Both SMBHs are in the truncated disc regime with r_ISCO < r_tr < r_ThD for the whole simulation. For the primary SMBH, the system is initially in the Bardeen–Petterson configuration with r_tr < r_warp and then switches from the solid-body precession regime back to Bardeen–Petterson and finally enters the solid-body regime again. This applies to both the uniadios and uniadios_winds model, although the uniadios_winds model spends longer in the Bardeen–Petterson regime due to the smaller truncation radius. Similarly, the secondary SMBH switches between the different precession and alignment regimes, however with the uniadios_winds model, there is only one solid body precession phase at late times, with the early precession phase being suppressed by the slower growth in r_tr. Note that the truncation radius of the secondary increases significantly towards late times (over one magnitude higher than the primary’s truncation radius) as with our binary set-up, we have preferential accretion onto the primary (see discussion in Bourne et al. 2023).

Finally, we note that for the (primary) BHs in the binary simulations, the truncation radius closely tracks the warp radius. This is because our binary SMBHs are in the Eddington fraction regime of f_Edd ≳ 10⁻³, which corresponds to the accretion regime where these two radii intersect for highly spinning BHs (see Fig. 2). The primary BHs closely maintain this accretion rate throughout the simulation as the external inflow rates are similar to the mass flow rate through the accretion disc so that the latter does not change its mass significantly. Whilst the single BH simulations also initially have r_tr ∼ r_warp, these two radii rapidly diverge as the accretion rate decreases due to the significantly lower external inflows, with r_warp having a much shallower dependence on f_Edd than r_tr.

Having established the different disc regimes that govern the mass and angular momentum evolution of the SMBHs, we move on to a more detailed analysis of the gas accretion rates and luminosities in Section 4.2 and of the disc angular momentum and SMBH spin evolution in Section 4.3.

4.2 Gas accretion rates and luminosities

4.2.1 Single SMBH simulations

We begin our analysis of the gas mass accretion rate by considering the single SMBH simulations. Fig. 7 shows the time evolution of the gas mass accretion rate at the ISCO (top panel), bolometric luminosity (middle panel), and SMBH spin (bottom panel). The colour-coding indicates the subgrid model employed, as listed in the legend. For comparison, we also show the gas mass inflow rates onto the accretion disc in light green as well as the mean inflow rate in bins of 0.1 Myr in dark green.

We note that due to the low cooling rate, the cavity only fills in relatively slowly (also see Appendix  A), so the initial disc mass plays an important role, leading to significant differences between the different disc initializations (see Table 2).

We first consider the reference thin-disc only run (solid brown line). Here the mass accretion rate and luminosity are relatively steady at |$\sim 10^{-5} \ \mathrm{{\rm M}_{\odot }}\, \mathrm{yr}^{-1}$| and |$\sim 10^{41} \ \mathrm{erg \, s^{-1}}$|⁠, albeit slightly declining as the disc mass is depleted at a much faster rate than it is being replenished by the inflows, i.e. |$\dot{M}_\mathrm{inflow} \lesssim 10^{-7} \ \mathrm{{\rm M}_{\odot }}\, \mathrm{yr}^{-1}$|⁠. The SMBH spin a is very slowly increasing as the equilibrium spin for the co-rotating thin disc tends towards maximally spinning solution, i.e. a_eq  ∼ 1.

With the unified accretion disc model, the mass accretion rate is virtually identical for the uniadios set-up (dash-dotted dark pink line) as the disc is in the truncated state (see Section 4.1) so that the mass flow rate is set by the outer thin disc. Note that this modelling choice results in the mass flow rates through the inner hot flow being significantly smoothed out. In principle, this boundary condition only needs to be true in a time averaged sense, and due to the significantly lower viscous time-scale in the ADAF regime, the inner hot flow may display significant variability, as seen in GRMHD simulations (e.g. Lalakos et al. 2024). Furthermore, the thin disc may display significantly more variability than is captured in our prescription, for example due to opacity changes (e.g. Jiang & Blaes 2020). This has important implications for AGN bursts and multimessenger signatures, as discussed in more detail in Section 4.2.2.

The uniadios_winds model (dash-dotted light pink line) also takes into account wind loss in the inner hot flow, so that the accretion rate at ISCO is significantly suppressed. Both the uniadios and the uniadios_winds model have suppressed luminosities compared to the thindisc case due to the much lower radiative efficiency as well as wind loss for the ADIOS model.

Considering the SMBH spin evolution, we find that the truncated accretion disc model behaves in the opposite way from the thindisc case, with the spin steadily decreasing. Again this decrease in the spin amplitude is extremely slow due to the low accretion rate (especially for the ADIOS case).

With the massive disc initialization (uniadios_massive, dash-dotted yellow line), the disc is in the thin-disc state for ∼1 Myr and then switches back to the truncated disc mode as the accretion rate rapidly decreases (also see Section 4.1). This also leads to a very bright initial luminosity peak of |$L_\mathrm{bol} \sim 10^{44} \ \mathrm{erg \, s^{-1}}$|⁠. The spin rapidly increases during the thin disc phase but then quickly decays again in the truncated disc phase.

Finally, we also tested a light disc set-up (not shown here), which is characterized by extremely low accretion rates and luminosities (⁠|$\dot{M} \sim 10^{-11} \ \mathrm{{\rm M}_{\odot }}\, \mathrm{yr^{-1}}$| and |$L_\mathrm{bol} \sim 10^{33} \ \mathrm{erg \, s^{-1}}$|⁠) so that, combined with the low inflow rates, there is virtually no change in the disc and SMBH properties.

4.2.2 Binary SMBH simulations

Having established the basic impact of the subgrid accretion disc models on the gas accretion rates and luminosities with the single SMBH simulations, we turn to the more complex and astrophysically interesting binary simulations. Fig. 8 shows the gas mass flow rates (top row), corresponding Eddington fractions (middle row), bolometric luminosities (bottom row). In the first column, we show the results from the equal-mass, aligned binary case, bin_q1i0, again just plotting the rates of the secondary SMBH which very closely follows the evolution of the primary SMBH due to the symmetric set-up. For the misaligned case, we separately plot the primary SMBH in the second column and the secondary SMBH in the third column.

We first focus on the gas mass flow rates in the top row. The blue lines indicate the mass flow rates from the CBD through the cavity calculated at r = 2a_orbit = 4 pc, with the solid blue lines indicating net inflows for a given timebin and the dotted blue line showing net outflows, both averaged over ∼5 binary orbits. We correct these CBD rates by a factor of 0.5, to account for the average flow rate towards the primary/secondary SMBH (though note that preferential accretion means that for the unequal-mass binaries, the SMBHs will not receive an equal share of the CBD inflows). The CBD mass flow rates fluctuate between inflows and outflows, however, the net mass flow rate over the whole simulation time is inflowing. Furthermore, we show the gas mass flow rates from the mini discs onto the accretion disc as a light green shaded region with the mean inflow rate, in bins of ∼5 binary orbits, plotted in dark green. We also show the mass flow rates through the subgrid disc for the different accretion disc models explored, including the thin-disc only reference run (thindisc, solid dark purple lines) as well as the unified accretion disc model without wind loss (uniadios, dash-dotted light purple lines) and with wind loss (uniadios_winds, dash-dotted pink lines).

The mean feeding rate from the CBD and the gas mass flow rates through the mini disc onto the accretion disc are similar overall, in particular for the aligned set-up and the primary SMBH in the misaligned set-up. For the secondary SMBH, the CBD feeding rate somewhat exceeds the mini disc rate. The latter significantly decreases at t ∼ 30 Myr and we have preferential accretion onto the primary BH (also see Bourne et al. 2023).

As noted in Section 4.2.1, the subgrid accretion disc model smooths out the inflow rates due to the relatively long viscous time-scale of the thin disc, which also sets the accretion disc time-scale in the truncated disc state (see Hogg & Reynolds 2018). In principle, the inner hot flow could experience significant variability due to disc turbulence and much shorter viscous time-scale. Capturing this variability is beyond the scope of this work, however, we note that there may be intermittent accretion bursts that our model does not capture.

Comparing the mass flow rates between the different subgrid models, we find similarly to the single SMBH case that the thindisc and uniadios models yield very similar mass flow rates. The slight positive offset of the uniadios rates is due to the extremely low radiative efficiency in the ADIOS regime so that the disc mass is being drained more slowly leading to somewhat higher accretion rates at late times. However, this is vastly outweighed by the wind loss in the uniadios_winds set-up, which has much lower accretion rates. Note that this represents somewhat of a lower limit as we do not have feedback in our simulations so the lost wind mass is not re-injected. If we included this effect, we may have a fountain effect resupplying mass to the discs, especially in the inner region (e.g. Sa̧dowski et al. 2013), though the winds may also act to reduce the mini disc mass flow rate.

The Eddington fractions of the inflow rates from the mini discs are generally ≲ 10⁻³, and with our set-ups approximately matching these initial inflow rates, the mini discs therefore set the truncated disc state.

For the luminosities, we separately plot the contributions from the outer truncated thin disc (dashed lines) and inner hot flow (dotted lines). This demonstrates that the inner hot flow significantly dominates the disc luminosity in the truncated state, which means that the spectrum of the binary would also be shifted towards harder X-rays with the truncated disc model (Giustini & Proga 2019). Energy extraction in the outer truncated thin disc is extremely inefficient for intermediate Eddington ratios as the gravitational potential energy released up to the truncation radius is only minimal, with most of the energy release in the full thin disc mode occurring at ISCO. Indeed even the inner hot flow in the uniadios model is one order of magnitude fainter than the thindisc model. For the uniadios_winds model, the low radiative efficiency and wind loss combined lead to a dramatic decrease in the inner hot flow luminosity.

This has important implications for multimessenger signatures, and in particular electromagnetic counterpart predictions. To investigate this further, we show the X-ray and radio luminosities of the binaries in Fig. 9.

To obtain the X-ray luminosities for the thin discs, we assume the bolometric correction from Shen et al. (2020), appropriate for a quasar template. For the ADIOS X-ray luminosities, we use the Eddington-ratio dependent bolometric correction from Vasudevan & Fabian (2007). For comparison, we also show the detectable luminosities for upcoming X-ray missions, with NewAthena plotted as brown-shaded regions and the proposed NASA APEX missions as grey-shaded regions. In both cases, the light-shaded regions indicate the detectable luminosities at z = 0.25, whilst the dark-shaded regions indicate detectable luminosities at z = 1.0. The NewAthena lower luminosity limit is based on the flux limit for a deep survey with the maximum proposed angular resolution (5 arcsec).⁶ Similarly, as a representative example of the NASA X-ray probe-class missions that are currently under review, we plot the flux limits for the proposed AXIS deep and wide surveys (Reynolds et al. 2023). Though note that (in addition to infrared missions), there are several other proposed X-ray probes for the APEX programme that will target dual SMBH and close pairs, e.g. HEX-P which would have the ability to spatially resolve dual AGN targeting the hard X-rays (Civano et al. 2024; Pfeifle et al. 2024) and LEM which could spectrally resolve dual AGN at high redshifts (Pfeifle et al. 2023).

This demonstrates the crucial need for accurate accretion disc modelling as the different disc models make significantly different predictions for the detectability of the binaries. In particular, at z = 1, even with the state-of-the-art NASA probe concepts, the binaries would not be detectable according to our unified accretion disc model – regardless of whether wind loss is included. For nearby galaxies within z = 0.25, which may harbour a significant SMBH binary population as revealed by the population modelling based on gravitational wave background constraints from nanoGRAV (see Agazie et al. 2023b, though note that strongest constraints are for more massive SMBHs than the ones considered here), X-ray detection prospects are more favourable, although the proposed probes may struggle to detect these objects if wind loss is significant. With the thin disc model, the predicted X-ray luminosities are much higher and may be detectable out to z = 1, with the proposed NASA X-ray Probe-class missions.

Radio interferometers will play a particularly important role in constraining wide binaries, such as the ones simulated in this work, as very long baseline interferometry (VLBI) allows us to directly resolve these systems, whilst future X-ray missions will still be limited to kpc scales (De Rosa et al. 2019). Current radio interferometers have already identified parsec scale binaries at z ∼ 0.1 and future facilities will push this to much higher redshifts with SKA resolution sufficient to at least z ∼ 1 and the ngVLA resolving wide SMBH binaries out to redshifts z ≳ 3 (Burke-Spolaor et al. 2018; Bogdanović, Miller & Blecha 2022).

To this end, we plot radio luminosity predictions for the different disc states in the bottom row of Fig. 9. Both analytical theory (Meier 2001) and GRMHD simulations (Tchekhovskoy, Narayan & McKinney 2010) predict that geometrically thick accretion discs should be more efficient at jet production and may even explain the radio loud/quiet AGN dichotomy. The spin-driven jet model from Meier (2001) is also applied in Barausse (2012) and more recently in Mangiagli et al. (2022), and we follow their disc-state-dependent formalism here for our binary electromagnetic counterpart predictions. We then convert the jet luminosities into radio luminosities using the scaling relationship between jet power and synchrotron luminosity from Cavagnolo et al. (2010). Following Mangiagli et al. (2022), we also take into account beaming effects assuming a typical Lorentz factor of Γ = 10 for AGN.

For comparison, we also again indicate detectable luminosities for future radio facilities based on typical sensitivities in the GHz band (where most of the jet emission is expected to occur) for SKA (∼200–2000 nJy, see Braun et al. 2019) and the ngVLA (∼40–400 nJy, see Plotkin & Reines 2018), as collated by Latif, Whalen & Mezcua (2024). As a note of comparison, the upper limit for detectable luminosities with SKA corresponds to the lower limit for LOFAR.

Based on current SMBH spin constraints (Reynolds 2019), we have initialized our SMBHs as highly spinning (a = 0.9) which favours strong jet production (Blandford & Znajek 1977; Blandford & Payne 1982). In the ADIOS regime, the jets are further boosted due to the smaller solid angle subtended by the base of the jet, increasing the jet power (Tchekhovskoy et al. 2010).

Hence we obtain very high radio jet luminosities in the ADIOS regime, even when considering wind loss, which exceed the thin disc predictions by several orders of magnitude. This has important implications for VLBI SMBH binary searches, which will mainly be limited by angular resolution for highly spinning SMBHs in the truncated or ADIOS state. For thin discs, however, these searches will also be fundamentally limited by sensitivity. In our case, whilst the ngVLA may be expected to resolve a parsec-scale SMBH binary out to z ∼ 3, with the thin disc model we would not expect to detect the jet luminosities of this system.

One important caveat is the role of magnetic fields in jet production. Indeed recent GRMHD simulations have demonstrated that thin discs can have much more efficient jets if the inner region of the disc is magnetically saturated (‘MAD state’) which may provide a potential explanation for radio-loud quasars (Liska et al. 2019; Ricarte, Narayan & Curd 2023). Incorporating these effects is beyond the scope of this work though we note that a stronger magnetic field would increase the spin-driven jet power for both disc types following the Blandford-Znajek mechanism (Blandford & Znajek 1977), so that the qualitative statement of stronger jets in the ADIOS regime would likely still hold.

As a final caveat, we note that these strong jets would spin down the SMBH, decreasing the jet power. As we are only estimating the jet power in post-processing we do not include this effect here, however, in the future we plan to couple our unified accretion disc model to a resolved jet feedback implementation (also see Talbot, Bourne & Sijacki 2021; Talbot, Sijacki & Bourne 2022, 2024) and include the jet spin down self-consistently.

Overall, we emphasize that the thin disc model and truncated disc model for the same inflow rate predict markedly different counterparts. For the former, we obtain relatively bright X-ray AGN but weak radio jets, whilst the latter produces strong radio jets and weak X-ray emission. Hence it is of critical importance to incorporate as realistic AGN accretion disc models as possible when making electromagnetic counterpart predictions for SMBH binaries.

4.3 Evolution of SMBH spin and accretion disc angular momentum

Here, we explore the evolution of the SMBH and accretion disc angular momentum, which is governed by the external inflows from the mini discs, the accretion of gas angular momentum from the accretion disc onto the SMBH and the mutual Lense–Thirring torques. The latter two are set by the assumed subgrid model for the accretion disc. In particular the Lense–Thirring torques, which arise from a misalignment between the SMBH and disc angular momentum versors, operate in two very different physically regimes depending on the disc model assumed. In the thin disc case, these differential torques are communicated via viscosity, leading to the so-called Bardeen–Petterson configuration, where the inner part of the thin disc is aligned with the SMBH versor up to a warp radius r_warp. In the ADIOS regime, the whole flow is strongly coupled via pressure waves and precesses as a solid body, only aligning very slowly with the total angular momentum versor. See Section 2.5.2 for a detailed description of the subgrid treatment of Lense–Thirring precession in the different disc states.

4.3.1 Single SMBH simulations: angular momentum versor evolution

We first consider the angular momentum evolution in the single SMBH simulations. Here the spin alignment is mainly driven by the Lense–Thirring torques as both the SMBH accretion rates and the external gas inflow rates are negligible (see Section 4.2.1). Fig. 10 shows polar projections of the SMBH and accretion disc versors’ time evolution for the thin disc model (single_thindisc, top row), unified accretion disc model without wind loss (single_uniadios, middle row), and unified accretion disc model with a massive initial disc and without wind loss (single_uniadios_massive, bottom row). The colour coding of the tracks indicates the time evolution and the black diamond and square markers indicate the initial and final versor orientation, respectively. The blue filled circle indicates the location of the total angular momentum versor j_tot. Note that the total versor does not evolve throughout the simulation due the extremely low gas inflow levels.

For the thin disc model, single_thindisc, both the SMBH and the accretion disc versor are slowly aligning with the total angular momentum versor due to the Bardeen-Petterson alignment. For the unified accretion disc model, single_uniadios, we instead obtain efficient precession of the inner hot flow. Note that the wind loss model (single_uniadios_winds, not shown), results in an identical versor evolution as the precession frequency only depends on the disc and SMBH state variables (i.e. the masses and angular momenta) and not the current accretion rate (see equation 33). Since the state variables are only changing negligibly due to the low accretion rates, the solid body precession pattern is therefore identical with or without wind loss in these simulations. Finally, we also consider the versor evolution for different initial disc masses. With the massive disc set-up, single_uniadios_massive, the Eddington fraction is initially much higher (see Section 4.2.1), so that we obtain extremely rapid realignment of the SMBH and disc versor with the total versor in the Bardeen–Petterson configuration. As the accretion rate drops, the disc switches to the truncated disc mode and starts precessing as a solid body, albeit with an extremely small precession angle due to the previous realignment. With the light disc model (single_uniadios_light, not shown) the precession time-scales are much longer than the simulation time and the alignment time-scale exceeds the Hubble time, so that for the duration of the simulation the versor orientation remains constant.

4.3.2 Single SMBH simulations: precession and alignment analysis

We investigate these different evolution scenarios more quantitatively in Fig. 11, which shows the polar angle of the total, SMBH and accretion disc versors as a function of time and Fig. 12, which shows the time evolution of the precession and alignment time-scales in the Bardeen–Petterson and solid body precession regimes, respectively.

The left-hand panel of Fig. 11 shows how in the single_thindisc simulation, the disc and the SMBH versor are torqued towards the total versor by the Bardeen–Petterson effect. However, this alignment process is very slow with time-scales T_{align, BP} ≳ 10³ Myr (see left-hand panel of Fig. 12). Meanwhile the precession time-scales are even longer, so that the versor evolution in the single_thindisc set-up is dominated by the alignment process, with precession only playing a negligible role.

For the single_uniadios model, the disc and SMBH versors are precessing in antiphase to one another in the solid body regime (see middle panel of Fig. 11), with the initial precession time-scale T_{prec, SB} ≳ 20 Myr, which then increases to T_{prec, SB} ∼ 10³ Myr at the end of the simulation (see Fig. 12), so that the BH-disc-system never completes a full precession period.

Finally, with the single_uniadios_massive model, we have a phase change where the disc goes through rapid Bardeen–Petterson alignment in the first Myr with T_{align, BP} ∼ 1 Myr and then enters the solid-body precession regime at t ∼ 10 Myr, undergoing small-scale precession due to previous efficient realignment with the period increasing from T_{prec, SB} ≲ 10 Myr to T_{prec, SB} ∼ 10³ Myr.

4.3.3 Binary SMBH simulations: angular momentum versor evolution

Following on from the simpler single SMBH case, we move on to analyse the angular momentum evolution of the SMBH and accretion disc in the binary simulations. Here the inflow rates from the mini discs can at times be significant, so that the versor evolution is not just driven by the Lense–Thirring torques but also by external inflows of gas angular momentum.

Fig. 13 shows polar projections of the SMBH and accretion disc versors’ time evolution for the aligned q1i0 (left panel) and misaligned q3i45 binary runs (right panel). The blue-purple tracks represent the versor evolution of the primary, whilst the orange-red tracks show the versor evolution of the secondary SMBH. The red and blue dotted lines show the evolution of the total angular momentum versorj_tot.

The evolution of the total angular momentum versors for a given binary configuration is virtually identical for all of the subgrid models explored as it is governed by the external inflows from the mini discs. For the aligned binary, bin_q1i0, we find efficient realignment of the total versors with the poles as we have coherent accretion from the mini discs, whilst for the misaligned binary, bin_q3i45, the accretion is less coherent as the binary realigns with the outer CBD (see Bourne et al. 2023, for details).

The first row shows the runs with the thin disc model. As with the single SMBH set-up, we have Bardeen–Petterson alignment of the SMBH and disc angular momentum versors with the total angular momentum versors. However, as the total angular momentum versors are being reoriented towards the mini disc angular momentum and the subgrid disc is initially assumed to be aligned with the mini discs (see Section 3.3.2), the alignment between the disc and total versors is accelerated by the external inflows. For the SMBH versors, on the other hand, the external inflows lead to less efficient alignment as the total versor recedes from the SMBH versor.

The second row shows the simulations performed with the unified accretion disc model without wind loss and with standard precession. For the bin_q1i0 set-up, after a short Bardeen–Petterson alignment period, the SMBH and disc versors are rapidly precessing. The SMBH versors’ precession angles remain roughly constant, whilst for the disc versor the precession angle is steadily decreasing as the disc and total angular momentum versor align due to accretion from the mini discs. For the bin_q3i45 simulation, the picture is more complex as the system alternates between Bardeen–Petterson alignment and solid-body precession, so that the precession angle decreases for both the SMBH and disc versors.

The simulations shown in the third row were also performed with the unified accretion disc model, but account for wind loss with the ADIOS model. Note that this slows down the SMBH mass growth rate, which decreases the truncation radius so that the inner hot flow is more likely to be smaller than the aligned region of the warped outer thin disc, so that the inner hot flow is also twisted into alignment. Indeed for both binary configurations, the SMBH–accretion disc system tends to spend more time in the Bardeen–Petterson alignment phase with the ADIOS wind loss being taken into account. The impact is particularly dramatic for the secondary SMBH in the bin_q1i0 set-up. Here the combination of solid body precession followed by Bardeen–Petterson alignment realigns the SMBH and disc versors even more efficiently than the thin disc only case.

The simulations presented in the fourth row were performed with a reduced solid body precession rate for the truncated disc case, following the results from Bollimpalli et al. (2023), who find that the precession rate of the inner hot flow is reduced by 95 per cent due to torques from the outer thin disc (see Section 2.5.2 for details on the solid body precession model). With this model, the system never completes a full precession cycle within the simulation time and there is only very slow alignment from the short Bardeen–Petterson phases. The equivalent set-ups with wind loss (uniadios_winds_slowprec, not shown here) are very similar. Since the precession proceeds very slowly, the duration of the Bardeen–Petterson alignment has a much less significant impact on the overall evolution.

4.3.4 Binary SMBH simulations: precession and alignment analysis

We now turn to investigate the complex precession behaviour of the binary SMBHs more quantitatively. Fig. 14 shows the time evolution of the versors’ polar angle in the left-hand panels. Due to efficient angular momentum accretion from the mini disc, the total angular momentum versor evolves significantly in the binary simulations. To control for this effect, we also consider separately the time evolution of the ‘offset angle’ between the BH/disc versors and the total angular momentum versor in the right-hand panels. Moreover, Fig. 15 shows the alignment and precession time-scales in the Bardeen–Petterson and solid body regimes as well as the characteristic external inflow time-scale (defined as |$T_\mathrm{inflow} = M_\mathrm{d}/\dot{M}_\mathrm{inflow}$|⁠) and the binary inspiral time-scale.

We first focus on the angular momentum evolution in the equal-mass, aligned binary case (bin_q1i0). Due to the inflows from the mini discs, the angular momentum versorj_d (and therefore the total versorj_tot) is driven towards alignment with the pole. As discussed in the previous section,j_tot may only change due to external inflows, and therefore thej_tot evolution follows the same pattern in all of the bin_q1i0 runs. From Fig. 15, we can see that the mean inflow time-scale varies between ∼20 and 100 Myr (light-blue and light-red solid lines) and hencej_tot aligns with the poles during the duration of the simulation.

The alignment and precession of the accretion disc and SMBH angular momentum versors, however, depends significantly on the accretion disc model assumed. With the thin disc model (first row), the disc is in the Bardeen–Petterson configuration, where the SMBH and spin versors may efficiently realign, both being driven towardsj_tot on an alignment time-scale of a few 100 Myr (see light-blue and light-red dotted lines in Fig. 15). Considering the evolution of the angle betweenj_tot and the disc and SMBH versors in the right-hand panel of Fig. 14, we can see that the disc efficiently aligns withj_tot as both the external inflows and the Bardeen–Petterson alignment work in tandem. For the SMBH, on the other hand, the inflows drivej_tot away from the SMBH versor so that there is a phase where the misalignment angle increases. Though eventually, as the inflow time-scale increases, Bardeen–Petterson alignment becomes more efficient. Overall the SMBH spin evolution is dominated by the alignment process, as the precession time-scale in the Bardeen–Petterson configuration is much longer (T_{prec, SB} ∼ 10³ Myr), exceeding the runtime of our idealized set-up.

With the unified accretion disc model (second to fourth rows), we get significantly different behaviour. As the disc is in the truncated state, the relevant Lense–Thirring precession regime is determined by the ratio between the truncation radius and the warp radius from the Bardeen–Petterson alignment process of the outer thin disc (shown as inset panels in Fig. 14). If the truncation radius is larger than the warp radius, then the outer disc is feeding misaligned material to the inner hot flow which therefore enters the solid-body-precession mode, characteristic of misaligned thick discs (see Section 2.5.2 for details).

With the uniadios model (second row), we have a short period of Bardeen–Petterson alignment followed by solid-body precession for the remainder of the simulation. In the case of the disc the precession angle is steadily decreasing as the total angular momentum versor is being driven towards the poles. For the SMBH versor, on the other hand, as the precession is in antiphase (close to the pole), the total angular momentum versor is moving towards the SMBH versor and the precession angle is decreasing. When the SMBH precession has completed a full cycle (far from the pole), the total angular momentum versor is moving away from the SMBH versor and the precession angle increases. These small oscillations can be seen in the third panel of the second row in Fig. 14, which as expected, are shifted from the precession motion by half a phase. Overall this means that the mean precession angle is approximately conserved, with the SMBH versor following the translational motion of the total versor, thereby also aligning with the poles and the disc angular momentum.

With the uniadios_winds model (third row), we also take wind-loss into account which leads to lower SMBH masses and therefore lower (auxiliary) Eddington fractions and smaller truncation radii. For the bin_q1i0 simulation, the truncation radius to warp radius ratio is relatively close to one (see insets in left-hand panels of Fig. 14), so that the system is highly sensitive to small changes in the truncation radius and the secondary SMBH dips back into the Bardeen–Petterson regime at t = 20 Myr for the uniadios_winds set-up. This change-over occurs when the SMBH is in antiphase (close to the poles), so that now the SMBH and the total versor are being driven towards one another by both the external inflows and the Bardeen–Petterson alignment. Meanwhile, whilst the disc is far away from the pole, the disc angular momentum is directly affected by the external inflows (with a higher fractional change than the total versor) and therefore also realigns rapidly with the poles. At t = 40 Myr, when the system re-enters the solid-body precession regime, both the disc and the SMBH versor are within 10° of the total versor, so that we only have small-scale precession for the remainder of the simulation. Remarkably, the interplay between the Bardeen–Petterson regime and solid-body precession leads to even more efficient alignment than the thin disc model.

In the fourth row of Fig. 14, we also show the time evolution of the polar and precession angles for the slow precession model for truncated discs following Bollimpalli et al. (2023). Due to the longer precession period, the SMBH versor remains in the phase space where the total versor is receding and therefore the offset angle for the SMBH is steadily increasing. Hence with this model, the final SMBH versor is the least aligned out of all set-ups considered. Note that including wind loss (bin_q1i0_uniadios_winds_slowprec, not shown here) leads to very similar results as in neither case a significant re-orientation of the SMBH versor is achieved.

Next, we analyse the precession and alignment patterns in the non-equal-mass, misaligned binary simulations bin_q3i45 shown in the second columns of the two panels in Fig. 14. We find that in this misaligned scenario the primary total versor is initially already in close alignment with the disc, leading to efficient realignment of the primary disc and SMBH versors with Δθ_offset < 20° for all accretion disc models explored (see right-hand panel of Fig. 14). For the secondary SMBH, the initial offset angle is much larger and the inflow time-scale increases dramatically at t = 40 Myr (see Fig. 15), so that alignment is much less efficient. We also note that the accretion from the mini discs for the misaligned simulations is highly stochastic and incoherent (Bourne et al. 2023), so that the evolution of the total versors (see left-hand panel of Fig. 14) slightly differs between the different simulations performed.

With the thin-disc model, the initial Bardeen–Petterson alignment time-scale for both SMBHs is a few 100 Myr. This remains relatively constant for the primary, however, for the secondary SMBH the alignment time-scale steadily increases from t ∼ 40 Myr onwards, significantly slowing down the disc alignment. Note that since the inflows for the misaligned binary case are much less coherent we do not observe the slowing down of the SMBH versor alignment that we saw with the equal-mass binary. As discussed in Bourne et al. (2023), there is some small precessional motion for both the disc and the total versor due to the incoherent inflows from the mini discs which are realigning with the outer CBD. It is worth noting that the realignment process slightly differs in our simulation from the set-up presented in Bourne et al. (2023) as they have a different randomly chosen initial SMBH versor orientation, which is much more significantly offset for the primary SMBH leading to less efficient alignment with the mini discs.

With the unified accretion disc model, we obtain several changeovers between the Bardeen–Petterson and solid-body regime due to the stronger variation of the Eddington fraction in the misaligned case (see Fig. 8). First, we focus on the set-up without wind loss, uniadios, in the second row of Fig. 14. For the primary and secondary SMBH, we initially have a period of Bardeen–Petterson alignment and then transfer into the solid-body precession at t = 5 Myr and t = 15 Myr, respectively. During the solid-body precession period, we again have the precession angle decreasing for the disc versor as the external inflows provide additional realignment, whilst the SMBH precession angle remains constant. The primary and secondary SMBH then switch back to the Bardeen–Petterson configuration at t ∼ 20 Myr, further aligning with the total versor until t = 40 Myr and t = 50 Myr, respectively. Having further realigned, the SMBHs then carry out solid-body precession at a reduced angle for the remainder of the simulation. Note that the disc precession angle only marginally decays in this second phase of solid body precession as the angular momentum inflow rates are less coherent.

With wind loss, uniadios_winds (third row), we obtain qualitatively similar behaviour for the primary SMBH, although the second Bardeen–Petterson alignment phase is significantly prolonged as the SMBH mass grows more slowly, leading to a smaller truncation radius. This leads to a reduced precession angle in the second solid-body precession phase. For the secondary SMBH, on the other hand, the smaller truncation radius means that the BH-disc system remains in the Bardeen–Petterson configuration until t ∼ 40 Myr and then transfers to the solid body precession regime with a significantly reduced precession angle.

Finally, we also consider the modified set-up, where the precession rates are reduced by 95 per cent in the truncated disc mode, uniadios_slowprec. As with the bin_q1i0 set-up, this means that the SMBH versors are lagging behind the total versor alignment due to external inflows throughout, so that the precession angle increases, leading to a larger final offset angle for the primary SMBH. The secondary SMBH experiences a similar albeit much smaller effect as it spends significantly less time in the solid body precession phase.

4.3.5 Summary and implications of angular momentum evolution

Figs 10 and 13 clearly show that the accretion disc subgrid model has a significant impact on the spin alignment and precession process. In the truncated disc regime, precession is much more prevalent with the SMBH spin versors covering a large solid angle. This has crucial implications for distributed feedback, in particular in the context of spin-driven jets (see Section 5.2.2 for a more detailed discussion of implications for AGN feedback). Also note that VLA and Merlin observations have uncovered precessing jets with inferred precession periods of 1–10 Myr (Krause et al. 2019). Whilst this could be attributed to geodetic precession, it would also be consistent with solid-body precession in the truncated disc regime within our framework and LOFAR is expected to shed more light on the origin of these variability patterns (e.g. Horton, Krause & Hardcastle 2023).

Furthermore, it is interesting to note that despite the alignment torques between the accretion disc and SMBH being extremely weak in the ADAF/ADIOS regime, the SMBH may still realign with the disc if the external inflow time-scale matches the precession time-scale. Furthermore, Bardeen–Petterson alignment may be accelerated if the solid body precession reorientates the SMBH versor towards the mini discs. On the other hand, if the precession process is slowed down significantly due to additional torques in the truncated disc state (Bollimpalli et al. 2023), then the spins will not align within the inspiral time-scale, with crucial implications for recoil predictions, influencing merger rates and the gravitational wave background.

5 DISCUSSION

We have developed a unified accretion disc model suitable for galaxy formation simulations that can resolve the outer edge of the accretion disc. Our work builds on the thin disc model prescription presented in Fiacconi et al. (2018) and extends this scheme by harnessing the results from high-resolution GR(R)MHD simulations of accretion discs. Our model spans a range of radiative regimes from the optically thin, geometrically thick discs to the optically thick, geometrically thin discs via a truncated disc configuration. We also track the SMBH spin evolution with a disc-state-dependent Lense–Thirring precession model. Our predictions for the SMBH spin alignment and precession as well as SMBH luminosities are highly sensitive to the disc state with important implications for future multimessenger observational campaigns which we outline in Section 5.1. We also discuss the caveats of our model and important future extensions in Section 5.2.

5.1 Implications for observations

Parsec-scale SMBH binaries have proved extremely challenging to observe, with only a handful of detections so far (Rodriguez et al. 2006; Bansal et al. 2017; Kharb et al. 2017), as these are in the regime where extremely high angular resolution is necessary to resolve the binary, however, the SMBHs are not close enough to display significant variability on observable time-scales. This makes radio interferometers the tool of choice for searching for these binaries and future facilities such as SKA or ngVLA are expected to significantly increase sample sizes (e.g. Burke-Spolaor et al. 2018), though see a cautionary note from D’Orazio & Charisi (2023).

To estimate radio luminosities, it is crucial to have reliable estimates for the mass accretion rates, spin orientations, and solid angles (e.g. Mangiagli et al. 2022). With our unified accretion disc model, we can estimate all of the relevant quantities and make predictions for the jet luminosities depending on the disc state which reveals a much more promising avenue for detecting SMBH binaries accreting at intermediate luminosities than the thin disc modelling (see Fig. 9).

According to analytical calculations, our binary system should have an inspiral time-scale of ≲ 10³ Myr (see Fig. 15), so this represents a precursor to the gravitational-wave-emitting binaries. Nevertheless it is crucial to make predictions for these precursor luminosities, in particular for X-ray confirmations (Piro et al. 2023).

Furthermore, accurate AGN modelling will be crucial for dual AGN (and close pairs) searches with JWST where recent tantalizing observations have revealed hints at a large fraction of candidate dual AGN (Maiolino et al. 2023; Übler et al. 2023), highlighting the importance of imaging spectroscopy in tracking SMBH mergers throughout cosmic time. Combined with results from gravitational wave interferometers, these JWST constraints can provide crucial constraints on the merging SMBH population (e.g. Padmanabhan & Loeb 2024).

The AGN lifetime and dynamical evolution of the binary in our simulations also provide important information for astrometric searches for SMBH binaries with Gaia (e.g. D’Orazio & Loeb 2019).

Finally, our spin evolution predictions play an important role in estimating the gravitational wave background signal as they determine the recoil velocity of the SMBH remnant – in particular large and misaligned SMBH spins may lead to recoil velocities of |$\sim 5000 \ \mathrm{km \, s^{-1}}$| (Lousto & Zlochower 2011, 2013; Lousto & Healy 2019) which would have a significant impact on SMBH occupation fractions and merger rates.

5.2 Caveats and future work

5.2.1 Accretion modelling

We have extended the thin disc model from Fiacconi et al. (2018) into the heavily sub-Eddington regime based on the ADAF/ADIOS models. However, we note that there are several other disc models that have been proposed for the radiatively inefficient accretion regime, including convection-dominated accretion flows (CDAFs; Narayan, Igumenshchev & Abramowicz 2000; Quataert & Gruzinov 2000) or more recently simple convective flows (Xu 2023) and strongly magnetized Bondi models (Cho et al. 2023). Li, Ostriker & Sunyaev (2013) carried out an extensive simulation suite, capturing the accretion flows from Bondi scale to the small scale and demonstrated how the inflow and outflow patterns will change as a function of cooling. Similarly, they show how strong outflows suppress the Bondi rate, however, they also demonstrate a strong dependence on cooling. Note that in principle, we could incorporate these accretion flow scenarios within our framework by varying the s parameter from the ADIOS wind loss model, which effectively sets the flow density profile. However, this is beyond the scope of this paper and we plan to investigate the impact of s parameter choices more systematically in the future.

Whilst this work has focused on the radiatively inefficient regime, in the future it will also be important to further consider the disc modelling in the radiatively efficient regime. Several recent studies have questioned and improved the standard thin α-disc model (e.g. Jiang et al. 2019) and this will require careful consideration to improve the predictions of our model in the quasar regime.

Misaligned discs experience radius-dependent Lense–Thirring torques. How these differential torques are communicated depends on the accretion disc regime. In our unified accretion disc model, we evolve the SMBH spin according to the Bardeen–Petterson alignment (Bardeen & Petterson 1975) in the thin disc state and following the solid body precession model from Fragile et al. (2007) in the thick disc state. We also include the possibility of slowed down solid body precession in the truncated disc regime (Bollimpalli et al. 2023). However, GR(R)MHD simulations paint an even more complex picture, including disc tearing (e.g. Liska et al. 2023b) which could suppress Bardeen–Petterson alignment (e.g. Steinle & Gerosa 2023) and induce phase lags between the inner hot flow and outer thin disc (e.g. Liska et al. 2023a). Furthermore, if the thick disc is in the magnetically arrested state, strong magneto-spin alignment could significantly suppress solid body precession (e.g. McKinney et al. 2013; Chatterjee et al. 2023).

Finally, another crucial future extension to our model will be to incorporate the super-Eddington regime. This accretion disc regime may be particularly important in the early Universe, where JWST is starting to uncover a population of overmassive AGN (compared to the local scaling relations) which could stem from sustained super-Eddington accretion episodes (e.g. Goulding et al. 2023; Maiolino et al. 2023). Furthermore, cosmological simulations show that there is likely enough gas present around the early SMBHs for super-Eddington accretion (e.g. Rennehan et al. 2023; Bennett et al. 2024; Juodžbalis et al. 2024). Note that many of the key physical properties of our ADIOS modelling in the sub-Eddington regime will also apply to the super-Eddington regime. The latter is advection-dominated due to the long radiative diffusion times (rather than long cooling times, see Yuan & Narayan 2014), so that the radiative efficiency drops in the super-Eddington regime and accretion (as well as jets) cause the SMBH to spin down (also see Ricarte et al. 2023; Lowell et al. 2024).

5.2.2 Feedback modelling

Within our current framework we only model gas accretion through a unified accretion disc model. For some model variations we take into account mass-loss through a wind, however, we do not yet re-inject the wind mass into the simulation. The obvious next step is to incorporate AGN feedback mechanisms into our model. The most straightforward extension here would be to incorporate ADIOS winds with the mass loading directly informed by the wind mass-loss in our model. The momentum and energy loading of the wind could then be inferred from the kinetic wind efficiencies obtained from GRMHD simulations (e.g. Tchekhovskoy et al. 2012; Sa̧dowski et al. 2013; Sa̧dowski & Gaspari 2017). We note that the thin disc component is also expected to drive winds, in particular at high Eddington fractions, and including these winds would be an important extension for future work (also see Sala et al. 2021), representing the ‘quasar mode’ feedback employed in cosmological simulations which is especially crucial in the early Universe. In addition to wind feedback, radiative feedback would be another important channel to consider. Here our unified accretion disc model provides the ideal starting point as we could base the AGN emission on the subgrid accretion disc state.

Moreover, AGN jets are an ubiquitous phenomenon and are likely powered by the SMBH spin energy extraction via the Blandford-Znajek (BZ) mechanism (Blandford & Znajek 1977). Talbot et al. (2021) have coupled spin-driven BZ jets self-consistently with the thin disc model from Fiacconi et al. (2018) and considered their evolution in Seyfert galaxies (Talbot et al. 2022) and a galaxy merger (Talbot et al. 2024). In future work, it would be timely to extend this scheme to include spin-driven jets across different accretion disc types (also see Huško et al. 2022, 2024). In particular, the solid body precession model could have a crucial impact here leading to more ‘distributed’ jet feedback. Note that the precession periods in the solid body regime are highly uncertain (see Fig. 4), and there are hints at Lense–Thirring precession in AGN on the scale of ∼10 yr (Cui et al. 2023; von Fellenberg et al. 2023). In terms of AGN subgrid models, this would imply that for time-scales relevant for galaxy formation the AGN feedback would be effectively distributed over a cone, which may significantly increase its impact and potentially lead to more direct interaction with the interstellar medium as demonstrated in idealized simulations of precessing jets by Su et al. (2021).

5.2.3 Porting the unified accretion disc model to cosmological simulations

In addition to improving the subgrid SMBH modelling, it will also be crucial to move towards more realistic environments and test different accretion regimes. We found that more chaotic accretion, such as in the misaligned binary case, can have a significant impact on luminosities and spin evolution so a natural next step would be to test this framework within galaxy merger simulations (also see Liao et al. 2023; Talbot et al. 2024). Furthermore, extending our model to cosmological zoom-in simulations will allow us to track the mass and angular momentum flows from the cosmological environment down to the scale of the accretion disc. For the zoom-in simulations, the super-Lagrangian refinement technique will permit us to reach the required resolution (also see Curtis & Sijacki 2016a; Bourne, Sijacki & Puchwein 2019; Anglés-Alcázar et al. 2021; Bourne & Sijacki 2021; Hopkins et al. 2024a, b). For large cosmological volumes, applying the super-Lagrangian refinement technique to the central regions of tens of thousands of galaxies, however, is computationally infeasible. One avenue to circumvent these issues could involve using machine-learning techniques such as generative adversarial networks or learned coarse models using convolutional neural networks (e.g. Li et al. 2021; Stachenfeld et al. 2021). This would require appropriate training data sets covering the vast range of relevant scales and different environments and could build on recent pilot projects pioneering this approach for subgrid models in planet formation (e.g. Pfeil et al. 2022).

Finally, our unified accretion disc model could be modified for lower resolution applications, in particular for the SMBH mass evolution. Here more coarse-grained prescriptions such as the Bondi model or the torque-based accretion model could provide an estimate of the inflow rate onto the disc. The accretion disc could then act as a gas reservoir that gets drained on the viscous time-scale which gets calculated in a disc-state dependent fashion. Similar approaches have already been introduced in the literature (e.g. Power et al. 2011; Wellons et al. 2023), however, these employ either a constant viscous time-scale or a viscous time-scale based on the thin α-disc model, whereas our approach would also be directly applicable in the radiatively inefficient regime and therefore appropriate for both the quasar and radio mode in cosmological simulations.

6 CONCLUSIONS

We have developed a novel unified accretion disc model designed to be embedded in galaxy formation simulations. Drawing inspiration from high-resolution GR(R)MHD simulations, our model tracks mass flows through three distinct disc states and evolves the SMBH spin based on disc-state-dependent Lense–Thirring precession models. This innovative accretion disc framework has been integrated into the moving mesh code arepo, and our work includes a comprehensive suite of validation simulations. These simulations encompass idealized scenarios featuring single SMBHs and binary SMBHs embedded in gaseous discs, with variations in mass ratio, inclination angle, and subgrid accretion disc models. Our key findings are summarized below.

• The ADIOS model predicts substantial wind loss in the thick disc regime, leading to a significant suppression of gas mass accretion rates onto SMBHs by an order of magnitude or more. Additionally, the decrease in radiative efficiency within the advection-dominated accretion regime results in luminosity estimates that are decreased by another order of magnitude. Consequently, predictions for the luminosities of electromagnetic counterparts of SMBH binaries exhibit significant sensitivity to the disc state.

• Hence, it is imperative to incorporate these realistic accretion disc models when forecasting outcomes for future missions both for single and binary SMBHs. Specifically, our unified accretion disc model reveals a drastic reduction in the redshift range for accretion disc luminosities compared to predictions based on the thin α-disc model, regardless of whether winds are accounted for. This has important implications for SMBH binary searches in the X-ray, optical, and infrared. Furthermore, the predicted spectrum undergoes significant shifts, with the luminosity being dominated by the hard inner hot flow.

• For the radio jet luminosities, on the other hand, the situation is reversed with the unified accretion disc model predictions much more favourable than the thin disc predictions, since the geometrically thick ADIOS flow leads to significantly more efficient jet production. We predict that future VLBI searches with SKA and the ngVLA should significantly increase the known samples of parsec-scale binaries.

• The evolution of the spin magnitude depends only mildly on the disc state, primarily due to relatively low external inflow rates in all cases considered here. We note however, that depending on the disc state adopted we predict for the same large-scale gas inflow, the SMBH spin may either increase or decrease. Clearly, in case of higher and more variable accretion rates, significant differences in the spin magnitude evolution are expected.

• Importantly, the SMBH spin direction’s evolution is markedly influenced by the assumed disc model, ranging from slow alignment in the Bardeen–Petterson configuration to fast solid-body precession in the truncated disc regime. This has important implications for the impact of jet-mode feedback on galaxies and the circumgalactic medium, as more rapid precession would result in distributing the jet energy over a larger volume.

• For binary systems, the variations in spin evolution could significantly impact recoil velocities, consequently influencing BH merger rates and the gravitational wave background.

This model provides a crucial advancement in AGN accretion modelling, and there are numerous promising avenues for extensions, including spatially resolved feedback. In the near future, next-generation cosmological simulations, incorporating this type of subgrid physics, will be able to give us unprecedented insights into the complex interplay of SMBHs and their host galaxies, and provide the scientific community with much needed robust predictions for the era of multimessenger science.

ACKNOWLEDGEMENTS

The authors would like to thank the referee for a thorough review which improved the quality of the manuscript. Furthermore, the authors are grateful to Alexander Tchekhovskoy, Beverly Lowell, Matthew Liska, Jenny Greene, Shy Genel, Angelo Ricarte, Doug Rennehan, and Samuel Turner for useful discussions during the development of this work. The authors would also like to thank Lucy Reading-Ikkanda at the Simons Foundation for creating the illustration of our model in Fig. 1. SK is supported by a Flatiron Research Fellowship at the Flatiron Institute, a division of the Simons Foundation and a Junior Research Fellowship from St Catharine’s College, Cambridge. SK, MAB, and DS acknowledge support by European Research Council Starting Grant 638707 ‘Black holes and their host galaxies: coevolution across cosmic time’. MAB and DS additionally acknowledge support from the Science and Technology Facilities Council (STFC). KP is supported by the Simons-NSBP Scholars Program. This work was supported by the Simons Collaboration on ‘Learning the Universe’. The computations in this work were, in part, run at facilities supported by the Scientific Computing Core at the Flatiron Institute, a division of the Simons Foundation. This work used the DiRAC Data Intensive service (CSD3) at the University of Cambridge, managed by the University of Cambridge University Information Services on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC component of CSD3 at Cambridge was funded by BEIS, UKRI, and STFC capital funding and STFC operations grants. This work used the DiRAC Memory Intensive service (Cosma8) at Durham University, managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC service at Durham was funded by BEIS, UKRI, and STFC capital funding, Durham University and STFC operations grants. DiRAC is part of the UKRI Digital Research Infrastructure.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX: CND RELAXATION FOR SINGLE SMBH SIMULATION

Here we present more detailed plots illustrating the evolution of the cavity for the single SMBH simulations during the refinement relaxation and once the accretion disc subgrid model is switched on. Note that we deliberately do not let the cavity fully fill with gas before activating our model as we are mainly interested in the low-accretion-rate regime.

Fig. A1 shows gas density projections of the cavity during the refinement relaxation (panels 1–4) and the production run with the subgrid disc model (panels 4–5). After the refinement model is activated, several instabilities develop, amplified by cooling from the CND onto the central region (also see Bennett & Sijacki 2020; Anglés-Alcázar et al. 2021). This leads to the formation of spiral arms and a dense core which then propagates outwards (see third panel). By the time that the subgrid model is activated (fourth panel), this overdensity has dissipated. However, during the simulation, the cavity continues to fill with gas leading to a mild increase in inflow rates.

Fig. A2 shows this more quantitatively, plotting the gas density, cell sizes, sound speed, and azimuthal velocity as a function of radius during the refinement relaxation phase (green lines) and the subgrid model simulation phase (purple lines).