© J. W. Isbell et al. 2022

1. Introduction

Active galactic nuclei (AGN) are thought to play a crucial role in the formation and evolution of its host galaxy. Moreover, understanding the dust in the vicinity of supermassive black holes is key to understanding how AGN are fed and how they interact with their hosts. The dust traces dense molecular gas which feeds the AGN. Large, obscuring dusty structures are thought to be responsible for both funneling material toward the central engine, and for distinguishing between Seyfert 1 and Seyfert 2 AGN. In the original Unified Model of AGN (Antonucci 1993; Urry & Padovani 1995; Netzer 2015), a central obscuring torus of dust is oriented such that the broad-line region of the AGN is directly visible (Seyfert 1) or such that its observation is blocked by the torus (Seyfert 2; hereafter Sy2). So in order to fully understand the accretion process and the life cycle of an AGN, one must understand the parsec-scale dust structures surrounding it.

The so-called torus is comprised of several key features which vary in temperature from < 200 K to 1500 K and scale from tenths of a parsec to hundreds of parsecs. The inner edge is the radius at which radiation from the accretion disk (AD) causes the dust to sublimate. The sublimation radius is dependent on both the luminosity of the AD and the dust composition, but typically ∼0.1 pc for dust evaporating at 1500 K, for a L ∼ 1 × 10¹⁰ L_⊙ AGN. Beyond the sublimation zone, it is thought that a dense disk or torus of material is responsible for both hiding the broad line region (BLR) in Sy2 AGN and for feeding the AD. Previous mid-infrared (MIR) interferometric studies revealed that many “tori” have an additional component in the form of a polar extension (see, e.g., Hönig et al. 2012; Burtscher et al. 2013; López-Gonzaga et al. 2016; Leftley et al. 2018), the Circinus Galaxy’s chief among them (Tristram et al. 2007, 2014). The polar component is thought to be a radiation-driven outflow (e.g., Wada 2012; Wada et al. 2016), and it can represent a key mechanism of AGN feedback. This is called the fountain model, and it was shown by Schartmann et al. (2014) to reproduce the MIR polar extension and dusty hollow cone in the Circinus Galaxy (hereafter Circinus). A key finding of spectral energy distribution (SED) fits to nearby AGN as well as comparisons to radiative transfer models is that the dust in the central structures (and particularly in the wind) must be clumpy, allowing dust to reach high temperatures and exhibit “blue” spectra even at large distances from the AD (Krolik & Begelman 1988; Nenkova et al. 2008; Hönig & Kishimoto 2017; Martínez-Paredes et al. 2020; Isbell et al. 2021). The exact nature of these components and how they are connected to each other and to the host galaxy remains an open question. A holistic model of the central dust distribution is shown in Izumi et al. (2018), but only the resolution offered by infrared interferometry can probe the subparsec details of the dust near the active nucleus.

The Multi AperTure mid-Infrared Spectro-Scopic Experiment (MATISSE) is the second-generation MIR interferometer on the Very Large Telescope Interferometer (VLTI) at the European Southern Observatory (ESO) Paranal site (Lopez et al. 2014, 2022). MATISSE combines the light from four unit telescopes (UTs) or four auxiliary telescopes (ATs) measuring six baselines in the L-, M-, and N-bands simultaneously. MATISSE furthermore introduces closure phases to MIR inteferometry. The combination of the phase measurements on any three baselines ϕ_ijk ≡ ϕ_ij + ϕ_jk − ϕ_ik is called the closure phase; this summation cancels out any atmospheric or baseline-dependent phase errors (Jennison 1958; Monnier 2003). Closure phases are crucial for imaging because they probe the spatial distribution of target flux and because they are unaffected by atmospheric turbulence. Recent imaging studies of NGC 1068 with GRAVITY (GRAVITY Collaboration 2020) and MATISSE (Gámez Rosas et al. 2022) have illustrated the power of this approach in revealing new morphological details and spatially resolved temperature measurements of the circumnuclear dust. Until this work, NGC 1068 was the only AGN to have been imaged with MATISSE.

Circinus is of particular interest as it is one of the closest Sy2 galaxies (at a distance of 4.2 Mpc Freeman et al. 1977; Tully et al. 2009) and the second brightest in the MIR (only fainter than NGC 1068). Circinus is a prototypical Sy2 galaxy, exhibiting narrow emission lines (Oliva et al. 1994; Moorwood et al. 1996) and an obscured broad-line region (BLR; Oliva et al. 1998), as well as bipolar radio lobes (Elmouttie et al. 1998) and an optical ionization cone (Marconi et al. 1994; Maiolino et al. 2000; Wilson et al. 2000; Mingozzi et al. 2019). Additionally, Circinus exhibits a Compton thick nucleus and a reflection component in X-rays (Matt et al. 1996; Smith & Wilson 2001; Soldi et al. 2005; Yang et al. 2009). Finally, in- and outflows and spiral arms have been observed in CO down to ∼5 pc scales (Curran et al. 1998; Izumi et al. 2018; Tristram et al. 2022), further indicating the complexity of the central structures.

Circinus was observed extensively with the first generation MIR interferometer, the MID-infrared Interferometric instrument (MIDI; Leinert et al. 2003), in the N-band (e.g., Tristram et al. 2007, 2014, hereafter T14). These observations showed a warm (∼300 K) dust disk roughly aligned with the water maser emission (Greenhill et al. 2003), but the flux was dominated by large scale (≳​100 mas) emission roughly orthogonal to the disk. The orientation of the large scale emission’s major axis was found to differ significantly from the optical ionization cone central angle (PA_opt. = −45° vs. PA_dust = −73°), and follow-up modeling work by Stalevski et al. (2017, 2019) has indicated that the polar-extended dust emission may come from an edge-brightened outflow cone.

The proximity and declination of Circinus (at around −60°) make it an ideal target for imaging with MATISSE, as it provides high spatial resolution (10 mas = 0.2 pc) and because its nearly circular uv-tracks aid in the production of high fidelity reconstructions. MATISSE provides the first MIR measurements of the closure phase, which sample the (a) symmetry of a source and are crucial for image reconstruction. Previous analysis relied on Gaussian model fitting, which is a smooth, simplified representation of the source emission; but interferometric image reconstruction has the potential to build on these results through model-independent sampling of the source structure. Herein we present the first image reconstructions of the N-band circumnuclear dust in Circinus.

This paper is organized as such: in Sect. 2 we present the observations entering this work as well as the data reduction methods. In Sect. 3 we lay out the interferometric image reconstruction process and final image reconstruction parameters. We also compute image errors and assess the morphology of the resulting structure. In Sect. 4 we measure the temperature distribution of the dust in the central structure via blackbody fitting. In Sect. 5 we analyze the various components of the central dust structure in Circinus and discuss their implications. Finally, we conclude and summarize in Sect. 6.

2. Observations and data reduction

2.1. MATISSE observations

The MATISSE observations of Circinus were carried out on 13−14 March 2020, 27 Feb. 2021, and 31 May 2021 as part of guaranteed time observations. Data were taken with low spectral resolution in both the LM- (3−5 μm) and N-bands (8−13 μm). The observations were taken using the unit telescope (UT) configuration, with physical baselines ranging from 30 m to 140 m. At 12 μm this corresponds to angular resolutions between 9 and 41 mas with a “primary beam” of 153 mas. Each observation sequence consists of two sky exposures, a number of exposure cycles, N_cycles, consisting each of four 1 min interferometric exposures with different configurations of the beam commuting device (BCD) of MATISSE, as well as optional photometric exposures while chopping (for details see Lopez et al. 2022). Near the end of the night of 14 March 2020, we opted to repeat more exposure cycles to reduce the overhead time of re-acquisition on the target. The exact number of exposure cycles, along with the atmospheric conditions at the start of each observation, are given in Table 1. The observing conditions on 14 March 2020 were excellent, while on 13 March 2020 high-altitude cirrus negatively impacted acquisition, guiding, and adaptive optics in several individual exposures; we note that the final correlated flux error estimates on this night are higher. Observations on 28 Feb. and 01 Jun. 2021 were unaffected by such issues. We show the combined uv-coverage of all the observations in Fig. 1.

Fig. 1. MATISSE uv-coverage from all 25 exposure cycles. Squares denote observations taken in 2020, while circles represent observations from 2021. The mean correlated flux between 11.5 and 12.5 μm is used as the color scale. The discrete color binning is done in 0.2 Jy intervals, based on the measured correlated flux uncertainties.

On each night, we observed the calibration star HD 120404 (F_12 μm = 13 Jy) directly before and/or after the Circinus observations. The atmospheric conditions at the start of each calibrator observation are given in Table 1. This star serves a spectral calibrator, an instrumental phase calibrator, and an instrumental visibility calibrator. It has a MIR spectrum given by van Boekel (2004), and its diameter is given as 2.958 mas in Cruzalèbes et al. (2019). During the Feb. and May 2021 observations, we observed secondary calibrators, HD 120913 (F_12 μm = 5.7 Jy) and HD 119164 (F_12 μm = 1.2 Jy) in order to perform cross-calibration and closure phases accuracy checks.

We focus hereafter solely on the N-band observations. While LM-band data were recorded simultaneously, the low total flux of the core (458.16 ± 39.18 mJy; Isbell et al. 2021) results in very faint correlated fluxes for even a marginally resolved source. We leave the analysis of the LM data for a future paper, awaiting improvements in low-signal-to-noise calibration.

2.2. MATISSE data reduction and calibration

The N-band data were initially reduced using the MATISSE data reduction software¹ (DRS) version 1.5.1. We used the coherent reduction flags corrFlux = TRUE and coherentAlgo = 2 in order to produce correlated fluxes using the coherent integration algorithm as employed in the MIDI Expert Work Station (EWS; Jaffe 2004) and used in T14. We also use spectral binning 21 px (=1 μm) and the default values for all other parameters.

The correlated flux, F(u, v, λ) was then calibrated in the standard way:

(1)

where is the raw flux (in counts) of the calibrator, is the catalog flux of the calibrator, and is the raw flux of the target. This assumes the calibrator is unresolved; for the selected calibrators with diameter < 3 mas this is the case.

Within an observing cycle, individual exposures are taken minutes apart. The standard deviation of these correlated flux measurements is used as an uncertainty estimate, typically 0.2 Jy at all baselines. The uncertainties measured in this way broadly agree with the DRS-estimated values. The squared visibilities are finally calculated as , where the “zero-baseline” flux is is the arithmetic mean of the photometric flux spectra measured by each of the 8.2 m UTs.

