Discovering AGN-driven winds through their infrared emission: II. Mass outflow rate and energetics
Dalya Baron Hagai Netzer



Obscured AGN at 1.5 < z < 3.0 from the zCOSMOS-deep Survey
   The  AGN catalog with electronic data is available at the CDS via anonymous ftp to 
    cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/
    M. Mignoli1
          A. Feltre2,3,4
          A. Bongiorno5
          F. Calura1
          R. Gilli1
          C. Vignali6,1
          G. Zamorani1
          S.J. Lilly7
          O. Le Fèvre8
          S. Bardelli1
          M. Bolzonella1
          R. Bordoloi9,10
          V. Le Brun8
          K.I. Caputi11
          A. Cimatti6,12
          C. Diener13
          B. Garilli14
          A.M. Koekemoer15
          C. Maier16
          V. Mainieri17
          Y, Peng18
          E. Pérez Montero19
          J.D. Silverman20
          E. Zucca1

    Received ; accepted 
============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






The global influence of AGN-driven outflows remains uncertain, due to a lack of large samples with accurately-determined outflow properties. In the second paper of this series, we determine the mass and energetics of ionized outflows is 234 type II AGN, the largest such sample to date, by combining the infrared emission of the dust in the wind (paper I) with the emission line properties. We provide new general expressions for the properties of the outflowing gas, which depend on the ionization state of the gas. We also present a novel method to estimate the electron density in the outflow, based on optical line ratios and on the known location of the wind. The inferred electron densities, n_e∼ 10^4.5 cm^-3, are two orders of magnitude larger than typically found in most other cases of ionized outflows. We argue that the discrepancy is due to the fact that the commonly-used [SII]-based method underestimates the true density by a large factor. As a result, the inferred mass outflow rates and kinetic coupling efficiencies are Ṁ_out∼ 10^-2 (M_⊙/yr) and ϵ = Ė_kin/L_bol∼ 10^-5 respectively, 1–2 orders of magnitude lower than previous estimates. Our analysis suggests the existence of a significant amount of neutral atomic gas at the back of the outflowing ionized gas clouds, with mass that is a factor of a few larger than the observed ionized gas mass. This has significant implications for the estimated mass and energetics of such flows.





galaxies: general – galaxies: interactions – galaxies: evolution – galaxies: active – galaxies: supermassive black holes –  galaxies: star formation









§ INTRODUCTION


The discovery of correlations between the masses of super massive black holes (SMBHs) and several properties of their host galaxy (stellar velocity dispersion, bulge mass, and bulge luminosity; e.g. ) have led to a suggestion that the growth of the SMBH is linked to the stellar mass growth in its host galaxy (e.g., ). Active galactic nuclei (AGN) feedback, in the form of powerful outflows, is invoked as a way to couple the energy released by the accreting SMBH with the ISM of its host galaxy, providing a possible explanation for the observed correlations (e.g., ). A number of such models have successfully reproduced these correlations, by requiring that a significant amount of the accretion energy of the AGN will be mechanically-coupled to the ISM of the host galaxy (∼5–10% L_bol; e.g., ). 

Powerfull winds are routinely detected in the host galaxies of AGN, spanning a large range of evolutionary stages, luminosities, and outflow gas phases (e.g., ). In order to assess the effect of such flows on their host galaxy evolution, it is necessary to accurately determine their mass and energetics, and compare them to those required by feedback models (see however ). A major source of uncertainty is related to the relative contributions of AGN and supernovae to the observed outflows (see introduction in ). It is normally assumed that the source that photoionizes the stationary and outflowing gas is also the main driver of the observed winds (e.g., ). Thus, outflows that are observed in systems with emission line ratios that are consistent with AGN photoionization, are considered to be driven by the AGN. 

A large number of AGN-driven winds are observed in the warm ionized phase, traced by rest-frame optical emission lines (e.g., ). In order to properly quantify the mass and energetics of these outflows, it is necessary to determine their kinematics (e.g., outflow velocity, v_out), spatial extent (r_out), and electron density (n_e), all of which are subjected to various uncertainties and systematics (see e.g., ). For example, the electron densities in such outflows are claimed to be rather low, in the range ∼ 10^2 - 10^3 cm^-3, based on the [SII] λλ 6717,6731Å  doublet ratio (e.g., ). However, recent studies suggest that the densities are higher, around ∼ 10^3 - 10^5 cm^-3, a regime where the [SII]-based method cannot be used <cit.>. 

So far, only spatially-resolved spectroscopic observations, either by long-spit observations or using integral field units (IFUs), could be used to determine the outflow extent, r_out, and thus determine the mass and energetics of outflows (e.g., ). While such observations offer a more detailed view of the outflows, they are observationally demanding. Typically, only the most extreme outflow cases are followed-up with IFUs (e.g., ), resulting in biased samples that might not represent more typical outflow cases. Furthermore, IFU-based observations suffer from various uncertainties and systematics, such as projection effects and beam smearing, that can significantly affect the derived kinematics and extents of the outflows (e.g., ).

In <cit.>, we argued that the outflowing gas in active galaxies contains dust, which is heated by the central AGN, and emits in mid-infrared wavelengths. We examined the infrared spectral energy distribution (SED) of thousands of type II AGN, and showed that emission by such dust is detected in many systems that host ionized gas outflows. Our SED fitting provides the properties of the dust, including its distance from the central source. Since this dust is mixed with the outflowing gas, our method provides the luminosity-weighted (and therefore the mass-weighted) location of the outflow. This infrared emission is not subjected to various systematics affecting ground-based optical IFU observations, such as beam smearing, projection effects, and dust extinction that affects the receding part of the outflow. The estimated location of the dust, r_dust, can be combined with 1D spectroscopic observations to determine the mass and energetics of ionized outflows in hundreds of systems. 

In this work we focus on a subset of the sample presented in Paper I, in which ionized outflows are detected in many optical emission lines ([OIII], Hβ, [OI], [NII], Hα, and [SII]), and there is a clear detection of a dusty wind component at mid-infrared wavelengths. This combination allows us to accurately determine key properties of the gas in the wind, such as dust reddening, ionization state, and electron density, and thus constrain the mass and energetics of the winds. We describe our sample in section <ref>, and derive the spectral properties of the objects in section <ref>. We present new general expressions to derive the ionization parameter, electron density, and line emissivity of the ionized gas in section <ref>. We then estimate the outflowing gas mass, mass outflow rate, and kinetic energy of the winds in section <ref>. We discuss our results in section <ref>, and conclude in section <ref>.



§ SAMPLE SELECTION


The sample analyzed in this paper has been described in detail in paper I. It consists of type II AGN with spectroscopically-detected outflows, for which a dusty wind component is detected in mid-infrared wavelengths.  For clarity, we summarize here the main properties of the sample.

