© The Authors 2023

1. Introduction

Active galactic nuclei (AGNs) are thought to play a crucial role in the formation and evolution of their host galaxies. Moreover, understanding the dust in the vicinity of supermassive black holes is key to understanding how AGNs are fed and how they interact with their hosts. The dust traces dense molecular gas, which feeds the accretion disk (AD). Large, obscuring dusty structures are thought to be responsible for both funneling material toward the AD and for causing the apparent differences between Seyfert 1 and Seyfert 2 galaxies. In the original unified model of AGNs (Antonucci 1993; Urry & Padovani 1995; Netzer 2015), a central obscuring torus of dust is oriented such that the broad-line region (BLR) 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 comprises several key features that vary in temperature from < 100 K to 1500 K and scale from tenths of a parsec to tens of parsecs. The inner edge is the radius at which radiation from the AD causes the dust to sublimate. The sublimation radius is dependent on both the luminosity of the AD and the dust composition, typically ∼0.1 pc for a L ∼ 1 × 10¹⁰ L_⊙ AGN with dust evaporating at 1500 K (see, e.g., Barvainis 1987). Beyond the sublimation zone, it is thought that a dense disk or torus of material is responsible for hiding the BLR in Sy2 AGNs, for feeding the AD, and for reflecting X-rays. 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 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 that it can potentially explain the MIR polar extension (Schartmann et al. 2014). Stalevski et al. (2017, 2019) used a model combining a dusty hollow cone and a thin disk to reproduce the spectral energy distribution (SED) and morphology of the torus in the Circinus Galaxy (hereafter Circinus). A key finding of SED fits to nearby AGNs as well as comparisons to radiative transfer (RT) 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; Schartmann 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 interferometry. 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 VLTI/GRAVITY (GRAVITY Collaboration 2020) and VLTI/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.

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 BLR (Oliva et al. 1998; Ramos Almeida et al. 2016), 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; Kakkad et al. 2023). Additionally, Circinus exhibits a Compton-thick nucleus and a reflection component in X-rays (Matt et al. 1996; Soldi et al. 2005; Yang et al. 2009; Arévalo et al. 2014). Finally, inflows, 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 recently imaged for the first time with MATISSE in the N band (Isbell et al. 2022). These images revealed a dust disk roughly aligned with the water maser emission (Greenhill et al. 2003), as well as warm (∼250 K) large-scale (≳100 mas) emission roughly orthogonal to the disk, similar to previous results with the first-generation MIR interferometer, MIDI (Tristram et al. 2014). 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 the MATISSE images revealed flux enhancements along the position angle (PA) of the optical ionization cone (PA_opt. = −45°). Previous modeling work by Stalevski et al. (2017, 2019) has indicated that this enhanced 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 analyses relied on Gaussian model fitting, which is a smooth, simplified representation of the source emission; interferometric image reconstruction has the potential to build on these results through model-independent sampling of the source structure. In Isbell et al. (2022, hereafter Paper I) we presented the first image reconstructions of Circinus in the N band. In this work we extend the analysis to the LM bands and further consider RT models that represent the source.

This paper is structured as follows. We present the observations used in this work in Sect. 2 along with the necessary data reduction steps. In Sect. 3 we explain our modeling and imaging methods and show the results of each. In Sect. 4 we consider the spatially resolved temperatures of the recovered images and models. Finally, we discuss the implications of this work in Sect. 5 and conclude 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 February 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 UT configuration, with physical baselines ranging from 30 to 140 m. The observations are described in more detail in Paper I. 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, and circles represent observations from 2021. The mean correlated flux between 4.5 and 4.9 μm is used as the color scale. North is up, and east is to the left.

On each night, we observed the calibration star HD120404 (F_12 μm = 13 Jy) directly before and/or after the Circinus observations. This star serves as a spectral calibrator, as an instrumental phase calibrator, and as 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 February and May 2021 observations, we observed secondary calibrators, HD120913 (F_12 μm = 5.7 Jy) and HD119164 (F_12 μm = 1.2 Jy) in order to perform cross-calibration and closure phases accuracy checks. All Circinus and calibrator observations entering this analysis are given in Table 1 of Paper I.

2.2. MATISSE LM-band data reduction and calibration

The LM-band data for Circinus and the calibrators were reduced using both the official data reduction software (DRS) version 1.5.1 and custom scripts. We reduced the data both coherently and incoherently. For the coherent reduction we used the flags corrFlux=TRUE and coherentAlgo=2 in order to produce correlated fluxes. For the incoherent reduction we used the flag corrFlux=FALSE. In both cases we used spectral binning 11 px (= 0.5 μm) and the default values for all other parameters. The DRS is not optimized for coherent reduction of the LM bands, but coherent integration is necessary for faint sources (see, e.g., the flux limits given in Lopez et al. 2022). The above settings (both coherent and incoherent) resulted in strange spectra with (1) no M-band flux and (2) sharp emission features at 3.7 μm; but analysis of the intermediate products (specifically, the cleaned interferogram) found neither of these features. It was found that a bias floor was present in the DRS-reduced data. This is likely due to the fringe search being optimized for 3.6 μm, but the Circinus spectrum is very “red” and the L-band flux is very low. Instead, a fringe extraction using the M band was necessary.

The data are then re-reduced using a custom python script¹. The custom pipeline uses the intermediate products of the DRS, specifically the complex cleaned interferograms (files called OBJ_CORR_FLUX). Using the 4.6 μm flux, the six fringes are identified and extracted in each exposure snapshot and each beam commuting device (BCD) configuration. Additionally, a bias “fringe” per frame is extracted far from the science fringes. The extracted fringes for each baseline and the extracted bispectra for each closure triplet are then bias corrected and temporally averaged incoherently (over the exposure cycle). More details are provided in Appendix A.

