Ionization Distributions in Outflows of Active Galaxies: Universal Trends and Prospect of Future XRISM Observations

The physics behind the ionization structure of outflows from black holes is yet to be fully understood. Using archival observations with the Chandra\HETG gratings over the past two decades, we measured an absorption measure distribution for a sample of outflows in nine active galaxies (AGNs), namely the dependence of outflow column density, $N_H$, on the ionization parameter, $\xi$. The slope of $\log(N_H)$ vs. $\log\xi$ is found to be between 0.00 and 0.72. We find an anti-correlation between the log of total column density of the outflow and the log of AGN luminosity, and none with the black hole mass and accretion efficiency. A major improvement in the diagnostics of AGN outflows will potentially occur with the launch of the XRISM/Resolve spectrometer. We study the ability of Resolve to reveal the outflow ionization structure by constructing the absorption measure distribution from simulated Resolve spectra, utilizing its superior resolution and effective area. Resolve constrains the column density as well as HETG, but with much shorter observations.


INTRODUCTION
Outflows are prevalent in Active Galactic Nuclei (AGNs). In the X-rays, one can observationally categorize three types of outflows. First are the relatively slow, ∼ 100 km s −1 , outflows detected through their numerous absorption lines in the soft X-ray spectra of Seyfert galaxies (e.g., Kaastra et al. 2000). Second, are the ultra-fast outflows (∼ 0.1c) detected mostly via highly blue-shifted Fe-K lines (Tombesi et al. 2010;Nardini et al. 2015). Third, are the high column-density outflows that almost completely obscure the soft X-ray AGN emission (Kaastra et al. 2014a). The present work will focus on the Seyfert outflows and on their unique property of a broad distribution of ionization, which has not been identified in the other two types.
Seyfert outflows are known for their broad ionization distribution, from the most highly charge species to neutral and near neutral species Behar et al. 2001). The ionization of a steady-state photo-ionized outflow can be described by its ionization parameter, which represents the balance between ionizing photons and recombining electrons. In the X-rays, one usually uses where L is the ionizing luminosity between 1 -1000 Ryd, n denotes the hydrogen number density, and r is the distance of the absorber from the ionizing source at the AGN center. ξ therefore has dimensions, in cgs units that is erg s −1 cm. Seyfert outflows span a few orders of magnitude in ξ; henceforth we will refer to −1 < log ξ < 4 in these units. Previous works found several discrete ionization components in Seyfert outflows and discussed whether they represent a continuous distribution. Blustin et al. (2005) published a survey of 23 AGNs to obtain insights into the absorbers and found that most outflows have multiple ionization phases, but there was no conclusion regarding whether the ionization structure is discrete or continuous. The average number of ionization components in their sample is two, but two AGNs with high quality spectra were modeled with three components. Detmers et al. (2011) studied the outflow of Mrk 509 with a model that spans four orders of magnitude in ξ, but preferred discrete ξ components over a continuous distribution. McKernan et al. (2007) analyzed X-ray spectra of 10 Seyfert outflows. Their best fits had two or more ξ-components, and five of the sources were modeled with three or more components. Laha et al. (2014) fitted XMM-Newton/RGS spectra of 17 Seyfert outflows using 1 -3 ξ-components for each one, to obtain an ionization distribution common to all Seyfert outflows.
The column density N H = ndr of the outflow distributed over different ξ values is termed the Absorption Measure Distribution (AMD, Holczer et al. 2007): Since the AMD involves n and r, its measurement can be invaluable for obtaining the density profile n(r) of the outflow (Behar 2009), which is not otherwise accessible via absorption measurements. This profile, in turn, may hold the crucial hint to the launching mechanism of Seyfert outflows, whether magnetic, radiative, or thermal (or a combination thereof). Different launching mechanisms could produce different AMDs, to be compared with the observed ones. For example by parameterizing the AMD as a power-law, with a slope a. Magnetohydrodynamic (MHD) accretion disk winds predict an AMD slope that depends on the scaling of the magnetic field with radius (along the line of sight) B ∼ r q−2 , leading to or q = (2 + a)/(2 + 2a) (Fukumura et al. 2010). Another physical effect is radiation compression of a hydro-statically held slab of photo-ionized gas (Stern et al. 2014) which results in a broad and flat (a ≈ 0) AMD. As the radiation is absorbed in the slab, its pressure drops, thus the gas pressure increases. Since the temperature also drops (less radiative heating) the density rises sharply, hence compression. Finally, numerical models of thermally driven winds have been able to produce AMDs that are flat, at least above log ξ = 2 (Waters et al. 2021, Fig. 13 therein).
In the present work, we model the spectra of a sample of nine Seyfert galaxies, using the Chandra/HETG spectrometer, seeking patterns in their outflow AMDs, and aiming to find a universal AMD profile. The specific goal is to study AMDs with commonly available spectral analysis tools, unlike the custom method of (Holczer et al. 2007). Our secondary objective is to assess the abilities of the future calorimeter spectrometer XRISM/Resolve (Tashiro et al. 2020) for studying outflow AMDs. We examine advantages of the superior sensitivity of XRISM/Resolve compared to Chandra/HETG. A comprehensive view of AGN outflows with the Hitomi calorimeter, similar to that of XRISM/Resolve was published by Kaastra et al. (2014b).

