© ESO, 2011

1. Introduction

An optically- and geometrically-thick, circumnuclear dusty region (= “dust torus”) is one of the key ingredients of the unification scheme of active galactic nuclei (AGN; Antonucci 1993). It is held responsible for angle-dependent obscuration of the central accretion disk and broad-line region, explaining the apparent difference between type 1 (unobscured line-of-sight) and type 2 (obscured line-of-sight) objects. Its dust content absorbs part of the accretion disk’s UV/optical emission and reemits in the infrared (IR). Because the torus resides on parsec scales, direct observations have been challenging because its angular size is much smaller than the resolution power of single telescopes in the IR.

With the advent of IR long-baseline interferometry and reverberation mapping, it is now possible to determine the sizes of the dusty region around nearby AGN (e.g., Jaffe et al. 2004; Tristram et al. 2007; Beckert et al. 2008; Kishimoto et al. 2009b; Burtscher et al. 2009; Pott et al. 2010). Of special interest are observations of type 1 AGN where the torus is presumably seen closer to face-on than in type 2 AGN without too many complications caused by obscuration effects (Kishimoto et al. 2009a, 2011b). While sizes are the primary parameter to extract from interferometry data, reverberation mapping uses the fact that any variability in the accretion disk emission causes a delayed response in the near-infrared continuum emission from the hottest, inner region of the dusty torus because of the light-travel time from the disk to the torus. Both interferometry and reverberation mapping found that the time lags and sizes generally agree with theoretical expectation of the sublimation radius based on thermal equilibrium of large graphite grains or anisotropic accretion disk emission (e.g., Kishimoto et al. 2007). These results were confirmed by directly measuring near-IR K-band ring radii using the Keck interferometer (Kishimoto et al. 2009b, 2011a). It should be noted that variability of the accretion disk occurs on multiple time scales (from intraday variability up to many years) and amplitude ranges. The torus has been observed to react on variability in the range of few days to as long as decades (e.g., Glass 1992, 1997, 2004; Oknyanskij et al. 1999; Gallimore et al. 2001; Minezaki et al. 2004; Suganuma et al. 2006).

Infrared interferometry has proven as a tool to access more detailed information about the dust torus than just sizes. Based on mid-IR interferometric observations, Jaffe et al. (2004) argued that the dust in the torus is confined to small clumps rather than smoothly distributed, as theoretically and observationally suggested in early studies (e.g., Krolik & Begelman 1988; Tacconi et al. 1994). Subsequently clumpy torus models have been successful in reproducing spectral (e.g., Nenkova et al. 2002; Polletta et al. 2008; Ramos Almeida et al. 2011; Alonso Herrero et al. 2011) and spatial information (e.g., Hönig et al. 2006, 2008; Schartmann et al. 2008; Hönig et al. 2010). Moreover, we recently showed that the dependency of visibility on wavelength and spatial frequency can be used to constrain the radial brightness distribution of the dust (Kishimoto et al. 2009b; Hönig et al. 2010; Hönig & Kishimoto 2010; Kishimoto et al. 2011b). Using such an analysis, it has been found that the brightness distribution and the slope of the mid-IR spectra are strongly correlated. Both parameters are probably tracers for the radial distribution of the dust, so that combining interferometry and photometry provides a tool to characterize the properties of the dust torus in different objects.

At present IR interferometry has some tight limitations regarding brightness and spatial resolution, so that only the brightest nearby sources can be studies in detail. Reverberation mapping, on the other hand, may serve as a complimentary tool for objects that are not accessible by interferometry. In this paper we present a correspondence in reverberation mapping to the interferometric method of determining the brightness distribution in type 1 AGN. For that we use a simplified version of a clumpy torus model, which has been successfully applied to interferometric data (Kishimoto et al. 2009b, 2011b). In Sect. 2 we describe a theoretical approach to handle temperature variations in a dusty medium, generalizing the seminal work of Barvainis (1992). Based on this theory, we outline a simple model for type 1 AGN in Sect. 3 and present model light curves and cross-correlation functions at IR wavelength that show the dependence on the brightness distribution of the torus. As a proof of our concept we use our model to reproduce the K-band light curve of NGC 4151 in Sect. 4. The results are summarized in Sect. 5.

2. Luminosity and temperature variations in a dusty medium

In this section we describe theoretically how variability in the luminosity of the incident radiation changes the temperature of the dust. The new dust temperatures can then be used to calculate the changed emission of the dust. Under the assumption that heating and cooling times are negligible, no iterations are required.