The resulting correlated fluxes and closure phases are computed for each BCD independently, and the final values are taken as the mean of the four BCD configurations. The final errors are the standard deviations of the four BCD configurations. In the L band the correlated flux errors on individual baselines are typically 1.8 mJy, and in the M band the correlated flux errors are typically 5–10 mJy, in comparison to the typical L-band flux of 2–10 mJy and the typical M-band flux of 10–40 mJy. In both bands the closure phase errors are in general quite large, ≳90°, and their use is limited. This process was done for both the calibrators and Circinus, and the resulting observables were calibrated as usual. Reductions using the DRS and the custom pipeline for both a calibration star and for Circinus are shown in Fig. A.1 for comparison and validation of the approach.

The correlated flux, F(u, v, λ) is then calibrated in the same way as in Paper I. The squared visibilities are finally calculated as , where the total flux comes from the shortest baseline correlated fluxes; in this case we used the azimuthal maximum of the 30–35 m correlated fluxes as an estimate. The 30 m total flux and sample of all the correlated fluxes are shown in Fig. 2 (the remaining correlated fluxes are shown in Fig. B.1). We use this 30 m correlated flux rather than the “zero-baseline” correlated flux because squared visibilities can be scaled somewhat arbitrarily; the relative changes give substructure and the absolute changes are indicative of over-resolved emission. The 8.2 m total flux from ISAAC observations (F_nuc, L = 458.16 ± 39.18 mJy, F_nuc, M = 676.41 ± 44.58 mJy; Isbell et al. 2021) would result in extremely small squared visibilities (∼3 × 10⁻⁴) and cause numerical issues with little gain in understanding of the source.

Fig. 2. Excerpt sample of calibrated correlated fluxes for Circinus in the LM band. The top-left panel shows the 33 m “total flux” used to compute squared visibilities. The colors indicate the observing date: March 2020 (blue), February 2021 (yellow), and May 2021 (red). Black squares and black circles are the predicted correlated fluxes from the image reconstructions and Gaussian model fits, respectively. Continued in Fig. B.1.

The closure phases are calibrated as in Paper I, but they instead use the results of the custom pipeline. Due to the low signal-to-noise ratio (S/N) for the visibilities on many baselines, only four closure phase triangles give reasonable values. The others are all dominated by noise, as mentioned above. The closure phases are shown in Fig. B.2.

2.3. Single-dish observations

This section describes single-dish observations of Circinus in the L and M bands from NACO and VISIR at the Very Large Telescope (VLT). These observations provide large-scale context for the interferometric images, linking the subparsec scales to the galactic scale. Additionally, we looked for signatures of flux variability.

2.3.1. NACO observations

We observed the Circinus nucleus with VLT/NACO in the Lp and Mp filters (3.4–4.2 μm and 4.4–5.2 μm, respectively) in burst mode in 2018 as part of program 0101.B-0446. A point spread function (PSF) and flux calibrator, IRAS 14480-5828, was observed in concatenation. It is the only target within 10° that has roughly similar color and brightness in the near-infrared and is similar or brighter at L and M (W1 ∼ 6 mag, W2 ∼ 4.5 mag). It is a basically unknown object, but is classified as a star in the SIMBAD database. It has no apparent optical counterpart.

We fit a two-dimensional Gaussian to each PSF-limited observation. This yields the full width at half maximum (FWHM) of the major and minor axes as well as the PA of the major axis for the target and for the calibrators². The calibrator IRAS 14480-5828 seems unresolved in Mp with FWHM 0.13″ × 0.13″, PA = 28°. However, it turned out to be extended in Lp with 0.20″ × 0.19″, PA = 112°). From SED fitting of the Spitzer/IRAC and WISE data, we obtain the following fluxes for IRAS 14480-5828 in Lp and Mp of 1501.6 mJy and 2313.0 mJy, respectively. Because IRAS 14480-5828 is resolved in Lp, we used an alternative calibrator, 2MASS J16232835, with an Lp flux of 89.43 mJy. This star was observed as part of 0101.C-0924 (PI: Christiaens; unpublished) in the same night as our Lp observation. The star is unresolved with 0.12″ × 0.11″ and PA = 167°. The nucleus of Circinus appears to be extended with 0.17″ × 0.14″, PA = 114°, and 0.19″ × 0.14″, PA = 108° in Lp and Mp, respectively. The achieved angular resolution and depth is significantly better than in the archival NACO and ISAAC observations.

The calibrated Lp and Mp fluxes for the Circinus nucleus are listed in Table 1. Depending on the choice of the flux reference we obtain fluxes that differ by a factor of ∼2. The likely reason for this is that at least one of our flux references is variable. We see some evidence for this in the SED of IRAS 14480-5828 where there are factor of 1.5 differences between the IRAC and WISE measurements. Flux calibrated images are shown in Fig. C.2 and fluxes extracted in 0.4″and 4.0″apertures are reported in Table 1. We adopted the values obtained by calibrating the L band and M band with 2MASS J16232835 and IRAS 14480-5828, respectively.

The 0.4″Lp filter flux measurement presented herein is ∼50 mJy larger than the L-band unresolved flux measured with ISAAC reported by Isbell et al. (2021). This difference is likely due to the use of aperture extraction in this work, while those authors did PSF extraction via Gaussian fitting of multiple components; their method attempted to separate resolved and unresolved flux.

2.3.2. VISIR M-band observations

We observed the nucleus of Circinus with VLT/VISIR in the M-band filter in burst mode in 2017 as part of program 099.B-0235 (PI: Hoenig). During the observations, the VISIR M-band filter suffered from a red leak, but we were able to roughly correct for this by subtracting a scaled PAH2_2 filter image. Nevertheless, all the flux values reported here should be regarded as upper limits, including the typical 10% systematic uncertainty on the flux of the calibrator.

The nucleus of Circinus appears extended in both direct and PSF-subtracted images. The direct image is shown in Fig. C.3 and has a size of 0.33″ × 0.30″ with PA = 98°, while the calibrator, HD138538, has a size of 0.22″ × 0.18″ with PA = 108°, and was observed directly before Circinus in a similar direction of the sky. We measured the flux of Circinus using apertures of two different sizes: 0.4″ and 4.0″. The resulting fluxes are given in Table 2.

