Outflows in the Disks of Active Galaxies

In the last decade, a wealth of observations has provided an increasingly detailed characterization of galaxy-scale outflows driven by active galactic nuclei (AGNs). The early observations of ultrafast AGN-driven winds on small scales (from the accretion disk scale up to the dusty torus) with velocities ≈0.1c (see King & Pounds 2015 for a review) through blueshifted absorption lines in the X-ray spectra in a substantial fraction (≈40%) of AGNs (e.g., Piconcelli et al. 2005; Tombesi et al. 2010; Gofford et al. 2013) have been recently complemented by a wide set of measurements of fast (velocities of the order of 1000 km s⁻¹), massive flows of ionized, neutral, and molecular gas, extended on kiloparsec scales. These have been performed through deep optical/NIR spectroscopy (Nesvadba et al. 2006, 2008; Alexander et al. 2010; Riffel & Storchi-Bergmann 2011; Rupke & Veilleux 2011; Cano-Diaz et al. 2012; Greene et al. 2012; Harrison et al. 2012, 2014; Cimatti et al. 2013; Liu et al. 2013a, 2013b; Genzel et al. 2014; Tadhunter et al. 2014; Brusa et al. 2015a, 2015b; Cresci et al. 2015; Perna et al. 2015a; Carniani et al. 2015; Perna et al. 2015b; Zakamska et al. 2016; Bischetti et al. 2017), through interferometric observations in the (sub)millimeter domain (e.g., Feruglio et al. 2010, 2013, 2015; Alatalo et al. 2011; Krips et al. 2011; Aalto et al. 2012; Cicone et al. 2012, 2014, 2015; Maiolino et al. 2012; Combes et al. 2013; Morganti et al. 2013a, 2013b; Garcia-Burillo et al. 2014), and through far-infrared spectroscopy from Herschel (e.g., Fischer et al. 2010; Sturm et al. 2011; Spoon et al. 2013; Veilleux et al. 2013; Stone et al. 2016; Gonzalez-Alfonso et al. 2017). These observations have enabled the detailed physical properties of the outflows (velocities, mass outflow rate, kinetic energy rate) for a number of sources with different AGN luminosity and host galaxy properties (gas mass, circular velocity) to be determined. Recent works by Cicone et al. (2014) and Fiore et al. (2017) have allowed samples with more than a hundred outflow measurements with detected massive winds at different scales (subparsecs to kiloparsecs) and with different molecular/ion compositions to be assembled. For several molecular outflows, the complementary measurement of the host galaxy circular velocity and gas mass has been used to constrain the relationships among the wind parameters, AGN parameters, and host galaxy parameters.

Parallel theoretical works (Silk & Rees 1998; King 2003; Granato et al. 2004; Lapi et al. 2005; Silk & Nusser 2010; King et al. 2011; Faucher-Giguère & Quataert 2012) have focused on capturing the main features of the outflows and on pinning down their main expansion and cooling properties, mainly through the implementation of models based on shocks expanding into the interstellar medium (ISM) approximated as a sphere with a power-law density profile $\rho \sim {R}^{-\alpha }$ (for the extension to exponential disks, see Hartwig et al. 2018). Within the large uncertainties and approximations, energy-conserving shock models are consistent with present measurements that indicate that AGN-driven, galaxy-scale outflows may commonly have momentum fluxes $10\,{L}_{\mathrm{AGN}}/c$ (in terms of the AGN bolometric luminosity L_AGN). These models allowed the derivation of scaling laws for the run of the shock velocity V_s, for the associated mass outflow rate ${\dot{M}}_{S,\theta }$, and for their dependence on the AGN luminosity, e.g., for the case of an isothermal sphere (α = −2), models including cooling predict mass outflow rates ${\dot{M}}_{s}\sim {L}_{\mathrm{AGN}}^{1/3}$, while the energy-conserving model by Lapi et al. (2005) yields a slightly steeper dependence ${\dot{M}}_{s}\sim {L}_{\mathrm{AGN}}^{1/2}$. The observed steeper dependencies ${\dot{M}}_{s}\sim {L}_{\mathrm{AGN}}^{0.8}$ for molecular winds and ${\dot{M}}_{s}\sim {L}_{\mathrm{AGN}}^{1.3}$ for ionized winds (see Fiore et al. 2017) then point toward a medium where the density profile is flatter than in the isothermal case (Faucher-Giguère & Quataert 2012), although such conclusions may be affected by biases in the observational results.

Although refined, recent shock models have started to compare with the distribution of observational outflow measurements (see, e.g., Richings & Faucher-Giguère 2018b), the increasing wealth of data concerning the physical properties of a large number of AGN-driven outflows calls for a more detailed and quantitative comparison with models, starting with a "one-by-one" comparison of model predictions with the measured outflow properties in well-studied objects residing in galaxies with different measured gas and dynamical mass.

Toward this aim, we extend the shock model for AGN outflows to include realistic exponential density profiles for the ISM, where the normalization of the gas density is related to the global gas content of the host galaxy disk, and the disk scale radius is related to the total host galaxy mass. This allows us to compare the shock model results with the most recent compilations of data concerning AGN outflows with different AGN luminosities and host galaxy gas and dark matter (DM) masses, measured at different distances from the host galaxy center. Our goal is to incorporate most previous advances into a single yet manageable analytic framework so as to describe the expansion of AGN-driven shocks for realistic exponential density profiles for the ISM in a disk geometry, and to derive solutions in terms of the global properties of the host galaxy and of the central AGN. This allows us to perform a direct comparison with existing observations in a variety of galactic hosts with measured properties, and out to a distance of ∼10 kpc from the center. The goal is to provide an observationally based test ground for the current description of AGN-driven shocks in realistic galactic hosts and to assess to what extent the present status of the modeling is consistent with the existing observational distribution of the expansion and mass outflow rates.

While the main observables we compare with, i.e., expansion and mass outflow rates, can be reliably computed using the analytical formalism we adopt (as found in earlier works comparing simulations with analytical computations in well-studied cases, e.g., Richings & Faucher-Giguère 2018b), our analytical approach is complementary to numerical simulations. For example, while a precise description of the position-dependent molecular, ionization, and chemical properties of the shocked shell requires numerical simulations to account for the effects of Rayleigh–Taylor instabilities (Richings & Faucher-Giguère 2018a; Richings & Faucher-Giguère 2018b; see also Zubovas & King 2014), our treatment effectively follows the expansion velocity of the shock and the mass outflow rate out to large radii where the assumption of power-law gas density profiles (adopted in such simulations) fails, and where the disk is the dominant component with respect to the rapidly declining bulge component (the shock expansion in this case is treated, e.g., in King et al. 2011). Also, our analytical model allows us to easily explore the dependence of the AGN-driven outflows on a variety of quantities (including the AGN luminosity, the gas mass fraction, and the total mass of the host galaxy) over a huge range of values. In addition, our computation allows us to describe the two-dimensional structure of the outflow, as opposed to the isotropic situations considered in most simulations (for simulations of outflows in a nonspherical, elliptical distribution of gas, see Zubovas & Nayakshin 2014). In this sense, our approach is similar to that adopted in Hartwig et al. (2018) but focused on the exploration of a wide set of properties of galactic hosts (including the gas mass fraction) and on the systematic comparison with the most recent compilation of observational data encompassing a wide range of properties of the host galaxy and of the central AGN.

The plan of the paper is as follows. We provide the basic equations governing the expansion of AGN-driven shocks in Section 2. In Section 3, we first derive solutions for an isotropic distribution of gas with a power-law density profile (Section 3.1) to compare with previous studies. Then, we derive solutions in the case of exponential gas density profiles along the plane of the disk (Section 3.2.1) and in the other directions (Section 3.2.2). Section 4 is devoted to a detailed comparison with existing samples of observed AGN-driven outflows. In Section 4.1, we compare with outflows in galaxies with measured gas and total mass, allowing for a one-by-one quantitative comparison with the model predictions, while in Section 4.2 we consider a large sample of observed outflows in galaxies where measurements of gas and total mass are not available, so that the comparison has to be performed assuming observational scaling laws for the observed host galaxy properties. In Section 5, we reconsider the AGN wind scaling laws empirically derived by many authors in light of the results from our updated models. Section 6 is devoted to a discussion and our conclusion.

We adopt the standard shell approximation (see Cavaliere & Messina 1976; Ostriker & McKee 1988; King 2003; Lapi et al. 2005; King 2010; Faucher-Giguère & Quataert 2012; Ishibashi & Fabian 2014, 2015; King & Pounds 2015; Hartwig et al. 2018) for the expansion of shocks into the ambient ISM of the host galaxy. We assume that nuclear winds with velocities V_in ≈ 3 × 10⁴ km s⁻¹ generated by the central AGN accelerate a forward shock expanding into the ambient medium. In the general two-dimensional case (see Figure 1), the shock radius ${R}_{S,\theta }$ and velocity ${V}_{S,\theta }={\dot{R}}_{S,\theta }$ depend not only on time, but also on the angle θ between the direction of expansion and the plane of the disk. The shock expansion results in the formation of a shell of swept-up material defined by $R\approx {R}_{S,\theta }$, with a mass outflow rate ${\dot{M}}_{S,\theta }$ in the considered direction, enclosing a bubble of hot shocked medium. We compute the expansion of the bubble and the properties of the shocked shell, in turn.

Zoom In Zoom Out Reset image size

2.1. The Expansion of the Shock