The photometric flux was initially reduced via incoherent processing in the DRS (using corrFlux = FALSE). This mode extracts the photometric flux passing through the 2λ/D pinhole in each UT (0.61″ at 12 μm). This is not computed for each observing block, as the N-band photometry cycle adds 10 min to each observation, but once per epoch we record the photometry. The N-band photometry we obtain from the DRS is a factor ∼3 larger than expected from the MIDI and VISIR observations in T14, 36 ± 4 Jy vs. 12 ± 1 Jy at 12 μm. We doubt temporal flux variations in the source, as none of the correlated fluxes at any spatial scale exhibit a similar change since 2008 (see Sect. 2.5). When using EWS (Jaffe 2004), which was used previously for the MIDI observations, we extract a photometric flux of 12.4 ± 0.5 Jy at 12 μm for the same set of observations. This indicates that the photometric flux only exhibited a change due to the spatial filter used in each software; EWS employs a narrow Gaussian filter while DRS employs a wider top-hat. To compare consistently to the MIDI data, the EWS value is used.

We assume the calibration stars are symmetric and have zero closure phase on all phase triangles – any deviations from zero represent instrumental phase errors. As a first step in closure phase calibration, deviations from zero phase in the calibrator δϕ_⋆, ijk(λ) are subtracted from the target phase: . A typical MATISSE observation cycle includes 4 configurations of the BCD which serve to calibrate the closure phase. The varied BCD configurations (called out-out, in-in, in-out, out-in) should be identical save for sign flips on individual closure loops (as ϕ_ijk = −ϕ_ikj). We then average the star-calibrated closure phases. We first calculate the temporal mean value for each individual BCD configuration, as they are each repeated a number N_cycles times. Finally the mean of the four BCD configurations serves as the closure phase value at each wavelength, and the standard deviation is used as an estimate of our closure phase uncertainty (on the order of 15° for Circinus, on the order of 1° for HD 120404). We note, however, that all closure loops which include the ∼130 m baseline, UT1-UT4, have systematically higher uncertainties due to the low signal-to-noise ratio (S/N) correlated flux on this baseline. The uncertainty is on the order of 50° for Circinus, which means that only the closure triangles UT1-UT2-UT3 and UT2-UT3-UT4 provide high-precision phase information.

We have measurements in a total of 25 MATISSE exposure cycles, corresponding to 150 correlated flux measurements and 100 closure phase measurements. We define a position angle in the uv-plane as tanψ = v/u; we have sampled essentially all ψ between 0 and 110°, although the sampling is not uniform. This becomes especially noticeable on the longer baselines (> ​100 m). Two long-baseline regions at ψ ≈ [10, 40]° and ψ ≈ [80, 110]° are highly sampled, while a more sparse region is present between ψ = 45° and 60°. On the shorter baselines, no such gaps are present.

2.3. MIDI observations

We include short-baseline MIDI observations from T14. These short baselines provide the small spatial frequencies necessary for imaging or modeling of the large-scale structure in Circinus. These data were reduced using the MIDI Expert Work Station (EWS; Jaffe 2004). The exact procedure is given in T14. These data contain the correlated flux, the visibility amplitude, and the wavelength differential phase. We calculate the squared visibility as V² = (F_corr/F_tot)². Both the MATISSE and MIDI data have been calibrated with the same calibration star, HD 120404. The MIDI data do not provide closure phases, so we select only a small number of AT baselines rather than fully incorporating the MIDI uv-coverage. We selected the baselines to have (i) a projected baseline < 35 m; and (ii) u, v spacing of at least 8.2 m (the UT-diameter). This leaves us with 18 baselines from the small configuration. In the MATISSE OIFITS format, these 18 baselines correspond to 12 closure phase loops, which we give as 0 ± 180° such that these nonexistent closure phases have no weight on imaging. This assumption is supported by the closure phase measurements of the VLT spectrometer and imager for the mid-infrared (VISIR) sparse-aperture-masking data.

2.4. VISIR sparse-aperture-masking data

Circinus was observed with the sparse-aperture-masking (SAM) mode of VISIR. The observations were taken in the N-band (λ₀ = 11.3 μm; Filter Name = 11_3_SAM) on 02 June 2017 (099.B-0484A). The data consisted on five observing blocks on the science with interwoven observations with the calibrator star HD 125687. Each data set in the sequence SCI-CAL was observed with a DIT = 142 milliseconds and NDIT = 6 exposures. The data reduction consisted of two parts. The first one uses the ESOREX data reduction pipeline offered by ESO². It allowed us to correct for (i) the background, (ii) the bad pixels, (iii) to extract the interferograms from the chopping sequence and (iv) to center each frame on a 256 × 256 pixel grid. Frames with low signal-to-noise or with bad cosmetics were discarded manually from the data. Once the interferograms were cleaned, we extracted the interferometric observables from them.

To obtain the squared visibilities and closure phases from the data, we used the CASSINI-SAMPip³ software (see e.g., Sanchez-Bermudez et al. 2020). This algorithm fits the interferogram directly on the image plane, methods with similar performance based on fringe fitting are described by Greenbaum et al. (2015) and Lacour et al. (2011). The code uses a Single Value Decomposition method to obtain the interferometric observables. The algorithm works with monochromatic data and uses a sinc-filter for compensating the wavelength smearing of the broad-band VISIR filter. Each frame in the data was fitted independently. The uncertainties in the observables were obtained by averaging the observables of the six frames in each data set of science and calibrator, respectively. With the seven pin-holes mask available on VISIR, 21 squared visibilities and 35 closure phases were obtained per data set. The minimum baseline produced with the VISIR non-redundant mask has a length of 1.67 m (λ₀/2B_min = 600 mas) and the maximum one a length of 6.28 m (λ₀/2B_max = 184 mas), respectively. Figure C.1 shows, as example, one snapshot of the recorded interferogram of the science target and the uv-coverage obtained with our observations. Once the raw observables were extracted, the data were calibrated by dividing the squared visibilites of the target over the ones of the calibrator star; the closure phases were calibrated by subtracting the closure phases of the calibrator from the ones of the target. Figure C.2 shows the calibrated observables versus spatial frequency.

2.5. Correlated flux stability

Combining the MIDI and MATISSE datasets taken ≥10 years apart depends on the assumption that both the structure and photometry of Circinus are stable in the same period. T14 reported possible flux variation of Circinus between 2008 and 2009. Moreover, there may be instrumental biases which are not properly calibrated. Therefore, we compare the correlated flux values taken using MATISSE in 2020 and 2021 with those at similar u, v coordinates reported in T14. We identify and compare 30 baselines from MIDI and MATISSE which are within 4 m in u, v distance of each other; these are shown in Fig. 2. We find excellent agreement between the two epochs, with > 90% of baselines consistent within the 1σ_Fcorr ≈ 0.2 Jy calibrated correlated flux errors. Only two baselines are discrepant at 12 μm by > 1σ_Fcorr, but agree within 2σ_Fcorr. We find that there are no significant changes in correlated or total flux over the last ≥10 years.

Fig. 2. Comparison of MIDI and MATISSE correlated flux values on baselines cross-matched within 4 m. Color scale is difference in σFcorr from the MATISSE observations. Only two uv-points are 1σFcorr < ΔFcorr < 2σFcorr discrepant between the MIDI and MATISSE observations spaced more than 10 years apart.

2.6. Combination of MIDI and MATISSE data

Information at a large range of spatial frequencies is necessary for robust imaging of a source. The MATISSE UT observations have a shortest baseline of ∼30 m, which causes structures larger than 82.5 mas to be resolved out at 12 μm. Without MATISSE AT observations, we lack constraints on the large-scale structure. We know, however, that there is large-scale structure out to ≥600 mas from MIDI, VISIR-SAM, and VISIR data (T14, this work, and Asmus et al. 2014, respectively). In order to (a) avoid resolving out structure which is shown to be present, and (b) constrain the locations of small-scale structures, we perform the image reconstruction using a combination of MIDI and MATISSE data. The practice of including small spatial frequency data via modeling or data supplementation is common in imaging (e.g., Cotton 2017). The MIDI data do not include closure phase measurements, so as stated above we set the values to 0 ± 180° during imaging.

We claim that such a data supplementation is valid in the case of Circinus for the following reasons. First, both the MIDI modeling and VISIR-SAM data show that closure phases are small on these scales (ϕ ≤ 10° in the T14 modeling; ϕ = 0.1 ± 2.5° in the VISIR-SAM data). Secondly, it is safe to combine the squared visibility measurements directly, as we show in Sect. 2.5 that the fluxes are stable on all scales over the last 17 years. Finally, we can combine the AT and UT data despite their different inherent spatial filtering because at 30 m baselines, the MIDI AT and MATISSE UT data give consistent correlated flux values, indicating that they probe the same structure. We finally note that the 18 included MIDI AT baselines represent only a small fraction of the imaging data, and serve primarily as a spatial constraint. The results of imaging both with and without the AT data are described further in Appendix E, but in summary the primary small-scale features remain stable in either approach.

3. Image reconstruction

The primary advantage of MATISSE over MIDI is that the availability of closure phases makes it possible to reconstruct high-fidelity images. We employ the image reconstruction software, IRBis (Image Reconstruction software using the Bispectrum; Hofmann et al. 2014, 2016), which was designed for MATISSE and is incorporated into the standard data reduction package. IRBis includes six regularization functions, two minimization engines, and myriad fine-tuning parameters such as the pixel scale, hyperparameter, and object mask. For the VISIR-SAM data, a completely independent image reconstruction process was carried out. We kept the image reconstruction for the MATISSE+MIDI data and that for the VISIR-SAM data separate due to dynamical range concerns, the different wavelength ranges, and because the spacial scales they measure are completely independent. We first focus on the image reconstruction of the MATISSE+MIDI data, with MIDI closure phases assumed to be 0 ± 180° (see Sect. 2.6). The image reconstruction for the VISIR-SAM data will be discussed in Sect. 3.2.

3.1. MATISSE image reconstruction

We select seven wavelength bins in which to produce independent images: 8.5 ± 0.2 μm, 8.9 ± 0.2 μm, 9.7 ± 0.2 μm, 10.5 ± 0.3 μm, 11.3 ± 0.3 μm, 12.0 ± 0.2 μm, and 12.7 ± 0.2 μm. Any spectral information within each bin is averaged, producing a series of “gray” images. Each bin was imaged with a range of regularization functions and hyperparameters (hereafter μ; essentially a scaling on the amount of regularization), with the best selected via a modified χ² function:

(2)

with α and β serving as weights on either squared visibilities or the closure phases. In IRBis, there are three “cost functions” which vary the relative weighting of the closure phases and squared visibilities during the image reconstruction process. For cost function 1, α = β = 1; for cost function 2, α = 0 and β = 1. Cost function 3 is more complex, using the χ² coming from the sum of the bispectrum phasors and the squared visibilities (Hofmann et al. 2022); in essence replacing the closure phases in the second term of Eq. (2) with the bispectrum. We employed cost function 1 for the quality assessment of best-fitting images.

