SPITZER SPACE TELESCOPE MEASUREMENTS OF DUST REVERBERATION LAGS IN THE SEYFERT 1 GALAXY NGC 6418

In the active galactic nucleus (AGN) unification paradigm, direct observation of the nucleus is blocked by a toroidal structure of dusty molecular gas for a range of viewing angles (e.g., Antonucci 1993). As this dust absorbs UV–optical radiation from the accretion disk and re-emits in the infrared (IR), this structure is also thought to be the dominant source of IR radiation in most AGNs. Understanding this obscuration of the central engine is therefore important to understanding the physical processes operating in AGNs and, more generally, their role in galaxy evolution.

The observational evidence (Antonucci 1993; Jaffe et al. 2004; Tristram et al. 2007), indicates that the obscuring structure is geometrically and optically thick, although a warped thin disk that extends throughout the host galaxy has also been proposed (Sanders et al. 1989). The conventional picture is that of a compact, but geometrically thick, torus of optically thick molecular clouds with a size of a few parsecs (Antonucci & Miller 1985; Krolik & Begelman 1988; Pier & Krolik 1992). Models in which the vertical thickness is supported by large random velocities due to elastic collisions between clouds (Krolik & Begelman 1988), or by IR radiation pressure (Pier & Krolik 1992; Krolik 2007), or by turbulence induced by supernovae (Wada & Norman 2002; Schartmann et al. 2009) have been explored. In an alternative class of models, the dusty material is not part of an essentially static torus, but is rather embedded in an outflowing hydromagnetic wind launched from the accretion disk (e.g., Blandford & Payne 1982; Emmering et al. 1992; Bottorff et al. 1997; Elitzur & Shlosman 2006; Dorodnitsyn et al. 2012).

Dust radiative transfer models for the torus broadly reproduce the IR spectral energy distribution (SED) of AGNs. Of necessity, early radiative transfer models assumed smooth density distributions (e.g., Pier & Krolik 1993; Granato & Danese 1994; Efstathiou & Rowan-Robinson 1995), but more recently, models for clumpy dust distributions have been developed (e.g., Nenkova et al. 2002; Dullemond & van Bemmel 2005; Hönig et al. 2006; Nenkova et al. 2008a, 2008b; Schartmann et al. 2008). These "clumpy torus" models are more successful in reproducing certain details of the SED such as, for example, the strength of the 10 μm silicate feature (Nikutta et al. 2009; Nenkova et al. 2008b).

The torus is too small to be directly imaged by any existing single telescope. Some constraints on its size and structure can be inferred from SED-fitting using radiative transfer models (e.g., Nenkova et al. 2008b; Mor et al. 2009; Hönig & Kishimoto 2010; Ramos Almeida et al. 2011; Alonso-Herrero et al. 2011), but there are many theoretical and observational uncertainties which obfuscate the results. Other methods are therefore required, the two most important being reverberation mapping and, for relatively close objects, IR interferometry.

Following the seminal work of Blandford & McKee (1982), the reverberation mapping technique has been well developed and extensively applied to studies of the broad emission line region (BLR). Time series analysis of the response of the broad emission lines to variations in the UV or optical continuum (as proxies for the AGN ionizing continuum) has revealed the characteristic size of the BLR in about 50 AGNs, enabling estimates of black hole masses and Eddington ratios (Peterson 1993, 2006; Gaskell 2009; Galianni & Horne 2013; Du et al. 2014, and references therein). It has also been determined that the BLR follows a size–luminosity relationship of the form $R\propto {{L}^{1/2}}$ (Peterson et al. 2004; Greene et al. 2010; Bentz et al. 2013),

Near-infrared (NIR) (K-band) versus optical (V-band) reverberation lags have been measured for around 20 Seyfert galaxies (Oknyanskij & Horne 2001; Minezaki et al. 2004; Suganuma et al. 2006; Koshida et al. 2009, 2014). As dust grains emitting in the K-band have temperatures close to the sublimation temperature (∼1200–1800 K, depending on grain composition), these lags are thought to represent the inner radius of the torus. The K-band reverberation lags are found to be larger than those of the BLR, while following a similar $R\propto {{L}^{1/2}}$ size–luminosity relation, implying that the BLR is bounded by the dust distribution, consistent with the central idea of the AGN unification scheme.

The inner regions of several bright, nearby Seyfert galaxies have been directly studied using NIR (K-band) interferometry (Swain et al. 2003; Kishimoto et al. 2009, 2011; Pott et al. 2010; Weigelt et al. 2012). The effective ring radii derived from the observed visibilities scale approximately as ${{L}^{1/2}}$, and are comparable with or slightly larger than the radii derived from reverberation lags (Kishimoto et al. 2011). Since Jaffe et al. (2004)'s pioneering study of the archetypal Seyfert 2 galaxy, NGC 1068, mid-infrared (mid-IR) (8–12 μm) interferometric observations have also been obtained for $\approx 20$ AGNs (e.g., Tristram et al. 2007; Burtscher et al. 2009; Kishimoto et al. 2009; Tristram et al. 2009; Hönig et al. 2013) In a recent analysis of the available data, Burtscher et al. (2013) find that while the mid-IR source size scales with luminosity in a manner similar to that seen in the NIR, the inferred size is more than an order of magnitude larger than the measured K-band size and the scatter is quite large.

