© ESO, 2017

1. Introduction

The processes which regulate the co-evolution of supermassive black-holes (SMBH) with their host galaxies are still not well understood, as well as the activity of SMBH over cosmological time scales. However, SMBH have been successfully detected both directly and indirectly at the centres of the galaxies and their masses span from millions to billions of solar masses (Fan et al. 2004, 2006; Mortlock et al. 2011; Venemans et al. 2012; Willott et al. 2003). A widely accepted view suggests that the growth of SMBH occurs due to rapid accretion and merging. However, it is still unclear what physical mechanisms drive SMBH to reach such high masses and what gives rise to the correlation with the physical properties of their host galaxies.

Substantial efforts were undertaken and numerous models have been proposed to explain the origin and formation of SMBH through different mechanisms (Rees 1984). However, to explain the formation of SMBH at higher redshift requires study of the first black-hole (BH) seeds. So far, there are three main scenarios which include the formation of BH seeds through the core collapse of Population III stars (Trenti & Stiavelli 2007), the collapse of dense stellar clusters (Devecchi et al. 2010, 2012) and the collapse of protogalactic metal free gas clouds in primodial halos (see e.g. Latif & Ferrara 2016). The latter mechanism provides the potentially largest seed masses and may thus favor the formation of supermassive black-holes at z> 6 (e.g. Lodato & Natarajan 2007; Wise et al. 2008; Begelman & Shlosman 2009; Schleicher et al. 2010a, 2013; Latif et al. 2013a,b; Ferrara et al. 2014). These high-redshift black-holes can potentially be detected with ALMA using CO and H2 line emission (Spaans & Meijerink 2008; Schleicher et al. 2010b).

Several studies have found tight correlations between black-hole mass and galaxy bulge properties in the nearby universe. Such correlations have indicated scaling relations of black-hole mass with bulge mass as well as black-hole mass with the velocity dispersions of their host (Magorrian et al. 1998; Haring & Rix 2004; Marconi & Hunt 2003; Merritt 2000; Gebhardt et al. 2000; Ferrarese & Merritt 2000). Correlations between star formation rate and SMBH growth have also been found by Diamond-Stanic & Rieke (2012). Their work includes a complete sample of Seyfert galaxies and they found a strong correlation between nuclear (r< 1 kpc) star formation rates and the accretion rates of the SMBH with a scaling relation $\mathit{Ṁ}{}_{\mathrm{\ast }}\mathrm{~}\mathit{Ṁ}\begin{array}{c}\mathrm{0.8}\\ \mathrm{BH}\end{array}$. Measurements on scales r> 1 kpc with extended star formation in the host galaxy showed only a weak correlation and scaling as ${\mathit{M}}_{\mathrm{\ast }}\mathrm{~}\mathit{Ṁ}\begin{array}{c}\mathrm{0.6}\\ \mathrm{BH}\end{array}$. Their results suggest a connection between gas on sub-kiloparsec scales that is forming stars and the gas on sub-parsec scales that goes into black-hole accretion. In this case, the transport of angular momentum in the galactic disk would be a natural mechanism which plays a significant role in feeding the central engine and triggering star formation. Additionally, prominent circumnuclear disks with star formation rings have been observed in the inner regions of some galaxies from parsec scale up to a few kilo parsecs from the galactic centre (Lenc & Tingay 2009; Hsieh et al. 2011; van der Laan et al. 2013; Falcon-Barroso et al. 2014).

In recent years, for several active galaxies gas mass inflows to the nucleus have been detected through optical and near-infrared spectroscopy with estimated inflow rates of the order of 10^-2M_⊙ yr^-1 (Schnorr-Müller et al. 2016, 2014a,b; Lena et al. 2015; Storchi-Bergmann et al. 2009, 2010). However, it is difficult to determine the actual kinematics of the inflowing gas to the nucleus from the current observations (Storchi-Bergmann et al. 2010). Another signature of mass inflows onto the SMBH from the CND has been recently reported by Izumi et al. (2016). In their work, through the use of the HCN(1–0) emission line and 2–10 keV X-ray observations of a sample of Seyfert galaxies they have found a positive correlation between dense molecular gas mass in the CND and black-hole mass accretion rate in ten Seyfert galaxies. Such study supports the scenario that rich gas environments such as the CND trigger star formation, which drives mass accretion onto the central SMBH.

Co-evolution models and numerical simulations are necessary for explaining and describing the processes that regulate and control the accretion onto SMBH as well as the gas inflow from the host galaxy which triggers nuclear starbursts from kilo-parsec scale down to parsec scale in active galaxies. Kawakatu & Wada (2008) and Wutschik et al. (2013) have proposed models of co-evolution, where the transport of angular momentum from the outer border of the CND to the SMBH is driven by turbulence injected via supernovae (SN) that lead to different accretion phases in the disk around a SMBH. An important aspect studied by Wutschik et al. (2013) is the functional form of the star formation rate, as different dependencies have been suggested in the literature with respect to the surface density of the gas as well as the turbulent velocity. The models explored included both linear and non-linear dependencies on the gas surface density, as well as models independent of the turbulent velocity v_turb and models where the star formation rate scales as ${\mathit{v}}_{\mathrm{turb}}^{-1}$. Specific models adopted from the literature included in particular the model employed by Kawakatu & Wada (2008) and by Elmegreen & Burkert (2010; hereafter model EB10) .

In this work, we further extend the framework of Wutschik et al. (2013) and apply it to the formation of both nearby AGN and high-redshift quasars. For this purpose, we have included the star formation recipe by Krumholz & McKee (2005; hereafter model KM05), including both non-linear dependencies on gas surface density and turbulent velocity, from which the Kennicutt-Schmidt relation can be derived in a systematic manner. We also explored whether our assumptions regarding the gravitational dynamics have a strong impact on the results. With this new model, we will show that the non-linear dependence in the star formation model, along with feedback from SN, can give rise to star formation histories that are similar to the observed ones.

In Sect. 2, we describe the general implementation along with the new star formation model considered here. Section 3 shows the results of semi-analytical calculations and simulations which describe the time evolution of the black-hole mass, stellar mass and disk radius, exploring different assumptions on the dominant component in the gravitational field. We also show the time evolution of the Mach numbers and the expected evolution of observable quantities as the luminosities. In Sect. 4, the model is then applied to specific objects including local AGN and high-redshift quasars. Finally Sect. 5 presents our discussion and conclusions.

2. Evolution of the SMBH and star formation in a circumnuclear disk

2.1. Disk model

The co-evolution model of Kawakatu & Wada (2008) considers the following system: an accreting SMBH surrounded by a CND, all embedded in the galaxy. The host galaxy supplies matter (dusty gas) with the rate Ṁ_sup to the CND. The disk becomes gravitationally unstable at a critical radius r_c, where the Toomre Q parameter is equal to one. This separates the disk into an outer disk (self-gravitating region) and into an inner disk (stable region). In the gravitationally unstable outer disk stars are forming and the turbulent pressure produced by SN drives gas to the inner disk to feed the SMBH.