For dust grains in radiative equilibrium of size a at distance r from a radiating source with spectrum L_ν, the received power per dust grain, ${\mathit{L}}_{\mathrm{rec}}\mathrm{=}{}^{\mathrm{\int }}\mathit{L}_{\mathit{\nu }}\mathit{/}\mathrm{\left(}\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}\mathrm{\right)}{\mathit{Q}}_{\mathrm{abs;}\hspace{0.17em}\hspace{0.17em}\mathit{\nu }}\mathit{\pi }{\mathit{a}}^{\mathrm{2}}\mathrm{d}\mathit{\nu }\mathrm{=}{\mathit{L}}_{\mathrm{abs}}{\mathit{a}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathrm{4}{\mathit{r}}^{\mathrm{2}}\mathrm{\right)}$ equals the emitted power L_em = ^∫4πa²  Q_abs;ν  πB_ν(T)dν, where Q_{abs;      ν} is the absorption efficiency of the dust and B_ν(T) denotes the Planck function of a black-body at temperature T. Equating L_rec and L_em leads to ${\mathit{L}}_{\mathrm{abs}}\mathrm{=}\mathrm{16}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}{\mathit{Q}}_{\mathrm{abs;}\hspace{0.17em}\hspace{0.17em}\mathrm{P}}\mathrm{\left(}\mathit{T}\mathrm{\right)}{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{T}}^{\mathrm{4}}\mathit{,}$(1)where ${}^{\mathrm{\int }}\mathit{Q}_{\mathrm{abs;}\mathit{\nu }}\mathit{\pi }{\mathit{B}}_{\mathit{\nu }}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{d}\mathit{\nu }$ has been replaced by the Planck mean absorption efficiency Q_abs;P(T)·σ_SBT⁴. If the emission source is varying with time, i.e. dL/dt ≠ 0, the temperature of the dust grain will change. From Eq. (1) we find $\frac{\mathrm{d}{\mathit{L}}_{\mathrm{abs}}}{\mathrm{d}\mathit{T}}\mathrm{=}\frac{\mathrm{4}{\mathit{L}}_{\mathrm{abs}}}{\mathit{T}}\mathrm{+}{\mathit{L}}_{\mathrm{abs}}\frac{{\mathit{Q}}_{\mathrm{abs;P}}^{\mathrm{\prime }}}{{\mathit{Q}}_{\mathrm{abs;P}}}\mathrm{·}$(2)Here, we define ${\mathit{Q}}_{\mathrm{abs;P}}^{\mathrm{\prime }}\mathrm{\equiv }\mathrm{d}{\mathit{Q}}_{\mathrm{abs;P}}\mathit{/}\mathrm{d}\mathit{T}$ as the T-derivative of the Planck mean absorption efficiency of the dust. Both the Planck mean absorption efficiency and its derivative can be calculated from the optical properties of dust and only depend on the dust temperature.

