Free Access
Issue
A&A
Volume 578, June 2015
Article Number A57
Number of page(s) 18
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/201525733
Published online 04 June 2015

© ESO, 2015

1. Introduction

To explain the apparent dichotomy between Type 1 and Type 2 active galactic nuclei (AGNs), the unified AGN model postulates a so-called dust torus (Antonucci 1993; Urry & Padovani 1995) around the broad-line region (BLR) of AGNs. Depending on viewing angle, this optically thick structure will, or will not, obscure the central region, resulting in different spectral energy distributions (SEDs) for Type 1 and Type 2 AGNs. Direct and indirect observational evidence (Barvainis 1987; Heisler et al. 1997; Osterbrock & Ferland 2006; Jaffe et al. 2004; Elitzur 2006) has been accumulated over past decades, confirming the existence of the assumed dusty torus around the BLR. Nevertheless, our knowledge of the long-term stability of the dust around AGNs and of its origin remains vague.

Self-consistent AGN torus models (Krolik & Begelman 1988; Schartmann et al. 2010) assume that dust is ejected into the interstellar medium by asymptotic giant branch (AGB) stars and brought from outside to the central nuclear region. Unless the luminosity of the central source has decreased on a timescale shorter than the dynamical timescale, the innermost dust for these models is always expected at the sublimation radius. A different scenario, which suggests a location or formation of dust farther out than the sublimation radius, has been proposed by Elvis et al. (2002). Here, the dust formation results from BLR clouds that are assumed to be part of an outflowing accretion disk (AD) wind1. In the course of a subsequent expansion and cooling process, the conditions of dust condensation are fulfilled and dust is formed – by the AGN itself – at a distance farther out than the sublimation radius.

To be able to distinguish between competing dust formation scenarios, it seems essential to determine the stability and evolution of the innermost dust radius. Monitoring the hot dust temperature is a further indicator of the immediate dust state, i.e. whether the hot dust is at the sublimation temperature. Observations, even if performed on one and the same AGN, NGC 4151, have so far given seemingly inconclusive indications:

  • VK reverberation measurements indicate strongly varying time lags between 30 and 70 days over the years 2001–2006 (Koshida et al. 2009), suggesting a strong variation in the location of the innermost dust over the observed period, possibly by means of dust destruction and fast re-formation. A reanalysis of these data by Hönig & Kishimoto (2011) with a simplified clumpy torus model resulted in an apparently constant time lag of τ = 43.8 ± 8.5 days for that same period.

  • Interferometric K-band observations in 2010 (Pott et al. 2010) have shown that the hot circumnuclear dust torus around NGC 4151 apparently does not expand owing to dust sublimation when the AD brightens, but found indications of dust survival in bright times.

  • From our 2010 z to K photometry of NGC 4151, we found no signatures of dust destruction following increased AD emission (Schnülle et al. 2013, hereafter Paper I). Instead, we measured a significant increase in the hot dust temperature, indicating that the hot dust in NGC 4151 was apparently outside the current sublimation radius in that period.

  • Recently, with K-band interferometric observations (Oct. 2010–May 2012), Kishimoto et al. (2013) have found indications of an expanding K-band dust location for NGC4151 caused by increased AD emission and a correlation of the hot dust radius with the averaged AD flux over the past six years.

While generally the dust location around AGNs is known to correlate with the average AD luminosity as (Minezaki et al. 2006; Kishimoto et al. 2011), the above results imply that it is not a tight function of the instantaneous luminosity state.

We present here the updated results from our monitoring program of NGC 4151 during 2010–2014, in which we used dust reverberation to study the evolution of the temperature and reverberation lag of the hot dust around this archetypical Sy1 AGN. Compared to Paper I, we substantially extended the data set by using a significantly longer time series and a broader wavelength coverage. We also extended our reduction and analysis techniques by performing precise photometry with an image subtraction method and employing a multiepoch, multiwavelength Markov chain Monte Carlo (MCMC) fit. NGC 4151 is part of our AGN hot dust reverberation campaign, in which we monitor roughly 25 nearby, bright, and variable Type 1 AGNs from the optical to the near-infrared (NIR).

Observations, data reduction and photometry are described in Sect. 2, followed by an introduction of our data analysis methods in Sect. 3. We present our updated results on NGC 4151 in Sect. 4, and discuss these in Sect. 5. A summary of our results, as well as the outlook for this project, is given in Sect. 6.

2. Observations and data processing

2.1. Near-infrared data

From 2010 January–June, we obtained NIR multiband photometric observations of the Sy1 NGC 4151 (see Paper I), using the Omega 2000 (Kovács et al. 2004) NIR wide-field camera mounted on the 3.5m telescope in Calar Alto, Spain2. Omega 2000 has a field of view , and its 2048 × 2048 pixel detector has a scale of pixel-1.

We continued our monitoring of NGC 4151 from 2012 February–2014 June, with roughly two to four weeks of sampling (apart from data gaps in 2010 July–2012 February and 2012 May–November). In each of the 29 epochs, we obtained broadband photometry of NGC 4151 and of three calibration stars in the same field of view in the z, Y, J, H, and K spectral bands. These five filters were chosen to differentiate between variations in the AD and hot dust emission. The weather conditions ranged from clear to photometric, and the seeing was .

Following the data reduction with IRAF (see Paper I), we extracted the nuclear fluxes of NGC 4151 at various epochs and filters, using two different approaches. First, we performed a PSF-bulge-disk decomposition with GALFIT (Peng et al. 2002), which allows for the simultaneous fitting of different galaxy components to the objects in an image. For NGC 4151, the fit model consisted of the point spread function (PSF) – constructed with the help of the IRAF DAOPHOT package for each frame – plus a Sersic profile and an exponential disk profile. This strategy led to small statistical errors in the fit (typically 0.01 mag for the PSF and Sersic components, and 0.02 mag for the disk component) as well as sufficiently smooth residuals for our purpose. The Sersic half-light radii Re, Sersic indices n, and disk scale lengths h for the different filters and epochs were found to be temporally stable. In a final fit, Re, n, and h were fixed to their average values determined from the previous fits, i.e. , n = 3.0, . These parameter values for NGC 4151 are consistent with the literature, for example Bentz et al. (2009). The residuals deviate significantly from zero only in the PSF regime, while they are very smooth on larger scales. Substantial residuals in the PSF domain are likely to be caused by the complex PSF (due to defocusing to avoid saturation, plus seeing), which is variable over the detector. Aperture photometry on these residuals yielded low fluxes, on the order of 2–3% of the PSF flux of the AGN. The fluxes of the reference stars were determined by performing a simultaneous PSF model fit in GALFIT.

Absolute flux was calibrated in the JHK bands with the known fluxes of our calibrators, as derived from the 2MASS PSC (Skrutskie et al. 2006). For the zY photometry, we used a database containing model spectra of main sequence stars (kindly provided by van Boekel, priv. com.). For each calibrator, we determined the best fit spectrum to the BVJHK fluxes as given by SIMBAD and calculated the zY photometry from that spectrum.

Table 1

NIR fluxes of the nucleus of NGC 4151, derived with ISIS.

As an alternative to GALFIT, we used the ISIS image subtraction package (Alard & Lupton 1998; Alard 2000), which enables precise measurement of flux variations in variable objects. Following the procedures described by Shappee & Stanek (2011), the images in each band are first aligned using the program Sexterp (Siverd et al. 2012). Then a reference image has to be chosen or constructed – in our case, we chose the image with the best seeing. A spatially variable convolution kernel K is determined that minimizes the discrepancy (1)i.e. the kernel optimally transforms the reference image R to the seeing and flux level of a given individual image I. Each frame is then subtracted from the convolved reference image. Light curves, which are difference fluxes relative to the reference image, are then extracted from the subtracted images by performing aperture photometry on the residual nuclear flux. The error in the difference flux is estimated by measuring the residual flux of non-variable stars in the field (which should ideally be zero). Absolute calibration of the difference flux of each individual frame was performed with the same three calibration stars as in our GALFIT analysis. The absolute nuclear flux level in the reference image was also adopted from our GALFIT results. As can be seen in Fig. 1 for the K band, our GALFIT and ISIS fluxes agree well within the given errors. For targets at substantially farther distances than NGC 4151 (as are present in our broad AGN sample), a GALFIT decomposition might be more difficult owing to increased degeneracies between PSF and the bulge, making ISIS the method of choice for these objects. For this reason and overall consistency, we use the ISIS photometry (listed in Table 1) throughout this paper.

