VARIATION OF INNER RADIUS OF DUST TORUS IN NGC4151

The unified model of active galactic nuclei (AGNs; Antonucci 1993; Urry & Padovali 1995) assumes the existence of a dust torus that surrounds the central hierarchical structure consisting of a supermassive black hole, an accretion disk, and a broad-line region (BLR). The dust grains in the torus absorb the UV/optical continuum emission from the accretion disk and re-radiate in the near-infrared (IR) wavelength region with some lag time corresponding to the light travel time from the accretion disk to the inner radius of the dust torus. Since the heated dust eventually sublimates at a constant temperature of 1500–1800 K, a more luminous AGN should have a larger dust torus and hence a larger lag time, yielding a correlation of Δでるたt ∝ L^0.5 between lag time Δでるたt and AGN luminosity L (Barvainis 1987, 1992).

In fact, based on the long-term multicolor monitoring data from the MAGNUM project (Yoshii 2002; Yoshii et al. 2003) and available archival data, Suganuma et al. (2006) recently reported such a correlation from a sample of AGNs spanning a wide range of absolute V-magnitude from M_V = −15 to −24. The optical luminosity is a good indicator of UV luminosity, because their variations are well synchronized with each other (Edelson et al. 1996). Therefore, the reported correlation in the V band indicates that the inner radius of the dust torus scales with the UV luminosity as well, which provides a strong piece of evidence for dust reverberation proceeding in the central region of AGNs.

Given the above result of scaling among AGNs of different luminosities, it is possible to examine how Δでるたt would scale with changing L in an individual AGN. This is currently the only diagnosis on its inner structure, which is too small to be resolved by direct imaging. However, since the range of flux variation in an AGN is as small as <0.1 mag, detection of a change in Δでるたt is indeed a challenge and needs an elaborate method that enables precise estimation of Δでるたt.

We here attempt to detect a change, if any, of Δでるたt for the case of a nearby Seyfert 1 AGN of NGC 4151, using the six year long accurate and frequent monitoring data for bright AGNs in the program of MAGNUM.

In Section 2, we show details of the observation of NGC4151. In Sections 3 and 4 we formulate estimations of Δでるたt and Ls by making use of model light curves. In Section 5, we apply the formulation to the case of NGC 4151 and derive its time-changing track in the Δでるたt versus L diagram. In Section 6, we discuss the timescale of dust formation in the central region of NGC 4151.

NGC4151 is one of the bright AGNs most frequently observed since 2001 by the MAGNUM telescope. The monitoring period analyzed in this Letter covers a length of 2000 d from MJD 52,000 to 54,000 (2001–2006), and photometric images were obtained for about 220 nights in the V and K bands. The observations, photometry, and data analysis are mostly in common with those in Minezaki et al. (2004), and we do not repeat them except for some important basics below.

Each night we took images of the AGN and the reference star alternately using the multicolor imaging photometer (MIP) with a viewing angle of 15, and performed relative photometry of the AGN with an aperture size of ϕ = 83. This size was chosen to minimize the photometric error and the flux variation in the aperture by the seeing effect.

The flux in the aperture contains the flux from the host galaxy, especially from its bulge. The host component was estimated by the model fitting to its surface brightness as $f_{V,\rm{gal}}=17.95\ \mathrm{mJy}$ (Sakata et al. 2009) and $f_{K,\rm{gal}}=44.22\ \mathrm{mJy}$ (Minezaki et al. 2004), and we subtracted these contributions from the V and K light curves, respectively. The narrow line flux was also significant in the V band for NGC4151, so we subtracted 10.41 mJy following Sakata et al. (2009). Furthermore, the IR flux contains the flux from the accretion disk, which makes Δでるたt systematically shorter than the actual lag. According to Tomita et al. (2006), the IR disk component should vary to synchronize with the optical variation that is originated genuinely from the accretion disk. This component in the K band was then estimated as $f_{K,\rm{acc}}=(\nu _{V}/\nu _{K})^{\alpha } f_{V}$, and we also subtracted this contribution from the K light curve by assuming αあるふぁ = 0 for the disk emission between the V and K bands.