In the extension of the model by Wutschik et al. (2013; hereafter model WS13), the stable and the unstable part of the CND are treated separately, assuming a power-law surface density distribution where the power-law index is obtained from normalization constraints. The mass exchange between both parts of the disk is driven by direct mass exchange due to accretion from the unstable outer disk to the inner disk, as well as through geometric effects, that is, if the critical radius moves towards the inner or outer part of the disk. The model also includes the cases where the entire CND is stable or unstable where the critical radius is equal to the inner or outer radius of the CND (see also Kawakatu & Wada 2008). Wutschik et al. (2013) and Kawakatu & Wada (2008) show that depending on the stability of the CND, one can distinguish between the following phases:

• a)

  High accretion phase (constant supply from the hostgalaxy): Due to higher matter density or lower turbulent velocity, stable matter in the disk becomes unstable , which results into the critical radius moving inward: r_c(t_o) >r_c(t) or r_c ≤ r_in. In consequence, gas ΔM_g from the inner gas reservoir flows to the outer reservoir. Hence, the gas mass in the outer region M_g,o is defined as M_g,o(t) = M_g,o(t_o) −ΔM_g. Thus, the disk becomes geometrically thick, stars form and their feedback through SN leads to the turbulent supersonic phase. In this phase, the disk thickness is regulated by the turbulent pressure of the gas P_g which is defined as ${\mathit{P}}_{\mathrm{g}}\mathrm{=}{\mathit{\rho }}_{\mathrm{g}}{\mathit{v}}_{\mathrm{turb}}^{\mathrm{2}}$, where ρ_g is the gas density and the turbulent velocity v_turb is generally larger than the sound speed c_s, that is, v_turb>c_s. The turbulent energy in the disk is ${\mathit{E}}_{\mathrm{turb}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{g}}{\mathit{v}}_{\mathrm{turb}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$, where Σ_g is the surface density of the gas. The total energy injected by SN is given as $\mathit{Ė}{}_{\mathrm{inj}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathrm{SN}}{\mathit{\eta }}_{\mathrm{SN}}{\mathit{E}}_{\mathrm{SN}}\hat{\mathrm{\Sigma }}{}_{\mathrm{\ast }}\mathrm{(}\mathit{r,t}\mathrm{-}{\mathit{t}}_{\mathrm{SN}}\mathrm{)}$, with t_SN is the average life-time of massive stars which then explode into SN, E_SN the energy injected into the ISM via SN, η_SN the heating efficiency. and f_SN is a constant parameter. $\hat{\mathrm{\Sigma }}{}_{\mathrm{\ast }}\mathrm{(}\mathit{r,t}\mathrm{-}{\mathit{t}}_{\mathrm{SN}}\mathrm{)}$ is the local star formation rate per unit area. We assume here that the energy not going into turbulence is dissipated (E_dis) through other processes, and due to the energy balance, we then have Ė_turb(r) + Ė_dis(r) = Ė_inj. As a typical lifetime for massive stars, we adopt 10⁶ yr, but we will also explore variations of this parameter in Sect. 4. The feedback from SN is the primary mechanism which drives the gas to the inner disk to feed the black-hole. In summary, during this phase, accretion becomes highly efficient.

  Table 1

  Model parameters.

• b)

  Low accretion phase (after termination of supply from host galaxy): Gravitationally unstable matter in the disk (in a state of active star formation) becomes stable which results into the critical radius moving outward: r_c(t_o) <r_c(t) or r_c>r_in. So the gas reservoir of the inner stable disk therefore increases, as a larger part of the CND is becoming stable over time. Once stable, the disk is vertically supported only by thermal gas pressure, ${\mathit{P}}_{\mathrm{g}}\mathrm{=}{\mathit{\rho }}_{\mathrm{g}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$, where v_turb<c_s. The scale height of the disk is determined by the balance between gravity and thermal gas pressure. In general, accretion is inefficient in this phase.

In our model we employ the input parameters for the outer and inner radius of the disk as ${\mathit{r}}_{\mathrm{out}}\mathrm{=}\sqrt{{\mathit{M}}_{\mathrm{BH}}\mathit{/}\mathit{\pi }{\mathrm{\Sigma }}_{\mathrm{host}}}$ and ${\mathit{r}}_{\mathrm{in}}\mathrm{=}\mathrm{3}\mathrm{pc}\sqrt{{\mathit{M}}_{\mathrm{BH}}\mathit{/}{\mathrm{10}}^{\mathrm{8}}{\mathit{M}}_{\mathrm{\odot }}}$ respectively (Wutschik et al. 2013; Kawakatu & Wada 2008). Table 1 summarises the fixed and free parameters of our model. The values for the time of supply t_sup, the surface density of the host galaxy Σ_host, the power law exponent for the surface gas density γ, the parameters related to the energy injected by SN (t_SN, f_SN, η_SN and E_SN) and the thermal sound speed c_s are taken from previous works (Wutschik et al. 2013; Kawakatu & Wada 2008). Considering mass conservation, the mass balance in the system is M_sup = M_BH + M_g + M_∗, where M_sup is the mass supplied to the CND, M_BH is the black-hole mass, M_g is the gas mass and M_∗ is the stellar mass. The time evolution of the gas mass in the disk is given as ${\mathit{M}}_{\mathrm{g}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathrm{\left[}\mathit{Ṁ}\sup \mathrm{\right(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{)}\mathrm{-}\mathit{Ṁ}\mathrm{\ast }\mathrm{(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{)}\mathrm{-}\mathit{Ṁ}\mathrm{BH}\mathrm{(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\left)}\mathrm{\right]}\hspace{0.17em}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\mathit{,}$(1)considering the difference between the mass that was supplied and the mass going into stars and the central black-hole. Ṁ_sup is the mass supply rate, Ṁ_BH is the black-hole mass accretion rate and Ṁ_∗ is the star formation rate in the disk. The latter is given by the expression: $\mathit{Ṁ}\mathrm{\ast }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{{\mathit{r}}_{\min }}^{{\mathit{r}}_{\mathrm{out}}}\mathrm{2}\mathit{\pi r}\mathrm{\Sigma ̇}\mathrm{\ast }\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r,}$(2)where the value of r_min depends also on the disk stability and is given by r_min = max(r_in,r_c). To estimate the black-hole mass accretion rate in a viscous accretion disk, we refer to the expression given by Pringle (1981) as $\mathit{Ṁ}\mathrm{BH}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\pi \nu }{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\frac{\mathrm{d}\ln \mathrm{\Omega }\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{d}\ln \mathit{r}}\mathit{,}$(3)where ν = αv(r)h(r) denotes the viscosity, with α the viscous coefficient, v(r) the gas velocity and h(r) the scale height of the disk. Here Ω represents the angular velocity of the gas, which is obtained from the balance between the centrifugal force and the gravitational force. The latter depends on the dominant gravitational component, which can be the black-hole, the CND itself and the host galaxy. In the general case, it is thus given as ${\mathrm{\Omega }}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathit{G}{\mathit{M}}_{\mathrm{BH}}}{{\mathit{r}}^{\mathrm{3}}}\mathrm{+}\frac{\mathit{\pi G}}{\mathit{r}}\mathrm{\left(}{\mathrm{\Sigma }}_{\mathrm{disk}}\mathrm{\right(}\mathit{r}\mathrm{)}\mathrm{+}{\mathrm{\Sigma }}_{\mathrm{host}}\mathrm{)}\mathit{,}$(4)where Σ_host is the host galaxy surface density, Σ_disk = Σ_g + Σ_stars is the disk surface density. It follows a power law function of the radius, represented as Σ_disk ≈ Σ_g = Σ₀(r/r_c)^{− γ}, where γ is a free parameter (Wutschik et al. 2013). The evolution of the black-hole mass M_BH from an initial seed mass M_BH,seed is given as ${\mathit{M}}_{\mathrm{BH}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{M}}_{\mathrm{BH}\mathit{,}\mathrm{seed}}\mathrm{+}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{t}}\mathit{Ṁ}\mathrm{BH}\mathrm{\left(}{\mathit{t}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{t}}^{\mathrm{\prime }}\mathit{.}$(5)

2.2. Star formation rate model

As a new ingredient in our model, we consider here the star formation model of KM05 due to its effective applicability not only to the Milky Way’s star formation rate but also to a vast variety of galaxies and its consistency with observations and simulations. The model particularly allows to derive the two forms of the Kennicutt-Schmidt relation and successfully predicts star formation rates based on observable quantities over a large range of conditions, spanning from normal galactic disks, circumnuclear starburts and ULIRGs to nearby star-forming galaxies such as Markarian 273, Arp 193 and Arp 220.

In the analytical description of KM05, the main premise is that the gas collapses in sub-regions of molecular clouds in a supersonically turbulent environment to form stars. The total star formation rate is formulated in terms of the observable properties of the gas in the galactic disk as follows $\mathrm{\Sigma ̇}\mathrm{\ast }\mathrm{\approx }\mathrm{0.073}{\mathit{M}}^{-0.32}{\mathit{\phi }}_{\mathit{\rho }}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{Q}}_{\mathrm{1.5}}^{-1}{\mathit{f}}_{\mathrm{GMC}}\mathrm{\Omega }{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{[}{\mathit{M}}_{\mathrm{\odot }}\hspace{0.17em}{\mathrm{p}{\mathrm{c}}^{-2}\hspace{0.17em}\mathrm{y}{\mathrm{r}}^{-1}}^{\mathrm{]}}\mathit{,}$(6)where M is the mach number which is the ratio of the turbulent velocity v_turb to the local sound speed c_s. The observable parameter φ_ρ is given by the ratio of the mean molecular density to the mean midplane density. This ratio has been approximated as , where (Krumholz & McKee 2005). The latter can be expressed as ${\mathit{\phi }}_{\mathit{\rho }}\mathrm{\approx }\mathrm{5}\mathrm{[}\frac{\mathrm{10}\mathrm{-}\mathrm{8}{\mathit{f}}_{\mathrm{GMC}}}{\mathrm{6}}{\mathrm{]}}^{\mathrm{3}\mathit{/}\mathrm{4}}$, where f_GMC is the molecular gas fraction. f_GMC = 0.5 (it implies φ_ρ = 0.6) is considered as an universal parameter which applies from starbursts to normal galaxies, as stated in the Appendix B of Krumholz & McKee (2005). The parameter Q represents the Toomre parameter (Toomre 1964; Rice 2016) in Eq. (6), Q_1.5 = Q/ 1.5.

To evaluate the effect of the star formation model on the BH mass accretion rate and on the star formation rate evolution in the CND, we have re-written Eq. (6) by emphasizing the dependence of the star formation rate on the turbulent velocity and angular velocity: $\mathrm{\Sigma ̇}\mathrm{\ast }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{\Psi \Omega }\mathrm{\left(}\mathit{r}\mathrm{\right)}{\left(\frac{{\mathit{v}}_{\mathrm{turb}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{c}}_{\mathrm{s}}}\right)}^{-0.32}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}$(7)We have written the fix parameter Ψ in terms of Q and f_GMC as $\mathrm{\Psi }\mathrm{=}\mathrm{0.073}\sqrt{\mathrm{5}}{{}^{\mathrm{\right(}}\frac{\mathrm{10}\mathrm{-}\mathrm{8}{\mathit{f}}_{\mathrm{GMC}}}{\mathrm{6}}^{\mathrm{\left)}}}^{\mathrm{3}\mathit{/}\mathrm{8}}\frac{\mathrm{1.5}}{\mathit{Q}}{\mathit{f}}_{\mathrm{GMC}}$. Now that we have condensed Eq. (6), we employ in the next section the formulation given in Eq. (7) for studying the evolution of the circumnuclear disk. If (part of) the CND is gravitationally unstable, the turbulent velocity in this region will be calculated based on the energy input from SN, leading to supersonic turbulence. As in the WS13 model, we approximate the turbulent velocity here as a power law via ${\mathit{v}}_{\mathrm{turb}}\mathrm{=}{\mathit{v}}_{\mathrm{o}}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{out}}}\right)}^{\mathit{\lambda }}\mathit{,}$(8)where ${\mathit{v}}_{\mathrm{o}}\mathrm{=}\sqrt{{\mathit{E}}_{\mathrm{turb}}\mathrm{\left(}{\mathit{r}}_{\mathrm{c}}\mathrm{\right)}\mathit{/}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}{\mathit{r}}_{\mathrm{c}}\mathrm{\right)}}$ and λ is the power law index expressed as $\mathit{\lambda }\mathrm{=}\frac{\mathrm{1.5}\mathrm{-}\mathit{\gamma }\mathrm{(}\mathit{\theta }\mathrm{-}\mathrm{1}\mathrm{)}}{\mathrm{(}\mathrm{2}\mathrm{+}\mathit{ϵ}\mathrm{)}}$, which is 0 in the sub-sonic turbulence phase. The parameters θ and ϵ come from the parametrization of the total star formation rate as $\mathrm{\Sigma ̇}{}_{\mathrm{\ast }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{\Psi }\mathrm{·}{{}^{\mathrm{(}}\mathrm{\Sigma }_{\mathrm{g}}\mathrm{(}\mathit{r}{\mathrm{\left)}}^{\mathrm{\right)}}}^{\mathit{\theta }}\mathrm{·}\mathrm{(}{\mathit{v}}_{\mathrm{turb}}{\mathrm{)}}^{\mathrm{-}\mathit{ϵ}}$ (Wutschik et al. 2013) . In the KM05 model the latter two parameters have the values of 1 and 0.32 respectively (see Eq. (7)), hence, we used them for calculating the value of λ, which is approximately equal to 0.65.

2.3. Dependence of accretion on dominant gravitational source

As described in Sect. 2.2, the total star formation rate depends non-linearly on the turbulent velocity and linearly on the surface density of gas. However, since in Eq. (4) the angular velocity depends on different gravitational source terms, we have inserted this relation into Eq. (7) to write the total star formation rate as $\begin{array}{ccc}\mathrm{\Sigma ̇}\mathrm{\ast }\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \mathrm{\Psi }\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathrm{BH}}}{{\mathit{r}}^{\mathrm{3}}}\mathrm{+}\frac{\mathit{\pi G}}{\mathit{r}}\left({\mathrm{\Sigma }}_{\mathrm{disk}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}{\mathrm{\Sigma }}_{\mathrm{host}}\right)}\\ & & \mathrm{\times }{\left(\frac{{\mathit{v}}_{\mathrm{turb}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{c}}_{\mathrm{s}}}\right)}^{-0.32}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}\end{array}$(9)To evaluate the evolution of the star formation rate Ṁ_∗(t) in the CND, we insert Eq. (9) into (2). To calculate Eq. (2) analytically, we have first considered three different limiting cases of , depending on the dominant source of gravity in the system: the central black-hole, the CND itself or the gravitational field from the host galaxy. As a result, we obtained expressions for the star formation rate Ṁ_∗(t) and the energy injection rate Ė_inj in different regimes.