In the shell approximation, we consider the motion of a mass element ${\rm{\Delta }}{M}_{S,\theta }$ of the swept-out gas at a shock radius ${R}_{s,\theta }$ propagating with velocity ${V}_{S,\theta }={\dot{R}}_{S,\theta }$ in the direction defined by the angle θ with respect to the plane of the disk. Within a solid angle ${\rm{\Delta }}{\rm{\Omega }}=2\,\pi \,\cos (\theta ){\rm{\Delta }}\theta$, the mass element ${\rm{\Delta }}{M}_{S,\theta }\,={\rm{\Delta }}{\rm{\Omega }}\,{\int }_{0}^{{R}_{S,\theta }}\,{{dR}}^{2}\,{R}^{2}\,\rho (R,\theta )$ is given by the initial galaxy gas mass (distributed according to a density profile $\rho (R,\theta )$) enclosed within the shock radius ${R}_{S,\theta }$ in the considered direction. We define the mass swept out per unit solid angle ${M}_{S,\theta }^{{\prime} }\equiv {{dM}}_{S,\theta }/d{\rm{\Omega }}$, so that the total mass of the swept-out gas is ${M}_{S}=\int \,d{\rm{\Omega }}{M}_{S,\theta }^{{\prime} }$. Then, the expansion of the shock in the solid angle ΔΩ along the considered direction θ is given by

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d}{{dt}}\left({M}_{s,\theta }^{{\prime} }{\rm{\Delta }}{\rm{\Omega }}\,{V}_{s,\theta }\right) & = & {\rm{\Delta }}{\rm{\Omega }}\,{R}_{s,\theta }^{2}\left({P}_{b}-{P}_{0}\right)\\ & & -\,\displaystyle \frac{G\,{M}_{s,\theta }^{{\prime} }\,{\rm{\Delta }}{\rm{\Omega }}\,M(\lt {R}_{S,\theta })}{{R}_{s,\theta }^{2}}.\end{array}\end{eqnarray} \tag{ 1 }$$

The above equation accounts for the balance between the pressure term (fueling the expansion) acting on the surface element corresponding to the solid angle in the considered direction, θ, and the counteracting gravitational term determined by the total mass M (contributed by the DM and by the central black hole) within the shock radius ${R}_{s,\theta }$. Here, P₀ is the initial ambient pressure, while the pressure of the hot gas in the bubble is ${P}_{b}={E}_{b}/(3/2)Q(t)$, where Q(t) is the volume enclosed by the bubble at the considered time t. In turn, the thermal energy E_b of the bubble evolves according to

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}{E}_{b}=\epsilon \,{L}_{\mathrm{AGN}}-\,{P}_{b}\int d{\rm{\Omega }}\displaystyle \frac{{dS}}{d{\rm{\Omega }}}\,{V}_{S,\theta }-{L}_{\mathrm{cool}}.\end{eqnarray} \tag{ 2 }$$

Here, S(t) is the surface enclosing the hot bubble at the considered time t, $\epsilon \approx {v}_{\mathrm{in}}/c$ is the efficiency of the AGN radiation transferring energy to the ISM medium, while L_cool is the cooling rate of the bubble ${L}_{\mathrm{cool}}=Q(t)({{\rm{\Lambda }}}_{\mathrm{IC}}+{{\rm{\Lambda }}}_{\mathrm{ff}})$, in turn related to the cooling functions for inverse Compton (Λ_IC) and free–free emission (Λ_ff). Notice that we have assumed the inner boundary of the shocked wind bubble to be much smaller than its outer boundary ${R}_{s,\theta }$, an approximation that is known to impact the results by less than 5% (see Richings & Faucher-Giguère 2018b). In the following, we will also ignore the initial pressure P₀ because it is found to be much smaller that the pressure P_b in light of the small initial temperature ∼10⁴ K of the unperturbed medium when compared to the temperatures ${T}_{b}\sim {10}^{9}\mbox{--}{10}^{11}$ K of the bubble (see also Richings & Faucher-Giguère 2018b).

We then define the geometrical factor $C({R}_{S,\theta })\,\equiv Q(t)/(4\pi /3\,{R}_{S,\theta }^{3})$ that expresses the deviation of the volume from the spherical case. With the above notation and within the approximations for the bubble volume and for the initial pressure of the ambient ISM discussed above, we can recast Equation (1) as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\,{V}_{S,\theta }}{{dt}} & = & \displaystyle \frac{2}{C({R}_{S,\theta })}\,\displaystyle \frac{{E}_{b}}{4\,\pi \,{M}_{S,\theta }^{{\prime} }\,{R}_{S,\theta }}\\ & & -\,\displaystyle \frac{M(\lt {R}_{S,\theta })}{M}\,\displaystyle \frac{{R}_{v}}{{R}_{S,\theta }}\,\displaystyle \frac{{V}_{c}^{2}}{{R}_{S,\theta }}-{V}_{S,\theta }\,\displaystyle \frac{{\dot{M}}_{S,\theta }^{{\prime} }}{{M}_{S,\theta }^{{\prime} }}.\end{array}\end{eqnarray} \tag{ 3 }$$

The last term on the right-hand side can be readily computed for a specific assumed density profile. We start from the general form $\rho \equiv {\rho }_{0}\,g(R/{R}_{v},\theta )$, where R_v is the virial radius of the host galaxy. Defining the rescaled radius $x\equiv R/{R}_{v}$, the normalization ρ₀ is related to the total gas content of the galaxy by the relation ${\rho }_{0}={M}_{\mathrm{gas}}/{R}_{v}^{3}\,I$, where the form factor I is obtained by integrating the density profile $g(x,\theta )$ over the volume occupied by the galactic gas. Thus, $I=4\,\pi {\int }_{0}^{1}\,g(x){x}^{2}{dx}$ in the case of a spherical distribution, while $I=4\,\pi {\int }_{0}^{{arctg}(h/x\,{R}_{v})}{\int }_{0}^{1}\,g(x){x}^{2}\,\cos \theta \,d\theta \,{dx}$ in the case of an isotropic distribution inside a disk, with a sharp cutoff at a distance h in the direction perpendicular to the disk. With such a notation, we get

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\,{V}_{S,\theta }}{{dt}} & = & \displaystyle \frac{2}{C({R}_{S,\theta })}\,\displaystyle \frac{{E}_{b}}{4\,\pi \,{M}_{S,\theta }^{{\prime} }\,{R}_{S,\theta }}-\displaystyle \frac{M(\lt {R}_{S,\theta })}{M}\\ & & \times \,\displaystyle \frac{{R}_{v}}{{R}_{S,\theta }}\,\displaystyle \frac{{V}_{c}^{2}}{{R}_{S,\theta }}-\displaystyle \frac{g({R}_{S,\theta }/{R}_{V},\theta )}{{\displaystyle \int }_{0}^{{R}_{S,\theta }/{R}_{v}}g(x,\theta ){x}^{2}\,{dx}}{\left(\displaystyle \frac{{R}_{S,\theta }}{{R}_{v}}\right)}^{3}\,\displaystyle \frac{{V}_{S,\theta }^{2}}{{R}_{S,\theta }},\end{array}\end{eqnarray} \tag{ 4 }$$

where we have defined the total mass (mainly contributed by DM) within the virial radius $M\equiv M(\lt {R}_{v})$, and we have expressed the ratio ${GM}(\lt {R}_{S,\theta })/{R}_{v}={V}_{c}^{2}\,[M(\lt {R}_{S,\theta })/M]$ in terms of the host galaxy circular galaxy velocity ${V}_{c}={GM}/{R}_{v}$. In the following, we will express distances in units of ${R}_{0}=1\,\mathrm{kpc}$, velocities in units of V₀ = 1000 km s⁻¹, and masses in units of ${M}_{0}={10}^{12}\,{M}_{\odot }$. Correspondingly, energies are expressed in units of ${E}_{0}={M}_{0}\,{V}_{0}^{2}$ and time in units pf ${t}_{0}={R}_{0}/{V}_{0}$. After defining $r\equiv R/{R}_{0}$, ${r}_{v}\equiv {R}_{v}/{R}_{0}$, $v=V/{V}_{0}$, $m\equiv M/{M}_{0}$, $e\equiv E/{E}_{0}$, and $\tilde{t}\equiv t/{t}_{0}$, the set of equations defining the expansion of the shock into the ambient ISM is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\,{v}_{S,\theta }}{d\tilde{t}} & = & \displaystyle \frac{2}{C({r}_{S,\theta })}\,\displaystyle \frac{{e}_{b}}{4\,\pi \,{m}_{S,\theta }^{{\prime} }\,{r}_{S,\theta }}-\displaystyle \frac{m(\lt {r}_{S,\theta })}{m}\\ & & \times \,\displaystyle \frac{{r}_{v}}{{r}_{S,\theta }}\,\displaystyle \frac{{v}_{c}^{2}}{{r}_{S,\theta }}-\displaystyle \frac{g({r}_{S,\theta }/{r}_{V},\theta )}{{\displaystyle \int }_{0}^{{r}_{S,\theta }/{r}_{v}}g(x,\theta ){x}^{2}\,{dx}}{\left(\displaystyle \frac{{r}_{S,\theta }}{{r}_{v}}\right)}^{3}\,\displaystyle \frac{{v}_{S,\theta }^{2}}{{r}_{S,\theta }},\end{array}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\,{e}_{b}}{d\tilde{t}} & = & 1.5\,{10}^{-4}\epsilon \,{l}_{\mathrm{AGN}}-\displaystyle \frac{2\,{e}_{b}}{3}\\ & & \times \,\displaystyle \frac{1}{Q(t)}\,\displaystyle \int \,d{\rm{\Omega }}\displaystyle \frac{{dS}(t)}{d{\rm{\Omega }}}\,{v}_{S,\theta }-{l}_{\mathrm{cool}}({e}_{b},{l}_{\mathrm{AGN}}),\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{r}_{S,\theta }={\int }_{0}^{\tilde{t}}d\tilde{t}^{\prime} \,{v}_{S,\theta }(\tilde{t}^{\prime} ),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{m}_{S,\theta }^{{\prime} }={m}_{\mathrm{gas}}{\int }_{0}^{{r}_{S,\theta }/{r}_{v}}\,g(x){x}^{2}\,{dx}/I,\end{eqnarray} \tag{ 8 }$$

