Physical properties and the variability mechanism of the He i outflow in NGC 4151

The six spectra were observed at optical and near-IR wavelengths in order to span the $2^{3}S$ He i emission and absorption features resulting from two of the components, at 3889 Å and 10 830 Å. In both wavebands small additive wavelength shifts were made to account for differences in wavelength calibration across the various epochs by aligning the spectra to the wavelengths of strong skylines prior to correction for redshift. Before any analysis of the spectra were undertaken, small variations in the effective resolution were removed by identifying the spectrum across both sets of observations with the broadest narrow features. All other spectra were then convolved with a Gaussian of unit area and appropriate FWHM to match the resolution of this spectrum ($\sim$400 km s^-1) while preserving the area of emission and absorption features.

Reconstruction of the unabsorbed spectrum is required in order to measure the depth of the absorption features. This process is complicated by the overlap of emission from several different sources, namely the AGN power-law continuum, the broad line region (BLR) including a continuum of blended broad Fe ii emission, the NLR, the torus and the host galaxy. To deal with the relatively featureless and continuous aspects of this emission, the aforementioned Fe ii continuum together with the host galaxy starlight and AGN power-law contributions were fitted simultaneously to the same relatively line-free regions of the spectra at each epoch using the iraf specfit package, which was also used in all subsequent spectral fitting described in this paper. For the near-IR spectra, an additional contribution from the dusty torus was modelled using a blackbody. The Fe ii model was scaled from the template obtained from I Zw 1 by Véron-Cetty et al. (2004), which was also Gaussian-blurred to approximate the width of the broad emission lines in the spectrum. The host galaxy contribution was built from the Sa and Sb-type spiral galaxy templates listed in the SWIRE library (Polletta et al., 2007). Epochs 2, 3 and 4 have estimated galaxy fluxes, derived from HST observations, listed in Landt et al. (2011). Galaxy templates were therefore normalised to match these fluxes prior to running the specfit fitting routine, minimising the problem of degenerate fitting. Epochs without such observations used an average value, corrected when necessary for differences in slit width. An example of the resulting power-law, galaxy starlight, dust torus and Fe ii emission models fitted to the observed spectra is illustrated in Fig. 1.

The simplest method of measuring outflow column densities is to assume that the absorber is unsaturated and covers 100 per cent of the emission region. This requires a simple integration across the absorber profile (Savage & Sembach, 1991) given by

Using the unabsorbed reconstructions in the spectral region of the He i $\lambda$20 587 line shown in Fig. 4, a value of the column density in this state can be estimated from the normalised spectral profile of the absorbed region (excluding the narrow emission line as the NLR is probably uncovered by the absorber). Since the He i transition at 20 587 Å is a singlet, only a lower limit for the column density can be obtained using direct integration of the profile, which will be an underestimate if partial coverage and saturation are present. However, we would expect these limits to be close to the true value given the shallowness of the absorber (Leighly et al., 2011). As mentioned in §3.2.2, column densities also need to be lower limits due to a lack of separation of the various emission sources at the wavelengths spanned by the singlet state absorber. The normalised absorption profiles are shown in Fig. 9. The derived column densities for both $2^{3}S$ (from the 3889 Å line and the 10 830 Å line) and $2^{1}S$ state helium are recorded in Table 3.

The unreliability of the optical measurements (due to the quality of the spectrum, the weakness of the continuum and the weakness of any 3889 Å He i absorber) is clear from these values for all but epoch 5, which is the only observation having a column density of greater than 2$\sigma$ significance (4.8$\sigma$). The strength of the near-IR absorber allows more reliable measurements to be made of absorber variability, which will be discussed further in §5. The $2^{1}S$ column density values are mutually consistent across the epochs, so no variability in this absorber is identifiable. This absorber is weak and the position of underlying emission is uncertain, so the lack of variability is not considered indicative.

