© ESO, 2011

1. Introduction

A key step in understanding the accretion process in active galactic nuclei (AGN) is to determine the role and structure of the putative obscuring torus. The obscuration here essentially refers to the radiative absorption by dust grains in the torus, which hides the UV/optical emission of the central engine and the broad line emission (see Antonucci 1993 for a review). In type 2 objects, which do not display these spectral components, our line of sight towards the nucleus is thought to be relatively close to the system’s equatorial direction along which the torus obscuration becomes effective. Type 1 objects are believed to be rather face-on and permit us to directly observe the central components unobscured. The torus dust grains thermally re-emit the absorbed energy, and the innermost dust sublimates at  ~1500 K, setting the inner boundary of the dust distribution. This dust sublimation region thus radiates mostly in the near-infrared (near-IR; wavelength λ  ~ a few μm), with outer dust grains emitting at lower temperatures in the mid-IR (~10 μm). Without much spatial information, the torus structure has long been inferred from the observed spectral energy distributions (SEDs) mainly in the near- and mid-IR, as well as theoretical arguments.

The SEDs were first modeled by the tori with smooth density distributions (Pier & Krolik 1993; Granato & Danese 1994; Efstathiou & Rowan-Robinson 1995). However, the case for the material being in a clumpy medium has theoretically been argued from an early stage (Krolik & Begelman 1988). Accordingly, there have been many efforts to develop clumpy torus models in the past decade (Nenkova et al. 2002; Dullemond & van Bemmel 2005; Hönig et al. 2006; Schartmann et al. 2008; Nenkova et al. 2008a; Hönig & Kishimoto 2010), and the torus structure has been inferred in the context of these clumpy torus models (Nenkova et al. 2008b; Ramos Almeida et al. 2009; Mor et al. 2009; Hönig et al. 2010; Ramos Almeida et al. 2011; Alonso-Herrero et al. 2011). However, we are finally beginning to have far more direct information with the long-baseline interferometry in both the near- and mid-IR by spatially resolving the inner structure.

1.1. Early type 2 investigations

Detailed interferometric studies in the mid-IR have targeted a few bright type 2 AGNs, with the first ones being NGC 1068 and Circinus (Jaffe et al. 2004; Tristram et al. 2007), using the mid-infrared interferometric instrument (MIDI) at the Very Large Telescope Interferometer (VLTI). These have achieved the first determination of the overall size and projected orientation of the mid-IR emission region with two-temperature elliptical Gaussian models. The AGN NGC 1068 was also successfully observed with the near-IR (2.2 μm) long-baseline interferometric instrument VINCI at VLTI (Wittkowski et al. 2004), and followed up in the mid-IR with more detailed MIDI observations (Raban et al. 2009). One hard aspect is that, since the nucleus is heavily obscured in type 2 objects, one needs to disentangle the effects of radiative transfer in order to study the intrinsic structure and distribution of the material. Another type 2 object, Cen A, was the first radio-loud object to be studied with the IR long-baseline interferometry (Meisenheimer et al. 2007), which was followed by more detailed observations by Burtscher et al. (2010). Interestingly, the corresponding mid-IR source is found to be dominated by an unresolved source, interpreted to be a synchrotron core, which might, however, also complicate the study of the structure.

1.2. Type 1 study in the near-IR

Studies of type 1 objects are largely unaffected by the inclination/obscuration. We can directly scrutinize the inner intrinsic structure in both the near-IR and mid-IR. In the near-IR, Swain et al. (2003) made the first measurement of the K-band (2.2 μm) visibility of the brightest type 1 AGN NGC 4151 with the Keck interferometer (KI), where they obtained quite high visibility (squared visibility V² ~ 0.8). This result was recently confirmed by Kishimoto et al. (2009a) and also Pott et al. (2010). Furthermore, Kishimoto et al. (2009a, 2011) found similarly high visibilities for in total eight type 1 AGNs. They argued that this is an indication of the marginal and partial resolution of the dust sublimation region in these type 1 AGNs. Firstly, they marginally detected a visibility decrease over the baseline length for NGC 4151, which would be robust evidence for a partially resolved structure. Secondly, its implied size calculated as a geometric thin-ring radius for NGC 4151 as well as those for the other type 1s are found to roughly match independent radius measurements from near-IR reverberations, which are thought to be probing the dust sublimation radius.

The near-IR reverberation radius R_τK is the light-crossing distance over the time lag between the variabilities in the UV/optical and near-IR (2.2 μm). The former emission is thought to originate from the central accretion disk (AD) and the latter from the hot dust grains in the sublimation region heated by the AD (e.g. Glass 2004; Minezaki et al. 2004; Suganuma et al. 2006; Koshida et al. 2009). The radius R_τK has been shown to be proportional to L^1/2 where L is the AGN luminosity (Suganuma et al. 2006; our specific definition of L is given in Sect. 3.1), thus, the interferometric ring radii at K-band described above also roughly scale with L^1/2. This is what we would expect for the dust sublimation radii (e.g. Barvainis 1987). In addtion, the measured ring radii have been shown to be almost unaffected by the near-IR AD component (Kishimoto et al. 2007, 2009a). Therefore, the case for the near-IR visibilities probing the dust sublimation region is strongly made.

Furthermore, Kishimoto et al. (2009a, 2011) pointed out that the ring radii derived from the observed visibilities are either roughly equal to or slightly larger than the L^1/2-fit to the reverberation radii as a function of L. Since reverberation measurements generally place weight on the small responding radii (e.g. Koratkar & Gaskell 1991), R_τK is expected to represent the radius very close to the inner boundary of the dust distribution. On the other hand, the interferometric ring radii rather correspond to brightness-weighted effective radii. In this sense, the L^1/2-fit reverberation radius gives a good normalization for the radial extent, taking out the simple L^1/2 scaling. This normalization is used extensively in this paper.

1.3. Type 1 study in the mid-IR

In the mid-IR, Beckert et al. (2008) presented the first long-baseline mid-IR interferometry for a type 1 AGN, NGC 3783. Kishimoto et al. (2009b) published two-baseline data for the same galaxy, revising the results shown in Beckert et al. (2008) with dedicated software (see more details below) and also adding one more data set taken at a different baseline. They proposed a generic way of implementing a uniform comparison of the interferometric data from different targets with different luminosities and at different distances, by using the inner radius R_in described above for normalizing the spatial scale probed by the inteferometer. Mainly using the data of NGC 3783, but supplementing them with those from other type 1s and the type 2 AGN NGC 1068, they argued that the radial brightness distribution in the mid-IR of type 1 AGNs seems to follow a relatively shallow power-law. This was a first attempt to directly map the radial distribution in type 1 AGN tori.

Subsequently, Burtscher et al. (2009) also presented two-baseline MIDI data for another type 1 AGN, NGC 4151, and measured the mid-IR emitting region size using a model of a Gaussian plus point-source, which gives a visibility curve relatively similar to that of a power-law at low spatial frequencies. Tristram et al. (2009) presented the results of a MIDI snap-shot survey for AGNs of both types 1 and 2, and discussed the mid-IR emission size as a function of luminosity, which has recently be further evaluated with additional data by Tristram & Schartmann (2011).

Here we present a first systematic investigation of type 1 AGNs with the mid-IR interferometer, expanding the study of Kishimoto et al. (2009b). Combining our mid-IR data with the near-IR interferometry and photometry, we aim to comprehensively map the radial structure of the thermally emitting dusty region in the mid- and near-IR. We aim to derive physical constraints as directly as possible from the observations, with the ultimate goal of finding whether the radial structure is systematically dependent on the properties of the central engine.

2. Observations and data reductions

2.1. MIDI