Here we report initial results from a mid-IR (3.6 and 4.5 μm) reverberation-mapping campaign using the Spitzer Space Telescope in its "warm mission." Our motivation is to probe the dust distribution at spatial scales intermediate between the innermost regions probed by the K-band observations and the outer, cooler regions probed by mid-IR interferometry. Furthermore, variability at 3.6 and 4.5 μm should be less susceptible than the 2.2 μm K-band to complicating effects such as dust sublimation (Minezaki et al. 2004; Kishimoto et al. 2013), or contamination by variable accretion disk emission (Tomita et al. 2006; Kishimoto et al. 2007). During a two-year campaign, we monitored a sample of 12 Seyfert 1 AGNs at cadences of 3 and 30 days during the first and second year, respectively. We selected our targets based on their proximity ($z\lt 0.4$) and their location near one of Spitzer's continuous viewing zones. We obtained B and V images of the targets over the same period using the Liverpool Telescope (LT), the Faulkes Telescope North (FTN) and the Southwestern University (SU) 0.4 m telescope.

In this work we describe our analysis of the first 17 months of measurements of the Seyfert 1 NGC 6418 (Véron-Cetty & Véron 2006), a Hubble-type Sab galaxy (Nair & Abraham 2010) with an apparent magnitude g = 14.87 at a redshift of z = 0.0285 (Ahn et al. 2014). It is classified spectroscopically as a Seyfert 1 on the basis of a strong, broad Hα emission line, but it is otherwise dominated by the stellar continuum (see Remillard et al. 1993, who described it as an "embedded" AGN). Nevertheless, it is also an X-ray source with a 0.1–2.4 keV luminosity of ${{L}_{X}}={{10}^{42.26}}$ erg s⁻¹ (Anderson et al. 2007). We selected NGC 6418 out of our sample due to its larger than average variations in the Spitzer channels for the first year of data; the results of the analysis of the other targets will be presented in a future publication.

We present our observations and describe our methods for measuring the light curves in Section 2. In Section 3 we describe the time series analysis technique that was used to extract the time lags between the 3.6, 4.5 μm and optical light curves. We discuss the implications of our results in Section 4, and present our conclusions in Section 5. Details of our photometric measurements and a comparison of two methods for determining time lags can be found in the appendices.

We will discuss the mid-IR and optical observations separately. See Appendix A for a detailed discussion of our photometric analysis.

2.1. Mid-infrared

2.2. Optical

Contemporaneous optical monitoring was performed in the B and V bands with three ground-based telescopes: Bessel B images were obtained with the 2 m LT on La Palma and the 2 m FTN on Maui; Johnson-Cousins B and V images were obtained with the 0.4 m telescope at SU Fountainwood Observatory (see Table 1). It was not possible to coordinate these observations with each other or with the Spitzer observations, but together they approximately span the time period covered by the Spitzer campaign except during 2011 November, when NGC 6418 was unobservable from the ground. The start and end dates of the observations with each telescope are given in Table 1.

Table 1.  Observations

Telescope	Start Date	End Date	# Obs	Instrument	Filter	Aperture
Liverpool Telescope	2011 Aug 06	2012 Oct 21	64	RATCam	Bessel B	15
Faulkes Telescope North	2011 Aug 10	2012 Sep 30	60	fs02	Bessel B	12
SU Fountainwood 0.4 m	2012 May 19	2012 Dec 18	48	SBIG ST-8300	Johnson-Cousins B/V	35
Spitzer	2011 Aug 01	2013 Jan 04	170	IRAC	ch1/ch2	18

Download table as:  ASCIITypeset image

The exposure times for the optical observations range from 60 to 180 s. Dark/bias subtraction and flat-field division of all images from SU were performed using the XVISTA software package (Treffers & Richmond 1989). Images from the RATCam instrument at the LT were bias subtracted and flat fielded by an automatic pipeline (Steele et al. 2004), as were images taken by the FTN. When more than a single exposure per night was available from LT and FTN, we stacked and registered the images using MATCH, an implementation of the star matching algorithm of Tabur (2007), and the XVISTA package. We then extracted photometry from the stacked image.

Hereafter, we refer to the light curve compiled from the LT and FTN observatories as the combined optical light curve. The SU observations are used to determine the AGN/host ratio. The mean flux densities measured within the aperture in Table 1 for all bands are tabulated in Table 2. These flux densities are not host subtracted.

Table 2.  Mean Flux Density

Name	Mean Flux Density
3.6 μm	3.62 mJy
4.5 μm	3.54 mJy
SU B band	0.53 mJy
SU V band	1.48 mJy
LT B band	0.50 mJy
FTN B band	0.54 mJy

Download table as:  ASCIITypeset image

2.3. Photometry

The photometric analysis proceeds in two stages for the SU dataset: in the first we measure instrumental magnitudes for each object (the target plus comparison stars) in all exposures; in the second the measurements from all exposures in a given passband are combined and the measured instrumental magnitudes are subjected to inhomogeneous ensemble photometry (Honeycutt 1992). For a detailed discussion of these steps see Appendix A. The LT and FTN datasets are reduced using image differencing (Alard 2000 and references therein). The combined optical and Spitzer light curves are shown in flux density, normalized to the mean, in Figures 1 and 2. The light curves are also shown after applying a shift equal to the time lag computed by the cross-correlation analysis (Section 3). In Figures 1 and 2 the time lag shifts applied to the combined optical light curves were 37.2 and 47.1 days, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The optical and infrared curves all show clear variations with similar features on timescales of ∼100 days, but with the variations in the infrared lagging behind those in the optical.