In order to locate the epoch 5 triplet state absorber in ionization parameter, number density and column density space, grids of photoionization simulations were calculated. The photoionization simulations calculated for this paper were performed using cloudy version c13.02 (Ferland et al., 2013), for which the setup was similar to that described in paragraph 1 of §3.1 in Leighly et al. (2011). Simulations were available which utilised a triplet state He i column density of $n_{{\rm He}}(2^{3}S)$=3.67$\times$10¹⁴ cm^-2, well within the errors on the epoch 5 absorber determined from the PPC analysis. These simulations are those discussed in this section.

The ionization parameter ($U$) is a useful quantity for measuring the ratio of Lyman continuum photons to hydrogen atoms as follows

Balmer absorption was detected in NGC 4151 by Hutchings et al. (2002), arising from the H$\beta$ transition. In Leighly et al. (2011), the lack of Balmer absorption in the quasar FBQS J1151+3822 was used to set an upper limit on the n=2 hydrogen state column density. The spectral region of H$\beta$ at epoch 5 was examined to determine whether or not H$\beta$ absorption was present. This employed the multiple Gaussian reconstruction used previously (see Fig. 5) as the unabsorbed reconstruction. No absorber is apparent from a visual inspection of the spectrum and the reconstruction. This lack of absorption can be used to set limits on the allowed ionization parameter and total hydrogen column density of the absorber as described in §4.3.

Our calculations require two runs of the code over each grid. The first run is used to find the thickness of the absorbing cloud to be integrated through to reach the spectroscopically-determined column density of helium in the $2^{3}S$ metastable state at epoch 5. This allows the total hydrogen ($N_{{\rm H}}$) and $2^{1}S$ state helium column densities to be calculated. The second run integrates through to the $N_{{\rm H}}$ values determined from the first run to find the hydrogen level populations, which can then be used to determine if a particular gridpoint exceeds the allowed Balmer absorption. Spectroscopic constraints on the column density of helium in the $2^{1}S$ state are lower limits, meaning that cloudy-derived values must exceed the observed value (1.76$\times$10¹² cm^-2 from Table 3) to be considered valid.

The cloudy output provides the fraction of helium in the sum of all singlet states at each principal quantum state rather than the individual singlet states at different total azimuthal quantum numbers. Helium in the n=2 singlet state exists in both ${}^{1}S$ and ${}^{1}P$ configurations, meaning the latter of these must be taken into account to derive the fraction in the $2^{1}S$ state. This state is formed by the combined effect of He⁺ + e^- recombinations and collisional excitation from the $2^{3}S$ metastable state. Approximately one third of helium recombinations to a singlet state of n=2 or higher ultimately results in population of the $2^{1}S$ state, while the rest lead to the ground state. Therefore, using the recombination and collisional rate coefficients derived from Smits (1996); Benjamin et al. (1999) and Bray et al. (2000) respectively as listed in Osterbrock & Ferland (2006), the fraction of helium in the $2^{1}S$ state at each cloudy output zone can be found as follows

It was convenient to initially generate a “coarse” grid measuring 8$\times$7 using a range in ionization parameter $U$ of $-$2.5$\leq$log$U$$\leq$1.0 in steps of 0.5 dex and a range in hydrogen number density $n_{{\rm H}}$ of 3.0$\leq$log($n_{{\rm H}}$/cm^-3)$\leq$9.0 in steps of 1.0 dex. Having found the gridpoints which satisfy the conditions derived from the spectra, a “fine” grid could then be calculated to narrow down the allowed values in the parameter space neighbouring the allowed coarse gridpoints. The gas hydrogen column density at each gridpoint was set equal to log$N_{{\rm H}}$=24.5+log$U$ to allow the code to integrate through the relatively high-ionization region occupied by $2^{3}S$ He i. The grid was generated by performing the following steps: (i) The He i $2^{3}S$-state column density output from cloudy is compared to the epoch 5 spectroscopically-derived value. Gridpoints which produce values lower than this are ruled out, (ii) for those gridpoints that satisfy condition (i), the $2^{3}S$ He i number density in each output zone generated by the code is calculated by multiplying the fraction of helium in this state by the helium abundance and $n_{{\rm H}}$, (iii) the $2^{1}S$ He i number density in each output zone is found by calculating $F(2^{1}S)$ in each zone using Equation 7 and multiplying the result by the helium abundance and $n_{{\rm H}}$, (iv) the absorber thickness is found by performing the integral $\int$$n_{{\rm He}}(2^{3}S)$dr through each zone from the irradiated face of the cloud inward until the epoch 5 column is reached, with the radius at which this value is reached giving the total thickness, (v) this integration is performed over the same zones as in (iv) to find the $2^{1}S$ column density using $\int$$n_{{\rm He}}(2^{1}S)$dr, and (vi) the total thickness is multiplied by $n_{{\rm H}}$ for the gridpoint in question to find the hydrogen column density $N_{{\rm H}}$.