where all luminosities ${l}_{\mathrm{AGN}}\equiv L/{L}_{0}$ and ${l}_{\mathrm{cool}}\equiv {l}_{\mathrm{cool}}/{L}_{0}$ are expressed in units of ${L}_{0}={10}^{45}$ erg s⁻¹, and the volume integral in Equation (7) is performed over the regions where the gas is initially distributed and extends up to a rescaled radius ${r}_{S,\theta }/{r}_{v}$. To compute the fraction of total mass within the shock radius $m(\lt {r}_{S,\theta })/m$, we include in $m(\lt {r}_{S,\theta })$ the contributions from both the central black hole mass m_BH and the DM mass m_DM. For the latter, we assume a Navarro et al. (1997) form, ${m}_{\mathrm{DM}}(r)$/${m}_{\mathrm{DM}}({r}_{v})\,=[{ln}(1+{cx})-{cx}$/$(1+{cx})]$/$[{ln}(1+c)-c$/$(1\,+c)]$, where c is the concentration for which we adopt the expression given in Macciò et al. (2008), and $x=r/{r}_{v}$ . The computation of the cooling term requires the bubble temperature T_b and densities n_b. These are computed after Equations (3.7) and (3.8) in Richings & Faucher-Giguère (2018b) assuming a fully ionized plasma with mean molecular weight μ = 14/23 to get ${T}_{b}\,=28/69({m}_{p}\,{E}_{b}\,{V}_{\mathrm{in}}^{2}/{k}_{b}\,\epsilon \,{L}_{\mathrm{AGN}}\,t)$ and ${n}_{b}={L}_{\mathrm{AGN}}\,t\,{X}_{{\rm{H}}}/Q(t){V}_{\mathrm{in}}^{2}\,{m}_{p}$ (here k_b is the Boltzmann constant, m_p is the proton mass, and X_H = 0.7 is the hydrogen mass fraction); the associated electron density is taken to be ${n}_{e}=1.2\,{n}_{b}$. From these quantities, the inverse Compton and free–free cooling functions entering the computation of l_cool are computed using Equations (3.4), (3.5), (B7), and (B8) in Richings & Faucher-Giguère (2018b).

2.2. The Properties of the Shocked ISM Shell

To compute the properties of the shocked ISM shell, we adopt the approach in Richings & Faucher-Giguère (2018b). We first compute the evolution of the energy in the layer as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d}{{dt}}{E}_{s} & = & {P}_{b}\,\int \,d{\rm{\Omega }}\displaystyle \frac{{dS}}{d{\rm{\Omega }}}\,{V}_{S,\theta }-\int d{\rm{\Omega }}\displaystyle \frac{{GM}(\lt {R}_{S,\theta }){M}_{S,\theta }^{{\prime} }}{{R}_{S,\theta }^{2}}\,{V}_{S,\theta }\\ & & -\,{L}_{\mathrm{cool},\mathrm{shell}}({E}_{s},{n}_{s}).\end{array}\end{eqnarray} \tag{ 9 }$$

This accounts for the balance between the work done on the shocked ISM layer by the shocked wind bubble pressure, and the effects of the gravitational potential and of the cooling. The thermal energy of the shell is ${E}_{s,\mathrm{th}}={E}_{s}\,-(1/2)\int \,d{\rm{\Omega }}{M}_{S,\theta }^{{\prime} }\,{V}_{S,\theta }^{2}$ and the associated temperature is ${T}_{s}\,=(28/69){E}_{s,\mathrm{th}}\,{m}_{p}/{M}_{S}\,{k}_{{\rm{B}}};$ it evolves according to Equation (8), where the bubble energy E_b satisfies Equation (5). The associated hydrogen number density in the shocked shell is ${n}_{s}=(3/2){P}_{b}\,{M}_{S}\,{X}_{{\rm{H}}}/{m}_{p}\,{E}_{s,\mathrm{th}}$ (Equation (3.17) in Richings & Faucher-Giguère 2018b), which can be recast (in our usual units defined in Section 2.1) as ${n}_{s}\approx 0.5({e}_{b}/{e}_{s,\mathrm{th}})[{m}_{S}/Q(t)]\,{10}^{4}$ cm⁻³, and the associated electron number density is 1.2 n_s. The above values of temperature T_s and density n_s are used to compute the cooling rate ${L}_{\mathrm{cool},\mathrm{shell}}({E}_{s},{n}_{s})$ from free–free and line emission following Equations (3.11) and (3.12) in Richings & Faucher-Giguère (2018b).

In sum, the following set of equations describes the energy evolution in the shocked ISM shell, in terms of our rescaled variables:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{d\,{e}_{s}}{d\tilde{t}}={e}_{b}\,\displaystyle \frac{2}{3\,Q(t)}\,\displaystyle \int d{\rm{\Omega }}\,\displaystyle \frac{{dS}}{d{\rm{\Omega }}}\,{v}_{s,\theta }\\ \quad -\,\displaystyle \int d{\rm{\Omega }}\displaystyle \frac{m(\lt {r}_{S,\theta })}{m}\,\displaystyle \frac{{r}_{v}}{{r}_{S,\theta }}\,\displaystyle \frac{{v}_{c}^{2}}{{r}_{S,\theta }}\,{v}_{S,\theta }\,{m}_{S,\theta }^{{\prime} }-{l}_{\mathrm{cool},\mathrm{shell}}({e}_{s},{n}_{s}),\end{array}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{e}_{s,{th}}={e}_{s}-(1/2)\int d{\rm{\Omega }}\,{m}_{S,\theta }^{{\prime} }\,{v}_{S,\theta }^{2},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{n}_{s}=0.5({e}_{b}/{e}_{s,{th}})[{m}_{S}/Q(t)]\,{10}^{4}\,{\mathrm{cm}}^{-3},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{T}_{s}=({e}_{s,{th}}/{m}_{S})0.5\,{10}^{8}\,{\rm{K}},\end{eqnarray} \tag{ 13 }$$

where the volume Q(t) is expressed in units of ${R}_{0}^{3}=1$ kpc³. We note that this set of equations is coupled with those describing the expansion of the shock (Equations (4)–(7)) through the bubble energy e_b, the shock position ${r}_{S,\theta }$, and the shock velocity ${v}_{S,\theta }$.

The solutions of Equations (4)–(7) and (9)–(12) depend on the assumed initial density distribution $\rho ={\rho }_{0}\,g(x,\theta )$, where $x=R/{R}_{v}$ and the normalization ${\rho }_{0}$ is related to the total gas content M_gas of the host galaxy (see Section 2). Although we will focus on exponential density profiles in a nonisotropic disk geometry, we first derive solutions for a spherical initial density distribution with a scale-free, power-law dependence on the radius R, as in this case analytical solutions exist in the limit of energy-conserving shock. This allows us to test the reliability of our numerical solutions against analytical results.

3.1. Testing the Numerical Solutions: The Case of Power-law Density Profiles

To test our solutions, we first consider a spherically symmetric initial gas density distribution with a power-law profile $g(x)={x}^{-\alpha }$, and the proper form factor I entering Equation (7) is simply $I=4\pi \,{\int }_{0}^{1}\,{x}^{2-\alpha }{dx}$. In Figures 1(a)–(c), we show our results for different values of the power-law index α, the AGN bolometric luminosity L_AGN, and gas mass M_gas for a host galaxy with DM mass $M={10}^{12}\,{M}_{\odot }$. Due to the spherical symmetry, in this case we have shock solutions that are independent of the inclination θ, i.e., in all equations in Section 2.1, we have ${R}_{S,\theta }={R}_{s}$, ${V}_{S,\theta }={V}_{S}$, C = 1, $Q(t)\,=(4\,\pi /3){R}_{s}^{3}(t)$, and ${dS}/d{\rm{\Omega }}={R}_{S}^{2}(t)$. Notice that in this case the equations can be written in terms of the total swept-up mass ${M}_{S}=4\,\pi \,{M}_{S}^{{\prime} }$.

We explore the dependence of our numerical solutions on the assumed value of α in the left panels of Figure 2. We note that decreasing α corresponds to a faster decline of the velocity V_S and to a steeper increase of the mass outflow ${\dot{M}}_{S}$ as a function of the shock position. This behavior was already found by earlier numerical and analytical works (see Faucher-Giguère & Quataert 2012). Indeed, in the case of a scale-free, power-law density profile, self-similar analytical scalings can be derived for ${R}_{S}\sim {t}^{3/5-\alpha }$ (and hence for ${V}_{S}={\dot{R}}_{S}$) in the limit of negligible cooling (energy-conserving outflows). These self-similar solutions are shown as dashed lines and provide a test for our numerical solutions. The match between numerical and self-similar solutions also indicate the approximate energy-conserving behavior of the expanding bubble (as already found by Faucher-Giguère & Quataert 2012 and Richings & Faucher-Giguère 2018b), although in the center the large gas densities achieved in the α = 2 case result in efficient cooling, yielding the slower shock velocity V_S visible in the central panels of Figure 1(a) compared to the self-similar solution.

Zoom In Zoom Out Reset image size

As a further test for our numerical solutions, we present in the central panels of Figure 1 the scaling of R_S, V_S, and ${\dot{M}}_{S,}$ with the AGN bolometric luminosity L_AGN. The solutions are characterized by increasing the normalization for all such quantities for increasing L_AGN, due to the larger energy injection powering the bubble expansion. Again, we can test our numerical results against self-similar analytical solutions, yielding $\mathrm{log}{R}_{S}\,\sim (1/3)\mathrm{log}{L}_{\mathrm{AGN}}$ (see Faucher-Giguère & Quataert 2012; see also Lapi et al. 2005) for the normalization of the shock expansion, again finding an excellent agreement.

Finally, we study the dependence of our solutions on the total host galaxy mass M_gas (right panels of Figure 1), corresponding to varying the normalization ${\rho }_{0}$ of our assumed density profile. We find that increasing M_gas results in faster shock velocities and smaller mass outflow rates, in agreement with previous works. In this case, self-similar solutions yield ${R}_{S}\sim {\rho }_{0}^{1/(5-\alpha )}\sim {M}_{\mathrm{gas}}^{1/(5-\alpha )}$, again in excellent agreement with our numerical results.

3.2. Solutions for Exponential Gas Density Profiles

Having tested the reliability of our numerical solution, we can proceed toward a detailed comparison between the properties of outflows observed in different galaxies and the predictions of shock models. Toward this aim, we consider a disk geometry (see Figure 2) for the distribution of galactic gas, with a gas density depending only on the galactocentric distance $x=R/{R}_{v}$, but confined within vertical boundaries corresponding to a disk scale height h. This is assumed to be constant with radius for a given galaxy (although for ${M}_{\mathrm{BH}}\lt {10}^{8}\,{M}_{\odot }$, models predict the gravitational bending of the interstellar gas below 100 pc, due to the black hole gravitational field; see Lamastra et al. 2006), and to increase with the galaxy circular velocity according to the observed average relation $h=0.45({V}_{c}/100\,\mathrm{km}\,{{\rm{s}}}^{-1})-0.14\,\mathrm{kpc}$ (see van der Kruit & Freeman 2011 and references therein). Inside the disk (where the vertical distance from the plane of the disk Y is smaller than the scale height h), we adopt an exponential density profile $\rho (R)={\rho }_{0}\,\exp (-R/{R}_{d})$ depending only on the galactocentric distance R and on scale length R_d. Outside the disk (i.e., for $Y\geqslant h$), the density is assumed to drop rapidly to zero. We assume that the processes occurring in the regions reached by the expanding shell (white regions in Figure 2) do not affect the expansion of the shock in the other regions interior to the disk (the red region in Figure 2), a reasonable assumption for supersonic shocks.