In order to produce images, we performed a grid search of the IRBis parameters, varying the field of view (FOV), the pixel number, the object mask, the regularization function, the hyperparameter μ, the cost function, and the reduction engine (ASA-CG or L-BFGS-B, see Hofmann et al. 2016 for more details). We use uniform weighting in the uv-plane, corresponding to weighting = 0 in IRBis. An initial best image is selected in each wavelength bin using Eq. (2), and a follow-up round of imaging using the best regularization function and pixel scale is performed. Regularization is a crucial component of ill-posed problems such as image reconstruction where the number of free parameters () is much larger than the number of data points. Regularization is the enforcement of an a priori constraint (e.g., smoothness, compactness, edginess, etc.) to prevent overfitting, but the strength of enforcement is set by the hyperparameter. Starting from the initial images, we finely vary the hyperparameter to construct L-curves – diagnostic comparisons between the amount of regularization and the residuals of the reconstruction (first applied by Lawson & Hanson 1995). One identifies the “elbow” of the curve as the image with optimal regularization parameters. This selection is necessary to strike a balance between over-regularization and allowing too many image artifacts to manifest. We give the final parameters for the reconstructions in Table 2. We note that different regularization functions in the same wavelength bin often result in very similar morphology, implying that the result is robust and simply not a consequence of regularized noise. Furthermore, the cost function has little effect on the final morphology or image quality and primarily aids convergence. We show the reconstructed images in Fig. 3, separating the continuum images from those inside the Si absorption feature. We also show the flux-weighted mean of the continuum images in Fig. 4 which represents an N-band image. Finally, we show the fit quality of each image in Figs. A.1 and B.1; we simulate the correlated fluxes and closure phases represented by each image at each uv-point and compare to the observed data. We see that overall the images trace the closure phase and correlated flux spectra well, although specific wavelengths at a handful of uv-coordinates are discrepant.

Fig. 3. Compilation of MATISSE images reconstructed independently in each wavelength bin. Pixel scale and field of view are matched in all panels. The FWHM of the beam is shown in the bottom right corner of each panel. The top row consists of images of the continuum emission, and the bottom row holds the images within the Si absorption feature. The bottom row also includes the reconstructed VISIR-SAM image (panel h), which has the field of view of the MATISSE images overlayed. All images are scaled to the power of 0.6. Contours are drawn at 5× the mean image error in each wavelength channel (see Sect. 3.1.1).

3.1.1. Image error analysis

We use the values in Table 2, which represent the “best” reconstruction parameters, to estimate the image-plane uncertainties. We do this through delete-d jackknife resampling of our uv-coverage (the method is developed in Shao & Wu 1989). In each Monte Carlo realization we randomly discard 10% of each the squared visibilities and the closure phases (i.e., 15 squared visibilities and 10 closure phases). This choice satisfies the criterion for being asymptotically unbiased: , where n is the sample size and d is the number of deleted elements. We then perform the image reconstruction at each wavelength using the parameters given in Table 2 and save the results. After 100 realizations, we calculate the median and standard deviation in each pixel of each image. The median image at each wavelength is used as the final image shown in Fig. 3. The standard deviation image serves as an error map with which we calculate the S/N at each pixel. The error maps and S/N maps are given in Appendix F. We use the median image at each wavelength for our morphological analysis. The patchiness of the extended structure at 12.0 and 12.7 μm is moreover confirmed through measurement of flux within a 14 mas aperture at several points in the polar emission. Taking the image errors into account, the differences between adjacent bright and dark regions (e.g., at [(51.6, 23.4), (51.6, 46.9)] mas and [(18.8, 37.5), (32.8, 37.5)] mas from the image center) are ≥2σ.

3.1.2. Morphology

In the final, independently reconstructed images (shown in Fig. 3), we find several consistent and key features. We discuss each below and have labeled them in Fig. 4 for reference.

Fig. 4. N-band continuum image and component labels. Left panel: flux-weighted mean MATISSE image – a proxy for an N-band continuum image. Contours are drawn at 5× the flux-weighted mean of the individual image errors. The cyan dashed ellipses represent the FWHM of the Gaussians fitted to MIDI observations of Circinus by T14. Center panel: same image as a contour map with levels at [5, 10, 20, 40, 80]× the 5× the flux-weighted mean of the individual image errors. Key morphological features are labeled: the 1.9 pc disk with i ≳ 83°, the polar emission, and the polar flux enhancement. Right panel: same image with the Greenhill et al. (2003) masers overplotted. The black dashed line represents the direction of the radio jet (Elmouttie et al. 1998). The cyan lines show the central PA and opening angle of the optical ionization cone (Fischer et al. 2013). Pixel scale and field of view are matched in both panels. All images and contours are scaled to the power of 0.65. The FWHM of the beam is shown in the bottom right corner of each panel. North is up and east is to the left.

A central disk-like component. This component is resolved in the NE-SW direction (≈1.9 pc across), but unresolved at all wavelengths in the NW-SE direction. Its orientation varies slightly in the different wavelength channels: along PA_disk ≈ 45° in the 8.5 and 8.9 μm images; and along PA_disk ≈ 30° in the images red-ward of the Si feature.

Central, unresolved flux. It is ≈10 mas (=0.2 pc) to the NE of the center of the disk in the 12 μm image. This is the brightest feature of the image at all wavelengths.

Significant extended emission in the polar direction (PA ∼ 295°), perpendicular to the maser emission and roughly aligned with the radio jet (see Fig. 4 for the orientations). The large-scale emission is more prominent at longer wavelengths. In the 11.3, 12.0, and 12.7 μm images, the extended emission is approximately symmetric about the photo-center, and it is roughly 4 × 1.5 pc. This emission is notably not smooth, and shows patchiness far above the noise level.

Two bright components, forming a rough line with the photo-center at PA_E-W  =   − 80° and superimposed on the polar emission, are observed for the first time. These substructures become more prominent at longer wavelengths, but are nonetheless present in all channels. They each extend to ∼65 mas (=1.2 pc) from the center and are both roughly 30 mas across at 12 μm.

At each wavelength we measure the flux contributions of the unresolved component, the disk, and the polar emission. These values are the total flux inside elliptical apertures placed at the image center with dimensions (10 × 10 mas), (100 × 10 mas) with PA = 45°, and (220 × 120 mas) with PA = −65°, respectively. The contributions from the disk and unresolved component are subtracted from the largest aperture. Similarly, the contribution from the unresolved component is subtracted from the disk aperture. These values are presented in Table 3. The fractional contribution of the unresolved component increases to shorter wavelengths, indicating that it contains relatively hot dust; conversely, the polar emission contribution increases at longer wavelengths because it is cooler.

There are several features which, while containing a large amount of flux, we consider to be artifacts of the image reconstruction process. A first, simple criterion is to take a S/N cut of σ_image ≥ 3, using the errors derived in Sect. 3.1.1. This simple cut agrees well with the following, more detailed considerations. If structures increase their size or radial distance from the photo-center linearly with wavelength, it is likely that they are artifacts of the uv-coverage. This is complicated, however, as different wavelengths probe different temperatures in the thermal infrared, and thus real structures may become “larger” at longer wavelengths where cooler dust is observed. An example of an artifact which varies with wavelength is the pair of arc-like emission features ∼100 mas to the NE and SW of the photocenter in the 9.7 μm image. These appear to correspond to the secondary peaks of the dirty beam (Fig. D.1). Finally, we assume that structures in the continuum should vary smoothly between adjacent imaging bins, and so we only consider those structures which are present in both the 8.5 and 8.9 μm images and those structures which are in all of the 11.3, 12.0, and 12.7 μm images to be real. We take the flux-weighted average of these five continuum images to produce a proxy for an N-band image. This is shown in Fig. 4. The continuum-average image emphasizes consistent features of the images while suppressing artifacts.

3.1.3. Effects of uv-coverage on image morphology

In this section we check what effects the attainable uv-coverage could have on our final images. On the two longest baselines, we have no uv-coverage for ψ > 110°, and even more notably, we have no uv-measurements on any baseline for 135°< ψ < 180°. These uv-holes are currently unavoidable due to VLTI delay line shadowing on the UT1–UT4 and UT2–UT3 baselines for Circinus at Dec = −65 : 20 : 21. While MATISSE can be used in a two-baseline configuration, we would not be able to measure the closure phases necessary for imaging. This uv-region has been shown by T14 to be important nonetheless, as the disk-like structure present in their modeling is primarily constrained by long baselines in this direction. T14 reports MIDI measurements of the UT1–UT3 baseline (∼90 m) in the uv-region we cannot currently measure. In order to test the effects of including measurements at these ψ, we performed a second round of imaging, incorporating MIDI baselines with BL ∈ [30, 100] m, ψ_MIDI ∈ [100, 180]°, and which are separated by at least 4 m in uv-coordinates. These criteria resulted in 18 additional baselines with correlated flux measurements from MIDI. As there were no closure phases for these baselines, we use the same procedure as when including the MIDI AT measurements, setting ϕ_T3, MIDI = 0 ± 180°.

Adding these 18 MIDI UT correlated fluxes to the 150 MATISSE measurements and the 18 MIDI AT measurements, we produced independent images at 8.9 and 12.0 μm. At 8.9 μm the resulting image is essentially unchanged, indiscernible by eye from the image shown in Fig. 3; an explicit comparison is shown in Fig. E.1. At 12.0 μm, however, the disk becomes more prominent and changes position angle slightly, while all other features remain constant. The disk-like structure in the initial imaging lies along PA_disk ∼ 35°, while after the addition of MIDI UT baselines the same disk-like structure lies along ∼40°. The latter value more closely resembles the 46 ± 3° given in T14. However, given the size of the beam at 12 μm, 9 mas, the disk orientation could quite easily vary in the image by ∼3°. The overall differences in the image plane are small when we include these baselines, so we proceed in our analysis without the MIDI UT measurements. It is clear from T14, however, that these baselines are important to understand the size and orientation of the disk-like structure in Circinus, and the planned doubling of the VLTI delay lines will make closure phase measurements including these baselines possible.

3.2. VISIR-SAM image reconstruction

The VISIR-SAM data allowed us to reconstruct an image of the target’s large scale at 11.3 μm. We reconstruct the VISIR-SAM image and the MATISSE images separately, rather than combining the uv-coverage for two reasons: first, the longest SAM baselines (≈7 m) are shorter than the shortest MIDI-AT baselines, meaning there is no overlap in measured visibilities; second, the SAM data exhibit squared visibilities > 0.4, which are much larger than the MATISSE values on longer baselines and would result in dynamical range issues during imaging. The data are nonetheless useful as a way to contextualize the MATISSE images and to measure the true extent of the polar emission. We used BSMEM (Buscher 1994; Lawson et al. 2004) to reconstruct the VISIR-SAM image. This code uses a regularized minimization algorithm to recover an image from infrared interferometric data. The regularized optimization engine uses a trust-region gradient-descent method with entropy (i.e., the sum of the logarithm of the pixel values in the image grid) as regularization function. Images were reconstructed using squared visibilities and closure phases simultaneously. The reconstructed image used a pixel scale of 15 mas in a pixel grid of 512 × 512 pixels. The code converges to a χ² close to unity. Figure 3 shows our VISIR-SAM image.