The absence of H$\beta$ absorption can be used to obtain an upper limit on the column density of hydrogen in the n=2 state and therefore rule out gridpoints where the cloudy models suggest this upper limit is exceeded. To do this, cloudy is run for a second time over the coarse grid, using a method very similar to that described in paragraph 3 of §3.1 in Leighly et al. (2011). The normalised epoch 5 profile shape for the He i $\lambda$3889 absorber is assumed to be the shape that any hydrogen absorption would follow. Balmer absorption can arise from hydrogen atoms in either of the two n=2 states (${}^{2}S$ or ${}^{2}P$), so to calculate the correct column density the relative populations of atoms in each state must be known. Using the $N_{{\rm H}}$ values obtained from the first run, cloudy was set to integrate through to these values and then run for a second time to determine the hydrogen n=2 population values.

The minimum H i n=2 column density at which a measurement would become significant is estimated by replacing the f$\lambda$ value in Equation 5 by a weighted mean value which depends upon the proportion of atoms in each of the two hydrogen n=2 states and using the same $C(v)$ profile as was calculated for $2^{3}S$ He i. The depth at every velocity bin of this profile is then scaled by factors of between 0.001 and 1 in steps of 0.001, with the column density calculated at each step using Equation 5. The first value of the column density at which the 1$\sigma$ uncertainty in the calculation is exceeded is assumed to be the maximum column density of hydrogen in the n=2 state. The constraints on the allowed coarse gridpoints obtained from the metastable Helium, n=2 singlet-state helium and maximum Balmer absorption are illustrated in Fig. 12

A second set of cloudy simulations using the same two-run method described in §4.2 and §4.3 was undertaken to further constrain the absorber properties at a finer resolution. A new grid measuring 24$\times$19 spanning the ranges $-$1.9$\leq$log$U$$\leq$0.4 and 7.1$\leq$log($n_{{\rm H}}$/cm^-3)$\leq$8.9, both in step sizes of 0.1 dex, was generated. Given the minimum value of $N_{{\rm limHe}}(2^{1}S)$=1.76$\times$10¹² cm^-2 listed in Table 3, all log$n_{{\rm H}}$ values in the find grid are allowed. No attempt is made to rule out gridpoints in log$U$ space where the minimum singlet helium column density is not reached as the $2^{1}S$ He i column density is a weak function of log$U$ at a given log$n_{{\rm H}}$ in the grid. Further use is made of the fine grid in the Discussion section (§5).

The presence of turbulence or differential velocity in the outflowing gas can reduce the Ly$\alpha$ optical depth and therefore lead to a higher gas density limit than that derived in these cloudy models. This is due to reduced pumping into the H i n=2 state. Given the FWHM of the absorber ($\sim$600 km s^-1) and the spectral resolution ($\sim$400 km s^-1) the turbulence or differential velocity of the line cannot be larger than $\sim$500 km s^-1. However it is likely to be lower than this, given the saturated nature of the near-IR line (Arav, 2004). Figure 22 (right panel) in Appendix C of Leighly et al. (2011) shows that, for a $2^{3}S$ He i column density of the same order-of-magnitude, the upper limit on log($n_{{\rm H}}$/cm^-3) increases by about 1.5 dex at a turbulent velocity $v_{turb}$=500 km s^-1. At $v_{turb}$$<$500 km s^-1, the upper limit on log$n_{H}$ rapidly declines with decreasing $v_{turb}$.

