Testing the Lamp-post and Wind Reverberation Models with XMM-Newton Observations of NGC 5506

The long-term variations in the spectrum of NGC 5506 are characterized by a roughly constant shape that varies in flux. This is illustrated in Figure 1, which shows the spectra from all nine EPIC-PN exposures; the effective area curve has been factored out by unfolding the spectra to a constant. The plot shows that the shape is roughly constant over the years while the flux changes by ∼30% above 2 keV. As pointed out by Guainazzi et al. (2010), the absorption does not appear to change with time as indicated by the same factor change at 2 and 10 keV, suggesting that the absorption is distant and likely produced in the large dust lane seen across the optical image of the galaxy (Baillard et al. 2011). The soft part of the spectrum (<2 keV), which is dominated by photoionized gas emission originating at scales of hundred parsecs from the nucleus, shows very little variability. The lower panel of Figure 1 shows the long-term fractional variability amplitude (e.g., Vaughan et al. 2003a) measured between observations. It has a constant shape between 2 and 10 keV implying a constant spectral shape. It drops below 1 keV because the spectrum is not variable in that band. It also shows a dip at 6.4 keV, corresponding to a less variable narrow Fe Kα line, with less clear structure between 6 and 7 keV.

Zoom In Zoom Out Reset image size

We focus the spectral modeling on the nuclear emission dominating the 2–10 keV band. The spectrum below 2 keV is not directly emitted in the nuclear region (Bianchi et al. 2003). A basic model would include an absorbed power law to model the primary coronal emission (modeled with ztbabs∗powerlaw), and a model for the strong narrow Fe Kα emission line. The latter component can originate in Compton-thin material such as the broad or narrow line regions (BLR or NLR), or Compton-thick material like the outermost regions of the disk or a torus. The two can be distinguished by the Compton reflection (CR). Matt et al. (2001) used the edge of the CR to imply that the line-emitting region is Compton-thick with optical depth larger than unity. We model this line using the xillver model (García et al. 2014). The starting base model therefore has the form tbabs∗(ztbabs∗powerlaw + xillver), where the tbabs model accounts for the Galactic absorption column density, fixed at 4.1 × 10²⁰ cm⁻². Here, we make the inconsequential assumption that the reflection is outside the local absorber. Including the reflection inside the absorber only changes the inferred intrinsic flux of the narrow Fe Kα line. Fitting the nine spectra separately, we find that this model accounts for the general shape of the spectrum, but leaves strong residuals peaking at 6.7 keV.

Given the detector resolution, it is not clear whether these residuals are due to a single broadened ionized disk emission, or a blend of narrow emission lines (H-like and He-like iron). These latter lines can be produced by fluorescence and resonant scattering in photoionized matter (Matt et al. 1996), and likely produced in the region producing the spectrum below 1 keV.

To quantify the significance of the residuals in each observation, we show them in Figure 2 as contours of significance levels in the energy–intensity space, produced by adding a narrow Gaussian function to the base model with a grid of energy and intensity values. The significance was calculated following the prescription in Zoghbi et al. (2015), where we find the Δχ² improvement at grid point, we then fake 10,000 spectra assuming model parameters from the base model and their uncertainties, and find the distribution of Δχ² that are due to the observational statistics. These Δχ² values are considered at any observed energy, so the number of trials is properly accounted for. The observed Δχ² is then compared to this simulation distribution to obtain a significance. The nine spectra are assessed independently.

Zoom In Zoom Out Reset image size

Figure 2 shows that the 6.7 keV feature is detected at more than 3σ confidence in 5 out of 9 observations. The weaker detection in the other four is likely because of their low exposure and low count rate. The apparent absorption line at 6.4 keV in the last three observations is an artifact of the base model trying to model the strong unmodeled residuals at 6.7 keV.

The residuals can be modeled as a sum of two narrow unresolved Gaussian functions at 6.7 and 6.9 keV (i.e., their width is consistent with the detector resolution) or a single broad line. The narrow lines can be due to recombination lines from Fe xxv and Fe xxvi respectively (Guainazzi et al. 2010). When fitted with two narrow lines, the 6.9 keV line is consistent with having a constant flux, while the constant flux hypothesis for the 6.7 keV line is ruled out at more than 99.99% confidence (this can also be seen from Figure 2 by comparing observations 8 and 9 for instance).