The reverberation lag, τ, between the driving optical continuum variations and those of the responding IR emission gives the characteristic size of the IR emitting region. The lag can be determined by cross-correlating the two light curves. The application of this technique to the broad emission line variability of AGNs ("reverberation mapping") is well developed (Gaskell & Sparke 1986; Gaskell & Peterson 1987; Edelson & Krolik 1988; Maoz & Netzer 1989; Koratkar & Gaskell 1991) and has been widely used to measure the size of the BLR (e.g., Peterson et al. 2004; see Peterson 2001 for a tutorial). As already noted, it has also been applied to optical and K-band light curves in order to determine the inner radius of the torus (Oknyanskij et al. 2006; Suganuma et al. 2006; Koshida et al. 2009, 2014).

We performed cross-correlation analyses for three pairs of data sets: 3.6 μm versus combined optical, 4.5 μm versus combined optical, and 4.5 μm versus 3.6 μm. The time series analysis was performed between the dates of MJD 55900 (2011 May 12) and MJD 56300 (2013 August 1). This time span was selected due to the significant optical and IR variations of the light curves and because there were no large gaps in coverage. For a comprehensive and detailed analysis of individual datasets see Appendix B. For each pair, the cross-correlation function (CCF) was computed using a lag step size of 1 day. The optical observations were not synchronized with the Spitzer observations and are typically separated by irregular intervals. On the other hand, the Spitzer light curves are for the most part more evenly and densely sampled than the optical measurements. Therefore, in order to compute the IR–optical CCFs, we generate IR data points corresponding to the optical observations by interpolating within the Spitzer light curves. For examples of the CCFs computed for the three pairs of light curves see Appendix C.

The maximum of the CCF yields the lag, τ, between the two light curves. However, the maximum is not always well defined, since computed CCFs typically exhibit a broad peak (see Appendix C) and structure in the wings (at large positive or negative lags), which can influence the calculation of the centroid or mean. A common approach is to calculate the centroid of the CCF using a subset of points whose correlation coefficients exceed a certain value; for example, 80% of the peak value (Peterson 2001). Here, we use a different method in which we fit a cubic spline to the CCF and use it to set a threshold for the minimum correlation coefficient. This minimum correlation coefficient is defined as:

$$\begin{eqnarray}\begin{array}{rcl} {\rm C}{{{\rm C}}_{{\rm min} }} & = & {\rm max} ({\rm CCF}(\tau ))-2\sigma ({\rm CC}{{{\rm F}}_{{\rm fit}}}(\tau ) \\ {} & {} & -{\rm CCF}(\tau )) \\ \end{array}\end{eqnarray} \tag{ 1 }$$

where CC_min is the minimum correlation coefficient, $CCF(\tau )$ is the CCF and $\sigma ({\rm CC}{{{\rm F}}_{{\rm fit}}}(\tau )-{\rm CCF}(\tau ))$ is the standard deviation of the difference between the fitted and actual CCF value. The CCF centroid is computed using only values exceeding CC_min. For more details see Appendix C.

To estimate the uncertainty on the CCF lags, we used the cross-correlation centroid distribution (CCCD) method (Gaskell & Peterson 1987; Maoz & Netzer 1989; Peterson et al. 1998), generating 1000 random realizations of the light curves. The CCCDs for the three pairs of data sets are shown in Figure 3 and the derived lags are listed in Table 3 and in Appendix B. The lag is taken to be the median of the distribution and the uncertainty is given by the interquartile range. The CCCDs for the 3.6 μm versus optical, 4.5 μm versus optical and 3.6 μm versus 4.5 μm light curves yield lags of $37.2_{-2.2}^{+2.4}$ days ($31.2_{-1.9}^{+2.0}\times {{10}^{-3}}\ {\rm pc}$), $47.1_{-3.1}^{+3.1}$ days ($39.5_{-2.6}^{+2.6}\times {{10}^{-3}}\ {\rm pc}$), and $13.9_{-3.8}^{+3.7}$ days ($11.7_{-0.1}^{+0.4}\times {{10}^{-3}}\ {\rm pc}$), respectively.

Zoom In Zoom Out Reset image size

Table 3.  Comparison of Cross-correlation Methods

Name	3.6 μm-optical	4.5 μm-optical	3.6–4.5 μm
	(lag(day) ± δ)	(lag(day) ± δ)	(lag(day) ± δ)
Peterson et al.	36.7 ± 3.4	48.6 ± 3.7	14.6 ± 6.0
Zu et al.	$40.4_{-6.5}^{+0.7}$	$49.5_{-4.7}^{+1.2}$	$13.2_{-2.9}^{+5.8}$
Vazquez et al.	$37.2_{-2.2}^{+2.4}$	$47.1_{-3.1}^{+3.1}$	$13.9_{-0.1}^{+0.5}$

Download table as:  ASCIITypeset image

For comparison, we also analyzed our data following the slightly different cross-correlation methods described by Peterson et al. (2004) and Zu et al. (2011). The results are compared in Table 3. We find that all methods yield results which are consistent within the uncertainties for all three pairs of light curves.