In §4, the calculations were performed exclusively for epoch 5 due to the relatively narrow constraints obtained on the column density for that absorber, made possible by the detection of metastable helium absorption in both the optical (3889 Å) and near-IR (10 830 Å) wavebands. All other epochs only detected a near-IR absorber, resulting in the calculation of lower column density limits only. Although absorption variability cannot be detected in the optical waveband, there is significant variability in the near-IR profile as indicated in Table 3. There are two main possibilities for this variability mechanism, either (a) changes in the ionization state of the gas under the influence of a varying incident ionizing continuum, or (b) the fraction of the emission source covered by the absorber changes over time due to the absorber’s motion across the line-of-sight. These possibilities are now examined in detail. We only consider velocity space occupied by the epoch 5 absorber due to the availability of the simulation results and the lack of constraints on the (possibly physically separate) absorber at higher outflow (more negative) velocities.

Many AGN are known to have variable continuum emission in the optical and UV wavebands. This is true of NGC 4151 which shows variability across X-ray, UV and optical wavebands on timescales ranging from hours/days to multiple years e.g. Edelson et al. (1996); Crenshaw et al. (1996); Kong et al. (2006). Ionization of gas in the AGN vicinity is mainly controlled by the incident extreme ultraviolet (EUV) flux originating from the continuum emission region. If the gas is sufficiently dense, resulting in short recombination timescales, changes in the strength of associated emission/absorption lines should be correlated with changes in the unobserved EUV continuum. It can be assumed that the optical-UV power-law continuum extends into the EUV, meaning measurements of the behaviour of the optical continuum can be used as a proxy for (ionizing) EUV variations (Kriss et al., 1999). This is the basis of reverberation mapping, where it is known that changes in the broad emission lines correlate with changes in the optical continuum, e.g. Peterson et al. (1998). However, this correlation may not be true for all AGN absorption lines. Studies showing evidence for ionizing continuum driven variability in AGN absorption lines include Wildy et al. (2015); Wang et al. (2015).

To test the plausibility of ionization-driven variability, we measured the relative changes in the flux of the power-law continuum model. The photometric accuracy of the spectra is not sufficient to directly compare continuum strengths across the observational epochs. Instead we scale the flux using the forbidden narrow line [Ne iii] $\lambda$3868, which is not expected to vary between epochs due to the spatially extended nature of the NLR and the low responsivity of forbidden lines to ionizing flux changes. The following procedure is therefore undertaken. At each epoch we multiply the power-law model by the flux of the [Ne iii] $\lambda$3868 template before dividing by the flux of the epoch 6 template and finally dividing the result by the power-law model at epoch 6 (template construction was described in §3.2.3). As noted, no template was generated for epoch 1, so instead at this epoch we use the same flux ratio but over only the unblended (blue side) of the line after subtracting a local linear continuum. The results are flux-corrected continuum models normalised to epoch 6. The values recorded at 4000 Å for each epoch are displayed in Table 4. The continuum peak–trough periodicity is consistent with V-band photometry of NGC 4151 recorded in Guo et al. (2014), which overlaps our epochs 3–5.