Absorption efficiencies of astronomical dust species (graphite, silicates) and compositions typically peak in the UV/optical in about the same ν-range where many astronomical sources radiate (e.g., stars or AGN accretion disks). Therefore we assume for practical purposes¹ dL_abs/L_abs ≈ dL/L (with $\mathit{L}\mathrm{=}{}^{\mathrm{\int }}\mathit{L}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }$) and obtain the variability function $\mathrm{d}\mathit{T}\mathrm{=}\frac{\mathrm{d}\mathit{L}}{\mathit{L}}{\left(\frac{\mathrm{4}}{\mathit{T}}\mathrm{+}\frac{{\mathit{Q}}_{\mathrm{abs;P}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{Q}}_{\mathrm{abs;P}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}\right)}^{-1}\mathit{.}$(3)For a given temperature T of the dust (e.g., as obtained by equilibrium calculations) the change of temperature dT can be analytically calculated for any change in incident radiation. However, the change in temperature depends on the temperature itself. But this additional T-dependence is minor. In Fig. 1 we show the dependence of the fractional change of temperature dT/T on the dust temperature T for a constant dL/L = 1. The black-solid line shows this dependence for standard ISM dust with 53% silicates, 47% graphite grains (optical constants from Draine 2003), and a grain size distribution according to Mathis et al. (1977, MRN). For temperatures below  ~100 K and above 1000 K the fraction is nearly constant at about 0.16–0.19. Between 100 K and 1000 K the Planck mean opacity has a “knee” that is reflected in dT/T as a bump up to 0.24 at 400 K. The blue-dashed line in Fig. 1 shows large graphite grains with sizes 0.07   μm < a < 1   μm (MRN size distribution). Below 100 K the graphite grain dT/T wiggles around the same mean value of 0.17 as the standard ISM dust. From 100 K to the sublimation temperature (~1500–2000 K) it increases to about 0.2. Overall, these changes are not very big and may be approximated by a constant factor. Indeed, if we assume the black-body limit (${\mathit{Q}}_{\mathrm{abs;P}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}$) we obtain $\frac{\mathrm{d}\mathit{T}}{\mathit{T}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{4}}\frac{\mathrm{d}\mathit{L}}{\mathit{L}}\mathit{,}$(4)where the relative change in temperature is directly linked to the relative change of luminosity without any extra T-dependence, and the constant factor is 0.25. This is shown as a red-dashed line in Fig. 1

Fig. 1 Dependence of the relative temperature change dT/T on the actual dust temperature for a variability in luminosity of dL/L = 1. The black solid line shows standard ISM (see text for details) while the blue dashed line is for large graphite grains. The red dotted line indicates the black-body limit.

It should be emphasized that this treatment is not restricted to a special class of objects (e.g., AGN, see Sect. 3) but can be used for dusty media in general. The presented treatment of temperature variations implicitly assumes that the reemission occurs in the optically thin wavelength regime. Otherwise an additional “self-heating term” has to be included on the left side of Eq. (1) and the treatment becomes iterative. However, even for a clumpy dust torus, where the primary emitting regions in the near- and mid-IR are the directly illuminated, optically thin layers of otherwise optically thick dust clouds, the temperature and emission profile is dominated by direct heating from the central source (Hönig & Kishimoto 2010) and the presented treatment should be applicable.

3. IR variability in AGN

3.1. Type 1 variability model

One of the applications of a variability model can be reverberation mapping measurements of AGN. Fluctuations in the incident radiation from the accretion disk, which primarily emits at UV/optical wavelengths, causes the re-emission of the dust in the torus to change over time (e.g., Glass 2004). Because the dust is located at some distance from the accretion disk, the variability in the IR is delayed with respect to the optical emission by the light-travel time from the disk to the dust. This method has been exploited to measure the near-IR radius in nearby AGN (e.g., Glass 2004; Suganuma et al. 2006).

Reverberation mapping is usually performed in type 1 AGN where obscuration effects are, in general, significantly lower than in type 2 AGN. In these objects the dust-/brightness distribution is projected onto a “disk” where the hottest dust emits at the smallest radii and the cooler dust at longer distances. Kishimoto et al. (2009a) presented a simple analytic model for the IR emission, which has been shown to be a good representation of a clumpy dust torus in type 1 AGN (Hönig & Kishimoto 2010). In this concept the dust emission S_ν(r) (dominated by the optically thin layer of an otherwise optically thick dust cloud; see Sect. 2) at a distance r from the AGN is characterized by the corresponding radiative equilibrium temperature T. The total emission L(tor)_ν of the torus is then calculated by integrating from the sublimation radius to a (well-selected; see Hönig & Kishimoto 2010) outer radius r_out, and accounting for a radial dust distribution that is parametrized as a power law with index α: $\mathit{L}\mathrm{(}\mathrm{tor}{\mathrm{)}}_{\mathit{\nu }}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{sub}}}^{{\mathit{r}}_{\mathrm{out}}}\mathit{\pi }{\mathit{B}}_{\mathit{\nu }}\mathrm{(}\mathit{T}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{)}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{sub}}}\right)}^{\mathit{\alpha }\mathrm{+}\mathrm{1}}\mathrm{d}\mathit{r}\mathit{.}$(5)In this concept the term (r/r_sub)^α is directly related to the surface filling factor of the torus (for details see the discussion in Hönig & Kishimoto 2010, Sect. 2.2). The dust equilibrium temperature can be obtained either by approximating a power-law (Barvainis 1987) or self-consistently from Eq. (1) (and using the sublimation radius and temperature as reference values) by solving the following equation: $\frac{{\mathit{Q}}_{\mathrm{abs;P}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{Q}}_{\mathrm{abs;P}}\mathrm{\left(}{\mathit{T}}_{\mathrm{sub}}\mathrm{\right)}}{\left(\frac{\mathit{T}}{{\mathit{T}}_{\mathrm{sub}}}\right)}^{\mathrm{4}}\mathrm{=}{\left(\frac{{\mathit{r}}_{\mathrm{sub}}}{\mathit{r}}\right)}^{\mathrm{2}}\mathrm{·}$(6)For practical purposes this equation can be solved using a pre-calculated look-up table for the left term and interpolating for the respective r.

With these prerequisites it is possible to calculate wavelength-dependent light curves of type 1 AGN based on Eq. (3), considering the traveling time of any variability through the torus.

3.2. Model results

The model presented in Sect. 3.1 has only one parameter, α, that represents the radial brightness distribution of the torus. We have shown that this parameter is well-correlated with the actual radial distribution of the dust at least in type 1 AGN even in more complicated models (Hönig & Kishimoto 2010). In consequence the radial dust distribution can be probed by observables such as the mid-IR spectral index or the comparison between mid-IR and near-IR interferometric size (Hönig et al. 2010; Kishimoto et al. 2011b). Here we analyze the effect of radial brightness distribution (and, in turn, radial dust distribution) on IR light curves of AGN.

3.2.1. Light curves

Fig. 2 Model light curves at 2.2   μm (red) and 8.5   μm (blue) for different values of the radial brightness distribution power law index α. From dark to light the distribution changes from very extended to very compact. The black dashed line illustrates the variability signal from the accretion disk.

In Fig. 2 we show light curves Δf(t)/f₀ (where f₀ = f(t₀) is the flux at time t₀, the starting time of the monitoring, and Δf(t) = f(t) − f(t₀) is the flux difference between t₀ and t) at 2.2   μm and 8.5   μm for various radial brightness distributions, parametrized by the power-law index α. The initial AGN variability signal is a step function with a width of 0.5   r_sub/c and an amplitude of 0.5   ΔL/L, and the light curves essentially show the transfer function of this signal. The lighter colors correspond to more compact distributions. In these cases most of the light is coming from a region close to the sublimation radius, i.e. the dust is concentrated in the inner torus. The darker colors have more extended brightness distributions. Here a significant fraction of the total dust mass is located at larger radii and directly heated by the AGN.