thumbnail Fig. 1

Nuclear fluxes of NGC 4151 derived with ISIS (green) and GALFIT (magenta), shown here for the K band. Within the photometric errors, we observe excellent agreement between the two methods. The agreement is on a similar level for all other bands, suggesting that both methods manage to separate the nuclear flux of interest from the bulge and disk.

Open with DEXTER

To remove emission line contributions in our zYJHK fluxes, a NIR spectrum of NGC 4151 was obtained with the SpeX spectrograph (Rayner et al. 2003) at the NASA Infrared Telescope Facility (IRTF) on 2010 February 26, as already described in Paper I. For each filter, we applied a correction factor to our data (8% in zY, 15% in J, 2% in HK), derived from the ratio of the synthetic IRTF photometry to the continuum flux level at the respective effective wavelength.

2.2. Mid-infrared data

Table 2

WISE W1, W2, W3 photometric fluxes of NGC 4151, derived with GALFIT.

Furthermore, we used mid-infrared (MIR) images of NGC 4151 in the WISE W1, W2, and W3 bands, downloaded from the NASA/IPAC Infrared Science Archive (IRSA), and extracted the nuclear fluxes with GALFIT. The downloaded images in each band are a stack of 31 individual frames, taken on 2010 May 31 and 2010 December 8. Although images are only available in stacked form, the total magnitude per band is given in the ALLWISE multiepoch catalog table for each of the 31 individual observations. From a weighted magnitude difference between 2010 May and 2010 December (that we attribute to the variability of the nucleus), we derived a correction term to scale our nuclear fluxes determined with GALFIT to the epoch 2010 May 31 (MJD 55 347, then almost coinciding with Epoch 5 of our NIR observations). Our derived AGN and host fluxes are listed in Table 2. The total flux in each band agrees very well with the average total flux for NGC 4151 given in the ALLWISE catalog.

2.3. Optical data

Furthermore, optical spectra were obtained between 2012 January and April as part of a reverberation mapping campaign described elsewhere (De Rosa et al., in prep.). Observations were made on the 1.3-m McGraw-Hill Telescope at the MDM Observatory on Kitt Peak with the Boller & Chivens CCDS spectrograph. A 350 lines mm-1 grating yielded a dispersion of 1.33 Å pixel-1. The entrance slit oriented north–south (PA = 0o) with a projected width of . This configuration yielded spectra covering the range 4400–5850 Å with a spectral resolution of 7.9 Å. We used an extraction window of in the cross-dispersion direction.

Spectra were scaled to a common [O iii]λ4959 flux of 3.76 × 10-12 erg s-1cm-2. The starlight contribution to each spectrum was estimated from host-galaxy surface brightness models based on Hubble Space Telescope images as 18.4( ± 1.8) × 10-15 erg s-1 cm-2 Å-1 (Bentz et al. 2013).

The continuum fluxes at 5100 Å (listed in Table 3) are then derived by averaging the spectrum over a 20 Å range around the lowest point, which is located between the [O iii]λ49595007 lines and the Fe ii blends.

Table 3

Optical continuum fluxes of the nucleus of NGC 4151 at 5100 Å.

3. Methods

3.1. Dust reverberation

The hot dust around the central nuclear region absorbs radiation emitted by the AD and re-emits it in the infrared. This dust emission can be roughly approximated by a blackbody function. Owing to the absorption and re-emission of incident radiation by the dust, variations in the AD flux can also be detected in the dust emission with a characteristic lag time of τ = Rdust/c, where Rdust is the radial distance of the hot dust from the AD and c the speed of light. With dust reverberation, that is to say, by monitoring the hot dust signal as it responds to the AD signal, we determine the reverberation lag time and the temperature evolution of the hot dust around NGC 4151, to infer the stability of the hot dust and to conclude whether it is close to sublimation.

We note that the dust sublimation temperature depends strongly on the assumed dust species, composition and grain sizes, all of which might change significantly with radial distance from the source. Without detailed knowledge of these dust properties, it is therefore not possible to determine whether the innermost dust is close to sublimation from measuring the dust temperature in one single epoch. Monitoring the hot dust temperature for an enhanced period, however, allows us to draw more reliable conclusions about the immediate dust state. If we observe rising temperatures following increased AD emission, we can conclude that the hot dust has not yet reached its sublimation temperature. Clearly, the temperature evolution is a more robust indicator than the absolute dust temperature in one single epoch.

To be able to measure the dust response to varying AD emission, we decompose the total fluxes in each epoch and band as the sum of an AD and a blackbody component that we attribute to the dust. In contrast to the single-epoch SED fits (Eq. (5)) used in Paper I, we perform a multiepoch, multiwavelength fit on our complete data set, of the form (2)Here, x denotes our complete vector of model parameters. The first term on the righthand side of Eq. (2) refers to the AD emission, usually described as a power law in terms of λ, with a power-law index α which we assume does not change significantly over the observed flux range. For completeness, however, we also test for a varying α and its influence on our results, see Sects. 4.2 and 5.3. C1 is a proportionality constant. In the second term on the righthand side of Eq. (2), BB(.) denotes the blackbody function (3)with T the blackbody temperature of the hot dust and C2 the blackbody constant (consisting of a mixture of emissivity, solid angle, and surface filling factor). In Eq. (2), we assume that the AD signal has the shape of the interpolated z-band signal, which is justified because the z band is known to be dominated by AD emission (Riffel et al. 2009). Starting from an initial temperature T0, the blackbody temperature T evolves by construction, for t>τ, in response to the AD variations as (4)The variability factor v accounts for the possibility that the z-band variability is not completely reprocessed by the hot dust (Hönig & Kishimoto 2011), which would lead to values of v< 1. Furthermore, values of v ≠ 1, might indicate that the variability of the heating signal is underestimated or overestimated by the z-band variability. The hot dust around AGN is heated mainly by the UV part of the AD radiation. We take the z-band flux as a proxy for this UV heating; however, the actual UV radiation is probably even more variable. Our full parameter vector thus comprises x = (C1,α,C2,T0,τ,v). Our new multiepoch approach Eq. (2)–(4), in particular the temperature evolution given by Eq. (4), is justified by the measured hot dust temperature evolution obtained from preceding single-epoch SED fits of the form (5)i.e. AD plus blackbody contribution. These single-epoch fits were used in Paper I, and were also performed for each epoch of the complete 2010–2014 data set, prior to the multiepoch approach described by Eqs. (2)–(4). As can be seen in Fig. 7, the blackbody temperature obtained from the single-epoch SED fits tracks the AD signal closely, with a stable delay of roughly 30 days, and the temperature variation is significant, thus justifying our new model. Furthermore, behavior as in Eq. (4) is expected according to Hönig & Kishimoto (2011). Compared to Eq. (5), the data are temporally connected in our new approach. Using the temperature evolution described by Eq. (4) drastically reduces the number of parameters of the multiepoch fit, since we only need to estimate the initial blackbody temperature T0 instead of a temperature for each single epoch. This makes the algorithm more robust against degeneracies between T and C2.

To check whether the distribution of the innermost hot dust is fairly compact or, in contrast, substantially radially extended (as we conjectured in Paper I), we tested for the possibility of a further blackbody component contributing to the hot dust signal. Therefore, in addition to the 1BB-model of Eq. (2), we fitted a 2BB-model of similar type to our data: (6)with the temperature evolution of the second blackbody given by an expression as in Eq. (4), though with different starting temperature, time lag, and variability factor. This 2BB-model has four additional parameters, and x = (C1,α,C2,T0,1,τ1,v1,C3,T0,2,τ2,v2).

3.2. Interpolation of the AD signal