The dusty torus absorbs UV/optical radiation from the accretion disk and re-emits it as infrared radiation (Telesco et al. 1984; Sanders et al. 1989). Variability in the accretion disk emission results in corresponding variations in the dust IR emission, but with a delay due to differing light travel times between the source, various points in the torus and the observer. The lags between the optical continuum light curve and the IR light curves can therefore be interpreted as measures of the distance from the source to the dust clouds that predominantly emit the 3.6 and 4.5 μm radiation. Our results indicate these clouds are located at a distance $\approx 1$ light-month ($\approx 0.03\;{\rm pc}$) from the source of the AGN UV–optical continuum. However, the two Spitzer bands have significantly different lags, with the 4.5 μm-optical lag being longer by $9.9\pm \;3.9$ days. The lag between the 4.5 and 3.6 μm light curves is 13.9 ± 0.5 days and is consistent with this difference. This implies that the clouds producing the bulk of the 4.5 μm emission are about 10 light-days (∼27%) further from the UV–optical continuum source.

In most models, the innermost radius of the torus is taken to be the dust sublimation radius which, for a typical ISM dust composition with silicate grains of average size, is (Barvainis 1987; Nenkova et al. 2008b)

$$\begin{eqnarray}&&{{R}_{d,{\rm Si}}}\simeq 1.3{{\left( \frac{{{L}_{{\rm bol}}}}{{{10}^{46}}\ {\rm erg}\;{{{\rm s}}^{-1}}} \right)}^{1/2}}{{\left( \frac{1500\ {\rm K}}{{{T}_{{\rm sub}}}} \right)}^{2.6}}\ {\rm pc}\end{eqnarray} \tag{ 2 }$$

where L_bol is the bolometric luminosity of the AGN and T_sub is the dust sublimation temperature.

However, many broad-line AGNs exhibit a distinct near-infrared "bump," peaking around 2–4 μm, which has a blackbody temperature $T\gtrsim 1000\;{\rm K}$ (e.g., Edelson & Malkan 1986; Barvainis 1987; Rodríguez-Ardila & Mazzalay 2006; Riffel et al. 2009a, 2009b). This feature often dominates the NIR and it has been found that it cannot be reproduced by torus models alone in fits to the infrared SED; instead, one must add a separate hot ($T\sim 1400$ K) blackbody component. The latter has been attributed to hot pure graphite dust located within the torus (Mor et al. 2009; Mor & Trakhtenbrot 2011), and Mor & Netzer (2012) have modeled this component as dust embedded in the outermost BLR, between the sublimation radius for pure-graphite grains,

$$\begin{eqnarray}&&{{R}_{d,{\rm C}}}\simeq 0.5{{\left( \frac{{{L}_{{\rm bol}}}}{{{10}^{46}}\ {\rm erg}\ {{{\rm s}}^{-1}}} \right)}^{1/2}}{{\left( \frac{1800\;{\rm K}}{{{T}_{{\rm sub}}}} \right)}^{2.8}}\ {\rm pc}\end{eqnarray} \tag{ 3 }$$

and the torus inner radius as given by Equation (2). The hot dust spectrum computed by Mor & Netzer (2012)suggests that this hot graphite dust contributes significant luminosity at 3.6 and 4.5 μm.

In order to estimate the sublimation radii given by Equations (2) and (3), it is necessary to determine L_bol. However, this is difficult to determine accurately for NGC 6418, as the optical spectrum is dominated by the stellar continuum and the AGN itself is evidently heavily reddened (Remillard et al. 1993).

The Sloan Digital Sky Survey (SDSS) optical spectrum of NGC 6418 (Ahn et al. 2012) (Figure 4) shows broad Hα and narrow lines of [O iii]$\lambda 5007$, Hα, [N ii]$\lambda 6548,6583$ and [S ii]$\lambda 6717,6731$, but the continuum is dominated by an evolved stellar population. The fact that the broad Hβ line is not evident in the spectrum indicates a steep broad-line Balmer decrement and suggests classification as a Seyfert Type 1.9 (Sy1.9; Osterbrock 1977, 1981). However, NGC 6418 is unusual in that the narrow Hβ emission is also very weak (in fact, this line appears in absorption) and the [O iii]$\lambda 4959,5007$ lines are much weaker relative to the stellar continuum than is typical in Seyferts, even Sy 1.9 s. Interestingly, these lines are not obviously visible in the earlier (1989) spectrum obtained by Remillard et al. (1993), even though the broad Hα line is clearly much stronger relative to the narrow Hα and [N ii] lines than in the SDSS spectrum. Evidently, the strong stellar continuum, the foreground reddening and the variable broad emission lines make the classification of this source somewhat ambiguous.

Zoom In Zoom Out Reset image size

To determine the bolometric luminosity of the AGN, we used the relationship established between the broad Hα luminosity (${{L}_{b{\rm H}\alpha }}$) and the bolometric AGN luminosity (L_AGN) in a large sample of quasars and Sy1 (Richards et al. 2006; Stern & Laor 2012).

$$\begin{eqnarray}&&{{L}_{{\rm bol}}}=130_{\div 2.4}^{\times 2.4}\times {{L}_{b{\rm H}\alpha }}.\end{eqnarray} \tag{ 4 }$$