2.3.3. Checking for variability

Because of the 0.4″aperture observations presented in this work, we can check the stability of the LM-band fluxes in Circinus over 18 years. We show the collected flux measurements in Fig. 3, bringing together data from Prieto et al. (2004), Stalevski et al. (2017), Isbell et al. (2021), and this work. We find that Circinus’s LM flux has remained fairly constant over the past decades. However, the L-band measurement from Stalevski et al. (2017) is ∼2× higher than the other three L-band fluxes; there are no apparent issues with the calibrator, HD129858 according to the MIR calibrator catalog of Cruzalèbes et al. (2019). However, Stalevski et al. (2017) used the MIDI observations of Tristram et al. (2014) to measure in a 0.4″ aperture, while Isbell et al. (2021) uses single-dish observations in which the PSF was distorted or significantly larger than 0.4″. The apparent difference could be due to this and/or the different method of measurement (PSF-fitting vs. aperture extraction).

Fig. 3. Circinus 0.4″ fluxes in the L and M bands between 2000 and 2018. We compile the measurements of Prieto et al. (2004), Stalevski et al. (2017), Isbell et al. (2021), and this work. L-band values are plotted in blue, and M-band values are plotted in red. The mean flux for each band is plotted as a dashed line, with the standard deviation of the fluxes given as a shaded region.

3. Modeling and image reconstruction

Both image reconstruction and model fitting rely on minimization of a cost function, q, which measures the flux distribution’s similarity to a weighted combination of the closure phases and squared visibilities. As in Paper I and Hofmann et al. (2016, 2022), we use a modified χ² function:

(1)

In the case of image reconstruction, a regularization term is added, yielding the cost function for an image, I, sampled at the vector of uv coordinates, x,

(2)

where q represents Eq. (1), and ϕ_model are the squared visibilities and closure phases of the image flux distribution sampled at x, μ is the so-called “hyperparameter” that sets the amount of regularization applied, and R the regularization function to be applied.

We first attempted model-independent image reconstruction, as these are the first interferometric observations of Circinus in the LM bands, and the source flux distribution was unknown. Following the relatively simple results of the image reconstruction, we performed follow-up modeling of the flux using elongated Gaussian components.

3.1. Image reconstruction

Following the procedure in Paper I, we use the package included in the DRS, Image Reconstruction Software using the Bispectrum (IRBis; Hofmann et al. 2014, 2022) to reconstruct images in the LM band. We selected two wavelength bins in which to produce independent images: 3.7 ± 0.1 μm and 4.7 ± 0.1 μm. Any spectral information within each bin is averaged, producing a series of “gray” images. Each bin is imaged with a range of regularization functions and hyperparameters (hereafter μ; essentially a scaling on the amount of regularization), with the best selected via Eq. (1). We performed a grid search of the IRBis parameters, varying the field of view, the pixel number, the object mask scale, the regularization function, and the hyperparameter μ. We used nonuniform weighting in the uv plane, setting weighting=0.5 in IRBis to de-emphasize the sparsely sampled and low-S/N points on baselines longer than ∼60 m.

An initial best image is selected in each wavelength bin using minimization of Eq. (1), and a follow-up round of imaging using the best regularization function and pixel scale is performed. We give the final parameters for the reconstructions in Table 3. The images are shown in Fig. 4 and the predicted correlated fluxes and closure phases from the reconstruction are overplotted in Figs. 2 and B.1–B.2; here it is apparent that the sparse uv coverage and low S/N of the correlated fluxes have resulted in significant image artifacts, particularly in the M band. We also see that the highest correlated flux values are not always matched in the M band, usually when low correlated flux values are found at adjacent uv points. Image errors are estimated as in Paper I, using delete-d jackknifing.

Fig. 4. Images and models of the Circinus LM-band dust. Top: Image reconstructions at 3.7 and 4.7 μm. The contours are at 10 × σIm, estimated from the error maps produced by delete-d jackknifing the uv coverage. The FWHM of the estimated dirty beam is given in the bottom-right corner of each panel. Also included is the N-band continuum image from Paper I for reference. Bottom: Gaussian model fits to the same wavelengths, specifically the 1+ model in the L band and the two-Gaussian model in the M band. The dashed ellipse in the L-band model image represents the 1σ PA uncertainty of the fitted model. The colored circles in the bottom-right panel illustrate the extraction apertures for the measured fluxes.

Despite the artifacts, the results of (1) a point source in the L band and (2) a highly elongated component in the M band at ∼45° PA are robust above the noise (i.e., > 10σ, where σ is estimated on a pixel-by-pixel basis using the delete-d jackknifing method described in Paper I). The point-like source is present in the L band, but it is not present in the M-band image reconstruction. There is, nonetheless, a peak in flux in the M-band image at the same pixel position as the point source. The disk-like component is roughly 4 mas in FWHM. This is comparable to the ideal resolution at 4.7 μm (3.7 mas), so the disk width is only marginally resolved. The structure extends 20.2 mas in FWHM along PA≈46°. There are several secondary features that cross the M-band image along PA ∼40; these are certainly artifacts due to their regular spacing, symmetry, and flux falloff with radius. Peaks along these lines indicate potential real polar flux, but they are at too low significance to analyze robustly.

3.2. Gaussian modeling

The image reconstructions result in simple components: a point source at 3.7 μm and a disk-like structure at 4.7 μm. These simple structures can be modeled using elongated Gaussian components, which serves to minimize the effect of the dirty beam and image artifacts due to sparse uv coverage and a small number of closure phases.

The fitted Gaussian components are fixed to the center of the image and only their major axis Θ, their minor-to-major axis ratio r ≡ θ/Θ, the PA of the major axis ϕ, and their relative flux f are allowed to vary. We fit a number of Gaussian components N_comp to each wavelength channel, but in general favor models with fewer components. The relative flux f of one Gaussian is fixed to 1. The fitted Gaussians’ parameters and error estimates are obtained through Markov-chain Monte Carlo likelihood maximization. We sampled the parameter space using the package emcee (Foreman-Mackey et al. 2013). The log-probability function to be maximized is given by the typical Bayesian formulation,