We observed six type 1 AGNs listed in Table 1 on 8−12 May 2009 (UT) using MIDI, which is a two-beam combiner working in the mid-IR wavelengths (Leinert et al. 2003). The observation log for each fringe track is listed in Table 2. We also collected a limited number of data sets for two of the six targets on 15 Mar. and 5 Aug. 2009 (Table 2). The adaptive optics system MACAO at each of the 8.2 m Unit Telescopes was locked on the nucleus of each target AGN. A single baseline was allocated for each night, and we used four different baselines with their projected length spanning from 27 m to 129 m. For one of the targets, NGC 3783, we also used in our analysis MIDI data taken in 2005 published by Beckert et al. (2008) and Kishimoto et al. (2009b).

The data were reduced using the software EWS (version 1.7.1; Jaffe 2004), but with modifications using our own IDL codes. The old NGC 3783 data were also re-reduced with the same procedure. The details of the process are described in the Appendix A. Briefly, we firstly implemented an additional background subtraction for the photometry frames using adjacent sky strips in the 2D spectrum. Secondly, to reduce errors in group delay determinations, we smoothed delay tracks over typically  ~20−40 frames. To avoid possible positive biases in the correlated flux, we also averaged  ~20−40 frames to determine the phase offsets. The system visibility was obtained from the observations of visibility calibrators taken right after or before each target observation with a similar airmass. These data were also reduced with similar smoothings above to calibrate out the effect of the time averaging. These calibrators are also the photometric standards found in the list by Cohen et al. (1999) and/or by R. van Boekel (priv. communication). We obtained the correlated flux, which was corrected for the system visibility, and total flux separately. The errors for the total flux and correlated flux were firstly estimated from the fluctuation of the measurements over time. The total flux spectrum for each object was then derived as the weighted mean of many measurements over different nights (for NGC 3783, the spectrum was derived separately for 2005 and 2009 data; see below). For the correlated flux measurements at adjacent uv points, we also took weighted means. Finally, for the correlated flux, we assigned a systematic uncertainty as a 5% of the total flux based on the analysis presented in Appendix A.

For NGC 3783, the mid-IR total flux apparently increased by  ~20% from 2005 to 2009 observations. To facilitate a direct comparison, we scaled the total flux and correlated flux from the 2005 observations by the same factor. This means that we assumed no visibility change in the mid-IR over this period. This is conceivable based on a comparison of near-IR interferometry data for NGC 4151 that found no significant visibility change over a one-year timescale (Pott et al. 2010; Kishimoto et al. 2011).

To test the accuracy of our MIDI reduction software specifically for targets at sub-Jy flux levels, we observed stars with flux  ~0.3−0.8 Jy at 12 μm during our May 2009 run, and also in Mar. 2009. Details are described in Appendix A. One significant issue is that, if the frames are not averaged over a sufficiently long time-interval, the correlated flux is affected by a positive bias apparently in the wavelength region for which there are lower (intrinsic) counts. Since MIDI count spectra often have a peak at shorter wavelengths, the resulting biased spectrum is often redder than it should be. In addition, coherence loss toward shorter wavelengths might not be fully compensated for by calibrator frames, leading also to redder correlated-flux spectra. With the data for the sub-Jy unresolved stars, we demonstrate that the dedicated codes can correctly estimate the correlated flux, and quantify the possible systematic uncertainty.

In Table 2, the quality of each fringe track is indicated by the maximum signal-to-noise ratio (S/N) per smoothed frame (see Sect. A.3). The analysis of the sub-Jy unresolved stars described above currently covers this S/N down to  ~2.5, thus in this paper we decided not to use data with a S/N lower than 2.5. Accordingly, five fringe tracks in Table 2 were excluded from our analysis.

2.2. The near-IR data

We combined the mid-IR MIDI data above with data in the near-IR to comprehensively map the radial structure of the thermally emitting dusty region. For the photometry in the near-IR, we have contemporaneous imaging data at five broad-band wavelengths from 0.9 to 2.2 μm for NGC 4151, IRAS13349, and H0557-385, taken with the Wide Field Camera (WFCAM) on the United Kingdom Infrared Telecope (UKIRT). The data for the former two objects were published in Kishimoto et al. (2009a), but, together with the new data for H0557-385, we present here the data for all three AGNs in Table 1. The wide field-of-view of WFCAM provides simultaneous point-spread-function (PSF) measurements from stars in the same field. We used them to implement accurate PSF + host galaxy fits, and measured the point-source-only SEDs in the near-IR (see Kishimoto et al. 2009a for more details). For the three other targets in the sample, we measured the point-source flux by applying the same procedure to 2MASS images.

We have the near-IR interferometry data taken with the Keck interferometer (KI) at 2.2 μm for two of the six MIDI targets, NGC 4151 and IRAS13349. The KI data were taken on 15 May 2009 (Kishimoto et al. 2009a), less than a week from our MIDI observations of these targets. For the general comparison over the sample between the mid- and near-IR, we also use the sample of in total eight type 1 AGNs, taken with the KI in the near-IR at 2.2 μm (Kishimoto et al. 2009a, 2011).

2.3. VISIR

To obtain accurate photometry and check the accuracy of the total flux spectra from MIDI, we took mid-IR images of each object in two passbands, namely at 8.6 μm and 11.9 μm with VISIR on UT3. These images were mostly taken within about four weeks of the MIDI observations, to minimize the effect of variability. The observation log and the results are shown in Table 3. For ESO323-G77, a VISIR spectrum was taken by Hönig et al. (2010) within several days of our MIDI observations.

All the VISIR images were taken in “perpendicular” mode with each dataset consisting of four image tiles. Using the data from the ESO pipeline (version 3.7.2), we implemented additional subtraction of the background determined at annuli of radii  ~1.5−3.0′′ from each source on each of the four tiles. After checking the stability of the flux calibration as a function of aperture radii, we determined the flux conversion factor at a radius  ~1′′ where its uncertainty was estimated from the scatter in the results from the four tiles.

2.4. Additional photometric data in the mid-IR and near-IR

When fitting the models described in the next section to the SEDs, we only used our MIDI total flux spectra and WFCAM (or 2MASS) point-source fluxes described above. However, we gathered in addition IRS spectra and IRAC imaging data from the Spitzer Heritage Archive, listed in Table 4, for all the targets. For NGC 3783, we also have VLT/ISAAC data taken on 2 Jan. 2011 (Hönig et al. in prep.; Table 4). These extra data points from Spitzer and ISAAC fall approximately on our model SEDs presented below, with possibly some additional flux from the host galaxy, thus are viewed as approximate upper limits, and not used in the spectral fitting. We summarize the usage of our data in Table 5.

2.5. Reddening corrections

We corrected the SED for Galactic reddening using values from Schlegel et al. (1998) listed in Table 1. We used the reddening curves of Cardelli et al. (1989), supplemented with that of Chiar & Tielens (2006) at long wavelengths. We also corrected the SED for the possible reddening within the host galaxy for IRAS13349, ESO323-G77, and H0557-385, using the A_V values listed in Table 1. For H0557-385, in terms of the optical continuum slope, Balmer decrement, and host galaxy inclination (Winkler 1992; Fairall et al. 1982), the reddening is likely to be between those of ESO323-G77 and IC 4329A, which is another Seyfert 1 galaxy known to be strongly affected by host reddening. For the latter two objects, the estimates of A_V from the different methods exist including that derived from the gradients of flux variation (Winkler et al. 1992), and are  ~1 and  ~2, respectively (Winkler et al. 1992). Without any additional solid information about the reddening value for H0557-385, we roughly infer that A_V = 1.5 ± 0.5.

3. Results

3.1. The whole data sets and prerequisites