We plot the light curves as a function of the intrinsic time scale t_sub = r_sub/c, which corresponds to the light-travel time from the source to the sublimation radius r_sub. In this way the variability signal from the AGN can be considered a tomographic device to map out the brightness distribution in the torus, and elapsed time and distance from the AGN become equivalent independent of the object’s luminosity (t_sub = r_sub/c ∝ L^1/2, see Sect. 2).

The near-IR (2.2   μm) light curves in Fig. 2 show that the more concentrated distributions are stronger peaked, while the distributions with α closer to 0 have longer tails. In the very compact objects almost all of the near-IR light comes from the peak black-body² emission of the hottest dust. When the dust distribution becomes more extended, then the Wien tails of the cooler dust emission start to contribute relatively more because the relative amount of hot dust decreases.

Fig. 3 Model cross-correlation function (CCF) for light curves of AGN variability at 2.2   μm (red) and 8.5   μm (blue). From dark to light the distribution changes from very extended (α = 0.0) to very compact (α =  −2.0). The dashed line marks a lag time of rsub/c.

The behavior in the near-IR has its correspondence at longer wavelength. The mid-IR (8.5   μm) light curves with a compact brightness profile are peaked in the inner torus where the hot dust is located. Over time the brightness decays quickly, but not as fast as in the near-IR. As the brightness distribution becomes shallower (more extended), the initial peak at r_sub/c vanishes, the light curves get flatter, and for α >  −0.5 they keep on increasing out to  ≥ 10 × r_sub/c. The difference between the mid-IR and the near-IR is that in the mid-IR the emission comes predominantly from either the Rayleigh-Jeans tail of hot dust emission in the inner torus or from the peak black-body emission of cooler dust at longer distances. In the near-IR the flux originates predominantly in the black-body peak of hot dust with little contribution from the steeply dropping Wien tail of cooler dust. For steep brightness profiles, the hot-dust Rayleigh-Jeans tail dominates the mid-IR emission because of the lack of extended dust. Around α ≈  − 1, the contribution of the black-body peak from extended dust becomes about equal to the hot dust Rayleigh-Jeans tail, which takes over for even shallower brightness profiles.

3.2.2. Lag times