If gravity is dominated by the central black-hole, we assumed that $\mathrm{\Sigma ̇}{}_{\mathrm{\ast }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{\Psi }\sqrt{{}^{\mathrm{(}}\frac{\mathit{G}{\mathit{M}}_{\mathrm{BH}}}{{\mathit{r}}^{\mathrm{3}}}^{\mathrm{)}}}{{}^{\mathrm{\right(}}\frac{{\mathit{v}}_{\mathrm{turb}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{\left)}}}^{-0.32}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ and the analytic solution for Ṁ_∗(t) and Ė_inj becomes $\begin{array}{ccc}\mathit{Ṁ}\mathrm{\ast }\mathrm{\left(}\mathit{t}\mathrm{\right)}& & \mathrm{=}\frac{\mathrm{2}\mathit{\pi }\mathrm{\Psi }\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{BH}}}}{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }}{\mathrm{\Sigma }}_{\mathrm{o}}^{\mathit{\theta }}{\left(\frac{{\mathit{v}}_{\mathrm{o}}}{{\mathit{c}}_{\mathrm{s}}}\right)}^{\mathrm{-}\mathit{ϵ}}\\ & & \mathrm{\times }{\mathit{r}}_{\mathrm{c}}^{\mathit{\theta \gamma }\mathrm{+}\mathit{ϵ\lambda }}\left[{\mathit{r}}_{\mathrm{out}}^{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }}\mathrm{-}{\mathit{r}}_{\mathrm{c}}^{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }}\right]\mathit{,}\mathrm{and}\end{array}$(10)$\begin{array}{ccc}\mathit{Ė}\mathrm{inj}& \mathrm{=}& \mathrm{7.9}\mathrm{\times }{\mathrm{10}}^{-3}{\mathit{\eta }}_{\mathrm{SN}}{\mathit{E}}_{\mathrm{SN}}\frac{\mathit{Ṁ}\mathrm{\ast }}{\mathrm{2}\mathit{\pi }}\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }\right)\\ & & \mathrm{\times }{\mathit{r}}^{\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}}{\left[{\mathit{r}}_{\mathrm{out}}^{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }}\mathrm{-}{\mathit{r}}_{\mathrm{in}}^{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }}\right]}^{-1}\mathit{.}\end{array}$(11)If the gravitational force is dominated by the CND itself, then $\mathrm{\Sigma ̇}{}_{\mathrm{\ast }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{\Psi }\sqrt{{}^{\mathrm{(}}\frac{\mathit{\pi G}{\mathrm{\Sigma }}_{\mathrm{disk}}}{\mathit{r}}^{\mathrm{)}}}{{}^{\mathrm{\right(}}\frac{{\mathit{v}}_{\mathrm{turb}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{\left)}}}^{-0.32}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$, and the analytic solution for Ṁ_∗(t) and Ė_inj follows as $\begin{array}{ccc}\mathit{Ṁ}\mathrm{\ast }\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{2}\mathit{\pi }\mathrm{\Psi }\sqrt{\mathit{G\pi }{\mathrm{\Sigma }}_{\mathrm{o}}}}{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\gamma }\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}{\mathit{\theta }}^{\mathrm{)}}\mathrm{-}\mathit{ϵ\lambda }}{\mathrm{\Sigma }}_{\mathrm{o}}^{\mathit{\theta }}{\left(\frac{{\mathit{v}}_{\mathrm{o}}}{{\mathit{c}}_{\mathrm{s}}}\right)}^{\mathrm{-}\mathit{ϵ}}{\mathit{r}}_{\mathrm{c}}^{\mathit{\theta \gamma }\mathrm{+}\mathit{ϵ\lambda }\mathrm{+}\frac{\mathit{\gamma }}{\mathrm{2}}}\\ & & \mathrm{\times }\left[{\mathit{r}}_{\mathrm{out}}^{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\gamma }\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\mathit{\theta }\right)\mathrm{-}\mathit{ϵ\lambda }}\mathrm{-}{\mathit{r}}_{\mathrm{c}}^{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\gamma }\mathrm{(}\frac{\mathrm{3}}{\mathrm{2}}\mathrm{+}\mathit{\theta }\mathrm{)}\mathrm{-}\mathit{ϵ\lambda }}\right]\mathit{,}\mathrm{and}\\ \mathit{Ė}\mathrm{inj}& \mathrm{=}& \mathrm{7.9}\mathrm{\times }{\mathrm{10}}^{-3}{\mathit{\eta }}_{\mathrm{SN}}{\mathit{E}}_{\mathrm{SN}}\frac{\mathit{Ṁ}\mathrm{\ast }}{\mathrm{2}\mathit{\pi }}\left(\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\gamma }\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\mathit{\theta }\right)\mathrm{-}\mathit{ϵ\lambda }\right)\\ & & \mathrm{\times }{\mathit{r}}^{\mathrm{-}\mathit{\gamma }\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\mathit{\theta }\right)\mathrm{-}\mathit{ϵ\lambda }\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}{\left[{\mathit{r}}_{\mathrm{out}}^{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}\mathrm{+}\mathit{\theta }\right)\mathrm{-}\mathit{ϵ\lambda }}\mathrm{-}{\mathit{r}}_{\mathrm{in}}^{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\gamma }\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\mathit{\theta }\right)\mathrm{-}\mathit{ϵ\lambda }}\right]}^{-1}\mathit{.}\end{array}$Finally, if gravity is dominated by the host galaxy, the star formation rate is expressed as $\mathrm{\Sigma ̇}{}_{\mathrm{\ast }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{\Psi }\sqrt{{}^{\mathrm{(}}\frac{\mathit{\pi G}{\mathrm{\Sigma }}_{\mathrm{host}}}{\mathit{r}}^{\mathrm{)}}}{{}^{\mathrm{\right(}}\frac{{\mathit{v}}_{\mathrm{turb}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{c}}_{\mathrm{s}}}^{\mathrm{\left)}}}^{-0.32}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$, and the analytical solution for Ṁ_∗(t) and Ė_inj becomes $\begin{array}{ccc}\mathit{Ṁ}\mathrm{\ast }\mathrm{\left(}\mathit{t}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{2}\mathit{\pi }\mathrm{\Psi }\sqrt{\mathit{G\pi }{\mathrm{\Sigma }}_{\mathrm{host}}}}{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }}{\mathrm{\Sigma }}_{\mathrm{o}}^{\mathit{\theta }}{\left(\frac{{\mathit{v}}_{\mathrm{o}}}{{\mathit{c}}_{\mathrm{s}}}\right)}^{\mathrm{-}\mathit{ϵ}}{\mathit{r}}_{\mathrm{c}}^{\mathit{\theta \gamma }\mathrm{+}\mathit{ϵ\lambda }}\\ & & \mathrm{\times }\left({\mathit{r}}_{\mathrm{out}}^{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }}\mathrm{-}{\mathit{r}}_{\mathrm{c}}^{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }}\right)\mathit{,}\mathrm{and}\\ \mathit{Ė}\mathrm{inj}& \mathrm{=}& \mathrm{7.9}\mathrm{\times }{\mathrm{10}}^{-3}{\mathit{\eta }}_{\mathrm{SN}}{\mathit{E}}_{\mathrm{SN}}\frac{\mathit{Ṁ}\mathrm{\ast }}{\mathrm{2}\mathit{\pi }}\left(\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }\right)\\ & & \mathrm{\times }{\mathit{r}}^{\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}}{\left[{\mathit{r}}_{\mathrm{out}}^{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }}\mathrm{-}{\mathit{r}}_{\mathrm{in}}^{\frac{\mathrm{3}}{\mathrm{2}}\mathrm{-}\mathit{\theta \gamma }\mathrm{-}\mathit{ϵ\lambda }}\right]}^{-1}\mathit{.}\end{array}$In the most general case, different regimes (where different terms in the gravitational field may become dominant in different parts of the CND) and their evolution may further be time-dependent. As we aim here to work with a still simplified model, we will thus show in the following sections that the results are not highly sensitive to the choice of the gravitational force term.