Within the above framework, we numerically solve Equations (5)–(8) for the expansion of the shock, and Equations (10)–(13) describing the evolution of the shock temperature T_S and density n_S with the shocked gas shell. We consider a grid of 20 equally spaced values of $0\leqslant \theta \leqslant \pi /2$ to derive at each time step the shock radius ${R}_{S,\theta }(t)$ at different inclinations θ, and update the corresponding values of the surface S(t) and volume Q(t) of the hot bubble and of the associated geometrical factor $C({R}_{S,\theta })$ entering Equations (5)–(8) and (10)–(13). Notice that until the shock radius becomes larger than the disk size h, the evolution is isotropic, due to our assumed isotropic form g(x) for the gas density distribution inside the disk, so that C = 1. When the shock breaks out of the disk, and the bubble expansion no longer retains an isotropic shape, the values of Q(t), C, and S(t) are computed numerically.

We first derive our solutions for the expansion of the shock in the plane of the galactic disk. Then, we use the properties of such solutions to understand the full two-dimensional structure of the outflows.

3.2.1. Solutions for Shock Expansion on the Plane of the Disk

We consider a galactic gas density profile $g(x)=\exp (-x/\xi )$, where $\xi \equiv {R}_{d}/{R}_{v}$ is the ratio between the disk scale length R_d and the virial radius R_v, in turn related to the circular velocity ${V}_{c}=10\,H(z){R}_{v}$ (Mo et al. 1998). We take $\xi =1/60$, a value consistent with determinations from both detailed disk models (Mo et al. 1998; assuming a DM angular momentum parameter λ = 0.05) and existing observations (Courteau 1996, 1997; for a recent review, see Sofue 2018). The form factor I entering the normalization of the density profile ${\rho }_{0}={M}_{\mathrm{gas}}/4\,\pi \,{R}_{v}^{3}\,I$ is obtained by performing the volume integral $I=4\,\pi {\int }_{0}^{1}\,g(x){x}^{2}{\int }_{0}^{{arctg}(h/x\,{R}_{v})}\,\cos \theta \,d\theta \,{dx}$ over the regions interior to the disk (see Section 2).

The results are illustrated in Figure 3 for different values of the AGN luminosity L_AGN and galactic gas mass M_gas, assuming a fixed fiducial value for DM mass $M=2\,{10}^{12}\,{M}_{\odot }$ (corresponding to a circular velocity V_c = 200 km s⁻¹).

Zoom In Zoom Out Reset image size

The qualitative behavior of the shock velocity and outflow mass is similar to that obtained for a power-law density profile, with a direct dependence of ${M}_{S}={M}_{S,\theta =0}$ on both AGN luminosity and galactic gas content, while ${V}_{S}={V}_{S,\theta =0}$ increases with increasing AGN luminosity and decreases with gas mass M_gas. We note that in this case the outflow rate ${\dot{M}}_{S}$ reaches a maximum value and declines at large radii. This is due to two factors: the first is the rapid drop in the gas density ρ, while the second is related to the large volume encompassed by the bubble expansion when it breaks out of the disk and is expressed by the quantity $C({R}_{S,\theta })$ in Equations (4) and (5). This retains the initial value $C({R}_{S,\theta })=1$ until the bubble reaches the disk boundary in the vertical direction. The subsequent rapid expansion of the bubble in the vertical direction (due to the low density encountered in this direction; see Section 3.2.2) reduces the pressure exerted on the portion of the bubble surface contained in the disk, thus reducing the expansion and mass outflow in the direction parallel to the disk. In fact, in such direction, the volume of the sphere with radius ${R}_{S,\theta }$ increasingly becomes smaller than the actual volume of the bubble Q(t), yielding progressively larger values for the quantity $C({R}_{S,\theta })$ and resulting in a smaller efficiency of the propulsive effect of the pressure term in Equations (1), (3), (4), and (5).

As for the properties of the shocked shell, the initial decline of the temperature T_S is followed by a sharp drop, due to fast radiative cooling when ${T}_{S}\approx {10}^{6.5}$ K, and T_S drops rapidly to 10⁴ K. Correspondingly, the density increases sharply; for a given total gas mass M_gas, the density reached by the shocked gas shell depends on the position of the shock at the moment of gas cooling and can easily reach values as large as ${10}^{4}\,{\mathrm{cm}}^{-3}$ at ${R}_{S}\lesssim 1\,\mathrm{kpc}$, typical of observed molecular outflows. Of course, the temperature drop (and hence the increase in density) of the shocked gas shell is delayed for increasing values of the AGN luminosity (due to the heating term in Equation (9)), while it is favored by increasing values of the total gas mass M_gas, which also correspond to increasing shocked gas densities n_S.

3.2.2. The Two-dimensional Structure of the Outflows

We now proceed to compute the full two-dimensional structure of the outflows. In this case, we compute the expansion of the shock in all directions by solving Equations (4)–(7), assuming the usual exponential density profile $\rho (R)={\rho }_{0}\,\exp (-R/{R}_{d})$ until the shock position reaches the disk boundary and a vanishing density in the regions external to the disk, i.e., where the vertical distance Y from the plane of the disk is larger than the scale height h. Our approach is similar to that adopted by Hartwig et al. (2018), although the latter authors adopt a scale height h that depends on the position along the disk.

The solutions for the outflow velocity ${V}_{S,\theta }$, mass outflow rate ${\dot{M}}_{S,\theta }$, and shock position ${R}_{S,\theta }$ as a function of time are plotted in Figure 4 for a reference galaxy with DM mass $M={10}^{12}\,{M}_{\odot }$, gas mass ${M}_{\mathrm{gas}}={10}^{10}\,{M}_{\odot }$, and AGN bolometric luminosity $L\,={10}^{45}$ erg s⁻¹. The X coordinate represents the distance from the center in the direction parallel to the plane of the disk, while the Y coordinate corresponds to the distance in the (vertical) direction perpendicular to the disk. Along the plane of the disk, the velocity ${V}_{S,\theta }$ rapidly decreases with increasing radius (top-left panel), while in the vertical direction the shock decelerates until it reaches the disk boundary h, but it rapidly accelerates afterward, due to the drop in the gas density outside the disk. The opposite is true for the mass outflow rate (top-right panel), which instead grows appreciably only along the plane of the disk, where the larger densities allow values ${\dot{M}}_{S,\theta =0}\sim {10}^{3}$ M_⊙ yr⁻¹ to be reached, as we have seen in the previous section.

Zoom In Zoom Out Reset image size

As for the shock expansion radius, this follows the paths of least resistance (see the bottom panel of Figure 4), yielding an elongated shock front in the vertical direction. For example, inspection of Figure 4 (bottom panel) shows that while in the direction perpendicular to the disk the outflows reaches a distance of 20 kpc in approximately 10⁷ yr, it takes about 10⁸ yr to reach the same distance in the plane of the disk. This has important implications for studies of AGN feedback in galaxy formation models—e.g., for an AGN lifetime of $\sim {10}^{8}$ yr, this would result in null gas expulsion along the plane of the disk.

We notice that such a behavior does not depend on the particular choice of the cutoff in the initial density distribution outside the disk (i.e., for $Y\geqslant h$). Indeed, Hartwig et al. (2018) find similar results for a radius-dependent scale length and with a different functional form for the cutoff. Thus, although the shock expansion follows the paths of least resistance (see the bottom-left panel of Figure 4), yielding an elongated shock front in the vertical direction, it is only in directions close to the plane of the disk that massive outflows (${\dot{M}}_{S}={10}^{2}\mbox{--}{10}^{3}\,{M}_{\odot }$ yr⁻¹) can be generated. This is shown in detail in the bottom-right panel of Figure 4, where the expansion of the shock position ${R}_{S,\theta }$ is shown as a function of time for both the vertical and the horizontal directions, along with the mass ${M}_{S,\theta }$ swept out by the outflows in the considered directions. While in the direction perpendicular to the disk the shock expands rapidly to reach 10 kpc in a short timescale, $\approx 2\times {10}^{7}$ yr, the denser medium encountered by the shock in the direction parallel to the disk results in a slower expansion (a distance of 10 kpc is reached only after $t\approx {10}^{8}$ yr). However, the mass swept out by the outflow in the vertical direction saturates to a small value ${M}_{s\perp }\approx {10}^{7}\,{M}_{\odot }$, while a much larger value ${M}_{s\parallel }\approx {10}^{9}\,{M}_{\odot }$ is attained along the direction parallel to the disk.

While outflows in the vertical directions can have a significant impact on the expulsion of gas from the disk (due to the large velocities attained on a short timescale), on the escape fraction of UV photons from the galactic center, and possibly on the formation of Fermi bubbles like those detected in our Galaxy (see Su et al. 2010), they are of minor importance in determining the observed massive molecular and ionized outflows in galaxies.

The above numerical solutions allow us to perform a detailed comparison with available data concerning observed molecular and ionized outflows. When key properties of the host galaxy mass are measured, we can use the observed AGN luminosity L_AGN, the total gas mass M_gas, and the dynamical mass M (or equivalently, the circular velocity V_c) as inputs for the model. This allows us to compute, for each observed galaxy, the expected expansion properties of the shock and to compare the results with the observed properties of the outflow in the considered galaxies. To perform a fair comparison, we must take into account how the properties of the outflows are derived from observations. First, present observations are not able to resolve the angular dependence of the outflow quantities (i.e., R_S, V_S, and ${\dot{M}}_{S}$). Thus, when comparing with observations, we first compute the full two-dimensional solutions in all directions (i.e., ${R}_{S,\theta }$, ${V}_{S,\theta },$ and ${\dot{M}}_{S,\theta }$), and then we derive their mass-weighted average over the directions θ, which are then compared with the observed values. In the following, for the sake of simplicity, we simply denote such averaged quantities by R_S, V_S, and ${\dot{M}}_{S}$. A second consideration concerns the measurement of the mass outflow rate, which is observationally derived as ${\dot{M}}_{S}={M}_{s}\,{V}_{S}/{R}_{S}$ (see the appendix of Fiore et al. 2017 for a discussion). Such a definition is not identical to the (θ-averaged) mass outflow rate ${\dot{M}}_{S}$ derived from the time derivative of Equation (8). When comparing with data, we will present the model predictions for both definitions, and we show that, in the regions usually covered by observations, they are basically equivalent.