To calculate the blackbody temperature T described by Eq. (4) of our new approach for any arbitrary value of τ in the range of the data, we first need to interpolate our input accretion disk signal, i.e. the z-band signal. This is done with the method for interpolation, realization, and reconstruction of noisy, irregularly sampled data by Rybicki & Press (1992) and Press et al. (1992), which was previously used for reverberation studies of the BLR, for example by Zu et al. (2011) or Hernitschek et al. (2015). The AD signal can be described as a stochastic process characterized by two structure function parameters, which are estimated from the data. With this we can predict the signal for unmeasured times.

Our M = 29 measurements y = { yi } ,i = 1,...,M of the z-band flux can be written as (7)where s represents the intrinsic variability signal, E is the vector (1,1,1...,1)T, is an appropriate estimator of the mean of the data, and n represents the measurement noise. The intrinsic signal s can be described as a stochastic process with a certain correlation structure, and we can calculate an estimate of the signal at any (measured or unmeasured) point using the already obtained M measurements: (8)Here, are linear coefficients that depend on the particular point to be estimated, and x is the discrepancy between the estimated and the true value. Minimizing the discrepancy with respect to the linear coefficients then yields the solution for the coefficients, and the least variance estimate of our signal (9)Here, S is the covariance matrix of the data, S is the covariance vector between the measured data points and the new data point, and N is the noise covariance matrix. In Eqs. (8) and (9), the mean is subtracted from the data before determining the coefficients3 and corresponds to the value that minimizes and is added back to the estimate afterward. The variance of the true value s about the best estimate is given by (10)To apply Eq. (9), we need to estimate the covariance structure of the underlying stochastic process from our data. Assuming stationarity, so that SijS(titj) ≡ S(τ), the covariance matrix can be replaced by the population mean square and the structure function V(τ): (11)Following the method described by Press et al. (1992), we calculate a lag τij = | titj | and an estimate of the structure function for each pair (i,j), where ni is the measurement error of data point yi. The M(M − 1)/2 = 406 pairs are binned into 20 bins, equally spaced in τ. In Fig. 3, is plotted against τij on a logarithmic scale. We fit a straight line (12)to these data, representing a power-law model of the form (13)where A is the average flux variability on a one-year timescale, and γ is the gradient of this variability. For the structure function parameters, we obtain best fit values of A = 0.011 ± 0.001 Jy (resp. A = 0.423 ± 0.019 W m-2μm-1) and γ = 0.881 ± 0.067. With these values, we calculate the covariance matrix and vector, and interpolate our z-band light curve according to Eq. (9). The interpolated z-band light curve is shown in Fig. 4. Since our measurement errors are generally small (at the 5% level), the interpolated light curve runs almost perfectly through our data points. In Fig. 3, a potential substructure seems to be evident around τ = 2.5 yr. However, this apparent feature might be merely an artifact, caused by insufficient sampling in that range of τ. Indeed, the bins around τ = 2.5 yr contain by far the least number of data points, obviously owing to the large data gap between June 2010 and February 2012. An apparent feature seen in Fig. 3 seems to be common in the empirical structure functions of single targets (see e.g. figures in Press et al. 1992; Schmidt et al. 2010; Morganson et al. 2014), while it is averaged out in the observed structure function of samples of AGNs (see e.g. Schmidt et al. 2010).

Table 4

Upper and lower limits a,b for uniform pdfs of the model parameters, for the 1BB model (with constant and varying α) and the 2BB model.

Our derived structure function parameters agree with values found in the literature. While for quasars, typical values of the power-law slope are γ ≈ 0.3–0.4 (Bauer et al. 2009; Schmidt et al. 2010), the structure function of Seyferts is generally found to be steeper (Hawkins 2002), except for very short timescales, which are dominated by noise. Indeed, specifically for NGC 4151, Czerny et al. (2003) report slopes of 0.65–1.0 for the V band structure function on different timescales and in different epochs.

An alternative approach of the structure function to describe AGN variability is the damped random walk (DRW) model: (14)where σ2 is the long-term variance of the process, and τd the damping timescale. While for a substantial τ range the variability increases with increasing time lag following a power law V(τ) ∝ τγ, the DRW model considers that there is a plateau where the variability saturates, at time lags that are longer than the longest correlation timescale of the stochastic process. When we fit this model to our observed structure function for the NGC 4151 z band signal, we obtain values σ = 0.148 ± 0.505 Jy and τd = 230.0 ± 576.1 yr. Clearly, the fitted damping timescale, typically between 0.1 and 3 years (MacLeod et al. 2010), seems physically unreasonable. When looking at the observed structure function data plotted in Fig. 3, we can see that our time series is obviously not long enough to sufficiently constrain the damping timescale for NGC 4151, as also evident from its large error. Obviously, no plateau is reached, but the variability increases throughout the observed time range. This is consistent with the results of Czerny et al. (2003), who find an unusually long damping timescale of roughly ten years for this object. We performed the fit again with fixed τd = 10 yr, leading to σ = 0.032 ± 0.002 Jy.

The obtained structure functions for the power law and the DRW model are plotted in Fig. 3. To test the influence of the particular choice of structure function on the robustness of our results, we performed all our following analyses (i.e., the interpolation of the AD signal as well as the various fits presented in Sects. 4.2 and 5) with all three models shown in Fig. 3.

3.3. Model fitting using a DE-MC algorithm

3.3.1. Differential evolution Markov chain

Compared to Paper I, in which we used the routine mpfit to fit the model of interest to our data, we now infer the best model parameters using a differential evolution Markov chain (DE-MC) algorithm, which was proposed by Ter Braak (2006). DE-MC is a population MCMC algorithm, with multiple chains running in parallel. Compared to random walk metropolis (RWM) type MCMC algorithms, the main advantage of DE-MC is that the problem of choosing a jumping distribution, hence an appropriate scale and direction of the jumps, is solved by generating jumps as a fixed multiple of two randomly selected parameter vectors xR1 and xR2 that are currently in the population: (15)where γ> 0, and the vector e of small random numbers is added to ensure that the whole parameter space can be reached. These random numbers are drawn from a symmetric (e.g. Gaussian) distribution that has small variance compared to the target variance. In the proposal scheme of Eq. (15), the chains learn from one another, hence the name “differential evolution”.

3.3.2. Implementation of the algorithm

We ran the DE-MC algorithm for the six-parameter 1BB model of Eq. (2), with 15 simultaneous Markov chains and 25 simultaneous chains for the case of the ten-parameter 2BB model of Eq. (6). We found to be a good value for achieving the desired acceptance rate of roughly 25%. In each iteration step and for each chain, the algorithm evaluates – for the current parameter vector xi and for the proposal vector xp – the posterior probability density function (pdf) (16)which is proportional to the prior pdf times the likelihood function according to Bayes’ theorem. Here, x represents the parameter vector, and y are the data. The notations mk and yk respectively refer to the values of the model and the data for the kth of the N = 211 NIR plus optical data points (the WISE data were not included in the model fit), and σk is the related photometric error. To not overweight the low flux epochs around MJD 56 000, where all of the 5100 Å measurements are assembled, we decided to use a weighted . Thus, we obtain an effective number of N = 160 data points.

For each parameter x, we chose a uniform prior pdf (17)within the limits a,b as given in Table 4, in order to exclude unphysical parameter ranges. The posterior pdf is then directly proportional to the likelihood, and the maximum of the posterior occurs when χ2 is smallest. The upper limit for C1 is given by the fact that for the z band (λ ≈ 0.9 μm), C1·λα must not exceed 1. In the 2BB model, the maximum temperature for the second blackbody component is given by , corresponding to the case that all radiation from the AD can reach the second blackbody, i.e. the covering factor of the first blackbody is negligible. We chose τ2 = 100 days as an upper limit for the reverberation lag of the second blackbody component, because at that distance (and resulting low temperatures), its contribution to the NIR fluxes would become negligible, unless the blackbody constant C3 would be several orders of magnitude higher than C2. We further set p(τ2) = 0 for τ2<τ1.

As initial distribution for the parameters, we took a uniform distribution within the limits of our priors, which is thus sufficiently overdispersed.