3. Results

3.1. Evolution of the SMBH accretion rate, the star formation rate and masses

The model described in the previous sections allows us to explore the evolution of different physical quantities of the system that include the evolution of the black-hole growth from an initial seed, the evolution of the star formation rate occurring on the outer disk, the evolution of the gas velocity as well as the AGN luminosity. In general, in this model we adopt a marginally stable self-gravitating disk Q ≈ 1 (consistent with Wutschik et al. 2013) and the universal fraction of molecular clouds (f_GMC = 0.5 as shown by Krumholz & McKee 2005).

As mentioned in the previous section, we examine the impact of the star formation model KM05 on the evolution of the system. In Fig. 1 we compare the final black-hole mass and final stellar mass as a function of the mass supply rate M_⊙ under different assumptions for the dominant gravitational source, using simulations with the input parameters displayed in Table 1. It can be seen that the final stellar masses are independent on the assumptions regarding the gravitational field. The black-hole mass depends somewhat on the underlying gravitational field, but the difference amounts to only a factor of a few. Additionally, in Fig. 1 we compare our results produced with the KM05 model to the results produced with the EB10 model.

We note here that by setting an initial black-hole seed of 10³M_⊙ and having a Ṁ_sup ≥ 100M_⊙ yr^-1 from the galaxy to the circumnuclear disk, both models KM05 and EB10 enable the growth of supermassive black-hole mass to the order of 10⁹M_⊙. It is interesting to note that the final black-hole masses converge to similar values both in the KM05 and EB10 models, here explored in the limit where the CND dominates the gravitational field.