There are indications in the data that the total flux from these residuals (the sum of two narrow lines or a single broad line) is correlated with the total continuum flux. The Spearman's rank correlation coefficients are 0.57 and 0.65 when using the fluxes from the narrow lines and the broad line compared to the 7–10 keV flux, with the no-correlation hypothesis rejected at the 95% and 98% confidence, respectively.

Distinguishing between narrow line versus broad line models is not possible based on spectroscopy alone, and this is the origin of the discrepancy in interpreting the spectra of NGC 5506 in the literature (e.g., Bianchi et al. 2003; Guainazzi et al. 2010; Sun et al. 2018). The suggestion that the line flux is correlated with the continuum flux may be an indication that the broad line origin is more likely. In Section 4, we include additional information from the fast variability (see the light curve in Figure 3). For the completeness of the spectral analysis, we also model the spectra using a full relativistic model and present the model parameters next.

3.2.1. Relativistic Model

We start by estimating the power spectrum density (PSD) using the periodogram following the standard procedures (e.g., Vaughan et al. 2003a). Using all the observations, we find that the periodogram strongly favors models that include a break around 0.1 mHz. Fitting a bending power-law model with the lower index fixed at 1 gives a break frequency of 0.19 ± 0.04 mHz, and a high frequency slope of a₂ = 3.2 ± 0.4, both consistent with published work (González-Martín & Vaughan 2012, using observations 7 and 8 only). Fitting the same model to the high and low flux intervals separately gives parameters that are consistent with those found using all the data, therefore, significant nonstationary behavior is ruled out.

We also calculate the covariance spectra, which measure the shape of the correlated variability at some Fourier frequency (e.g., Uttley et al. 2014). We use two frequency bins 0.008–0.1 and 0.1–0.26 mHz and refer to them as the low and high frequency bins. The justification of the bin choice is based on the lag measurements discussed in Section 4.3. We find that the covariance spectra can be described by an absorbed power-law model for both frequency bins. The fit parameters provide evidence that the continuum is significantly harder at higher frequencies: Γ = 1.91 ± 0.03 versus 1.75 ± 0.05 for the low and high frequency bins respectively. This is consistent with observations of other sources that show flatter PSDs at higher energies (McHardy et al. 2007; Zoghbi et al. 2011; González-Martín & Vaughan 2012), and it is an indication that the variability is driven by the primary continuum, not by other processes such as intervening absorption.

Zoom In Zoom Out Reset image size

Adding a Gaussian line to the covariance spectra around 6.7 keV did not provide significant fit improvement. This applies to the whole data set and also when considering the high and low flux data separately. The 90% confidence upper limits on the reflection fraction⁵ at the line peak (assuming the energy and width are similar to the time-averaged spectra) are (6%, 14%) for the low and high frequency bins respectively. For comparison, the fraction of the line contribution, also at the peak, from the time-averaged spectrum is 5.4 ± 0.4%. The main conclusion here is that the covariance spectrum provides only upper limits on the fractional variability of the 6.7 keV feature, and these limits are higher than its fractional flux contribution to the time-averaged spectrum. Assuming that the line varies only as a response to continuum variability (i.e., it cannot vary more than the continuum), the reflection fraction in the time-averaged spectrum itself can be used as an upper limit on the variable reflection fraction.

Time lags generally depend on both frequency and energy. For the frequency dependence, light curves at three energy bands, 2–6, 6–7, and 7–10 keV, are extracted. The first and the last bins are generally dominated by the continuum, and the second bin is dominated here by the emission from the 6.7 keV feature. The light curves are tapered with a Hanning function (Bendat & Piersol 2000) to reduce red noise leak effects before calculating time lags following the standard procedure (Uttley et al. 2014). The frequency axis was binned with a geometric factor of 1.7, requiring that each frequency bin contains at least 15 independent Fourier frequencies. The result is shown in Figure 5. Positive lags indicate that the hard band lags the softer band. Figure 5 shows that below ∼0.3 mHz, where the signal is not dominated by Poisson noise, the iron band, 6–7 keV, appears to be delayed with respect to the 2–6 keV band between 0.1 and 0.2 mHz, while the 7–10 keV is the one that lags the 2–6 keV at lower frequencies.