For monitoring convergence of the DE-MC algorithm, we used the convergence criteria proposed by Gelman & Rubin (1992) and Brooks & Gelman (1998), i.e., the value, and visual inspection of W and , as explained in Sect. A. The 1BB model converged within 10 000–15 000 iterations on average. We considered as close enough to 1, Ter Braak (2006) uses . Figure 5 shows the evolution of the parameter T0,1 (which took longest to converge) for one exemplary run, as well as the evolution of , and W for this parameter. Visual inspection shows that W and are not evolving anymore and convergence was diagnosed correctly. The alternative ten-parameter 2BB model took 50 000–100 000 iterations on average to converge.

thumbnail Fig. 2

Calibrated 2010–2013 photometry of NGC 4151. The optical data are plotted in orange for the epochs that coincide with our NIR epochs, and in black otherwise. Photometry of the NIR data was done with ISIS, and is highly consistent with our previously determined fluxes with GALFIT (also see Fig. 1).

Open with DEXTER

thumbnail Fig. 3

Estimated (vij)1/2 (red dots) from the z band data versus τij, shown on logarithmic scale. The three different structure function models that we fit to these data (see text of Sect. 3.2) are represented by the green, magenta, and blue lines.

Open with DEXTER

4. Results

4.1. Constant power-law index

All results presented in this section were obtained using the power-law structure function model and a constant AD power-law index α. For results under the use of different structure function and α models, see Sects. 4.2 and 5.3. The temporal evolution of the photometry of the nucleus of NGC 4151 is shown in Fig. 2. The flux variations in the other bands show a time lag behind those of the AD-dominated z band – except for the 5100 Å fluxes, which seem to be concurrent with the z-band fluxes – and this lag increases with wavelength. This band-dependent behavior can be explained by the reverberation delay of the hot dust and by the increasing dust contribution and decreasing AD contribution for longer wavelengths.

thumbnail Fig. 4

z-band light curve, interpolated with the method of Rybicki & Press (1992) and Press et al. (1992). We interpolated using the three models described in Sect. 3.2. For clarity, only the interpolation resulting from the power-law structure function is shown here. This interpolated signal is the least-variance estimate of the stochastic process and therefore very smooth. The error of this prediction is indicated by the shaded gray region. Single realizations of the process have much more structure on short timescales and will also make excursions outside of the estimated error regions.

Open with DEXTER

thumbnail Fig. 5

Evolution of the parameter T0,1 (for one chain and the population mean) of the 1BB model, for one exemplary run. Also shown are the evolution of the pooled within chain variance W (Eq. (A.4)), maximum variance estimate (see Eq. (A.1)) and the resulting convergence parameter, i.e., the PSRF (Eq. (A.1)). Convergence was diagnosed at iteration step n = 13 599, where was reached for all parameters.

Open with DEXTER

thumbnail Fig. 6

Our data and the 1BB model fit for the epochs 3,5,9, and 10. Overplotted is the AD model, the 1BB model, and the sum of both components. From the data, we observe a clear rise in the blackbody flux, and a clear shift of the blackbody emission peak to shorter wavelengths until epoch 5. This emerging bump is confirmed by the blackbody temperature maximum in epoch 5 of our single-epoch fit (see Fig. 7).

Open with DEXTER

Figure 7 shows a single-epoch decomposition of the fluxes into AD and hot dust according to Eq. (5). The hot dust temperature closely follows the overplotted AD z-band flux with a delay of 30 days, which justifies the model described by Eqs. (2)(4). The minimum that is seen in the temperature around MJD 56700 (and also in the corresponding JHK flux minima) is not covered within the sampling of the z-band data (see Fig. 8 and Fig. 2), therefore we observe a substantial discrepancy here. Comparing epochs 1–9 with epochs 10–29, one can see that the ratio FAD,z/T is higher in the first epochs, already indicating our later finding that in these early epochs the reverberation lag might be larger than in the later epochs, so that the AD flux is more diluted and the dust heated less in epochs 1–9.

The results of our multiepoch, multiwavelength fit is shown in a temporal plot in Fig. 8 for the 1BB model and in Fig. 9 for the 2BB model. For the 1BB model, the best-fit lag is τ = 31.0 ± 1.6 days, and the best-fit value for the power-law slope is α = 1.63 ± 0.04 (cf. Table 5 for values of the other parameters), nicely matching the Fνν1/3 law expected from a standard Shakura-Sunyaev accretion disk. The initial blackbody temperature has a best-fit value of T0 = 1436 ± 15 K, and then evolves (by construction) as in Eq. (4). Epochs 1 and 2 are not included in the evaluation of the fit, as in our model, T evolves only for tτ. We find for this fit, indicating that we have probably slightly underestimated our measurement errors or that the model is too simple.

For the 2BB model, we find best-fit values of τ1 = 29.3 ± 1.9 days, T0,1 = 1479 ± 27 K, τ2 = 67.2 ± 8.3 days, T0,2 = 698 ± 64 K, α = 1.65 ± 0.04, and a reduced χ2 value of . Here, epochs 1–3 are not included into the fit, because T2 evolves only for tτ2. The best-fit values for the complete set of parameters are given in Table 5. For both the 1BB and the 2BB model, the data can be fit well with a stable time lag.

In Fig. 6, we show the results from our 1BB fit in the spectral domain, for the epochs 3, 5, 9, and 10. It is clear how the observed SED changes from epochs 3 to 5, and we can see a hot dust bump emerging, with the peak of the blackbody emission shifted to lower wavelengths (λ ≈ 1.9 μm), thus indicating our detected temperature increase. Though degeneracies between the blackbody constant and the blackbody temperature are non-negligible (in the single-epoch fits, where the temperature can evolve freely, the correlation coefficient is as high as − 0.75), visual inspection of the data and the systematic changes in the NIR color, correlating with the delayed AD brightness, underline the actual temperature increase. The hot dust peak is shifted to much longer wavelengths (λ ≈ 2.5 μm) in epoch 9, and until epoch 10, we again observe rising temperatures.

thumbnail Fig. 7

Hot dust temperature evolution derived from the single-epoch decomposition. The hot-dust temperature changes follow the relative AD z-band flux changes with a delay of roughly 30 days.

Open with DEXTER

thumbnail Fig. 8

Results of our 1BB model fit. Plotted are the K,H,J,Y,z, and 5100 Å data (green) over time. The bigger dots with errors bars in the 5100 Å panel mark the data points that coincide with our NIR epochs. The red line represents the blackbody contribution for each band, blue the AD contribution, and the black line is the sum of both. Epochs 1 and 2 are not included in the evaluation of the fit (as t<τ).

Open with DEXTER

thumbnail Fig. 9

Results of our 2BB model fit. Plotted are the K,H,J,Y,z, and 5100 Å data (green) over time. The bigger dots with error bars in the 5100 Å panel denote the data points that coincide with our NIR epochs. The red line represents the contribution of the inner blackbody component (at 29 light-day distance) for each band, the brown line is the second blackbody (at 67 light-days), blue the AD contribution, and the black line is the sum of all three components. Epochs 1, 2, and 3 are not included in the evaluation of the fit (as t<τ2). The resulting 2BB model fit only differs significantly from the 1BB model fit in the high-flux epochs 4-6, in the K band, as indicated by a red circle.

Open with DEXTER

thumbnail Fig. 10

Marginalized posterior probability distributions for the four parameters C2,T0, and α of the 1BB model fit, with the 1-dimensional projections shown along the diagonal, and the 2-dimensional projections in the other panels. Contours mark the 10%, 25%, 50%, 85% and 99% confidence intervals.

Open with DEXTER

Table 5

Global mean and errors for the parameters of our various models, and reduced χ2 of the fit.

In Fig. 10, we show the marginalized posterior probability distributions for the parameters C2,T0, and α. As can be seen, C2 and T0 are strongly anti-correlated. However, any increase or decrease in C2, hence decrease or increase in T0, will only shift the resulting temperature curve in the vertical direction, owing to our approach Eq. (4). The significance of the temperature variations justifying this approach has already been shown with our single-epoch fits (Fig. 7) and can also be seen from the SEDs shown in Fig. 6. The reverberation delay τ does not seem to show significant correlations with any other parameter, while α is slightly anti-correlated resp. correlated with C2 resp. T0. Interestingly, we observe a multimodal pdf of the reverberation lag τ. This is on the one hand caused by slightly different reverberation delays in 2010 and 2013–2014 (see Sect. 5.3), but mainly by a strong bimodality in τ in the 2012–2014 part of the data set (epochs 7–29, see Fig. 11). Because that second period dominates the fit due to a higher amount of data points and lower photometric errors, the bimodal pdf is also visible in the fit of the complete data set. This bimodality is discussed in more detail in Sect. 5.3.