Sample and Data
In order to study the AMD of AGN outflows, we selected all of the observations in the Chandra/HETG archive with visibly detectable absorption lines that can be ascribed to outflows. The targets and observations used in the present paper are detailed in Table 1. For each object, we used only observations from the same year in order to minimize variability of the absorber. We checked one such example, the NGC 5548 HETG observations from 2000, 2002 and 2019, all with relatively high photon counts. All AGNs in the sample have a Chandra/HETG exposure of at least 190 ks and at least 160,000 source counts. The rightmost column in Table 1 lists previously published works with multi-ξ models for these observations, and for each target. Despite the variety of approaches, the models share some common attributes, as follows. They include 2 − 6 ξ-components (e.g., Netzer et al. 2003). Most of the models identify at least two kinematic components (e.g., Steenbrugge et al. 2005a); half of the objects have a slow outflow with velocity 100 − 500 km s −1 and a fast one with velocity > 1000 km s −1 (e.g., Silva et al. 2016). Out of these four objects with fast outflows, three are only apparent in the high-ξ components. Two objects have no high-ξ components nor fast outflows (Gupta et al. 2013;Detmers et al. 2011). In contrast with this variety of models, the present analysis uses a fixed set of ionization components, in order to obtain a uniform AMD structure. Table 2 lists the physical parameters of the AGNs. The unabsorbed X-ray luminosity L X in the 0.5 -10 keV band is measured from the present data, while the distance and black-hole mass M BH are taken from the literature. It can be seen that all AGNs are at low z ≤ 0.034, yet M BH spans two orders of magnitude between 1.6 × 10 6 M − 1.4 × 10 8 M , as does the accretion efficiency 0.003 ≤ L X /L Edd ≤ 0.3. Therefore, if the outflow AMD depends on any of these, we expect to be able to find hints to that dependence in our sample. The simultaneous fitting of all spectra of each target assumes the absorber did not vary (or it yields results for a mean outflow of that target) during the temporal period of observations. Since the overwhelming majority of observations for a given target were obtained within a year, or even much less, this seems to be a reasonable assumption for Seyfert outflows, also supported by the previous analyses. The unresolved (narrow) absorption lines suggest the absorber is relatively far from the nucleus, and would not vary over such short time scales. However, there are reports of absorber variability (Krongold et al. 2005(Krongold et al. , 2007.

Spectral Model
In each observation, first diffraction order (±1) of both the HEG and the MEG gratings were used. The total exposure time for all targets is at least 190 ks, and the total photon counts between 1.5 and 20Å are greater than 166320 for all objects. We prepared different XSTAR (Bautista & Kallman 2001) tables for each target, with the appropriate power-law spectral index Γ, a soft excess component when needed, and using solar abundances. We use a turbulent velocity of 100 km s −1 , except for the broad lines of NGC 3516 (Holczer & Behar 2012) and the log ξ = 4 component of MCG -6-30-15 (Holczer et al. 2010). The grid includes 7 logarithmic steps between logN H = 10 18 − 10 24 , and 9 logarithmic steps in log ξ from -1 to 4.
We fitted the spectra using the the Xspec implementation of C-statistic (Cash 1976), since some of the bins in the data contain only few photons. The above tables were sufficient to well-fit (C-stat/dof ∼ 1) all spectra with components describing the continuum, and several absorption ionization components N H (ξ).
We fitted the spectra of each target between 1.5-20Å using Xspec (Arnaud 1996). The continua are modeled as a power-law with a soft excess (when needed) that rises above the power-law around 15Å. This soft-excess is commonly modeled as a black-body component, which is satisfactory for our purpose of characterizing the absorber. The neutral galactic absorption was also taken into account with a fixed absorption component (Table 2). We modeled the absorber with six components at fixed log ξ i values of -1, 0, 1, 2, 3, and 4. The exception is NGC 3783, where we added two more components at 0.5 and 2.5, required by its superb spectrum.
The velocity of each ξ i -component was fitted individually. The outflows of NGC 3516 and MCG -6-30-15 have more complex kinematic structure. NGC 3516 requires at least two kinematic components. One between -650− -50 km s −1 and the second at −1650 ± 50. Both kinematic components have all six of the ξ i -components. The turbulent velocity in each outflow was kept fixed for all ξ i -components. The outflow of MCG -6-30-15 also has two kinematic components; the log ξ = 4 component is at −1800 ± 80 km s −1 , while the lower-ξ-components are between −350 − −130 km s −1 .
Finally we allowed the fit to include the following narrow emission lines -1.78Å (Fe Lyα), 1.87Å (Fe +24 f ) , 1.94Å (Fe Kα), 13.45Å (Ne +8 r), 13.70Å (Ne +8 f ). NGC 4151 included the following emission features as well: 7.13Å (Si Kα), 9.17Å (Mg +10 r), 9.31Å (Mg +10 f ), 12.13Å (Ne +9 Lyα), 16.01Å (O +7 Lyβ), 16.8Å (O +6 RRC). The best-fit column densities N H (ξ i ) of each AGN are subsequently used to build the AMD, by plotting log N H as a function of log ξ i . For each target, we obtain the AMD slope by fitting a linear regression in log space to log N H vs. log ξ i , taking into account N H (ξ i ) uncertainties, which are calculated by Xspec with the standard 90% confidence. For NGC 3516 we built the AMD based on the values of the slower outflow velocity. The uncertainties on each N H (ξ i ) for this purpose are taken as a mean of the lower and upper statistical uncertainties extracted from the spectral fit. The total column density, N tot H = Σ i N H (ξ i ), is obtained for each AGN. The reported uncertainty of N tot H reflects the maximal lower and upper limits.

AMD Reconstruction
We subsequently take the best-fit Chandra/HETG model to simulate XRISM/Resolve spectra using its anticipated response matrix (Ishisaki et al. 2018). We simulate a standard exposure time of 100 ks for all targets. For Ark 564, Mrk 509, IC 4329A, NGC 4051, NGC 4151 and NGC 5548, the highest ξ component, log ξ = 4 is not well constrained with Chandra/HETG. Since the simulation randomly draws photons, it is meaningless to simulate such low column densities that the spectrometer cannot detect. Thus, in the XRISM/Resolve simulations of these targets we take the HETG upper limit. In order to get an idea of the ability of XRISM/Resolve to constrain the AMD, we then fitted the same model to the simulated spectra, followed by constructing an AMD for each object. Next, we carried out the linear regression and calculated N tot H to be compared with the Chandra/HETG results.  Figure 2 shows the same part of the spectrum in the XRISM/Resolve simulation, plotted in energy instead of wavelength to best demonstrate the performance of the future spectrometer.

Column Densities
The best-fit column densities for each ξ i -component of each target, as well as N tot H , are listed in Table 3. The best-fit continuum parameters and cstat/dof values are listed in Table 4. As expected, the best-fit column densities of Chandra/HETG observations and XRISM/Resolve simulations agree to within 60%, and are consistent within the uncertainties. However, the fractional uncertainties of the XRISM/Resolve fits are generally smaller by up to a factor 4, and especially in the highest ξ-components. The improvement of XRISM/Resolve is most noticeable in the log ξ = 4 component of NGC 4151 and NGC 5548. We also examined the difference between Resolve simulations when taking the best-fit HETG value, versus using the upper limit for log ξ = 4. The results are very similar for Ark 564 and NGC 4051, where the column density of the log ξ = 4 component is not constrained neither in HETG nor in Resolve. For Mrk 509 and IC 4329A the column of the log ξ = 4 component is again not constrained when the original column density value is used, as opposed to the upper limit value.
For the most part, our results of N tot H and outflow velocities approximately agree with previous works, apart from the two outstanding following exceptions. Gupta et al. (2013) fitted the outflow of Ark 564 with two ξ-components, both at a velocity of ∼ −100 km s −1 , at log ξ = 0.39 ± 0.03 and at −0.99 ± 0.13. Their total columns at log N H (cm −2 ) = 20.94 and 20.11, respectively, are much lower than what we find, likely because Gupta et al. (2013)     spectra above ∼ 9Å, therefore finding no high-ξ components. The XMM-Newton/RGS spectrum of the outflow of NGC 4051 was fitted by Silva et al. (2016) with four ionization components: log ξ = 0.37 ± 0.03, 2.60 ± 0.10, 2.99 ± 0.03, and 3.70 ± 0.04 corresponding to outflow velocities of, respectively, −340 ± 10, −530 ± 10, −4260 ± 60, and − 5770 ± 30 km s −1 . We do not identify the three high-velocity components in the HETG spectra. For NGC 3783 we found two velocity ranges between the ξ i -components, within the range of ionic velocities reported in Kaspi et al. (2002).

AMDs
Figure 3 presents all the AMDs and their fitted slopes, for both spectrometers, and including N H (ξ i ) uncertainties. It can be seen from the figure and Table 3 that the highest column densities occur in the high-ξ i components. In some AGNs it is difficult to identify these components using the HETG observations, therefore Resolve spectra could have a meaningful advantage. The values of AMD slopes are detailed in Table 5. The HETG slopes range between 0.00 -0.72, and between -0.07 -0.85 including the uncertainties. Apparently, a slowly increasing AMD is a common property of Seyfert outflows. Apart from Mrk 509, all slopes are tightly constrained to ∼ 0.1 or better. The anomalously high N tot H in NGC 4151 is due to its intermittent obscuration (see Section 1). Its high N H (ξ i ) > 10 22 cm −2 values at log ξ = −1 and log ξ = 1, which none of the other AGNs features (see Table 3) likely have little to do with the outflow.
The consistency between AMD slope values of HETG and Resolve (Table 5) implies one can accurately constrain N H (ξ i ) and reconstruct the AMD based on XRISM/Resolve, with observation times half as long as those of Chandra/HETG, or shorter. Our results in Table 3 show that XRISM/Resolve has the sensitivity to measure N H (ξ i = 10 4 ) down to ∼ 10 21 cm −2 , with a 100 ks exposure. Column densities of lower-ξ i -components can even be constrained as low as ∼ 10 20 cm −2 .  In order to demonstrate the effect of the different ξ i -components on the XSTAR based model, we plot each of them individually in Figure 4, for the model of NGC 3783 which has an abundance of resolved features. The Figure reveals several attributes of the model. The lower ξ i -components absorb mainly the continuum, without many absorption lines. This is most evident for log ξ = −1. Higher ξ values absorb the continuum less and less; log ξ = 4 absorbs virtually no continuum. The column density of each ξ i -component is predominantly determined by the imprint of its continuum slope (Figure 4). The Fe-M UTA, at 16 − 17Å appears mainly in the log ξ = 1 component. Apparently, this XSTAR model does not fit the conspicuous UTA of NGC 3783 properly, see residuals in Figure 1. Since there is overlap in the lines between components, there is much freedom for the fit to lower or raise column densities of adjacent ξ i -components, which is reflected in the uncertainties. The main driver of C-stat minimization is therefore the continuum.

AMD and other AGN Parameters
It remains to be understood what drives the outflow properties. We examine the connection between the AGN physical properties and N tot H , meaning the sum of column densities of all log ξ i components. In Figure 5 we show a relation between logN tot H and log L X , which appear to be anti-correlated. The Pearson's correlation coefficient r and Spearman's rank coefficient r s are -0.22 and -0.38, and p-values 0.56 and 0.31, respectively. However, by removing NGC 4051, the anti-correlation improves dramatically (-0.83 and -0.83, with p-values 0.01 and 0.01, respectively). Note that the two main groups of AGNs in Fig. 5 differ by their log ξ = 4 column density (high N tot H low L X ) or lack thereof (low N tot H high L X ). This separation may turn out to be a smooth transition, once XRISM/Resolve better constrains this component. Conversely, there seems to be no clear relation between M BH and N tot H , as can be seen in Figure 6. The coefficients there are r = 0.17, r s = 0.22, with p-values 0.66 and 0.58, respectively. Therefore, there is also no clear relation between N tot H and L X /L Edd (∼ L X / M BH ), see Figure 7. The coefficients there are r = −0.40, r s = 0.36, with p-values 0.29 and 0.36, respectively. We also examine the connection between the AMD slope and the above AGN parameters. Since the slopes span a narrow range, we do not expect to find a strong relation. Indeed, in both cases we find no significant relation between the AMD slope and these AGN parameters.

DISCUSSION
Using archival Chandra/HETG grating observations of nine AGNs we constructed the AMDs of their outflows using at least 6 pre-defined ξ i -components ranging over −1 < log ξ i < 4. Mrk 509 requires three components, while the other three only provide upper limits, while NGC 3783 requires seven components and one upper limit. This is a somewhat broader range than reported by McKernan et al. (2007). The reason is that we assumed all pre-defined components, while McKernan et al. (2007) sought the minimal number of components that provided a satisfactory fit.
The best-fit slopes of the various AMDs in the Chandra/HETG spectra span a range of 0.00 -0.72 (-0.07 -0.85 with uncertainties), which is consistent with the range of 0.0 -0.4 reported in Behar (2009). Steenbrugge et al. (2005b) found a slope of 0.40 ± 0.05 for the outflow of NGC 5548, which is inconsistent with our results. This may be due to the fact we used later observations of the outflow as well as the ones used in Steenbrugge et al. (2005b). In a comprehensive linear regression for the ξ-components of 17 different Seyfert outflows (including all of the present ones except NGC 4151), Laha et al. (2014) find an AMD slope of 0.31 ± 0.06, which is consistent with the present slopes. This suggests a ubiquitous AMD shape for AGN outflows of a shallow positive slope (< 0.72).
The present slope range of 0.00 -0.72, corresponds in the MHD self similar solutions of B ∼ r q−2 (Fukumura et al. 2010), to q = 0.79 − 1.0, or approximately B ∼ 1/r in all outflows. Our results are marginally consistent with those of Stern et al. (2014), a ∼ 0.03. XRISM/Resolve spectra will allow us to measure the high-ξ components with smaller  Table 2. uncertainties, providing a more definitive AMD slope value. The MHD outflow model is scalable with q, but other models will be confronted with these refined AMDs.
The anomalous high column densities of NGC 4151 (N H (ξ i ) > 10 22 cm −2 ) in low-ξ components are greater than those of any other AGN (Fig. 3). The soft X-rays of NGC 4151 are often heavily absorbed (George et al. 1998;Kraemer et al. 2005) by K-edges of light elements that are likely not related to the steady outflow. We analyzed the 2002 spectra where NGC 4151 was in a high flux state. Nevertheless, residual continuum absorption results in these high N H (ξ i ) values for low-ξ i components, although there are no clear absorption lines above 10Å (see also Kraemer et al. 2005).
Previous works (Holczer et al. 2007;Laha et al. 2014) find a gap in the AMD, between log ξ = 0.5 − 1.5, and suggest this could be a universal feature due to thermal instability. Waters et al. (2021) show that in thermally driven winds, the buoyancy of gas clumps, and their disintegration within thermally unstable regions can remove this gap from the AMD. The present method of a rigid ξ-grid does not provide unambiguous evidence for this gap, although marginal evidence can be seen in the AMDs of NGC 3516, NGC 3783, NGC 4151 (Fig. 3).
Since the AMD slopes are relatively similar between AGNs, they point to a basic physical attribute of the outflows which is universal. On the other hand, the large dispersion in N tot H , allows us to correlate it with the AGN fundamental properties. We find that the logN tot H plausibly anti-correlates with log L X . A similar anti-correlation is found in the SUBWAY quasar sample between N HI and L bol (Mehdipour et al., in preparation). Conversely, there is no correlation between N tot H and M BH or L X /L Edd . Radiatively driven winds are actually expected to drive more mass with luminosity. The anti-correlation with L X thus might suggest that the X-ray flux moves gas out of the line-  Table 2. of-sight, leading to lower column densities. An alternative explanation is that high-L X AGNs totally ionize the wind, thus hiding its most ionized components. However, the similar AMDs of all AGNs and specifically the lack of increasing AMD slope with luminosity, suggests that outflows of low-L sources are as ionized as those of high-L ones. Blustin et al. (2005) found a possible, weak correlation between N tot H and the bolometric luminosity (cf. Fig. 5 in Laha et al. 2016), which hangs on four luminous quasars. Two of them, PG0844+349 and PG1211+143, have ultra-fast velocities, quite different from the Seyfert outflows (e.g., Laha et al. 2014). The two other sources, IRAS 3349+2438 (L X ∼ 6 × 10 44 erg s −1 and N tot H = 1.2 ± 0.3 × 10 22 cm −2 , Holczer et al. 2007), and MR 2251-178 (L X = 1.7 − 5.2 × 10 44 erg s −1 and N tot H = 3.2 − 6.3 × 10 21 cm −2 , Kaspi et al. 2004), would strengthen the anticorrelation of Fig. 5.

CONCLUSIONS
Following the uniform analysis of a sample of nine Seyfert outflows we reach the following conclusions: • The AMD slope, a proxy of the ionization distribution in the outflow is relatively flat. This slope is found to be a universal characteristic of the outflows, indicating a common wind-launching mechanism, or micro-physics that is not related to global properties of the AGN.
• The log of the total column density in the outflow N tot H anti-correlates with the log of the X-ray luminosity L X , perhaps indicating that high-L X sources clear absorbing gas from the line of sight.