4.1. Comparison with Single Objects

For molecular outflows, and for a single ionized outflow, observations in the literature have recently led to the assembly of a sizable sample of objects for which the gas mass and the rotation velocity of the host galaxies have been measured. In the following, we compare with the data sample summarized in Table 1, which extends the sample of molecular outflows in Fiore et al. (2017) to include those in M51, Circinus, XID2028, zC400528, APM08279, and 3C 298 (references are given in the caption). In addition, the data for I11119 have been updated using the recent results by Veilleux et al. (2017). We use only AGNs for which there is not only a measurement of the physical properties of the outflow (the physical size R_S, the velocity V_S, and the mass outflow rate ${\dot{M}}_{S}$) but also an estimate of the gas mass within R_S, the projected rotation velocity ${V}_{\mathrm{rot}}(\lt {R}_{S})\sin \,i$ within R_S, and the inclination angle i; for some objects, the asymptotic rotation velocity V_asympt is also available.

Table 1.  Sample of Observed Outflows

Object	Redshift	LAGN	MBH	log ${M}_{\mathrm{gas}}(\lt {R}_{S})$	${V}_{\mathrm{rot}}(\lt {R}_{S})\sin (i)$	i	Vasympt	RS	VS	${\dot{M}}_{S}$	References
		(1045 erg s−1)	$({10}^{9}\,{M}_{\odot }$)	(M⊙)	(km s−1)	(deg)	(km s−1)	(kpc)	(km s−1)	(M⊙ yr−1)
mark231	0.042	5	0.087	8.9	77	36	340	0.3	750	1000	1, 2, 3, 23
mark231	0.042	5	0.087	9.3	77	36	340	1	850	700	1,2, 4, 5, 23
n6240	0.025	0.63	0.1	9.3	230	70	⋯	0.6	500	500	6, 7 8, 9
n6240	0.025	0.63	0.1	9.8	188	70	⋯	5	400	120	6, 7, 8, 9
I08572	0.06	4.6	⋯	9.1	100	75	⋯	1	1200	1200	9, 10
I10565	0.04	0.65	0.02	9.3	75	20	250	1.1	600	300	9, 10, 11, 12
I23060	0.17	11.5	⋯	10.4	175	75	⋯	4	1100	1100	10
I23365	0.06	0.47	0.037	9.47	130	30	260	1.2	600	170	9, 10, 11, 12
J1356	0.12	1.25	0.3	8.5	200	45	⋯	0.3	500	350	13
NGC 1068	0.03	0.087	0.01	7.8	52	41	270	0.1	200	120	5, 14, 15
IC 5063	0.01	0.1	0.055	7.7	166	74	⋯	0.5	400	22	16, 17, 18, 19
NGC 1266	0.01	0.02	0.003	8.6	110	34	⋯	0.45	360	13	20, 21
I17208	0.04	1.3	0.05	11.13	130	30	⋯	1	370	65	12, 22, 23
I11119	0.19	15	0.016	9.95	142	30	⋯	7	1000	800	24
M51	0.002	0.1	0.001	9.8	⋯	22	200	0.04	100–200	11.6	25, 26
Circinus	0.001	0.04	0.0017	8.46	⋯	65	220	0.45	150	3.1	27
XID2028	1.6	20	⋯	10	210	30	350	10	700	350	28, 29
zC400528	2.3	1.7	⋯	11	250	37	⋯	4.2	450	768	30, 31
APM08279	3.9	280	10	11.15	550	30	⋯	0.27	1340	1000	32, 33
3C 298	1.43	70	3.2	9.81	190	54	⋯	1.6	400	2300	34

References. 1 = Feruglio et al. (2015); 2 = Lonsdale et al. (2003); 3 = Davies et al. (2004), 4 = Veilleux et al. (2009), 5 = Davies et al. (2007); 6 = Feruglio et al. (2013); 7 = Tacconi et al. (1999); 8 = Engel et al. (2010); 9 = Howell et al. (2010); 10 = Cicone et al. (2014); 11 = Dasyra et al. (2006); 12 = Downes & Solomon (1998); 13 = Sun et al. (2014); 14 = Garcia-Burillo et al. (2014); 15 = Krips et al. (2011); 16 = Morganti et al. (1998); 17 = Morganti et al. (2013); 18 = Woo & Urry (2002); 19 = Malizia et al. (2007) 20 = Alatalo et al. (2011); 21 = Alatalo et al. (2014); 22 = Veilleux et al. (2013); 23 = Xia et al. (2012); 24 = Veilleux et al. (2017); 25 = Querejeta et al. (2016); 26 = Shetty et al. (2007); 27 = Zschaechner et al. (2016); 28 = Perna et al. (2015a); 29 = Brusa et al. (2015b); 30 = Genzel et al. (2014); 31 = Herrera-Camus et al. (2019); 32 = Feruglio et al. (2017); 33 = Riechers et al. (2009); 34 = Vayner et al. (2017).

Download table as:  ASCIITypeset image

For each object, the observed AGN and host galaxy properties summarized in the left side of Table 1 (left of the vertical line) are used to obtain the input quantities for the model. The most uncertain quantity is the host circular velocity V_c. For most objects, we derive a lower limit from the rotation velocity within R_S (corrected for the inclination) and explore the effect of changing the assumed value of V_c. For objects where the asymptotic rotation velocity has robust estimates, we adopt such a value as an upper limit (for the two objects where only the asymptotic is available, we explore the effect of assuming lower values). The input value for the total gas mass M_gas is derived by extrapolating the observed value ${M}_{\mathrm{gas}}(R\lt {R}_{S})$ out to the virial radius using an exponential profile with disk scale radius R_d related to the circular velocity as explained in Section 3.2. For each object, the input quantities L_AGN, V_c, and M_gas derived as above allow us to compute the corresponding predicted values of the outflow velocity and outflow rate. These are compared with the observed values of V_S and M_S (shown to the right of the vertical line in Table 1).