The resulting V and K light curves of NGC 4151 from 2001 to 2006 are shown in Figure 1. The V light curve shown here is consistent with the optical continuum light curves reported by Shapovalova et al. (2008). We see that the optical flux of NGC 4151 reached a minimum around MJD 52,000 and then increased until reaching a maximum at MJD 52,700. Thereafter, it turned to decrease again to a minimum at MJD 53,500. Apparently, the flux variation in the K band followed that in the V band with clear time delay.

Zoom In Zoom Out Reset image size

The light curve shown in Figure 1 is complex, with many peaks and valleys. From it we need to extract a table of Δでるたt as a function of continuum flux which has sufficient precision for us to test the hypothesis that Δでるたt is proportional to f^0.5_V as the flux varies, and/or to understand any difference from that expected result.

The techniques used previously are not adequate. In our series of papers on dust reverberation in AGNs (Minezaki et al. 2004; Suganuma et al. 2006; Tomita et al. 2006), we estimated Δでるたt and its error based on the Cross-Correlation Function (CCF) analysis by making use of the structure function (SF; White & Peterson 1994) for the interpolation of optical flux variability during the monitoring period. In the equal sampling scheme of interpolation which we used for the CCF analysis (Suganuma et al. 2006), the number of data points generated at equal intervals of time can exceed that actually observed. In such a case, the simulated light curve shows larger flux variation compared to the overall trend seen from the observed light curve, especially when the data points are generated far from the location of observed epochs. This necessarily gives rise to too large an error of Δでるたt.

We here adopt the bi-directional interpolation (BI) scheme, where all V-band flux data are paired with K-band flux data and at least either the V or the K flux in each pair is measured by actual observation. While the BI scheme overcomes the problem of generating too many data points for CCF calculations, the CCF still shows an unreal peak with increasing number of data pairs. This drawback is remedied by estimating Δでるたt from the CCF centroid around the peak which is interpreted as a luminosity-weighted distance of the dust torus from the AGN central energy source (Koratkar & Gaskell 1991). As a practice of determining the CCF range for centroid estimation, the threshold is usually set to be 0.8 times the CCF peak value. If the threshold is larger than 0.9, Δでるたt is influenced by an unreal peak of the CCF. On the other hand, if it is smaller than 0.7, Δでるたt is necessarily influenced by sub-features such as a flared foot of the CCF.

We have so far used a simple mean of observed V-band fluxes as an indicator of optical luminosity that corresponds to Δでるたt, but this could be seriously biased by two factors. First, the flux mean depends on the sampling frequency of the light curve during the monitoring period. When an AGN is often observed in the bright phase rather than in the faint phase, the simple flux mean preferentially gives a brighter flux value. Second, the flux mean depends on the shape of the light curve. When an AGN undergoes an abrupt burst on a flat light curve, the simple flux mean gives a value which is more or less the same as the value of the flat part of the light curve, although the flat part contributes very little to the value of the cross-correlation (CC) coefficient for the estimation of Δでるたt.

Considering the above, we take a mean of the V-band fluxes generated at equal intervals of time, after excluding the data which statistically show no flux variation in the time bins. To pick up such data, we use the equally interpolated optical light curve for which the size of the time bins is taken as the mean interval of monitoring observations. The variability in each bin is tested by χかい² with a significance level of 0.99, and the data which are judged as showing no flux variation are excluded in calculating the mean of the V-band fluxes. The distribution of this flux mean for all simulated light curves is then used to derive an average optical flux f_V and its error.