4.2. Further structure function models and the time-variable power-law index

As described in Sect. 3.2, we tested for the robustness of all our performed multiepoch, multiwavelength fits using three different structure function models. It turned out that the results are highly stable under the exchange of the structure function model or its particular parameters. In Table 5, the results are only listed for the power-law model, while Table 6 shows the influence of the structure function on a subset of our results.

Table 6

Global mean and errors for the parameters in the first period (epochs 1–9) and in the second period (epochs 7–29) for the 1BB model.

Table 7

Global mean and errors for the parameters in the first period (epochs 1–9) and in the second period (epochs 7–29) for the 1BB model.

We also tested the influence of the AD power-law slope α on our results. As an alternative approach to keeping α fixed over the whole time and flux range, we allowed for a varying α. Since the continuum emission from the AD is found to get harder as the AD brightens (see e.g. Trèvese et al. 2001 and references therein), we performed an alternative 1BB fit allowing for a varying power-law slope of the form αvar(t) = α + α1·Lz(t)/ ⟨ Lz(t) ⟩. We obtained best fit values α = 1.56 ± 0.07 and α1 = 0.10 ± 0.07, while all other parameters stayed nearly the same (Table 5). At this point, we note that unfortunately the power-law slope is only well-determined in those epochs where optical measurements were also available (epochs 7–9). Thus, our applied approach of fitting a varying α to our data might be insufficient to determine the true variability range of the AD power-law slope. Therefore, we alternatively used published empirical relations between α and the AD luminosity. According to Trèvese et al. (2001), who analyzed multiepoch data of a sample of quasars, there is a relation Δα = a + bΔ(log Fν) between the change of the spectral index and the logarithmic optical continuum flux change of the AD, with a = (−8.49 ± 5.50) × 10-2 and b = 2.55 ± 0.75. Specifically for NGC 4151, Fanti et al. (1984) report the relation α = b·log Fν with b ≈ 4. Making use of the inferred α values in epochs 7–9 and the flux differences of the z band signal with respect to the z flux in these three epochs, we applied the cited two models for deriving an alternative, more representative evolution of α over the whole time and flux range. Thus, instead of fitting the evolution of α to our data, which may problematic because of the absence of optical data in most epochs, we now use those two models as input for our fit. We find our results to be qualitatively robust under the use of these different models (also see Sect. 5.3), even though values as high as α ≈ 3.3 are reached in the high-flux epochs. The influence of the particular choice of α for a subset of our results is shown in Table 7.

5. Discussion

5.1. Single-blackbody model

As for our 2010 data of NGC 4151, we observe a significant rise of the emission and temperature of the innermost hot dust, following states of increased AD brightness. The hot dust temperature follows the AD flux with a time delay of roughly 31 days. This indicates that the hot dust in NGC 4151 currently observed is simply heated up by increased AD irradiation and not destroyed owing to sublimation. In the case of significant dust destruction, one would expect an increase in the reverberation delay, which is not seen in our data. Obviously, the major part of the hot dust in the nucleus of NGC 4151 is not located at its current sublimation radius, but is cooler than sublimation temperature.

There are strong indications that the hot circumnuclear dust around AGNs mainly consists of large graphite grains (0.2 μm grain size), with sublimation temperatures 1500 K (Gaskell et al. 2004; Kishimoto et al. 2007, 2011). Since our inferred dust temperatures do not reach 1500 K (see Figs. 7 and 12), it is to be expected that we see no dust sublimation in our data.

Limitations of our method are given by our z-band sampling, and thus the interpolation of our input AD signal. Because, for example, the minimum observed in JHK around MJD 56 700 is missed in the z-band observations (therefore in the interpolated AD signal), the resulting model fits the JHK fluxes around MJD 56 700 very poorly. Apart from this mismatch, the 1BB model already fits the data remarkably well within the errors. Only in the epochs 4–6 are the modeled K-band fluxes systematically lower than the actually measured K-band peak fluxes.

5.2. Two-blackbody model

In the 2BB model, epochs 4–6 are fit better, indicating that a second BB component might actually contribute significantly to the observed K-band flux. However, we must point out that the goodness of the fit is not improved when using the 2BB model () instead of the 1BB model ()). This is not surprising, since the weight of the few 2010 HK data points is low compared to the global data set.

In our Paper I, we already argued that the hot, NIR dust around NGC 4151 might be substantially radially extended and better represented by more than one blackbody component.

However, when including data at longer wavelengths, i.e. the WISE W1-W3 bands, we must see that our 2BB model has to be rejected4. From Figs. 14 and 15, it becomes apparent that the fluxes are fit better with the 1BB model.

This finding agrees with the data of Burtscher et al. (2009). The unresolved point source that they measure with N-band interferometry and which they attribute to the inner rim of the hot dust torus (at 0.04–0.05 pc 50 light-days) is apparently the same source that constitutes our observed NIR fluxes, because the observed Burtscher PS fluxes are consistent with the SED of our 1BB model (Fig. 14). In Figs. 14 and 15, we went on to overplot the fluxes of an extended Gaussian source measured by Burtscher et al. (2009), and the sum of both components (PS + Gaussian). The extended source is interpreted as the warm component of the clumpy torus located farther out, at 2.0 ± 0.4 pc 240 light-days. Because this extended source is not resolved in the WISE photometry, it contributes to the WISE PSF flux of the W3 band (which matches the Burtscher total flux at 12.5 μm), while the W1 and W2 bands obviously show no significant warm dust contribution. A clear excess of the observed PS + Gaussian fluxes over the modeled 1BB fluxes can only be seen from the N band (λ ≥ 8 μm) on. Here, the single blackbody approximation is no longer sufficient, but a second blackbody component contributes to, or even dominates, the measured fluxes. This is, however, not the 698 ± 64 K blackbody component from our 2BB fit, but a blackbody of lower temperature () that is located farther out (Burtscher et al. 2009). A second blackbody component within a 100 light-day distance from the source (which we set as an upper limit in the prior for τ2, see caption of Table 4) is thereby clearly ruled out.

5.3. Variable time lag

thumbnail Fig. 11

Marginalized posterior probability distributions for the four parameters C2,T0, and α of the 1BB model fit, for the first part of the data set (epochs 1-9) in green and the second part (epochs 7–29) in red. The 1-dimensional projections are shown along the diagonal, and the 2-dimensional projections in the other panels.

Open with DEXTER

Besides a second blackbody component contributing to the NIR fluxes, there might be another explanation for the missed 2010 K-band peak by our fit. The resulting lag of τ1 = 31.0 ± 1.6 days found for the 1BB model fits our global data very well. This holds for all data except for the K-band peak in epochs 4–6. Here, it looks as if the actual lag might still be a bit higher. If we split the data set into two subsets (set 1: epoch 1–9 (2010 – beginning of 2012 data), set 2: epochs 7–29 (2012–2014 data)5), we find that the fitted hot dust reverberation delay is τ2010 = 42.5 ± 4.0 days if we include only epochs 1–9 (see Table 5 for the other parameter values), consistent with Hönig & Kishimoto (2011) and Kishimoto et al. (2011a). For the epoch 7–29 fit, we get a best-fit delay of τ2013 = 29.6 ± 1.7 days (also see Table 5). We excluded that this decrease in τ is merely an effect of the improved sampling in the second period. The blackbody constant also changes between the two different fits, from to , and the initial temperature from to . However, refers to the initial temperature in epoch 1 and thus has no real physical meaning for the second-period fit, since only epochs 7–29 are included. Figure 12 shows the resulting temperature evolution of this fit. Interestingly, while all other parameters stay nearly the same, the variability factor also changes between the two fits, from v2010 = 0.96±0.04 to v2013=0.81±0.03, which might indicate that the AD illumination is less efficiently reprocessed by the hot dust in the second period.