Zoom In Zoom Out Reset image size

Given these indications, we calculate the lag-energy dependence at two frequency bins: 0.008–0.1 and 0.1–0.26 mHz. The orange points in Figure 6 show the results of calculating the lags using eight energy bins between 2 and 10 keV in log-steps. The top panels show the result when using all observations while the middle and lower panels are for the high and low flux intervals respectively. The shaded green regions show the lags produced when using 16 energy bins, and explores the effect of the energy binning choice.

Zoom In Zoom Out Reset image size

The plots in Figure 6 show that there are clear deviations from zero lag. The significance of these measurements against a constant and a linear (versus log-energy) null models are shown in Table 2. The constant model tests for the presence of any inter-band delays. The log-linear model tests the data against a featureless trend.

Table 2.  Results of the Hypothesis Test for the Lag Spectra in Figure 6

Obs	constant		linear
	Bin-1	Bin-2	Bin-1	bin-2
All obs.	${\bf{3}}({{\bf{10}}}^{-8})$	${\bf{2}}({{\bf{10}}}^{-4})$	0.40	${\bf{4}}({{\bf{10}}}^{-3})$
High flux	${\bf{7}}({{\bf{10}}}^{-10})$	${\bf{1}}({{\bf{10}}}^{-6})$	0.16	${\bf{5}}({{\bf{10}}}^{-6})$
Low flux	${\bf{1}}({{\bf{10}}}^{-4})$	${\bf{6}}({{\bf{10}}}^{-3})$	$1.8({10}^{-2})$	$5.7({10}^{-2})$

Note. The table shows the null hypothesis p value from each test for a constant and a linear models. Values in bold are those below p = 0.01, corresponding to a 99% confidence in rejecting the null hypothesis. The tests are done using the 16 bin data. Other binning schemes give comparable results.

Download table as:  ASCIITypeset image

Table 2 shows that the constant model is rejected at more than 99% confidence in all cases. The log-linear model is only rejected in the high frequency bin of the combined data set and the high flux subset. The featureless lag at the low frequency bin appears to be similar to the continuum lags that often dominate at low frequencies in many sources (Miyamoto & Kitamoto 1989; Kotov et al. 2001; Arévalo & Uttley 2006; Zoghbi et al. 2011).

The only case therefore where complexity beyond the two simple models is required is the high frequency bin of the high flux observations (in addition to the combined observations, whose lags are driven by the high flux data). Modeling the additional complexity with a Gaussian function gives a centroid energy of E = 6.5 ± 0.1 keV. This energy suggests that these lags are related to the 6.7 keV spectral feature found in the time-averaged spectra discussed in Section 3.2.

Interpreting the delay values requires that the observed inter-band lags are converted to component delays by accounting for their relative contribution to the variable spectrum at different energies (i.e., the reflection fraction). Before using a lamp-post model to fit the lag data, we first consider a model-independent dilution correction.

Based on the discussion in Section 4.2, we use the upper limit on the reflection fraction from the time-averaged spectrum. We use the best-fitting models for the time-averaged spectra to calculate the 99.7% confidence upper limit on the fractional contribution of the reflection component at each energy band f_r(E). Two cases are considered for modeling the reflection line: a Gaussian function, and a full relxill model (see Section 3.2). Then a model of the form a + b × f_r(E) is used to fit the lag spectra, where a and b are model parameters. a is a shift parameter that accounts for the reference band in the lag calculations, and b measures the delay between the reflection and direct components. This modeling assumes there are no inter-band delays within the continuum and the reflection spectra themselves.

Fitting this model to the high frequency high flux lag data, we obtain lower limits (99.7% confidence) on the dilution-corrected lag of b > 2755 s when using a Gaussian function, and b > 2106 s when using the full relxill model. The limit is smaller when using the full reflection model because reflection contributes to the whole reference band, not just where the line is present, and that dilutes the delays. Because we are discussing lower limits on the lag (given the upper limits on the reflection fraction), we consider the smaller value from the relxill model, giving a 99.7% confidence lower limit on the dilution-corrected lag of τ > 2106 s.