We next use model light curves to verify that our analysis procedure reliably finds a characteristic lag Δでるたt and its associated flux f_V for each of the individual variability events that combine to make the complex light curve shown in Figure 1. In order to quantify the above estimation of a characteristic lag Δでるたt and its associated flux f_V as well as their errors, we model the light curves to calculate CCFs. Given a V light curve f_V(t), we construct a corresponding K light curve as f_K(t) = f_V(t − Δでるたt), where Δでるたt/50 = (f_V/200)^γがんま for γがんま = 0.5 ± 0.3 which roughly reproduces the lags seen in Figure 1. Here, the units are days in time and arbitrary in flux. Then we generate 30 data points along each of the V and K light curves, and disperse the points according to a Gaussian distribution of photometric errors taken from the actual observations by the MAGNUM telescope.

The models include three characteristic shapes of (1) linear bottom, (2) quadratic bottom, and (3) bursting Gaussian peak, as shown in order from top to bottom in the left column of Figure 2. The models 1 to 3 are shown in Figure 2 to examine the dependence on shape of the time-symmetric incident V light curve used for CCF calculation. In addition to all these models with a standard index of γがんま = 0.5, we consider two cases of γがんま = 0.2 and 0.8 for each of light curve shapes 1 and 2. We note, however, that γがんま>0.5 is considered only for academic purposes and contradicts the geometrical 1/r² dilution of the incident flux from the central energy source.

Zoom In Zoom Out Reset image size

Following Suganuma et al. (2006), we applied the BI method to each of the model light curves to calculate the CCF and also the CC centroid distribution (histogram) by Monte Carlo simulations, as shown in the middle column of Figure 2 for γがんま = 0.5. The results of the characteristic Δでるたt and its error are given on the panels of this column and in Table 1. The results of Δでるたt for models 1 to 3 give values around 80 d. This indicates that Δでるたt could be estimated independently of the shape of the light curve.

Now we derive an average flux f_V for each of the models 1 to 3, following the method described in Section 3. The horizontal arrows in the left column of Figure 2 specify the periods in which the data are rejected in the estimation of flux mean due to their invariability. We calculate the distribution of this flux mean for the light curves generated by Monte Carlo simulations, and present it as a histogram from which the characteristic f_V and its error are derived, as shown in the right column of Figure 2 for γがんま = 0.5. The results of f_V and its error are given on the panels of this column and in Table 1.

To be consistent with the measured value of the characteristic Δでるたt, we compare our estimate of the characteristic f_V with that from the relation Δでるたt/50 = (f_V/200)^γがんま adopted in the CCF analysis. The relative difference of δでるたf_V/f_V ≡ (f_V − f_V(Δでるたt))/f_V(Δでるたt) is shown in Table 1. We see that |δでるたf_V/f_V| for γがんま = 0.5 is less than 10% for all the models 1 to 3, which is considerably improved over our previous estimate of 40% by using the simple flux mean.

The results of f_V for different values of γがんま = 0.2 and 0.8 tabulated in Table 1 show little difference from that for γがんま = 0.5, which indicates that the estimated f_V is almost independent of this parameter. In addition, various other light curves in asymmetric shapes also give essentially the same result.

Given that the consistent estimation of Δでるたt and f_V has been confirmed as above by using the model light curves, we apply it to the real data for NGC 4151. We divided the whole monitoring period into eight separate periods, each of which contains a single feature of flux variation, as above, with a sufficient number of more than 40 flux data pairs for the CCF analysis. The middle date of each period is taken as an epoch that represents each period. Following the procedure described in Section 3, the Monte Carlo simulation was repeated 500 times for each period to form the distribution of CC centroid lag from which Δでるたt was obtained as the median and its error as ± 34.1 percentile from the median. Similarly, we calculated the distribution of flux mean from which f_V and its error were obtained. The results for eight separate periods are given in Table 2 as well as in Figure 3. We here note that the K-band data of NGC 4151 were well sampled and the observed K light curves were smooth enough. Therefore, in order to avoid additional IR flux variation which consequently overestimates the errors in Δでるたt, use of the SF for the interpolation was turned off and the simple linear interpolation was adopted for the K-band data in making pairs with the V-band data.

Zoom In Zoom Out Reset image size