If continuum changes are driving the variability, it is important to understand in what sense the absorber should respond. For the epoch 5 absorber, Fig. 13 shows a representative curve of absorber column density varying as a function of log$U$ taken at example physical conditions allowed by the cloudy fine-grid simulation (log($n_{{\rm H}}$/cm^-3)=8.2 cm^-3, log($N_{{\rm H}}$/cm^-2)=21.5). Since log$n_{{\rm H}}$ is held constant, variations in log$U$ must correspond to changes in ionizing photon flux. From epoch 1, the continuum drops from a high state to a low trough at epochs 3 and 4, before reaching to its highest value of all 6 observations at epoch 5. It then returns to a lower state at epoch 6. Given the shape of the column density response in Fig. 13, this pattern of continuum variability would have to explain the absence of the epoch 5 absorber at epochs 1 and 2 (possibly also absent at epochs 3 and 4) and its subsequent disappearance at epoch 6. This is consistent with a situation where log$U$ is initially sufficiently low that any $2^{3}S$ He i absorption is undetectable (left of the peak in Fig. 13). The subsequent highly significant increase in continuum flux at epoch 5 boosts the metastable hydrogen column density to within the limits allowed by the partial coverage model, with the situation reversing between epochs 5 and 6.

From the example in Fig. 13, a reduction in log$U$ of at least $\sim$0.6 dex is required to remove the epoch 5 absorber from detectability at epochs 1 and 6. Given the fact that variability at ionizing wavebands can be many times greater than that in the optical e.g. Edelson et al. (1996), this is achievable, as the 4000 Å flux ratios at these epochs relative to epoch 5 are $\sim$0.5 and $\sim$0.1 respectively. The observations therefore support a situation where changes in the ionizing continuum incident on the absorber generate its measured variability.

The plausibility of changes in gas covering fraction as an explanation for the observed variability can be tested using a method similar to the “crossing disc” model applied in Capellupo et al. (2013); Capellupo et al. (2014). In this model a disc-shaped absorber moves across a disc-shaped emission source, changing the fraction of the emission covered and hence the apparent strength of the absorber. Given the values for the near-IR absorber recorded in Table 3, a plausible scenario that occurred in NGC 4151 over the course of the observations is that of a cloud starting to eclipse the AGN emission region as viewed from Earth at some point between epochs 2 and 3, reaching a peak coverage near epoch 5 before moving out of the line-of-sight prior to the epoch 6 observation. For the purposes of this study the model relates the velocity of gas across the line-of-sight to the crossing time as follows

Covering fraction changes ($\Delta$A) are computed across time intervals beginning or ending with epoch 5. Epoch 1 is excluded as the sense of any change between epochs 1 and 2 is completely unknown. Therefore, calculations for the intervals between epochs 2–5, 3–5, 4–5 and 5–6 are performed. Checking the consistency of the crossing-clouds model requires calculating the value of $\sqrt{\Delta{}A}$/$\Delta$t and its associated error at each epoch interval used. These are listed in Table 6.

The velocity values calculated at the two epoch intervals showing significant absorption at both endpoints (3–5 and 4–5) are 420 km s^-1 and 570 km s^-1 respectively. The best value is therefore taken as the mean of these ($v_{c}$=495 km s^-1) under the assumption that crossing-clouds model is correct. One possible explanation for movement of absorbing gas across the line-of-sight to an AGN is that the gas is in orbit around the black hole with a lateral Keplerian velocity. The radial distance in such a scenario is $r_{{\rm kep}}$=GM/$v_{c}$². Using the estimated mass of the black hole in NGC 4151 of $M_{{\rm BH}}$=4.57$\times$10⁷ $M_{\odot}$ from Bentz et al. (2006), a value of $r_{{\rm kep}}$=2.5$\times$10¹⁸ cm (0.81 pc) is obtained.

If a crossing-clouds model were correct, consistency would be expected between values of $\sqrt{\Delta{}A}$/$\Delta$t recorded at intervals 3–5 and 4–5. Additionally, intervals 2–5 and 5–6 should each show values consistent with or less than those at intervals 3–5 and 4–5, since the time the absorber spends outside the line-of-sight between epochs 2 and 3 and epochs 5 and 6 is unknown. From Table 6 it is clear that these conditions are met. Although this leads to self-consistency within the crossing clouds scenario, it does not provide conclusive evidence of its reality. Such evidence would require: (a) tightly constrained non-zero covering fractions across several epochs, and (b) knowledge that the same absorber is likely recorded at each epoch. It also seems coincidental within this model that the peak covering fraction should occur at maximum continuum strength. Therefore ionization change driven by the continuum is favoured as the variability mechanism.