The flux in the broad Hα line was measured from the SDSS spectrum using Gaussian profiles to fit and deblend the [N ii], [S ii] and Hα lines. In the fit, the wavelengths of the components representing the narrow lines were fixed at the values determined by the SDSS spectroscopic measurement pipeline (spec1d; Bolton et al. 2012). The amplitudes and widths were free parameters, with the exception of [N ii]$\lambda 6548$, which is constrained so as to preserve its fixed 1:3 intensity ratio with [N ii]$\lambda 6583$. The variances provided by the SDSS spectroscopic data reduction pipeline (spec2d; Stoughton et al. 2002) were used to assign weights to each data point; in addition, we assigned a 10% systematic error to the derived fluxes (Bolton et al. 2012). The resulting fit is shown in Figure 5 and the parameters derived from the fit are summarized in Table 4. The broad Hα component has a flux of (2563 ± 120) $\times {{10}^{-17}}\;{\rm erg}\;{{{\rm s}}^{-1}}\;{\rm c}{{{\rm m}}^{-2}}$. Using this Hα flux and assuming a distance of 122 Mpc (Mould et al. 2000), we calculate the observed Hα broad line luminosity to be $L_{{\rm H}\alpha }^{{\rm obs}}$ = (4.56 ± 0.85) ×10⁴⁰ erg ${{{\rm s}}^{-1}}$.

Zoom In Zoom Out Reset image size

Table 4.  Emission Line Fit Parameters

Line	λa	Fluxb	FWHMa
	(Å)	(${\rm erg}\ {{{\rm s}}^{-1}}\;{\rm c}{{{\rm m}}^{-2}}\times {{10}^{-17}}$)	(Å)
${{{\rm H}}_{\alpha \;{\rm broad}}}$	6742 ± 2	2563 ± 120	156 ± 4
${{{\rm H}}_{\alpha \;{\rm narrow}}}$	6753	41 ± 20	5 ± 1
[N ii]	6739	90 ± 25	8 ± 1
[N ii]	6775	269 ± 26	8 ± 1
[S ii]	6912	159 ± 28	10 ± 1
[S ii]	6926	151 ± 27	9 ± 1

Download table as:  ASCIITypeset image

It is clear, however, that a large extinction correction needs to be applied in order to obtain the intrinsic Hα luminosity. From the SDSS spectrum we estimate a lower limit to the broad line Balmer decrement of Hα/Hβ ≥ 6. We used the mean ${\rm H}\alpha /{\rm H}\beta$ from Dong et al. (2005) and their expression to correct for reddening:

$$\begin{eqnarray}\begin{array}{rcl} {\rm log} L_{{\rm H}\alpha }^{\operatorname{int}} & = & {\rm log} L_{{\rm H}\alpha }^{{\rm obs}}+1.87({\rm log} ({\rm H}\alpha /{\rm H}\beta ) \\ {} & {} & -{\rm log} (2.97)) \\ \end{array}\end{eqnarray} \tag{ 5 }$$

which yields a lower limit to the intrinsic broad Hα luminosity of $L_{{\rm H}\alpha }^{\operatorname{int}}\geqslant (1.70\pm 0.32)\times {{10}^{41}}\;{\rm erg}\;\ {{{\rm s}}^{-1}}$.

With this lower limit, Equation (4) yields a lower limit to the bolometric luminosity of the AGN in NGC 6418 of ${{L}_{{\rm bol}}}\geqslant$ ($2.21_{-1.29}^{+3.09}$) × 10⁴³ erg$\ {{{\rm s}}^{-1}}$. Using Kaspi et al. (2000)'s relation ${{L}_{{\rm bol}}}\sim 9\lambda {{L}_{\lambda }}(5100)$ and assuming ${{L}_{\lambda }}(5500)\sim {{L}_{\lambda }}(5100)$ we obtain a lower limit to the AGN V-band luminosity of ${{L}_{V}}\;\geqslant$ ($2.46_{-1.44}^{+3.43}$) × 10⁴² erg$\ {{{\rm s}}^{-1}}$. For comparison, we used the flux variation gradient (FVG) method (Choloniewski 1981; Sakata et al. 2010; Haas et al. 2011, Pozo Nuñez et al. 2012, 2014) to estimate the (constant) host galaxy contribution within our photometric aperture as illustrated in Figure 6. Using the B and V fluxes obtained from the SU observations (35 aperture), we find an AGN/host ratio of 1.55, yielding an an estimate for the AGN contribution to the V-band luminosity of $1.54\pm 0.53\times \;{{10}^{42}}\ {\rm erg}\ {{{\rm s}}^{-1}}$ (as reddening corrections have not been applied to the B and V fluxes, this value should be regarded as a lower limit.) Thus, within the admittedly large uncertainties, the AGN V-band luminosity estimated from the FVG method is consistent with that determined from the ${\rm H}\alpha$ luminosity.

Zoom In Zoom Out Reset image size

Having determined the lower limit on the bolometric luminosity of the AGN, we can determine the dust sublimation radii given by Equations (2) and (3). For silicate dust with a sublimation temperature $\approx 1500$ K, we find ${{R}_{d,{\rm Si}}}\;\geqslant$ $60_{-21}^{+33}\times {{10}^{-3}}\ {\rm pc}$ ($71_{-25}^{+39}$ light days), whereas for pure graphite dust with sublimation temperature ∼1800 K, we find ${{R}_{d,{\rm C}}}\;\geqslant$ $24_{-8}^{+13}\times {{10}^{-3}}\ {\rm pc}$ ($28_{-10}^{+15}$ light days).