As discussed in Sect. 4, we tested the influence of different structure-function models and different evolutions of the AD power-law index α on our results. It seems particularly important to test the robustness of the apparent decrease in the reverberation lag under the different models. While the parameter values are very insensitive to the different structure function and α models for the global 1BB model, which includes all epochs, the parameters inferred for the epoch 1–9 fit and epoch 7–29 fit do undergo certain changes when applying the different structure functions and α models. The influence of the various structure function models on our results are shown in Table 6, while the influence of α is shown in Table 7. We do see a change in the derived absolute parameter values, especially the time lag seems to be sensitive to the applied variability model of α. Nevertheless, relative to each other, the inferred lags for the epoch 1–9 and epoch 7–29 fits and the decrease in the delay remain unchanged. While for a constant α (and different structure functions), the lag decreases from τ2010 ≈ 43 days to τ2013 ≈ 30 days, we infer τ2010 ≈ 37 days, τ2013 ≈ 28 days using the model according to Trèvese et al. (2001), and τ2010 ≈ 34 days, τ2013 ≈ 26 using the one presented by Fanti et al. (1984). Thus, our results are qualitatively unaltered: we see a significant decrease in the reverberation delay from 2010 to 2013–2014.

thumbnail Fig. 12

Hot dust temperature Tsingle derived from the single-epoch decomposition, where T is a free parameter in each epoch, versus the temperature Tmulti from the multiepoch fits, in which T evolves according to Eq. (4) and only the initial temperature T0 is a free parameter. For epochs 1–9, the temperature resulting from the epoch 1–9 fit is plotted, for epochs 10–29 the temperature from the epoch 7–29 fit. The minimum that is seen in Tsingle around MJD 56 700 (and also in the corresponding JHK flux minima) is not covered within the sampling of the z-band data (see Figs. 8 and 2), therefore we observe a substantial discrepancy between Tsingle and Tmulti around that date.

Open with DEXTER

In Fig. 7, the temperature evolution resulting from the epoch 1–9 and epoch 7–29 fits is plotted versus the temperature from the single-epoch decomposition (Eq. 5). Naively one would expect that due to the decreased AD radiation by 50% on average from 2010 to end of 2012–2014 (Lz2013/ ⟨ Lz2010 ≈ 0.49), the blackbody temperature would have decreased by a factor 0.51/4 ≈ 0.84 following LT4. This temperature decrease due to the dimmer AD would be roughly balanced by a temperature increase due to a new dust location farther inside by a factor of . Thus, one would expect the blackbody temperature for 2010 to be on the same average level as at the end of 2012–2014. Nevertheless, we observe a decreased average temperature in the second part of the data set (T2013/ ⟨ T2010 ≈ 0.86) in the single-epoch fits. This is confirmed by the temperature curve of our multiepoch fit for epochs 7–29 and achieved through a lower variability factor with v2013/v2010 ≈ 0.84. Obviously, in the second part of the data set, the dust is heated less efficiently by the AD radiation estimated from the z band flux. The physical reasons for this could be an increase in dust grain size or a geometrical cause. In a model proposed by Czerny & Hryniewicz (2011), the dust clouds resulting from an AD wind only become exposed to the AD radiation in the process of moving farther outside, because dust that is still located farther inside tends to rise only a short distance above the disk, whereas the dust farther outside has a greater height. The innermost dust would then be heated less efficiently assuming an anisotropic radiation characteristic of the AD. The discussed parameter changes between the two periods might point to a changed dust distribution.

thumbnail Fig. 13

Marginalized posterior probability distributions for the parameters T0 and τ of the epoch 7–29 1BB model fit, obtained by fitting the J, H, and K bands separately (also see text of this section). The J band fit is shown in blue, H band in green, and K band in red. The 1-dimensional projections of the pdfs are shown along the diagonal, and the 2-dimensional projections in the other panel.

Open with DEXTER

In Fig. 11, we show the marginalized posterior probability distributions of the parameters C2,T0, and α for the epoch 1–9 fit and the epoch 7–29 fit. Clearly, we see a significantly decreased reverberation lag from τ2010 = 42.5 ± 4.0 in 2010 to τ2013 = 29.6 ± 1.7 in 2012–2014. As already mentioned in Sect. 4, the pdf of τ2013 is bimodal. Interestingly, no bimodality is observed for τ2010. Here, the pdf seems perfectly unimodal. Possibly, the bimodality is washed out by the higher photometric errors in the first period. Increasing the errors in the second period to the level of the epoch 1–9 errors, however, only slightly washes out the bimodality but does not make it vanish completely. We furthermore excluded that the detected bimodality in epochs 1–9 in contrast to the missing bimodality in epochs 7–29 is merely due to the improved sampling in the second period. This supports our supposition that we see a changed dust distribution in the second period.

The resolved bimodality in τ2013 is caused by slightly different delays for the bands JHK (constituting the blackbody emission), as confirmed by performing separate fits for each band6, i.e. , and , pointing to a slight radial extent of the dust within a narrow region around τ2013. Interestingly, what we observe is not just a smoothly extended structure, but at least two separate blackbodies at two discrete, nearby, but different radii. This feature is visible in Fig. 13, especially in the H band. While the two-dimensional projections are hard to disentangle by the eye, in the right panel of Fig. 13 one can nicely see an “inflection point” of the two different reverberation delays around τ ≈ 30 days. The pdf peaks at two different delays, τinner ≈ 26 days and τouter ≈ 31 days. In the K band, τouter dominates the pdf, while τinner is slightly visible as well. In the J band, the pdf is clearly dominated by the shorter lag, while the longer lag can only be marginally “resolved” due to a slight asymmetry in the pdf. The intersection point between the pdfs of the two different delays always occurs at τ ≈ 30 days, for all three bands, and H and K even exhibit a minimum there. As expected, the fitted blackbody temperatures of the three bands decrease with wavelength, from K, K, and K. Due to missing information on the color in the single-band fits, the degeneracies between T0 and C2 are extremely high, so the fitted temperatures should not be taken too seriously. The detected 2BB structure within a narrow range around τ2013 does not contradict the result of our previous 2BB model fit, which was mainly performed to test for a significantly radially extended dust distribution rather than a compact distribution. The two distinct blackbodies resolved within the epoch 7–29 fit still represent a fairly compact dust distribution, since they are very close, and are therefore not in conflict with our finding from the 2BB model fit.

Our results seem to indicate a decreased time delay from 2010 (MJD 55 300–55 400) to 2012–2013 (MJD 56 000–56 300). This stands in strong contrast to interferometric observations by Kishimoto et al. (2013), who measured an increased delay, from roughly 40–50 days in 2010 to 70–80 days in 2012 – a delay that is clearly not consistent with our data.

thumbnail Fig. 14

Resulting fluxes from the 1BB model fit in the high temperature epoch 5, plotted over wavelength. Overplotted are our derived NIR and MIR fluxes (zYJHK and W1, W2, W3). The 1BB model roughly matches the MIR point source fluxes, given by Burtscher et al. (2009). Since the model was fit to our zYJHK data alone (also see Sect. 5.2), these are marked with green circles, while the other data are represented by triangles.

Open with DEXTER

thumbnail Fig. 15

Resulting fluxes from the 2BB model fit in the high temperature epoch 5, plotted over wavelength. Overplotted are our derived NIR and MIR fluxes (zYJHK and W1, W2, W3). While matching the observed HK fluxes better than the 1BB model, the 2BB model total fluxes by far exceed the WISE W1, W2 measurements. A second blackbody component inside of 100 days can be clearly rejected. Again, the model was only fit to our zYJHK data (marked with green circles, while the other data are represented by triangles).

Open with DEXTER