Figures 1 and 2 present all the data for each of the six targets. The total and correlated flux spectra are shown in the top panel, the visibility as a function of spatial frequency in the middle, and the uv coverage in the bottom. In the top panel, all the additional flux data described in the previous section are also plotted. The MIDI data were binned by seven pixels (Δλ ~ 0.4   μm) with the error in each bin taken as the median of the errors over the binned spectral channels (since errors over adjacent channels are correlated).

For the visibilities, the spatial frequency is shown in units of cycles per R_in. As described in Sect. 1.2, as the most empirical measure of the inner radius R_in of the dust distribution, we adopt the L^1/2 fit to the K-band reverberation radius R_τK given by Suganuma et al. (2006). Here L is the UV luminosity of the central engine, which is defined as a scaled optical luminosity, L = 6νL_ν(5500   Å) (Kishimoto et al. 2007), assuming a generic SED over the UV/optical wavelengths. Here ν designates the frequency in the rest frame, and L_ν is the isotropic luminosity per frequency. The formula for the L^1/2 fit by Suganuma et al. (2006) is given in Kishimoto et al. (2007) as (1)If the dust temperature structure is primarily determined by the illumination from the central source, the radial temperature structure scales with the sublimation radius R_in. In this case, by normalizing a length scale with R_in, we can place various interferometric data on a uniform scale (Kishimoto et al. 2009b). A given configuration of an interferometer probes a certain spatial frequency, or a certain spatial scale represented by a spatial wavelength, which is the reciprocal of the spatial frequency. Here we normalize this spatial wavelength with R_in (see the upper axis of the middle panel of Figs. 1 and 2). This corresponds to writing the spatial frequency in units of cycles per R_in. With this normalization, we can eliminate the simple luminosity scaling originating from R_in being proportional to L^1/2, as well as a simple distance scaling. Different objects with different luminosities and at various distances have different angular sizes for R_in, but the data can be uniformly viewed when the probed spatial scale is normalized with R_in. We can further directly determine whether there is a luminosity dependence in the structure beyond the simple L^1/2 scaling.

Following the same idea, the total and correlated fluxes can be normalized using . They are shown in the top panel of Figs. 1 and 2 as and have the same dimension as the surface brightness in the rest frame of each target. This is equal to , where ν_obs, f_νobs, and θ_in are the frequency, observed flux, and angular size of R_in in the observed frame, respectively. The correlated flux in the top panel is color-coded with the spatial wavelength, while each visibility data point in the middle panel is color-coded with the (radiation) wavelength in the rest frame, as indicated by the color bars. The total flux, correlated flux, and visibility spectra are shown in Appendix B using more conventional units to facilitate direct comparisons with the data in the literature.

3.2. MIDI total flux spectra

With MIDI, total fluxes are measured separately and independently of correlated flux measurements. We derived a total flux spectrum for each object as a weighted mean of many measurements over different nights. The fluxes from our VISIR photometry (or VISIR spectrum in the case of ESO323-G77) were found to be either in excellent agreement (within a few percent) with or slightly larger (~10−20%; see Appendix B) than the mean MIDI total flux spectrum for each object. Since the latter is from the beams corrected with MACAO and from a smaller extraction window (adaptive mask in EWS), it is conceivable that it has less contribution from the extended emission, e.g. from the host galaxy. In the case of NGC 4151, which seems to have quite an extended mid-IR emission (see the low visibilities in Fig. 1a), some of the extended emission in the VISIR flux is probably excluded in the MIDI flux. All these comparisons with VISIR fluxes confirm the accuracy of the mean MIDI total flux, and its error might be slightly over-estimated in some cases.

In the middle panel for visibilities in Figs. 1 and 2 (and also in Fig. 3 discussed in the next section), we show the visibility uncertainties originating from each correlated flux measurement without the contribution of total flux errors, but show the total flux errors separately as data points of visibility unity at the lowest spatial frequency. In this way, we can clearly show the relative uncertainties between each correlated flux measurement and the total flux measurement. In addition, these total flux data points are shown at the spatial frequency with the baseline length equal to the single telescope diameter. As long as the source remains unresolved by a single telescope, the total flux gives the correlated flux at this spatial frequency, and the visibility is unity at this spatial frequency. Thus, the data points illustrate the constraint in such a case, which is expected to hold at least approximately for our sample of type 1 objects. We note however that these are only for the data presentation. For the model fitting in Sects. 3.5 and 4, we include the total flux error in the visibility error budget in a normal way.

3.3. Mapping radial structure

Each panel of Figs. 1 and 2 shares a common range of axes over different objects so that we can directly compare them. We can also place together the data for all the targets on a single plane. Fig. 3 is such a figure, showing the observed mid-IR visibilities of all the targets, each averaged from 10 to 12 μm in the rest frame of each object. In Appendix C, we also show the same figure with the spatial frequency axis on a linear scale.

Our mid-IR observations cover the range of several tens to hundreds of R_in in spatial wavelengths. Over this range, the observed visibilities are relatively concentrated in one locus, from 0.3 to 0.9. If the visibilities of different objects followed the same curve, this would mean that the size scales only with L^1/2 and objects share the same radial structure as inferred by Kishimoto et al. (2009b). However, this is not the case.

In Fig. 3, each object is color-coded with its UV luminosity as defined in Sect. 3.1. It is quite clear from the figure that within the small sample, which nevertheless spans over  ~2.5 orders of magnitudes in luminosity, the radial structure does depend on luminosity, i.e. the structure in units of R_in looks more compact in higher luminosity objects. This implies that the increase in the mid-IR physical size (e.g. in pc) with luminosity is more gradual than L^1/2. We discuss this luminosity dependence in Sect. 4.1.

While the results above are essentially model-independent, we attempt in the following to describe the radial surface-brightness distribution more specifically with simple models. The dotted curves in Fig. 3 are visibility functions for power-law brightness distributions, which we argue below approximately depict the observed brightness distributions (see more in Sect. 3.5).

3.4. Single Gaussian/ring versus power-law

We first consider the results for ESO323-G77, where we have the data with the least ambiguity. The total flux spectrum provides with an excellent match to the simultaneous VISIR spectrum (see Sect. 2.3; Fig. B.1) and we have the two datasets at almost exactly the same PA of the baselines (~38°) with good S/N. The visibility data are shown in Fig. 4 as a function of spatial frequency.

As in the middle panel of Figs. 1 and 2, the visibilities in Fig. 4 are presented in different colors for different wavelengths as indicated by the color bar. For a given baseline, shorter wavelengths probe the source at higher spatial resolutions, and if the overall source size is the same when seen at different wavelengths, then the visibility is expected to become lower, i.e. more resolved, toward shorter wavelengths. However, the observed visibility spectra are roughly flat, or increase toward shorter wavelengths, which means that the overall size of the source is smaller at shorter wavelengths. This is also generally seen in other objects, except at some objects’ shortest wavelengths where visibilities might be slightly underestimated owing to a coherence loss (see Appendix A.4).

On the other hand, when we consider the visibilities at each wavelength as a function of baselines (or spatial frequencies), we can easily see that a simple Gaussian or ring (which we treat here almost equally, since they give almost identical visibility curves at low spatial frequencies well before the first null; see Sect. 4.2) does not fit them well, as shown by the dotted and dashed Gaussian curves in the top panel of Fig. 4. If we derive a Gaussian HWHM or ring radius for each baseline measurement at a fixed wavelength, we obtain a smaller size for longer baselines as we show in the bottom panel of Fig. 4. We appear to observe here a brightness distribution where the corresponding Gaussian/ring sizes become progressively smaller at longer baselines. This probably means that there is a continuous distribution of brightness spread over a wide range of continuous spatial scales. This can be understood if the brightness distribution resembles a power-law distribution at a given wavelength, with its index generally becoming steeper at shorter wavelengths. This is what we advocated in Kishimoto et al. (2009b) (see their Fig. 1). This power law is at least partially expected physically, since the temperature of the directly illuminated dust grains is expected to follow a power-law with the radius from the illumination source, simply based on thermal equilibrium (e.g. Barvainis 1987).