Because of the contribution from different emission regions changes as a function of wavelength and brightness distribution, we can expect delays between the peak of the emission in light-curves. Classical reverberation mapping tries to quantify the lag time between the AGN variation and its correspondence in the near-IR by means of the cross correlation function (CCF) $\mathit{CCF}\mathrm{=}\mathrm{\int }{\mathit{f}}_{\mathit{X}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}{\mathit{f}}_{\mathit{Y}}\mathrm{(}\mathit{t}\mathrm{-}\mathrm{\Delta }\mathit{t}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{t,}$(7)where f_X(t) and f_Y(t) are light curves at wavebands X and Y, and Δt is the introduced time lag between both light curves. This method can be applied, in principle, also to the mid-IR. We note that in the framework of this paper we define lag times as the peak in the CCF.

In Fig. 3 we show the CCF of the AGN variability signal with the 2.2   μm and 8.5   μm light curves, respectively. As expected, the time lag between AGN variability and near-IR emission is close to r_sub/c and the CCF is quite narrow. That it is actually slightly longer than unity is in part a result of the used variability function (step function with width 1/2t_sub) and the smoothened-out response by the dust. Evidently, the lag time depends slightly on α, so that extended emission profiles show slightly longer lag times (~1.5 × r_sub/c) while the most compact brightness distributions are at  ~ 1.1 × r_sub/c.

The mid-IR light curves show a different behavior. If the brightness distribution is compact, the peak in the CCF is still well-defined and close to r_sub/c. However, when the distribution becomes more extended, the CCF becomes very broad. For α ≤  − 0.5 the peak is no longer well-defined and shifts to very long time lags. This is expected from the light curves and reflects the change in contribution from hot and cooler dust to the mid-IR emission (see Sect. 3.2.1).

The change in CCF (or light curve) with α may be used to distinguish between compact and extended dust distributions. Interestingly, because the near-IR CCFs and light curves are peaked around r_sub/c for all α, it is not necessary to take the AGN variability signal as a reference, but instead use the near-IR as the CCF-reference signal for the mid-IR. While the actual size of r_sub cannot be determined in that way, the dust distribution relative to r_sub is still accessible. In Fig. 4 we show the CCF between near-IR and mid-IR light curves. For each α the near-IR light curve was used as reference for the corresponding mid-IR light curve. As expected, the compact brightness distributions have no time lag between near- and mid-IR, illustrating that the mid-IR emission originates from the Rayleigh-Jeans tail of the hot dust emission. With α <  −1 the CCF curves become significantly wider and the lag shifts to longer time lags. For 0.0 < α <  −0.5 the time lags reach 5 to 20 × r_sub/c. It is therefore possible to exploit IR light curves to determine the radial brightness distributions and, in turn, the radial dust distribution as an alternative to interferometry. We note that the actual shape of the CCFs also depends on the auto correlation function of the input signal which does not change the overall effects and conclusions, however.

3.3. Complications when dealing with AGN light curves

One of the complications when trying to model optical and IR light curves of AGN is a possible dependence of the variability on frequency. While the model uses a frequency-integrated version of dL/L, observations are usually in one specific waveband. Meusinger et al. (2011) report Δf_λ ∝ λ^-2, converting to Δf_ν ∝ ν⁰ (i.e. independent of frequency), so that Δf_ν/f_ν (and correspondingly dL_ν/L_ν) depends on the slope of f_ν. In the near-IR the slopes are assumed f_ν ∝ ν^1/3 and probably redder toward the optical and UV (e.g., Zheng et al. 1997; Vanden Berk et al. 2001; Kishimoto et al. 2008). This means that Δf_ν/f_ν in the UV is larger than in the optical, and because more energy is radiated in the UV, any light curve observed in the optical potentially underestimates the variability amplitude of the absorbed emission.³ On the other hand, in the process of radiative transfer the whole accretion disk spectrum is integrated over all frequencies, convolved with the dust absorption efficiency, a process that probably smoothens out the extreme variability amplitudes. We can introduce a constant factor w_x ≡ (dL/L)/(dL_x/L_x), the variability efficiency factor at waveband x, as a free parameter to compensate for these complications without requiring additional assumptions on the wavelength-dependence of the variability. In this way w_x will also include any effects that potentially alter the relative variability (e.g., non-variable light from the host), but are not included in this simple model.

Fig. 4 2.2   μm–8.5   μm model cross-correlation function (CCF). From dark blue to light blue the distribution changes from very extended (α = 0.0) to very compact (α =  −2.0).

An additional complication is a possible change in dust composition with distance from the AGN. Several studies found evidence that the inner torus is dominated by large graphite grains while the outer torus contains a mix of silicates and graphites. This leads to comparably small sublimation radii and a change of emissivity from the near-IR to mid-IR resulting in a “bump” at around 3 μm in many type 1 SEDs (e.g., Kishimoto et al. 2007; Mor et al. 2009; Kishimoto et al. 2011b). Therefore, when comparing observed light curves in the near- and mid-IR it may become necessary to adjust the absorption efficiencies in Eqs. (3) and (6) accordingly.

4. Modeling the near-IR light curve of NGC 4151

A comparison between near-IR and mid-IR light curves seems to be a promising tool to constrain the brightness profiles of AGN dust tori. As of yet, mid-IR light curves are very sparse. However, as described in Sect. 3.2.1, near-IR light curves have different peak intensities and tails for compact or extended distributions with respect to the AGN variability signal. We may therefore hope that comparing light curves in the optical and near-IR can also help decide if the brightness distribution in the torus is compact or extended, at least in the inner part of the torus.

One of the best-monitored AGN in the IR is the bright nearby Seyfert 1 galaxy NGC 4151. Koshida et al. (2009) presented optical and near-IR light curves with a very high temporal coverage, and the data presented below were extracted from their work. They were taken in the course of the MAGNUM project (Yoshi et al. 2003) over about 2000 days from 2001 to 2006. We use the reduced and host- and accretion-disk-subtracted data as shown in Fig. 1 of Koshida et al. (2009). As noted in Minezaki et al. (2004), the typical photometric errors for NGC 4151 data are 0.01 mag (~1% in flux; but see below for more details). Koshida et al. (2009) report that the K-band emission at 2.2 μm varies similarly to the V-band emission with a notable delay that reflects the light-travel time from the accretion disk to the torus. Near-IR interferometry of the same object showed that the K-band emission is probably coming from a confined region at the dust sublimation radius Kishimoto et al. (i.e. consistent with a thin ring; 2009b, 2011a). Therefore it seems that the object is well-suited to apply our model.

The AGN variability model has essentially two parameters: the brightness distribution power law index α (see Sect. 3.2) and the variability efficiency factor w_V (see Sect. 3.3). For the simulations we first determined the V-to-K time lag by Monte Carlo simulations of light curves based on the observed V- and K-fluxes, following the method outlined in Suganuma et al. (2006). For that we first calculated the structure functions $\mathit{s}\mathrm{(}{\mathit{t}}_{\mathit{i}}\mathrm{-}{\mathit{t}}_{\mathit{j}}\mathrm{)}\mathrm{=}{}^{\mathrm{(}}\sum _{\mathit{i}\mathit{<}\mathit{j}}\mathrm{\left[}\mathit{f}\mathrm{\right(}{\mathit{t}}_{\mathit{i}}\mathrm{)}\mathrm{-}\mathit{f}\mathrm{(}{\mathit{t}}_{\mathit{j}}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}{}^{\mathrm{)}}\mathit{/}\mathit{N}\mathrm{(}\mathit{i}\mathit{<}\mathit{j}\mathrm{)}$ of the V- and K-band light curves, where f(t_i) is the flux at epoch t_i and N(i < j) is the number of (i,j) pairs. From s we can determine the typical change of flux for a given time interval between observations. Based on s we performed Monte-Carlo simulations to interpolate the gaps in the observed light curves as follows: first we defined a regular t-sampling (T₁,T₂,...) that is finer than the observed one. Then linearly interpolated fluxes for this t array were calculated based on the observed fluxes. We pick a random epoch T_i from the pre-defined t-sampling, calculated the distance T_i − t_obs to the nearest observed epoch t_obs, and determined the corresponding flux change Δf_s(T_i − t_obs) from the structure function s. We then calculated a random flux change Δf from the linearly interpolated flux at T_i based on a Gaussian distribution with standard deviation Δf_s. The corresponding epoch and flux was added to the array of observed fluxes and the procedure repeated until fluxes were calculated for all T_i in the t-sampling. For a proper comparison we used the same t-sampling for the V- and K-band light curves.

Fig. 5 Observed and interpolated V-band light curve of NGC 4151. The blue filled circles are the observed V data points. Error bars are plotted but smaller than the symbols in most cases. Based on these observations, 50 interpolated light curves have been calculated and are shown as gray dotted lines, reflecting the uncertainty of the “true” light curve. The red dashed line indicates the mean of all these models, which has been used as the input signal to our model simulations.

The Monte Carlo simulations were repeated 50 times to obtain an idea of the uncertainties in the light curves introduced by the finite sampling of the observations. These simulated light curves can then be used to analyze the time lag between the V- and K-band emission. We found an average time lag measured over the full light curve, excluding the strongest (and widest) peak, of 43.8 ± 8.5 days, which corresponds to a (sublimation) radius of r_sub = 0.037 ± 0.007 pc, in agreement with near-IR interferometry (Kishimoto et al. 2009b; Pott et al. 2010; Kishimoto et al. 2011a). When using only individual features of the light curves, the error is generally larger but the resulting time lag is consistent with the average one. We will discuss consequences of this seemingly constant time lag below.

For the modeling we used the average time lag as an offset between the V and K light curves for the model simulations. In a next step the simulated V-band light curves have bee used to calculate a mean V-band light curve. The observed and simulated V-band light curves, as well as the mean light curve, is shown in Fig. 5. This mean V-band light curve was used as the input signal for our torus variability model.

Fig. 6 Observed and modeled K-band light curve of NGC 4151. The blue filled circles are the observed data points. Error bars are plotted but smaller than the symbols in most cases. Based on these observations, 50 interpolated light curves have been calculated and shown as gray dotted lines, reflecting the uncertainty of the “true” light curve. The model light curve is overplotted as a red dashed line.

In Figs. 6 and 7 we compare the observed K-band light curve (blue-filled circles) and the CCFs based on observations (gray-dotted lines), respectively, to the best-fitting model. The blue-filled circles in Fig. 6 are individual photometric observations extracted from Koshida et al. (2009). Based on the previously described Monte Carlo interpolation scheme, we simulated 50 (nearly continuous) light curves accounting for the gaps in the observed light curve and the resulting uncertainty. They are shown as gray-dotted lines, as light curves in Fig. 6 and as CCFs in Fig. 7. The range of these simulated curves reflects the range of curves if the temporal coverage were infinite. The simulated CCFs based on the observations have been normalized so that the peak value is unity. Overplotted (red-dashed lines) is the best-fitting K-band model in both plots. The formal fitting was made with the light-curve only. To asses the quality of the fit, we used two different versions of the error estimates of the observation. Although Koshida et al. (2009) do not provide typical errors, an initial fraction of the NGC 4151 K-band photometry using the same instrumentation was published in Minezaki et al. (2004) and a photometric error of about 0.01 mag (or about 1% in flux) was reported. When using this value for the full data set, we obtain a reduced ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\mathrm{=}\mathrm{4.2}$. An alternative approach to estimate the error in the data comes from the structure function that we used in the Monte Carlo interpolation of the light curves. The structure function of an AGN is supposed to follow a power law from short to long intervals t_i − t_j until a break at long intervals (e.g., Suganuma et al. 2006). At short time intervals the intrinsic variation should approach 0 for t_i − t_j → 0. In the K-band data we see, however, a flattening of the structure function for t_i − t_j ≲ 6 days, i.e. the (squared) variation of the fluxes becomes independent of the observed interval. This is more typical for uncertainty in measurements than intrinsic variability. When interpreted in this way, the measurement error is 3.5% and the best fit has ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\mathrm{=}\mathrm{1.2}$.

Despite the simplicity of our model, it reproduces the overall peaks and dips of the observations and interpolations (gray region) quite well. There are, however, some notable deviations, reflected in the moderate ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}$ when assuming 1% photometric error. In particular the observed light curve drops more steeply after the major outburst at 20 r_sub/c, leading to a lower general Δf(t)/f₀ at the peak around 33 r_sub/c although the shape of this peak is adequately reproduced. This general flux level is only recovered at the last peak after 40 r_sub/c.

The nominal best-fit model has a fairly compact brightness distribution with α =  −1.75. However, the probability distribution in α parameters is very broad, ranging from about  −0.5 to  −2.0 (and probably beyond) with very similar ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}$ values (assuming 1% photometric error), meaning that this parameter is poorly constrained. A combination of K-band with longer wavelength light curves as discussed in Sect. 3.2 is needed to better constrain the dust distribution of the hot-dust/sublimation zone emission. The efficiency factor w_V is much better constrained. In Fig. 8 we show the probability distribution for α and w_V for both 1% and 3.5% photometric error. In the range of acceptable α-values around the best fit we find ${\mathit{w}}_{\mathit{V}}\mathrm{\approx }\mathrm{0.}{\mathrm{4}}_{-0.1}^{\mathrm{+}\mathrm{0.2}}$, meaning that only about half of the power in the V-band variability is converted into variability in the K-band. Based on the discussion in Sect. 3.3 one may have expected that the V-band probably underestimates the amplitude of the (relative) variability of the integrated AGN emission and, therefore, leading to w_V > 1. However, a number of other factors can play a role in the value of w_V, e.g., uncertainties in host and accretion disk subtraction in either of the wavebands, or additional light close to the nucleus that does not participate in the variability, at least not on the covered time scales (see below for such a scenario). A more extended study using a sample of objects and/or additional UV/optical wavebands could help in better understanding the energy conversion from the accretion disk emission into the dust emission and its observational caveats.