It is possible that this apparent shift in the time lag only represents the complicated and dynamic dust morphology with the dust being clearly extended and with clouds moving turbulently within this dust distribution. It seems possible that even slight changes in the dust distribution through motion of single clouds could cause noticeable changes in the measured torus response, e.g. by means of shielding. However, the shift in time lag could also mean that the dust radius has indeed decreased by inflow of dust from outside due to the low luminosity state of the AD and thereby decreased dust sublimation radius. Alternatively, BEL clouds could have been launched from farther inside due to decreased AD luminosity. Dust condensation only occurs when the clouds have expanded to roughly three times their initial radius, which happens at a distance d ≈ 1000·d0, with d0 the initial distance of the clouds from the source (Elvis et al. 2002). Dust condensation would set in after roughly three to nine years from the initial launching of the clouds. Interestingly, our measured decrease in time lag seems to continue the decline from τ ≈ 60 days observed end of 2008 (Pott et al. 2010) to τ ≈ 50–45 days in 2009/2010 (Kishimoto et al. 2009, 2011a). This total decrease in τ from 2008–2012 seems to track a decline in the AD flux from 2003 until 2006 with a delay of roughly five to six years (see Fig. 3 in Kishimoto et al. (2013), thus perfectly matching the timescale expected within the accretion wind scenario.

An upper limit for the infall velocity of dust produced by stars, which is moving in from outside, is the free fall velocity onto the central black hole (MBH ≈ 4.6 × 107MSun, Bentz et al. 2006). We obtain an infall velocity of vfree ≈ 0.01c at the location of the hot dust torus. It would thus take 1300 days to decrease the dust radius from 43 to 30 days. This time span seems only marginally possible within our estimate. Although these qualitative arguments seem to favor the accretion wind scenario over the inflow model, we note that detailed modeling of the different scenarios is needed to reliably conclude on the dust formation mechanism in AGNs.

Our measured 2010–2012 and 2012–2014 dust radii are inferred with one and the same method; however, the decrease over time of these measured radii will not fit the generally expected relation, because it cannot be reasonably related to a decay in the momentary AD brightness delayed by τ (as also visible in Fig. 3 in Kishimoto et al. (2013) for other observations of NGC 4151). A variable time lag (that possibly traces the averaged AD signal of several years before, but not the momentary AD signal) might help explain the intrinsic scatter found in the radius-luminosity relation of dust around AGNs.

6. Conclusions

We presented updated results from our 2010–2014 NGC 4151 photometric and spectroscopic data, which are part of our AGN hot dust reverberation project to monitor the evolution of the hot dust temperature and reverberation lag around AGNs. Our findings for NGC 4151 are:

  • Dust sublimation: if dust sublimation occurred in responseto increased AD flux, it would happen at the in-ner edge of the dust distribution, thereby increasing the timedelay. Although the AD brightness increased substantiallywithin the observed time range, we see no signatures of any dustsublimation traced by our data. In particular, after detectedAD flux increases by more than 30% in 2010, andby roughly 100% from March 2012 to November 2012, we ob-served no increase in the hot dust reverberation delay, thus rulingout significant dust destruction for the time range 2010–2014. Itseems that the hot dust in this galaxy is currently located beyond itssublimation radius.

  • Dust temperature: the hot dust emission, and moreover, the host dust temperature in NGC 4151 closely tracks the AD flux variations on a short-term response timescale, which is roughly one month on average. We measured maximum hot dust temperatures lower than 1500 K throughout 2010–2014. The large graphite dust grains that are typically assumed for the hot circumnuclear dust around AGNs have a sublimation temperature of Tsub ≳ 1500 K. It thus seems perfectly consistent that we saw no dust sublimation in the observed period.

  • Dust distribution radius: on a long-term timescale, we measured a change of the hot dust reverberation delay of 13 days in two years. Detailed comparison of our results with seemingly comparable interferometric experiments reveals different variations in the derived dust radii. These apparent observational inconsistencies could imply that the real dust distribution in NGC 4151 is fairly complex, so that the comparison of hot dust radii measured with interferometry and reverberation is not straightforward.

  • Dust distribution: In our 2010 data, we saw slight indications of a second blackbody component of 700 K, which could not be confirmed without broader wavelength coverage. Here, we presented a new analysis, now including WISE photometry, which rules out significant thermal radiation at 700K from radii smaller than 100 days.

To sum up, we found a decreased reverberation radius of the hot, circumnuclear dust in NGC 4151. Dust destruction in the observed epoch seems highly improbable from our data, but a slight change of the dust morphology seems likely. While the observed decrease in the reverberation delay matches the timescale expected within the accretion wind scenario perfectly, a radius decrease due to the inward motion of dust from outside seems only marginally possible, as estimated from an upper limit for the infall velocity of dust produced by stars (see Sect. 5.3). From our analysis, new dust formation in a cooling BLR wind appears to be more likely than a radius decrease due to inflow. We emphasize, however, that detailed modeling is indispensable to reliably distinguish the different dust formation scenarios in AGNs. In our AGN host dust reverberation project, we monitor roughly 25 additional bright and variable Seyfert 1 AGNs with the GROND camera (grizJHK bands, ESO La Silla) and in the optical with the All-Sky Automated Survey for Supernovae (Shappee et al. 2014). It will be interesting to see whether our results for NGC 4151 apply to other AGNs as well.


1

It is not yet understood well (Peterson & Horne 2006) whether the BLR clouds are infalling, outflowing or in rotation. Recent observations indicate that the signatures of low-ionization broad emission lines are consistent with gas that is infalling or in a rotating Keplerian disk (Grier et al. 2013; Peterson & Horne 2006), whereas the higher ionized lines seem to trace gas that is part of an outflowing AD wind (Crenshaw et al. 1999; Peterson 2006; Kollatschny & Zetzl 2013).

2

Based on observations collected at the Centro Astronómico Hispano Alemán (CAHA), operated jointly by the Max-Planck Institut für Astronomie and the Instituto de Astrofisica de Andalucia (CSIC).

3

The appropriate mean is calculated by when subtracted from the data.

4

The WISE data were not included in the χ2 of the fit, because the aim was to test for an extended structure in the NIR regime up to the K band. Nevertheless, any resulting 2BB model would have to match the WISE fluxes as well.

5

We split whole data set into these two overlapping data sets because, while it was obvious that the delay in epochs 1–6 would be longer than the delay in epochs 10–29, a priori the delay of epochs 7–9 was not apparent, because fitting only 3 epochs is rather ambiguous.

6

We fixed C1 and v to the best fit values of the epoch 7–29 fit for these single-band fits, to avoid high degeneracies between the parameters.

7

The DE-MC proposal scheme presented in Eq. (15) might at first glance seem to violate one of the basic assumptions for monitoring convergence with the -statistic of Gelman & Rubin (1992), Brooks & Gelman (1998), namely the assumtion that the m individual chains are independent of each other. However, Ter Braak (2006) demonstrates that the conditional stationary pdf of one individual chain does not depend on the states of the other chains and is identical for all chains, so that the joint stationary pdf p(x1,...,xm) of the m chains is given simply by the product p(x1) × ... × p(xm). Thus, the states of the individual chains x1,...,xm are independent of each other at any iteration step after the burn-in phase (Ter Braak 2006; Mengersen & Robert 2003), i.e., after the algorithm has become independent of the starting distribution. Convergence of the DE-MC algorithm can therefore be monitored using Gelman’s convergence criterion.

8

Only the second half is evaluated in order to lower effects of the influence of the starting distribution, i.e., the burn-in phase is discarded.

Acknowledgments

We thank the anonymous referee for valuable comments and suggestions that helped improve this manuscript. We are very thankful to R. Andrae and M. Fouesneau for helpful discussions and comments on this work. We thank R. van Boekel for providing a database of main sequence star atmospheres. K.S. acknowledges support by “IMPRS for Astronomy & Cosmic Physics at the University of Heidelberg”. B.M.P. and G.D.R. are grateful for the support of the US NSF through grant AST-1008882. B.S. is a Hubble, Carnegie-Princeton Fellow, and is supported by NASA through Hubble Fellowship grant HF-51348.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555.

References

Appendix A: Gelman convergence diagnostics

We check for convergence using the Gelman convergence diagnostics (Gelman & Rubin 1992; Brooks & Gelman 1998). Given m independent7 chains and 2n iterations of the sampler, we can calculate for the second half of the iterations8 in each step and for each parameter estimand x the potential scale reduction factor (PSRF) (A.1)which serves as an (over-) estimate of the true scale reduction factor . Here, (A.2)is an estimator of the true target variance σ2 and is calculated by a weighted mean of the between-chain variance B/n(A.3)i.e. the variance between the m chain means (where denotes the global mean over all chains and iterations), and the pooled within-chain variance W: (A.4)The true target mean μ is estimated by the global sample mean of x, i.e. , and the term B/ (mn) in Eq. (A.1) refers to the sampling variability of . As W always underestimates the true variance (Wσ2 for n → ∞) for any finite n, the PSRF will always overestimate the true scale reduction factor and can be used as convergence diagnostic:

High values of , i.e. values significantly above 1, indicate that further simulations may lead to an improved inference of the target distribution – either because the variance estimates in the numerator of Eq. (A.1) can be decreased further by running more iterations or because W will increase with continuing iterations, because the chains have not yet covered the complete target distribution. If is close to 1, it can be assumed that W has nearly converged to the target variance and that each of the m chains of n iterations is sufficiently close to the target distribution.

The above holds if the sampler is started with a sufficiently overdispersed starting distribution. Otherwise values close to 1 might falsely diagnose convergence (Brooks & Gelman 1998). Therefore, it is strongly recommended to graphically inspect the evolution of W and , to make certain that they are not still evolving.

All Tables

Table 1

NIR fluxes of the nucleus of NGC 4151, derived with ISIS.

Table 2

WISE W1, W2, W3 photometric fluxes of NGC 4151, derived with GALFIT.

Table 3

Optical continuum fluxes of the nucleus of NGC 4151 at 5100 Å.

Table 4

Upper and lower limits a,b for uniform pdfs of the model parameters, for the 1BB model (with constant and varying α) and the 2BB model.

Table 5

Global mean and errors for the parameters of our various models, and reduced χ2 of the fit.

Table 6

Global mean and errors for the parameters in the first period (epochs 1–9) and in the second period (epochs 7–29) for the 1BB model.

Table 7

Global mean and errors for the parameters in the first period (epochs 1–9) and in the second period (epochs 7–29) for the 1BB model.

All Figures

thumbnail Fig. 1

Nuclear fluxes of NGC 4151 derived with ISIS (green) and GALFIT (magenta), shown here for the K band. Within the photometric errors, we observe excellent agreement between the two methods. The agreement is on a similar level for all other bands, suggesting that both methods manage to separate the nuclear flux of interest from the bulge and disk.

Open with DEXTER
In the text
thumbnail Fig. 2

Calibrated 2010–2013 photometry of NGC 4151. The optical data are plotted in orange for the epochs that coincide with our NIR epochs, and in black otherwise. Photometry of the NIR data was done with ISIS, and is highly consistent with our previously determined fluxes with GALFIT (also see Fig. 1).

Open with DEXTER
In the text
thumbnail Fig. 3

Estimated (vij)1/2 (red dots) from the z band data versus τij, shown on logarithmic scale. The three different structure function models that we fit to these data (see text of Sect. 3.2) are represented by the green, magenta, and blue lines.

Open with DEXTER
In the text
thumbnail Fig. 4

z-band light curve, interpolated with the method of Rybicki & Press (1992) and Press et al. (1992). We interpolated using the three models described in Sect. 3.2. For clarity, only the interpolation resulting from the power-law structure function is shown here. This interpolated signal is the least-variance estimate of the stochastic process and therefore very smooth. The error of this prediction is indicated by the shaded gray region. Single realizations of the process have much more structure on short timescales and will also make excursions outside of the estimated error regions.

Open with DEXTER
In the text
thumbnail Fig. 5

Evolution of the parameter T0,1 (for one chain and the population mean) of the 1BB model, for one exemplary run. Also shown are the evolution of the pooled within chain variance W (Eq. (A.4)), maximum variance estimate (see Eq. (A.1)) and the resulting convergence parameter, i.e., the PSRF (Eq. (A.1)). Convergence was diagnosed at iteration step n = 13 599, where was reached for all parameters.

Open with DEXTER
In the text
thumbnail Fig. 6

Our data and the 1BB model fit for the epochs 3,5,9, and 10. Overplotted is the AD model, the 1BB model, and the sum of both components. From the data, we observe a clear rise in the blackbody flux, and a clear shift of the blackbody emission peak to shorter wavelengths until epoch 5. This emerging bump is confirmed by the blackbody temperature maximum in epoch 5 of our single-epoch fit (see Fig. 7).

Open with DEXTER
In the text
thumbnail Fig. 7

Hot dust temperature evolution derived from the single-epoch decomposition. The hot-dust temperature changes follow the relative AD z-band flux changes with a delay of roughly 30 days.

Open with DEXTER
In the text
thumbnail Fig. 8

Results of our 1BB model fit. Plotted are the K,H,J,Y,z, and 5100 Å data (green) over time. The bigger dots with errors bars in the 5100 Å panel mark the data points that coincide with our NIR epochs. The red line represents the blackbody contribution for each band, blue the AD contribution, and the black line is the sum of both. Epochs 1 and 2 are not included in the evaluation of the fit (as t<τ).

Open with DEXTER
In the text
thumbnail Fig. 9

Results of our 2BB model fit. Plotted are the K,H,J,Y,z, and 5100 Å data (green) over time. The bigger dots with error bars in the 5100 Å panel denote the data points that coincide with our NIR epochs. The red line represents the contribution of the inner blackbody component (at 29 light-day distance) for each band, the brown line is the second blackbody (at 67 light-days), blue the AD contribution, and the black line is the sum of all three components. Epochs 1, 2, and 3 are not included in the evaluation of the fit (as t<τ2). The resulting 2BB model fit only differs significantly from the 1BB model fit in the high-flux epochs 4-6, in the K band, as indicated by a red circle.

Open with DEXTER
In the text
thumbnail Fig. 10

Marginalized posterior probability distributions for the four parameters C2,T0, and α of the 1BB model fit, with the 1-dimensional projections shown along the diagonal, and the 2-dimensional projections in the other panels. Contours mark the 10%, 25%, 50%, 85% and 99% confidence intervals.

Open with DEXTER
In the text
thumbnail Fig. 11

Marginalized posterior probability distributions for the four parameters C2,T0, and α of the 1BB model fit, for the first part of the data set (epochs 1-9) in green and the second part (epochs 7–29) in red. The 1-dimensional projections are shown along the diagonal, and the 2-dimensional projections in the other panels.

Open with DEXTER
In the text
thumbnail Fig. 12

Hot dust temperature Tsingle derived from the single-epoch decomposition, where T is a free parameter in each epoch, versus the temperature Tmulti from the multiepoch fits, in which T evolves according to Eq. (4) and only the initial temperature T0 is a free parameter. For epochs 1–9, the temperature resulting from the epoch 1–9 fit is plotted, for epochs 10–29 the temperature from the epoch 7–29 fit. The minimum that is seen in Tsingle around MJD 56 700 (and also in the corresponding JHK flux minima) is not covered within the sampling of the z-band data (see Figs. 8 and 2), therefore we observe a substantial discrepancy between Tsingle and Tmulti around that date.

Open with DEXTER
In the text
thumbnail Fig. 13

Marginalized posterior probability distributions for the parameters T0 and τ of the epoch 7–29 1BB model fit, obtained by fitting the J, H, and K bands separately (also see text of this section). The J band fit is shown in blue, H band in green, and K band in red. The 1-dimensional projections of the pdfs are shown along the diagonal, and the 2-dimensional projections in the other panel.

Open with DEXTER
In the text
thumbnail Fig. 14

Resulting fluxes from the 1BB model fit in the high temperature epoch 5, plotted over wavelength. Overplotted are our derived NIR and MIR fluxes (zYJHK and W1, W2, W3). The 1BB model roughly matches the MIR point source fluxes, given by Burtscher et al. (2009). Since the model was fit to our zYJHK data alone (also see Sect. 5.2), these are marked with green circles, while the other data are represented by triangles.

Open with DEXTER
In the text
thumbnail Fig. 15

Resulting fluxes from the 2BB model fit in the high temperature epoch 5, plotted over wavelength. Overplotted are our derived NIR and MIR fluxes (zYJHK and W1, W2, W3). While matching the observed HK fluxes better than the 1BB model, the 2BB model total fluxes by far exceed the WISE W1, W2 measurements. A second blackbody component inside of 100 days can be clearly rejected. Again, the model was only fit to our zYJHK data (marked with green circles, while the other data are represented by triangles).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.