4. Measuring the dust temperature distribution

The images produced above supply not only morphological information, but also information about the temperature and optical depth of the dust in different regions. In this section, we fit one-temperature blackbody models with extinction to a series of apertures. As shown in Gámez Rosas et al. (2022), Gaussian modeling, point-source fitting, and image reconstruction all resulted in similar extracted SEDs in NGC 1068; therefore we can with confidence use the extracted apertures from our reconstructed image to undertake a temperature analysis of Circinus. We first convolve the images to the beam of the lowest resolution reconstruction (12.7 μm, corresponding to 10.1 mas). The individual images are aligned using cross-correlation before SED extraction, but effectively the photo-centers are simply matched. Then we define 13 apertures (shown in Fig. 5) which are 5 px (23.4 mas) in diameter and which do not overlap; their exact locations were chosen by hand to cover key features of the disk and polar emission. This is > 2× larger than the lowest resolution “beam” in our images, and in this way we do not make claims based on any hyper-resolved features. We extract the mean flux from each aperture in each image and estimate the flux error from the same apertures on the error maps estimated in Sect. 3.1.1. Finally, we add the calibration error of the total photometric flux at each wavelength in quadrature to the extracted flux error.

Fig. 5. Left panel: continuum image contours on top of the 8.9 μm image. Overplotted on these are the 13 apertures we used for SED extraction. The contours are at [5, 10, 20, 40, 80] times the continuum image error (see Fig. 4). All images and contours are scaled to the power of 0.65. Right panel: extracted mean flux at each wavelength in each aperture as well as the best-fitting blackbody functions to each SED. These fits had uniform priors on T and Gaussian priors on AV with μ = 28.5 mag and σ = 8.1 mag. The shaded regions represent 1σ uncertainties as estimated from the posterior probability distributions.

We fit a single blackbody (BB) curve with absorption to each aperture-extracted spectrum with the form⁴

(3)

where τ(λ)/τ_v = κ(λ)/κ_v, and we use the standard interstellar medium κ(λ) profile from Schartmann et al. (2005) which is based on the standard interstellar medium profile of Mathis et al. (1977). Fitting of T and A_V was done in two iterations using Markov chain Monte Carlo sampling with the package emcee (Foreman-Mackey et al. 2013). Final values in each iteration are the median of the marginalized posterior probability distribution. The 16th and 84th percentiles of the resulting temperature and extinction distributions are used as the 1σ fit uncertainties⁵.

In the first iteration, we use uniform prior probability distributions with T ∈ [100, 600] K and A_V ∈ [0, 100] mag. We find that the extinction does not vary significantly across the field. The central aperture, fit with the smallest uncertainties, shows mag which is using the mass-extinction profile from Schartmann et al. (2005). The other apertures have nominally higher values, but large uncertainties which make the differences insignificant. Only D_S40, W65, and E65 show differences > 1σ from the central value.

In the second iteration, we use again a uniform prior for temperature (T ∈ [100, 600] K), but a Gaussian prior for A_V with μ = 28.5 mag and σ = 8.1 mag based on the fit to D0 in the first iteration. The central aperture, D0, serves as a good estimate of the overall extinction because it (a) has the highest S/N and (b) has significant flux on both sides of the Si absorption feature. The resulting temperatures are consistent with the unconstrained case but are typically lower. The qualitative behavior of the temperature distribution is unchanged, but the fitted uncertainties are greatly diminished due to the degenerate nature of A_V and T for a fit to the Rayleigh-Jeans tail of the Planck function. In the following discussion, we therefore use the values from the second iteration. We show the best fit parameters for each aperture in each iteration in Table 4.

We do not find evidence of an extinction gradient across the field, indicating that there is a relatively uniform screen of foreground absorption. In the first fitting iteration, with A_V allowed to vary, the mean extinction values are similar to the east and west. In the second iteration, we restricted A_V around 28.5 mag, to get better constraints on the temperature. Based on HubbleK-band imaging, Wilson et al. (2000) estimated an extinction of A_V = 28 ± 7 toward a compact (< 2 pc) nucleus. Burtscher et al. (2016) measured a value of A_V = 27.2 ± 3 using SINFONI in the K-band. Roche et al. (2006) found 2.2 ≤ τ_9.7 ≤ 3.5 using T-ReCS on Gemini-South. Previous measurements are nearly identical to the fitted value in D0, mag, and furthermore consistent with the rest of the field. Uniform absorption, however, is in contrast to the Δτ = 27 arcsec⁻¹ gradient across the polar emission measured by T14. This discrepancy is puzzling, but we recognize there are major differences between our approach and that of T14. Specifically, T14 used differential phases and Gaussian modeling due to the lack of closure phase data. Their differential phases were measured on the UT and AT baselines, and thus probe larger-scale material than the MATISSE UT closure phases alone. On the other hand, we use no differential phases and had to assume an unconstrained AT closure phase value of 0 ± 180°. However, we note that on the UT baselines (probing ≲1 pc scales), the 9.7 μm closure phases are well matched by our images without an extinction gradient. The phase signals are instead produced by small-scale structure that was smoothed out in the Gaussian modeling approach of T14. The two approaches emphasize different aspects of the data, but differential phases could be included with the closure phases in future work in chromatic image reconstructions. Future closure phase measurements at 9.7 μm are required on shorter baselines (e.g., with MATISSE AT observations) to directly measure the Si absorption across the large-scale component.

We separate the apertures into two rough categories based on their locations. Those oriented NE and SW from the photocenter at PA ∼ 30° are labeled as “disk” apertures, based on the presence of a thin disk-like structure in both our reconstructed images and in the Gaussian modeling of T14. The other points, extending NW and SE from the photocenter are labeled “outflow” apertures, as they lie in the direction of the polar extension. The extracted spectra and the fitted blackbody curves (with uncertainty estimates as shaded regions) are shown in Fig. 5. The two-dimenstional temperature distribution based on the fits is shown in Fig. 6. We find that on average the “disk” apertures show a much steeper temperature falloff with projected distance than the “outflow” apertures.

Fig. 6. Two-dimensional temperature distribution as fitted in each of 13 apertures. Colors at the edge of each circle match those given by the aperture definitions in Fig. 5. Temperatures within the polar extension remain high even at large distances from the center when compared to the “disk” apertures. The contours are at [5, 10, 20, 40, 80] times the continuum image error and scaled to the power of 0.65 (see Fig. 4).

4.1. Temperature gradient analysis

In current modeling, the dust in the outflow is anisotropically illuminated by a face-on accretion disk. We can use the temperature profile of the outflow to characterize the dust environment.

We begin with a comparison to the simple analytic model of Barvainis (1987):

(4)

where L_acc is the luminosity of the accretion disk in L_⊙, r is the distance from the accretion disk in parsec, and τ_uv is the optical depth to the ultraviolet continuum. Here we use L_acc = 6 × 10⁹ L_⊙, which is the lower bound on estimates of the accretion disk luminosity in Circinus (6 × 10⁹ − 7 × 10¹⁰ L_⊙), as inferred from X-ray (Arévalo et al. 2014; Ricci et al. 2015), IR (T14) and optical (Oliva et al. 1999) observations. We also use A_V = 40 mag → τ_uv = 7.2, set roughly by the mean of the first-iteration fitted extinction values in Table 4 and converted using the dust extinction curve of Schartmann et al. (2005), but we note that the best-fitting value, mag, would result in even higher temperatures at a given radial distance. With these assumed values, we compare the radial temperature profile of the optically thin, continuous dust environment described by Barvainis (1987) to the fitted SED temperatures of the “outflow” in Fig. 7. At all radii, the Eq. (4) temperatures are larger than the measured Circinus temperatures by a factor of ∼2. This is not completely unexpected, as the inefficient re-radiation by the dust grains in the Barvainis (1987) model leads to higher temperatures at large radii; this model should be considered an upper limit on the dust temperatures at a given radius. Moreover, the Barvainis (1987) model does not take the anisotropy of the radiation into account. We also plot the expected temperature profile arising due to the simplest case of radiation equilibrium for perfectly efficient blackbodies as given in Tristram et al. (2007) for comparison.

Fig. 7. Radial temperature profile of the “torus” in Circinus based on 13 extracted apertures. Measurements along the disk-like structure in the center of image are marked with squares; measurements in the polar extensions are marked with diamonds. Colors correspond to the apertures given in Fig. 5. B+1987 is the radial profile from Barvainis (1987) and is given as the solid black line. The radial profile arising from simple radiation equilibrium is given as the dashed red line. In both analytic profiles the luminosity of the Circinus accretion disk is assumed to be Lacc = 6 × 109 L⊙. The shaded blue region labeled S+2008 represents the range of temperatures of the dust clumps at each radius in the standard clumpy torus model shown in Schartmann et al. (2008), Fig. 3. The hashed region represents the range of temperatures of dust cells in the disk+hyp model. The boxes/diamonds represent fitted blackbody temperatures in our 13 apertures applied to disk+hyp models based on Stalevski et al. (2019); each box (for disk apertures) or diamond (for polar apertures) center is the median temperature T and the height of each represents the range of temperatures fitted to the disk+hyp models.

4.2. Comparisons to radiative transfer models

4.2.1. Clumpy torus models

Modern AGN “torus” modeling takes the clumpiness of the dust, as implied from infrared interferometry, as well as anisotropic illumination from the accretion disk into account. We compare the temperatures at different radial distances in the standard clumpy torus model of Schartmann et al. (2008) to those fit in Circinus. These models consist of a wedge-shaped torus filled with randomly placed, optically-thick spheres of dust. The clump density falls off with radius, r, from the anisotropically illuminating source as ρ ∝ r^−0.5 and the clump size increases as a ∝ r^1.0. These models consist only of a puffy “disk” with half-opening angle θ = 45°, as they predate the observations of polar dust in Circinus. The clumpy torus models produce a range of dust temperatures as a function of radius which serve as a theoretical bound on the temperature distribution in the central few parsec. The temperatures found by our blackbody fits are clearly within these theoretical bounds of the model, cf. Fig. 7. A similar result was already found by Tristram et al. (2007, 2014).

4.2.2. Disk + wind models