Now we consider the frequency at which the lags are measured. A time lag τ cannot be measured at frequencies higher than the wrap frequency ν_w = 1/2τ, because at ν_w, a lag of τ is the same as a lead of −τ, as the phase delay ϕ = 2πν_wτ = π, and a phase lag of π is the same as a lead of −π. In fact, this is the case for all frequencies ν_n = n/2τ where n ≥ 1 is an integer. Also, τ here is the dilution-corrected value, not the observed inter-band lag.

Given the lower limit on the dilution-corrected lag of τ > 2106 s, the first (n = 1) wrap frequency is ν_w < 0.24 mHz. This is comparable to the frequencies in which the lag is observed (0.1–0.26 mHz), implying that direct interpretation of the lags requires careful considerations.

First, we consider the possibility that the lag we measure is above the wrap frequency. This would suggest that the observed lag is in the oscillatory part of the delay transfer function, and the actual intrinsic lag is much larger than the observed lags. This would also mean that the lag of the iron band relative to the continuum instead of a lead is just a coincidence. A prediction of this possibility is that the lag should flip sign when we probe lower frequencies, the iron band will be leading instead of lagging due to the oscillatory nature of the transfer function.

This prediction can be tested given that we have data in one lower frequency bin, namely 0.008–0.1 mHz, so we can test for the presence of the 6.7 keV feature and for the change of the lag sign. Unlike the high frequency lag data, the low frequency has a trend that increases with energy. Accounting for the trend, we find that the presence of a feature at 6.7 keV in the low frequency band that is similar to the one at high frequencies cannot be ruled out. The contribution in the low frequency bin is constrained to be a factor f = 0.8 ± 0.5 of the lag at the high frequency bin. In other words, we do not see any significant change in the lag of the 6.7 keV feature in the low frequency bin, and a flip in the sign (f = −1) can be ruled out at more than 99.9% confidence.

This implies that the lag of the 6.7 keV feature is present at least down to 0.008 mHz, or a factor of ∼30 in frequency (see also Figure 5), implying that we are not observing the lag in the oscillatory part of the transfer function, at least for simple reverberation models (see a similar discussion in Zoghbi et al. 2011), because the lag is expected to switch sign at multiples of ν_w, and the smaller its value, the more sign changes should be observed. The question of whether there is a sign flip at frequencies lower than 0.008 mH remains open, and future observations that target timescales >100 ks will be able to address it.

In this section we attempt to model the iron band lag, taking the information from the time-averaged spectrum into consideration. In Section 3.2.1, we showed that the 6.7 feature in the time-averaged spectrum can be modeled with weakly relativistic disk reflection in a lamp-post geometry. We use the reverberation model kynxilrev,⁶ which computes the time dependent reflection spectra of the disk as a response to a flash of primary power-law radiation from a point source located on the axis of the black hole accretion disk (Dovčiak et al. 2004).

We start by fixing the model parameters to those found in the spectral modeling in Section 3.2.1. The only free parameter is the black hole mass (M₈ in units of 10⁸ M_⊙), which converts the scale from r_g to physical units. The black hole mass in NGC 5506 is uncertain, with stellar dispersion measurements suggesting M₈ = 0.88 (Papadakis 2004), and X-ray variability timescale suggesting M₈ = 0.03 (Nikołajuk et al. 2009). We model the high frequency high flux lag spectra, and the best-fit model is shown as a solid line in the right panel of Figure 7. We find χ² = 45.6 for 15 dof, and ${M}_{8}=({6}_{-2}^{+52})\times {10}^{-3}$, closer to the X-ray variability mass estimate. This model is rejected by the lag data with high confidence (p = 6 × 10⁻⁵). Although it is not unreasonable to think that increasing M₈ increases the lags for the same reflection fraction, doing so also reduces the frequency over which those lags are observable, and that is what puts an upper limit on the mass. Fundamentally, the reflection fraction in the model that fits the time-averaged spectra is too small to account for the lags.

Zoom In Zoom Out Reset image size

To improve the fit to the lag data, we tested allowing several parameters to change. The most significant changes were found by allowing the inner radius R_in (Δχ² = 7.9 for 1 dof) and inclination θ (Δχ² = 15.1 for 1 dof) to change. The best-fit model that allows for both these parameters to change is shown as the dotted line in the right panel of Figure 7. The parameters of the model are ${R}_{\mathrm{in}}={4}_{-3}^{+22}{r}_{g}$, $\theta ={({15}_{-15}^{+9})}^{^\circ }$, and M₈ = 0.013 ± 0.003. This model describes the lag data very well (χ² = 19.9 for 13 dof, p = 0.1).