3.5. Half-light radius for power-law brightness distribution

For given inner and outer radii, a power-law index specifies a power-law brightness distribution, and a convenient way of characterizing the corresponding size is to use the half-light radius R_1/2, the radius within which half of the total integrated light at a given wavelength is contained. This can be directly compared with the HWHM of a Gaussian brightness distribution, since R_1/2 is equal to the HWHM in a Gaussian case. When the brightness distribution is adequately described by a power law, the values of R_1/2 derived from visibilities at different baselines should all be the same for a given wavelength. The bottom panel of Fig. 4 presents the deduced R_1/2 for each visibility observed. We can see that the deduced R_1/2 is quite independent of the baseline lengths, in contrast to the Gaussian HWHM. This reassures that the power-law description in this case is quite adequate. We can also see that the size becomes smaller toward shorter wavelengths, as discussed above.

The disadvantage of adopting a power-law model and its corresponding R_1/2 is that we need to define both the inner and outer radii, which are set here to be 1R_in and 100R_in, respectively. However, we can also define another, very similar, more model-independent radius, the half-visibility radius R_V0.5 as follows. A visibility at a given low spatial frequency (within the first lobe of the visibility function before the first null) roughly corresponds to the fraction of the flux that is contained within the spatial scale being probed and remains unresolved by the interferometer’s configuration. Here the configuration is set by the observing wavelength λ and baseline length B. In this sense, R_1/2 roughly corresponds to the spatial scale probed by the interferometer for which the visibility becomes 0.5. We can refer to this spatial scale using the spatial wavelength Λ ≡ λ/B, which is the reciprocal of the spatial frequency for the configuration, and corresponds to the spatial resolution of the configuration (Kishimoto et al. 2009b; see Sect. 3.1). It is straightforward to show that the spatial wavelength Λ_0.5 at which visibility becomes 0.5 for a Gaussian is related to R_1/2 by (2)If we define the half-visibility radius R_V0.5 as (3)then this is equal to R_1/2 for a Gaussian, and can generally be referred to as a good approximation of the half-light radius model-independently. This radius can approximately be deduced by interpolating observations when observed visibilities cross over 0.5 at the longest available baseline. We can also say more generically that the visibility at a given spatial wavelength Λ within the first lobe gives the approximate fraction of flux originating from the region within the radius  ~Λ/4.5. We note again that for a Gaussian brightness distribution, all three quantities, namely HWHM, half-light radius R_1/2 and half-visibility radius R_V0.5, are identical.

In the top panel of Fig. 5, we compare the visibility functions of Gaussian and power-law distributions with various R_1/2. The curves are shown as a function of spatial frequency in units of cycles per R_1/2 (or spatial wavelength in units of R_1/2), so that we can directly compare different R_1/2 power-law cases with the same Gaussian curve. We can see that power-law curves with smaller values of R_1/2 become asymptotically more similar to the Gaussian/ring curve. The bottom panel shows the ratio of Gaussian HWHM, derived from each visibility, to R_1/2 as a function of the same spatial frequency. For extended power-law distributions, the dependence of the Gaussian HWHM on the baseline length is quite strong, but for steep distributions, the dependence on baseline becomes quite small. All curves cross the visibility V of 0.5 at spatial wavelength Λ close to 4.5 R_1/2, meaning that the half-visibility radius R_V0.5 closely represents the half-light radius R_1/2 in all cases.

The dotted curves in Fig. 3 are the visibility functions for the power-law brightness distributions with three different steepnesses, which are quantified by the normalized half-light radius R_1/2/R_in = 3, 10, and 30 (corresponding to power-law indices of −2.6, −2.0, and −1.5, respectively). These seem to closely describe the mid-IR surface brightness distribution over the sample. For the more compact structures seen in the higher luminosity objects, it is more difficult to differentiate the power-law from a Gaussian/ring for a given measurement uncertainty. However, based on relatively resolved cases such as in ESO323-G77, we infer that the power-law description with different values of R_1/2 over the sample is adequate at least approximately. We can then translate the observed visibility curves on the 2D plane (Fig. 3, or the middle panels of Figs. 1 and 2) into a value of R_1/2 at each wavelength, and then disscuss the relation between this size information and other observables.

We derived the best-fit normalized half-light radius R_1/2/R_in for each object by fitting the visibility curve of a power-law brightness distribution to the observed visibility data, separately for four wavelength bins, namely 8.2−10, 10−11, 11−12, and 12−13 μm in the observed frame. To ensure a uniform comparison at the same rest wavelength over the sample, we fitted a power law as a function of wavelength over the derived R_1/2/R_in (and its error) at the four wavelength bins, and used this fit to derive values at 8.5 and 13 μm in the rest frame. The results are summarized in Table 6. The uncertainties are formal 1-σ uncertainties for the fits, except for NGC 4151 and NGC 3783 where the targets are well-resolved (V < 0.5 at long baselines) and the uncertainty in R_1/2 was estimated as the dispersion in those derived from different baseline data. Since the dispersion in R_1/2 about its best-fit value was more symmetric in log space, the uncertainties are shown in dex.

We note that the fits were done with the inner and outer radius fixed at 1 and 100 R_in. The latter can affect the radius estimation for the extended cases of NGC 4151 and NGC 3783. However, as long as each source is at least approximately unresolved by the single telescope, the half-visibility radius R_V0.5 is quite well-constrained as can be seen in the middle panels of Figs. 1a and b as well as in Fig. 3 (see Sect. 3.2). Therefore, our R_1/2 estimation would also be quite robust for the two targets.

4. Discussion

4.1. Sizes at various wavelengths as a function of UV luminosity

Using the half-light radius R_1/2 derived in the previous section, we discuss below the overall relations of the size to observing wavelengths and flux-related quantities such as spectral shape, total flux, and luminosity.

Figure 6 shows R_1/2/R_in as a function of rest-frame wavelengths. For the near-IR, we have plotted the thin-ring radii from the KI data for the two objects in our present MIDI sample, NGC 4151 and IRAS13349 (Kishimoto et al. 2009a). The thin-ring approximation should be fairly robust in the near-IR, where the radial brightness distribution is expected to be much steeper and the power-law visibility curve becomes very close to that of a Gaussian/ring (Fig. 5). A more detailed estimation of R_1/2 would require knowledge of the shape of the innermost dust distribution, which is not yet well constrained. Thus, we simply adopt here the thin-ring approximation for the near-IR data. We would infer that the other four targets, for which no KI data are available, probably have near-IR ring radii of  ~1−2 R_in based on the KI-observed sample of type 1 AGNs (see Fig. 7; Kishimoto et al. 2011). From the mid-IR to the near-IR, the normalized emission size would then become certainly smaller with wavelength, very roughly following some power law for each object, but with different indices for different objects.

As we demonstrated in Fig. 3, the mid-IR emission size seems to have a luminosity dependence. This is now quantitatively shown in Fig. 7, where R_1/2/R_in is plotted against UV luminosity L. The data points are color-coded for the observing wavelengths, and we include the 2.2 μm thin-ring radii of all the KI-observed sample from Kishimoto et al. (2011). As we can see, the ratio of half-light radius to R_in becomes smaller with increasing UV luminosity. This implies that, for a given radial temperature distribution, the higher luminosity objects have a steeper radial density structure, at least in the mid-IR emitting radii, as we discuss in greater detail below. We also plot power-law fits to the R_1/2/R_in at 8.5 and 13 μm in Fig. 7.

