pAGN: the one-stop solution for AGN disc modelling

ABSTRACT

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 and Thompson et al. 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.

accretion, accretion discs, black hole physics, galaxies: active

1 INTRODUCTION

Active galactic nuclei (AGNs) are compact regions at the centre 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, Mainieri & Padovani 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⁵ K and surface densities of ∼10⁵ g cm⁻². 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 (Guo et al. 2022a; Guo, Barth & Wang 2022b; Jha et al. 2022).

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. 2011, 2012; Secunda et al. 2019; Yang et al. 2019; Fabj et al. 2020; Tagawa, Haiman & Kocsis 2020; Trani, Quaini & Colpi 2023). 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; Vaccaro et al. 2023; Whitehead et al. 2023). The occurrence of hierarchical mergers in AGN discs crucially depends on the presence of the so-called migration traps, namely locations in the disc where the migration torque changes sign, which is still an open issue in AGN-disc modelling (Bellovary et al. 2016; Tagawa et al. 2020; Grishin, Gilbaum & Stone 2023).

Early models of AGNs discs consist of 1D, 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, Smarr & Wilson 1984; Falle 1991; Schartmann et al. 2005, 2008; Wada, Papadopoulos & Spaans 2009; Huško & Lacey 2023). Nevertheless, 1D 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 1D 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, Quataert & Murray (2005) (but see also Goodman 2003; Levin 2003; Haiman, Kocsis & Menou 2009; Cantiello, Jermyn & Lin 2021; Gilbaum & Stone 2022; Grishin et al. 2023; 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 1D 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-analyse 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 Section 2, we lay out the equations for the Sirko & Goodman (2003) and Thompson et al. (2005) models. In Section 3, we explore some of the input parameter space for both models. In Section 4, we showcase our implementation, looking in particular at the occurrence of migration traps in either of the two disc models. In Section 5, we present the public code pAGN. In Section 6, we draw our conclusions and present prospects for future work.

2 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) varies because it directly takes into account the mass lost to star formation.

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.

Figure 1.

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.

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.

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 |$\dot{M}$| is a fixed parameter for the SG03 disc but a profile parameter for the TQM05 disc.

Symbol	Definition	I/F/P
M	Mass of the central BH	I
Rs	BH Schwarzschild radius	F
LEdd	Eddington luminosity	F
|$\dot{M}_\mathrm{Edd}$|	Eddington accretion rate	F
X	Hydrogen abundance in disc	I
κes	Electron scattering opacity	F
r	Radial distance from the central BH	P
|$\dot{M}$|	Mass accretion rate	F or P
rmin	Inner edge of the disc	I
T	Midplane temperature	P
Teff	Midplane effective temperature	P
ρ	Midplane density	P
h	Height of disc from the midplane	P
Σg	Midplane surface density	P
Σtot	Midplane total dynamical density	P
fg	Gas fraction	P
τv	Midplane optical depth	P
κ	Midplane opacity	P
cs	Midpane sound speed	P
pgas	Gas pressure	P
prad	Radiation pressure	P
Sirko & Goodman (2003) parameters
Symbol	Definition	I/F/P
α	Shakura-Sunyaev viscosity parameter	I
rmax	Outer edge of the disc	I
lE	Disc Eddington ratio	I
ϵS	Radiative efficiency	I
b	Switch for viscosity-pressure relation	I
QS	Toomre stability parameter	P
ν	Disc viscosity	P
ΩS	Rotational velocity	P
β	Gas pressure to total pressure ratio	P
Thompson et al. (2005) parameters
Symbol	Definition	I/F/P
σ	Stellar dispersion velocity	I
rout	Effective outer edge of the disc	I
|$\dot{M}_{\rm out}$|	Accretion rate at rout	I
mT	Global torque efficiency	I
ϵT	Star formation efficiency	I
ξ	Supernova radiative efficiency	I
|$\dot{\Sigma }_\mathrm{*}$|	Star formation rate	P
η	Star formation efficiency fraction	P
QT	Toomre stability parameter	P
ΩT	Rotational velocity	P

2.1 Sirko & Goodman (2003)

2.1.1 Modelling 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 fuelling the hotter inner regions all under a constant gas accretion rate |$\dot{M}$|⁠.

The stability of the disc is encoded by the parameter first defined by Toomre (1964) for circular Keplerian orbits

$$\begin{eqnarray} Q_\mathrm{S} \equiv \frac{c_\mathrm{s} \Omega _\mathrm{S}}{\pi G \Sigma _\mathrm{g}} \approx \frac{\Omega _\mathrm{S}^2}{2 \pi G \rho } \, , \end{eqnarray}$$

(1)