Although both the KM05 and EB10 models produce similar final masses, model KM05 produces more stable accretion during the high accretion phase whereas model EB10 produces over-accretion during that phase (see Fig. 2). Stable accretion is reflected as less sharp fluctuations on the accretion rate and on the star formation rate in the period of 10⁶ to 10⁷ yr. This effect is also reflected in the turbulent Mach number, which is driven by the feedback from star formation. Compared to the EB10 model, the dependence on the turbulent velocity is however less severe, scaling as M^-0.32 rather than M^-1, thus reducing the impact of the negative feedback and damping the fluctuations within the system. In addition, as already described by Wutschik et al. (2013), there are hints of anticyclic variations between the black-hole mass accretion rate and the star formation rate, which may be interpreted as competition between gas consuming processes in the CND.

We have also investigated the influence of the black-hole seed mass on the evolution of the system. For this purpose, Table 2 compares the final black-hole masses at adopting different seeds and under different assumptions regarding the central source of gravity. In all cases, it is clearly visible that the final black-hole mass has an extremely weak dependence on the mass of the initial seed, even when considering seed masses of 10⁶M_⊙, which are at the upper limit of what is suggested in the literature. This effect has been reported already by Armijo & Pacheco (2011) and Wutschik et al. (2013), the latter of which adopts different assumptions on star formation. However, as shown in Fig. 1, a large mass supply rate to the system indeed enhances the final black-hole masses. For instance, for a supply rate of 100 M_⊙ yr^-1, the final masses are of the order of 10⁹M_⊙ which is in agreement with the results given by Volonteri et al. (2015). While for moderate mass supply rates of 1M_⊙ yr^-1 and 10 M_⊙ yr^-1, the final black-hole masses only weakly depend on our assumptions concerning the gravitational acceleration, they differ by roughly a factor of three when considering large supply rates of 100 M_⊙ yr^-1.

3.2. Evolution of the Mach number

Figures 33–5 show the time evolution of the Mach number for different mass supply rates under different assumptions of the central gravitational field. We note in general that the turbulent velocity of the gas becomes larger compared to the local sound speed during the high accretion phase (starting at 10⁶ yr), which leads to large Mach numbers. When comparing the Mach number evolution in the different dominant gravitational scenarios, we notice that it is highest when the gravitational field is dominated by the central black-hole, with M_max ≈ 65 at Ṁ_sup =1 M_⊙ yr^-1. In this scenario, the black-hole accretes efficiently with Ṁ_BH ≈ 1 M_⊙ yr^-1 (see also Fig.11).