We can also analyze the luminosity dependence of unnormalized R_1/2 in a physical scale. Figure 8 shows the half-light radius R_1/2 in pc plotted against the same UV luminosity. In the physical scale, the radius in the mid-IR increases with luminosity much slower than L^1/2. Power-law fits to R_1/2 in pc at 8.5 and 13 μm as a function of luminosities are (4)and (5)where we see that R_1/2 at 13 μm is consistent with being constant.

These results are contrary to those of previous studies of the mid-IR size as a function of luminosity by Tristram et al. (2009) and Tristram & Schartmann (2011), who concluded that the mid-IR size is consistent with being proportional to L^1/2. They used a single Gaussian for the size estimation, and this can lead to a large systematic error owing to the baseline dependency. When the radial brightness distribution has a relatively shallow power-law form, a Gaussian size would represent the emission size that is resolved by a given interferometer configuration, leading to the strong dependence on the sampled uv points. Thus, in this case, the simple Gaussian size is quite inadequate for representing the emission size for the whole distribution. Here we account for the effect of the uv sampling more properly, based on the inference that the brightness distribution is of roughly a power-law form. This enables us to perform a more accurate, uniform comparison over the sample, as shown in Fig. 8. To gain a greater accuracy and model-independency, we can always return to the uniform view of the visibility curves over the sample shown in Fig. 3.

The half-light radius in units of R_in has a quite close relationship with the mid-IR spectral shape. Figure 9 shows R_1/2/R_in plotted against spectral index a in f_ν ∝ ν^a in the mid-IR. The shape becomes bluer when the size is smaller. In terms of the radial structure, this is what we would expect, as the steeper radial structure looks bluer as the contribution from hotter dust becomes more dominant. We elaborate on this point later using simple models.

We can also calculate an average surface brightness in the rest frame using R_1/2 as , which is equivalent to where θ_1/2 is the angular size of R_1/2. Figure 10 shows this as a function of rest-frame wavelengths. This can be directly compared with the Planck function as shown in the same figure, where the comparison gives the brightness temperature¹ at a given wavelength. In the near-IR, all the objects show the brightness temperature  ~1400−1000 K, while in the mid-IR the brightness temperature ranges from  ~1000 K to  ~200 K. In the near-IR, the brightness temperature is quite close to the observed color temperature of  ~1400 K, meaning that the average surface filling factor at the innermost dusty region must be close to unity. Figure 11 shows the surface brightness as a function of UV luminosity, and the corresponding brightness temperatures at 2.2 and 13 μm are shown in the second and third y-axis, respectively. We can again clearly see the luminosity dependence.

The close relationship between R_1/2/R_in and the mid-IR spectral index (Fig. 9) and the relationship between R_1/2/R_in and UV luminosity (Fig. 7) implies that there should also be a correlation between the mid-IR spectral index and luminosity. In fact, the luminosity dependence of the radial dust distribution was already inferred by Hönig et al. (2010), based on the IR spectral comparison between Seyfert galaxies and high-luminosity type 2 QSOs. Here we directly show the distribution change by spatially resolving it.

4.2. Multiple-temperature Gaussian fit

While the power-law brightness description at each wavelength gives a relatively model-independent measurement of the emission size, it does not physically describe the change in the overall size with wavelength and the flux that we observe. A simple, though very approximate way of describing the power-law brightness with changing steepness over different wavelenths is to consider two or more rings/Gaussians of different size at different temperatures. At each wavelength, the multiple components with different sizes contribute differently to the visibility function, giving a wider range in spatial frequency than the range covered by a single ring/Gaussian. The contribution ratio of the multiple temperature components changes with wavelength, causing the overall apparent size to change with wavelength.

We adopted a Gaussian geometry for each component, though rings might provide a more physically accurate description of in particular the innermost dusty region. In terms of the visibility function, a ring and Gaussian essentially give the same curve at the low spatial frequencies before the first null (Fig. 5; see e.g. Millour 2008). Both a thin-ring and a Gaussian geometry have only one parameter, but a Gaussian permits us to more simply quantify the average surface brightness of the component, and its visibility curve (which is also a Gaussian) is the simplest (no nulls, no lobes).

We describe each Gaussian component as a brightness distribution S_ν at a radius r from the center given as (6)where B_ν is the Planck function at frequency ν and temperature T, and R_1/2 is the physical size of the Gaussian in HWHM, where HWHM = R_1/2 for Gaussian. The isotropic luminosity L_ν of this component is (7)The factor f₀ gives the average surface brightness with respect to the blackbody when all the emission is contained within HWHM (see more on emissivity and surface filling factor described in Sect. 4.4). In addition to these two Gaussians, we also include the accretion disk component, which is assumed to remain unresolved and have a spectrum of f_ν ∝ ν ^+ 1/3 at IR wavelengths longward of 0.8 μm (Kishimoto et al. 2009a, 2008).

Table 7 shows the results of two-Gaussian-component fits. We assumed that HWHM = 1 R_in for the  ~1400 K near-IR component here, except for the two objects for which we have KI data. We see that the  ~300 K mid-IR component has a HWHH of several to a few tens of R_in, which implies an effective temperature radial gradient of β from  ~−0.8 to  ~−0.4 where T ∝ r^β. This effective index includes the effect of the radial gradient of surface density distribution. In Sect. 4.4, we attempt to separate the temperature and density gradients using more physical but simple models.

4.3. PA dependence

The PA coverage of our baselines for each object is very limited. Since our targets are all type 1 AGNs, we do not expect to see a large PA dependence in visibilities. Nevertheless, we seem to see some PA dependence at least in one object, NGC 3783. This has already been pointed out by Hönig et al. (2010). The brightness distribution along the SE-NW direction seems more extended than that in the NE-SW direction.

We can complement our inner PA dependency study with optical polarization measurements. The projected direction of the system polar axis can be inferred to be either parallel or perpendicular to the PA of the optical continuum polarization, depending on the line polarization properties (Smith et al. 2004; Kishimoto et al. 2004). If a linear radio structure is detected, this direction also indicates that of the system polar axis. All the targets except IRAS09149 have optical polarization measurements (Brindle et al. 1990; Wills et al. 1992; Martel 1998; Schmid et al. 2003), and a clear radio axis is also available for NGC 4151 (Mundell et al. 2003). The resulting inferred PAs are summarized in Table 5, and shown in the bottom panel of the uv coverage for each object in Figs. 1 and 2 with polar and equatorial directions in dashed and solid lines, respectively.

To evaluate the elongation of the mid-IR emission in NGC 3783, we attempted to fit an elliptical Gaussian for the mid-IR component with the major axis fixed at the polar axis PA 136°, and the results are listed in Table 7. The emission seems to be a factor of  ~2 elongated along this polar direction, suggesting that some of the mid-IR emission originates from the polar region. More extensive data and analyses will be presented for this object by Hönig et al. (in prep.).

4.4. Power-law brightness distribution: temperature/density-gradient model

4.4.1. Model description