(3)

with measurements at uv coordinates x, parameters θ, and error estimates σ scaled by some constant c. For maximum likelihood estimation, the log likelihood function for the models is written as

(4)

where q (Eq. (1)) is the cost function for the squared visibilities and closure phases produced by the model with parameters θ, c represents the underestimation of the variance by some fractional amount, and . We estimate the best-fit value as the median of each marginalized posterior distribution and the 1σ errors from the values at 16th and 84^th percentiles. Because the closure phases are very low S/N, we fix all components to the center of the image, and we set β = 0 in Eq. (1) to fit the squared visibilities alone.

We first fit one Gaussian (i.e., N_comp = 1) based on the simplicity of the image reconstruction, but Gaussian fits with more components were attempted as well. We favor models with fewer parameters based on the Akaike information criterion (AIC; Akaike 1981). The AIC for a model with k parameters and maximum likelihood L is

(5)

The model with the minimum value of a is considered the “best” because it is a sufficient representation of the data without overfitting. For each fitted model with Gaussians fixed at the center, we fit up to k = 4N_comp − 1 parameters (because the flux f of one component is fixed). The fit results for 3.7 μm with N_comp = 1 and 4.7 μm with N_comp ∈ {1, 2, 3} are given in Table 4 and the best fitting models are shown in Fig. 4.

In the L band, a true single-Gaussian model does not perform well. It becomes large in order to produce low visibilities at short baselines, but then the long-baseline visibilities are far too low. An augmented single-Gaussian model was then fit, wherein a second, large component was added. This second component has fixed size (128 × 128 mas) and orientation, but its flux is allowed to vary. It plays the role of over-resolved flux. This model (called 1+) gives a marked improvement in AIC over the single-Gaussian model with only one additional parameter. It results in a marginally extended source with FWHM 7.1 × 2.8 mas. A two-component Gaussian model produces a similarly extended source (6.5 × 2.0 mas) with a much less extended second component. The two-component model produces the same AIC value as the 1+ model; for the rest of this work, we used the simpler 1+ model.

In the M band, all three modeling results include a disk-like component with a similar size and orientation; this disk-like component reproduces the image reconstruction’s morphology. In the two- and three-component fits, a point-like source is introduced on top of the disk. In the three-component model, a diffuse, roughly polar extended source is added. While its orientation is suggestive (based on the N-band polar emission), this additional component is disfavored by the AIC. If there is a signature of the polar dust in the M band, it is at low significance, and more observations would be necessary to confirm it. The model with two Gaussian components is preferred at 4.7 μm, and it will be used for the rest of this work. We note, however, that the flux in the central aperture (see the following section) differs by only 1% between the one- and two-component models. The selected model is marginally resolved with a width of 4.2 mas (the 4.7 μm resolution is 3.7 mas) and has a major axis with mas and PA = . The increased size of the disk in the modeling when compared to the images likely compensates for the large-scale flux that is allowed in the imaging. This scenario is supported by the observation that the PA of the disk component decreases with an increasing number of fitted components.

The M-band model and image each show maximum extension nearly perpendicular to the beam rather than in the polar direction. This is not unprecedented when comparing to the N-band images of Paper I (and to some extent to the L-band image or model). In the N band, while the maximum extension is in the polar direction, a significant amount of flux is perpendicular to the beam in a disk-like component. Similarly, the L-band model has uncertainties in PA that mean it could also be aligned with the M-band extension, and the L-band image shows secondary peaks along that same direction. We would argue that indeed all images and models in all bands either show or allow for flux in this direction. The apparent lack of polar dust in the M band is most likely due to a combination of low dynamic range in image reconstructions and the fact that large-scale structures (> 16 mas = 0.3 pc) are resolved out by the interferometer.

4. Measuring component temperatures

In Paper I, circular apertures with diameter 23.4 mas were used to extract the flux from each N-band image at a number of locations. We made the assumption here that the LM models can be astrometrically matched with the N-band images by aligning the photocenters. This was done for NGC 1068 in Gámez Rosas et al. (2022), and their cross-correlation matching done in the N band was in the end equivalent to photocenter matching. Therefore, we used the same aperture diameter and distribution as Paper I.

In the L band, although there is a point source within the central aperture, D0, much of the flux is contained in the background component. This is true in both the imaging and in the 1+ model. Accordingly, only 86% (in the Gaussian model) and 37% (in the image) of the total L-band flux is contained in the aperture D0. In the M band, the majority of the flux is found in the central aperture with only a minuscule amount falling in the disk apertures D_S40 and D_N40. In the imaging and modeling of neither band is there significantly measured flux in the polar direction. For these regions we present an upper limit from the “sky” background F_upperlim. ≤ 2σ_sky in the image reconstructions. All apertures from the N-band analysis that are not listed in Table 5 are considered to have only upper limits: F_3.7 μm ≤ 1.2 mJy and F_4.7 μm ≤ 2.2 mJy.

The Gaussian models and image reconstructions give slightly different morphologies, and therefore yield different flux measurements. We considered both sets of results independently. The 3.7 and 4.7 μm fluxes for the apertures from each imaging method are given in Table 5. The uncertainties in each case come from either the image or model uncertainty (each is described above) and the uncertainty on the total correlated flux (at 30 m). The total flux uncertainty (∼5 mJy at 4.7 μm) dominates in both imaging and modeling.

4.1. Blackbody fitting

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

(6)

where η is an absolute flux scaling due to the filling factor of the dust in the aperture, τ(λ)/τ_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 to the LMN SEDs is done in two iterations using Markov chain Monte Carlo sampling with the package emcee (Foreman-Mackey et al. 2013), similar to the approach in Sect. 3.2. 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.