The comparison between model predictions and observations is shown in Figure 5 for each object in Table 1. The values of V_c, M_gas, and L_AGN that have been adopted as inputs for the model are shown in the labels for each objects. The model seems to capture the basic dependence of the outflow velocity and mass outflow rate on the AGN luminosity L_AGN and host galaxy gas mass M_gas for a relatively large range of input parameters ${10}^{44}\lesssim {L}_{\mathrm{AGN}}/\mathrm{erg}\,{{\rm{s}}}^{-1}\leqslant {10}^{46}$ and ${10}^{8}\lesssim {M}_{\mathrm{gas}}/{M}_{\odot }\leqslant {10}^{10}$ covered by the data in Table 1. For all objects, the density n_S of the shocked gas at the observed outflow positions is close to the critical threshold for emission from the rotational transition of CO (corresponding to 2700 cm⁻³), although in a few cases the predicted densities are slightly below the critical threshold (for I11119, I10565, NGC 1266, J356). However, as noted in the introduction, a detailed treatment of the position-dependent ionization properties of the shocked gas and its molecular content requires numerical simulations (Richings & Faucher-Giguère 2018a, 2018b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The evolution of the outflow velocity in the models is characterized by an upturn. This is related to the two-dimensional properties of the shocks discussed above. Although the mass-weighted average V_S is initially dominated by the component θ = 0 aligned with the plane of the disk (due to the large mass involved in such a direction; see Figure 4) when the shock in the plane of the disk reaches a standstill, the average V_S is mainly contributed by components not aligned with the plane of the disk, characterized by an increasing expansion velocity (see the top-left panel of Figure 4) related to the low gas density encountered in such a direction. The same effect is responsible for the downturn of the average mass outflow rate ${\dot{M}}_{s}$. In fact, this is largely contributed by the gas mass in the disk; however, in the disk direction, the combined effect of large densities and of the drop in the pressure term associated with the bubble expansion as it breaks out of the disk (see Section 3.2.1) leads to the drop in the mass outflow rate.

We note that our results are not sensitive to the uncertainties affecting V_c (and hence the extrapolated M_gas) defining the properties of the host galaxy. This results from the balance between the effect of changing V_c and the correction that relates the observed values ${M}_{\mathrm{gas}}(\lt {R}_{S})$ in Table 1 to the overall gas content M_gas. For example, increasing V_c tends to shift the predicted curves on the right along the x-axis; however, such effect is balanced by the larger value of M_gas corresponding to the observed ${M}_{\mathrm{gas}}(\lt {R}_{S})$.

In the last panel of Figure 5, we present the case of SDSS J1148 at z = 6.4, where Maiolino et al. (2012) and Cicone et al. (2015) reported the detection of a massive [C ii] outflow (i.e., associated with the cold gas phase of the ISM) powered by a high-luminosity QSO for which the properties of the observed host galaxy (i.e., M_gas, V_rot) are known. Even for this object, characterized by a huge AGN luminosity, a large R_S, and an extreme mass outflow rate, the agreement with the model predictions is excellent. In this case, we also obtain a predicted density n_S of the shocked gas shell that is lower than the critical density for the molecular CO emission, as expected in the case of ionized outflows.

4.2. Comparison with High-luminosity, Ionized Outflows

For all other ionized outflows present in the literature, detailed measurements of the host galaxy M_gas and V_rot are not available. In Table 2, we report the values of V_S and ${\dot{M}}_{S}$ for a large sample of objects taken from Fiore et al. (2017), where objects are listed in order of increasing bolometric luminosity L_AGN. From the sample in Fiore et al. (2017), we have excluded I10565, due to the ambiguous interpretation of the outflow (possible earlier bubbles due to previous ejection episodes; see Rupke & Veilleux 2013), Mrk 231 (the analysis of Rupke & Veilleux 2013 excludes a part of the nuclear emission), and J1339 (its identification with an AGN is uncertain; see Harrison et al. 2014).

Table 2.  Sample of Observed Ionized Outflows: Objects Sorted by Increasing AGN Luminosity

Object	Redshift	LAGN (1045 erg s−1)	RS (kpc)	VS (km s−1)	${\dot{M}}_{S}$ (M⊙ yr−1)	References
SDSS J0958	0.10	45.0	2.6	866	1.1	1
SDSS J1356	0.12	45.1	3.1	1049	1.6	1
SDSS J1130	0.13	45.1	2.8	616	0.3	1
SDSS J1125	0.17	45.2	2.9	1547	0.75	1
SDSS J1430	0.08	45.3	1.8	999	1.7	1
SDSS J1316	0.15	45.4	3.1	1216	1.48	1
SDSS J0945	0.13	45.5	2.7	1511	1.62	1
SDSS J10100	0.10	45.6	1.6	1267	1.46	1
GS3-19791	2.22	45.6	1.3	530	3.23	2
SDSS J1000	0.15	45.7	4.3	761	1.16	1
SDSS J1355	0.15	45.7	3.5	797	0.57	1
GS3-28008	2.29	45.9	1.3	300	2.34	2
XID 5395	1.47	45.9	4.3	1600	2.65	4
SDSS J10101	0.20	46	3.9	1523	1.82	1
SDSS J1100	0.10	46	1.9	1192	1.65	1
SDSS J0210	0.54	46.1	7.5	560	2.62	5, 6, 7, 8
SDSS J1040	0.49	46.2	7.6	1821	3.16	5, 6, 7, 8
COS 11363	2.10	46.2	1.3	1240	2.83	2
SMM J1636	2.38	46.3	7	1054	1.44	9
MRC 0406	2.44	46.3	9.3	960	3.82	10
XID 5321	1.47	46.3	11	1950	1.84	11, 12
RG J0302	2.24	46.3	8	1234	1.48	9
SDSS J0319	0.62	46.4	7.5	934	2.32	5, 6, 7, 8
SDSS J0321	0.64	46.5	11	946	2.30	5, 6, 7, 8
SDSS J0841	0.64	46.5	6.4	675	2.60	5, 6, 7, 8
MIRO 20581	2.45	46.6	4.8	1900	2.29	13
MRC 1138	2.20	46.6	20	800	2.39	14
MRC 0828	2.57	46.6	9	800	3.87	10
SMM J1237	2.06	46.7	7	1200	1.48	9
SMM J0943	3.35	46.7	15	1124	1.57	9
SDSS J0842	0.56	46.8	9.	522	2.59	5, 6, 7, 8
HB 8905	2.48	46.8	1.3	500	2.65	9
SDSS J1039	0.58	46.9	5.8	1046	2.81	5, 6, 7, 8
SDSS J0149	0.57	46.9	4.1	1191	2.60	5, 6, 7, 8
SDSS J0858	0.45	47.2	5.6	939	2.79	5, 6, 7, 8
HB 8903	2.44	47.3	1.9	1450	1.76	13
SDSS J0759	0.65	47.3	7.5	1250	2.87	5, 6, 7, 8
HE 0109	2.40	47.4	0.4	900	3.14	5, 6, 7, 8
LBQS0109	2.35	47.4	0.4	1850	2.84	13
SDSS J1326	3.30	47.6	7	2160	3.81	14, 15
SDSS J1201	3.51	47.7	7	1850	3.50	14, 15
SDSS J1549	3.30	47.8	7	1380	3.42	14, 15
SDSS J0900	3.30	47.9	7	2380	3.52	14, 15
SDSS J0745	3.22	48.0	7	1890	3.76	14, 15

References. 1 = Harrison et al. (2014); 2 = Genzel et al. (2014), assuming ${H}_{\alpha }/{H}_{\beta }$ = 2.9, extinction corrected; 3 = Cicone et al. (2014), 4 = Brusa et al. (2016); 5 = Liu et al. (2013a); 6 = Liu et al. (2013b), extinction corrected; 7 = Wylezalek et al. (2016); 8 = Reyes et al. (2008); 9 = Harrison et al. (2012); 10 = Nesvadba et al. (2008); 11 = Brusa et al. (2015a); 12 = Perna et al. (2015a); 13 = Perna et al. (2015b); 14 = Nesvadba et al. (2006); 13 = Carniani et al. (2015); 14 = Bischetti et al. (2017); 15 = Duras et al. (2017).

Download table as:  ASCIITypeset image

Because the properties of the observed host galaxies are not available, we cannot perform a detailed on-by-one comparison with the model as we did for molecular outflows. Thus, we have divided the observed AGN luminosity range into different bins. For each luminosity bin, the input quantities for the model, i.e., the total mass M (or the circular velocity V_c) and the gas mass M_gas of the host galaxy, are related to L_AGN through available average scaling relations. We adopt the total mass, M, from the central value of the relation $\mathrm{log}({M}_{\mathrm{BH}}/{M}_{\odot })\,=(1.55\pm 0.02)\mathrm{log}(M/{M}_{\odot })-(11.26\pm 0.20)$ found from the analysis of the Illustris N-body simulations by Mutlu-Pakdil et al. (2018); such a relation is consistent with a wide set of observational data (see references in the above paper). An average $M\mbox{--}{L}_{\mathrm{AGN}}$ relation is then derived after converting the black hole mass to bolometric luminosity assuming Eddington emission; assuming an Eddington ratio peaked at 0.3 (as indicated by some observations; see, e.g., Kauffmann & Heckman 2009; Shankar et al. 2013) does not appreciably change our results. Although, observationally, the large scatter of the relation at small galaxy masses makes the correlation weak (see, e.g., Kormendy & Bender 2011; Sabra et al. 2015), the scatter reduces appreciably for large halos with ${V}_{c}\gtrsim 200$ km s⁻¹ and large black hole masses ${M}_{\mathrm{BH}}\gtrsim {10}^{8}\,{M}_{\odot }$, such as those corresponding to the objects in Table 2. To derive the other input quantity M_gas, we use the relation $\mathrm{log}({M}_{\mathrm{gas}}/M)\,\approx -2-[[\mathrm{log}(M/{M}_{\odot })-12]$ approximating the scaling found from the abundance matching by Popping et al. (2015) for galaxies with $12\leqslant \mathrm{log}M/{M}_{\odot }\leqslant 13$, the range covered by the masses corresponding to the luminosities in Table 2. With such an approximation, the gas mass stays constant at ${M}_{\mathrm{gas}}\approx {10}^{10}\,{M}_{\odot }$ over the entire interval of interest,$M={10}^{11}\mbox{--}{10}^{13}\,{M}_{\odot }$ (corresponding to ${V}_{c}=150\mbox{--}400$ km s⁻¹).

With the above approximations for the input values of M and M_gas, we computed the model predictions for the outflow velocity V_S, outflow mass rate ${\dot{M}}_{S}$, and shocked gas density n_S for different bins of L_AGN, and compare them with the observed values taken from Table 2 for each L_AGN bin. The results are shown in Figure 6 for AGN luminosities ranging from 10⁴⁵ to 10⁴⁸ erg s⁻¹. To account for uncertainties in the input values of M_gas derived from the scaling law in Popping et al. (2015), in each bin we show the effect of assuming an input value of M_gas differing from the average relation by a factor of 3 above and below the mean value. In all panels, the model predictions are computed at z = 0. Because we cannot perform a one-by-one comparison with data points, the model predictions are all computed at z = 0. However, objects at high redshifts are more compact (see, e.g., Mo et al. 1998), and all sizes are expected to scale accordingly. Thus, to compare with data corresponding to objects with different redshifts in the same plot, we have rescaled all observed sizes to z = 0 according to the expected evolution of the disk radius ${r}_{d}{(z)/{r}_{d}{(z=0)=({{\rm{\Omega }}}_{{\rm{\Lambda }}}+{{\rm{\Omega }}}_{0}(1+z)}^{3})}^{-0.5}$ (Mo et al. 1998; with density parameters ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.7$ and ${{\rm{\Omega }}}_{0}=0.3$ for the dark energy and matter, respectively).

Zoom In Zoom Out Reset image size

Within the unavoidable uncertainties due to the derivation of the input quantities from the above average relations, the model predictions are in general in agreement with the observations, the agreement becoming excellent for the highest luminosity bins. Also, for all objects, the predicted shocked gas density is below the value required for the CO emission, as appropriate for ionized outflows. Notice that in the vast majority of cases (although not in all of them), the shocked gas has reached the cooling radius at the observed shock position, so the expected temperatures for the shocked shell are ${T}_{S}\sim {10}^{4}$ K. Nevertheless, the large AGN luminosities push the cooling radius to the outer regions, where the lower gas densities yield values for n_S smaller than the threshold for CO emission. The huge range spanned by the model predictions when L_AGN is changed is also noticeable, with mass outflow rates ranging from ${\dot{M}}_{S}\,\sim 10\,{M}_{\odot }$ yr⁻¹ for the lowest luminosity bin to ${\dot{M}}_{S}\sim {10}^{3}\,{M}_{\odot }$ yr⁻¹ for ${L}_{\mathrm{AGN}}\sim {10}^{48}$ erg s⁻¹. The agreement of model predictions with observations over such a large range of input and output quantities provides a strong support to the model predictions. On the other hand, in the lowest luminosity bin, $45\leqslant {{\rm{l}}{\rm{o}}{\rm{g}}{L}}_{{\rm{A}}{\rm{G}}{\rm{N}}}/{\rm{e}}{\rm{r}}{\rm{g}}\,{{\rm{s}}}^{-1}\leqslant 45.4$, the model overpredicts the mass outflow rates. However, as discussed in Fiore et al. (2017), the measured ionized mass outflow in low-luminosity objects is likely to represent only a fraction of the total mass outflow. In particular, from the comparison with model predictions, we expect a correction factor of ∼10–50 in this luminosity range. Thus, observational determinations of the fraction of mass outflow in ionized winds in low-luminosity objects will constitute an important consistency check for the model predictions in this regime.

Finally, we focus on the scaling properties of the outflow quantities V_S and ${\dot{M}}_{S}$. In Figure 7, we show the dependence of both quantities on the AGN luminosity for observed outflows at small (${R}_{S}\leqslant 1$ kpc) and large (${R}_{S}\geqslant 1$ kpc) distances from the galaxy center. The observational data points are compared with the model predictions for V_S and ${\dot{M}}_{S}$ at different luminosities computed at ${R}_{S}=0.5\,\mathrm{kpc}$ (top panels) and ${R}_{S}=7\,\mathrm{kpc}$ (bottom panels). In all cases, the assumed value for V_c has been computed following the procedure described in Section 5, while we considered three equally spaced values for M_gas in the range $(0.3\mbox{--}3)\times {10}^{10}\,{M}_{\odot }$ as done in Section 5 and in Figure 6. Thus, we do not perform a one-by-one comparison between the data and the model predictions as the latter are computed only at particular values of R_S and M_gas.

Zoom In Zoom Out Reset image size

The predicted V_S scales as ${V}_{S}\sim {L}_{\mathrm{AGN}}^{0.35}$ at small radii and as ${V}_{S}\sim {L}_{\mathrm{AGN}}^{0.37}$ at larger radii ${R}_{S}\geqslant 1\,\mathrm{kpc}$. The correlation is consistent with the best fit to the data, ${V}_{S}\sim {L}_{\mathrm{AGN}}^{0.3\pm 0.4}$, given in Fiore et al. (2017), although we stress that the model slope is computed at fixed M_gas and R_S, while the observed points in the ${V}_{S}\mbox{--}{L}_{\mathrm{AGN}}$ plane are characterized by different values of gas mass and shock position.

The corresponding predicted scaling of the mass outflow rate is ${\dot{M}}_{S}\sim {L}_{\mathrm{AGN}}^{0.3}$ for ${R}_{S}=0.5\,\mathrm{kpc}$, and ${\dot{M}}_{S}\sim {L}_{\mathrm{AGN}}^{0.35}$ for ${R}_{S}=7\,\mathrm{kpc}$. In this case, the observed overall correlations show a sensibly stronger dependence (${\dot{M}}_{S}\sim {L}_{\mathrm{AGN}}^{0.76\pm 0.06}$ and ${\dot{M}}_{S}\sim {L}_{\mathrm{AGN}}^{1.29\pm 0.38}$ for molecular and ionized outflows, respectively), although with some variance depending on the observational sample (see, e.g., Bischetti et al. 2019) However, besides the scatter in the R_S and M_gas of the observed points, the observed steep correlation for ionized winds is largely determined by the points with small ${\dot{M}}_{S}$ at large radii ${R}_{S}\geqslant 1\,\mathrm{kpc}$. Indeed, most of the ionized outflows in the lower-right panel of Figure 7 are below the value expected by our model at low ${L}_{\mathrm{bol}}\lesssim {10}^{46}$ erg s⁻¹. This is due to the fact that at low L_bol, ionized outflows represent only a small fraction of the total mass outflow rate. At high L_bol, the mass outflow rates of ionized and molecular outflows are about similar, meaning that both can be used as relatively good tracers of the total mass outflow rate (see Fiore et al. 2017).

In any case, the above comparison is largely affected by the fact that the model predictions are computed at particular values of R_S and M_gas, while the observed points correspond to objects with a wide range of R_S and M_gas. Thus, the steeper logarithmic slope of the observed correlation can be due to either a true inadequacy of the model in describing ionized outflows in low-luminosity objects, or large biases affecting the determination of the total mass outflow rate from the observation of ionized outflows in low-luminosity AGNs, or to the intrinsic scatter in R_S and M_gas of the observational data points.

Although present data are too sparse and incomplete to allow for detailed determination of the ${\dot{M}}_{S}\mbox{--}{L}_{\mathrm{AGN}}$ and ${V}_{S}\mbox{--}{L}_{\mathrm{AGN}}$ relations in different bins of M_gas and R_S, a step forward to test the predicted correlations can be performed by scaling the observed values of V_S and ${\dot{M}}_{S}$ (corresponding to different outflows with different AGN luminosity and gas mass and circular velocity) to a reference value of V_c, M_gas, and L_AGN using the predictions of the model. The run of the rescaled observed quantities with R_S can then be compared with the model predictions for ${V}_{S}({R}_{S})$ and ${\dot{M}}_{S}({R}_{S})$ computed at the reference values for V_c, M_gas, and L_AGN.

To this aim, we first express the model predictions in a compact form. In fact, the scaling properties of our solutions with the main input quantities of the model at different shock radii R_S can be approximated by the following power- law relations:

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}\left(\displaystyle \frac{{V}_{S}}{\mathrm{km}\,{{\rm{s}}}^{-1}}\right) & = & 2.65+0.35\,\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{AGN}}}{{10}^{45}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)\\ & & -\,\displaystyle \frac{1}{3}\mathrm{log}\left(\displaystyle \frac{{M}_{\mathrm{gas}}}{{10}^{10}\,{M}_{\odot }}\right)+\mathrm{log}\left(\displaystyle \frac{{V}_{c}}{200\,\mathrm{km}\ {{\rm{s}}}^{-1}}\right)\\ & & -\,0.65\,\mathrm{log}\left(\displaystyle \frac{{R}_{S}}{1\mathrm{kpc}}\right),\end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}\left(\dot{\displaystyle \frac{{M}_{S}}{{M}_{\odot }{\mathrm{yr}}^{-1}}}\right)=1.9+0.35\,\mathrm{log}\left(\displaystyle \frac{{L}_{\mathrm{AGN}}}{{10}^{45}\,\mathrm{erg}\ {{\rm{s}}}^{-1}}\right)\\ \ \ +\,0.65\,\mathrm{log}\left(\displaystyle \frac{{M}_{\mathrm{gas}}}{{10}^{10}\,{M}_{\odot }}\right)-1.7\,\mathrm{log}\left(\displaystyle \frac{{V}_{c}}{200\mathrm{km}\ {{\rm{s}}}^{-1}}\right)\\ \qquad +\,1.3\,\mathrm{log}\left(\displaystyle \frac{{R}_{S}}{1\,\mathrm{kpc}}\right).\end{array}\end{eqnarray} \tag{ 15 }$$