We now discuss the direct physical constraints on the radial structure of the inner dusty region thermally emitting in the IR. To this end, we use a simple but generic description of the radial surface brightness distribution with as few parameters as possible. To physically describe a face-on radial brightness distribution S_ν(ν,r) at radius r and observing frequency ν, we would need at least two separate distributions, namely the temperature and surface density distributions. Here we slightly generalize the surface brightness description in Kishimoto et al. (2009b) as (8)The idea behind this description is that, at each radius, the face-on surface brightness is dominated by the emission from the heated dust grains that have a maximum temperature at that radius (e.g. Hönig & Kishimoto 2010). These grains are probably at around the surface of the torus, and most likely directly illuminated by the central source, if there are such directly illuminated grains at that radius. Under this maximum-temperature approximation, we refer to the dust at T_max simply as heated dust grains/clouds. We note that our observations are insensitive to much colder, interior dust grains that may exist deeper along our line of sight at that radius². We write this maximum temperature T_max(r) at radius r as (9)We mainly consider a discrete, clumpy, or inhomogeneous cloud distribution. The surface brightness of each heated cloud at radius r would approximately scale with r in proportion to B_ν(T_max(r)). The additional factor f_ε(r) above, which we write in the form of a power law as (10)describes the radial run of the surface density (i.e. per unit area) distribution of these heated discrete clouds (i.e. at T_max), but also includes the emissivity of each heated cloud surface. More specifically, the dimensionless factor f_ε can be regarded as a surface filling factor multiplied by the emissivity. In the case of optically thick illuminated clouds where their emissivity does not depend sensitively on the radial distance from the illuminating source or the observing frequency (e.g. Hönig & Kishimoto 2010, see their Fig. 3), the factor f_ε(r) is roughly proportional to the radial surface density distribution of the heated clouds, or heated dust grains. Thus, as long as the emissivity is roughly constant over the radius and frequency, we can regard f_ε(r) as a normalized surface density of heated dust grains (i.e. again those at T_max). In this way, our description is generic enough to ensure that we can obtain as direct constraints as possible from the observations.

The model would be called a temperature gradient model, if the density distribution is set to be uniform, i.e. α = 0, and the brightness described only by the temperature index β. In this case, the index would be an effective temperature index, which would include the effect of the density gradient. Here we try to separate the temperature and density gradient. At least in the case of the direct illumination of a dust grain, we would expect that β =  −0.5 for large grains and β ~−0.36 for standard-size ISM grains (e.g. Barvainis 1987). By referring to these values, we attempt here to obtain constraints on the surface density distribution.

4.4.2. Mid-IR spectral slope and SED-size relation

In Fig. 9, we show a grid for the two indices α and β on the plane of the mid-IR spectral index and half-light radius. The upper-most grid with dotted lines is for T_in = 1400 K, which corresponds to the observed color temperature in the near-IR and is thought to be roughly equal to the dust sublimation temperature. The range of density index α from −2 to 0 for a given temperature index reproduces the observed range of normalized half-light radii R_1/2/R_in. However, the expected mid-IR spectral shape as shown by the dotted grid is systematically bluer than observed. In fact, the observed range of the mid-IR spectral index, as well as that of R_1/2/R_in, is reproduced well with the power-law distribution if T_in is  ~700 K, where the corresponding grid is shown in dashed lines in Fig. 9. We need to accommodate this size-color constraint.

The total flux spectra, including the mid-IR spectral shape, and the size information for all the objects can be summarized on the plane of the SED, normalized at 2.2 μm where all objects show similar brightness temperature (Fig. 10), plotted against R_1/2/R_in as shown in Fig. 12a. For objects without Keck near-IR interferometric data, we set their R_1/2 at 2.2 μm to have the expected range of values from 1 to 2 R_in (see Fig. 7; Kishimoto et al. 2011). In the mid-IR, the objects show a relatively uniform behavior. From the high to low luminosity objects, as the normalized size R_1/2/R_in becomes larger, the mid-IR spectral shape becomes gradually redder (which is already seen in Fig. 9), and the mid-IR total flux level gradually increases with respect to the near-IR flux, although two objects, IRAS13349 and H0557-385, seem to have slightly (a factor of a few) higher mid-IR flux level than the other four objects (see more below).

Overplotted in dashed lines in Fig. 12b on the same data set is the SED-size relation at three wavelengths expected from the power-law distribution with the temperature index β =  −0.36 over the range of density indices α =  −1.6 to −0.4 (the case with slightly steeper β is similar, but with a slightly different density index range). The total flux mismatch can be clearly seen across this plane. It is quite difficult to reproduce the wide range of R_1/2/R_in and the overall near- to mid-IR SED which is relatively flat in νf_ν, while keeping the relatively red mid-IR spectral shape, for this single power-law structure.

4.4.3. Mapping the radial structure with a bright near-IR rim

However, we would be able to roughly understand both the total flux and size information if, in addition to the power-law structure, there is a central brightness concentration on the scale of  ~R_in emitting nearly at the sublimation temperature of  ~1400 K. This concentration would correspond to a dust sublimation region, and could be a bright rim of the innermost dust distribution, such as the puffed-up rim that is often discussed for the inner region of young stellar objects (see Dullemond & Monnier 2010 for a review). The exact structure of this central concentration is not yet well constrained with our data. Figure 12c shows a grid of models with the power-law plus the central concentration represented simply by a 1400 K ring at r = R_in, which seems to closely reproduce the overall SED-size relation with a range of density indices. In this case, the two objects, IRAS13349 and H0557-385, simply have a slightly higher flux level in the power-law component with respect to the near-IR ring.

In Figs. 1 and 2, we also show the best-fit results with this model for each object and Table 8 summarizes the parameters. The normalization of the ring luminosity is parameterized in the same way as in the case of Gaussian given in Eq. (7), where R_1/2 is replaced by the radius of the ring. Figure 13 shows the resulting surface brightness profiles of the models for each object. The power-law component in higher luminosity objects becomes steeper, as expected. Figure 14 shows the corresponding normalized surface density distribution of the heated dust grains (i.e. with a maximum temperature at each radius). It becomes steeper, and also slightly denser, for higher luminosity objects. The density indices are from  ~0 for lower luminosity to  ~−1 for higher luminosity objects in our sample.

We note that the central temperature T_in of the mid-IR power-law distribution has to be as low as  ~700 K, to reproduce the relatively red mid-IR spectral shape. (We fixed T_in to be 700 K in the fits implemented above, and we would need visibility data that fill the gap between near-IR and mid-IR wavelengths to constrain T_in more tightly.) This would mean that there is a compact region of such a low temperature quite close to the radius R_in, which could be the shadowed region just behind the bright sublimating rim.

4.4.4. NGC 4151

In the case of NGC 4151, this low-temperature compact region is probably what we see in the mid-IR visibility curves, which become relatively flat at spatial wavelengths smaller than  ~100 R_in as shown in the mid-panel of Fig. 1a. In our interpretation, this is caused by both the high spatial frequency part of the low-T_in power-law component and the hot central brightness concentration. Here the effect of the latter on visibility is indicated in this panel (also in the other mid-panels in Figs. 1 and 2 for other objects), where the dotted curves show the visibilities without the bright near-IR rim.

These can also be understood in terms of the correlated flux observed at long baselines (see the top-panel of Fig. 1a), which provides an approximate spectrum of the unresolved part of the source. Its color temperature is  ~400 K (or  ~500 K if we exclude slightly less reliable data at wavelengths shortward of  ~10 μm), which is much higher than that of the total flux (~280 K) and much lower than the sublimation temperature.

The deviation of our model fit from the data might indicate that the low-T_in core has even more flux while the outer radiating source is even more extended. This could suggest a model of three components (including the near-IR component), which was implied by Burtscher et al. (2009), who modeled their two-baseline mid-IR data with a Gaussian and an unresolved source, and suggested another near-IR component. We note, however, that, puzzlingly, the correlated flux in their two-baseline data did not show the relatively blue color temperature, thus gave almost no indication of the size becoming smaller toward shorter wavelengths within the mid-IR.

4.4.5. Warm mid-IR and hot near-IR parts

The composite brightness structure discussed above, consisting of emission components of a warm mid-IR dominating part and a hot near-IR dominating part, also seems to be discernible in the SED, particularly in those of relatively red objects such as NGC 3783 and NGC 4151 (Fig. 1). The two-component explanation of the near to mid-IR SEDs has been advocated by different authors from SED analyses (e.g. Mor et al. 2009). Historically, the bump around 3 μm in the SEDs of QSOs, which could correspond to the hot near-IR part, has been known for decades (e.g. see Sect. V in Neugebauer et al. 1979).