For the first BB component, we use uniform prior probability distributions with T₁ ∈ (100, 500] K, η₁ = 1, and A_V, 1 ∈ [20, 37] mag. These priors are based on the N-band fit results, particularly mag for D0 (see Paper I). The second component is forced to be strictly hotter and smaller than the first component, resulting in the uniform priors T₂ ∈ (500, 1500] K, η₂ ∈ (10⁻³, 0.1], and A_V, 2 ∈ (0, 700] mag. For the central aperture, D0, which should cover the sublimation zone, another set of priors is also used. They come from the assumption that dust is the sublimation temperature is indeed present but can be heavily obscured. We estimated a representative sublimation radius for silicate dust at 1500 K using the formula from Barvainis (1987):

(7)

As Moorwood et al. (1996) report L_UV = 5 × 10⁹ L_⊙, we estimate r_sub = 0.05 pc = 2.8 mas at a distance of 4.2 Mpc. This gives an upper limit on η₂ ≤ (2.8/11.7)² = 0.06, where 11.7 mas is the aperture radius. Therefore, we define the priors of the second component: T₂ = 1500 K, η₂ ≤ 0.06, and A_V, 2 ∈ (0, 700] mag. This fit gives a rough estimate on the minimum amount of extinction necessary to hide dust at the sublimation temperature. It is important, however, to note that for different dust compositions and grain sizes, the sublimation radius can vary from 0.05 to 0.2 pc. This is compounded by uncertainties in the luminosity of the AGN. The value we adopted is meant only as a fiducial, representative value that is roughly the mean radius of the sublimation zone.

The recovered temperatures for the one- and two-component BB fits to the image reconstruction fluxes and to the Gaussian model fluxes are given in Table 6. The fitted SEDs for the Gaussian models using one and two BB fits are shown in Figs. 5 and C.1, respectively. The nuclear K-band flux from Burtscher et al. (2015) for Circinus is shown in Figs. 5 and C.1, and it serves as an upper limit on the near-infrared flux; all extrapolated BB flux values are far below this limit.

Fig. 5. One-BB fits for the aperture-extracted Circinus LMN fluxes. The colors are the same as in Fig. 4, with D0 in red, DN40 in yellow, and DS40 in red. The top-right plot compares the D0 fluxes to the various RT model modifications. The fits using the N-band data alone are included for comparison. In aperture D0 the K-band measurement from Burtscher et al. (2015) is included as an upper limit for the near-infrared flux.

For the aperture D0, the image reconstructions’ fluxes result in lower fitted temperatures than the Gaussians’. However, the fitted T₁ values in both the images and the Gaussian models are consistent within the uncertainties to the temperature inferred from the N band alone ( K). The fitted “cool” component temperature T₁ has essentially the same value both with and without the additional hot BB. Importantly, in neither the images nor the models is an additional, hot component necessary. Large extinction values (∼450 mag) are preferred, and the hot component makes up only 0.2% of the 3.7 μm flux and 1.3% of the 4.7 μm flux. In fact, looking at the posterior probably distributions for all apertures, only second component extinction values ≳250 mag are allowed in any of the fits.

For the disk apertures, D_N40 and D_S40, the BB-fit results are consistent from both the modeling and image reconstruction fluxes. The results, furthermore, are consistent with the N-band fitted temperatures and extinctions. A second component is once again disfavored due to the large fitted extinction values and the fact that the fit is not markedly improved with the addition of this component.

4.2. Comparisons to radiative transfer models

To describe the circumnuclear dust, a disk+wind model has become the consensus in both RT modeling (e.g., Stalevski et al. 2017, 2019; Hönig & Kishimoto 2017) and hydrodynamical modeling (Williamson et al. 2020). Specifically, Stalevski et al. (2017, 2019) undertook RT 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 RT 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. This model was shown in Paper I to broadly agree with the N-band MATISSE data, though new constraints on the clump number density were reported.

In this section we explore the applicability of those RT models and their SEDs to the combined LMN-band results. We used the same procedure for the RT model images as for the observations: we extracted fluxes in each of the 23.4 mas apertures and fit temperatures and extinctions using Eq. (6) to the disk+hyp model grid at λ ∈ [3.7, 4.7, 8.5, 8.9, 9.3, 9.7, 10.1, 10.6, 11.0, 11.5, 12.0, 12.5] μm. The total flux of each RT model image is scaled such that the central 500 × 500 mas flux matches the total MATISSE-observed flux at that wavelength.

The most striking difference between the RT models and the observations is found in the central aperture. Both the extracted flux and the fitted temperature are found to be low in the disk+hyp models compared to the observations ( K vs. 290 K). We modified the disk+hyp models in several different ways to determine where additional warm dust could plausibly be found:

1. The inner radius of the fiducial disk+hyp model is 0.16 pc. We decreased this to 0.1, which increases the flux impinging the inner hyperboloid walls. This inner radius change produced only a small difference in fitted temperature (∼4 K)

2. Diffuse dust with mass fractions of (0.7 to 0.97) was distributed between the clouds of the fiducial disk+hyp model. Diffuse dust actually decreased the fitted temperature in some cases, likely due to more self-shielding of the dust structures.

3. Finally, we added a cluster of clouds positioned just above the disk and inside the hyperboloid at a distance of 0.15 pc (∼3r_sub). This distance was set by the temperature versus radius estimates of Eq. (7); much closer in and the clouds would both be too hot (we see no evidence of T ∼ 1000 K dust) and/or blocked by the disk. In a standalone model, they were able to reach temperatures ∼330 K. When added to the existing disk+hyp model, this effect is reduced, but changing the optical depth or clump density of the cone wall would make it possible to glimpse these clouds more directly.

In each case, we measure the temperature in the central aperture, D0. We show the predicted temperatures and SEDs in Fig. 5. We find that neither decreasing the inner radius nor adding diffuse dust made significant changes to the measured temperature. On the other hand, adding clouds above the disk but inside the hyperboloid greatly increased the measured temperatures, bringing the models into much better agreement with the observations. We show a schematic of the resulting RT model in Fig. 6. We did not constrain the amount or exact location of these clouds because they are found on scales that remain spatially unresolved in our data. We instead emphasize that this simple modification of the nominal disk+hyp model at the base of the cone causes a large qualitative improvement in the agreement between the (modified) disk+hyp models and the observations.