Fig. 7 V/K-band cross-correlation functions of NGC 4151. Based on the observed V- and K-band light curves, 50 interpolated light curves for both V- and K-band have been simulated and CCFs calculated accordingly, shown as gray dotted lines. These reflect the observations and the uncertainty in the determination of the time lag owing to the gaps and non-uniform sampling of the data. The CCF based on the light-curve model is overplotted as a red dashed line. All CCFs have been normalized to their peak value.

Fig. 8 Probability distribution of our model grid for the radial brightness distribution power law index α and the V-band efficiency factor wV. The gray shaded regions with red dashed boundaries show the probability distribution assuming an error of the observations of 1% from Minezaki et al. (2004), leading to ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\mathrm{=}\mathrm{4.2}$. The blue dotted lines assume an error of 3.5% based on the analysis of the structure function with a ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\mathrm{=}\mathrm{1.2}$.

Interestingly, the observed change in flux in the V-band does not seem to have shown then naively expected change of the time lag. Koshida et al. (2009) already note that a potential change in time lag does not scale with (ΔL)^1/2. With a V-band flux change of about a factor of 20 from peak to valley⁴, we could have naively expected that the observed time lag/sublimation radius changed by a factor of 20^1/2 = 4.5 (see Eq. (1)). Koshida et al. (2009) conclude that the time lag changes from about 60–70 days in the first 2/3 of the light curve to about 35–50 days in the last 1/3, although at different scaling than (ΔL)^1/2. If such a change were significant, we should have noticed a shift in the peaks when comparing observed and modeled light curves because our models assume a constant time lag. Over about 600 days (~13.7   r_sub/c) the shift between observed and modeled light curve would be  ~6   r_sub/c. However, such a shift is not seen in Fig. 6. All the peaks and valleys in the modeled and observed light curve overlap within less than 1   r_sub/c. This would also have affected the comparison of observed and modeled CCF, but both are consistent within errors (see Fig. 7).