The un-normalised SED used in the cloudy calculations can be appropriately scaled using the value of the bolometric luminosity given by $L_{{\rm bol}}$=7.3$\times$10⁴³ erg s^-1 (Kaspi et al., 2005). Using this scaled SED, we can estimate the distance from the ionizing continuum source to the incident surface of the epoch 5 absorber at each of the fine gridpoints described in §4.4, as shown in Fig. 14. Radial distances are within 0.49 pc of the continuum source for 90 per cent of all valid gridpoints, with an inner limit of 0.03 pc.

At a given log$U$ gridpoint, total hydrogen column densities are a very weak function of log$n_{{\rm H}}$. Therefore only one set of points corresponding to log($n_{{\rm H}}$/cm^-3)=7.9, representing log$N_{{\rm H}}$ as a function of log$U$, is indicated in Fig. 15. Allowed limits are 21.2$\leq$log($N_{{\rm H}}$/cm^-2)$\leq$23.3. Hydrogen column densities are required to calculate mass outflow rates, can be estimated from total hydrogen column densities using the method of Crenshaw et al. (2009) as follows

Since cases of absorption from metastable triplet state He i are effectively high ionization transitions, the observed absorption may also be visible in high ionization UV transitions resulting from species such as C iv and P v. Although no simultaneous UV spectra are available, the centroid of the epoch 5 absorber is similar to the ’D+E’ feature described in Kraemer et al. (2005, 2006) and Crenshaw & Kraemer (2007). The values listed in Crenshaw & Kraemer (2012) for the ionization parameter (log$U$=$-$1.08), radial distance ($\sim$0.1 pc) and mass outflow rate (0.014 $M_{\odot}$ yr^-1) found for the UV sub-component labelled ’D+Ec’ lie within the range allowed in Fig. 14. Therefore the $2^{3}S$ state-helium absorber studied in this paper is consistent with the physical properties of ’D+E’ without the high cloumn density components labelled ’D+Ea’ and ’Xhigh’ contributing. The lack of high column absorbers in our study may not be totally surprising given the lack of detection of H$\beta$ absorption, in contrast to that found in the June 2000 HST observation reported in Crenshaw & Kraemer (2007). It is likely that our epoch 5 observations detect an outflow between the inner narrow line region and the BLR, possibly within an intermediate line region where self-absorption occurs (Crenshaw & Kraemer, 2007). A hydrogen number density range of 7.1$\leq$log($n_{{\rm H}}$/cm^-3)$\leq$8.8 is intermediate to the BLR and NLR densities and is therefore consistent with the location found in the aforementioned papers. If this UV absorption results from the same obscurer as our epoch 5 observation, it rules out the crossing cloud scenario since, at the velocities listed in §5.3, the absorber would have completely crossed the emission region at least $\sim$3 years before epoch 5.

Given the range of possible radial distances (see Fig. 14), the $2^{3}S$ helium epoch 5 absorber likely lies at or beyond the radial distance of the inner edge of the dusty torus. The torus is thought to surround the central regions of AGN including NGC 4151, with a structure extending from the inner face out to a few pc (Riffel et al., 2009). In Koshida et al. (2014) the observed time lag between the V-band and K-band light-curve behaviour in NGC 4151 suggested an inner radius of the dust torus of approximately 30 to 80 lt-days. The earlier study of Minezaki et al. (2004) was consistent with this finding, having calculated a hot inner edge of the torus in NGC 4151 located at $\approx$0.04 pc from the central ionising source. It is therefore possible that the outflow consists of material originating in the torus, as has been proposed for a broad absorption line in Leighly et al. (2015).