More recently, Stalevski et al. (2017, 2019) undertook radiative transfer modeling of VISIR imaging data, the MIR SED, and MIDI interferometric data of Circinus. Their best-fitting model (presented in Stalevski et al. 2019) consists of a compact, dusty disk and a hollow hyperbolic cone extending in the polar direction (hereafter disk+hyp). In this modeling, a parameter grid for the radiative transfer models was searched such that the overall SED as well as the interferometric observables were well reproduced. This was not a model fit, but rather an exploration of the parameter space. For comparison with the MATISSE data, we started from the best model of Stalevski et al. (2019) and varied its parameters with finer sampling of the parameter space. We significantly expanded the explored range of the parameters which define the clumpiness: the number of clumps (i.e. filling factor) and different random realizations of the clumps’ positions (set by the “seeds” for the random number generator). Using the MATISSE uv-coverage, we simulate the squared visibilities and closure phases of each model image and compute the χ² to the data (the comparisons and resulting model are shown in Fig. G.1). This comparison placed constraints on the system inclination (i ∼ 85°), the hyperboloid opening angle (θ_OA ∼ 30°), the disk Si feature depth (τ_Si, DSK ∼ 15), and the outer radius of the disk (r_out ∼ 3 pc). The closure phase comparison favored a small number of clumps. We then performed a finer parameter search based on these constraints, focusing on the filling factors of the disk and the hyperboloid. After comparing with the MATISSE data, the parameter values defining the boundaries of the model geometry remain unchanged (the dusty disk outer radius, angular width, optical depth; hyperboloid shell position, width and optical depth). However, our modeling converged on a smaller number of clumps (30% less than in the MIDI model) and found that random positions of the clumps have a significant impact on the quality of the fits. The selected model exhibits a sky covering fraction of 78% due to the dust clumps at 0.55 μm. We show in Fig. 7 the average dust temperature as a function radius; these are indicative temperatures obtained by averaging the local thermal equilibrium temperatures over the dust species and grain sizes. We finally extract fluxes in each of the 13 apertures and fit blackbody temperatures using Eq. (3) to the disk+hyp model grid at λ ∈ [8.53, 8.91, 9.29, 9.70, 10.12, 10.56, 11.02, 11.50, 12.00, 12.52] μm.

We compare the extracted model spectrum in each aperture to the observed spectra. We quantify this through the χ² but do not perform any model fitting. These comparisons are shown in Fig. G.2. In the polar extension, the model and image extracted fluxes and temperatures agree well. The preferred model of Stalevski et al. (2019) includes the polar dust flux-enhancements E-W of the center. Along the disk, and particularly in the central aperture, D0, we see significant discrepancies. The central aperture temperature is ∼100 K less in the models than in the observations and the extracted flux is ≲10% of the observations. These discrepancies may indicate that the model disk is perhaps too dense. The disk apertures D_S77 and D_N77 also show much lower observed temperatures than the model predicts, indicating that the model disk can be further improved. Given that the outer radius, angular width, average edge-on optical depth and inclination of the disk appear to be well constrained, it is likely that the disk is actually inhomogeneous, or perhaps with a gap, thus allowing more warm emission to escape. LM-band images at ∼3 mas resolution are required to further constrain the disk component in modeling of Circinus. The MATISSE observations and imaging of Circinus agree very well with clumpy modeling, but it is beyond the scope of this work to place constraints on the specifics of a clumpy medium.

5. Discussion

MIR interferometry of Circinus has revealed several major components of the thermal dust: a disk-like central emission, large-scale polar emission, and a central point source along the disk. Image reconstruction has recovered these features in unprecedented detail and brought forward new substructures. The morphological features are labeled in Fig. 5. In the following, we examine each of these features separately. After the discussion of the individual features, from the smallest scales to the largest, we discuss the overall morphology.

The orientations of the central structures in Circinus are compared to those of the optical ionization cone and of the warped maser emission in the center. The well-studied optical ionization cone has a central axis along PA_opt = −52° and a projected half-opening angle between 36° and 41° (see, e.g., Marconi et al. 1994; Maiolino et al. 2000; Wilson et al. 2000; Fischer et al. 2013; Mingozzi et al. 2019). The observed ionized emission is not symmetric; it only extends toward the NW with no optical counterpart seen to the south, though a southern counterpart can be seen in the NIR (Prieto et al. 2004). The ionization cone is thought to coincide with an outflow of dense material, driven by radiation pressure and fed by a gaseous nuclear bar (Maiolino et al. 2000; Packham et al. 2005). Notably, the O[III] and Hα emission in the ionization cone is much brighter along its southern edge (PA ∼ −90°). The ionization cone is observed out to ∼40 pc from the nucleus (Wilson et al. 2000).

The warped H₂O maser disk was separated by Greenhill et al. (2003) into 3 components: the blueshifted emission (0.11 < r ≲ 0.4 pc; PA_{maser, blue} = 56 ± 6°), the central emission (0 < r < 0.11 ± 0.02 pc; PA_{maser, central} = 29 ± 3 deg), and the redshifted emission (0.11 < r ≲ 0.4 pc; PA_maser, red = 56 ± 6°). The central maser emission, which may trace the orientation of the accretion disk and the dense material around it, is nearly perpendicular to the radio jet axis (PA_jet = 115 and 295 ± 5°; Elmouttie et al. 1998), which is not aligned with the central axis of the optical ionization cone. These orientation markers are shown in Fig. 4 for comparison to the MATISSE images.

5.1. Unresolved central flux

We find a central, bright component unresolved at all wavelengths (≤6.7 mas at 8.5 μm – ≤10.1 mas at 12.7 μm). It is in the same position relative to the image photocenter in each reconstruction, and is therefore likely the same unresolved object present throughout. This point source is consistently found ≈​10 mas to the NE of the photocenter of the disk. Our central aperture, D0 (Fig. 5, Table 4), is centered on this point source and the extracted fluxes are brighter than the surrounding features by more than a factor of 2. We find that the fitted blackbody is relatively hot, K. While this source was well-fit by a single blackbody, we note that this is difficult to motivate physically and only serves as an estimate.

These results are similar to those of T14, who found a central unresolved component lying along the disk-like component. Their point source was shifted 14 mas to the NE of the disk-center, similar to the 10 mas which we find. T14 measured the temperature of this component to be 317 ± 22 K, which is ∼2σ lower than our measured temperature. The temperature difference is perhaps a result of the overlapping contributions of the three Gaussian components in T14, while we fit an isolated mean temperature at each extraction location. Nonetheless, no directly visible hot (≳900 K) dust is found by either T14 or this work.

The central aperture is almost certainly probing a column of much cooler dust along the line of sight and it may indeed reach dust at the sublimation temperature. A large range of spatial scales and temperatures are being merged into one aperture because of projection effects. It is thus difficult to draw any strong conclusions about the temperature in this feature without the LM-bands which should be more sensitive to hot dust. The L- and M-band fluxes measured using VLT/ISAAC by Isbell et al. (2021) represent the AGN flux within 630 mas, and are certainly an upper limit on the LM flux within the central aperture. However, if we perform a two-blackbody fit to the ISAAC LM measurements in addition to our central aperture fluxes, we see that a very compact and extincted 1500 K blackbody in addition to a larger 310 K blackbody fit the data very well. So, it is possible that the central aperture contains dust at the sublimation temperature, but we cannot draw any strong conclusions without fully analyzing the LM MATISSE data and spatially resolving the flux inside this central ∼10 mas region. This will be done in a subsequent work.

5.2. Central disk

We find a thin disk-like structure along ≈30°. It is present in all wavelength channels, but most prominent at longer wavelengths. In the continuum image, this disk is almost 1.9 pc in diameter and is unresolved in width. The extent of the disk is set by the 5σ contours at 12 and 12.7 μm. Considering the dirty-beam (Appendix D), we expect artifacts in the form of secondary lobes at ∼100 mas from the center along the disk PA, which indeed manifest themselves as low surface-brightness features near the edges of the images. Nonetheless, the central part of the disk in our images has a high flux density and is robustly detected at S/N ≥ 5.

Evidence that the dust in this disk is relatively dense comes from the blackbody fits performed on the “disk apertures”. Here we see that the temperature falls quickly as one moves farther from the photo center; indeed, the lowest fitted temperatures in the image occur in the disk at a projected distance of only 0.7 pc from the center. Taken together, the disk apertures (D0, D_S40, D_N40, D_S77, D_N72) exhibit a much steeper radial temperature gradient than apertures in the polar direction. The temperature profile of the disk is shown in both Fig. 7 and in the images themselves, as the disk becomes much less prominent at short wavelengths, indicating that the dust is relatively cool and the emission drops off significantly below 9 μm. The steep temperature gradient is possibly indicative of a dense environment wherein only the innermost dust has a direct view toward the accretion disk, and the outer clouds are heated only through re-radiation and photon scattering (e.g., Krolik 2007).

The disk component places constraints on the inclination of the system. Assuming that the disk is both thin and axisymmetric with diameter 1.9 pc, the fact that we do not resolve the width of the disk (≤9.5 mas = 0.18 pc at 12 μm) indicates an inclination i ≳ 83°. This is in agreement with the best disk+hyp model with i ∼ 85° matched to the closure phases and squared visibilities. This estimate can be considered a lower limit, as a more-realistic “puffed-up” disk would be thicker. ≳83° is closer to edge-on than on the galactic scale (∼65° Freeman et al. 1977; Elmouttie et al. 1998), the ALMA CO(3−2) tilted ring estimation (i ≳ 70° Izumi et al. 2018), and the T14 estimate for the MIR disk (i > 75°). However, Izumi et al. (2018) note that from 10 pc inward the inclination seems to increase, eventually reaching i ∼ 90° for the warped H₂O maser disk (Greenhill et al. 2003). The relatively dense dust of the MIR disk is consistent with the above and may lie in the same plane as the maser emission. The above assumes that the disk-like emission is the edge of the disk, rather than from reflected light on the top of a disk-like structure (see e.g., T14; Stalevski et al. 2019). We make this assumption because we observe no absorbing band on either side of the disk-emission.

The disk is aligned very well with the inner position angle of the H₂O maser emission (29 ± 3°; Greenhill et al. 2003) as well as with the compact nuclear disk (CND) at 10 s of pc found in ALMA CO(3−2) and [C I](1−0) (32 ± 1.9°; Izumi et al. 2018). The entirety of the warped maser disk, moreover, fits within our ≤9.5 mas = 0.18 pc thick dust disk. It is for this reason that we place the maser emission in the center of our disk; we do not have absolute astrometry from MATISSE, and so we must base the correspondence on the coincidence of PA and scale. Through Gaussian modeling, T14 also found a thin disk oriented along 46 ± 3° and with a FWHM of 1.1 ± 0.3 pc. The size of the disk in the T14 modeling is similar to what we measure. The T14 disk orientation differs slightly from that of our imaged disk, but they (a) used differential phases rather than closure phases in their modeling; and (b) used Gaussian modeling which simplifies the structure and may combine components. In Sect. 3.1.3 we found that with the T14 uv-coverage, our image disk could be oriented along ≈40°.

This dense disk of dust may play the role of the classical “torus”, obscuring a direct view toward the BLR. However, we find two competing phenomena. First, we see in the central aperture that hot dust is present, and depending on the exact distribution of the dust in the LM bands, we may even have a direct line of sight to dust at/near the sublimation temperature. This is, however, somewhat at odds with the steep temperature gradient we see across the disk. The thin disk must somehow be dense enough to shield some or most of the dust from directly seeing the sublimation zone or the central engine, but clumpy or low-density enough that we can see evidence of hot dust at or near the sublimation temperature. Authors such as Kishimoto et al. (2011) and Hönig et al. (2012) hypothesize that a “puffed-up” inner region (a few sublimation radii in size) may act as the classical obscuring torus.