To test for even small effects of possible dust destruction and reformation, we removed the assumption of a constant time lag/sublimation radius and account for dust sublimation once it heats over the sublimation temperature (set as a free parameter) and instant or delayed reformation of the dust after cooling. This changes the sublimation radius and time lag as the variability progresses through the torus. But we found that the smallest ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}$ is always found for a model with constant sublimation radius/time lag, i.e. disfavoring significant dust destruction or reformation over the observed time span. This implies that the dust can either strongly overheat before it is destroyed, or that the dust is efficiently self-shielded. The survival of a dust grain depends on the balance between gas pressure and vapor pressure. If partial gas pressure dominates, then a dust grain is stable; otherwise it evaporates. The vapor pressure of dust, p_vap ∝ exp( − 1/T), is a strong function of the temperature T, while the partial pressure, p_gas ∝ T, depends only linearly on T. Therefore, we would expect that the lifetime of individual dust grains is short once they are heated over the sublimation temperature. As a consequence the observed behavior would favor shielding in a locally-dense environment instead of overheating of dust grains. This is consistent with the idea of a clumpy torus where the dust is confined in optically-thick clouds. A change in luminosity first acts on individual clouds, which may sublimate part of their dust content but can resist longer overall at the same location than smoothly distributed dust that is not shielded locally. This scenario may also explain the low w_V value: If individual clouds close to the sublimation radius are heated up to the sublimation temperature (and above in corresponding equilibrium temperature) and their dust becomes only gradually sublimated from the surface (e.g., as a “melting snowball”), their actual peak temperature would remain essentially constant, leading to a more or less constant K-band flux over some time. Hence only cooler clouds would contribute to the variability on the same time scales as the incident radiation.

5. Summary and conclusions