This model, as Figure 7 shows, has a higher reflection fraction, so larger lags can be produced for a given black hole mass, allowing it to model the lags properly. The increase in reflection fraction, however, now overpredicts the strength of the line in the time-averaged spectrum, which is better described by a lower reflection fraction (left panel of Figure 7).

It is also possible to model the lags with a model where the ratio of the primary to the reflected normalization are not restricted to the static lamp-post case. As expected, the fit improves significantly (Δχ² = 19 for 1 less degree of freedom), and the primary to reflected ratio is ${3.9}_{-0.4}^{+2.6} \%$ (parameter Np:Nr in kynxilrev). In other words, a significant fraction of the photons from the illuminating source need to hit the disk rather than reach the observer, which could be an indication of beaming. Such a model, however, also overpredicts the strength of reflection in the time-averaged spectrum.

In summary, the lag spectra requires the reflection fraction to be larger than that observed in the time-averaged spectrum. Because of this, we find that a static lamp-post model cannot reproduce both the lag and time-averaged spectra at the same time.

NGC 5506 shows a strong, resolved residual feature at 6.7 keV that has been attributed to either a blend of at least two narrow emission lines from Compton-thin plasma, or to a weakly relativistic disk line. The spectral modeling of a new observation, along with old observations, cannot distinguish between the models, primarily because of the limited energy resolution of the detectors. The data suggest that the flux from the line(s) is correlated with the primary continuum flux on timescales of years (Section 3.2.1), indicating that the line is responding to changes in the continuum.

On timescales of ∼10 ks, the data allows only for an upper limit on the fractional variability of the 6.7 keV feature (Section 4.2), which is consistent with its fractional flux contribution, an indication that the feature is likely varying with the continuum at these timescales.

The lag analysis in Section 4.3 shows that the 6.7 keV feature is delayed with respect to the continuum, and that the delay depends on the flux level. The simplest interpretation, again, is that the 6.7 spectral feature is due to a weakly relativistic disk line that is delayed with respect to the continuum.

Using a disk reverberation model that assume a static lamp-post geometry (Section 4.5) with parameters from the time-averaged spectrum fails to fit the lag spectra, underpredicting the lags. A model with a lower disk inner radius and inclination has a higher reflection fraction, and is able to reproduce the lags, but now overpredicts the contribution of reflection to the time-averaged spectrum (Figure 7). In other words, the lag and energy spectra, taken together, are inconsistent with a lamp-post geometry where a static point source illuminates the accretion disk. The main reason is that the reflection fraction in the time-averaged spectra is too small to account for the lags.

The failure of the lamp-post reverberation model results is similar to several other studies that attempted to model both the spectral and timing data at the same time. This was manifested for instance as wavy residuals in the lag-frequency data Caballero-García et al. (2018) or a 3 keV dip in the lag-energy spectra (Kara et al. 2013a; Chainakun et al. 2016; Wilkins et al. 2016; Ingram et al. 2019). We note that the 3 keV dip may not be a dip, but a consequence of the fact that the reflection fraction required to model the Fe K lags is higher than that inferred from the time-averaged spectrum (e.g., Figure 10 in Chainakun et al. 2016 and Figure 16 in Ingram et al. 2019), similar to what we find here in NGC 5506. We note also that the statistical significance of the 3 keV dip has not always been addressed explicitly, so its existence as a separate feature may not be significant in all cases discussed in the literature (e.g., Wilkins et al. 2017).

The main observable discrepancy between the lag and time-averaged spectra when modeled with a static lamp-post model is the enhanced reflection fraction in the former. The case of NGC 5506 presents a clear manifestation of this discrepancy.

It is conceivable that the lag spectra, through the frequency filtering, selects regions with enhanced reflection, whose effect is averaged out in the time-average spectra. Sheared hot spots in the disk may produce such an effect. However, these are expected to be random, transient, and likely present throughout the disk, and it seems unlikely that they produce a consistent picture in which reflection is enhanced when selecting some timescales and not others.