5.3. Polar extension

We find a large-scale structure oriented in the same direction in all wavelength channels. This structure is referred to as a “polar extension” because its primary axis lies perpendicular to the AGN orientation and along the radio jet. The polar extension in our imaging is a large (∼4 × 1.5 pc) structure made up of warm (> 200 K) dust with major axis along ≈ − 60°. This larger envelope contains significant substructure: most prominently enhanced brightness directly E and W of the disk center. The polar emission exhibits “patchiness” at a significance ≥3σ on scales similar to the beam size, most prominently in the 12.0 and 12.7 μm images. Patchiness in the image could arise from clumpy dust emission, though it is unlikely that we resolve individual clumps at this scale (10 mas = 0.19 pc). Nonetheless, these images provide direct evidence that the polar emission is not a smooth, continuous structure.

We find that the substructures of the polar emission exhibit spatial variation in temperature. At a similar projected distance, apertures E65 and W65 are marginally hotter than W64 and E64 (∼270 K vs. ∼230 K). Additionally, the dust comprising the polar extended regions remains warm (∼200 K) out to a projected distance of ∼1.5 pc from the center of the structure. This behavior is significantly different than the dust temperature gradients along the disk, indicating differences in environment and density. The dust in the polar direction is likely less dense and/or more clumpy, as high temperatures at large distances require a relatively unobscured line of sight to the accretion disk. As shown in Fig. 7, the temperatures in the polar emission are entirely consistent with predictions from radiative transfer modeling of clumpy media (e.g., Schartmann et al. 2008; Stalevski et al. 2019), however only the latter reproduce the interferometric observables. At much lower resolution, the MIR SEDs of nearby AGN have shown that clumpy formalism is necessary to reproduce the relatively “blue” spectra indicative of an abundance of warm dust (e.g., Nenkova et al. 2008; Stalevski et al. 2016; Hönig & Kishimoto 2017). At the parsec scale, our results support the predictions of clumpy models.

5.3.1. E-W flux enhancements

The morphology we recover is in accordance with previous single-dish N-band estimates of the polar dust, and with the MIDI results of T14. The primary position angle of the polar extension was estimated from VLT/VISIR observations to be −80 ± 10° (Asmus et al. 2016). Similarly, the modeling done by T14 resulted in a mas FWHM (≈2 pc) Gaussian component with K and with a major axis along −73 ± 8 deg. Both the single-dish PA and that of the large Gaussian component in T14 are directed more closely to E-W orientation than our imaging suggests. This is likely explained by the lack of resolution and the simplicity of the Gaussian modeling; the large structure in our imaging shows significant nonuniformities. Namely, enhancements in flux directly to the E and to the W of the image photocenter. If one considers a flux-weighted mean of the polar emission in our imaging, it would certainly be more similar to the PA = 75 ± 8° as seen in T14. Indeed, the analysis by Stalevski et al. (2017, 2019) claims that the T14 large component is a simplified representation of an edge-brightened outflow cone, and they use this hypothesis to explain the discrepancy between the orientation of their polar outflow and the true pole of Circinus.

We present two possible explanations for the bright E-W substructure of the polar emission. The first is that the accretion disk in Circinus is tilted with respect to the central dust structures. If one considers that the central maser emission traces the orientation of the accretion disk (supported by the agreement with the radio jet position angle, assuming the jets originate in the central region), then one can relatively simply explain the asymmetric illumination of the polar extension. We show in Fig. 4 a line tracing the radio jet orientation. This line touches both the E and W flux-enhanced regions of the image. Due to the anisotropic nature of accretion disk emission, any dust some angle θ away from the “face” of the accretion disk is illuminated by a factor ∝ cos θ(2cosθ + 1) less than the dust which does see the “face” (Netzer 1987). The features we observe end more abruptly than this function suggests, but this could be due to patchiness or clumpiness of the dust. The idea of an accretion disk tilted with respect to the large-scale structures in Circinus is not new. Greenhill et al. (2003) suggests that the orientation of the accretion disk should only be “weakly coupled via gravity to the surrounding large-scale dynamical structures” because the central engine has a sphere of influence with a radius of only a few pc (Curran et al. 1998). Using VISIR images, MIDI observations, and the SED of Circinus, T14 as well as Stalevski et al. (2017, 2019) hypothesized that a warped or tilted accretion disk (as described by e.g., Petterson 1977; Nayakshin 2005) was required to asymmetrically illuminate the polar dust in their modeling. Hydrodynamic modeling of the central structures by Wada (2012) predicts that symmetric radiation-driven outflow cones should form perpendicular to the accretion disk. So while our observations suggest that the illumination of the polar dust is asymmetric – possibly from a tilted accretion disk – we cannot at this time explain why or how such a tilt occurred. The second possibility is that there is simply more material along the E-W direction; indeed Greenhill et al. (2003) speculated that the warped accretion disk could channel material in the nuclear outflow. This hypothesis is in better agreement with the Wada (2012) modeling, as in this case the polar outflows would be symmetric w.r.t. the accretion disk. The higher temperatures of the E-W flux-enhancements with respect to the apertures at the same projected distance argue in favor of the direct-illumination hypothesis. An overdensity of material should exhibit cooler temperatures due to dust self-shielding (as seen in the disk). This is merely a qualitative agreement, and in order to distinguish between these two hypotheses, detailed modeling of the formation of the outflow cones in the presence of a warped accretion disk will be crucial.

5.3.2. Connection to larger scales

It is clear in the MATISSE imaging that the majority of thermal dust emission in the center of Circinus comes from the polar extension, but its full extent is poorly constrained. In our imaging, any structure larger than those probed by the shortest baselines is resolved out; this means for imaging using the MIDI AT baselines we are not sensitive to structures larger than 688 mas at 12 μm. This is strictly an upper limit, however, and in the image reconstruction process we (a) limit the FOV to 600 mas, and (b) apply an object mask with a radius 160 mas. The object masking heavily suppresses any structure which falls outside of the specified radius. We can, nonetheless, confidently state that there is N-band emission out to ∼1.5 pc from the center to both the NW and SE, and that the emission shows a flux enhancement to the E and W of the image center.

The VISIR-SAM data were fit in the image plane with a Gaussian having FWHM 3.3 × 2.2 pc and major axis along PA = 72°. This is larger than either the MATISSE images or the MIDI modeling (with FWHM = 2 pc), indicating that the MATISSE images do not capture the true extent of the structure. The position angle of the SAM data matches the T14 result, though it is likely also flux-biased toward the South due to the E-W flux enhancements. Continuing to lower resolution, the N-band VLT/VISIR images in Asmus et al. (2014, 2016) show that in Circinus roughly 60% of the flux is extended farther than 5.24 pc and at PA = 100 ± 10°. It is clear that the polar structures we see in our images extend continuously outward past 5 pc.

5.4. Overall morphology

We present the first model-independent image of the circumnuclear dust in Circinus. The recovered combination of a geometrically thin disk and large-scale polar emission supports previous MIR interferometric findings, but newly imaged substructures hint at complexity unmatched in existing modeling. In particular, we find that the disk is simultaneously dense and yet allows emission from hot dust to radiate through; we find that an unresolved component lies 10 mas NE of the photocenter along the disk; and we find significant flux enhancements in the polar emission E and W of the disk-center.

The size of circumnuclear dust structures has been shown to vary with AGN luminosity (e.g., Kishimoto et al. 2011; Burtscher et al. 2013). The scales measured herein of the circumnuclear structures in Circinus – namely a thin disk with diameter 1.9 pc and ≳4 pc polar emission – with L_AGN = 6 × 10⁹ − 7 × 10¹⁰ L_⊙ (Arévalo et al. 2014; Ricci et al. 2015; Tristram et al. 2014; Oliva et al. 1999) place a constraint on the luminosity-dependent scaling of the dust structures in AGN. Leftley et al. (2019) showed that the ratio of extended flux to unresolved flux increased with Eddington ratio (ϵ_Edd), claiming that this implied the presence of more dust in a radiation-driven wind for a higher ϵ_Edd. Circinus, with ϵ_Edd ∼ 0.2 (Greenhill et al. 2003), is dominated by polar dust emission. We measure the flux of the unresolved component to be F_{pt, 12 μm} = 0.77 ± 0.04 Jy, which is 6.3 ± 0.3% of the total flux at 12 μm. At 8.9 μm we measure F_{pt, 8.9 μm} = 0.39 ± 0.01 Jy, which is 9.0 ± 0.2% of the total flux. The fraction at 12 μm is significantly smaller than previously reported (20% and 10% at 12 μm in Leftley et al. 2019; López-Gonzaga et al. 2016, respectively), but they relied on simple two-Gaussian modeling of MIDI data.

Disk+wind radiative-transfer models (Hönig & Kishimoto 2017; Stalevski et al. 2019) have recently been invoked to explain the polar emission found in a number of nearby AGN (e.g., Tristram et al. 2007, 2014; Burtscher et al. 2013; López-Gonzaga et al. 2014, 2016; Leftley et al. 2018). Fits to the NIR and MIR SEDs of nearby AGN have shown that disk+wind models provide the best match to the overall SED, reproducing the MIR flux through large-scale emission and NIR flux via reflected light from the accretion disk in the windy outflow (e.g., Martínez-Paredes et al. 2020; Isbell et al. 2021). The disk+wind morphology in the radiative-transfer models is supported by hydrodynamical and radiation-hydro modeling (Wada 2012; Wada et al. 2016; Williamson et al. 2020; Venanzi et al. 2020), but has had few direct observational constraints. The images presented in this work, with a thin disk (1.9 pc × ≤ 0.18 pc) and polar emission (∼4 × 1.5 pc) perpendicular to it, resemble the disk+wind models only in broad strokes. Modifications to the disk+wind model in Stalevski et al. (2019) explain the E-W flux enhancements in the polar emission via a tilted accretion disk, but the dynamical stability of such a shift in radiation pressure remains untested. Hydrodynamical models produce structures symmetric about the accretion disk (Wada 2012; Venanzi et al. 2020), so tilting with respect to the dusty structures may play a larger role. Whether this is specific to Circinus or a more general feature remains to be explored.

Only one other AGN has been imaged with MATISSE so far: NGC 1068. Imaging work by Gámez Rosas et al. (2022) has revealed a quite different circumnuclear dust morphology than we recover. In NGC 1068 at 12 μm, they find a disk-like structure ∼2 pc in diameter with emission extending nearly perpendicular to it, similar to what we see with the disk and E-W polar flux enhancements. However, the NGC 1068 and Circinus morphologies differ significantly at other wavelengths. At 8.5 μm and in the LM-bands, the NGC 1068 emission is resolved into a ring-like structure with 720 K dust embedded within. We have shown that hot dust can make up the unresolved flux in Circinus, and the LM data can help clarify the situation, as they probe the ≳500 K dust morphology. Finally, Jaffe et al. (2004), López-Gonzaga et al. (2014), and Gámez Rosas et al. (2022) showed that in NGC 1068, the standard ISM dust we use does not reproduce the observed SEDs. The effects of varying dust composition will be explored in future work.