It is apparent that Δでるたt is not constant in time. For the most extreme case, the difference of Δでるたt between the earlier periods 1 through 4 and the later periods 6 through 8 is larger than the measurement error by many σしぐまs. Since Δでるたt corresponds to a light-travel distance to the surrounding dust from the AGN center, this should be the first convincing evidence that the inner radius of dust torus did change in an individual AGN (cf. Oknyanskij et al. 1999; Minezaki et al. 2004).

A change of inner radius of the dust torus should occur reflecting the variation of the incident UV/optical flux from the central energy source. Some theoretical models suggest the existence of a sublimation radius inside which dust grains are sublimated (Barvainis 1987, 1992). The sublimation radius expands when the UV/optical continuum becomes bright, then it retreats when the continuum becomes faint. If dust grains were sublimated or replenished immediately after the UV/optical flux variation, like the model we assumed in Section 3, the expected changes of Δでるたt and f_V should trace the simple relation Δでるたt ∝ L^0.5 which has the slope indicated by the thick line in Figure 3. We see from this figure that the observed change of Δでるたt follows the variation of f_V in overall trend as expected. These results indicate that the dust torus is not a distinct, separate physical structure but is part of a continuous structure starting from the BLR component out to the dust component with its inner dust edge set by the sublimation radius that changes with the UV/optical flux variation (Suganuma et al. 2006; Nenkova et al. 2008).

On the other hand, however, the observed time-changing track does not exactly trace the simple relation Δでるたt ∝ L^0.5. The lag Δでるたt at period 1 is well above that expected for M_V, and does not change so much with brightening of about 1 mag from period 1 to 3. Then, with subsequent fading from period 3 to 6, Δでるたt still remains constant until period 4, and starts to rapidly decrease thereafter. These deviant behaviors suggest that the inner dust torus did not reach an equilibrium state immediately after the UV/optical flux variation.

In the case of period 1, the light curves in Figure 1 indicate that there was a brighter period before period 1. The expanded inner radius of the dust torus at the brighter period could still remain larger than the radius that would be expected at period 1. In this context, the behavior of almost constant Δでるたt followed by its rapid decrease in the fading period places a constraint such that the timescale of dust replenishment in the central region of an AGN should be as long as the elapsed time of about one year from period 4 through 6. A promising mechanism of such replenishment is re-formation of dust grains rather than their infall from the outer region. This is understood by considering that the decrease of Δでるたt from period 4 to 6 is 24.5 d or a light-travel distance of 6.35 × 10¹¹ km in the interval of 309 d between the middle dates of these periods. If such a decrease of Δでるたt would occur with redistribution of dust grains supplied from the outer region of the dust torus, the infall velocity would be 2.4 × 10⁴ km s⁻¹ or 7.9% of the light velocity. This seems highly unlikely, when we compare this infall velocity with a few ×10³ km s⁻¹ for the velocity dispersion of BLR clouds which exist just near the inner dust torus. Consequently, we conclude that re-formation of dust grains did occur in the central region where they had been sublimated.

After the dust re-formation period, Δでるたt turned almost constant again from period 6 through 8. This indicates that the dust sublimation radius became larger than the actual inner radius of the dust torus and dust sublimation would start. The length of the period of constant Δでるたt suggests that the timescale of dust destruction could be longer than about one year. Exact measurement of this timescale might have been possible if the light curves had been obtained well after period 8.

The formulation for improving the accuracy of lag and luminosity measurements was applied to the light curves for eight separate periods containing single features. This formulation, when applied to the light curves for the whole monitoring period containing several features, gives an almost average of their respective Δでるたt estimates derived from single features, and is still useful in discussing an overall lag–luminosity relationship for a sample of many AGNs with a wide range of luminosity.

We thank the staff at the Haleakala Observatories for their help with facility maintenance. This research has been supported partly by the Grants-in-Aid of Scientific Research (10041110, 10304014, 11740120, 12640233, 14047206, 14253001, 14540223, and 16740106) and the COE Research (07CE2002) of the Ministry of Education, Science, Culture and Sports of Japan.