Fig. 6. Schematic of the dust within the RT models that best reproduce the observed dust temperatures. The basis of the model comes from Stalevski et al. (2019), and a rough parameter range is given in Paper I. Here we have added a cluster of dust clouds above the disk and behind the hyperbolic cone, at a position that in projection corresponds to the central aperture. These clouds could represent, for example, a puffed-up sublimation zone, freshly launched winds, or a smoother boundary for the hyperbolic cone. The spatial resolution of these observations is not sufficient to distinguish between these scenarios.

5. Discussion

Both the image reconstructions and the Gaussian models of the LM-band data are remarkably consistent with the structures imaged in the N band. Specifically, in the M band we recover a thin disk at PA = . This disk is mas FWHM across – comparable to the ideal 3.73 mas MATISSE resolution at 4.7 μm (note that the actual resolution is likely lower due to sparse and elongated uv coverage). The M-band disk width is thus at best marginally resolved, and it is therefore consistent with the edge-on measurements of Paper I and the general trend of increasing inclination toward the center of Circinus (Izumi et al. 2018). Additionally, in the M-band image reconstruction (Fig. 4), we see relatively bright artifacts in the polar direction. This suggests the presence of over-resolved flux at these scales, which we also see in single-dish ISAAC and NACO images at ∼400 mas resolution (this work and Isbell et al. 2021).

Taken together, the single-dish and interferometric data suggest that the observed LM emission relates more to large-scale dust structures than to the scales or dust temperatures observed with dust-reverberation mapping (e.g., Koshida et al. 2014). We see little evidence of temporal changes in the single-dish flux over the last decades (Fig. 3). This lack of variability indicates that the emission originates not in directly illuminated, fast-moving clouds near the AD, but rather from a more distributed, large-scale dust structure. This conclusion is supported by the interferometric images and SED fits. We show that the LMN emission can be simply explained by a single moderately warm dust structure missing any direct sightlines to very hot dust clouds. Specifically, aperture D0 is dominated by ∼367 K dust emission, and any hotter dust emission is shown to be hidden by A_V ≳ 250 mag of extinction. Furthermore, the morphological continuity between the LMN bands suggests that the emission traces the same structure from ∼0.1 pc to > 1 pc scale.

It is tempting to conclude that the M-band image reveals a warp in the central disk. There is much evidence of a tilted (with respect to the ionization cone) AD in Circinus, for example the warped H₂O maser emission (Greenhill et al. 2003), fits of RT models to MIDI data (Stalevski et al. 2019), and recent polarization mapping data (Stalevski et al. 2023). In short, it seems that the warped maser emission does indeed trace a warped or tilted AD, which asymmetrically illuminates the circumnuclear dust. This was a possible explanation of the flux- and temperature-enhanced regions of the polar dust imaged in the N band in Paper I. Radiation transfer modeling of dust around a warped maser disk by Jud et al. (2017) showed that the warp would indeed be visible in the warm, innermost circumnuclear dust. Despite the evidence of a warp from other observations, we caution that the M-band closure phases used to reconstruct our image have very large (> 100°) errors in most closure triangles. So, while the correlated fluxes strongly constrain the spatial scale of the emitting M-band flux, the specific morphology relies on a handful of significant measurements. Either repeated observations under excellent conditions or improvements in data processing would be necessary to confirm the presence of a true warp in the M-band-emitting dust.

Both the image reconstructions and Gaussian models give very similar SED and fitted temperature results. We extracted the flux in three apertures with diameters of 23.4 mas. These are a subset of the 13 apertures studied in Paper I, allowing us to study the spatially resolved LMN-band SEDs. In all apertures, a single BB temperature is sufficient to describe the measured fluxes. There is no compelling indication of the expected hot, approximately sublimation temperature dust in the central aperture. This begs the question whether it is simply not there, or whether it is present and being obscured by the disk component. If the hot dust is simply not there, this raises serious questions and is difficult to reconcile with the latest radiation transfer and hydrodynamical modeling (e.g., Wada et al. 2016). When we assume the dust is there, we find that very large extinction values (≳250 mag) are necessary to match the observed flux in the central aperture. This value, however, is not incredible. Using the X-ray hydrogen column density to extinction relation of Predehl & Schmitt (1995),

(8)

and the value N_H = 10²⁴ cm⁻² from Matt et al. (1999), we find A_V = 558.6 mag. Recently, ALMA observations by Izumi et al. (2023) indicate extinction toward the AD of A_V = 210 − 440 mag. These estimates are fully in line with our fitted value. Our fitted value is also in broad agreement with predictions from the disk+hyp RT model; the A_V of the model at the given inclination, averaged over the azimuthal angle, is approx A_V = 211 mag (Stalevski et al. 2019); however, this might be a lower limit, since the model grid did not include optical depth values larger than τ_9.7 = 15 for the disk. All in all, it seems that the optically thick, geometrically thin disk plays the role of the dusty obscurer in the classic dusty torus model.