Future N-band observations with the MATISSE ATs will yield the first closure phase measurements of the ≳1 pc dust, further improving our imaging capabilites beyond the MIDI data. In a subsequent paper, we will utilize the LM-band MATISSE data, with ∼3 mas resolution, to probe the hotter dust at small scales both within the disk and perhaps at the origins of the polar extension. We show that the central aperture can contain dust near the sublimation temperature, and a detailed study of the LM-data can give insights into the density, thickness, and perhaps the clumpiness of the disk. Circinus and NGC 1068 are laboratories in which to study the circumnuclear dust at extremely small physical scales, but the ongoing the MATISSE AGN Programme aims toward a statistical understanding of the central dust scaling and relation to the SMBH.

6. Conclusions

In this work we present the first images of the circumnuclear dust in the Seyfert 2 galaxy Circinus. These images were reconstructed with IRBis using 150 correlated fluxes and 100 closure phases in the N-band from VLTI/MATISSE. Closure phase measurements of Circinus are reported here, and their novel inclusion in MATISSE observations makes imaging possible for the first time. The above results are largely in agreement with previous observations from MIDI (Tristram et al. 2007, 2014) and VISIR (Asmus et al. 2014). But our images, moreover, are model-independent and show new substructure which can be used to further constrain physical modeling of circumnuclear dust in AGN. Through analysis of the interferometric observables and the images reconstructed in seven independent wavelength channels we

1. Show that correlated flux measurements on individual baselines have not changed over the last 17 years, implying that the underlying structures remain unchanged from the MIDI observations obtained between 2004 and 2011.

2. Find significant substructure in the circumnuclear dust. The circumnuclear dust can be separated into several components: central, unresolved flux; a thin disk 1.9 pc in diameter; polar emission (∼4 × 1.5 pc) extending orthogonal to the disk and exhibiting patchiness; and flux enhancements E and W of the disk embedded within the polar dust.

3. Report that the polar dust makes up ∼60% of the total flux, increasing toward longer wavelengths. The unresolved flux makes up ≲10%, increasing toward shorter wavelengths and further hinting at the presence of hot dust.

4. Measure SEDs in 13 apertures across the structures and fit temperature and extinction values to blackbodies in those apertures. We fit hotter dust temperatures ( K) in the central aperture along with warm dust (T ≳ 200 K) 1.5 pc from the center, indicating a clumpy circumnuclear medium. We clearly distinguish the radial temperature profiles of the disk and the polar extension: the disk shows a steeper temperature gradient indicating dense material; the polar emission shows a much flatter temperature profile with warm temperatures out to 2 pc from the center.

5. Recover a remarkably symmetric object, in terms of both flux and temperature distributions. We fit mag, consistent with the galactic-scale value (A_V = 28 ± 7 Wilson et al. 2000). We find no evidence of an absorption gradient across the field, in contrast to previous results (i.e., Tristram et al. 2014). Our new results indicate the presence of a foreground dust screen with very little local variation.

6. See that on large scales, the recovered morphology of the N-band dust in Circinus resembles the results of disk+wind modeling (e.g., Wada et al. 2016; Stalevski et al. 2019), but new questions are raised because the subparsec dust is imaged here for the first time. We find that the temperature distribution is well-reproduced by the clumpy torus models of Schartmann et al. (2008) and Stalevski et al. (2019). The Schartmann et al. (2008) models do not, however, match the imaged morphology. The disk+hyp models better match the structure, but discrepancies are found in the central and disk apertures, indicating modifications to the disk component are necessary in the models. Using a suite of disk+hyp models based on Stalevski et al. (2019), we find that a large range of clump densities and disk filling-factors can match the data within the uncertainties of the images and interferometric observables.

7. Discover inhomogeneities in the polar dust emission: namely significant patchiness on scales of the resolution element; and flux enhancements directly to the E and W of the disk. The here-discovered patchiness is the first direct evidence that the polar dust is not a smooth, continuous structure but is rather clumpy. The E-W flux enhancement raises questions about the relation of the accretion disk to the larger dust structures.

The imaged substructures and temperature distributions presented herein serve as a direct constraint on future physical modeling of the circumnuclear dust. In the disk+wind model, the thin disk we image is related to inflowing material, while the polar emission represents a radiation-driven outflow. How these components relate to large-scale (≳10 pc) structures and furthermore to the host galaxy can be tested in both hydrodynamical modeling and future observations, specifically with the MATISSE ATs. It is clear that the classic geometrically-thick torus is not present in our imaging, but the (nearly) rotationally symmetric structures we recover can play much the same role; yet the detailed implications for AGN Unification remain to be explored through modeling and the MATISSE AGN Programme.

Acknowledgments

We thank the VLTI Paranal staff for all their help and good company during observations, the entire MATISSE Consortium for their support, and the anonymous referee for their helpful comments which have improved the manuscript. J.W.I. thanks Felix Bosco and Thomas Jackson for insightful discussions. J.S.B. acknowledges the full support from the UNAM PAPIIT project IA 101220 and from the CONACyT “Ciencia de Frontera” project 263975. M.S. acknowledges support by the Ministry of Education, Science and Technological Development of the Republic of Serbia through the contract no. 451-03-9/2022-14/200002 and by the Science Fund of the Republic of Serbia, PROMIS 6060916, BOWIE. This research has made extensive use of NASA’s Astrophysics Data System; the SIMBAD database and VizieR catalogue access tool, operated at CDS, Strasbourg, France; and the python packages astropy, emcee, scipy, and matplotlib.

References

Appendix A: MATISSE correlated fluxes

In Fig. A.1 we present the correlated flux for each uv-point, reduced and calibrated as described in Sect. 2. The total photometric flux (“the zero-baseline flux”) is included in the first panel of Fig. A.1.

Fig. A.1. Correlated flux data for Circinus from March 2020 (blue), February 2021 (yellow), and May 2021 (red). The black points are simulated values extracted from the final images, with errors estimated using the 1σ error maps (described in Sect. 3.1.1). The total photometric flux is included in the first panel. Presented errors come from both the calibrator flux uncertainty and the statistical variation of the observables within a set of observing cycles. Near 8 and 13 μm one can see flux variations due to the edges of the atmospheric window.

Fig. A.1. continued.

Fig. A.1. continued.

Fig. A.1. continued.

Appendix B: MATISSE closure phases

In Fig. B.1 we present the closure phase for each closure triangle, reduced and calibrated as described in Sect. 2.

Fig. B.1. Closure phase data for Circinus from March 2020 (blue), February 2021 (yellow), and May 2021 (red). Presented errors come from both the calibrator phase uncertainty and the statistical variation of the observables within a set of observing cycles. The black points are simulated values extracted from the final images, with errors estimated using the 1σ error maps (described in Sect. 3.1.1). The panels are sorted by length of the longest projected baseline in the closure triangle.

Fig. B.1. continued.

Fig. B.1. continued.

Fig. B.1. continued.

Appendix C: VISIR-SAM data

In Figs. C.1 and C.2 we show the measured VISIR-SAM data from the reduction as described in Sect. 2.4.

Fig. C.1. u-v coverage obtained with the VISIR-SAM data (left), the different colors indicate the 21 different baselines in the data. Snapshot of the Circinus interferogram obtained with the VISIR-SAM data (right).

Fig. C.2. VISIR-SAM data used in the diagnosis of large-scale closure phases. In the (left) panel we show squared visibilities. The individual observations are in gray, and the mean value over the cycles at a given baseline is in red with 1σ error bars coming from the standard deviation of the cycles. In the (right) panel we show the same for closure phase, with mean values in green.

Appendix D: Dirty beam

We estimate the dirty beam in the typical way in order to identify image artifacts. In the uv-plane, we set the squared visibility at each uv-point we observed (±4.05 m, the UT radius) to 1 and the surrounding points to 0. We set the phase to 0 deg across the uv-plane. We finally take the inverse Fourier transform of this complex array to obtain an estimate of the dirty beam (shown in Fig. D.1).

Fig. D.1. Dirty beam estimated for the combined MIDI AT and MATISSE UT uv-coverage. On the (left) we show the final uv-coverage, and on the (right) we show the resulting dirty beam with square-root scaling inside a 600 mas window.

Appendix E: Imaging with and without ATs

In Fig. E.1 we show the effects of imaging with and without the MIDI AT baselines. We stress that the MIDI AT baseline inclusion is necessary due to the resolved nature of this AGN, as shown in both MIR interferometric and single-dish observations. The MIDI AT baselines require the synthesis of closure phase triangles in order to match the IRBis formatting. We set the closure phases involving these baselines to 0 ± 180°, such that they do not bias the imaging. We justify the inclusion of these baselines through the following arguments: First, the correlated flux values for all 30 MIDI uv-points within 4m of a MATISSE point show < 2σ variation over 10 years (Sect. 2.5). Second, the AT baselines from MIDI transition continuously to the MATISSE UT baselines around 30m (i.e., variations within the 0.2 Jy correlated flux uncertainties). Finally, VISIR-SAM imaging of Circinus shows 0.1 ± 2.5° closure phases on ≤6.3 m baselines (Sect. 2.4). This agrees with T14’s Gaussian modeling of the MIDI data which gives ≈0° closure phases for baselines ≤30 m.

Fig. E.1. Comparison of image reconstruction using different uv-samplings. In the (leftmost panel) we show the image resulting from MATISSE UT uv-coverage alone. In the (second panel) we show the MIDI UT+AT Gaussian model from Tristram et al. (2014). In the (third panel) we show the image reconstruction resulting from the combination of MATISSE UT and MIDI AT uv-coverage. In the (rightmost panel) we show the results of imaging using the MATISSE UT, the MIDI AT, and the MIDI UT data, with closure phases in the MIDI data set by the T14 Gaussian model. The interior structures (a disk, an unresolved source, and bright E-W flux enhancements) are present in all image reconstructions, implying their fidelity.

Nonetheless, it is instructive to see which structures arise as a result of the MATISSE-only imaging. We show the 12 μm UT-only reconstruction in Fig. E.1 alongside the Gaussian model of T14 which used both UTs and ATs from MIDI (without closure phases) and the 12 μm image reconstruction as detailed in Sect. 3. We see that the central ≈1 pc is nearly identical in the two images, and notably the bright features E-W of the center remain prominent in both setups. The disk-like component is perhaps even more obvious in the UT-only image, given the same color scaling. The largest difference between the images is the lack of large-scale extended flux, but this is expected as the UTs shortest baseline corresponds to ≈40 mas at 12 μm, and structures larger than this are suppressed.