These sublimation radii bracket the radii derived from the lags at 3.6 μm (${{R}_{\tau ,3.6}}=31.2_{-1.9}^{+2.0}\times {{10}^{-3}}\ {\rm pc}$) and 4.5 μm (${{R}_{\tau ,4.5}}=39.5_{-2.6}^{+2.6}\times {{10}^{-3}}\ {\rm pc}$). As ${{R}_{d,{\rm Si}}}$ and ${{R}_{d,{\rm C}}}$ are lower limits, this suggests that the bulk of the 3.6 and 4.5 μm emission comes from the region bounded by the graphite and silicate sublimation radii, respectively, and is conceivably emitted by the same graphite dust that is thought to be responsible for the NIR bump. As already noted, the model graphite dust emission spectrum computed by Mor & Netzer (2012), while peaking in the 2–3 μm range, also emits strongly in the 3.6–4.5 μm range. Nevertheless, the longer lag exhibited by the 4.5 μm emission implies the presence of a temperature gradient in the emitting region.

In K-band reverberation mapping studies of Seyfert 1 galaxies it has been found that the reverberation radius derived from the time lag is quite tightly correlated with $L_{{\rm opt}}^{0.5}$, where L_opt is the AGN optical luminosity (Suganuma et al. 2006; Koshida et al. 2009, 2014). This is consistent with the $R\propto {{L}^{0.5}}$ relation expected for dust in radiative equilibrium. However, Kishimoto et al. (2007) found that the K-band reverberation radii are a factor ∼3 smaller than the sublimation radii as predicted by Equation (2). One possible explanation is that the NIR dust emission is dominated by graphite grains; sublimation radii predicted by Equation (3) are a factor ∼3 smaller than the Silicate radii and thus much closer to the K-band reverberation measurements (see Figure 7). Several other explanations have been advanced for the apparent discrepancy between the measured dust radii and the sublimation radii predicted for the standard ISM dust composition. For example, the dust may include larger grains than the typical size (a ≈ 0.05 μm) assumed in Equation (2) (Kishimoto et al. 2007). Kawaguchi & Mori (2010) investigated the effect of anisotropic illumination of the torus inner wall by the accretion disk, which permits a smaller torus inner radius close to the disk plane. Another possibility, proposed by Pozo Nuñez et al. (2014), is that the torus is very optically thick in the NIR so that only emission from the facing rim of the torus inner wall is seen, leading to a "foreshortened" lag. Modeling of the time-dependence of the optical–NIR SED of NGC 4151 by Schnülle et al. (2013) suggests that the innermost dust is well below the sublimation temperature. This implies that the dust is located beyond the sublimation radius, suggesting anisotropic illumination or geometrical foreshortening, as envisaged by Pozo Nuñez et al. (2014).

Zoom In Zoom Out Reset image size

In Figure 7 we plot reverberation radii versus V-band luminosity ($\lambda {{L}_{\lambda }}(V)$) for both the 3.6 and 4.5 μm lags reported here and K-band results taken from Clavel et al. (1989), Suganuma et al. (2006) and Koshida et al. (2014). For this purpose, we use the lower limit to the AGN V-band luminosity of NGC 6418 inferred from $L_{{\rm H}\alpha }^{{\rm obs}}$, as described above.

We also plot Kishimoto et al. (2007)'s fit to the K-band lag data points,

$$\begin{eqnarray}&&{{R}_{\tau ,K}}=0.47{{\left( \frac{6\lambda {{L}_{\lambda }}({\rm V})}{{{10}^{46}}\ {\rm erg}\;{{{\rm s}}^{-1}}} \right)}^{1/2}}{\rm pc}.\end{eqnarray} \tag{ 6 }$$

With the caveat that the NGC 6418 points represent lower limits in luminosity, it can be seen that the mid-IR reverberation radii are located above the trend defined by the K-band lag times, as expected if the 3.6 and 4.5 μm emission is dominated by cooler dust located somewhat deeper in the torus. Equation (6) predicts ${{R}_{\tau ,K}}\gtrsim 11.6\;\times \;{{10}^{-3}}\ {\rm pc}$ for NGC 6418, given our lower limit on the V luminosity, implying that ${{R}_{\tau ,3.6}}\lesssim 2.7{{R}_{\tau ,K}}$ and ${{R}_{\tau ,4.5}}\lesssim 3.4{{R}_{\tau ,K}}$, respectively.

For dust grains in radiative equilibrium, the radius at which grains have a temperature T is approximately,

$$\begin{eqnarray}&&\frac{{{R}_{d}}}{{{R}_{{\rm sub}}}}\simeq {{\left( \frac{T}{{{T}_{{\rm sub}}}} \right)}^{\alpha }},\end{eqnarray} \tag{ 7 }$$

where R_sub is the sublimation radius and $\alpha \approx 2-2.8$ depends on the dust composition. In combination with Wien's Law, Equation (7) provides a rough estimate of the largest radius at which the dust contributes to the torus emission at a specific wavelength. For the typical ISM composition of Equation (2) ($\alpha =2.6$), we find ${{R}_{3.6}}/{{R}_{K}}\simeq 3.6$, ${{R}_{4.5}}/{{R}_{K}}\simeq 6.4$ and ${{R}_{4.5}}/{{R}_{3.6}}\simeq 1.8$. The values for ${{R}_{3.6}}/{{R}_{K}}$ and ${{R}_{3.6}}/{{R}_{K}}$ exceed the empirical upper limits determined from reverberation mapping, while the value of ${{R}_{4.5}}/{{R}_{3.6}}$ agrees with the ratio of the reverberation lags (${{R}_{\tau ,4.5}}/{{R}_{\tau ,3.6}}=1.3\pm 0.7$).