We test several ad hoc modifications to the disk+hyp models of Stalevski et al. (2017, 2019) in order to better reproduce the observed fluxes and temperatures of the central aperture, D0. We find that a cluster of clouds positioned above the disk and inside of the hyperboloid at a distance of 0.15 pc (∼3r_sub) matches the observations most closely. While we did not constrain the exact scale or origin of these clouds, they indicate a thicker hyperboloid wall in this region and/or a puffed-up disk near the center. In both cases, they indicate (dynamic) dust structures in this region strongly affected by the asymmetric radiation pressure of the tilted AD (see, e.g., Paper I Stalevski et al. 2019, 2023, for evidence of this tilt). A puffed-up inner region has been previously proposed by Hönig et al. (2012). Within a zone of 5–10 sublimation radii (0.25–0.5 pc or 12–24 mas), the predicted puffed-up region described in Hönig (2019; based on 3–5 μm SEDs) should be resolved with our MATISSE observations. As we imaged this inner region for the first time, we placed constraints on the scale of this puffed-up zone in Circinus, roughly 3 r_sub. If the puffed-up region indeed extends farther (to ∼10 r_sub), it could contribute to the A_V ∼ 500 mag obscuration of the “missing hot dust” mentioned above. Future high-resolution MIR observations of other Seyfert galaxies can explore how typical the circumnuclear dust in Circinus is. Additionally, the dust in this regime could be related to the radiation-driven fountain suggested in Wada (2012). Whether the new clumps represent failed flows or freshly launched dusty outflows is unclear. The latter would represent the onset of feedback with the host galaxy. Future, very high spatial resolution infrared spectra in this region could give an idea of the dust dynamics in Circinus.

6. Conclusion

Following up on the N-band analysis of Paper I, we present the first-ever L- and M-band interferometric observations of Circinus. These observations allowed us to reconstruct images and fit Gaussian models to the L- and M-band data. Using these images and models:

1. We find a thin disk whose width is marginally resolved (0.08 pc = 4.23 mas). This disk is shown to be the spectral continuation of the disk imaged in the N band to shorter wavelengths, as the measured fluxes correspond to the fitted N-band temperatures. In addition to this thin disk, there is a point-like source found in the L and M bands that was identified with the N-band point source based on the LMN-band SED fit.

2. We show that there is no trace of hot dust (T ∼ 1500 K) in the circumnuclear dust structure of Circinus. By assuming the dust is there, we find that obscuration within the disk of A_V ≳ 250 mag is necessary to reproduce the measured fluxes. With dust extinction this high, the imaged disk could then play the role of the obscuring “torus” in the unified scheme of AGNs.

We also explore the parameter space of the disk+hyp models from Stalevski et al. (2019), identifying a simple modification that better matches the observations toward the center. We added a cluster of dust clouds above the disk but inside the radius of the base of the hyperbolic cone, at a position that in projection corresponds to the central aperture. These clouds could represent, for example, a puffed-up sublimation zone, freshly launched winds, or a smoother, radially wider boundary for the base of the hyperbolic cone. The spatial resolution of the presented observations is not sufficient to distinguish between these scenarios, but the presence of these clouds provides a much better match to the observed temperature distribution of the circumnuclear dust in Circinus than the standard disk+hyp models. In the future, we plan to test the applicability of this inner dust structure to other MATISSE-observed AGNs.

Acknowledgments

The authors would like to thank the late Prof. Dr. Klaus Meisenheimer for his valuable contributions to both this work and to their growth as scientists, through his inspiring supervision, thoughtful questions, and ambitious vision for the VLTI. We thank the VLTI Paranal staff for all their help and good company during observations. M.S. is supported by the Ministry of Science, Technological Development and Innovations of the Republic of Serbia through the contract No. 451-03-9/2023-14/200002. 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. The authors thank the anonymous referee for their helpful comments which greatly improved this paper.

References

Appendix A: Details on custom data reduction

The MATISSE DRS has been optimized for sources with an L-band flux F_3.6 > 0.5 Jy or an N-band flux F₁₂ ≳ 1 Jy. Moreover, it uses 3.6 μm as the reference wavelength for the estimation of the optical path delay (OPD) in the L band. This works well for bright sources or for blue sources (such as young stellar objects). For red objects, such as Seyfert 2 galaxies, OPD corrections in the L band are severely impacted by low S/N and so the M band should be used. At the time of writing, there is no option to change the OPD reference wavelength. For the data reduction of Circinus, it was necessary to manually extract the fringes using the M band as a reference.

Fortunately, the Charge-Coupled Device (CCD) calibration steps and the production of the clean interferogram could be taken directly from the DRS. These steps resulted in a time series of two-dimensional Fourier transforms of the interferogram (called OBJ_CORR_FLUX_*.fits). In each frame, the positions of the fringes could be measured and shifted (in other words, phase-corrected) before the correlated flux and closure phase could be extracted using the following methods; the fundamental equations we used can be found in the MATISSE instrument paper, Lopez et al. (2022).

A.1. Squared (incoherent) correlated flux

The squared correlated flux³ for a baseline B_ij is given as

(A.1)

where ⟨...⟩_t is the time average, u is the spatial frequency integrated between (B_ij − D)/λ and (B_ij + D)/λ, D is the pupil diameter, and β is the bias present in each fringe peak (estimated from the value between fringe peaks). The pupil here refers to the spatial-filter pinhole with diameter 1.5λ/d in the L band and 2λ/d in the N band with d representing the telescope diameter.

A.2. (Coherent) correlated flux

The correlated flux is integrated coherently in the form

(A.2)

where once again u is the spatial frequency integrated between (B_ij − D)/λ and (B_ij + D)/λ and D is the pupil diameter.

A.3. Closure phase

Recalling that interferometric observations are complex numbers that have both an amplitude and a phase (the argument of the complex number), one can compute the closure phase for a telescope triplet ijk, ϕ_ijk in the following way:

(A.3)

The coordinates u₁, u₂, and u₁ + u₂ correspond to the baselines B_ij, B_jk, and B_ik, respectively, integrated over the pupil width as for the squared correlated flux. The parameter γ represents the photon bias in the bispectrum. This photon bias contains both an additive and a multiplicative term due to the combination of the fringes. The estimation of this bias is nontrivial, and the method used in the DRS is given in, for example, Gordon & Buscher (2012). However, in the limit where read-noise is negligible and photon noise dominates, γ ≈ |β_ij|² + |β_jk|² + |β_ik|² − 2N where β_ij is the photon bias used in Eq. A.1 for a baseline B_ij and N is the mean number of photons in the interferogram (Wirnitzer 1985). For bright sources, the γ term is negligible.