A direct size constraint on the hot central brightness concentration can be derived from near-IR interferometry, and its structure might be a little distinct from the warm power-law distribution. For example, the near-IR visibility of NGC 4151 suggests that it has a quite compact structure with R_1/2 quite close to R_in, while most of its mid-IR emission is quite extended. In that case, another parameter, in addition to the heating luminosity, may be responsible for forming the structure of the inner region. We plan to explore these regions with further measurements in both the near-IR and mid-IR.

5. Conclusions

We have presented mid-IR long-baseline interferometry of six type 1 AGNs. Within the small sample, which nevertheless spans over  ~2.5 orders of magnitudes in the UV luminosity L of the central engine, we have shown that the radial structure of the dust distribution changes systematically with luminosity. We have argued that the mid-IR radial brightness distribution can approximately be described by a power-law with the inner boundary set by the dust sublimation radius R_in. In this case, the power-law steepness can be quantified with a half-light radius R_1/2. The mid-IR half-light radius R_1/2 in units of R_in becomes smaller (i.e. the power-law distribution becomes steeper) as the luminosity increases, with R_1/2 at 13 μm ranging from a few tens of R_in down to a few R_in over the range of L from  ~10^43.5 to  ~10⁴⁶ erg/s. This means that, in contrast to the results of previous studies, the physical mid-IR radii in pc is not proportional to L^1/2, but increases with L much more slowly. Our current estimate is that the radii are  ∝ L^{+0.21 ± 0.05} at 8.5 μm, and nearly constant at 13 μm.

This mid-IR size normalized by R_in is correlated with the mid-IR spectral shape, where the normalized size is smaller for bluer spectral shapes. The mid- to near-IR SED as well as the measured size in mid- and near-IR wavelengths can roughly be explained in a picture where there is a central brightness concentration at the scale of  ~R_in emitting nearly at the dust sublimation temperature of  ~1400 K, over an underlying power-law-like structure extending outwards. The latter structure is inferred to have an innermost temperature as low as  ~700 K, and a radial surface density distribution of the heated dust that becomes steeper from lower to higher luminosity objects, ranging from  ~r⁰ to  ~r^-1. The radial temperature run is consistent with T ∝ r^β of β ~−0.35 ± 0.05. The exact structure of the innermost concentration will further be constrained with the planned near-IR interferometry.

---

Acknowledgments

This research is primarily based on observations made with the European Southern Observatory telescopes obtained from the ESO/ST-ECF Science Archive Facility under the program 083.B-0452. Additional data sets were collected under the programs 082.B-0330 and 083.B-0288. The authors are grateful to Roy van Boekel for providing the spectral database of MIDI calibrators. A part of the data presented here was obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology (Caltech), the University of California and the National Aeronautics and Space Administration (NASA). The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The Keck Interferometer is funded by NASA as part of its Exoplanet Exploration program. We are also grateful for the data obtained under the service observing program of the United Kingdom Infrared Telescope, which is operated by the Joint Astronomy Centre on behalf of the Science and Technology Facilities Council of the UK. This work is also based in part on observations made with the Spitzer Space Telescope, obtained from the NASA/ IPAC Infrared Science Archive, both of which are operated by the Jet Propulsion Laboratory, Caltech under a contract with NASA. The research is also based in part on observations with AKARI, a JAXA project with the participation of ESA. S.H. acknowledges support by Deutsche Forschungsgemeinschaft (DFG) in the framework of a research fellowship (“Auslandsstipendium”). This work has made use of services produced by the NASA Exoplanet Science Institute at Caltech, and the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, Caltech, under contract with NASA. This research has also made use of the SIMBAD database, operated at CDS, Strasbourg, France.

References

Online material

Appendix A: Reduction and calibration of MIDI data for faint targets

We used a part of the software EWS (Jaffe 2004; version 1.7.1) to reduce the MIDI data described in this paper. The mid-IR correlated flux of all the targets here is  ≲ 0.5 Jy, which is much fainter than normally handled with the software. Therefore, we implemented several additional procedures to reduce and calibrate the data as described below, using our own IDL codes. In short, we simply averaged or smoothed over a relatively large number of frames, designated as w here (typically w ~ 20−40) to determine both the group-delay and phase-offset tracks, and applied the same averaging to the calibrator frames to compensate for the side effects of the averaging.

We also observed several calibration stars, listed in Table A.1, which are fainter than  ~0.8 Jy at 12 μm, have known spectral type, and are unresolved with VLTI/MIDI. The deduced correlated flux is directly compared with the total flux, which is obtained by scaling the template spectra of Cohen et al. 1999 (not from less accurate total flux frames obtained with MIDI). This facilitates a reliable evaluation of the accuracy of our reduction and calibration procedures. The scaling factors are derived from both IRAS and AKARI measurements. On the basis of the close match between the two, we estimate that the total flux accuracy is  ~4%.

A.1. Large gsmooth parameter

When the group delay is determined from the delay function (the Fourier transform of each fringe spectrum from the frequency domain to the delay domain), several frames are averaged or frames are smoothed over time direction as specified by gsmooth parameter in EWS to suppress instrumental delay peaks. The default is four, but we used a larger number (typically 10−20) to increase S/N in delay peak determinations, as normally done for faint targets.

A.2. Significant averaging over time for phase-offset determination

Figure A.1(a) compares the calibrated correlated flux of the unresolved star HD107485 with the known total flux spectrum (i.e. the scaled spectrum described above), and Fig. A.1(b) shows the deduced visibility (using the known total flux). The raw count spectrum of the correlated flux is shown in Fig. A.1(c). The star has a peak count rate of  ~800 cts/s in the PRISM-mode spectrum in this particular observation. In these figures, the dotted line shows the spectrum from a default reduction (though with the same large gsmooth parameter), which determines phase offsets with a small averaging width. In this case, the spectrum suffers a red bias, i.e. extra counts at the longer wavelength side.

This seems to be at least partly due to the determination of phase offsets at too low S/N. This can be visualized on the complex plane of the broad-band correlated flux counts integrated over the N-band of  ~7 to 13 μm. The phase offsets are determined from this integrated broad-band flux. When a large enough number of frames are averaged, the distribution of the broad-band correlated flux measurements on the complex plane shows a (partial) circle as shown in Fig. A.1(d), where its radius represents the amplitude of the correlated flux with some phase offset slowly changing over time.

A noisy phase-offset determination (and subsequent rotation using the erroneous phase) would lead to a positive bias in the correlated flux (e.g. it would detect non-zero correlated flux even if the real correlated flux is zero). The red bias is probably caused, since the MIDI count spectra usually have a peak at shorter wavelengths, by the phase-offset determination having a larger weight at shorter wavelengths (or more accurately count peak), thus leading to the longer wavelength side having a stronger positive bias. We note also that, in general, correlated flux spectra would also tend to be redder if larger coherence loss at shorter wavelengths is not correctly calibrated.

Therefore, we chose a large w for the phase offset track determination, and a boxcar smoothing was applied to the group delay track for the same number of frames. By applying the same averaging to the delay and phase tracks of calibrator frames, we estimated the low effective system visibility arising from this averaging process, and calibrated it out from the target frames. As we show below, we chose w to maximize the S/N of the phase-offset determination and the result is insensitive to the exact value of the chosen w.

A.3. Optimum averaging width and evaluation of atmospheric conditions

Here we quantify the effect of this averaging process as a function of w, and describe how we can choose an optimal value of w.

