pAGN : the one-stop solution for AGN disc modeling

Models of accretion discs surrounding active galactic nuclei (AGNs) find vast applications in high-energy astrophysics. The broad strategy is to parametrize some of the key disc properties such as gas density and temperature as a function of the radial coordinate from a given set of assumptions on the underlying physics. Two of the most popular approaches in this context were presented by Sirko & Goodman (2003) and Thompson et al. (2005). We present a critical reanalysis of these widely used models, detailing their assumptions and clarifying some steps in their derivation that were previously left unsaid. Our findings are implemented in the pAGN module for the Python programming language, which is the first public implementation of these accretion-disc models. We further apply pAGN to the evolution of stellar-mass black holes embedded in AGN discs, addressing the potential occurrence of migration traps.


INTRODUCTION
Active galactic nuclei (AGNs) are compact regions at the center of galaxies powered by gas accretion onto supermassive black holes (BHs) as opposed to solely the radiation from stars.The underlying theory describing the accretion disc of AGNs was first introduced by Zel'dovich (1964) and Salpeter (1964).AGNs have been studied at low and high redshifts across several electromagnetic wavelengths, capturing a wide range of astrophysical phenomena (see Netzer 2015;Padovani et al. 2017;Hickox & Alexander 2018;Bianchi et al. 2022 for broad reviews on the topic).Due to the deep gravitational well surrounding the central BH, the gas in the accretion disc is expected to reach temperatures of ∼10 5 K and surface densities of ∼10 5 g cm −2 .The accretion disc of the central BH extends to sub-pc scales and is surrounded by optically thick material which is coupled to the disc itself.These components are collectively referred to as the AGN disc, which is expected to extend to separations of 1-10 pc (Netzer 2015).Because of high obscurations and uncertainty in observations, the actual size of AGN discs is somewhat unclear but tends to be larger than what is expected from theoretical models (Jha et al. 2022;Guo et al. 2022a,b).
AGN discs are unique astrophysical environments with a rich phenomenology, including high-energy jets, dusty torii, and accreting BHs (Padovani et al. 2017).In the context of gravitational-wave observations, AGN discs are studied as host environments for compact-binary formation and mergers (McKernan et al. 2012(McKernan et al. , 2011;;Yang et al. 2019;Secunda et al. 2019;Fabj et al. 2020;Tagawa et al. 2020;Trani et al. 2024).The large escape velocity around a supermassive BH implies that objects are likely to be retained in the disc vicinities, potentially forming a large population of stellar-mass BHs that have a higher likelihood of interacting.The dense gas in the disc can facilitate binary formation, accelerate the inspiral, and induce chains of hierarchical BH mergers (Gerosa & Fishbach 2021;Santini et al. 2023;Whitehead et al. 2023;Vaccaro et al. 2023).The occurrence of hierarchical mergers in AGN discs crucially depends on the presence of the socalled migration traps, namely locations in the disc where the migration torque changes sign, which is still an open issue in AGN-disc modeling (Bellovary et al. 2016;Tagawa et al. 2020;Grishin et al. 2024).
Early models of AGNs discs consist of one-dimensional, steady-state, semi-analytic solutions utilizing parametric prescriptions.Subsequent computational advancements allowed for models capturing more complex physics, such as radiative transfer, gas phase transitions, magnetic fields, and general relativity (e.g.Wilson 1972;Hawley et al. 1984;Falle 1991;Schartmann et al. 2005Schartmann et al. , 2008;;Wada et al. 2009;Huško & Lacey 2023).Nevertheless, one-dimensional models remain highly valuable today due to their computational efficiency and insightful perspectives on the structure of AGN discs.This makes them particularly useful in the study of interactions between compact objects and BHs.The first of these one-dimensional approaches dates back to Shakura & Sunyaev (1973), who first model geometrically thin, optically thick discs around a BHs.Building on this seminal work, two models are most commonly used in the field, namely those by Sirko & Goodman (2003) and Thompson et al. (2005) (but see also Goodman 2003;Levin 2003;Haiman et al. 2009;Cantiello et al. 2021;Gilbaum & Stone 2022;Grishin et al. 2024;Hopkins et al. 2024a,b).Both these models assume some heating mechanisms in the disc that marginally support the outer regions from collapsing due to self-gravity and formulate one-dimensional sets of equations for the AGN-disc profile as a function of a number of parameters such as the mass of the central BH and the accretion rate.
The Sirko & Goodman (2003) and Thompson et al. (2005) models are widely used and underpin some of the key, qualitative results in the field of AGN-disc physics.Despite that, the underlying parameters and methods are often left unspecified.Achieving a stable numerical implementation of these disc solutions is not straightforward and codes in this area have not been released in the public domain.The goal of this paper is to critically re-analyze the AGN disc models by Sirko & Goodman (2003) and Thompson et al. (2005).In particular, we clarify the model equations one needs to solve (and crucially the order one needs to solve them), highlight the choices one has to make to obtain stable solutions, and provide a highly customizable implementation.Our software is made publicly available in a Python package called pAGN (short for "parametric AGNs", pronounced as "pagan").
This paper is organized as follows.In Sec. 2, we lay out the equations for the Sirko & Goodman (2003) and Thompson et al. (2005) models.In Sec. 3, we explore some of the input parameter space for both models.In Sec. 4, we showcase our implementation, looking in particular at the occurrence of migration traps in either of the two disc models.In Sec. 5, we present the public code pAGN.In Sec. 6, we draw our conclusions and present prospects for future work.

AGN DISC MODELS
We first summarize the AGN disc models by Sirko & Goodman (2003) and Thompson et al. (2005).We refer to the models as SG03 and TQM05, respectively.Both models consist of an inner, thin accretion disc extended to larger radii to explain observed AGN luminosities.In the outer regions, the disc needs to remain marginally stable against fragmentation.With respect to the Shakura & Sunyaev (1973) thin-disc solution, the SG03 model additionally assumes the existence of some heating mechanism generating radiation pressure that can support the outer parts of the disc against collapse.The TQM05 model further modifies the SG03 model, with the most notable change being that the mass advection is driven by non-local torques rather than local viscous stresses.Furthermore, the Sirko & Goodman (2003) accretion rate is constant across the disc while that of Thompson et al. (2005) Table 1.Key parameters entering the Sirko & Goodman (2003) and Thompson et al. (2005) AGN-disc models.The third column indicates whether the parameter is an input of the model (I), a fixed value for the entire disc (F), or a profile parameter obtained by running the model (P).The accretion rate Ṁ is a fixed parameter for the SG03 disc but a profile parameter for the TQM05 disc.varies because it directly takes into account the mass lost to star formation.

Symbol
We now introduce each model in closer detail and present the key equations one needs to solve to build the resulting disc profiles.For clarity, the parameters entering each model are reported in Table 1.A step-by-step guide on how the equations are solved is provided in Figs. 1 and 2.
Calculate Ω S to find T eff using Eq. ( 8).
Guess values of T eff , T .
Guess values of T eff , T .

Modeling strategy
In the inner regions, the SG03 model assumes a thin and viscous accretion disc to be the source of AGN luminosity (as proposed by Pringle (1981)), similar to the disc model by Shakura & Sunyaev (1973).Such a self-gravitating disc cannot be extended to large radii, where the gravitational pull in the vertical direction causes disc fragmentation and star formation, thus depleting the disc of gas to sufficiently fuel the inner regions.The SG03 model resolves this by assuming that some auxiliary heating (i.e.heating that does not come from orbital energy) lowers the density of the gas in the outer region gas, thus reducing the gravitational pressure.This heating is most likely sourced by star formation, but this is left unspecified in the SG03 model.The auxiliary heating process is prescribed so that gas supply from the marginally gravitationally stable outer regions keeps fueling the hotter inner regions all under a constant gas accretion rate Ṁ .
The stability of the disc is encoded by the parameter first defined by Toomre (1964) for circular Keplerian orbits where cs is the speed of sound, ΩS = GM/r 3 is the angular velocity of the disc, Σg = 2ρh is the midplane mass surface density, ρ is the midplane mass density, and h is the height from the midplane.The disc collapses and fragments whenever QS < 1.The SG03 model is made of two regimes.
In the inner region one has QS ≫ 1: the angular frequency and temperature are high and there is no risk of fragmentation.The outer region instead presents QS ∼ 1: the disc is only marginally stable and auxiliary heating sources become necessary to prevent vertical collapse.
The construction of the model proceeds from an inner boundary rmin and assumes a zero-torque boundary condition, see Fig. 1.Sirko & Goodman (2003) approximate the innermost stable circular orbit to be rmin = Rs/4ϵS, where Rs = 2GM/c 2 is the Schwarzschild radius of the BH and ϵS is the radiative efficiency of the BH, which is set to ϵS = 0.1.For each gas ring at a cylindrical radius r from the central BH, one first assumes that the ring is located in the inner regime where QS ≫ 1.The equations presented in Sec.2.1.2below are then solved to find ΩS and ρ.In turn these are used to evaluate QS from Eq. (1).If QS < 1, one switches to the QS = 1 regime and solves the equations from Sec. 2.1.3instead.This process is then repeated for every value of r until r = rmax.Unless specified, we set rmax to the minimum between 10 7 Rs and 1 pc.An unreasonably large value of rmax leads to a spectral energy distribution that does not match observations (cf.Sirko & Goodman 2003).
The accretion rate of the SG03 disc is parameterized by the Eddington ratio where L Edd is the Eddington luminosity and the normalization is set to the luminosity of a non-self gravitating disc.In turn, the Eddington luminosity is where κes = 0.2(1 + X) cm 2 g −1 is the electron scattering opacity for a fractional abundance of hydrogen which we assume to be X = 0.7.The SG03 model thus depends on the mass of the central BH M through both the angular velocity of the disc ΩS and the accretion rate Ṁ .The disc viscosity is prescribed using the Shakura & Sunyaev (1973) where is the fraction of gas pressure pgas to total pressure ptot; the latter contains contributions from both gas and radiation.
The parameter b = {0, 1} acts as a switch flag to determine how viscosity and pressure relate in the disc.The two cases are often referred to as α-disc (b = 0) and β-disc (b = 1), see e.g.Haiman et al. (2009).For the gas pressure, we use the ideal gas law where kB is the Boltzmann constant and mU is the atomicmass constant.The radiation pressure is given by which is constructed such that in the optically thick regime it recovers p rad = 4σSBT 4 /3c, but retains a dependency on τv in the optically thin regime (Sirko & Goodman 2003).The source of the radiation pressure is not made explicit by Sirko & Goodman (2003), but is assumed to come from stellar processes such as supernovae and nuclear fusion in stars.

Inner regime
For each value of r, the model first assumes that there is no star formation (QS ≥ 1).Each annulus is treated as a black-body with an effective temperature T eff .This is found by equating the locally radiated flux to the viscous heating rate per unit area (Shakura & Sunyaev 1973): where σSB is the Stefan-Boltzmann constant and we have defined Ṁ ′ = Ṁ (1 − rmin/r).Equation ( 8) assumes that all material below r = rmin falls into the BH and cannot energetically interact with the rest of the disc.Mass conservation relates the viscosity of the gas ring to the accretion rate (Sirko & Goodman 2003) which gives two families of solutions: b = 0 (where the viscosity is proportional to total pressure) and b = 1 (where the viscosity is proportional to the gas pressure only).The sound speed in the disc is defined as In this regime, for each value of r, we look for solutions in cs and T and rearrange all other parameters as functions of cs and T only.The midplane height can be expressed as a function of cs by assuming hydrostatic equilibrium The value of the density as a function of cs and T is then given by substituting Σg = 2ρh, Eq. ( 10) and Eq. ( 11) into Eq.( 9) for the b = 0 case, and combining them with the equation for the gas pressure [see Eq. ( 6)] for the b = 1 case.The temperature profile in the disc depends on the optical depth τv = κρh, where κ(ρ, T ) is the opacity.The latter is obtained using interpolated values by Semenov et al. (2003) when T < 10 4 K and Badnell et al. (2005) when T > 10 4 K, the set of which we refer to as the "combined" opacity.These are newer prescriptions for the opacity compared to those by Iglesias & Rogers (1996) and Alexander & Ferguson (1994) used by Sirko & Goodman (2003).The opacities in Semenov et al. (2003) are calculated for silicate grains; the effect of graphite that is important at temperatures of ∼2000 K and may be responsible for the broad line region in AGN observations (see Baskin & Laor 2018) is ignored.The inclusion of the effect of graphite in pAGN is left to future work.From the opacity and effective temperature, we look for solutions in T by assuming the disc ring is in radiative equilibrium: where the functional form of the equation was chosen to match the temperature dependence on T eff and τv across both the optically-thick and optically-thin regimes, cf.Sirko & Goodman (2003).Finally, one can look for solutions in cs and T by considering Eq. ( 10).The Toomre stability parameter QS is calculated from the second expression in Eq. ( 1).If this falls below 1, it is assumed that the ring is in the outer regime and a different set of equations is used, which we present next.

Outer regime
In the outer regions, the model expects the disc to be only marginally gravitationally stable, i.e.QS = 1.In this case, Eq. ( 8) no longer applies since there is additional auxiliary heating.The density is then given by which is a rearrangement of Eq. ( 1) where QS = 1.In the inner regime we look for solutions in cs and T ; in the outer regime we know the value of the density ρ and instead look for solutions in T and T eff .This is done by rearranging Eq. ( 9) and substituting in Eqs. ( 10), ( 11) and ( 13) to obtain an expression for cs as a function of T .In the b = 0 case, cs can be independently determined for a given value of r; in the b = 1 case, cs is a function of the temperature T [see Eq. ( 6)].To find values for T and T eff , we look for solutions that satisfy Eqs. ( 10) and ( 12) simultaneously.In fig. 2 in Sirko & Goodman (2003), there are three different solutions for the disc parameters from r ≳ 5×10 5 Rs.Our implementation recovers the same behavior, but we only accept the continuous, high-temperature, low-opacity solution.For this case, the midplane temperature of the disc remains above 10 3 K and the opacity drops to 10 −3 cm 2 g −1 in the outer regime, both of which affect the gas and radiation pressure profiles, as reflected in the parameter β.The transition between the inner and outer disc regimes takes place at r ≈ 10 3 Rs, which is consistent with the original results by Sirko & Goodman (2003).

Disc profiles
Figure 3 also compares the b = 0 case with the b = 1 case for the same AGN disc.The difference between the two is that the viscosity is assumed to be proportional to the total and gas pressure, respectively.The b = 1 case remains in the optically thick regime out to larger separations, thus also maintaining higher temperatures in those outer regions.For this set of input parameters, the aspect ratio of the disc becomes > 1 at separations as small as r ≈ 5 × 10 4 Rs, which breaks the thin-disc assumption.
In the outer regime, the density scales as ρ ∝ Ω 2 S ∝ r −3 [see Eq. ( 13)].Furthermore, the condition r ≫ rmin implies that Ṁ ′ is approximately constant.Depending on the value of b, one can use Eq. ( 9) to relate ρ, T and ΩS.By approximating the system as optically thick (τv ≫ 1), Eq. ( 12) gives T 4 ∝ τvT 4 eff .Using Eq. ( 10), one can then find simple power-law scalings for most parameters in the outer regime, as long as the opacity κ is kept constant.These are shown Fig. 3 for the 10 3 Rs ≲ r ≲ 10 5 Rs region.In particular, one has T ∝ r −3/4 for b = 0 and T ∝ r −1/2 for b = 1.From Eq. ( 10), we find that in the optically thick regime, cs is approximately constant when b = 0 and proportional to r 1/2 when b = 1.At r ≳ 10 5 Rs, one has κ ≪ 1 and both discs fall back to the optically thin regime.In this case, Eq. ( 12) gives T 4 ∝ T 4 eff /τv, which from Eq. ( 10) gives T 4 ∝ ρc 2 s /τ 2 v for both b = 0 and b = 1.If we assume κ to be constant, then T will also remain constant for all values of r.
Figure 3 shows the "most unstable mass" Mmost unst.≡ c 4 s /G 2 Σg at a given radius.This is the mass enclosed in protostellar clumps with a characteristic radius rc = c 2 s /GΣg (Toomre 1964) and corresponds to the maximum mass that can be present in local perturbations and is thus available for star formation.Figure 3 shows that, for both the b = 0 and b = 1 cases, Mmost unst.has a minimum at r ≈ 10 3 Rs, corresponding to high Σg values and low cs values.Below this radius, where Q > 1 and star formation ceases, the value of Mmost unst.no longer provides meaningful information.Thompson et al. (2005) 2.2.1 Modeling strategy Thompson et al. (2005) proposes an AGN model for which the outer areas of the disc are vertically supported against gravitational collapse by radiation pressure from star formation by-products, dominated in the optically thick regime by dust grains around massive stars.The angular momentum transport in the TQM05 disc is described by global torques instead of a local viscosity mechanism like in the

2.2
Radial profile in the SG03 disc model for the opacity κ, the mass density ρ, the pressure fraction β, the optical depth τv, the disc mass M disc , the sound speed cs, the temperature T , the most unstable mass Mmost unst., the Toomre stability parameter Q S , the effective temperature T eff , the surface mass density Σg and the half-thickness h.All input values are the same as in Fig. 2 by Sirko & Goodman (2003): M = 10 8 M ⊙ , α = 0.01, l E = 0.5, ϵ S = 0.1.Blue (orange) curves indicate the case where b = 0 (b = 1) and the viscosity is proportional to total (gas) pressure.
SG03 model, which provides rapid radial advection rates in the outer regions of the disc.
In the TQM05 model, the angular velocity which is taken from their full galaxy sample.We stress that both of these expressions were obtained for surveys of non-AGN galaxies, meaning that they do not appropriately account for selection biases (Barausse et al. 2017;Shankar et al. 2017;Menci et al. 2023).
The TQM05 model accounts for the star-formation rate per unit area: which is parametrized using the fraction η of the disc dynamical timescale.By means of Σ * , the TQM05 model explicitly tracks changes in the accretion rate Ṁ throughout the disc due to star formation.The gas accreted onto the central BH is supplied by material outside of a radius rout at a constant rate Ṁout.As Thompson et al. (2005) point out, the AGN disc for the TQM05 model does not have a clear outer boundary because the gas is expected to be fed to the central BH by the surrounding interstellar medium.Unlike rmax in the SG03 model, which is a chosen value after which the gas is expected to fragment into stars, here rout represents a transition point beyond which the accretion rate is constant and within which the accretion rate varies due to star formation.Opposite to the SG03 case, in the TQM05 model one integrates from the outer boundary of the AGN disc rout down to the inner edge of the disc, here set to rmin = 3Rs.
In the TQM05 model, the Toomre (1964) stability criterion is written as where κ 2 Ω = 4Ω 2 T + dΩ 2 T /d ln r is the epicyclic frequency.To first order in 1/r, Eq. ( 14) gives dΩT/dr ≈ −ΩT /r such that κΩ ≈ √ 2Ω.When QT ≫ 1, we expect conditions to be unfavorable to star formation so that Σ * and η are close to zero.In the outer area of the disc where QT ≈ 1, stellar feedback plays a key role in stabilizing the disc.
Much like the SG03 model, 1 the TQM05 one also has two regimes according to the value of QT, see Fig. 2. We initialize our numerical root finder at the outer boundary assuming that the disc is optically thick to its own infrared radiation and that QT = 1, thus obtaining initial values for T , ρ and η (see Appendix A).

Non star-forming regime
For every value of r under consideration, we first assume that there is no star formation and that QT > 1.In this case, the accretion rate is constant and thus the value of Ṁ is the same as that of the preceding separation, i.e.Ṁ (r) = Ṁ (r + ∆r), where ∆r is the numerical radial resolution.At r = rout, one has the boundary condition Ṁ (rout) = Ṁout.The gas ring at cylindrical radius r is assumed to radiate as a black body with effective temperature: which is the same as Eq. ( 8).The TQM05 model assumes that the angular momentum in the disc is transported by global torques, so that the radial velocity of the gas vr is a fraction mT of the sound speed cs.The resulting accretion rate is 1 For small values of r one has that Ω T is approximately Keplerian and Q T ≈ Q S .Since Q T is expected to be ≫ 1 near the BH, the factor √ 2 is negligible.
where we have assumed hydrostatic equilibrium, h = cs/ΩT [see Eq. ( 11)].Using Eq. ( 19), one can compute the disc half thickness h as a function of the accretion rate and density.We then interpolate the opacity tables of our choosing to find the κ(ρ, T ), which in turn gives us the optical depth τv = κρh as a function of T and ρ.Notably, Thompson et al. (2005) use the opacities by Semenov et al. (2003) which are provided for temperatures up to T ≃ 10 4 K and extrapolate them to higher temperatures by keeping κ(ρ, T ) constant; in the following we refer to this set as the "Semenov" opacities.
In pAGN, we instead use the combined set of opacities with values by Semenov et al. (2003) up to T = 10 4 K and then values by Badnell et al. (2005) for higher temperatures.
We look for solutions in T and ρ so that the gas ring is in radiative equilibrium and the sound speed is consistently defined.The condition for radiative equilibrium adopted by Thompson et al. (2005) is which is the same as Eq. ( 12) but doubled.The definition of the sound speed cs = ptot/ρ is almost identical to that given in Eq. ( 10) for the SG03 model.The sound speed definition assumes hydrostatic equilibrium and the pressure definitions pgas = ρkBT /mU, p rad = σSBτvT 4 eff /c.The additional factor of 2 in the TQM05 model's definition of p rad ensures that in the optically thick regime using Eq. ( 20) gives p rad ≈ 4σSBT 4 /3c.
Solutions for ρ and T are then found by balancing Eqs. ( 10) and (20).One can then compute QT once more using Eq. ( 17).If QT < 1, the ring at radius r is instead assumed to be in the outer star-forming regime.

Star-forming regime
In the case where there is star formation, the accretion rate is no longer constant.Instead, it is calculated numerically by taking the difference between the initial Ṁout and the integrated accretion rate from star formation down to the current ring: where the subscript j denotes that the given parameter is taken at r = rj.Like in the SG03 model, we assume marginal stability, i.e.QT = 1.Rearranging Eq. ( 17) for the mass density yields The two parameters we are solving for in this case are the temperature T and the star formation fraction η of the disc ring.One calculates h from Eq. ( 19), interpolates the value of κ(T ) and finds τv(T ), finally calculating T 4 eff (T ) using Eq. ( 20).
We now look for solutions in η and T that balance the radiated flux which now directly accounts for radiation from stars unlike Eqs. ( 8) and ( 18), while assuming hydrostatic equilibrium.One has where ϵT and ξ are free parameters describing the efficiency of star formation in the disc and the radiative efficiency of supernovae, respectively.In this regime, it is expected that the gas will be optically thin, and therefore radiation pressure from supernovae is included through the ξ parameter.Thompson et al. (2005) sets ϵT = 10 −3 and ξ = 1.We seek the values of η and T that simultaneously solve Eqs. ( 24) and (25).

Accretion criterion
Unlike Sirko & Goodman (2003), the Thompson et al. (2005) model presents an accretion rate Ṁ that changes as a function of the radial separation r.This naturally means that if the accretion rate at the outer boundary is too low, not enough gas is able to reach the central BH to maintain high temperatures and bright AGN luminosities, which are expected to be in the 10 −3 − 0.5 L Edd range, see Heckman et al. (2004); Kollmeier et al. (2006); Suh et al. (2015); Kong & Ho (2018).This introduces a minimum threshold for Ṁout.Thompson et al. (2005) argue that accretion rates of ∼1 − 10 M⊙yr −1 at the inner disc boundary rmin are sufficient to produce a bright AGN when the central BH mass is ∼10 9 M⊙.This is equivalent to a minimum BH accretion rate of Ṁ ∼ 0.2 ṀEdd = 0.2 L Edd /(0.1c 2 ) at r = rmin.Using Eq. ( 47) in Thompson et al. (2005), we find that over a wavelength range of [10 −8 m, 10 −3 m], setting 0.2 ṀEdd gives a disc bolometric luminosity of 2 × 10 −4 L Edd .
There is no general expression that relates the accretion rate Ṁout and outer radius rout to the BH mass M that would ensure a bright AGN disc.Nonetheless, we can attempt to find such a relationship by considering how the accretion rate at the outer boundary Ṁout scales with the size of the disc rout and the central BH mass.Thompson et al. (2005) proposes a critical value Ṁc, obtained by equating the star formation timescale τ * = 1/ηΩ with the advection timescale τ adv = r/vr to determine whether enough material reaches the central BH to form a luminous signal.From Eq. ( 19) we find Ṁc = 4πr 2 ρhηΩT . (26) Together with Eq. ( 26), this result can be used to introduce a dependence on rout and M to Ṁout.A BH with M = 10 8 M⊙ surrounded by a TQM05 disc that has rout = 95 pc, Ṁout = 320 M⊙yr −1 , and σ = 188 km s −1 satisfies Ṁout > Ṁc(r = rout) and has an accretion rate near the central BH of ≈1.93 M⊙yr −1 = 0.74 ṀEdd (giving a disc luminosity of 9.61 × 10 −4 L Edd ).We use these values to keep the ratio of Ṁout/ Ṁc constant.Using Eqs. ( 16), ( 14), ( 26), and assuming the optically thick regime (see Appendix A), one can show that Ṁc ∝ rσ 2 .Therefore, we scale Ṁout with the outer boundary of the disc and the dispersion velocity, i.e.

Ṁout = 320 M⊙yr
However, for high masses, Eq. ( 27) is not enough to fulfill the the Ṁout > Ṁc criterion.The inability of Ṁc to accurately predict whether a bright AGN disc is formed is not surprising, as it compares the timescales for only one value of r.In Thompson et al. (2005), it is stated that discs with Ṁout < Ṁc(r = rout) cannot form bright AGNs.We find that using Ṁc from Eq. ( 26) as a threshold is too stringent and often omits signals that produce bright AGN discs.In the following, we use Eq. ( 27) as an initial guess for Ṁout but then make adjustments if the accretion rate is not large enough to form a luminous AGN.Developing a full prescription is left to future work.As a precaution to avoid setting an Ṁout that is too high, TQM05 suggests a maximum limit for Ṁout equal to Ṁmax = 8πρhσ 2 r/ϵTc = LM/ϵTc 2 , where LM is the limiting Eddington-like luminosity below which a galaxy will not have momentum driven winds that are high enough to significantly reduce the gas in the disc (Murray et al. 2005).
Other authors have used different values for Ṁout and rout.For instance, Stone et al. (2017) scale the AGN disc down to a Milky-Way type galaxy, using M = 3 × 10 6 M⊙, Ṁout = 15 ṀEdd , and rout = 10 pc, where the latter was motivated by the radius of the dusty tori from AGN disc observations (Jaffe et al. 2004;Burtscher et al. 2013;García-Burillo et al. 2019;Sajina et al. 2022).Our implementation results in disc profiles that diverge when the disc is close to the central BH, around r ∼ 10 −3 pc in this case; this follows from the r → rmin limit in the definition of Ṁ ′ .We report good agreement between the two opacity implementations, with a noteworthy difference being the presence of the iron opacity bump (Jiang et al. 2016) at a radius of r ∼ 2 × 10 −4 pc, seen only for the combined opacities.For both sets of opacities, the disc profile presents a sharp feature at r ∼ 5 × 10 −1 pc where the temperature becomes high enough to leave the so-called opacity gap (the dip in κ for temperatures 10 3 K ≲ T ≲ 10 4 K, see Semenov et al. 2003;Thompson et al. 2005).Figure 4 shows that for this set of parameters the disc profile is not sensitive to the choice of opacity tables.

PARAMETER-SPACE EXPLORATION
We now present a brief exploration of the phenomenology predicted by the Sirko & Goodman (2003) and Thompson et al. (2005) disc models.

Mass dependency
We first investigate the behavior of both models as a function of the mass of the central BH. Figure 5  Radial profile in TQM05 disc model for the temperature T , the effective temperature T eff , the mass density ρ, the surface mass density Σg, the gas pressure pgas, the radiation pressure p rad , the gas pressure fraction β, the half-thickness of the disc h, the sound speed cs, the opacity κ and the optical depth τv.The input values have been chosen to reproduce Fig. 6 in Thompson et al. (2005): σ = 300 km/s, ϵ T = 10 −3 , m = 0.2, Ṁout = 320M ⊙ yr −1 , and rout = 200 pc.Models shown in blue use the opacities by Semenov et al. (2003), models shown in orange use the combined datasets from Semenov et al. (2003) and Badnell et al. (2005).and TQM05 discs profiles of four output parameters, namely the disc height from the midplane h, the mass density ρ, the optical depth τv, and the temperature T , for three central BH masses: M = 10 6 , 10 8 , and 10 10 M⊙.These five output quantities can be used to fully reconstruct an AGN disc for both models.Results are presented using the combined opacity datasets.
For the Sirko & Goodman (2003) model, we set α = 0.01, lE = 0.5, and only consider the α disc (i.e.b = 0).For each disc, we find the solution up to a radius of 10 7 Rs, with the M = 10 6 M⊙ case having a maximum extension of ∼1 pc, and the M = 10 10 M⊙ case extending to ∼1 kpc.
The temperature of the SG03 disc is higher at small separations for lower masses.In particular, one has r ∝ Rs ∝ M in Fig. 5, so that ΩS ∝ M −1/2 , and thus T ∝ M −3 in the inner region, cf.Eq. ( 12) for the optically-thick regime in the SG03 model.In the outer regions of the SG03 disc, all three models have the same temperature T ≈ 7.5 × 10 3 K, which is reached at the separation where the disc becomes optically thin (τv < 1).At large radii, if the disc is dominated by radiation pressure and the gas is optically thin (T 4 eff ∝ τvT 4 ), then from Eq. ( 10) we find that c 2 s ∝ τ 2 v T 4 /ρ.If κ is independent of r, then τv ∝ ρh, which in hydrostatic equilibrium gives a constant T independent of both r and M .
Figure 5 shows that the density ρ is lowest when the central BH mass is highest, with ρ ∝ M 2 in the inner region of the SG03 disc.The model with M = 10 10 M⊙ presents the thickest SG03 disc, reaching h/r > 1 at r ≳ 10 6 Rs; this is outside the regime of validity of our equations but only applies for large radii suggesting a diffuse envelope of gas around the AGN disc.
The Thompson et al. (2005) model shown in Fig. 5 uses mT = 0.2, ϵT = 10 −3 and ξ = 1.In Fig. 5, we linearly scale the outer boundary of the disc rout using the Schwarzschild radius so that rout = 10 7 Rs for all three BH masses.We calculate Ṁout using Eq. ( 27) for all three BH masses, but find that for the M = 10 10 M⊙ case the scaled Ṁout does not satisfy the Ṁout > Ṁc condition and the disc profile looks significantly different from the AGN discs with smaller masses (the height ratio h/r monotonically decreases and the temperature in the disc does not reach 10 4 K).Instead, we opt for Ṁout = 1.5×10 6 M⊙yr −1 when M = 10 10 M⊙ instead.Equation ( 27) gives Ṁout = 0.37 M⊙yr −1 when M = 10 6 M⊙ and Ṁout = 322 M⊙yr −1 when M = 10 8 M⊙.The AGN disc with M = 10 8 M⊙ has an outer boundary of 100 pc, which is about half the size of the model shown in Fig. 4.
The M = 10 6 M⊙ case in Fig. 5 shows an AGN disc with an outer boundary rout ≈ 1 pc and a BH accretion rate 0.37 M⊙yr −1 .Its accretion rate Ṁ is higher than both the star formation rate and Ṁc for all values of r, leading to temperatures as large as T ∼ 10 6 K at r = rmin and a disc luminosity of 2 × 10 −5 L Edd .The radiation pressure in such a high-temperature region leads to a thick disc, with h/r > 1 below r ∼ 5 × 10 2 Rs.At this aspect ratio, the thin-disc approximation no longer applies and caution must be applied when interpreting our results.In order to reduce h/r in the inner regime, one can decrease Ṁout or decrease mT.
We find that the model with M = 10 10 M⊙ also reaches h/r > 1 but at r > 10 5 Rs.This is due to a combination of low densities, a large optical depth, and a large accretion rate which all increase the radiation pressure at the outer boundary.The M = 10 10 M⊙ AGN disc extends out to 10 kpc and has an accretion rate of ∼10 M⊙yr −1 = 0.04 ṀEdd at r = rmin, giving a disc luminosity of 0.07 L Edd .For the TQM05 model with M = 10 10 M⊙, the optical depth τV shows oscillations at r ∼ 200Rs (see Fig. 5) which are due to the model switching between the inner and outer regimes back and forth when close to the QT = 1 boundary.

Input parameters
The SG03 model has five input parameters: the mass of the central BH M , the luminosity ratio lE (or alternatively the accretion rate Ṁ ), the disc viscosity α, the BH radiative efficiency ϵS, and the pressure flag b = 0, 1.We consider a fiducial model with M = 10 8 M⊙ ϵS = 0.1, α = 0.01, lE = 0.5 and b = 0. Of these parameters, Fig. 6 explores the effect of varying α and lE.
The density ρ in the outer regime is largely independent of α and lE.The Shakura & Sunyaev (1973) parameter α relates the viscosity to pressure and accretion, cf.Eq. ( 9).A larger α in the Sirko & Goodman (2003) model implies a lower density and lower temperature in the inner regime, cf.Fig. 6.In the outer regions, the density is independent of the viscosity and thus independent of α.
We vary the Eddington ratio from lE = 10 −3 to lE = 1, capturing the range of observed AGNs (Heckman et al. 2004;Kollmeier et al. 2006;Suh et al. 2015;Kong & Ho 2018).The Eddington ratio parameterizes the accretion rate, which plays a key role in the disc dynamics at all radial distances from the BH.Scaling relations in the optically thick regime far from the disc (see Sec. 2.1.4)indicate that the SG03 model maintains a constant temperature and density at r ≳ 10 5 Rs.
Higher accretion rates leads to higher effective temperatures [Eq.( 8)], higher disc temperatures overall [Eq.( 12)], and higher total pressure in the disc [Eq.( 9)], which also implies that h must be higher to maintain hydrostatic equilibrium.
The TQM05 model has six input parameters: the mass M of the SMBH from which we get the velocity dispersion σ using Eq. ( 15), the star formation efficiency ϵT, the efficiency of angular momentum transport mT in the disc, the supernovae radiative fraction ξ, the outer boundary of the disc rout, and the accretion rate at this outer boundary Ṁout.Figure 7 assumes a fiducial model with M = 10 8 M⊙, rout = 10 7 Rs, ϵT = 10 −3 , ξ = 1, mT = 0.2 and Ṁout ≈ 312 M⊙yr −1 from Eq. ( 27).Starting from this set of parameters, we explore how the disc profile changes when varying either Ṁout or mT.
We consider three values of the accretion rate: Ṁout = 15, 100, 300 M⊙yr −1 .The lowest accretion rate considered, Ṁout = 15 M⊙yr −1 , falls below the critical accretion rate Ṁc ≈ 21 M⊙yr −1 from Eq. ( 26) at r = rout.According to this criterion, this model should not produce an AGN that is sufficiently bright.At r = rmin, the accretion rate for the Ṁout = 15 M⊙yr −1 case is ∼ 0.58 M⊙yr −1 = 0.22 ṀEdd , which is below the 1 − 10 M⊙yr −1 threshold indicated by Thompson et al. (2005).For this case, the disc luminosity is 1.8 × 10 −4 L Edd , which still falls in the range of Eddington ratios one might expect for AGN discs.This further shows that Ṁc is too strict a criterion for determining whether a TQM05 disc forms an AGN.The disc with such a low accretion rate has a different structure compared to the other two cases, with temperatures that are typically lower.As illustrated in Fig. 7, these low temperatures lead to low radiation pressure that fails to effectively counteract the vertical collapse of the disc and thus lower h/r values.On the other hand, for cases where the outer accretion rates clear the Ṁc criterion, we find that the profiles become identical when in the inner, non-star forming regime, see the region left of the QT = 1 line in Fig. (7).For these cases, the advection timescale and star formation timescale reach an equilibrium at the opacity gap (τ adv = τ * when T ≈ 10 3 K).This leads to discs of the same temperature, density, aspect ratio and accretion rate ( Ṁ = 2.23 M⊙yr −1 at r = rmin for both discs).
The global torque efficiency parameter mT is strongly correlated to the behaviour of the disc for all radial distances.Much like α for the SG03 model, here mT parametrizes the relationship between the angular momentum transport and the accretion rate, cf.Eq. ( 19).In the outermost regions of the disc, where the accretion rate Ṁ ≈ Ṁout is roughly constant and h/r ∼ 1, the total pressure is inversely proportional to mT, see Eq. ( 19) and the definition of cs.This gives a thinner TQM05 disc, with lower values of h/r for higher values of mT at the outer boundary.The density is constant in the marginally stable outer region because of Eq. ( 17), but the low total pressure causes some temperature deviations at r ≈ 10 7 Rs for each disc we consider.These variations contribute to different initial conditions in τv for each value of mT.The TQM05 disc has similar behavior for all three mT values once the solutions enter the opacity gap at r ≈ 10 5 Rs, though differences in the optical depth impact the disc profiles at small values of r.In the innermost regions of the disc, we find that high mT values lead to thick, low density discs due to low radiation pressure (which is proportional to τv by definition).Model variations for the SG03 model, showing in particular the aspect ratio h, the midplane mass density ρ, the optical depth τv and the midplane temperature T .For both columns, we set M = 10 8 M ⊙ and b = 0.In the left column, we consider AGN discs with an Eddington fraction l E = 0.5 and vary the viscosity with α = 0.01 (blue), α = 0.05 (orange) and α = 0.1.In the right column, we consider AGN discs where α = 0.01, and vary the Eddington ratio l E = 0.001 (blue), l E = 0.01 (orange), l E = 0.1 (green) and l E = 1 (red).For each disc instance, the radius at which Q S = 1 is marked by a vertical line.
migration.Migration traps are an effective mechanism for merging stellar-mass BH binaries embedded in AGN discs, especially in a hierarchical manner (McKernan et al. 2012;Bellovary et al. 2016;Yang et al. 2019;Santini et al. 2023;Vaccaro et al. 2023, see Gerosa & Fishbach 2021 for a review).Earlier works by Bellovary et al. (2016) and Grishin et al. (2024) showed that the the location of migration traps does not depend on the properties of the migrating object.The location of these migration traps turns out to be a non-trivial function of the AGN disc parameters, ultimately set by the complex interplay of the gradients of the surface density Σg and temperature T .Migration traps are thus an ideal context to showcase our implementation of the Sirko & Goodman (2003) and Thompson et al. (2005) disc models.Table 2 summarizes all parameters used for this section.

Torque implementation
In particular, we apply our AGN disc models to the methods by Grishin et al. (2024), adopting their migration torque and thermal torque expressions.Grishin et al. (2024)  report the existence of migration traps.However, migration traps disappear when considering the updated migration torque formulas by Jiménez & Masset (2017).Grishin et al. (2024) then add a new type of migratory torque, namely the thermal torque by Masset (2017), and find that migration traps are able to form in their AGN disc model once more.We apply the same methodology and formulas to our more complex AGN models.
Migration induces two over-dense spiral arms in the disc.Each arm will produce a torque acting on the migrating object with a magnitude (Korycansky & Pollack 1993 where q ≡ mBH/M is the BH mass ratio and Ω is equal to either ΩS or ΩT depending on the AGN disc model.The net torque ΓI acting on the migrator in a locally isothermal limit is given by (Paardekooper et al. 2010 where γ = 5/3 is the adiabatic index.The parameter is the Lindblad torque where is a function that adds a dependence on the thermal diffusivity for the Lindblad torque and can be approximated to 1/γ in the case where the diffusivity is small (Masset & Casoli 2010).The thermal diffusivity of the disc is defined as The thermal torque Γ therm originates from the temperature build-up around the migrating object due to lack of heat release during its orbital evolution.If heat is trapped around the migrator, two cold and dense lobes are formed in the disc, which leads to inward migration (Lega et al. 2014).If the migrator is instead able to release heat back into the disc around it, two hot and under-dense lobes form, leading to outward migration (Benítez-Llambay et al. 2015).The total heating torque is (Masset 2017) where xc is the corotation radius of the migrating object, λ is the typical size of the lobes, and L is the luminosity generated by the migrator through thermal heating, and is the critical luminosity.If L = Lc, the hot and cold torques acting on the migrator balance out and Γ therm = 0. We approximate the luminosity of the migrator, L, to be its Eddington luminosity [see Eq. ( 3), replacing the mass M with the mass of the migrator mBH].The size of the lobes λ is given by (Grishin et al. 2024) and the corotation radius is (Grishin et al. 2024): We approximate d ln ptot/d ln r by combining the equation for vertical hydrodynamical equilibrium ptot ≈ ρh 2 Ω 2 with the definition of the sound speed c 2 s = h 2 Ω 2 , resulting in dptot/dr ≈ ρc 2 s /r.The thermal torque given by Eq. ( 34) is expected to diminish in optically thin discs.Following Grishin et al. (2024), we multiply Eq. ( 34) by a factor of 1 − exp{−λτv/h}.Additionally, when the mass of the migrator exceeds the thermal mass m th , the thermal torque will be reduced (Guilera et al. 2021).The thermal mass is defined by: where RB is half the Bondi radius In the regions where h < RB we use the disk height h in place of half the Bondi radius RB.To correct for the critical thermal mass, we split Eq. 34 into its heating component (the positive L/Lc term) and its cooling component (the negative term), which we label as Γ therm, hot and Γ therm, cold respectively.The total thermal torque is described by Eq. 34 unless µ th ≡ m th /mBH < 1.In regions of the disc where µ th < 1, the thermal torque is instead given by: which is an approximation of numerical fits detailed in Velasco Romero & Masset (2020) and used in Guilera et al. (2021); Grishin et al. (2024).

Migration traps
Figure 8 shows the migration-torque profiles for a M = 10 6 M⊙ SMBH in both the Sirko & Goodman (2003) and Thompson et al. (2005) model.We use ϵS = 0.1, α = 0.01, lE = 0.5 and b = 0 for the SG03 AGN disc and rout = 10 7 Rs, ϵT = 10 −3 , ξ = 1, mT = 0.2, and Ṁout = 1.5 × 10 −2 M⊙ yr −1 for the TQM05 model.The outer accretion rate Ṁout was set smaller than the value given by Eq. ( 27) in order to enforce h/r < 1 throughtout the disc, unlike the small mass case in Fig. 5. Identifying a migration trap corresponds to regions of the disc where the net migration torque is zero and goes from negative (i.e.inward migration) to positive (i.e.outward migration) as r increases.
The top panel shows the migration torque using both Eq. ( 29) by Paardekooper et al. (2010) and Eq. ( 30) by Jiménez & Masset (2017).When using the former, we find migration traps at r ≈ 22Rs and r ≈ 10 3 Rs for the SG03 model, which is in line with the results reported by both Bellovary et al. (2016) and Grishin et al. (2024).When using the updated migration torque values by Jiménez & Masset (2017) for the SG03 model, we find that the migration torque is always negative and thus the migrator moves across the disc without being trapped.This result is in agreement with those by Grishin et al. (2024).Once the thermal torque from Eq. ( 34) is added to the updated migration torque of Eq. ( 30), the bottom panel in Fig. 8 shows that we again obtain migration traps.In the SG03 AGN disc, we find two migration traps for a M = 10 6 M⊙ central BH and a 10M⊙ migrator occurring at r ≈ 1.4 × 10 3 Rs = 1.4 × 10 −4 pc and r ≈ 6.8 × 10 4 Rs = 6.5 × 10 −3 pc.
When considering the Thompson et al. (2005) model, we obtain a larger number of migration traps, irrespective of the torque prescriptions adopted and the inclusion of the thermal torque contribution to the net torque.For the migration torque by Jiménez & Masset (2017) (the top panel in Fig. 8), we find migration traps form in the TQM05 disc when the gradient d ln Σg/d ln r discretely changes values, as can be seen in the lower panels of Fig. 8 at r ≈ 2.5 × 10 3 Rs and r ≈ 8.3 × 10 3 Rs.When both migration and thermal torques are considered, we find traps at r ≈ 2.5 × 10 3 Rs, r ≈ 8.4 × 10 3 Rs and r ≈ 2.6 × 10 6 Rs for the Thompson et al. (2005) model.
Our implementation of both the SG03 and TQM05 models is released publicly in the pAGN module for the Python programming language.
pAGN is distributed under git version control at github.com/DariaGangardt/pAGN(code repository) The documentation is provided at In addition, the code distributions include opacity tables by Semenov et al. (2003) and Badnell et al. (2005) as well as an interpolation routine.External opacity tables can also be provided by the user.The overall solution strategy follows what is presented in this paper as illustrated in the flowcharts of Figs. 1 and 2.

CONCLUSIONS
This work presents a critical re-analysis of the AGN disc models by Sirko & Goodman (2003) and Thompson et al. (2005).Our findings are implemented in the public pAGN module for the Python programming language (Gangardt & Trani 2024).We presented the equations from the original papers and emphasized their solution strategy.Compared to the original model, our results consider updated opacity tables, relate some of the input parameter (most notably the scaling of the outer accretion rate with the central BH mass for TQM05 case), validate AGN discs through limits on the accretion rate at the disc boundaries, and investigate the limits of the thin-disc approximation.While the parameter exploration presented in this work provides valuable insights, there is room for further enhancement to fully explore the predictions of these models across the entire parameter space.An example of such research, Ballantyne (2008) presented an observation-motivated study of how the TQM05 input parameters affects the properties of AGN discs in Seyfert-like with a particular focus on the "starburst" disc regions with high star formation.
As a further example, in this paper we have applied our pAGN code to the disc-migration problem, reproducing the analysis by Grishin et al. (2024) with the more complex disc profiles by Sirko & Goodman (2003) and Thompson et al. (2005).While we largely confirm previous findings for the SG03 case, our TQM05 disc shows a large number of migration traps, with potential implications for the formation of hierarchical merging stellar-mass BH binaries detectable with current gravitational-wave detectors (Gerosa & Fishbach 2021).This is an interesting avenue for future work.
The AGN disc models by Sirko & Goodman (2003) and Thompson et al. (2005) are widely used in the literature.We hope our full, public implementation of these approaches, together with the details of the underlying evolutionary equations, might facilitate further advances in this area while clarifying their underlying limitations.Both the SG03 and TQM05 models can be applied to various problems and compared to newer AGN-disc modeling approaches.The goal of pAGN is precisely that to aid further research in the growing field of AGN and gravitational-wave science.

Figure 1 .
Figure1.Flowchart showing detailing our solution strategy for the SG03 model.Construction proceeds from the inner disc to the outer disc, with initial guesses on the stability parameter Q Q which are then checked a-posteriori.

Figure 2 .
Figure 2. Flowchart showing detailing our solution strategy for the TQM05 model.Construction proceeds from the outer disc to the inner disc, with initial guesses on the stability parameter Q T which are then checked a-posteriori.

Figure 3
Figure3shows the radial profile of some key disc parameters in the SG03 model, tailored to reproducing fig.2fromSirko 14) is only approximately Keplerian and includes the effect of the dispersion velocity σ.The dispersion and the central mass are related by the M − σ relation from observations.Thompson et al. (2005) used the expression by Tremaine et al. (2002), while we opted for an updated fit by Gültekin et al. (2009) log σ 200 km/s = 1 4.24 log M M⊙ − 8.12 ,

Figure 4
Figure4reproduces Fig.6inThompson et al. (2005), assuming either the Semenov opacities (as was done byThompson et al. 2005) or the combined opacities.The input parameters are σ = 300 km s −1 , M ≈ 10 9 M⊙ [instead of Eq. (15) we use the M -σ relation byTremaine et al. (2002) as was done byThompson et al. (2005)], Ṁout = 320M⊙ yr −1 , rout = 200 pc, m = 0.2, ϵT = 10 −3 , and ξ = 1.Our results shown in Fig.4are generally in agreement with those byThompson et al. (2005).Our implementation results in disc profiles that diverge when the disc is close to the central BH, around r ∼ 10 −3 pc in this case; this follows from the r → rmin limit in the definition of Ṁ ′ .We report good agreement between the two opacity implementations, with a noteworthy difference being the presence of the iron opacity bump(Jiang et al. 2016) at a radius of r ∼ 2 × 10 −4 pc, seen only for the combined opacities.For both sets of opacities, the disc profile presents a sharp feature at r ∼ 5 × 10 −1 pc where the temperature becomes high enough to leave the so-called opacity gap (the dip in κ for temperatures 10 3 K ≲ T ≲ 10 4 K, seeSemenov et al. 2003;Thompson et al. 2005).Figure4shows that for this set of parameters the disc profile is not sensitive to the choice of opacity tables.
Figure4.Radial profile in TQM05 disc model for the temperature T , the effective temperature T eff , the mass density ρ, the surface mass density Σg, the gas pressure pgas, the radiation pressure p rad , the gas pressure fraction β, the half-thickness of the disc h, the sound speed cs, the opacity κ and the optical depth τv.The input values have been chosen to reproduce Fig.6inThompson et al. (2005): σ = 300 km/s, ϵ T = 10 −3 , m = 0.2, Ṁout = 320M ⊙ yr −1 , and rout = 200 pc.Models shown in blue use the opacities bySemenov et al. (2003), models shown in orange use the combined datasets fromSemenov et al. (2003) andBadnell et al. (2005).
Figure6.Model variations for the SG03 model, showing in particular the aspect ratio h, the midplane mass density ρ, the optical depth τv and the midplane temperature T .For both columns, we set M = 10 8 M ⊙ and b = 0.In the left column, we consider AGN discs with an Eddington fraction l E = 0.5 and vary the viscosity with α = 0.01 (blue), α = 0.05 (orange) and α = 0.1.In the right column, we consider AGN discs where α = 0.01, and vary the Eddington ratio l E = 0.001 (blue), l E = 0.01 (orange), l E = 0.1 (green) and l E = 1 (red).For each disc instance, the radius at which Q S = 1 is marked by a vertical line.

Figure 7 .
Figure7.Model variations for the TQM05 model, showing in particular the aspect ratio h, the midplane mass density ρ, the optical depth τv and the midplane temperature.For both columns, we set M = 10 8 M ⊙ , ϵ T = 10 −3 , ξ = 1 and rout = 200 pc.In the left column, we consider AGN discs with a global torque efficiency of m T = 0.2 and vary the accretion rate Ṁout = 15 M ⊙ yr −1 (blue), Ṁout = 100 M ⊙ yr −1 (orange) and Ṁout = 300 M ⊙ yr −1 (green, dashed).The Ṁout = 300 M ⊙ yr −1 case is dashed to show that parameter profiles are identical to the those of the Ṁout = 100 M ⊙ yr −1 case close to the central BH.In the right column, we consider AGN discs with an outer accretion rate Ṁout ≃ 312 M ⊙ yr −1 and vary the global torque efficiency m T = 0.1 (blue), m T = 0.2 (orange), and m = 0.3 (green).For each disc instance, the radius at which Q T = 1 is marked by a vertical line.

Figure 8 .
Figure 8.The absolute values of the migration torques for a m BH = 10 M ⊙ BH orbiting a M = 10 6 M ⊙ central BH in AGN discs.The left panels show torque profiles for a SG03 disc with ϵ S = 0.1, α = 0.01, l E = 0.5 and b = 0.The right panels show torque profiles for a TQM05 disc with ϵ T = 10 −3 , ξ = 1, m = 0.2, rout = 10 7 Rs and Ṁout = 1.5 × 10 −2 M ⊙ yr −1 .Migration torques, thermal torques, and their combination are shown in first three rows from the top, respectively.The bottom panel shows the midplane surface density of the disc for each case.For the Type I migration torques considered in the top row, we show both results using prescriptions by both Paardekooper et al. (2010) (light curves) and Jiménez & Masset (2017) (heavy curves).Colors indicate the sign of the torque, with blue referring to inward migration (i.e.positive torques) and orange referring to outward migration (i.e.negative torques).Vertical grey lines indicate the migration traps for all torque prescriptions except for Paardekooper et al. (2010) in the top panel.
dariagangardt.github.io/pAGN(documentation)together with a set of minimal examples.Our pAGN module is available on the Python Package index.The code can be installed with pip install pagn Packages numpy, scipy, and matplotlib are specified as dependencies.The package is imported with import pagn and contains two main classes for the SG03 and TQM05 implementation, respectively: pagn.SirkoAGN pagn.ThompsonAGN

Table 2 .
Summary of the parameter entering our treatment of disc migration explored in Sec. 4.