However, in clumpy torus models (e.g., Nenkova et al. 2008a, 2008b), there is a wide range of dust temperature within a typical cloud, which therefore emits a broad IR spectrum. In the models of Nenkova et al. (2008b), the bulk of the emission at $\lambda \lesssim 5\;\mu {\rm m}$ emerges from clouds at no more than twice the inner radius (see Nenkova et al. 2008b their Figure 13). Thus, the relative sizes of the reverberation radii at 3.6, 4.5 μm and K-band seem consistent with at least some clumpy tori models.

It is well established, mainly from Balmer line reverberation mapping (Greene et al. 2010; Bentz et al. 2013 and references therein) that the broad emission line region follows a similar $R\sim {{L}^{1/2}}$ size–luminosity relationship. For a given AGN luminosity, the BLR reverberation radius is a factor 4–5 smaller than the K-band dust emission reverberation radius (Suganuma et al. 2006; Koshida et al. 2014), as expected in the AGN unification paradigm. Interestingly, radii derived from Mor & Netzer (2012)'s SED fits suggest that the NIR emission component attributed to hot graphite dust clouds occupies a region intermediate between the BLR and K-band reverberation radii (see Koshida et al. 2014, Figure 13), consistent with the idea that this dust resides in the outer BLR clouds. In their analysis of mid-IR (12 μm) interferometric observations, Burtscher et al. (2013) find that although source sizes scale in a similar way with luminosity, there is a much larger scatter, with mid-IR source radii ranging from $\lesssim 4$ to $20\times {{R}_{\tau ,K}}$. A clearer picture of the structure of the AGN emission regions beyond the accretion disk is therefore beginning to emerge. Placing our results in this context, the reverberation radii derived from the 3.6 and 4.5 μm light curves are consistent with the variable emission at these wavelengths arising in the inner clouds of the torus. However, we note as a caveat that NGC 6418 exhibits an atypical optical spectrum for a Seyfert 1, with a relatively strong, broad Hα line but with a steep Balmer decrement, relatively weak narrow lines (for instance, the equivalent width of [O iii]$\lambda 5007$ is only ∼3 Å, that of narrow H $\alpha \sim 0.5$ Å) and with stellar emission dominating the optical continuum. This indicates that the BLR and AGN UV–optical continuum are subject to heavy extinction along the line-of-sight, raising the possibility that the circum-nuclear dust distribution may be more quasi-spherical than toroidal in nature.

We have presented initial results from the first year of a two-year campaign of IR (3.6 and 4.6 μm) and optical (B and V) monitoring of a sample of 12 Seyfert 1 galaxies using the Spitzer Space Telescope supported by ground-based optical observations. In NGC 6418, we have found a lag between the mid-IR and optical light curves, with a time delay of $37.2_{-2.2}^{+2.4}$ days ($31.2_{-1.9}^{+2.0}\times {{10}^{-3}}\ {\rm pc}$) at 3.6 μm and $47.1_{-3.1}^{+3.1}$ days ($39.5_{-2.6}^{+2.6}\times {{10}^{-3}}\ {\rm pc}$) at 4.5 μm, respectively. The 3.6 μm emission leads the 4.5 μm emission by $13.9_{-0.1}^{+0.5}$ days ($11.7_{-0.1}^{+0.4}\times {{10}^{-3}}\ {\rm pc}$). These results indicate that the dust emitting the bulk of the 3.6 and 4.5 μm emission is located at a distance $\approx 1$ light-month ($\approx 0.03\;{\rm pc}$) from the source of the AGN UV–optical continuum.

The nucleus of NGC 6814 appears to be heavily reddened, with a broad line Balmer decrement of Hα/Hβ ≥ 6. For this reason, we can only determine a lower limit for the intrinsic luminosity of the AGN and hence lower limits on the dust sublimation radii. The reverberation radii are a factor ∼2 smaller than the sublimation radius lower limit for silicate grains (sublimation temperature $\approx 1500\;{\rm K}$; ${{R}_{d,{\rm Si}}}\geqslant 60_{-21}^{+33}{\mkern 1mu} \times {{10}^{-3}}\ {\rm pc}$), but consistent with that for pure-graphite grains (sublimation temperature $\approx 1800$ K; ${{R}_{d,{\rm C}}}\geqslant 24_{-8}^{+13}\times {{10}^{-3}}\ {\rm pc}$). Reverberation radii derived from K-band variability studies of other Seyferts are similarly a factor ∼3 smaller than the silicate sublimation radius. It seems possible that some of the emission in the 3.6–4.5 μm range comes from hot graphite dust located within the region bounded by ${{R}_{d,{\rm C}}}$ and ${{R}_{d,{\rm Si}}}$, whose presence is suggested by SED model-fitting.

The 3.6 and 4.5 μm reverberation radii fall above the extrapolated K-band size–luminosity relationship by factors $\lesssim 2.7$ and $\lesssim 3.4$, respectively, while the 4.5 μm reverberation radius is only 27% larger than the 3.6 μm radius. This indicates a steeper temperature gradient than expected for optically thin dust in radiative equilibrium but is consistent with clumpy torus models, in which individual optically thick clouds emit strongly over a broad wavelength range.

We dedicate this paper to the memory of our great friend, colleague and mentor David Axon, who initiated this project and brought the collaboration together. We would like to thank Davide Lena for improvements to the manuscript and discussions on spectroscopic data reduction. S.U. thanks Lex Shaw for a generous instrumentation donation. This work is based (in part) on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Support for this work was provided by NASA through an award issued by JPL/Caltech, GO-80120. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France, and of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. 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. The SDSS-III web site is http://sdss3.org/ 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 Astrofísica 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. The Liverpool Telescope is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias with financial support from the UK Science and Technology Facilities Council. The Faulkes Telescope Project is an educational and research arm of the Las Cumbres Observatory Global Telescope Network (LCOGTN). BMP is supported by NSF grant AST-1008882.