A set of broad-band correlated flux measurements for a given w as shown in the complex plane in Fig. A.1(d) gives a distribution of the correlated flux amplitude measurements. The corresponding histogram is shown in Fig. A.1(e). We can evaluate the distribution as a function of w by calculating a fractional dispersion σ/m (standard deviation σ divided by mean m) for each w. This statistical property would be valid as long as the effective number of correlated flux measurements is not low. Here, the effective number of correlated flux measurements is the number of non-rejected frames (see Sect. A.5 below) divided by w. The reciprocal m/σ corresponds to the S/N of the correlated flux measurements per smoothed frame.

Figure A.1(g) shows the change of σ/m as a function of w (or corresponding averaging time interval, “accumulation time”; see upper axis) for the target and its visibility calibrator. The curve is cut when the effective number of measurements becomes lower than 10. In the small w regime for the faint target, the dispersion naturally decreases with larger w as photon noise decreases. The photon noise here comprises a background Poisson noise and fluctuation due to background subtraction residuals. This can be estimated from the noise spectra (a few hundred frames with large OPD offsets, taken at the start of each fringe track with MIDI) and the fractional fluctuation of the delay function peaks. This is shown in gray squares in Fig. A.1(g), and closely matches the σ/m curve at small w range. The fractional dispersion σ/m starts to increase at some width where the effect of a large time-averaging width to smear and effectively reduce the correlated flux becomes more dominant over the photon noise decrease.

For the bright calibrator, σ/m stays relatively constant at small w where flux fluctuation slowly decreases (indicated by gray plus signs), while the smearing effect of averaging slowly increases. The dash-dotted curve is the relative correlated flux counts of the bright calibrator (normalized at the smallest w), showing a relative effective system visibility at each w. For larger w, σ/m increases quickly (system visibility declines quickly) when the blurring effect of the large-w averaging starts to dominate. This gives an indication of the coherence time in the mid-IR at the time of the observation (~1.5 s in the case of Fig. A.1(g)). In the large w regime, the fractional dispersion σ/m is thus determined predominantly by the atmosphere at the time of the observation, rather than the observed source.

The coherence time and thus the blurring effect is approximately the same between target and calibrator frames when the two observations take place close in time and on the sky. In this case, we can actually calculate the expected σ/m for the target frames as a function of w from the calibrator frames. We subtract the photon noise (gray plus signs) in quadrature from the σ/m of the calibrator (solid plus signs with dashed lines) to obtain the dispersion caused by averaging, and add it in quadrature to the photon noise in target frames. This expected σ/m curve for the target is shown as triangle symbols with a dotted line in Fig. A.1(g), and closely matches the observed σ/m, meaning that atmospheric conditions were very similar. In this case, we simply apply the same averaging width for the calibrator frame to obtain an appropriate estimate of the effective system visibility spectra, or transfer function over wavelengths, for the target frames processed with this particular averaging width.

In Fig. A.1(g), the green circles show the relative calibrated broad-band correlated flux (normalized at w giving minimum σ/m, or maximum S/N) as a function of w, and the green crosses show the relative power-law index of the calibrated correlated flux spectrum. The calibrated correlated flux is systematically higher, and often redder, at small w because of the positive bias from low S/N determination of phase offsets. We have chosen w to yield maximum S/N for the target. If w is chosen to be larger than this value, the result is quite insensitive to the exact value of w, although the effective S/N of correlated flux measurements is lower at larger w owing to the smearing effect of averaging.

A.4. Adjustment for the atmospheric difference

Figure A.2 shows a similar case to that shown in Fig. A.1, where the count rate is slightly lower and the coherence time is slightly shorter. The maximum S/N is thus slightly lower. However, the atmospheric conditions between target and calibrator observations seem to be well matched.

In Figs. A.3 and A.4, we show the cases where there seems to be a difference in atmospheric conditions between target and calibrator fringe tracks. In these cases, we try to evaluate the difference by comparing the expected σ/m curve for the target (triangles with dotted line in Figs. A.3(i) and A.4(i)), as derived from the calibrator frames, with the observed curve for the target (squares with solid lines). We find that we can fit these two curves if we simply scale the time length of the calibrator frames (a factor of  ~0.7 in Fig. A.3(g), and  ~2.5 in Fig. A.4(g)). The scaling factor probably corresponds to the difference in the coherence time between target and calibrator frames. Therefore, we infer that we can obtain appropriate calibrations of the target frames if we use this time scaling factor. This correction also stabilizes the calibrated correlated flux over different w.

Figure A.5 shows the calibrated visibility of the 12 low and intermediate S/N observations of unresolved stars as a function of the maximum S/N per smoothed frame. We calculated the visibility in three different wavelength bins, avoiding the range affected by atmospheric ozone features. We confirm that the calibration procedure above correctly calibrates the data, although we see a slight tendency that, toward shorter wavelengths, the correlated flux is underestimated as the maximum S/N decreases. Based on the mean deviation and dispersion of these calibrated visibilities, we assign a systematic error in the correlated flux as 5% of the total flux, when the maximum S/N is at least larger than  ~2.5. Since the cases at lower S/N remain untested, we did not use the target datasets under this level in this paper.

A.5. Refinement of group delay determination

As described in Sect. A.1, we use a large gsmooth parameter for faint targets to smooth over a large number of delay function frames and thus increase the S/N for delay determination and also strongly suppress the instrumental delay peaks. However, this also suppresses the atmospheric delay peaks at the time intervals where the delay is relatively quickly changing. As a result, the delay function image becomes quite “blobby”, and strong delay peaks are left only where the atmospheric delay did not change rapidly.

To partially compensate for this effect, we implemented a simple iteration to determine the delay track. Namely, we first interpolated (over time) between the strong delay peaks to derive an approximate delay track. Then we used this approximate atmospheric delay track (plus instrumental delay) to rotate the original fringe spectra in the complex plane, and re-derived the delay function smoothed with the same gsmooth parameter. This makes the resulting delay peaks line up straighter over time, and thus recovers a larger number of frames with strong delay peaks. We cycled over this process a few times. For calibrator frames, this resulted in almost all the frames having strong peaks (i.e. even with the large gsmooth parameter). For target frames, significant fraction of frames is recovered, resulting in a good representation of the delay track even in a relatively low S/N case as shown in Fig. A.3.

In low S/N cases, we flagged frames that show delay peaks at positions that deviate significantly from the overall track obtained above. This was done by taking the histogram of the difference between the determined group delay track with that smoothed over many frames (typically  ~20−40 frames, set as 2  ×  gmoooth). The distribution can roughly be descibed by a Gaussian, and frames with large deviations (>10σ) were excluded. In this rejection process, we did not select frames based on the strength of the delay peaks in order to avoid possible biasing.

A.6. A note on total flux frames

For the targets observed with UTs that have total flux of less than a few Jy, it is better to implement an additional background subtraction using the sky regions very close to the target position. The residual of the primary background subtraction using chopped frames can still dominate the target flux for these relatively faint targets. An optional way to further reduce the systematic error from the background subtraction is to iteratively reject frames with the sky region counts having large deviations from the average. This was implemented in a small number of cases when the time fluctuation over the photmetry frames was significantly reduced after the frame rejection. With this additional sky subtraction, we averaged many sets of data, even from different nights, to obtain reliable total flux spectra, as confirmed with our VISIR photometry and spectrum.

Appendix B: Mid-IR data plotted in conventional formats

In Figs. 1 and 2 we show our flux and visibility data with normalizations using the inner radius R_in. We show the same data in the mid-IR (after reddening corrections; see Sect. 2.5) using conventional units in Figs. B.1 and B.2, to facilitate direct comparisons with the data in the literature.

Appendix C: Radial structure plot in linear scale

In Fig. 3 we showed observed mid-IR visibilities of all the targets in a single plot to uniformly compare the radial structure of the objects in normalized units with spatial frequency in log scale. Here in Fig. C.1 we show exactly the same figure but with spatial frequency in linear scale, again to facilitate comparisons with the data in the literature.

All Tables

All Figures