We stress that the above equations constitute a valid approximation (within 10%) only in the inner region ($\lt 5\,\mathrm{kpc}$ for most cases) where the scaling of the solutions with R_S is approximately a power law. At larger radii, our solutions (as described in Section 3) cannot be fitted with a single power law, and are characterized by a turnover that depends on the input quantities V_c, M_gas, and L_AGN.

The above fitting formulas allow us to test the typical dependencies of our model on the input quantities V_c, M_gas, and L_AGN against observations. We first rescale the data (corresponding to different outflows with different AGN luminosity and gas mass and circular velocity) to a reference value of V_c, M_gas, and L_AGN through the dependencies in Equations (14) and (15). This give evidence of the dependence of the data points on the shock radius R_S, which can be compared with predictions. Such a test is performed in Figure 8, where all points (corresponding to objects in Table 1 for which measurements of M_gas, R_S, and V_S are available) have been rescaled to the same reference values ${L}_{\mathrm{AGN}}={10}^{45}$ erg s⁻¹, ${M}_{\mathrm{gas}}={10}^{10}\,{M}_{\odot }$, and V_c =200 km s⁻¹ after Equations (14) and (15), and the resulting dependence on R_S is compared with the model solutions for the reference input value of V_c, M_gas, and L_AGN.

Zoom In Zoom Out Reset image size

The correspondence of model predictions with the observed scalings of ${\dot{M}}_{S}$ is actually very good, despite the simple assumptions of the model and the large uncertainties affecting the data. Only for one object (NGC 1068, green points in the figure) does the model yield significant deviations from the observed values of V_S and ${\dot{M}}_{S}$ (see also Figure 5). However, this galaxy is characterized by radio jets not aligned with the line of sight (Crane & van der Hulst 1992; Gallimore et al. 1996), a property shared by Circinus and M51. Thus, in this case, the observed molecular outflow could be determined (or affected) by the interaction of the jet with the ISM (not described by our model).

While the comparison of model predictions with present data is encouraging, we stress that at present several uncertainties affect the data with which we are comparing. Measured mass outflow rates depend on the adopted conversion from CO luminosities into H₂ masses, and on the estimated size of the outflow. In particular, the conversion factor ${\alpha }_{{CO}}=0.8\,{M}_{\odot }$ (K km s⁻¹ pc⁻²)⁻¹ adopted in Fiore et al. (2017) can be a function of density, metallicity, and gas distribution (see the discussion in Fiore et al. 2017; see Bolatto et al. 2013 for a review), while the size of the outflow is based on the maximum radius up to which high-velocity gas is detected (baseline method, but alternative methods have been proposed in the literature; see Carniani et al. 2015). Both of these uncertainties will likely greatly be reduced by future higher resolution observations with the Atacama Large Millimeter Array (ALMA) and Northern Extended Millimeter Array (NOEMA).

We computed the two-dimensional expansion of outflows driven by AGNs in galactic disks as a function of the global properties of the host galaxy and of the luminosity of the central AGN. We derived the expansion rate, the mass outflow rate, and the density and temperature of the shocked shell in the case of an exponential profile for the disk gas, for different expansion directions θ with respect to the plane of the disk. Having expressed our model results in terms of global properties of the host galaxies, we compared our predictions to a large sample of 19 outflows (mostly molecular, except for one object) in galaxies with measured AGN luminosity and gas mass, and with estimated total mass. This allowed us to perform a detailed, one-by-one comparison with the model predictions, to assess to what extent the present status of the modeling is consistent with the existing observational distribution of outflow properties.

We find—in the vast majority of cases—an encouraging agreement for a wide range of gas masses and AGN bolometric luminosities (including the hyperluminous quasars with ${L}_{\mathrm{AGN}}\approx {10}^{47}$ erg s⁻¹; see Vietri et al. 2018). The model yields—for each considered galaxy and at the observed outflow radii—values of velocity and mass outflow rates (averaged over the directions θ) that are in good agreement with observations. The predicted densities of the shocked shell are consistent with the observed molecular emission of the outflows in the vast majority of cases. Significant deviations from the model predictions are found only for NGC 1068 (concerning both the shock velocity V_S and the mass outflow rate ${\dot{M}}_{S}$) and for IC 5063 and Circinus (concerning the gas density needed to produce the observed molecular outflow). However, these three objects are all characterized by weak radio jets not aligned with the line of sight (Crane & van der Hulst 1992; Gallimore et al. 1996; Elmouttie et al. 1998; Morganti et al. 2015), which could be the origin of (or contribute to) the observed outflow, a situation beyond the reach of our model.

We notice that some features characterizing the predicted expansion and mass outflow rates are specific to the exponential density profile and two-dimensional geometry that we consider. This includes the upturn of the expansion rate and the drop of the mass outflow rate at large radii $\approx 0.1\,{R}_{v}$ (see Figures 3 and 8). Interestingly, signatures of such a typical behavior at the expected distance from the galaxy center seem to be already present in the considered sample of measured outflows (Figure 8).