Fig. A.1. Comparison of DRS-extracted and custom-extracted correlated fluxes for Circinus (left) and a calibration star, HD120404 (right). From top to bottom, the panels show the cleaned fringe pattern (the Fourier transform of the cleaned interferogram), the six fringes extracted with the DRS, and the same six fringes extracted with our method. In the Circinus fringe pattern, there is significant flux in the M band that is missing in the DRS fringes but is present in our method.

Appendix B: Reduced data

B.1. Correlated fluxes

In Figs. 2 and B.1 we present the LM-band correlated flux for each baseline, reduced and calibrated as described in Sect. 2.1.

Fig. B.1. Continuation of Fig. 2.

Fig. B.1. continued.

Fig. B.1. continued.

B.2. MATISSE LM-band closure phases

In Fig. B.2 we present the LM-band closure phase spectrum for each closure triangle, reduced and calibrated as described in Sect. 2.1.

Fig. B.2. LM-band 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.

Fig. B.2. continued.

Fig. B.2. continued.

Fig. B.2. continued.

Appendix C: Two-BB fits

In Fig. C.1 we show the resulting two-BB fits for the Gaussian models discussed in Sect. 4.

Fig. C.1. Two-BB fits for the aperture-extracted Circinus LMN fluxes. The colors are the same as in Fig. 4, with D0 in red, DN40 in yellow, and DS40 in red. There is little discernible improvement in fit quality with the addition of a second component, and the second component is in all cases highly extincted. The fits using the N-band data alone are included for comparison. In aperture D0, the K-band measurement from Burtscher et al. (2015) is included as an upper limit for the near-infrared flux.

Fig. C.2. NACO flux-calibrated images for IRAS 14480-5828 and Circinus.

Fig. C.3. VISIR flux-calibrated M-band image of Circinus with 0.4", 1", 2", and 4" apertures overlaid. Here we clearly see the extended polar dust structures to the east and west of the nucleus.

All Tables

All Figures

	Fig. 1. MATISSE uv coverage from all 25 exposure cycles. Squares denote observations taken in 2020, and circles represent observations from 2021. The mean correlated flux between 4.5 and 4.9 μm is used as the color scale. North is up, and east is to the left.
In the text

Fig. 2. Excerpt sample of calibrated correlated fluxes for Circinus in the LM band. The top-left panel shows the 33 m “total flux” used to compute squared visibilities. The colors indicate the observing date: March 2020 (blue), February 2021 (yellow), and May 2021 (red). Black squares and black circles are the predicted correlated fluxes from the image reconstructions and Gaussian model fits, respectively. Continued in Fig. B.1.

In the text

	Fig. 3. Circinus 0.4″ fluxes in the L and M bands between 2000 and 2018. We compile the measurements of Prieto et al. (2004), Stalevski et al. (2017), Isbell et al. (2021), and this work. L-band values are plotted in blue, and M-band values are plotted in red. The mean flux for each band is plotted as a dashed line, with the standard deviation of the fluxes given as a shaded region.
In the text

Fig. 4. Images and models of the Circinus LM-band dust. Top: Image reconstructions at 3.7 and 4.7 μm. The contours are at 10 × σIm, estimated from the error maps produced by delete-d jackknifing the uv coverage. The FWHM of the estimated dirty beam is given in the bottom-right corner of each panel. Also included is the N-band continuum image from Paper I for reference. Bottom: Gaussian model fits to the same wavelengths, specifically the 1+ model in the L band and the two-Gaussian model in the M band. The dashed ellipse in the L-band model image represents the 1σ PA uncertainty of the fitted model. The colored circles in the bottom-right panel illustrate the extraction apertures for the measured fluxes.

In the text

Fig. 5. One-BB fits for the aperture-extracted Circinus LMN fluxes. The colors are the same as in Fig. 4, with D0 in red, DN40 in yellow, and DS40 in red. The top-right plot compares the D0 fluxes to the various RT model modifications. The fits using the N-band data alone are included for comparison. In aperture D0 the K-band measurement from Burtscher et al. (2015) is included as an upper limit for the near-infrared flux.

In the text

Fig. 6. Schematic of the dust within the RT models that best reproduce the observed dust temperatures. The basis of the model comes from Stalevski et al. (2019), and a rough parameter range is given in Paper I. Here we have added a cluster of dust clouds above the disk and behind the hyperbolic cone, at a position that in projection corresponds to the central aperture. These clouds could represent, for example, a puffed-up sublimation zone, freshly launched winds, or a smoother boundary for the hyperbolic cone. The spatial resolution of these observations is not sufficient to distinguish between these scenarios.

In the text

Fig. A.1. Comparison of DRS-extracted and custom-extracted correlated fluxes for Circinus (left) and a calibration star, HD120404 (right). From top to bottom, the panels show the cleaned fringe pattern (the Fourier transform of the cleaned interferogram), the six fringes extracted with the DRS, and the same six fringes extracted with our method. In the Circinus fringe pattern, there is significant flux in the M band that is missing in the DRS fringes but is present in our method.

In the text

	Fig. B.1. Continuation of Fig. 2.
In the text

	Fig. B.1. continued.
In the text

	Fig. B.1. continued.
In the text

	Fig. B.2. LM-band 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.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

	Fig. B.2. continued.
In the text

Fig. C.1. Two-BB fits for the aperture-extracted Circinus LMN fluxes. The colors are the same as in Fig. 4, with D0 in red, DN40 in yellow, and DS40 in red. There is little discernible improvement in fit quality with the addition of a second component, and the second component is in all cases highly extincted. The fits using the N-band data alone are included for comparison. In aperture D0, the K-band measurement from Burtscher et al. (2015) is included as an upper limit for the near-infrared flux.

In the text

	Fig. C.2. NACO flux-calibrated images for IRAS 14480-5828 and Circinus.
In the text

	Fig. C.3. VISIR flux-calibrated M-band image of Circinus with 0.4", 1", 2", and 4" apertures overlaid. Here we clearly see the extended polar dust structures to the east and west of the nucleus.
In the text