We start with the publicly-available catalog: AGN Line Profile And Kinematics Archive (ALPAKA; ), which provides emission line measurements for a sample of 24 264 optically-selected AGN from SDSS DR7 <cit.>. <cit.> performed a multi-component fitting to the emission lines [OIII]λ5007Å, [NII]λ6584Å, Hα, and Hβ, using continuum-subtracted spectra. The fit included one to three kinematic components to each of the emission lines, and the presence of an additional broader component in [OIII]λ5007Å  was interpreted as a spectroscopic signture of an ionized outflow. We selected systems that are classified in the ALPAKA catalog as type II AGN in the redshift range 0.05 ≤ z ≤ 0.15, with spectroscopically-detected outflows.

We cross-matched this sample with the MPA-JHU catalog, which provides stellar mass and star formation rate (SFR) measurements of SDSS galaxies <cit.>. We also used optical and infrared (IR) photometric data as follows. We used the riz optical photometry by the SDSS, the JHK near-infrared photometry from the 2MASS All-Sky Extended Source Catalog <cit.>, and the W1-W4 mid-infrared photometry from WISE <cit.>. We used the astrometric cross-matches provided by the SDSS[https://skyserver.sdss.org/CasJobs] to cross-match between systems observed by SDSS and systems observed by 2MASS and WISE.

In Paper I, we examined the infrared SEDs of 2 377 systems that show spectroscopic signatures of an ionized wind. The best-fitting SED included contributions from direct stellar light, which dominates the SED at optical and near-infrared wavelengths, torus and NLR dust emission, which dominate the SED at mid-infrared wavelengths, and dust in star forming (SF) regions which dominates the far-infrared wavelengths. The IR SED also included a new contribution from dust that is mixed with the outflow, is heated by the AGN, and emits at mid-infrared wavelengths. We found that 2 044 systems require the additional dusty wind component. Out of these, we were able to constrain the dust temperature, covering factor, and location in 1 696 systems. The initial sub-sample addressed in this paper consists of these 1 696 type II AGN. 





§ SPECTRAL PROPERTIES


The goal in paper II is to determine the level of ionization, the density, the mass outflow rate and the kinetic energy of the outflowing gas. To study the ionization state and density of the gas, we must detect the outflow in various emission lines: [OIII], [NII], [SII], [OI], Hα, and Hβ (see e.g., ). The outflowing gas is traced by broader emission lines that are blended with narrower emission lines coming from the stationary NLR. The broader emission lines are usually weaker than their corresponding narrow lines, and their characterization requires a careful emission line decomposition procedure. The ALPAKA catalog provides such a decomposition. However, <cit.> fit the continuum emission locally using a 5 degree polynomial, which they then subtract to obtain an emission line spectrum. Since they do not fit the spectra with a stellar population synthesis model, stellar absorption lines are not properly accounted for, and may affect the resulting emission line spectrum, particularly around the Hα and Hβ emission lines. The broader Hβ line is usually the weakest emission line in the spectrum, and its equivalent width is of the same order as the stellar Hβ absorption. 

To improve the accuracy of the broad Hβ emission line, we fit stellar population synthesis models to the 1 696 type II AGN sample, and performed our own emission line decomposition. We use Penalized Pixel-Fitting stellar kinematics extraction code (pPXF; ), which is a public code for extracting the stellar kinematics and stellar population from absorption-line spectra of galaxies <cit.>. It uses the MILES library, which contains single stellar population (SSP) synthesis models and covers the full range of the optical spectrum with a resolution of full width at half maximum (FWHM) of 2.3Å  <cit.>. We use SSP models with stellar ages that range from 0.03 to 14 Gyr, thus ensuring that we can properly describe spectra of systems with different star formation histories. The wide wavelength range of the SDSS spectra and the SSP models includes the Balmer absorption lines: Hζ, Hϵ, Hδ, Hγ, Hβ, and Hα, ensuring that the fitted stellar continuum around the Hβ absorption line is robust. The output of the code includes the relative weight of stars with different ages, the stellar velocity dispersion, the dust reddening towards the stars, and the best-fitting stellar model. Since we are not interested in the stellar properties of the galaxies in our sample, we only use the best-fitting stellar models and subtract them from the spectra.

The subtracted emission line spectra show varios emission lines: [OIII] λλ 4959,5007Å, Hα λ 6563Å, [NII] λλ 6548,6584Å, [OII] λλ 3725,3727Å, Hβ λ 4861Å, [OI] λλ 6300,6363Å, and [SII] λλ 6717,6731Å  (hereafter [OIII], Hα, [NII], [OII], Hβ, [OI], and [SII]). We model each emission line as a sum of two Gaussians - one which represents the narrow component and one that represents the broader, outflowing, component. We perform a joint fit to all the emission lines under several constrains: (1) we force the intensity ratio of the emission lines [OIII] λλ 4959,5007Å  and [NII] λλ 6548,6584Å  to the theoretical ratio of 1:3, (2) we tie the central wavelengths of all the narrow and broad lines so that the gas has the same systemic velocity, and (3) we force the widths of all the narrow lines to show the same velocity dispersion, and we do the same for the broader lines. 





Out of our initial sample, we detect a broad component with at least 2σ in [OIII], [NII], [SII], Hα, and Hβ, in 234 systems. We focus on these 234 systems throughout the paper. In figure <ref> we show the distribution of redshift, SFR, stellar mass, and dust-corrected narrow [OIII] luminosity for this subsample. The distributions are roughly similar to the distributions of these properties in the parent sample of 1 696 type II AGN. The two-sample KS test for the distributions give p-values of 0.53, 0.04, 0.50, and 0.03, for the redshift, SFR, stellar mass, and [OIII] luminosity respectively.

In figure <ref> we show the narrow and broad emission lines on a line-diagnostic diagram <cit.>. We show three separating criteria. The first is a theoretical upper limit which separates starbursts and AGN-dominated galaxies (, Ke01; black line). The second is a modified criterion which includes composite galaxies showing contributions from both SF and AGN (, Ka03; grey line). The third is a line from <cit.> to separate between LINERs and Seyferts (CF10; dashed black). Clearly, both the stationary (narrow lines) and the outflowing (broad lines) gas in most of the 234 systems are dominated by AGN photoionization, with the majority of the systems classified as Seyferts and some as LINERs. 

In the left panel of figure <ref> we show the FWHM of the broad emission lines. Since the fit requires all the broad lines to have the same velocity dispersion, the FWHM is the same for all the broad lines. In the right panel of figure <ref> we show the velocity shift of the broad emission lines with respect to the systemic velocity, defined by the centroid of the narrow emission lines. As in the left panel, this velocity shift is the same in all the lines.

Next, we use the measured Hα/Hβ flux ratios to estimate the dust reddening towards the two kinematic components. Assuming case-B recombination, a gas temperature of 10^4 K, a dusty screen, and the <cit.> extinction law, the colour excess is given by:

    E(B-V) = 2.33 × log [ (Hα/Hβ)_obs/2.85] ,

where (Hα/Hβ)_obs is the observed line ratio. In the left panel of figure <ref> we show the distribution of the colour excess for the narrow and broad lines respectively. For comparison, we show the distribution of the colour excess for the narrow emission lines in the ALPAKA catalog. Our narrow-line-based distribution of E(B-V), and the one from the ALPAKA catalog are roughly consistent, with somewhat larger number of sources with higher E(B-V) in the latter distribution. We attribute this difference to the different procedures used to fit and subtract the stellar continuum emission. For this reason, we do not compare the colour excess distributions in the broad emission lines.

In the right panel of figure <ref> we show the colour excess towards the broad lines versus the colour excess towards the narrow lines. Surprisingly, we find no correlation between the two. This is consistent with the findings of <cit.>, who analysed a sample of AGN-driven outflows in ultra luminous infrared galaxies. The figure clearly shows that the narrow and the broad lines suffer from different amounts of extinction (see also ). Various earlier studies assume the same amount of dust extinction, and use the estimated colour excess in the narrow lines to correct the broad line luminosities for dust extinction. Figure <ref> suggests that such a procedure is not justified. 



We derived dust-corrected line luminosities using the measured E(B-V). We propagate the uncertainties on the best-fitting parameters to obtain the uncertainties of the line luminosities, which are typically in the range 0.05–0.1 dex. We use the narrow emission line luminosities, which are based only on the narrow components from the multi-Gaussian fits, to estimate the bolometric luminosity of the AGN with two different methods. The first uses the narrow Hβ and [OIII] luminosities (equation 1 in ), and the second uses the narrow [OI] and [OIII] luminosities (equation 3 in ). In figure <ref> we show the distribution of AGN bolometric luminosity in our sample. The two methods give consistent estimates. One can see that our sample spans about two orders of magnitude in AGN bolometric luminosity, roughly from log L_bol = 43.5 to log L_bol = 45.5 erg/sec.





§ GAS PROPERTIES FROM PHOTOIONIZATION MODELLING


In the next section we estimate the mass outflow rates and energetics of the observed outflows, which depend on various properties of the emission line gas, such as its electron density and level of ionization. Studies usually assume constant fiducial values for the electron densities and line emissivities in all the objects in the sample (see e.g., ), which correspond to specific gas properties that might not be appropriate in a general case. In this section we provide more general expressions for the electron densities and line emissivities, which depend on the observed properties of individual sources. 

Many of the gas properties depend on the ionization parameter, U = Q(Lyman)/ 4 π r^2 n_H c, where Q(Lyman) is the number of the hydrogen-ionizing photons, n_H is the hydrogen density, and c is the speed of light. In section <ref> we show that the ionization parameter can be constrained using observed emission lines ratios. We then present a novel method to estimate the electron density in the gas, which is based on the estimated ionization parameter and on the known location of the gas, in section <ref>. Finally, we provide more general expressions of the Hα and [OIII] line emissivities, where the latter depends on the estimated ionization parameter, in section <ref>.



 §.§ Ionization parameter estimation


The expressions we provide are based on extensive photoionization modelling, which details are listed in appendix <ref>. We use a "single cloud" approach, where we calculate the properties of spherical shells of gas, with different density, pressure, column density, and distance from the central ionizing source. These models do not intend to replace more complex realistic models, that account for the geometrical distribution of clouds, with a range of densities and locations (see  for a more complete description of multi-component photoionization models). Nevertheless, they provide useful expressions that can be used in the more realistic scenarios. 

We focus on models with optically-thick geometrically-thin shells of gas, and do not model cases in which the gas is optically-thin and extends over regions that are as large, or even larger, than the distance from the central source (e.g., see ). The latter models require information on the gas spatial distribution, which is not available in this case. In addition, the outflowing gas can be a combination of matter-bounded (optically-thin) and radiation-bounded (optically-thick) clouds (see e.g., ). We did not consider such a combination since (1) the optically-thick models we examine provide a general and consistent description of the four emission lines studied here ([OIII], [NII], Hα, and Hβ), and (2) the contribution of matter-bounded clouds to Balmer line luminosities, and to the outflowing dust emission, must be small unless the covering factor of this component is considerably larger than the one found here for radiation-bounded gas. This is inconsistent with our measurements, and thus the presence of such a component must have a small effect on the expressions derived here.

We examined a large range in gas metallicity (from 0.5Z_⊙ to 2Z_⊙), ionization parameter (from log U = -3.8 to log U = -2), and several different slopes of the ionizing continuum (with mean energy of an ionizing photon of 2.56 Ryd to 4.17 Ryd). We compared the predictions of constant-density versus constant-pressure models, and found consistent results within the range of properties we considered. Our analysis suggests that the expressions we provide are general, and do not depend on specific assumptions made during the modelling.

The [OIII]/Hβ and [NII]/Hα emission line ratios are based on the strongest lines in the optical spectra of AGN. These ratios are insensitive to dust reddening and are widely used to distinguish between AGN and SF as the main source of photoionization <cit.>. The [OIII]/Hβ ratio is mostly sensitive to the ionization parameter in the gas, but it also depends on the gas metallicity. The [NII]/Hα ratio is mostly sensitive to the gas metallicity, but it also depends on the ionization parameter. For AGN-dominated systems, both ratios depend slightly on the shape of the ionizing continuum (e.g., see figure <ref> and ). 

We looked for an expression that ties the ionization parameter in the gas with the observed [OIII]/Hβ and [NII]/Hα line ratios. It is given by:

    log U = a_1 + a_2[log([OIII]/Hβ)] + a_3[log([OIII]/Hβ)]^2 + 
        + a_4[log([NII]/Hα)] + a_5[log([NII]/Hα)]^2,

where a_1, a_2, a_3, a_4, and a_5 are constants which are determined by a likelihood method. In appendix <ref> we provide details of this fitting process and derive best fitting values for the five constants. These are listed in table <ref>. Equation <ref> is only valid for ionization parameters in the range log U = -2 to log U = -3.8. The typical deviation in U, given the range of parameters listed in the appendix, and the two types of models (constant density and constant total pressure) is 0.11 dex. This includes a large range of metallicities and several SED shapes. 





 §.§ Electron density estimation


We use an independent method to estimate the electron density in the gas using optical emission lines. The method relies on the ionization parameter of the ionized gas, and the known location of the gas with respect to the central source. Using the definition of the ionization parameter, U, the electron density in the gas is given by:

    n_e≈ n_H = Q(Lyman)/4 π r^2 c U,

where Q(Lyman) can be estimated from the AGN bolometric luminosity and the assumed SED (see equation <ref> and related text in appendix <ref>). For the typical AGN in our sample, and the assumed SEDs, equation <ref> can be written as: 

    n_e≈ 3.2 (L_bol/10^45 erg/sec) ( r/1 kpc)^-2(1/U)  cm^-3,

where L_bol is the AGN bolometric luminosity and r is the known distance of the gas from the central source. The ionization parameter is estimated using the observed emission line ratios (equation <ref>). Equation <ref> can be used to estimate the electron density of both stationary and outflowing gas, as long as the AGN is the main source of ionizing radiation.



 §.§ Line emissivities estimation


In this section we provide more general expressions for the line emissivities in ionized gas. We focus on the Hα and [OIII] emissivities, γ_Hα and γ_[OIII], as these emission lines are typically used to estimate mass outflow rates of ionized outflows. Both emissivities depend on the electron temperature in the gas, with γ_Hα∝ T_e^-1 and γ_[OIII]∝ e^-1/T_e/√(T_e) (e.g., ). The electron temperature depends on the gas density, metallicity, and the ionization parameter. Our photoionization models suggest that the electron temperature in the ionized gas is in the range 10 000 K to 14 000 K, consistent with observations (e.g., ). We therefore select T_e = 12 000 K, resulting in Hα emissivity of γ = 3 × 10^-25 erg  cm^3  sec^-1. 

The [OIII] line emissivity is given by:

    γ_[OIII] = C_[OIII]× h ν_[OIII]×n(O^+2)/n(O)×n(O)/n(H),

where n_eC_[OIII] is the collisional excitation rate of [OIII], which depends on the electron temperature exponentially (C_[OIII]∝ e^-1/T_e/√(T_e); ), h ν_[OIII] is the photon energy, n(O^+2)/n(O) is the fraction of O^+2 ions, and n(O)/n(H) is the oxygen abundance relative to hydrogen. The emissivity depends on the temperature, primarily through C_[OIII], on the gas metallicity, which sets n(O)/n(H), and on the ionization parameter, which affects n(O^+2)/n(O). We note that different ionizing SEDs result in somewhat different n(O^+2)/n(O) ratios, but this effect is of secondary importance.

<cit.> provided an estimate for the [OIII] line emissivity which is based on equation <ref>. They assumed that n(O^+2)/n(O) ∼ 1, which is a reasonable assumption for their specific system, but does not apply in the general case. In particular, this assumption is not valid for most of the objects in our sample. Furthermore, they estimated the emissivity assuming solar abundance of oxygen and did not take into account that a non-negligible fraction of the oxygen atoms are depleted onto dust grains. Assuming electron temperature of 12 000 K, solar abundance of oxygen, n(O)/n(H) = 4.25× 10^-4, and a mean energy of an ionizing photon of 2.56 Ryd (model 2 SED; see appendix), the [OIII] line emissivity is given by:

    γ_[OIII]  (erg cm^3 sec^-1) ≈
       
    	2.06 × 10^-24×(n(O)/n(H)/4.25× 10^-4) × f(log U),    dusty gas
    
    	3.45 × 10^-24×(n(O)/n(H)/4.25× 10^-4) × f(log U),    dustless gas
    ,

where f(log U) represents the dependence on the ionization parameter, and is listed in table <ref>. The ionization parameter is estimated from observed line ratios in equation <ref>. The [OIII] emissivity for the dusty gas case is a factor of ∼ 0.6 smaller than the emissivity for the dustless gas case. This is a direct result of the depletion of oxygen onto dust grains, where only 60% of the oxygen atoms remain in the gas phase. In the next section, we use the [OIII] emissivity for the dusty gas case. Finally, we note that the constants in equation <ref> represent the average [OIII] emissivity in the models we considered, with the assumed ionizing SED and gas metallicity. While these give a general description of the ionized gas, they cannot replace full photoionization models, where the emissivities can be different by a factor of ∼1.5. 






§ MASS OUTFLOW RATE AND ENERGETICS


Having separated the stationary and outflowing gas components, we now turn to estimate the outflowing gas mass and its kinetic energy. The ionized gas mass in the outflow can be estimated either from the broad Hα luminosity, or from the broad [OIII] luminosity. The gas mass is given by (see e.g., ):

    M_out = μ m_H L_broad  line/γ_line n_e,

where μ is the mass per hydrogen atom, which we fix at 1.4, m_H is the hydrogen mass, L_broad line is the extinction-corrected broad line luminosity (L_Hα or L_[OIII]), n_e is the electron density in the outflowing gas, and γ_line is the effective line emissivity (γ_Hα or γ_[OIII]), which is given in section <ref>. The [OIII] line emissivity depends on the gas metallicity, which is typically unknown. Lacking additional information, we assume solar metallicity.

We also estimate the mass outflow rate, the kinetic power of the wind, and the kinetic coupling efficiency. The mass outflow rate is given by Ṁ_out = M_out/t_out, where t_out = r_out/v_out. The kinetic power in the wind is given by Ė_kin,  out = 1/2Ṁ_out v_out^2. The kinetic coupling efficiency relates the wind kinetic power to the AGN bolometric luminosity: ϵ = Ė_kin,  out/L_bol (see e.g., ). 

The above estimates require knowledge of the location of the wind, the electron density in the wind, and the effective outflow velocity. Since we use spatially-integrated (1D) spectroscopy and photometry, our estimates are based on luminosity-averaged properties, in particular the mass-weighted average properties of the wind. We describe in section <ref> our estimate of the wind location, which is based on the method presented in Paper I. We then describe various methods to estimate the electron density in the wind in section <ref>. Finally, in section <ref> we discuss our definition of outflow velocity.



 §.§ Effective wind location


In Paper I we argued that the outflowing gas in active galaxies contains dust, which is heated by the central AGN, and emits at mid-infrared wavelengths. We have shown that this emission component is detected in many type II AGN that host ionized gas outflows. This new component offers novel constrains on the outflow properties. Specifically, we used the dust temperature obtained from SED fitting and the AGN bolometric luminosity to estimate the distance of this dust from the central source. Our method provides luminosity-weighted, and therefore mass-weighted distance from the center of the galaxy. Since this dust is mixed with the outflowing gas, the method provides an estimate of the mass-weighted location of the outflow. In figure <ref> we show the distribution of outflow locations, traced by the dust infrared emission, for our sample. The distribution is similar to the distribution of our parent sample in Paper I. As discussed in Paper I, the uncertainty on the location is roughly 0.25 dex.





 §.§ Electron density


The electron density of the outflowing gas is a major source of uncertainty. The most robust estimates are obtained from full photoionization modelling of the outflowing gas (e.g., ). Such models require detailed, spatially-resolved spectroscopic observations, which are available only for a small number of objects, and hence they cannot be used to explore the entire range of population properties. Most earlier studies either assume a constant value for the electron density for all objects in the sample (e.g., ), or use the ratio of the weak [SII] emission lines to put limits on the density (e.g., ). Typical [SII]-based estimates are in the range 200–1000 cm^-3 (e.g, ). Here we discuss the merits and the limitations of this and other methods used to determine the gas density in the outflow



  §.§.§ [SII] line ratio




In the left panel of figure <ref> we show the [SII] density diagnostic diagram. The [SII]-based estimator is sensitive to electron densities which are close to the critical density of the [SII] λ 6717,6731Å  lines: about 1600 and 1.5 × 10^4 cm^-3 for a gas temperature of 10^4 K <cit.>. The [SII] ratio saturates at the low (n_e < 10^2 cm^-3) and high (n_e > 10^4 cm^-3) density regimes, and cannot be used as a density estimator in these ranges. We use the broad [SII] lines measured in section <ref> to estimate the electron density in the outflow, and show 20 randomly-selected measurements on the diagram. Due to the significant uncertainty in the emission line decomposition, the ratios are consistent with the entire range of possible ratios, from roughly 0.5 to 1.5. In particular, most measurements are consistent with the high density limit of about 10^4 cm^-3.

In the right panel of figure <ref> we show two distributions of the electron density based on the [SII] lines. The first is using the initial guess provided to the fitting function and without taking account the large uncertainties. This distribution peaks at roughly 10^2.5 cm^-3, similar to the estimates by others (e.g. ). Moreover, the exact location of the fit depends on the initial guess. We also show the distribution of [SII]-based densities after accounting for the measurement uncertainties. This, and the information provided in the left panel of the diagram, demonstrate that electron densities above about 10^3.5 cm^-3 cannot be constrained using the [SII] line ratio.

The [SII]-based electron density is subjected to an additional systematic. For the range of ionization parameters considered here, the [OIII] and Hα emission lines are emitted throughout most of the ionized cloud, while the [SII] lines are primarily emitted close to the ionization front in the cloud, where the electron density decreases with increasing depth into the cloud. Thus, the [OIII] and Hα lines tend to trace higher electron density regions, compared to the [SII] lines. The difference between the electron density traced by [OIII] and Hα and that traced by the [SII] becomes more significant for increasing ionization parameter, and for increasing hydrogen density in the cloud. Since the mass of the outflowing gas is estimated using the [OIII] or the Hα emission line luminosities (equation <ref>), the combination of [OIII] or Hα line luminosity with [SII]-based electron density may result in inconsistent gas mass estimates, with significant deviations in cases with high ionization parameters and large hydrogen densities.




  §.§.§ [OIII]/Hβ, [NII]/Hα, and r_dust


We now use the novel method presented in section <ref> to estimate the electron density in the outflowing gas. Given the coefficients in table <ref>, we can use equations <ref> and <ref> to estimate the electron density, where we define the distance of the outflowing gas to be r_dust, the mass-weighted location of the dust in the wind (see section <ref>). We propagate the uncertainties of the line ratios, the effective location of the dust, and the analytical expression of the ionization parameter, to obtain the uncertainty on the electron density, which is roughly 0.6 dex. In the right panel of figure <ref> we show the distribution of the electron densities in our sample using this method. The distribution is centered around n_e∼ 10^4.5 cm^-3, suggesting that the [SII]-based method, shown on the left side of the diagram, underestimates the true electron density in the outflowing gas by roughly 2 orders of magnitude. 



  §.§.§ Auroral and transauroral [OII] and [SII] emission line ratios


<cit.> discussed the uncertainties involved in estimating the electron density using the [SII] λ 6717,6731Å  emission lines. They pointed out that this estimator is not sensitive to n_e≳ 10^4 cm^-3. To estimate the electron density, they used the transauroral and auroral emission lines [SII] λ4068,4076Å  and [OII] λ7318,7319,7330,7331Å, which have critical densities around roughly 10^6 cm^-3 (e.g., ). Their method is based on the summed flux of each doublet, while the traditional [SII] method is based on the doublet flux ratio. As such, their estimated electron density does not suffer from the uncertainty involved in emission line decomposition, which is significant when using the traditional [SII] method. Their method is sensitive to a broad range of electron densities, from roughly 10^2 cm^-3 to 10^6 cm^-3. However, it is sensitive to reddening corrections.

<cit.> applied this method to study the AGN-driven outflows in 9 ultra luminous infrared galaxies (ULIRGs) observed with VLT/Xshooter. The exceptional quality of the VLT/Xshooter observations allowed them to perform emission line decomposition of the weak emission lines [SII] λ4068,4076Å  and [OII] λ7318,7319,7330,7331Å, and thus provide accurate estimates of the electron density in the outflowing gas. They found relatively high electron densities in the wind, from roughly log n_e(cm^-3) = 3.4 to log n_e(cm^-3) = 4.8  (see also ). <cit.> applied the method to an additional AGN, finding a remarkably high density of log n_e(cm^-3) ≃ 5.5. These estimates exceed the densities typically assumed or measured for ionized outflows in active galaxies. Our independent method confirms this in a much larger sample of non-ULIRG type II AGN. The densities we infer are somewhat larger than those found by <cit.>. A detailed comparison between the samples should be carried out with caution, since the two have very different size and consist of different types of objects, i.e. low luminosity type II AGN versus ULIRGs.

The two methods are sensitive to a large range of electron densities, and are expected to give consistent results. There are several differences between the method presented in <cit.> and our method. First, we assume that the distance of the gas from the central source is known, which is not required in the auroral and transauroral emission lines method. Furthermore, the method presented in <cit.> can be used to put constrains on the dust reddening in the gas, while we need to estimate it using the observed Hα/Hβ line ratio. However, the auroral and transauroral emission lines are very weak, and can only be detected in high-quality observations. In addition, this method requires a larger wavelength coverage, roughly from 3700Å  and 8000Å. Therefore, while less general, our method is expected to be more practical in many cases. 




 §.§ Velocity


There are several uncertainties associated with the estimate of the bulk velocity of the outflow (e.g., ). First, since we work with spatially-integrated spectra, the wind geometry is unknown and we cannot correct for projection effects. Second, dust-extinction may affect the observed emission line profiles, resulting in underestimation of the outflow extent and kinematics. 

Various studies define the bulk velocity in different ways (e.g., ). We follow the definition by <cit.>, since the properties of their sample, such as wind velocity, AGN bolometric luminosity, and redshift, closely match the properties of our sample. According to Karouzos et al. (2016b), v_out = √(v_shift^2 + σ^2 ), where v_shift is the velocity shift of the broad emission line centroid with respect to the narrow lines (see figure <ref>), and σ is the velocity dispersion of the broad emission lines. According to Karouzos et al. (2016b), using different definitions can result in mass outflow rates and kinetic powers which differ from the one used here by a factor of 1–3.

Having estimated the wind location and velocity, and the electron density in the outflow, we can now estimate the gas mass, mass outflow rate, and kinetic power of the wind. In figure <ref> we show the distribution of ionized gas mass, mass outflow rate, kinetic power, and the kinetic coupling efficiency for the objects in our sample. We show two estimates for these properties, one based on the dust-corrected broad Hα luminosity (black histogram), and one based on the dust-corrected broad [OIII] luminosity (grey histogram). The [OIII]-based estimates are consistent with the Hα-based estimates. This suggests that our assumption of solar metallicity is justified. 



Based on the broad Hα, we find that the median ionized gas mass in the outflow is roughly M_out∼ 10^4 M_⊙, the median mass outflow rate is Ṁ_out∼ 10^-2 M_⊙/yr, and the median kinetic coupling efficiency is roughly ϵ = Ė_kin/L_bol∼ 10^-5. These values are 1–2 orders of magnitude lower than typical estimates in type II AGN with similar bolometric luminosity (see e.g., ). These differences are the result of the much higher ionized gas density found in this work; roughly 2 orders of magnitude higher than in most previous studies. 




§ DISCUSSION


In paper I of this series we have shown that the IR SED of the dusty wind component in a large number of type-II AGN can be used to place constraints on the wind location. Here, in paper II, we analysed a sub-sample of 234 type II AGN with signatures of ionized gas outflows in their optical spectra, and with a detection of a dusty wind component in mid-infrared wavelengths, with well determined temperature, covering factor, and location. The detection of the outflow in multiple emission lines (Hβ, [OIII], [OI], [NII], Hα, and [SII]) allowed us to accurately determine key-properties of the ionized flow. In particular, we introduced a novel method to estimate the ionization parameter of the ionized gas. The combination of known ionization parameter and location allow us to use a novel method to estimate the electron density in the ionized wind, and thus to estimate its mass outflow rate and kinetic power. 

Our analysis of spatially-integrated spectroscopy and photometry provides the mass-weighted average properties of the ionized outflows. These estimates are not as accurate as estimates that are based on spatially-resolved IFU observations, but they provide the largest compilation of ionized outflow properties to date, and can be used to study the statistical properties of winds in type II AGN. Below, we discuss our results in the broader context of AGN feedback (section <ref>). We then elaborate on the neutral atomic gas which is associated with the observed ionized winds (section <ref>).















 §.§ AGN feedback




AGN feedback, in the form of powerful outflows, is invoked as a way to couple the energy released by the accreting SMBH with the ISM of its host galaxy, thus providing a possible explanation for the observed correlations between the masses of SMBHs and their host bulges (e.g., ). A number of such models have successfully reproduced the observed correlations, by requiring that a significant amount of the accretion energy of the AGN will be mechanically-coupled to the ISM of the host galaxy (∼5–10% L_bol; e.g., ). According to <cit.>, in a two-stage feedback scenario, the initial energy requirement may be 10 times lower. 

The coupling efficiency used in hydrodynamic cosmological simulations is often defined as the fraction of the total energy carried out by the nuclear wind close to its launching location. It remains unclear how much of this energy is transferred to the galactic-scale winds, due to various energy losses such as shocks, radiation, and work that is done against the gravitational potential of the host galaxy (e.g., ; see  for a comprehensive review). Furthermore, galactic-scale winds can be multi-phased, and it is often unclear what is the expected coupling efficiency of the different gas phases. Therefore, the coupling efficiencies required by numerical feedback models are not necessarily similar to kinetic coupling efficiencies derived from observations <cit.>. 

The typical mass outflow rate of the objects in our sample is roughly 10^-2 M_⊙/yr, corresponding to kinetic coupling efficiencies around ϵ = Ė_kin/L_bol∼ 10^-5, 3–4 orders of magnitude lower than the typical requirement in most numerical simulations. Our estimates are about 1–2 orders of magnitude lower than typical estimates in AGN host galaxies (e.g., ), which are based on electron densities in the range n_e∼ 10^2 - 10^3 cm^-3. Indeed, when we scale the mass outflow rates and coupling efficiencies found by <cit.> to have electron densities of n_e∼ 10^4.5 cm^-3, we find consistent results for AGN with the same bolometric luminosity. Furthermore, some of the earlier studies are based on IFU observations that target systems which show the most extreme outflow signatures. Our sample is less biased towards the most extreme outflow cases. 

In the left panel of figure <ref> we show the mass outflow rate versus the AGN bolometric luminosity for the objects in our sample. For comparison, we show the measurements by <cit.> for ionized outflows which, for the same L_bol, are much larger. Unlike <cit.>, we find no correlation between the mass outflow rate and the AGN bolometric luminosity in our sample. The middle panel shows the kinetic energy of the outflow versus L_bol, where we find no correlation as well. In the right panel of figure <ref> we show the mass outflow rate versus the star formation rate (SFR) in the host galaxy, where the latter is measured using the Dn4000 index (see details in Paper I). Surprisingly, the strongest ionized outflows in our sample are found in galaxies with high SFRs. This dependence might suggest that star formation is the main driver of the observed winds (see e.g., ). Since both the narrow and broad kinematic components in all the objects in our sample are classified as AGN using line diagnostic diagrams, it is possible that the process that photoionizes the stationary (NLR) and outflowing gas is not necessarily the process that produces the observed outflows. While the AGN in our sample populate only 3 out of the 4 quadrants indicated in the diagram, which might suggest that there are no active galaxies with low SFRs and strong winds, the small correlation coefficient (ρ = 0.16) and the large uncertainty on the mass outflow rates, prevent us from drawing more definite conclusions. 


The lack of correlation between Ṁ_out and L_bol might be related to the uncertainty in our mass outflow rate and bolometric luminosity measurements. Combining all the uncertainties discussed in section <ref>, we find the uncertainty on the mass outflow rate to be roughly 1 dex. The uncertainty on L_bol is 0.3–0.4 dex <cit.>. The combination of limited dynamical range (2 dex in L_bol) and the above uncertainties is not enough to explain the complete lack of correlation between Ṁ_out and L_bol. 

Finally, the lack of correlation with the AGN bolometric luminosity might be related to projection effects. If the AGN drives the observed winds through radiation pressure, we expect the gas to be accelerated within the AGN ionization cones, which are perpendicular to the line of sight in our type II AGN. Therefore, the measured wind velocity is underestimated by a factor of sinα, where α roughly represents the opening angle. Assuming that the scatter in the observed mass outflow rate is driven by scatter in the opening angle, such that logṀ_obs = logṀ_outflow + log(sinα), the scatter in the observed mass outflow rate is given by: Δ(logṀ_obs) = √(Δ(logṀ_outflow)^2 + Δ(log(sinα))^2). Under these assumptions, the lack of correlation between Ṁ_out and L_bol over 2 orders of magnitude implies that the scatter in sinα is roughly 1-2 dex, depending on the exact shape of the sinα distribution, and assuming logṀ_out/1   M_⊙/yr  = 1.3 logL_bol/10^45 erg/sec <cit.>. We find this very large range in sinα unlikely.

Given the above considerations, we suggest that the lack of correlation between Ṁ_out and L_bol is not entirely due to measurement uncertainties and/or projection effects, and is partially driven by the scatter in the intrinsic properties of the winds. This is at odds with the recent results by <cit.>, who found such a correlation assuming a constant value electron density in all the objects in their sample. Our new estimates of the electron density for individual sources both decrease the measured mass outflow rates and introduce a scatter. The lack of intrinsic correlation can be attributed to several different factors. First, different outflow histories may result in different spatial extents of the observed winds. Furthermore, since outflows are most likely to go through the path of least resistance, the observed winds may have very different geometries, depending on the mass distribution in their host galaxy. 



 §.§ Neutral atomic gas in the outflow


So far, we have focused on the warm ionized phase of the outflow, which is traced by strong optical emission lines. However, this gas can have large enough column density to be optically thick to the Lyman continuum radiation. For the ionization parameters found here, this requires a column in the range N_H = 10^19.5 cm^-2 to N_H = 10^20.5 cm^-2. The neutral gas at the back of such clouds is largely unconstrained by the observed optical lines, yet it may represent a significant fraction of the outflowing gas mass. The lower-ionization optical lines, such as [OI] λλ 6300,6363Å, [NII] λλ 6548,6584Å, and [SII] λλ 6717,6731Å, peak close to the ionization front. However, their emission drops significantly beyond it, and they cannot be used to place constraints on the fraction of neutral gas in the outflow. 

The NaID absorption line is commonly used to constrain neutral gas outflows in active galaxies (e.g., ). We examined the NaID absorption profile as a function of different outflow properties in our sample, and found no variation. However, the lack of variation does not necessarily imply that there is no neutral gas in the outflow, since the NaID absorption does not trace the same regions as the optical emission lines.

The dust mid-infrared emission can be used to constrain the neutral fraction of the outflows because the ratio of its infrared luminosity to the dust-corrected Hα luminosity, L(IR)/L(Hα), is sensitive to the ionization parameter and to the total column density (ionized and neutral). Since the dust is mixed with the outflowing gas, L(IR)/L(Hα) should not depend on the gas covering factor. We find that the observed L(IR)/L(Hα) ratios in our sample are broadly consistent with those expected from photoionization models, given the empirically-estimated ionization parameters, for hydrogen column densities of N_H∼ 10^20.5-10^21.5 cm^-2. This suggests significant amount of neutral atomic gas in most cases. It also supports the earlier conclusion (section <ref> that most of the outflowing gas is radiation bounded. The inferred hydrogen column density implies that the mass of the neutral gas is a factor of a few larger than the observed ionized phase, with important implications to the mass and energetics of the winds.

Unfortunately, the dust infrared luminosity, L(IR), was estimated in Paper I though SED fitting with significant uncertainties of ∼ 0.5 dex. This, and the uncertainty on U mentioned earlier, prevents us from reaching more definite conclusions about the exact amount of neutral atomic gas in individual cases. Furthermore, for a conical outflow, it is possible that the Hα emission originating in the receding side of outflow is completely obscured by dust in the ISM of the host galaxy. While the ISM of the host galaxy is optically-thick to optical radiation, it is optically-thin to mid-infrared radiation, and the observed L(IR) traces the full extent of the dust in the outflow. Therefore, the measured L(IR)/L(Hα) ratio is probably overestimated by a factor of ∼ 2.




§ SUMMARY AND CONCLUSIONS


In the first paper of this series <cit.>, we argued that the outflowing gas in active galaxies contains dust. The dust is heated by the central AGN and emits at mid-infrared wavelengths. We analysed the infrared spectral energy distribution of thousands of type II AGN and showed that this dust emission is detected in many systems that host ionized gas outflows. This infrared component offers novel constraints on the outflow properties, in particular, its location. In this work we focused on a subset of 234 galaxies from this sample, in which ionized outflows are detected in many optical emission lines ([OIII], Hβ, [OI], [NII], Hα, and [SII]), and there is a clear detection of a dusty wind component in mid-infrared wavelengths. The combination of the two allowed us to determine key properties of the gas in the wind, such as dust reddening, ionization state, and electron density, and thus constrain the mass and energetics of the winds. While our sample is based on spatially-integrated spectroscopy and photometry, and is thus less accurate than samples that are based on spatially-resolved observations, it forms the largest compilation of ionized outflow properties to date, and can be used to study the statistical properties of winds in type II AGN. Our results can be summarized as follows: 




  * We presented a novel method to estimate the ionization parameter in the outflowing gas. We also provided more general expressions for the Hα and [OIII] line emissivities, where the latter depends on the ionization state of the outflowing gas, and thus on the ionization parameter. The new line emissivities, which are then used to estimate the mass and energy of the ionized outflows, are expected to be more accurate than previously-used estimates.  


  * We presented a novel method to estimate the electron density in ionized gas, based on the strong optical line ratios [OIII]/Hβ and [NII]/Hα, and on the known location of the gas. We applied this method to estimate the electron density in the observed outflows and found electron densities of n_e∼ 10^4.5 cm^-3, which are about two orders of magnitude higher than most published estimates. We argue that the commonly-used method to estimate the electron density in the wind, which is based on the [SII] emission lines, underestimates the true densities by 1–2 orders of magnitude. Our estimates are more consistent with recent studies that are based on the density-sensitive auroral and transauroral [SII] λ4068,4076Å  and [OII] λ7318,7319,7330,7331Å  emission lines.


  * We found that the median ionized gas mass in the outflow is roughly M_out∼ 10^4 M_⊙, the median mass outflow rate is Ṁ_out∼ 10^-2 M_⊙/yr, and the median kinetic coupling efficiency is roughly ϵ = Ė_kin/L_bol∼ 10^-5. These values are 1–2 orders of magnitude lower than typical estimates in type II AGN with similar bolometric luminosity. This is a direct consequence of the high electron densities found in this study compared with the [SII]-based method. 


  * We found no correlation between the mass outflow rate of the ionized outflow and the AGN bolometric luminosity in our sample. This lack of correlation can be attributed to outflows having different locations, geometris, and different fractions of ionized gas in the outflowing clouds. Surprisingly, we found that the strongest outflows in our sample tend to occur in galaxies with high star formation rates, which might suggest that these are driven by supernovae in the host galaxy. Thus, the source that photoionizes the outflowing and stationary gas is not necessarily the source that drives the observed winds. 


  * Our study suggests the existence of a significant amount of neutral atomic gas at the back of the outflowing ionized gas clouds. We made an attempt to estimate the neutral gas fraction by combining the dust infrared luminosity with the Hα luminosity. We suggest that the neutral gas mass at the back of the outflowing clouds is a factor of a few larger than the observed ionized gas mass, which has significant implications for the estimated mass and energetics of such flows. The large uncertainties on the dust infrared luminosity prevent us from estimating the neutral gas mass in individual sources.






§ ACKNOWLEDGMENTS

We thank the referee for useful comments and suggestions that helped improve this manuscript.
This work is supported in part by the Israel Science Foundation grant 284/13 and by DFG/DIP grant STE 1869/2-1 GE 625/17-1.
The spectroscopic analysis was made using IPython <cit.>. We also used the following Python package: astropy[www.astropy.org/].

This work made use of SDSS-III[www.sdss3.org] data. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. 

mn2e








§ PHOTOIONIZATION MODELS


We model the central source using standard assumptions about AGN SEDs. The SED consists of a combination of an optical-UV continuum emitted by an optically-thick geometrically-thin accretion disk, and an additional X-ray power-law source that extends to 50 keV with a photon index of Γ = 1.9. The normalisation of the UV (2500Å) to X-ray (2 keV) flux is defined by α_OX, which we take to be 1.37.



The grid of models consists of a geometrically-thin, optically-thick, shells of dusty gas, with ISM-type grains and a density of n_H = 10^4 cm^-3. We used a range of gas metallicities: 0.5, 0.65, 0.8, 1, 1.5, and 2Z_⊙, and ionization parameter in the range log U = -3.8 to log U = -2. We consider four different shapes for the ionizing continuum SED, and we list their properties in table <ref>. Specifically, the SEDs have different slopes, with mean energies of ionizing photon of 2.56, 2.65, 3.15, and 4.17 Ryd. The ranges of metallicities, ionization parameters, and ionizing SEDs were chosen to fully cover the range of possible gas and AGN properties, and thus to provide a general description of the sources in our sample. For the typical AGN in our sample, we use SED 2 (mean energy of an ionizing photon of 2.56 Ryd), where the number of ionizing photons is given by:

    log Q(Lyman) = 10.06 + log L_bol,

where L_bol is the AGN bolometric luminosity, which can be estimated using narrow emission line luminosities (see ). Using the other SEDs (models 1, 3, and 4) result in a value which is different by roughly 10% from the value given in equation <ref>. Substituting equation <ref> into equation <ref> results in the expression given in equation <ref> in the main text.

We run a grid of 192 models with version 17.00 of  <cit.>. In figure <ref> we show the predicted [OIII]/Hβ versus [NII]/Hα by the different models. The colour represents the ionization parameter of the gas, and the size of the markers represents the metallicity, where the smallest markers correspond to 0.5Z_⊙, and the largest markers to 2Z_⊙. The shapes of the markers represent the ionizing SED, where the squares, stars, circles, and triangles correspond to SEDs with a mean energy of ionizing photon of 2.56, 2.65, 3.15, and 4.17 Ryd respectively. The predictions of these models are in agreement with predictions by other studies (e.g., ).

The models presented in figure <ref> are calculated under the assumption of constant density in the dusty gas. We also examined a similar grid of models, where we assume a constant total pressure (gas and radiation) in the gas. We find that the two sets of models give consistent optical line ratios for a given ionization parameter, and thus our conclusions do not depend on the assumption of constant density versus constant pressure in the gas. This is expected since we model gas clouds with ionization parameters of log U ≤ -2, where gas pressure is larger than radiation pressure. Therefore, these constant-pressure models do not represent a case of radiation pressure confinement (RPC), and are effectively constant-density models. RPC takes place when the radiation pressure is significantly larger than the gas pressure <cit.>, which is achieved for gas clouds with log U ≥ -2. The calculations of <cit.> show that such models can reproduce the spectrum of Seyfert galaxies with the most extreme [OIII]/Hβ line ratios (e.g., figure 10 in ), but cannot describe the bulk of the population. In particular, they do not cover a large region of the BPT diagram where many of our sources reside (figure 2). In such models, a range of ionization parameters (from log U ∼ 3 to log U ∼ -1) produce similar optical line ratios, and thus the line ratios cannot be used to place constrains on the average ionization parameter in the gas. In the rest of the section, we examine only constant-density models, since constant-pressure models give consistent results for our selected gas properties. 



Figure <ref> shows that the metallicity and the shape of the ionizing continuum are degenerate. That is, different combinations of gas metallicity and SED shape result in similar [OIII]/Hβ and [NII]/Hα line ratios. Therefore, these two properties cannot be determined from the optical line ratios alone, and require additional observables. On the other hand, the colour of the markers varies smoothly across the diagram, suggesting that the ionization parameter can be determined from the [OIII]/Hβ and [NII]/Hα line ratios, almost regardless of the gas metallicity or the shape of the ionizing continuum. 



The ionization parameter expression derived from our modeling is shown in equation <ref>. We find the coefficients a_1, a_2, a_3, a_4, and a_5, that minimize the residuals between the ionization parameter predicted by the photoionization models and the ionization parameter derived using equation <ref>. We list the best-fitting coefficients in table <ref>. In figure <ref> we show the predicted [OIII]/Hβ versus [NII]/Hα by the different models (same as figure <ref>), where the background colour represents the ionization parameter that is derived according to equation <ref>. Overall, there is a good agreement between the ionization parameters predicted by the models and those derived using equation <ref>. In figure <ref> we show the residuals between the two, which span the range of -0.2 to 0.2 dex. The standard deviation of the residuals is σ = 0.11 dex, which we adopt as the uncertainty of the analytic expression. 

To test the robustness of our analytical expression, and in particular the dependence of line ratios on the gas density, we considered a more complicated model, consisting of four geometrically-thin, optically thick, gas clouds. All gas clouds have solar metallicity and ISM-type grains. The central source is characterized by a mean energy of an ionizing photon of 2.56 Ryd. We varied the distance of each gas cloud and its hydrogen density, such that all models have the same ionization parameter of log U = -2.75. We set the hydrogen density in the clouds to be log n_H (cm^-3) = 5, log n_H (cm^-3) = 4, log n_H (cm^-3) = 3, and log n_H (cm^-3) = 2. 

As a first test, we examine the emission line ratios [OIII]/Hβ and [NII]/Hα of each individual cloud. Since the different clouds were defined to have the same ionization parameter, we expect their line ratios to be similar. The resulting line ratios differ from each other by less than 0.1 dex, supporting our suggestions that the hydrogen density does not affect the observed line ratios, for a constant ionization parameter. We used these line ratios and equation <ref> to predict the ionization parameter, and found log U of -2.68, -2.78, -2.87, -2.96 respectively. These values differ from the input ionization parameter, log U = -2.75, by at most ∼0.2 dex, consistent with the residuals shown in figure <ref>. 



Next, we consider the case of a combination of clouds with different densities and distances. We normalize the four models to have a similar covering factor. This is similar to using the mean of the four previous ionization parameters. We extracted the predicted emissivities of [OIII], [NII], Hα, and Hβ in the four different clouds, and summed them to obtain the total emissivity in each of the lines for the model. We then used the total [OIII]/Hβ and [NII]/Hα line ratios, and estimated the ionization parameter using the analytical expression from equation <ref>. The predicted ionization parameter is log U = -2.83. The difference between the input ionization parameter and the predicted is 0.08 dex, which is smaller than our adopted uncertainty. This suggests that the simple analytical form in equation <ref> is expected to hold for more complex scenarios, where the outflow is stratified, with a large range of density and location.