where c_s is the speed of sound, |$\Omega _\mathrm{S} = \sqrt{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 Q_S < 1. The SG03 model is made of two regimes. In the inner region one has Q_S ≫ 1: the angular frequency and temperature are high and there is no risk of fragmentation. The outer region instead presents Q_S ∼ 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 r_min and assumes a zero-torque boundary condition, see Fig. 1. Sirko & Goodman (2003) approximate the innermost stable circular orbit to be r_min = R_s/4ϵ_S, where R_s = 2GM/c² 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 Q_S ≫ 1. The equations presented in Section 2.1.2 below are then solved to find Ω_S and ρ. In turn these are used to evaluate Q_S from equation (1). If Q_S < 1, one switches to the Q_S = 1 regime and solves the equations from Section 2.1.3 instead. This process is then repeated for every value of r until r = r_max. Unless specified, we set r_max to the minimum between 10⁷R_s and 1 pc. An unreasonably large value of r_max leads to a spectral energy distribution that does not match observations (cf. Sirko & Goodman 2003).

The accretion rate of the SG03 disc is parametrized by the Eddington ratio

$$\begin{eqnarray} l_{\rm E} = \frac{\dot{M}\epsilon _{\rm S} c^2 }{ L_{\rm Edd}}\, , \end{eqnarray}$$

(2)

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

$$\begin{eqnarray} L_{\rm Edd} = \frac{4 \pi G M c}{\kappa _{\rm es}} \, , \end{eqnarray}$$

(3)

where |$\kappa _{\rm es} = 0.2(1 + X) \, {\rm cm}^2 {\rm 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 |$\dot{M}$|⁠.

The disc viscosity is prescribed using the Shakura & Sunyaev (1973) dimension-less parameter

$$\begin{eqnarray} \alpha = \frac{\nu }{ c_{\rm s} h \beta ^b} \, , \end{eqnarray}$$

(4)

where

$$\begin{eqnarray} \beta = \frac{p_{\rm gas}}{p_{\rm gas} +p_{\rm rad}} = \frac{p_{\rm gas}}{p_{\rm tot}} \end{eqnarray}$$

(5)

is the fraction of gas pressure p_gas to total pressure p_tot; 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

$$\begin{eqnarray} p_{\rm gas} = \frac{\rho k_{\rm B}T }{ m_\mathrm{U}} , \end{eqnarray}$$

(6)

where k_B is the Boltzmann constant and m_U is the atomic-mass constant. The radiation pressure is given by

$$\begin{eqnarray} p_{\rm rad} = \frac{ \sigma _{\rm SB} \tau _{\rm v}}{2 c}T_{\rm eff}^4\, , \end{eqnarray}$$

(7)

which is constructed such that in the optically thick regime it recovers p_rad = 4σ_SBT⁴/3c, but retains a dependence 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.

2.1.2 Inner regime

For each value of r, the model first assumes that there is no star formation (Q_S ≥ 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):

$$\begin{eqnarray} \sigma _{\rm SB} T_{\rm eff}^4 = \frac{3 \Omega _\mathrm{S}^2}{8 \pi } \dot{M}\left(1-\sqrt{\frac{r_{\rm min}}{r}} \right) = \frac{3\Omega _\mathrm{S}^2}{8\pi } \dot{M}^{\prime } \, , \end{eqnarray}$$

(8)

where σ_SB is the Stefan–Boltzmann constant and we have defined |$\dot{M}^{\prime } = \dot{M}(1-\sqrt{{r_{\rm min}}/{r}})$|⁠. Equation (8) assumes that all material below r = r_min 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)

$$\begin{eqnarray} \beta ^b c_{\rm s}^2 \Sigma _{\rm g} = \frac{\dot{M}^{\prime } \Omega _\mathrm{S} }{3 \pi \alpha } \, , \end{eqnarray}$$

(9)

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

$$\begin{eqnarray} c_\mathrm{s}^2 = \frac{p_{\rm tot}}{\rho } \, . \end{eqnarray}$$

(10)

In this regime, for each value of r, we look for solutions in c_s and T and rearrange all other parameters as functions of c_s and T only. The midplane height can be expressed as a function of c_s by assuming hydrostatic equilibrium

$$\begin{eqnarray} h = \frac{c_\mathrm{s} }{ \Omega _\mathrm{S}} \, . \end{eqnarray}$$

(11)

The value of the density as a function of c_s and T is then given by substituting Σ_g = 2ρh, equation (10) and equation (11) into equation (9) for the b = 0 case, and combining them with the equation for the gas pressure [see equation (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⁴ K and Badnell et al. (2005) when T > 10⁴ 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:

$$\begin{eqnarray} T^4 = \left(\frac{3}{8}\tau _{\rm v}+ \frac{1}{2} + \frac{1}{4\tau _{\rm v}}\right) T_{\rm eff}^4 \, , \end{eqnarray}$$

(12)

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 c_s and T by considering equation (10).

The Toomre stability parameter Q_S is calculated from the second expression in equation (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.

2.1.3 Outer regime

In the outer regions, the model expects the disc to be only marginally gravitationally stable, i.e. Q_S = 1. In this case, equation (8) no longer applies since there is additional auxiliary heating. The density is then given by

$$\begin{eqnarray} \rho = \frac{\Omega _\mathrm{S}^2}{2 \pi G}, \end{eqnarray}$$

(13)

which is a rearrangement of equation (1) where Q_S = 1. In the inner regime we look for solutions in c_s 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 equation (9) and substituting in equations (10), (11), and (13) to obtain an expression for c_s as a function of T. In the b = 0 case, c_s can be independently determined for a given value of r; in the b = 1 case, c_s is a function of the temperature T [see equation (6)]. To find values for T and T_eff, we look for solutions that satisfy equations (10) and (12) simultaneously.

2.1.4 Disc profiles

Fig. 3 shows the radial profile of some key disc parameters in the SG03 model, tailored to reproducing fig. 2 from Sirko & Goodman (2003). In particular, the figure shows a 10⁸ M_⊙ BH surrounded by a disc with α = 0.01, l_E = 0.5, and ϵ_S = 0.1, presenting both the b = 0 and b = 1 cases.

Figure 3.

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 c_s, the temperature T, the most unstable mass |$M_{\rm {most} \, \rm {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⁸ 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.

In fig. 2 in Sirko & Goodman (2003), there are three different solutions for the disc parameters from r ≳ 5 × 10⁵R_s. Our implementation recovers the same behaviour, but we only accept the continuous, high-temperature, low-opacity solution. For this case, the midplane temperature of the disc remains above 10³ K and the opacity drops to 10⁻³ cm²g⁻¹ 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³R_s, which is consistent with the original results by Sirko & Goodman (2003).

Fig. 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⁴R_s, which breaks the thin-disc assumption.

In the outer regime, the density scales as |$\rho \propto \Omega _\mathrm{S}^2 \propto r^{-3}$| [see equation (13)]. Furthermore, the condition r ≫ r_min implies that |$\dot{M}^{\prime }$| is approximately constant. Depending on the value of b, one can use equation (9) to relate ρ, T, and Ω_S. By approximating the system as optically thick (τ_v ≫ 1), equation (12) gives |$T^4 \propto \tau _{\rm v}T_{\rm eff}^4$|⁠. Using equation (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 in Fig. 3 for the 10³R_s ≲ r ≲ 10⁵R_s region. In particular, one has T ∝ r^−3/4 for b = 0 and T ∝ r^−1/2 for b = 1. From equation (10), we find that in the optically thick regime, c_s is approximately constant when b = 0 and proportional to r^1/2 when b = 1. At r ≳ 10⁵R_s, one has κ ≪ 1 and both discs fall back to the optically thin regime. In this case, equation (12) gives |$T^4 \propto T_{\rm eff}^4/\tau _{\rm v}$|⁠, which from equation (10) gives |$T^4 \propto \rho c_\mathrm{s}^2/\tau _{\rm v}^2$| for both b = 0 and b = 1. If we assume κ to be constant, then T will also remain constant for all values of r.

Fig. 3 shows the ‘most unstable mass’ |$M_{\rm most \, unst.} \equiv c_{\rm s}^4/G^2\Sigma _{\rm g}$| at a given radius. This is the mass enclosed in protostellar clumps with a characteristic radius |$r_{\rm c} = c_{\rm s}^2/G\Sigma _{\rm g}$| (Toomre 1964) and corresponds to the maximum mass that can be present in local perturbations and is thus available for star formation. Fig. 3 shows that, for both the b = 0 and b = 1 cases, |$M_{\rm most \, unst.}$| has a minimum at r ≈ 10³R_s, corresponding to high Σ_g values and low c_s values. Below this radius, where Q > 1 and star formation ceases, the value of |$M_{\rm most \, unst.}$| no longer provides meaningful information.

2.2 Thompson et al. (2005)

2.2.1 Modelling 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 SG03 model, which provides rapid radial advection rates in the outer regions of the disc.

In the TQM05 model, the angular velocity

$$\begin{eqnarray} \Omega _{\rm T} = \sqrt{\frac{GM}{r^3} + 2 \frac{\sigma ^2}{r^2}} \end{eqnarray}$$

(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)

$$\begin{eqnarray} \log {\frac{\sigma }{200 \, {\rm km\,s^{-1}}}} = \frac{1}{4.24} \left(\log {\frac{M}{ {\rm M}_{\odot }}}-8.12 \right) \, , \end{eqnarray}$$

(15)

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, Bernardi & Sheth 2017; Menci et al. 2023).

The TQM05 model accounts for the star-formation rate per unit area:

$$\begin{eqnarray} \dot{\Sigma }_* = \Sigma _\mathrm{g} \Omega _\mathrm{T} \eta \, , \end{eqnarray}$$

(16)

which is parametrized using the fraction η of the disc dynamical time-scale. By means of |$\dot{\Sigma }_*$|⁠, the TQM05 model explicitly tracks changes in the accretion rate |$\dot{M}$| throughout the disc due to star formation. The gas accreted onto the central BH is supplied by material outside of a radius r_out at a constant rate |$\dot{M}_{\rm 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 r_max in the SG03 model, which is a chosen value after which the gas is expected to fragment into stars, here r_out 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 r_out down to the inner edge of the disc, here set to r_min = 3R_s.

In the TQM05 model, the Toomre (1964) stability criterion is written as

$$\begin{eqnarray} Q_\mathrm{T} = \frac{\kappa _{\Omega } c_\mathrm{s}}{ \pi G \Sigma _\mathrm{g}} \approx \frac{\Omega _\mathrm{T}^2}{\sqrt{2} \pi G \rho } \, , \end{eqnarray}$$

(17)

where |$\kappa _{\Omega }^2 = 4 \Omega _\mathrm{T}^2 + \mathrm{d}\Omega _\mathrm{T}^2/\mathrm{d}\ln {r}$| is the epicyclic frequency. To first order in 1/r, equation (14) gives dΩ_T/dr ≈ −Ω_T/r such that |$\kappa _{\Omega } \approx \sqrt{2} \Omega$|⁠. When Q_T ≫ 1, we expect conditions to be unfavourable to star formation so that |$\dot{\Sigma }_*$| and η are close to zero. In the outer area of the disc where Q_T ≈ 1, stellar feedback plays a key role in stabilizing the disc.

Much like the SG03 model,¹ the TQM05 one also has two regimes according to the value of Q_T, 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 Q_T = 1, thus obtaining initial values for T, ρ, and η (see Appendix  A).

2.2.2 Non-star-forming regime

For every value of r under consideration, we first assume that there is no star formation and that Q_T > 1. In this case, the accretion rate is constant and thus the value of |$\dot{M}$| is the same as that of the preceding separation, i.e. |$\dot{M}(r) = \dot{M}(r + \Delta r)$|⁠, where Δr is the numerical radial resolution. At r = r_out, one has the boundary condition |$\dot{M}(r_{\rm out}) = \dot{M}_{\rm out}$|⁠. The gas ring at cylindrical radius r is assumed to radiate as a blackbody with effective temperature:

$$\begin{eqnarray} \sigma _{\rm SB} T_{\rm eff}^4 = \frac{3 \Omega _\mathrm{T}^2 }{8\pi } \dot{M}^{\prime } \, , \end{eqnarray}$$

(18)

which is the same as equation (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 v_r is a fraction m_T of the sound speed c_s. The resulting accretion rate is

$$\begin{eqnarray} \dot{M}= 4 \pi r \rho h v_{\rm r} = 4 \pi r \rho h m_\mathrm{T} c_\mathrm{s} = 4 \pi r \Omega _\mathrm{T} m_\mathrm{T} \rho h^2 \, , \end{eqnarray}$$

(19)

where we have assumed hydrostatic equilibrium, h = c_s/Ω_T [see equation (11)]. Using equation (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⁴ 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⁴ 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

$$\begin{eqnarray} T^4 = \left(\frac{3}{4} \tau _{\rm v}+ \frac{1}{2\tau _{\rm v}} + 1 \right) T_{\rm eff}^4 \, , \end{eqnarray}$$

(20)

which is the same as equation (12) but doubled. The definition of the sound speed c_s = p_tot/ρ is almost identical to that given in equation (10) for the SG03 model. The sound speed definition assumes hydrostatic equilibrium and the pressure definitions p_gas = ρk_BT/m_U, |$p_{\rm rad} = \sigma _{\rm SB} \tau _{\rm v}T_{\rm eff}^4/c$|⁠. The additional factor of 2 in the TQM05 model’s definition of p_rad ensures that in the optically thick regime using equation (20) gives p_rad ≈ 4σ_SBT⁴/3c.

Solutions for ρ and T are then found by balancing equations (10) and (20). One can then compute Q_T once more using equation (17). If Q_T < 1, the ring at radius r is instead assumed to be in the outer star-forming regime.

2.2.3 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 |$\dot{M}_{\rm out}$| and the integrated accretion rate from star formation down to the current ring:

$$\begin{eqnarray} \dot{M}(r) &= \dot{M}_{\rm out}-\int _{r}^{r_{\rm out}} 4 \pi r \rho h \Omega _\mathrm{T} \eta \mathrm{d} r \end{eqnarray}$$

(21)

$$\begin{eqnarray} &\approx \dot{M}_{\rm out}-\sum _{r_j = r}^{r_{\rm out}} 4 \pi r_j \rho _j h_j \Omega _{\mathrm{T,} j} \eta _j \Delta r_j \, , \end{eqnarray}$$

(22)

where the subscript j denotes that the given parameter is taken at r = r_j. Like in the SG03 model, we assume marginal stability, i.e. Q_T = 1. Rearranging equation (17) for the mass density yields

$$\begin{eqnarray} \rho = \frac{\Omega _\mathrm{T}^2}{\sqrt{2} \pi G} \, . \end{eqnarray}$$

(23)

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 equation (19), interpolates the value of κ(T) and finds τ_v(T), finally calculating |$T_{\rm eff}^4(T)$| using equation (20).

We now look for solutions in η and T that balance the radiated flux

$$\begin{eqnarray} \sigma _{\rm SB} T_{\rm eff}^4 = \rho h \Omega _\mathrm{T} \eta \epsilon _\mathrm{T} c^2 + \frac{3}{8\pi } \dot{M}^{\prime } \Omega _\mathrm{T}^2 \, , \end{eqnarray}$$

(24)

which now directly accounts for radiation from stars unlike equations (8) and (18), while assuming hydrostatic equilibrium. One has

$$\begin{eqnarray} \rho h^2 \Omega _\mathrm{T}^2 = \frac{\rho k_{\rm B}T}{m_\mathrm{U}} + 2 \rho h \eta \Omega _\mathrm{T} \epsilon _\mathrm{T} c \left(\frac{\tau _{\rm v}}{2} + \xi \right) \, , \end{eqnarray}$$

(25)

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⁻³ and ξ = 1. We seek the values of η and T that simultaneously solve equations (24) and (25).

2.2.4 Accretion criterion

Unlike Sirko & Goodman (2003), the Thompson et al. (2005) model presents an accretion rate |$\dot{M}$| 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_{\rm Edd}$| range, see Heckman et al. (2004); Kollmeier et al. (2006); Suh et al. (2015); Kong & Ho (2018). This introduces a minimum threshold for |$\dot{M}_{\rm out}$|⁠. Thompson et al. (2005) argue that accretion rates of |${\sim } 1 - 10 \rm \, M_\odot {\rm yr}^{-1}$| at the inner disc boundary r_min are sufficient to produce a bright AGN when the central BH mass is |${\sim } 10^9 \rm \, M_\odot$|⁠. This is equivalent to a minimum BH accretion rate of |$\dot{M} \sim 0.2~\dot{M}_{\rm Edd} = 0.2~L_{\rm Edd}/(0.1 c^2)$| at r = r_min. Using equation (47) in Thompson et al. (2005), we find that over a wavelength range of |$[10^{-8} \, \mathrm{m}, 10^{-3} \, \mathrm{m}]$|⁠, setting |$0.2~\dot{M}_{\rm Edd}$| gives a disc bolometric luminosity of |$2\times 10^{-4} \, L_{\rm Edd}$|⁠.

There is no general expression that relates the accretion rate |$\dot{M}_{\rm out}$| and outer radius r_out to the BH mass M that would ensure a bright AGN disc. None the less, we can attempt to find such a relationship by considering how the accretion rate at the outer boundary |$\dot{M}_{\rm out}$| scales with the size of the disc r_out and the central BH mass. Thompson et al. (2005) proposes a critical value |$\dot{M}_\mathrm{c}$|⁠, obtained by equating the star formation time-scale τ_* = 1/ηΩ with the advection time-scale τ_adv = r/v_r to determine whether enough material reaches the central BH to form a luminous signal. From equation (19) we find

$$\begin{eqnarray} \dot{M}_\mathrm{c} = 4 \pi r^2 \rho h \eta \Omega _\mathrm{T} \, . \end{eqnarray}$$

(26)

Together with equation (26), this result can be used to introduce a dependence on r_out and M to |$\dot{M}_{\rm out}$|⁠. A BH with |$M=10^8 \, \rm M_\odot$| surrounded by a TQM05 disc that has r_out = 95 pc, |$\dot{M}_{\rm out}= 320 \, \rm M_\odot {\rm yr}^{-1}$|⁠, and |$\sigma = 188~{\rm km\, s}^{-1}$| satisfies |$\dot{M}_{\rm out}\gt \dot{M}_\mathrm{c} (r = r_{\rm out})$| and has an accretion rate near the central BH of |${\approx } 1.93 \, {\rm M}_\odot {\rm yr}^{-1} = 0.74 {\dot{M}}_{\rm Edd}$| (giving a disc luminosity of |$9.61 \times 10^{-4} \, L_{\rm Edd}$|⁠). We use these values to keep the ratio of |$\dot{M}_{\rm out}/\dot{M}_{\rm c}$| constant. Using equations (16), (14), (26), and assuming the optically thick regime (see Appendix  A), one can show that |$\dot{M}_\mathrm{c} \propto r \sigma ^2$|⁠. Therefore, we scale |$\dot{M}_{\rm out}$| with the outer boundary of the disc and the dispersion velocity, i.e.

$$\begin{eqnarray} \dot{M}_{\rm out}= 320 \, {\rm M}_\odot {\rm yr}^{-1} \, \left(\frac{r_{\rm out}}{95 {\rm pc}}\right) \left(\frac{\sigma }{188 \,{\rm km\,s}^{-1}} \right)^{2} \, . \end{eqnarray}$$

(27)

However, for high masses, equation (27) is not enough to fulfill the the |$\dot{M}_{\rm out}\gt \dot{M}_{\rm c}$| criterion. The inability of |$\dot{M}_{\rm c}$| to accurately predict whether a bright AGN disc is formed is not surprising, as it compares the time-scales for only one value of r. In Thompson et al. (2005), it is stated that discs with |$\dot{M}_{\rm out}\lt \dot{M}_\mathrm{c} (r=r_{\rm out})$| cannot form bright AGNs.

We find that using |$\dot{M}_\mathrm{c}$| from equation (26) as a threshold is too stringent and often omits signals that produce bright AGN discs. In the following, we use equation (27) as an initial guess for |$\dot{M}_{\rm 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 |$\dot{M}_{\rm out}$| that is too high, TQM05 suggests a maximum limit for |$\dot{M}_{\rm out}$| equal to |$\dot{M}_{\rm max} = 8 \pi \rho h \sigma ^2 r/\epsilon _{\rm T} c = L_{\rm M}/\epsilon _{\rm T} c^2$|⁠, where L_M 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, Quataert & Thompson 2005).

Other authors have used different values for |$\dot{M}_{\rm out}$| and r_out. For instance, Stone, Metzger & Haiman (2017) scale the AGN disc down to a Milky-Way type galaxy, using |$M = 3 \times 10^6 \, \rm {\rm M}_{\odot }$|⁠, |$\dot{M}_{\rm out}= 15 \, \dot{M}_{\rm Edd}$|⁠, and |$r_{\rm out}= 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, Lacy & Pope 2022).

2.2.5 Disc profiles

Fig. 4 reproduces Fig.6 in Thompson et al. (2005), assuming either the Semenov opacities (as was done by Thompson et al. 2005) or the combined opacities. The input parameters are σ = 300 km s⁻¹, |$M \approx 10^9 \, {\rm M}_{\odot }$| [instead of equation (15) we use the M-σ relation by Tremaine et al. (2002) as was done by Thompson et al. (2005)], |$\dot{M}_{\rm out}= 320 {\rm M}_{\odot }$| yr⁻¹, r_out = 200 pc, m = 0.2, ϵ_T = 10⁻³, and ξ = 1. Our results shown in Fig. 4 are generally in agreement with those by Thompson et al. (2005).

Figure 4.

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 p_gas, the radiation pressure p_rad, the gas pressure fraction β, the half-thickness of the disc h, the sound speed c_s, 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⁻³, m = 0.2, |$\dot{M}_{\rm out}= 320 {\rm M}_{\odot }$| yr⁻¹, and r_out = 200 pc. Models shown in blue use the opacities by Semenov et al. (2003), models shown in orange use the combined data sets from Semenov et al. (2003) and Badnell et al. (2005).

Figure 5.

Aspect ratio h/r, mass density ρ, optical depth τ_v, and midplane temperature T as functions of cylindrical radius r for both the SG03 (left) and TQM05 (right) AGN disc models. We vary the central BH mass |$M=10^6 \, {\rm M}_\odot$| (blue), |$10^8 \, {\rm M}_\odot$| (orange), and |$10^{10} \, {\rm M}_\odot$| (green). For the SG03 case, we set α = 0.01, l_E = 0.5, and b = 0. For the TQM05 case, we set m = 0.2, ϵ_T = 10⁻³, and ξ = 1. The outer radius r_out and outer accretion rate |$\dot{M}_{\rm out}$| are both scaled with the central BH mass such that r_out = 95 pc and |$\dot{M}_{\rm out}= 320 \, {\rm M}_\odot {\rm yr}^{-1}$| when |$M = 10^8 \, {\rm M}_\odot$|⁠, except for the |$M = 10^{10} \, {\rm M}_\odot$| disc which has an outer accretion rate set to |$\dot{M}_{\rm out}= 1.5 \times 10^6 \, {\rm M}_\odot {\rm yr}^{-1}$|⁠.

Figure 6.

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 \, {\rm M}_\odot$| 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.

Our implementation results in disc profiles that diverge when the disc is close to the central BH, around r ∼ 10⁻³ pc in this case; this follows from the r → r_min limit in the definition of |$\dot{M}^{\prime }$|⁠. We report good agreement between the two opacity implementations, with a noteworthy difference being the presence of the iron opacity bump (Jiang, Davis & Stone 2016) at a radius of r ∼ 2 × 10⁻⁴ pc, seen only for the combined opacities. For both sets of opacities, the disc profile presents a sharp feature at r ∼ 5 × 10⁻¹ pc where the temperature becomes high enough to leave the so-called opacity gap (the dip in κ for temperatures |$10^3\, {\rm K} \lesssim T \lesssim 10^4\, {\rm K}$|⁠, see Semenov et al. 2003; Thompson et al. 2005). Fig. 4 shows that for this set of parameters the disc profile is not sensitive to the choice of opacity tables.

3 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.

3.1 Mass dependence

We first investigate the behaviour of both models as a function of the mass of the central BH. Fig. 5 compares the SG03 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⁶, 10⁸, and |$10^{10} \, {\rm M}_{\odot }$|⁠. These five output quantities can be used to fully reconstruct an AGN disc for both models. Results are presented using the combined opacity data sets.

For the Sirko & Goodman (2003) model, we set α = 0.01, l_E = 0.5, and only consider the α disc (i.e. b = 0). For each disc, we find the solution up to a radius of 10⁷R_s, with the |$M=10^6 \rm \, M_\odot$| case having a maximum extension of ∼1 pc, and the |$M = 10^{10} \rm \, M_\odot$| case extending to ∼1 kpc.

The temperature of the SG03 disc is higher at small separations for lower masses. In particular, one has r ∝ R_s ∝ M in Fig. 5, so that Ω_S ∝ M^−1/2, and thus T ∝ M⁻³ in the inner region, cf. Equation (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³ 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_{\rm eff}^4 \propto \tau _{\rm v}T^4$|⁠), then from equation (10) we find that |$c^2_\mathrm{s} \propto \tau _{\rm v}^2 T^4/\rho$|⁠. If κ is independent of r, then τ_v ∝ ρh, which in hydrostatic equilibrium gives a constant T independent of both r and M.

Fig. 5 shows that the density ρ is lowest when the central BH mass is highest, with ρ ∝ M² in the inner region of the SG03 disc. The model with |$M=10^{10} \rm \, {\rm M}_{\odot }$| presents the thickest SG03 disc, reaching h/r > 1 at r ≳ 10⁶R_s; 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 m_T = 0.2, ϵ_T = 10⁻³, and ξ = 1. In Fig. 5, we linearly scale the outer boundary of the disc r_out using the Schwarzschild radius so that r_out = 10⁷R_s for all three BH masses. We calculate |$\dot{M}_{\rm out}$| using equation (27) for all three BH masses, but find that for the |$M=10^{10} \, {\rm M}_\odot$| case the scaled |$\dot{M}_{\rm out}$| does not satisfy the |$\dot{M}_{\rm out}\gt \dot{M}_{\rm 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⁴ K). Instead, we opt for |$\dot{M}_{\rm out}= 1.5 \times 10^6 \, {\rm M}_\odot {\rm yr}^{-1}$| when |$M=10^{10} \, {\rm M}_\odot$| instead. Equation (27) gives |$\dot{M}_{\rm out}= 0.37 \, {\rm M}_\odot {\rm yr}^{-1}$| when |$M = 10^6 \, {\rm M}_\odot$| and |$\dot{M}_{\rm out}= 322 \, {\rm M}_\odot {\rm yr}^{-1}$| when |$M = 10^8 \, {\rm M}_\odot$|⁠. The AGN disc with |$M=10^8 \, {\rm M}_{\odot }$| has an outer boundary of |$100\,$| pc, which is about half the size of the model shown in Fig. 4.

The |$M=10^6 \, {\rm M}_{\odot }$| case in Fig. 5 shows an AGN disc with an outer boundary |$r_{\rm out}\approx 1\,$|pc and a BH accretion rate |$0.37 \, {\rm M}_\odot {\rm yr}^{-1}$|⁠. Its accretion rate |$\dot{M}$| is higher than both the star formation rate and |$\dot{M}_{\rm c}$| for all values of r, leading to temperatures as large as |$T \sim 10^6 \rm \, K$| at r = r_min and a disc luminosity of |$2 \times 10^{-5} \, L_{\rm Edd}$|⁠. The radiation pressure in such a high-temperature region leads to a thick disc, with h/r > 1 below r ∼ 5 × 10²R_s. 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 |$\dot{M}_{\rm out}$| or decrease m_T.

We find that the model with |$M=10^{10}\, \rm {\rm M}_{\odot }$| also reaches h/r > 1 but at |$r \gt 10^5\, R_{\rm s}$|⁠. 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} \, \rm {\rm M}_{\odot }$| AGN disc extends out to 10 kpc and has an accretion rate of |${\sim } 10 \, \rm {\rm M}_{\odot } {\rm yr}^{-1} = 0.04 ~\dot{M}_{\rm Edd}$| at r = r_min, giving a disc luminosity of |$0.07 \, L_{\rm Edd}$|⁠. For the TQM05 model with |$M=10^{10}\, \rm {\rm M}_{\odot }$|⁠, the optical depth τ_V shows oscillations at r ∼ 200R_s (see Fig. 5) which are due to the model switching between the inner and outer regimes back and forth when close to the Q_T = 1 boundary.

3.2 Input parameters

The SG03 model has five input parameters: the mass of the central BH M, the luminosity ratio l_E (or alternatively the accretion rate |$\dot{M}$|⁠), the disc viscosity α, the BH radiative efficiency ϵ_S, and the pressure flag b = 0, 1. We consider a fiducial model with |$M=10^8 \rm \, M_\mathrm{\odot }$| ϵ_S = 0.1, α = 0.01, l_E = 0.5 and b = 0. Of these parameters, Fig. 6 explores the effect of varying α and l_E.

The density ρ in the outer regime is largely independent of α and l_E. The Shakura & Sunyaev (1973) parameter α relates the viscosity to pressure and accretion, cf. equation (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 l_E = 10⁻³ to l_E = 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 parametrizes 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 Section 2.1.4) indicate that the SG03 model maintains a constant temperature and density at r ≳ 10⁵R_s. Higher accretion rates leads to higher effective temperatures [equation (8)], higher disc temperatures overall [equation (12)], and higher total pressure in the disc [equation (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 equation (15), the star formation efficiency ϵ_T, the efficiency of angular momentum transport m_T in the disc, the supernovae radiative fraction ξ, the outer boundary of the disc r_out, and the accretion rate at this outer boundary |$\dot{M}_{\rm out}$|⁠. Fig. 7 assumes a fiducial model with |$M=10^8 \, {\rm M}_\odot$|⁠, r_out = 10⁷R_s, ϵ_T = 10⁻³, ξ = 1, m_T = 0.2, and |$\dot{M}_{\rm out}\approx 312 \, {\rm M}_\odot {\rm yr}^{-1}$| from equation (27). Starting from this set of parameters, we explore how the disc profile changes when varying either |$\dot{M}_{\rm out}$| or m_T.

Figure 7.

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 \, {\rm M}_\odot$|⁠, ϵ_T = 10⁻³, ξ = 1, and r_out = 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 |$\dot{M}_{\rm out}= 15 \, {\rm M}_\odot {\rm yr}^{-1}$| (blue), |$\dot{M}_{\rm out}= 100 \, {\rm M}_\odot {\rm yr}^{-1}$| (orange), and |$\dot{M}_{\rm out}= 300 \, {\rm M}_\odot {\rm yr}^{-1}$| (green, dashed). The |$\dot{M}_{\rm out}= 300 \, {\rm M}_\odot {\rm yr}^{-1}$| case is dashed to show that parameter profiles are identical to the those of the |$\dot{M}_{\rm out}= 100 \, {\rm M}_\odot {\rm yr}^{-1}$| case close to the central BH. In the right column, we consider AGN discs with an outer accretion rate |$\dot{M}_{\rm out}\simeq 312 \, {\rm M}_\odot {\rm 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.

We consider three values of the accretion rate: |$\dot{M}_{\rm out}=15$|⁠, 100, |$300 \, {\rm M}_\odot {\rm yr}^{-1}$|⁠. The lowest accretion rate considered, |$\dot{M}_{\rm out}= 15 \, {\rm M}_\odot {\rm yr}^{-1}$|⁠, falls below the critical accretion rate |$\dot{M}_\mathrm{c} \approx 21 \, {\rm M}_\odot {\rm yr}^{-1}$| from equation (26) at r = r_out. According to this criterion, this model should not produce an AGN that is sufficiently bright. At r = r_min, the accretion rate for the |$\dot{M}_{\rm out}= 15 \, {\rm M}_\odot {\rm yr}^{-1}$| case is |$\sim 0.58 \, {\rm M}_\odot {\rm yr}^{-1} = 0.22 \, \dot{M}_{\rm Edd}$|⁠, which is below the |$1-10 \, {\rm M}_\odot {\rm yr}^{-1}$| threshold indicated by Thompson et al. (2005). For this case, the disc luminosity is |$1.8 \times 10^{-4} \, L_{\rm Edd}$|⁠, which still falls in the range of Eddington ratios one might expect for AGN discs. This further shows that |$\dot{M}_{\rm 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 |$\dot{M}_\mathrm{c}$| criterion, we find that the profiles become identical when in the inner, non-star forming regime, see the region left of the Q_T = 1 line in Fig. 7. For these cases, the advection time-scale and star formation time-scale reach an equilibrium at the opacity gap (τ_adv = τ_* when |$T \approx 10^3 \, {\rm K}$|⁠). This leads to discs of the same temperature, density, aspect ratio, and accretion rate (⁠|$\dot{M}= 2.23 \, {\rm M}_\odot {\rm yr}^{-1}$| at r = r_min for both discs).

The global torque efficiency parameter m_T is strongly correlated to the behaviour of the disc for all radial distances. Much like α for the SG03 model, here m_T parametrizes the relationship between the angular momentum transport and the accretion rate, cf. equation (19). In the outermost regions of the disc, where the accretion rate |$\dot{M}\approx \dot{M}_{\rm out}$| is roughly constant and h/r ∼ 1, the total pressure is inversely proportional to m_T, see equation (19) and the definition of c_s. This gives a thinner TQM05 disc, with lower values of h/r for higher values of m_T at the outer boundary. The density is constant in the marginally stable outer region because of equation (17), but the low total pressure causes some temperature deviations at r ≈ 10⁷R_s for each disc we consider. These variations contribute to different initial conditions in τ_v for each value of m_T. The TQM05 disc has similar behaviour for all three m_T values once the solutions enter the opacity gap at r ≈ 10⁵R_s, 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 m_T values lead to thick, low density discs due to low radiation pressure (which is proportional to τ_v by definition).

4 DISC MIGRATION

Migration in gas discs was first proposed based on formulae that considered how spiral-wave structures can be sustained in galaxies (Feldman & Lin 1973; Goldreich & Tremaine 1979; Lin & Papaloizou 1979). An object orbiting in a gas disc exchanges angular momentum with its surroundings, leading to changes in its orbit and thus a net radial migration in the disc (Armitage 2020). These migration torques were predicted by Goldreich & Tremaine (1979), described for planets in protoplanetary discs by Ward (1997), improved upon by Paardekooper & Mellema (2006); Paardekooper et al. (2010), studied for the case of planets by Nelson et al. (2000); Tanaka, Takeuchi & Ward (2002); Paardekooper & Mellema (2006); Kley & Crida (2008); Lyra, Paardekooper & Mac Low (2010); Paardekooper et al. (2010), and extended to the AGN case by Syer, Clarke & Rees (1991); Artymowicz, Lin & Wampler (1993); Levin (2003); McKernan et al. (2011, 2012); Bellovary et al. (2016); Derdzinski & Mayer (2023). A key phenomenon emerging from these studies is the potential occurrence of migration traps – locations in the gas disc where the net radial migration torque is zero. Depending on the mass ratio between the migrator, the central object or the disc, the migrator may clear a gap (referred to as Type II migration) or deplete material at the migration trap without clearing a gap (Type I migration) (Ward 1997). In this work, we only consider the case of Type I 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. (2023) showed that 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.

Table 2.

Summary of the parameter entering our treatment of disc migration explored in Section 4.

Symbol	Definition
mBH	Mass of the migrating object
q	Mass ratio between migrator and central BH
Γ0	Normalization migration torque
ΓI	Type I migration torque
γ	Adiabatic index
CL	Lindblad torque
χ	Thermal diffusivity of the disc
Γtherm	Thermal torque
xc	Corotation radius of the migrator
λ	Size of the thermal lobes
L	Migrator luminosity from thermal heating
Lc	Critical migrator luminosity

4.1 Torque implementation

In particular, we apply our AGN disc models to the methods by Grishin et al. (2023), adopting their migration torque and thermal torque expressions. Grishin et al. (2023) use a simpler AGN disc model where profiles are power laws in M, r, and accretion rate |$\dot{M}$|⁠. Their discs are relatively similar to the SG03 models with M = 10⁶M_⊙ and α = 0.01. When using migration torques by Paardekooper et al. (2010) which assume the disc is locally isothermal, Grishin et al. (2023) 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. (2023) 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 overdense spiral arms in the disc. Each arm will produce a torque acting on the migrating object with a magnitude (Korycansky & Pollack 1993)

$$\begin{eqnarray} \Gamma _0 = q^2 \Sigma _\mathrm{g} r^4 \Omega ^2 \left(\frac{h}{r}\right)^{-3} \, , \end{eqnarray}$$

(28)

where q ≡ m_BH/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):

$$\begin{eqnarray} \Gamma _{\rm I} = \left(-0.85 + 0.9 \frac{\mathrm{d}\ln \Sigma _\mathrm{g}}{\mathrm{d}\ln r} + \frac{\mathrm{d}\ln T}{\mathrm{d}\ln r} \right) \frac{h}{r} \, \Gamma _0 \, . \end{eqnarray}$$

(29)

Jiménez & Masset (2017) updates the migration torque formula to

$$\begin{eqnarray} \Gamma _{\rm I} = \left[C_\mathrm{L} + \left(0.46 + 0.96 \frac{\mathrm{d}\ln \Sigma _\mathrm{g}}{\mathrm{d}\ln r}-1.8 \frac{\mathrm{d}\ln T}{\mathrm{d}\ln r} \right)\gamma ^{-1} \right] \frac{h}{r}\, \Gamma _0 \, , \end{eqnarray}$$

(30)

where γ = 5/3 is the adiabatic index. The parameter

$$\begin{eqnarray} C_\mathrm{L} = \left(-2.34 -0.1 \frac{\mathrm{d}\ln \Sigma _\mathrm{g}}{\mathrm{d}\ln r} + 1.5 \frac{\mathrm{d}\ln T}{\mathrm{d}\ln r} \right) f_\gamma \left(\frac{\chi }{h^2 \Omega }\right) \, , \end{eqnarray}$$

(31)

is the Lindblad torque where

$$\begin{eqnarray} f_\gamma (x) = \frac{(x/2)^{1/2} + 1/\gamma }{(x/2)^{1/2} + 1} \, \end{eqnarray}$$

(32)

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

$$\begin{eqnarray} \chi = \frac{16 \gamma (\gamma -1)\sigma _\mathrm{SB} T^4}{3 \kappa \rho ^2 h^2 \Omega ^2} \, . \end{eqnarray}$$

(33)

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 underdense lobes form, leading to outward migration (Benítez-Llambay et al. 2015). The total heating torque is (Masset 2017)

$$\begin{eqnarray} \Gamma _\mathrm{therm} = 1.61 \frac{\gamma -1}{\gamma } \frac{x_{\rm c}}{\lambda } \left(\frac{L}{L_\mathrm{c}}-1 \right) \Gamma _0 , \end{eqnarray}$$

(34)

where x_c 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

$$\begin{eqnarray} L_\mathrm{c} = \frac{4 \pi G q M \rho }{\gamma } \chi \end{eqnarray}$$

(35)

is the critical luminosity. If L = L_c, 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 equation (3), replacing the mass M with the mass of the migrator m_BH]. The size of the lobes λ is given by (Grishin et al. 2023)

$$\begin{eqnarray} \lambda = \sqrt{\frac{2 \chi }{3 \gamma \Omega }} \, , \end{eqnarray}$$

(36)

and the corotation radius is (Grishin et al. 2023):

$$\begin{eqnarray} x_{\rm c} = - \frac{h^2}{3 \gamma r} \frac{\mathrm{d}\ln {p_{\rm tot}}}{\mathrm{d}\ln {r}}\, . \end{eqnarray}$$

(37)

We approximate dln p_tot/dln r by combining the equation for vertical hydrodynamical equilibrium p_tot ≈ ρh²Ω² with the definition of the sound speed |$c_{\rm s}^2 = h^2 \Omega ^2$|⁠, resulting in |$\mathrm{d}p_{\rm tot}/\mathrm{d}r \approx \rho c_{\rm s}^2/r$|⁠.

The thermal torque given by equation (34) is expected to diminish in optically thin discs. Following Grishin et al. (2023), we multiply equation (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:

$$\begin{eqnarray} \frac{m_{\rm th}}{m_{\rm BH}} = \frac{\chi }{c_{\rm s} R_{\rm B}} \, , \end{eqnarray}$$

(38)

where R_B is half the Bondi radius

$$\begin{eqnarray} R_{\rm B} = \frac{G m_{\rm BH}}{c_{\rm s}^2} \, . \end{eqnarray}$$

(39)

In the regions where h < R_B we use the disc height h in place of half the Bondi radius R_B. To correct for the critical thermal mass, we split equation (34) into its heating component (the positive L/L_c term) and its cooling component (the negative term), which we label as |$\Gamma _{\rm therm, \, hot }$| and |$\Gamma _{\rm therm, \, cold }$|⁠, respectively. The total thermal torque is described by equation (34) unless μ_th ≡ m_th/m_BH < 1. In regions of the disc where μ_th < 1, the thermal torque is instead given by:

$$\begin{eqnarray} \Gamma _{\rm therm} = \Gamma _{\rm therm, \, hot} \frac{4 \mu _{\rm th}}{1 + 4 \mu _{\rm th}} + \Gamma _{\rm therm, \, cold} \frac{2 \mu _{\rm th}}{1 + 2 \mu _{\rm th}} , \end{eqnarray}$$

(40)

which is an approximation of numerical fits detailed in Velasco Romero & Masset (2020) and used in Guilera et al. (2021); Grishin et al. (2023).

4.2 Migration traps

Fig. 8 shows the migration-torque profiles for a M = 10⁶ M_⊙ SMBH in both the Sirko & Goodman (2003) and Thompson et al. (2005) model. We use ϵ_S = 0.1, α = 0.01, l_E = 0.5, and b = 0 for the SG03 AGN disc and r_out = 10⁷R_s, ϵ_T = 10⁻³, ξ = 1, m_T = 0.2, and |$\dot{M}_{\rm out}= 1.5 \times 10^{-2} \, {\rm M}_{\odot } \, {\rm yr}^{-1}$| for the TQM05 model. The outer accretion rate |$\dot{M}_{\rm out}$| was set smaller than the value given by equation (27) in order to enforce h/r < 1 throughout 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.

Figure 8.

The absolute values of the migration torques for a |$m_{\rm BH} = 10 \, {\rm M}_\odot$| BH orbiting a M = 10⁶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⁻³, ξ = 1, m = 0.2, r_out = 10⁷R_s and |$\dot{M}_{\rm out}= 1.5\times 10^{-2} \, {\rm M}_{\odot } \, {\rm 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). Colours 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.

The top panel shows the migration torque using both equation (29) by Paardekooper et al. (2010) and equation (30) by Jiménez & Masset (2017). When using the former, we find migration traps at r ≈ 22R_s and r ≈ 10³R_s for the SG03 model, which is in line with the results reported by both Bellovary et al. (2016) and Grishin et al. (2023). 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. (2023). Once the thermal torque from equation (34) is added to the updated migration torque of equation (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 \, {\rm M}_\odot$| central BH and a 10 M_⊙ migrator occurring at r ≈ 1.4 × 10³R_s = 1.4 × 10⁻⁴ pc and r ≈ 6.8 × 10⁴R_s = 6.5 × 10⁻³ 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 dln Σ_g/dln r discretely changes values, as can be seen in the lower panels of Fig. 8 at r ≈ 2.5 × 10³R_s and r ≈ 8.3 × 10³R_s. When both migration and thermal torques are considered, we find traps at r ≈ 2.5 × 10³R_s, r ≈ 8.4 × 10³R_s, and r ≈ 2.6 × 10⁶R_s for the Thompson et al. (2005) model.

5 PUBLIC IMPLEMENTATION

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

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

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.

6 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. (2023) 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 modelling approaches. The goal of pAGN is precisely that to aid further research in the growing field of AGN and gravitational-wave science.

ACKNOWLEDGEMENTS

We thank Maria Paola Vaccaro, Andrea Derdzinski, Evgeni Grishin, Ari Laor, Hiromichi Tagawa, Shmuel Gilbaum, Sean McGee, Cressida Cleland for discussions. DG and DGer are supported by ERC Starting Grant No. 945155–GWmining, Cariplo Foundation Grant No. 2021-0555, MUR PRIN Grant No. 2022-Z9X4XS, and the ICSC National Research Centre funded by NextGenerationEU. AAT acknowledges support from JSPS KAKENHI Grant No. 21K13914 and from the European Union’s Horizon 2020 and Horizon Europe research and innovation programs under the Marie Skłodowska-Curie grant agreements No. 847523 and No. 101103134. DGer is supported by MSCA Fellowships No. 101064542–StochRewind and No. 101149270–ProtoBH. Computational work was performed at CINECA with allocations through INFN and Bicocca.

DATA AVAILABILITY

The data underlying this article are available at Gangardt & Trani (2024). Additional data will be shared on reasonable request to the corresponding author.

Footnotes

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 |$\sqrt{2}$| is negligible.

APPENDIX A: OPTICALLY THICK APPROXIMATION IN TQM05

For the TQM05 model, the case where the disc optically thick to its own infrared radiation can be approximated analytically. We assume that the gas is a constant fraction f_g ≡ Σ_g/Σ_tot of the total dynamical mass

$$\begin{eqnarray} \Sigma _\mathrm{tot} = \frac{\sigma ^2}{\pi G r} \, . \end{eqnarray}$$

(A1)

These assumption function best at large scales (i.e. r ≫ r_out) where the angular frequency is dominated by the velocity dispersion, so that |$\Omega _\mathrm{T} \approx \sqrt{2} \sigma /r$|⁠. The mass density from equation (23) reads

$$\begin{eqnarray} \rho \approx \frac{\sqrt{2} \sigma ^2}{\pi G Q_\mathrm{T} r^2} \, . \end{eqnarray}$$

(A2)

If f_g is constant, then the mid-scale height is given by

$$\begin{eqnarray} \frac{h}{r} = \frac{f_\mathrm{g} Q_\mathrm{T}}{2^{3/2}} \, , \end{eqnarray}$$

(A3)

and the sound speed is

$$\begin{eqnarray} \frac{c_\mathrm{s}}{\sigma } = \frac{f_\mathrm{g} Q_\mathrm{T}}{2} \, . \end{eqnarray}$$

(A4)

At large values of r, the disc is mostly radiation-pressure dominated, so that |$p_\mathrm{rad} = 4 \sigma _\mathrm{SB} T^4/3 c = \sigma _\mathrm{SB} \tau _{\rm v}T_{\rm eff}^4/c$|⁠. In the optically thick limit, the main contribution to T_eff that of star formation

$$\begin{eqnarray} \sigma _\mathrm{SB} T_{\rm eff}^4 = \frac{1}{2} \epsilon _\mathrm{T} \dot{\Sigma }_* c^2 \, . \end{eqnarray}$$

(A5)

Combining these equations, one can then find the temperature

$$\begin{eqnarray} T = \left(\frac{3c Q_\mathrm{T}}{2^{7/2} \pi G \sigma _\mathrm{SB}} \right)^{1/4} \left(\frac{f_\mathrm{g} \sigma ^2}{r} \right)^{1/2} \, , \end{eqnarray}$$

(A6)

and the star formation rate

$$\begin{eqnarray} \dot{\Sigma }_* = \frac{\sqrt{2}f_\mathrm{g} Q_\mathrm{T}}{\epsilon _\mathrm{T} \kappa c} \frac{\sigma ^2}{r} \, . \end{eqnarray}$$

(A7)