Modeling the lag data, therefore, clearly requires models beyond a static lamp-post geometry. Our modeling of the NGC 5506 suggests that the reflection fraction in the lag needs to be enhanced by a factor of 6 compared to that inferred from modeling the time-averaged spectrum.

A few modifications to the simple lamp-post geometry we used here have been considered in other works. Chainakun et al. (2016) model the additional complexity in the data by including an ionization gradient in the disk. Unlike the model we used here, where the whole disk is assumed to have the same ionization parameter, they find that if the innermost regions are highly ionized compared to the outer regions, the reverberation signatures produced through reflection from the inner parts (iron Kα photons redshifted to 3 keV) are less observable. The line is enhanced as it is dominated by reflection from colder reflection off the outer part. We test this model in the data by allowing the density (parameter density) and the density profile (den_prof) in the model to change. These provide no significant improvement in the fit (Δχ² = 0.5 for two degrees of freedom). This is likely a result of the fact that the inner radius of the disk in NGC 5506 is relatively large, and the impassivity index is flatter (q = 3) than the models in Chainakun et al. (2016), so the ionization gradient does not provide enough contrast between the inner and outer disks to provide the necessary reflection fraction enhancement.

Wilkins et al. (2016) considers a more complicated case of a combination of hard lags produced by propagation through an extended emitting region and reverberation from a compact, point-like corona. If luminosity fluctuations propagate upwards through a collimated, vertically extended corona, it is found that the interplay between light-travel delays and propagation delays can produce lag-energy spectra where photons at 3 keV are the earliest to arrive, producing a dip at 3 keV and an enhanced line feature. A slightly different, but comparable, geometry was considered by Chainakun & Young (2017), where two axial point sources illuminating an accretion disk to model the data of PG 1244+026 are used. The lags are produced from the light-travel delays between the two points and the disk and also from the assumption of different source responses for two X-ray sources. Unlike other source, strong relativistic effects are not present in NGC 5506, implying that the discrepancy between the lag, if attributed to relativistic reverberation, and the time-averaged spectra is not an effect of general relativity, but rather likely a geometrical effect.

It is worth noting that complexity beyond the lamp-post may well be required given the turbulent flow in accretion disks around black holes. However, the simplicity and predictive power in relativistic models is lost. Simple fits to the data to extract important parameters such as the black hole spin is no longer a simple task.

Given these implications for the static lamp-post relativistic model, we also consider alternative models to explain the spectral and timing observations in NGC 5506. While Miller et al. (2010a, 2010b) and Turner et al. (2017) showed that lags from a cloud at ∼100 r_g could explain the delays observed below 1 keV in several sources. Mizumoto et al. (2018, 2019) showed that lags in the Fe Kα band, similar to those observed, can be produced by a distant scattering wind that is out-flowing at ∼0.1c. The short lags compared to the light-travel time are explained by dilution effects, where the majority of the photons in the Fe Kα band are primary photons with zero delay, whereas the time-delayed reprocessed photons only make a subtle contribution.

To check if this model can explain both the spectra and the lags in NGC 5506, the observations are compared to spectral and timing products generated using the model in Mizumoto et al. (2018). This model is not designed to fit the data, so we focus here on products that are generated by only at few parameters that were found to roughly match the shape of the spectrum and lag scale. We find that the wind needs to be ionized to explain the 6.7 keV features. We also find that no simple model can model both the spectrum and the lag at the same time for the same reason the lamp-post model fails. The reflection fraction from the model that fits the time-average data is not large enough to explain the lags (see the resulting spectrum in the left panel of Figure 8 and the lag in the orange dotted line in the right panel of Figure 7). If the lags are made larger by increasing the size of the scattering clouds or by increasing the scattering fraction outside the line of sight, the frequencies at which those lags are produced decreases compared to the observations.

Zoom In Zoom Out Reset image size

If absorption in the line of sight from these same scattering clouds is also included in the model, we find that a larger lag for the same reflection fraction in the time-averaged spectra can be produced. This is a result of the reduced dilution in the primary continuum as the flux at the line energies is reduced due to the presence of absorption lines. This model, however, also predicts that several absorption lines other than those at 6.7 and 6.9 should be observed, and that means the model fails to reproduce the spectral data correctly above 7 keV as Figure 8 (right) shows. The spectral data themselves show no evidence for absorption lines as the model predicts.