Appendix F: Image error estimates

We performed delete-d jackknifing (Shao & Wu 1989) to estimate the errors present in our final images (see Sect. 3.1.1). We present the final images, the error maps, and the S/N maps in Figs. F.1–F.7. We use the S/N maps to determine which morphological features we trust. We perform an S/N cut of ≥3 on the final images to (a) define where valid apertures can be located, and (b) determine the extent of large features.

Fig. F.1. Final image reconstructions and error estimates. From left to right: the median image, error map, and S/N map for 8.9 μm as estimated from delete-d jackknifing. The white contour is S/N = 3 and the black contour shows S/N = 5.

Fig. F.2. As Fig. F.1, but for 8.9 μm.

Fig. F.3. As Fig. F.1, but for 9.7 μm.

Fig. F.4. As Fig. F.1, but for 10.5 μm.

Fig. F.5. As Fig. F.1, but for 11.3 μm.

Fig. F.6. As Fig. F.1 but for 12.0 μm.

Fig. F.7. As Fig. F.1, but for 12.7 μm.

Appendix G: SKIRT model parameter variation

In Fig. G.1 we compare the extracted squared visibilities and closure phases to the observed values for a range of model parameters. The simulated squared visibilities and closure phases use the uv-coverage of the MATISSE UTs. The parameters are the depth of the silicate feature in the disk and in the hyperboloid (τ_9.7), the outer radius of the disk (Disk R_out), the opening angle of the hyperboloid, the relative number of clumps (N_cl), and the inclination of the model (i where 90° is edge-on). The comparisons place constraints on the system inclination (i ≈ 85°), the hyperboloid opening angle (θ_OA ≈ 30°), the disk Si feature depth (τ_Si, DSK ≈ 14), and the outer radius of the disk (r_out ≈ 3 pc). The closure phases provide clearer constraints.

Fig. G.1. Comparisons via χ2 of simulated observables – squared visbilities and closure phases – to the MATISSE data for a range of disk+hyp model parameter values. In the top six panels we show the squared visibility comparisons. In the middle six panels we show the closure phase comparisons. In the bottom panels we show the model with parameters favored by the χ2 comparison at both its native and 10 mas resolution.

In Fig. G.2 we show a comparison via χ² between the model spectra and observed spectra in each of our 13 apertures defined in Sect. 4. We see that flux in the central apertures is under-represented in the models, indicating that modifications to the disk component (of e.g., clumpiness or thickness) may be necessary.

Fig. G.2. Comparisons of measured aperture-extracted spectra to those of disk+hyp models with disk filling factor varied. The displayed χ2 in each panel is the mean χ2 of all models and the ranges are given by the standard deviation of the model values. At large radii, the models agree well with observations, but the unresolved central flux is under-represented in the models.

All Tables

All Figures

	Fig. 1. MATISSE uv-coverage from all 25 exposure cycles. Squares denote observations taken in 2020, while circles represent observations from 2021. The mean correlated flux between 11.5 and 12.5 μm is used as the color scale. The discrete color binning is done in 0.2 Jy intervals, based on the measured correlated flux uncertainties.
In the text

	Fig. 2. Comparison of MIDI and MATISSE correlated flux values on baselines cross-matched within 4 m. Color scale is difference in σFcorr from the MATISSE observations. Only two uv-points are 1σFcorr < ΔFcorr < 2σFcorr discrepant between the MIDI and MATISSE observations spaced more than 10 years apart.
In the text

Fig. 3. Compilation of MATISSE images reconstructed independently in each wavelength bin. Pixel scale and field of view are matched in all panels. The FWHM of the beam is shown in the bottom right corner of each panel. The top row consists of images of the continuum emission, and the bottom row holds the images within the Si absorption feature. The bottom row also includes the reconstructed VISIR-SAM image (panel h), which has the field of view of the MATISSE images overlayed. All images are scaled to the power of 0.6. Contours are drawn at 5× the mean image error in each wavelength channel (see Sect. 3.1.1).

In the text

Fig. 4. N-band continuum image and component labels. Left panel: flux-weighted mean MATISSE image – a proxy for an N-band continuum image. Contours are drawn at 5× the flux-weighted mean of the individual image errors. The cyan dashed ellipses represent the FWHM of the Gaussians fitted to MIDI observations of Circinus by T14. Center panel: same image as a contour map with levels at [5, 10, 20, 40, 80]× the 5× the flux-weighted mean of the individual image errors. Key morphological features are labeled: the 1.9 pc disk with i ≳ 83°, the polar emission, and the polar flux enhancement. Right panel: same image with the Greenhill et al. (2003) masers overplotted. The black dashed line represents the direction of the radio jet (Elmouttie et al. 1998). The cyan lines show the central PA and opening angle of the optical ionization cone (Fischer et al. 2013). Pixel scale and field of view are matched in both panels. All images and contours are scaled to the power of 0.65. The FWHM of the beam is shown in the bottom right corner of each panel. North is up and east is to the left.

In the text

Fig. 5. Left panel: continuum image contours on top of the 8.9 μm image. Overplotted on these are the 13 apertures we used for SED extraction. The contours are at [5, 10, 20, 40, 80] times the continuum image error (see Fig. 4). All images and contours are scaled to the power of 0.65. Right panel: extracted mean flux at each wavelength in each aperture as well as the best-fitting blackbody functions to each SED. These fits had uniform priors on T and Gaussian priors on AV with μ = 28.5 mag and σ = 8.1 mag. The shaded regions represent 1σ uncertainties as estimated from the posterior probability distributions.

In the text

Fig. 6. Two-dimensional temperature distribution as fitted in each of 13 apertures. Colors at the edge of each circle match those given by the aperture definitions in Fig. 5. Temperatures within the polar extension remain high even at large distances from the center when compared to the “disk” apertures. The contours are at [5, 10, 20, 40, 80] times the continuum image error and scaled to the power of 0.65 (see Fig. 4).

In the text

Fig. 7. Radial temperature profile of the “torus” in Circinus based on 13 extracted apertures. Measurements along the disk-like structure in the center of image are marked with squares; measurements in the polar extensions are marked with diamonds. Colors correspond to the apertures given in Fig. 5. B+1987 is the radial profile from Barvainis (1987) and is given as the solid black line. The radial profile arising from simple radiation equilibrium is given as the dashed red line. In both analytic profiles the luminosity of the Circinus accretion disk is assumed to be Lacc = 6 × 109 L⊙. The shaded blue region labeled S+2008 represents the range of temperatures of the dust clumps at each radius in the standard clumpy torus model shown in Schartmann et al. (2008), Fig. 3. The hashed region represents the range of temperatures of dust cells in the disk+hyp model. The boxes/diamonds represent fitted blackbody temperatures in our 13 apertures applied to disk+hyp models based on Stalevski et al. (2019); each box (for disk apertures) or diamond (for polar apertures) center is the median temperature T and the height of each represents the range of temperatures fitted to the disk+hyp models.

In the text

Fig. A.1. Correlated flux data for Circinus from March 2020 (blue), February 2021 (yellow), and May 2021 (red). The black points are simulated values extracted from the final images, with errors estimated using the 1σ error maps (described in Sect. 3.1.1). The total photometric flux is included in the first panel. Presented errors come from both the calibrator flux uncertainty and the statistical variation of the observables within a set of observing cycles. Near 8 and 13 μm one can see flux variations due to the edges of the atmospheric window.

In the text

	Fig. A.1. continued.
In the text

	Fig. A.1. continued.
In the text

	Fig. A.1. continued.
In the text

Fig. B.1. Closure phase data for Circinus from March 2020 (blue), February 2021 (yellow), and May 2021 (red). Presented errors come from both the calibrator phase uncertainty and the statistical variation of the observables within a set of observing cycles. The black points are simulated values extracted from the final images, with errors estimated using the 1σ error maps (described in Sect. 3.1.1). The panels are sorted by length of the longest projected baseline in the closure triangle.

In the text

	Fig. B.1. continued.
In the text

	Fig. B.1. continued.
In the text

	Fig. B.1. continued.
In the text

	Fig. C.1. u-v coverage obtained with the VISIR-SAM data (left), the different colors indicate the 21 different baselines in the data. Snapshot of the Circinus interferogram obtained with the VISIR-SAM data (right).
In the text

	Fig. C.2. VISIR-SAM data used in the diagnosis of large-scale closure phases. In the (left) panel we show squared visibilities. The individual observations are in gray, and the mean value over the cycles at a given baseline is in red with 1σ error bars coming from the standard deviation of the cycles. In the (right) panel we show the same for closure phase, with mean values in green.
In the text

	Fig. D.1. Dirty beam estimated for the combined MIDI AT and MATISSE UT uv-coverage. On the (left) we show the final uv-coverage, and on the (right) we show the resulting dirty beam with square-root scaling inside a 600 mas window.
In the text

Fig. E.1. Comparison of image reconstruction using different uv-samplings. In the (leftmost panel) we show the image resulting from MATISSE UT uv-coverage alone. In the (second panel) we show the MIDI UT+AT Gaussian model from Tristram et al. (2014). In the (third panel) we show the image reconstruction resulting from the combination of MATISSE UT and MIDI AT uv-coverage. In the (rightmost panel) we show the results of imaging using the MATISSE UT, the MIDI AT, and the MIDI UT data, with closure phases in the MIDI data set by the T14 Gaussian model. The interior structures (a disk, an unresolved source, and bright E-W flux enhancements) are present in all image reconstructions, implying their fidelity.

In the text

	Fig. F.1. Final image reconstructions and error estimates. From left to right: the median image, error map, and S/N map for 8.9 μm as estimated from delete-d jackknifing. The white contour is S/N = 3 and the black contour shows S/N = 5.
In the text

	Fig. F.2. As Fig. F.1, but for 8.9 μm.
In the text

	Fig. F.3. As Fig. F.1, but for 9.7 μm.
In the text

	Fig. F.4. As Fig. F.1, but for 10.5 μm.
In the text

	Fig. F.5. As Fig. F.1, but for 11.3 μm.
In the text

	Fig. F.6. As Fig. F.1 but for 12.0 μm.
In the text

	Fig. F.7. As Fig. F.1, but for 12.7 μm.
In the text

Fig. G.1. Comparisons via χ2 of simulated observables – squared visbilities and closure phases – to the MATISSE data for a range of disk+hyp model parameter values. In the top six panels we show the squared visibility comparisons. In the middle six panels we show the closure phase comparisons. In the bottom panels we show the model with parameters favored by the χ2 comparison at both its native and 10 mas resolution.

In the text

	Fig. G.2. Comparisons of measured aperture-extracted spectra to those of disk+hyp models with disk filling factor varied. The displayed χ2 in each panel is the mean χ2 of all models and the ranges are given by the standard deviation of the model values. At large radii, the models agree well with observations, but the unresolved central flux is under-represented in the models.
In the text