When the gas supply ceases after 10⁸ yr, the gas accretion becomes less efficient, so the black-hole accretes afterwards at very low rates (see Fig. 11). At this timescale a very pronounced spike appears on the Mach number evolution at all field dominations. The latter is a result of delayed SN feedback at a time where the supply of new gas is inhibited, thus strengthening the effect of the feedback. Subsequently, on a timescale comparable to the lifetime of massive stars, the Mach number drops to a value of order unity, as turbulence is no longer driven by SN, but may result from hydrodynamical instabilities.

In the simulations presented here, the value of the viscosity coefficient α is of the order of unity in all our simulations. However, this quantity is still uncertain and it is worth mentioning that Kawakatu & Wada (2008) and Wutschik et al. (2013) have associated different values of α to the different accretion phases. In their assumptions, α has the order of unity in the turbulent supersonic phase whereas in the subsonic phase α ≈ 0.01–0.5 (due to magneto-rotational instabilities) where the turbulent velocity is comparable to the local sound speed. To study the dependence of our results on α, we performed simulations to compare different values of α with α = 0.01,0.5,1. Figure 6 shows the time evolution of the Mach number for the different values of α. It can be seen that small values of α< 1 have a weak impact on the accretion phases, i.e. the gas velocity remains comparable to the sound speed during the subsonic phase. Additionally, the Mach number evolution follows the same trend for all values of α except that it reduces very slightly the Mach numbers in the periods of one million and three million years. We find here that the time evolution of the BH mass accretion rate, star formation rate, final masses and radius are quite independent of the parameter α. Such independence on α shows that the accretion is dynamically regulated by the available mass reservoir and therefore remains efficient as long as sufficient matter will be available.

3.3. Evolution of the CND radii

Returning to the evolution of the accretion disk itself, we have also studied the evolution of the final inner and outer radius as a function of time and mass supply rate (see Figs. 7 to 10), exploring different assumptions regarding the gravitational field. The explored mass supply rates, ranging from 1 to 1000 M_⊙ yr^-1, allow the rapid growth of the disk radii from parsec to kiloparsec scale. We find no significant dependence of the final values of the outer radius on the assumption regarding the gravitational field. So the results converge to similar values after the gas supply ceases at 10⁸ yr. The final outer radius is comparable, for instance, with the size of some circumnuclear rings in galaxies such as NGC 1097 with 1.5 kpc of radius studied by Hsieh et al. (2011) and NGC 613 with a nuclear starburst of 700 pc wide found by Falcon-Barroso et al. (2014). The inner radius, on the other hand, always remains at similar values with distances from the centre spanning from 1 to approximately 30 pc.

3.4. Evolution of the luminosities

We also explore here the implications of our model on observable quantities, including the luminosity from the AGN (L_AGN), the nuclear starburst (L_NSB), the stellar luminosity (L_stellar) and for comparison the Eddington luminosity (L_edd). The AGN luminosity is calculated based on two well-known cases: the standard thin disk (Shakura & Sunyaev 1973) which is applied for the low accretion phase and the slim disks (Abramowicz et al. 1988; Abramowicz 2004) for the high accretion phase. Similarly to the method in Kawakatu & Wada (2008, 2009), Eqs. (16) and (17) show the dependence of the luminosities on the accretion rate. In the standard disk aproximation, L_AGN is given as ${\mathit{L}}_{\mathrm{AGN}}\mathrm{=}\left(\frac{\mathit{ṁ}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathrm{10}}\right){\mathit{L}}_{\mathrm{edd}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}\hspace{0.5em}\mathit{ṁ}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{<}\mathrm{20}\mathit{,}$(16)and in the slim disk approximation L_AGN is given as ${\mathit{L}}_{\mathrm{AGN}}\mathrm{=}\mathrm{2}\left(\mathrm{1}\mathrm{+}\ln \frac{\mathit{ṁ}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathrm{20}}\right){\mathit{L}}_{\mathrm{edd}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}\hspace{0.5em}\mathit{ṁ}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\ge }\mathrm{20}\mathit{,}$(17)where $\mathit{ṁ}\mathrm{=}\frac{\mathit{Ṁ}{\mathrm{BH}}_{}}{\mathit{Ṁ}{\mathrm{edd}}_{}}$, and $\mathit{Ṁ}{\mathrm{edd}}_{}\mathrm{=}\mathit{Ṁ}{\mathrm{crit}}_{}\mathrm{=}{\mathit{L}}_{\mathrm{edd}}\mathit{/}{\mathit{c}}^{\mathrm{2}}$ (Watarai et al. 2000). The Eddington luminosity L_edd is given by the following expression: ${\mathit{L}}_{\mathrm{edd}}\mathrm{=}\mathrm{1.25}\mathrm{\times }{\mathrm{10}}^{\mathrm{38}}\frac{{\mathit{M}}_{\mathrm{BH}}}{{\mathit{M}}_{\mathrm{\odot }}}\mathrm{[}\mathrm{erg}\mathit{/}\mathrm{s}\mathrm{]}\mathit{.}$(18)For calculating the nuclear starburst luminosity (L_NSB), we adopt the formula given by Kawakatu & Wada (2008), with ϵ = 0.007: ${\mathit{L}}_{\mathrm{NSB}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{0.14}\mathit{ϵ}\mathit{Ṁ}\mathrm{\ast }\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathit{c}}^{\mathrm{2}}\mathit{.}$(19)Since our model includes phases without active star formation, we estimate the total stellar luminosity by adopting the mass-to-light ratio (M/L) from the total stellar mass. Typical values range from two to ten solar masses per solar luminosity (Dickel 1978; Faber & Gallagher 1979). For our calculations we adopt a mass-to-light-ratio of 6 M_⊙/L_⊙.

The curves shown in Fig. 12 are the results of our simulations for an initial black-hole seed mass of 10³M_⊙, with a mass supply rate of 1 M_⊙ yr^-1. The rest of the input parameters are listed in Table 1. Thus Fig. 12 compares the time evolution of the AGN luminosity, the Eddington luminosity, the nuclear starburst luminosity and the total stellar mass using model KM05 and model EB10. The AGN luminosity is calculated assuming gravity to be dominated by the central black-hole. We note that in both models the time evolution of the Eddington luminosities are similar and by comparing them with AGN luminosites, we identify different timescales of the accretion phase: the high accretion phase (10⁶–10⁸ yr) and the low accretion phase (>10⁸ yr). We further note that the luminosities should more substantial fluctuations in the model EB10 compared to KM05, which is consistent with the results obtained for the star formation and accretion rates.