We begin by defining a region within each image containing NGC 6418 and several nearby reference stars. Next, we determine a background value for this region by fitting a Gaussian to the histogram of pixel values: the peak yields the background value and the width its uncertainty. Sources are detected using the STARS program of XVISTA, which employs an algorithm based on the FIND procedure within DAOPHOT (Stetson 1987). Candidate objects which survive cuts in several parameters such as FWHM, sharpness and roundness, are selected for aperture photometry. We measure the brightness of each object using the PHOT program in XVISTA, which sums all counts within a circular aperture, including weighted contributions from pixels that lie partially inside the aperture. PHOT also measures the median pixel value within an annulus around each object to determine a local sky value and subtracts this from the object counts. Finally, the remaining object counts are converted to an instrumental magnitude.

The second stage of the analysis subjects the measured instrumental magnitudes to inhomogeneous ensemble photometry (Honeycutt 1992). Small differences in sky brightness, transparency, exposure time, and other factors can cause all objects in some particular exposure to appear slightly brighter or dimmer than average; ensemble photometry is designed to identify these systematic changes and remove their effects.

Honeycutt (1992) defines the equation of condition as

$$\begin{eqnarray}&&m(e,s)=m0(s)+{\rm em}(e),\end{eqnarray} \tag{ A1 }$$

where $m(e,s)$ is the instrumental magnitude of star $s$ in exposure $e$ and $m0$ is the intrinsic instrumental magnitude of that star. The "exposure magnitude," ${\rm em}$, of an image accounts for variations in extinction, exposure time, background intensity and other effects that are common to all sources in an image. We note that even without the transparency issues that are typical of ground observations, the Spitzer IR data will have small variations due to changes in orientation and background illumination of the space telescope. The quantity that we want to minimize is

$$\begin{eqnarray}&&\beta =\mathop{\sum }\limits_{e=1}^{ee}\mathop{\sum }\limits_{s=1}^{ss}{{[m(e,s)-m0(s)-{\rm em}(e)]}^{2}}w(e,s),\end{eqnarray} \tag{ A2 }$$

where $w(e,s)$ is the weight of each instrumental magnitude; we take its value to be $\sigma {{(m(e,s))}^{-2}}$. This technique yields the best fit value of $m0(s)$ for each source, assuming no intrinsic variability, and an empirical estimate of the uncertainty. In an ideal experiment, the uncertainty would be equal to that derived from the quadrature sum of the shot noise of the source, the sky noise and the detector read noise. This empirical estimate of the uncertainty is valid only for constant sources, such as the reference stars, but not for sources which vary intrinsically from one image to the next.

In an effort to give a comprehensive picture of the results obtained by different software packages used to determine the lag time between light curves we have included this appendix with Table 5 of all results. The table contains the analyses of individual and combined optical datasets versus the infrared channels of the Spitzer Space Telescope. The first column indicates the Spitzer channel. The table has three sections, one for each of the software packages we used. Columns 2 through 6 are values obtained for our in-house cross-correlation package. Of those, columns 2–4 represent the difference from the median to the 25% value of the interquantile range (IQR), the median of the distribution and the difference from the median to the 75% value of the IQR, respectively. Columns 5 and 6 are the mean and the standard deviation. Column 7 is the mean and the standard deviation for Peterson's code (Peterson et al. 2004). Columns 8–10 are Zu's (Zu et al. 2011) corresponding to the low, mid and high values of the lag.

Our simulations employ 1000 realizations of the light curves. Each synthetic light curve is generated by replacing each flux density measurement with an artificial datum. This consists of the measured magnitude plus a random deviate drawn from a Gaussian distribution with a mean of zero and standard deviation equal to the uncertainty in the measured value. We compute the CCFs and the corresponding weighted mean lags for each set of synthetic light curves to construct a distribution of the CCF centroids, the CCCD.

The CCFs are often not symmetrical functions, and the skewness of these functions affects the calculation of their centroids. The question is—how to select the significant portion of each distribution, while discarding the uninteresting wings? Figure A1 shows representative single realizations of the CCFs; it is obvious that the centroid of each CCF will depend on the range of data chosen for further calculation. In this work, we have adopted an algorithm that uses properties of each distribution itself to select the subset of measurements for the centroid calculation. First, we fit a cubic spline to the distribution in each realization, and compute the standard deviation, σ, between the spline and the data. We adopt 2σ as a measure of the dispersion within the CCF. We set a threshold in correlation which is the peak of the CCF minus this dispersion: $K={\rm peak}-2\sigma$. All the CCF values greater than K are then used to calculate the centroid of that particular CCF. The fitted spline is shown together with the computed CCF(τ) data points. We found that for optical versus 3.6 μm, the top 24% of CCF data was used, for the optical versus 4.5 μm the top 23%, and for the 3.6 μm versus 4.5 μm the top 6%. The threshold clearly is dependent on the noise characteristics of the underlying light curves which explains why the Spitzer light curves have a smaller data percentage used in the centroid calculation.

Zoom In Zoom Out Reset image size

After calculating the centroid of each realization of the CCF in this manner, we then combine all the centroids to create the CCCD for that pair of light curves. We choose the median value in the CCCD as the time lag between the two light curves.