We presented model simulations of time-variable infrared emission from dust as a consequence of variability of the incident radiation. For that we first introduce a generalized treatment for temperature variations, which can be used for all kinds of dusty environments. We applied this scheme to a simplified model of a (clumpy) dusty torus around AGN and investigated how the variability of the accretion disk radiation influences the torus emission in the near- and mid-IR. The main parameter of this model is the radial brightness distribution of the torus that has previously been shown to be connected to the radial distribution of the dust in the torus. We showed that any variability signal in the optical is smoothened stronger if the brightness distribution is very extended. While this effect is true for both the near- and mid-IR, longer wavelengths show much wider transfer functions than short wavelengths. The time lags between the optical and near-IR emission is mostly representing the light travel time from accretion disk to the sublimation radius independent of the brightness distribution. For mid-IR wavelengths, however, time lags can become very long, up to 10s of r_sub/c. The effect that the brightness distribution influences the time lags seems to be much stronger than any similar effect from inclination (at least in type 1 AGN) or details of the shape of the inner torus, as recently presented by Kawaguchi & Mori (2011). This change of lag time from near-IR to mid-IR can be used to quantify the brightness distribution of the torus, either by comparing optical light curves to near-IR and mid-IR light curves, or by directly comparing near-IR to mid-IR light curves.

Finally, the model has been applied to the optical and near-IR light curves of the nearby Seyfert 1 galaxy NGC 4151. We showed that the simple model can reproduce the overall observed variability signal with some deviations in the details. For an even better match it may be necessary to use the variability scheme in the framework of more complex torus model, e.g. by incorporating it into CAT3D (Hönig & Kishimoto 2010). Nevertheless, we conclude from our modeling that a single-wavelength near-IR light curve is probably not enough to constrain the brightness distribution in the torus, requiring at least one other light curve at longer wavelengths as discussed before. We found, however, that about 40% of the energy in the variability signal in the V-band has been converted into K-band variability. This

may be explained by a scenario where individual clouds close to the sublimation radius gradually sublimate their dust from the surface inward (“melting snowball”), essentially keeping their temperature constant, and, therefore, do not contribute significantly to the variability at the measured time scales. We also note that our modeling does not support a change of time lag/sublimation radius over the observed light curve epoch in spite of a significant change in V-band emission.

Acknowledgments

S.F.H. acknowledges support by Deutsche Forschungsgemeinschaft (DFG) in the framework of a research fellowship (“Auslandsstipendium”). The authors are grateful for valuable comments and suggestions by the anonymous referee.

References

All Figures

	Fig. 1 Dependence of the relative temperature change dT/T on the actual dust temperature for a variability in luminosity of dL/L = 1. The black solid line shows standard ISM (see text for details) while the blue dashed line is for large graphite grains. The red dotted line indicates the black-body limit.
In the text

	Fig. 2 Model light curves at 2.2   μm (red) and 8.5   μm (blue) for different values of the radial brightness distribution power law index α. From dark to light the distribution changes from very extended to very compact. The black dashed line illustrates the variability signal from the accretion disk.
In the text

	Fig. 3 Model cross-correlation function (CCF) for light curves of AGN variability at 2.2   μm (red) and 8.5   μm (blue). From dark to light the distribution changes from very extended (α = 0.0) to very compact (α =  −2.0). The dashed line marks a lag time of rsub/c.
In the text

	Fig. 4 2.2   μm–8.5   μm model cross-correlation function (CCF). From dark blue to light blue the distribution changes from very extended (α = 0.0) to very compact (α =  −2.0).
In the text

Fig. 5 Observed and interpolated V-band light curve of NGC 4151. The blue filled circles are the observed V data points. Error bars are plotted but smaller than the symbols in most cases. Based on these observations, 50 interpolated light curves have been calculated and are shown as gray dotted lines, reflecting the uncertainty of the “true” light curve. The red dashed line indicates the mean of all these models, which has been used as the input signal to our model simulations.

In the text

Fig. 6 Observed and modeled K-band light curve of NGC 4151. The blue filled circles are the observed data points. Error bars are plotted but smaller than the symbols in most cases. Based on these observations, 50 interpolated light curves have been calculated and shown as gray dotted lines, reflecting the uncertainty of the “true” light curve. The model light curve is overplotted as a red dashed line.

In the text

Fig. 7 V/K-band cross-correlation functions of NGC 4151. Based on the observed V- and K-band light curves, 50 interpolated light curves for both V- and K-band have been simulated and CCFs calculated accordingly, shown as gray dotted lines. These reflect the observations and the uncertainty in the determination of the time lag owing to the gaps and non-uniform sampling of the data. The CCF based on the light-curve model is overplotted as a red dashed line. All CCFs have been normalized to their peak value.

In the text

Fig. 8 Probability distribution of our model grid for the radial brightness distribution power law index α and the V-band efficiency factor wV. The gray shaded regions with red dashed boundaries show the probability distribution assuming an error of the observations of 1% from Minezaki et al. (2004), leading to ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\mathrm{=}\mathrm{4.2}$. The blue dotted lines assume an error of 3.5% based on the analysis of the structure function with a ${\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\mathrm{=}\mathrm{1.2}$.

In the text