These results show that our model predicts upper limits of the luminosities: the highest AGN luminosity (~ 10⁴⁵ erg/s) and the lowest AGN luminosity (~ 10⁴² erg/s) as well as the highest nuclear starburst luminosity (~ 10⁴³ erg/s), which are in agreement with Kawakatu & Wada (2008). Furthermore, we note that our model predicts super-Eddington luminosities for the period between 10⁶ yr and 2 × 10⁷ yr. Sub-Eddington luminosities occur when the AGN luminosity starts to decrease in the period of 2 × 10⁷ yr to 10⁸ yr, then the AGN luminosity drops drastically to lower values after 10⁸ yr.

In general, we find no significant differences in the time evolution of the total stellar luminosities and Eddington luminosites when models KM05 and EB10 are applied. However, the time evolution of the AGN and nuclear starburst luminosities in both models differs during the high accretion phase. As we have described in Sect. 3.1 model EB10 produces over-accretion during the high accretion phase so in consequence these effect also produces the variabilites on the luminosities displayed in Fig. 12.

We want to emphasise that the typical timescale for activity in AGN depends on external parameters, and for instance, the external mass supply could last for a shorter time or longer time than 10⁸ yr. This would depend to some extend on the type of the AGN, the redshift and also the environment. For instance, for local AGNs, it may be possible that the accretion phases are typically shorter, while there might be rather long and extended periods of accretion necessary to explain high redshift quasars.

4. Application of the model to active galaxies and starburst rings

4.1. High redshift quasars

Mortlock et al. (2011) discovered the most distant luminous quasar ULAS J1120+0641, which resides at z = 7.085 (7.6 × 10⁸ yr after the big bang). The estimated mass of the supermassive black-hole hosted by that quasar is 2 billion solar masses which generates a very large luminosity of 2.4 × 10⁴⁷ erg/s. Another distant quasar at z = 6.42, J1148+5251 was discovered by Willott et al. (2003), with a supermassive black-hole mass of 4.9 × 10⁹M_⊙ and a luminosity of 5.2 × 10⁴⁷M_⊙ (Pezzulli et al.2016, with references therein). We adapt here our model to reproduce the observed AGN properties of these two high redshift quasars. First, we extended the time of mass supply to 8 × 10⁸ yr, and run simulations for mass supply rates of 30 M_⊙ yr^-1 and 100 M_⊙ yr^-1. The rest of the parameters correspond to the values given in Table 1.

As a result, our model produces final black-hole masses of 2.2 × 10⁹M_⊙ and 4.91 × 10⁹M_⊙ when the gas supply ceases after 8 × 10⁸ yr (see Fig. 13). Such large mass supply rates allow efficient accretion up to 10 M_⊙ yr^-1 (see Fig. 14). This large accretion rate lasts only for a short time in the high accretion phase and then produces an upper limit in the luminosity of almost ~ 10⁴⁷ erg/s at roughly 10⁷ yr (see left graph of Fig. 15). Shortly after the maximum luminosity is reached, the quasar emits continuous radiation at a rate of ~ 10⁴⁶ erg/s until the gas supply ceases.

We find again that our model predicts super-Eddington and sub-Eddington luminosities even when the time supply is extended up to one billion years. These results are displayed in Fig. 15, where we compare the time evolution of the luminosities using model KM05 and model EB10 in the same fashion as in Fig. 12. We note that the application of model EB10 for large accretion rates (>1 M_⊙ yr^-1) reduces the over-accretion, but it still produces variations of the AGN and nuclear starburst luminosities during the high accretion phase.

Figure 14 shows that the star formation rate is the preferred mode of gas consumption in the high accretion phase. This effect can be seen in the period of 10⁶ yr (SN feedback timescale) and 10⁹ yr, where the evolution of the star formation rate comprises larger and enhanced rates than the BH mass accretion rate. This behaviour can be attributed to the very large mass supply rate injected from the host galaxy to the outer disk of the CND, where the stars are forming. It is important to note that the injection of a large mass supply rate produces little gas competition between the BH mass accretion rate and the SFR, whereas for low supply rates, as it was shown in the previous section, a high competition for gas consumption and strong fluctuations develop as displayed in Fig. 11. Hence, our model predicts an upper limit of the SFR value, which is close to 100M_⊙ yr^-1 (see Fig. 14). The observed values of the SFR in the quasars at z = 7.085 and z = 6.42 are 160–440 M_⊙ yr^-1 and 2000 M_⊙ yr^-1 respectively, which were estimated both from far-infrared (FRI) luminosities reported by Venemans et al. (2012) and Pezzulli et al. (2016) (with references therein). While the comparison of the SFR values predicted by our model with the observed results shows an apparent discrepancy, one should note that the quoted values of the SFR refer to the entire host galaxy, while our model only follows the evolution in the CND. Higher-resolution observations may however help to constrain the SFR in the central region of the galaxy in the future.

4.2. The prominent starburst ring in NGC 6951

NGC 6951 is a Seyfert 2 galaxy (z = 0.004750) with a barred spiral morphology. It exhibits a characteristic starburst ring with radius of 580 pc (van der Laan et al. 2013) but other studies have claimed a radius of 480 pc (Sani et al. 2012). Van der Laan et al. (2013) has determined the stellar content on the circumnuclear rings that spans from young to old stars. The stellar clusters along the ring contain predominantly massive stars with intermediate ages. The total stellar mass estimated in the ring is 6 × 10⁷M_⊙ with an estimated age between 1 Gyr and 1.5 Gyr. No other determination of the age of the ring has been published so far.

The mass of the central black-hole ranges between 5.5 × 10⁶M_⊙ at 81° and 1.3 × 10⁷M_⊙ at 33° (Beifiori et al. 2012). As the angle under which we see the system is still uncertain, the unknown angle leads to an uncertainty in the dynamical quantities, in spite of the position of the AGN nucleus and its H_β luminosities having been accurately measured by Perez et al. (2000).

The presence of a bar in the galaxy can potentially affect the formation of a circumnuclear disk and partly influence its evolution. As previously shown by Athanassoula (1992) and Regan & Teuben (2003), the existence of a circumnuclear ring inside a large-scale bar depends on the orbit families that make up the circumnuclear region. While the x1 orbits are elongated along the asymmetric potential and make up the stellar bar, the x2 orbits are oriented perpendicular to the bar, and gas-orbits turn slowly from one into the other within a bar potential. While we cannot exclude that this has an impact on the disk, our previous analysis in Sect. 3 has shown that both the stellar mass and the black-hole mass are relatively insensitive to the details of the potential. With these caveats in mind, we therefore tentatively compare our model to NGC 6951, knowing that a description of the full dynamics should also consider the bar in more detail.