We then considered a larger sample of 48 (mostly) ionized outflows in galaxies with no reliable measurements of the gas and dynamical mass, and we perform an approximate comparison of the model predictions for different bins of AGN luminosity assuming different reference values for the gas mass and dynamical mass derived from average scaling relations. Within the unavoidable uncertainties due to the derivation of the input quantities from average relations, the model predictions are in general in agreement with the observed outflow properties, the agreement becoming excellent for the highest luminosity bins. Also, for all objects, the predicted shocked gas density is below the value required for the CO emission, as appropriate for ionized outflows. Notice that, in the lowest luminosity bin, the model overpredicts the mass outflow rates measured in ionized winds. However, as discussed in Fiore et al. (2017), the measured ionized mass outflows in low-luminosity objects are likely to represent only a fraction of the total mass outflow, while our model predictions concern the total mass outflow rate. When comparing with observations, we need to assess whether the observed quantities are good tracers of the total outflow rate. The much lower ionized outflow rates with respect to molecular rates found in the past, in particular at low bolometric luminosity, suggest that only the latter are good tracers of the total outflow rate, with the former probably a good tracer only at high L_bol (Carniani et al. 2015; Fiore et al. 2017). In particular, comparing our model predictions for the total mass outflow rates with the observed rates in ionized winds for low-luminosity AGNs, $45\,\leqslant \mathrm{log}{L}_{\mathrm{AGN}}/\mathrm{erg}\,{{\rm{s}}}^{-1}\leqslant 45.4$, we expect that ionized winds trace only a fraction of the correction factor, ∼0.1 of the total mass outflow rate. Observational determinations of the fraction of mass outflow in ionized winds in low-luminosity objects will constitute an important consistency check for the model predictions in this regime.

While the encouraging, quantitative agreement of the model predictions with a wide set of existing observations constitutes a baseline for the interpretation of forthcoming data and for a more detailed treatment of AGN feedback in galaxy formation models, we stress that the comparison with observations is still affected by large uncertainties related to the data. These mainly affect the estimates of the mass outflow rate ${\dot{M}}_{S}$ and the shock position R_S. In fact, the adoption of different approaches in the measurement of the outflow velocity V_S (the velocity peak of the broad emission lines versus the width of the emission lines at 80% of the line flux)—while resulting in a mild error of ≈4% on V_S itself—produces uncertainties of ∼35% in the associated estimates of ${\dot{M}}_{S}$. For molecular outflows, the latter is also affected by the uncertainties related to the conversion of CO luminosities into H₂ masses. While Fiore et al. (2017) adopted a constant conversion factor ${\alpha }_{\mathrm{CO}}=0.8$, this can actually be a function of density, metallicity, and gas distribution (see Bolatto et al. 2013 for a review). For ionized winds, estimates of the mass of outflowing gas depend linearly on the assumed gas temperature T and inversely on the density n. While the data we base on are derived for $T={10}^{4}$ K and n = 200 cm⁻³, uncertainties up to a factor 2 can affect both quantities (see Fiore et al. 2017 and references therein). As for the size of the outflow R_S, in most cases this is taken as the maximum radius up to which high-velocity gas is detected (baseline method). On the other hand, Carniani et al. (2015) evaluate a size of the ionized wind systematically lower than all other cases, because they adopt a different astrometric procedure. However, this uncertainty will likely be reduced by future higher resolution ALMA and NOEMA observations. While the above uncertainties in present data do not allow us to unambiguously determine the effectiveness of the model in providing a full description of the expansion of outflows, the typical dependence of the mass outflow rate ${\dot{M}}_{S}$ on R_S in Figure 8 (bottom panel) could be tested in detail to probe the position of the upturn when a larger sample of outflows at large ${R}_{S}\gtrsim 5\,\mathrm{kpc}$ will be detected. Also, testing dependencies on the properties of the host galaxy (summarized in Equations (13) and(14)) on a solid statistical ground will require a large sample of outflows with associated measurements of the galaxy properties. The above observational goals actually characterize the SUPER (Survey for Unveiling the Physics and Effect of Radiative feedback) ongoing ESO's Very Large Telescope (VLT)/SINFONI Large Programme (see Circosta et al. 2018). SUPER will perform the first systematic investigation of ionized outflows in a sizable and blindly selected sample of 39 X-ray AGNs at $z\approx 2$, linking the outflow properties to a number of AGN and host galaxy properties. The large sample of outflow in AGNs with high bolometric luminosities (up to ${L}_{\mathrm{AGN}}\approx {10}^{47}$ erg s⁻¹) will enable the evolution of the outflows to be traced up to large distances from the central AGN, providing a larger leverage to test the specific predictions of the model at large values of R_S.

On the model side, our results are in excellent agreement with previous studies in overlapping cases. When we test our model solutions for a spherical initial gas distribution with a power-law decline with radius, we find results that agree with existing analytical scaling laws in the limit of energy-driven winds (see, e.g., King et al. 2011; Faucher-Giguère & Quataert 2012) and with numerical solutions in Faucher-Giguère & Quataert (2012). When an exponential disk distribution is assumed for the unperturbed galactic gas, the two-dimensional structure that we obtain is similar to that obtained by Hartwig et al. (2018) and is characterized by a shock expansion that follows the paths of least resistance (see bottom panel of Figure 4) with an elongated shock front in the direction perpendicular to the disk. In such a direction, the velocity field is characterized by a decline with increasing distance Y from the disk, followed by a strong increase for large distances $Y\gt h$. However, massive outflows (${\dot{M}}_{S}={10}^{2}\mbox{--}{10}^{3}\,{M}_{\odot }$ yr⁻¹) can be generated only in the plane of the disk.

In the present paper, we have focused on the comparison with available observations of massive outflows, which do not resolve the spatial structure of the shock. Thus, when comparing with data, we did not fully exploit the full two-dimensional description of our model, because the predicted quantities (outflow radius, velocity, and mass rate) compared with the data have been averaged over the angular direction relative to the plane of the disk. However, two-dimensional spectroscopic maps obtained by present and upcoming integral field unit (IFU) facilities will allow more detailed comparison of the mass outflow rate and velocity maps with those predicted by the present model. For example, the ongoing MAGNUM (Measuring Active Galactic Nuclei Under MUSE Microscope, Venturi et al. 2018) survey by the MUSE instrument at the VLT, aimed at mapping the ionized outflows from local AGNs, and observations with ALMA and NOEMA will provide a crucial sample to test the two-dimensional picture described by our model. In the future, the NIRSPEC IFU facility at the James Webb Space Telescope will allow us to extend any comparison to higher redshifts. As for the molecular fraction of the AGN outflows, millimeter facilities (like ALMA) will increase both sensitivity and resolution of the comparison sample.

At the same time, our two-dimensional description of outflows can be applied to a number of different investigations. For example, the model can provide a detailed estimate of the escape fraction of ionizing photons in active galaxies, an issue relevant for studies on cosmic reionization. In fact, the fast motion of the outflow in the direction perpendicular to the disk can effectively sweep out the interstellar gas. Computing such an effect as a function of the AGN power and lifetime will enable the amount of H i ionizing photons escaping from the population of AGN host galaxies to be estimated, and provide a refined computation of the contribution AGN to the ionizing background, extending and improving the results derived from the analytic blast-wave model (Menci et al. 2008) adopted in Giallongo et al. (2012). Our model can also be applied to compute the acceleration of particles in the shock front of AGN outflows and the ensuing generation of gamma-rays and neutrinos, to compute in detail the contribution of AGNs to the extragalactic gamma-ray and neutrino backgrounds (Lamastra et al. 2017). Further, the two-dimensional description developed here can contribute to provide a quantitative description of the origin of Fermi bubbles (see Su et al. 2010). As suggested by earlier authors (see Zubovas et al. 2011; Zubovas & Nayakshin 2012; Lacki 2014), these are connected to the anticorrelation between the outflow speed and the gas density, which rapidly decreases in the direction perpendicular to the disk. The consideration of the azimuthal dependence of the gas density in our two-dimensional description will then allow for a detailed, quantitative comparison between the observed properties of the Fermi-LAT lobes and the prediction of the model: we plan to address this point in a subsequent paper.

Finally, the model can provide a refined description of the AGN feedback in galaxy formation models by enabling the ejected gas mass to be computed, i.e., the gas mass that passes the virial radius with a velocity larger than the escape velocity. Indeed, most analytic and numerical cosmological models of black hole growth have included mostly "thermal-like" AGN feedback recipes (e.g., Dubois et al. 2013), which are based on approximately isotropic injections of thermal energy into the gas surrounding the black hole. Barausse et al. (2017) showed that current implementations of AGN (quasar mode) feedback in comprehensive galaxy evolution models fall drastically short in reproducing the observed strong dependence of black hole mass with velocity dispersion at fixed stellar mass. A kinetic-like feedback, such as the one discussed in this work, may provide a stronger coupling to velocity dispersion in view of the possibly more efficient removal of gas in lower mass systems. In this case, the full two-dimensional description of the outflows plays a relevant role, because the larger velocities attained in the vertical direction (perpendicular to the disk) easily exceed the escape velocity at the virial radius in a short timescale, while the slower expansion of the shock in the plane of the disk can prevent the escape of gas in this direction within the lifetime of the AGN. Inspection of Figure 4 (bottom panel) shows that while in the direction perpendicular to the disk the outflows reach a distance of 20 kpc in approximately 10⁷ yr, it takes about 10⁸ yr to reach the same distance in the plane of the disk. For the AGN with lifetime ∼10⁸ yr, this would result into null gas expulsion along the plane of the disk. The description of the shock expansion that we provide in terms of global galactic quantities allows for a fast implementation of the above description in semianalytic models of galaxy formation. We plan to investigate the above issues in subsequent papers.

We acknowledge support from INAF under PRIN SKA/CTA FORECaST and PRIN SKA-CTA-INAF ASTRI/CTA Data Challenge. F.S. acknowledges partial support from the H2020 AHEAD project (grant agreement No. 654215) and a Leverhulme Trust Research Fellowship. E.P. acknowledges financial support from ASI and INAF under contract 2017-14-H.0 ASI-INAF. We thank the referee for the constructive comments that significantly helped to improve the paper.