Here we apply and adapt our model to reproduce the properties of the ring in NGC 6951 related to the stellar mass, the radius of the ring, the evolution of the star formation rate as well as the central SMBH mass. In our simulations, we study the evolution of the CND with an initial black-hole seed of 10³M_⊙. The time for mass supply from the host galaxy is extended up to a billion years. We have adopted a low surface density of the host galaxy as 200 M_⊙ pc^-2 and a low constant mass supply rate of 0.2 M_⊙ yr^-1. In addition, we have run the simulations for different timescales associated to SN explosions produced by nuclear starbursts (massive stars). Here we explore the effect of varying the average lifetime of the massive stars, which may be considered as revealing the effect of different stellar populations ending their life at different times. In particular, we adopt possible lifetimes of 5 × 10⁶ yr (21 M_⊙), 10⁷ yr (15 M_⊙) and 5 × 10⁷ yr (8 M_⊙). With these assumptions, our model is able to reach good agreement with the observed properties of the circumnuclear ring in NGC 6951. Table 3 shows the final results of the BH mass, mass of the ring and the ring radius using different assumptions for the central gravitational field. These results are in a good agreement with the observed values reported by van der Laan et al. (2013) and Beifiori et al. (2012).

Figure 16 displays the star formation rate history in the ring, where the stellar ages range from million to billion years. The simulated star formation rate versus time is shown in Fig. 17. For producing all different coloured curves illustrated in Fig. 17, we have employed the input parameters as mentioned above. The coloured curves correspond to the different lifetimes of massive stars exploding into SN. All curves have the same final results for the masses and for the ring radius when employing the same assumptions on gravity as displayed in Table 3. Due to the delayed SN feedback for larger lifetimes we can see that the SFR initially remains higher (>0.02 M_⊙yr^-1; see Fig. 17) until the feedback sets in. In general, such fluctuations in the SFR induced by SN feedback are in good agreement with the observed star formation rate variation (Fig. 16) by van der Laan et al. (2013), implying that our model can reproduce the observed star formation rate oscillating between 0 and 0.2 M_⊙yr^-1.

5. Discussion and conclusions

We have extended the co-evolution model of Wutschik et al. (2013) by implementing the turbulence-regulated star formation model of Krumholz & McKee (2005). We have additionally studied the time evolution of the Mach numbers as well as the AGN and stellar luminosities. The implementation of the KM05 model produces moderate fluctuations in the black-hole mass accretion rates and star formation rates during the accretion phase rather than the strong fluctuations produced by the model EB10. The latter predicts over-accretion, i.e. the sonic speed becomes extremely large compared to the local sound speed which produces, in consequence, extreme Mach numbers up to 10⁴.

We have studied the limiting cases of the total star formation rate under different assumptions for the gravitational field: dominated by the central black-hole, the CND or the host galaxy. For each of these cases, our model is able to predict the growth of a SMBH from 10 million solar masses up to a few billions solar masses, with low (1 M_⊙ yr^-1) and very large (1000 M_⊙ yr^-1) mass supply rates respectively, in agreement with earlier work by Wutschik et al. (2013), Kawakatu & Wada (2009) and Volonteri et al. (2015). On the other hand, though that model EB10 produces an over-accretion, the final black-hole masses are similar to the final values provided by model KM05 for gravity dominated by the CND. The final stellar masses generated by both models are also identical when compared at all field dominations.

An increased mass supply rate enhances black-hole mass accretion rate, star formation rate, disk radius final, black-hole masses, final stellar masses as well as luminosities. However, we note that large mass supply rates dampen the level of turbulence and the competition of gas accretion during high accretion phases. This effect becomes more evident on the Mach number evolution when gravity is dominated by the central black-hole. Large mass supply rates also significantly affect the disk radius; they lead to the growth of the inner and outer radii, from parsec scale to a few kilo-parsec. For instance, to produce an SMBH of a billion solar masses, our model predicts a final inner radius of 10 pc and outer radius of 500 pc. Moreover, our model shows that large mass supply rates have not only a significant impact on the final black-hole mass or black-hole accretion rates but also on the order of magnitude of AGN luminosities.

Since our model considers high and low accretion phases during the evolution of the system, we have also assumed that the time evolution of the AGN luminosities follows similar phases accordingly. Adopting the slim disk approximation during the high accretion phase our model predicts super-Eddington luminosities, whereas sub-Eddington luminosities are found during the low accretion phase when the standard thin disk approximation is applied. The latter is reflected in the decrease of the luminosity after the end of the mass supply. However, afterwards it seems that the central black-hole keeps accreting at very low rates of the order of 10^-3M_⊙ yr^-1.

As a result of the extended supply time and increased mass supply rate, our model is able to predict an upper limit on the AGN luminosity which is close to 10⁴⁷ erg/s. AGN luminosities lower than 10⁴⁷ erg/s are in agreement with Kawakatu & Wada (2008, 2009). However, the upper limit is reached at a different time than the cosmological age of the observed AGN by Mortlock et al. (2011) and Pezzulli et al. (2016) (with references therein). Nevertheless, we do not know when the mass supply started in the observed systems, and of course the real mass supply evolution may be more complex, i.e. time dependent. Rapid variability on larger scales caused by chaotic cold accretion (CCA; Gaspari et al. 2017; Gaspari & Sa¸dowski 2017) could also influence the mass supply rate. However, despite this possible discrepancy with respect to time scales, our model succeeds in predicting the observed final black-hole masses.

We find that turbulence is a main driver of angular momentum redistribution (Balbus & Hawley 1998; Klessen & Hennebelle 2010; Gaspari et al. 2017), and particularly high values of turbulence are found in the phase of high star formation rates and for the gravitational field dominated by the central black-hole. We note that the turbulent velocity of gas becomes larger not only during the high accretion phase but also when the gas supply is depleted.

While energy is injected by SN already at earlier stages of the accretion, the SN explosions occuring at the later stages are associated with the rise of the turbulent velocity due to a decreasing density, which is reflected in the peaks of the Mach number in the late evolutionary stages. We emphasize here that the timescale of SN feedback sets the timescale for fluctuations in both the black-hole mass accretion rate and star formation rate.

We emphasise that our model does not take into account any feedback from the AGN but if considered, it may indeed decrease the accretion rates and the jet emission may trigger star formation as investigated in MHD simulations by Gaibler et al. (2012). Overall, our model does neither assume a specific structure of the galaxy nor take into account the interaction with the bar in NGC 6951, some quantities like the black-hole mass and the stellar mass can still be reproduced. This suggests that these quantities may not strongly depend on the choice of the gravitational potential, though further investigations may be necessary in the future.

Acknowledgments

We thank the anonymous referee for insightful comments regarding our manuscript. We are grateful to Tessel van der Laan for permission to include Fig. 16 in our manuscript. D.R.G.S. thanks for funding through Fondecyt regular (project code 1161247), through the “Concurso Proyectos Internacionales de Investigación, Convocatoria 2015” (project code PII20150171), trough ALMA-Conicyt (project code 31160001) and from the Chilean BASAL Centro de Excelencia en Astrofísica y Tecnologías Afines (CATA) grant FB-06/2007. W.C. wishes to thank to M. Breuhaus for his encouragements and helpful discussions.

References

All Tables

All Figures