© The Authors 2024

1. Introduction

It is widely accepted that the spectral differences between type 1 (Sy1) and type 2 (Sy2) Seyfert active galactic nuclei (AGNs) can be attributed to an obscuring medium that blocks an observer’s line of sight (LOS) to the central engine and broad line region (BLR) in Sy2s, but not in Sy1s (Antonucci 1993). The obscurer is generally thought to be a dusty medium in the equatorial direction, such as a thick disk, which obscures the central region when the AGN is viewed edge-on, but not face-on. The medium is also thought to be clumpy in order to both explain how a thick disk can be vertically supported (Krolik & Begelman 1988) and observed short timescale variations of obscuration in X-rays (Risaliti et al. 2005). The medium then also acts as a material reservoir to fuel the central engine and emission regions (Hönig 2019; Hickox & Alexander 2018). Given the tremendous amounts of energy released by the central engine, some of the dusty material could be lifted out of the disk and into a dusty outflow through radiation pressure, opening up a pathway for AGN feedback to the host galaxy (e.g., Williamson et al. 2019; Wada 2012). Therefore, understanding this dusty structure is crucial to understanding AGN inner dynamics and interactions between AGNs and their host galaxy.

First-generation long baseline IR interferometry has offered a powerful new tool to AGN science by resolving the expected spacial scales of the obscuring dusty material at the wavelengths that it thermally emits. It has provided strong evidence for polar dust structures, associated with a dusty outflow (López-Gonzaga et al. 2016; Hönig et al. 2013; Burtscher et al. 2013) in the mid-IR (N band, 8 − 13 μm) by the Mid-Infrared Interferometer (MIDI, Leinert et al. 2003) and ring or disk-like structures (Weigelt et al. 2012; Kishimoto et al. 2009, 2011) in the near-IR (K-band, ∼2.2 μm) by the Astronomical Multi-BEam combineR (AMBER, Petrov et al. 2007) and Keck (Colavita et al. 2013). The discovered geometry was in agreement with (and helped develop upon) our understanding of the dusty material leading to an iteration of the unification scheme that contains an obscuring equatorial dust-gas disk and an outflowing dusty wind (Hönig et al. 2012).

Second-generation IR interferometers improved upon the first through increased uv coverage per observation and improved phase information. In particular, AMBER demonstrated the application of closure phases in IR interferometry, a lesson which GRAVITY (GRAVITY Collaboration 2017) and the Multi AperTure mid-Infrared SpectroScopic Experiment (MATISSE, Lopez et al. 2022) have well received. While first-generation instruments such as MIDI and AMBER have provided model derived geometry for AGNs, it is evident that MATISSE and GRAVITY can further allow for image reconstructions (e.g., for AGNs: Isbell et al. 2023, 2022; Gámez Rosas et al. 2022; GRAVITY Collaboration 2021, 2020b).

One such object is the Sy2 AGN, named NGC 1068. The GRAVITY data of NGC 1068 at 2.0 − 2.3 μm was reconstructed and a ring-like structure was found (GRAVITY Collaboration 2020b). The structure had multiple interpretations in the work; however, the two favoured geometries were either: the ring is dust near the sublimation radius or the ring is coincidental and consists of clumpy dust on the illuminated back wall of the obscuring disk or inner outflow. Alternatively, a radiation transfer (RT) modelling attempt of the GRAVITY data suggested a thin ring perpendicular to the accretion plane (Vermot et al. 2021), which has been interpreted as possibly a tidal disruption event unrelated to the dust structure of the unification scheme.

Gámez Rosas et al. (2022) used MATISSE to observe and reconstruct images for NGC 1068 from 3 to 13 μm and found a similar structure to that seen by GRAVITY at the shortest wavelengths of MATISSE. However, the wider wavelength range provided the possibility to extract thermal information from the images derived at different wavelengths. By assuming the MATISSE and GRAVITY images at each wavelength were aligned by their matching spatial flux distributions and brightest spots, an SED was extracted for different apertures aligned to the observed dust structures, and a two black body and dust obscuration model was fit. The image alignment had the further constraint that the extracted SEDs must be continuous. This modeling suggested that the dust in the ring feature observed with GRAVITY was too cool (≤600 K) to be sublimating (T_sub ≈ 1500 − 1800 K). The caveat of this interpretation is that the inter-wavelength image position registration was assumed by geometry and SED instead of known. The absolute image position would require absolute phase information which is destroyed by the atmosphere. To obtain the relative alignment between bands instead would require inter-band differential phases that are not currently available. The uncertainty leaves room for ambiguity over which interpretation is correct. A different registration could give rise to different temperatures and different positions of the dust within the AGN.

Furthermore, GRAVITY Collaboration (2020b) and Gámez Rosas et al. (2022) have provided multiple alignments between the image reconstructions and the H₂O maser emission of NGC 1068. The most recently published study has come from Gallimore & Impellizzeri (2023), who also provided additional alternate alignments. This has made it difficult to build a clear understanding of the coincidence between the maser and dust emission structures. It is clear that uncertainty in the image positions is a major limitation when studying structures at such high angular resolutions.

It is the aim of this work to provide an interpretation of the IR interferometric data that is free of a “visually based” image alignment via the modeling of MATISSE and GRAVITY data simultaneously with a chromatic model compatible with the unified scheme of AGN. Specifically, we have designed and used the model to: attempt to explain the complex IR interferometric data and reconstructed images of the MATISSE (Gámez Rosas et al. 2022) and GRAVITY (GRAVITY Collaboration 2020b) observations, determine whether the recovered model was compatible with the unification scheme of AGNs, and use the best fit models to produce a model-dependent image alignment and black hole location.

The derived geometry will be used for a full RT modelling of the dusty structure in NGC 1068 in future work. The model we present will also be applied to a larger sample of AGNs to create a more complete geometric description than had previously been possible before second-generation VLTI instrumentation.

Overall, RT models have been used to interpret the SEDs of AGNs successfully for many years (e.g., González-Martín et al. 2019; Stalevski et al. 2012; Hönig & Kishimoto 2010; Nenkova et al. 2008, 2002); they are flexible and robust tools for determining model dependent properties of dust in AGNs. Combining the SED fitting process with simultaneous fitting of IR interferometric data has proved to be a very powerful tool with greater constraints on the dust geometry than the SED alone (Hönig & Kishimoto 2017). SED modeling provides zero-order geometric information, the distance and distribution of the dust relative to the central engine can be derived but the absolute positions are ambiguous due to the unresolved nature of the information. By providing resolved information, distinctions can be made between otherwise well fitting models (e.g., Isbell et al. 2022; Leftley et al. 2019, 2018; Hönig & Kishimoto 2017).

The results of models such as CAT3D-WIND (Hönig & Kishimoto 2017) when simultaneously reproducing interferometric and SED data as well as the recent work of Isbell et al. (2022), Gámez Rosas et al. (2022) demonstrate the usefulness of combing spatially resolved observations with chromatic information when modelling dust structures. While SED modelling is well established, it is only after the first generation of VLTI IR interferometers were applied to AGNs was it possible to do this on milliarcecond scales in the IR where the dust thermally emits (e.g., Stalevski et al. 2019). With the second generation of instrumentation, it is now possible to apply these techniques to a much wider range of wavelengths and in greater detail due to improved uv coverage per observation and phase information. We take inspiration from the simultaneous RT and interferometric modeling of previous works and create a chromatic model for AGNs.

In Sect. 2, we describe the observations we used in the modeling. In Sect. 3, we provide details on how the model was constructed and how we created brightness distributions per wavelengths. In Sect. 4, we present a comparison of our model to the data and we determine the best set of parameters to reproduce the data. In Sect. 5, we provide the results of the modeling described previously with some discussion into what each set of results indicate about NGC 1068 alone. In Sect. 6, we discuss the results in the context of the wider literature. In Sect. 7, we discuss the different fitting methods and prospects for future work. Finally, in Sect. 8, we give our final remarks and summarize the main points of the work.

2. Observations

This work makes use of published observations of NGC 1068 made with GRAVITY and MATISSE on the VLTI (priv. comm.)¹. The detailed description of the observations and reduction for the GRAVITY and MATISSE data can be found in GRAVITY Collaboration (2020b) and Gámez Rosas et al. (2022), respectively. The raw data can be found on the ESO archive under the following IDs: MATISSE programme ID: 0103.C-0143, 0104.B-0322 GRAVITY programme ID: 0102.B-0667, 0102.C-0205, and 0102.C-0211.

We use the same data as the previous works to remove the possibility of any differences from data handling. The only difference we make is that we did not use the Auxiliary Telescope (AT) data given in GRAVITY Collaboration (2020b) because combining Unit Telescope (UT) and AT data is a non-trivial problem due to the different FoVs between the telescopes.

3. Model

To interpret the observational data, we created a semi-RT model of AGN dust morphology. The morphology we chose is based on the disk+wind interpretation of the unified model (e.g., Hönig et al. 2012) given its prior success in reproducing IR interferometric observations of AGNs such as with CAT3D-WIND (Hönig & Kishimoto 2017). CAT3D-WIND is a successful radiative transfer implementation of the disk+wind model. In their work, they used CAT3D-WIND to model the IR SED of many AGNs, including those that CAT3D (a clumpy torus model) was not able to. It has also been used to model MIDI observations of AGNs successfully (e.g., Leftley et al. 2018). More recently, it was used in one of the interpretations of the GRAVITY data we use in this work (GRAVITY Collaboration 2020b, model 4). Outside of CAT3D-WIND, the geometry has also been used to explain Circinus interferometric observations with RT models (Isbell et al. 2022; Stalevski et al. 2019). Furthermore, the disk+wind model is backed up by hydrodynamic and numerical simulations that find a dusty wind can be easily formed from an accreting dust disk around a central engine (Williamson et al. 2020; Venanzi et al. 2020; Wada et al. 2016; Wada 2015, 2012). IR radiation from the dust disk itself then pushes the wind vertically through radiation to form a hollow hyperbolic cone. Therefore, we decided to use this geometry in our model.

3.1. The geometry

We interpret the disk+wind model in terms of two components: a disk with an opening angle and inner radius, and two hyperbolic hollow cones that are perpendicular to the disk. We define our initial geometry, such that z is perpendicular to the plane of the disk, namely, the polar axis of the AGN and (x, y) are in the disk plane. For convenience, we also define z as north and x as east.

We define the disk to have an inner radius (r_in) and an opening angle 2ϕ where ϕ is a free variable (the half-opening angle); the disk opens from the inner radius. In a 2D plane, the surface of the disk can be defined as:

$$\begin{array}{cc}\hfill z& =\pm (\sqrt{x^2+y^2}-r_{\mathrm{in}})tan\varphi ,\hfill \\ \hfill & =\pm (r-r_{\mathrm{in}})tan\varphi ,\mathrm{for}r>r_{\mathrm{in}},\hfill \end{array}$$(1)

where $r=\sqrt{x^2+y^2}$ and the half opening angle ϕ is defined from the disk plane (x, y).

We define the wind to be a hyperbolic cone described by:

$$\begin{array}{cc}\hfill z& =\pm a_{\mathrm{w}}(\sqrt{{\left(\frac{r}{b}\right)}^2+1}-\sqrt{{\left(\frac{r_{\mathrm{in}}}{b}\right)}^2+1}),\hfill \\ \hfill & \mathrm{for}r>r_{\mathrm{in}}\hfill \end{array}$$(2)

where a_w and b are fitted variables that define the hyperbolic cone. These two components are the base geometry of our model. A schematic of the disk and wind can be found in Fig. 1.

Fig. 1. Schematic of the disk described in Eq. (1) and wind defined in Eq. (2). We provide the example aw and b used to make the wind schematic.

3.2. Applying clumpiness

We want our structure to be clumpy so for each component we distribute a number of clumps. We define two variables for the number of clumps. N_p we give as the total number of clumps in the system and f_w is the log₁₀ fraction of N_p that is attributed to the wind. Therefore, we can derive the number of clumps in the disk as N_d = (1−10^f_w)N_p and the number of clumps in the wind as N_w = 10^f_wN_p.

To distribute clumps in the disk, we do so randomly with a probability based on the geometry in Eq. (1). The probability in each dimension is:

$$\begin{array}{cc}\hfill P(x,y),P\left(z\right)& =0\mathrm{for}r<r_{\mathrm{in}},\hfill \\ \hfill P(x,y)& =Aexp-\frac{x^2}{2{\beta }^2}exp-\frac{y^2}{2{\beta }^2},\mathrm{for}r>r_{\mathrm{in}}\hfill \\ \hfill P\left(z\right)& =\frac{3}{\sqrt{2\pi }(r-r_{\mathrm{in}})tan\varphi }exp-\frac{{\left(3z\right)}^2}{2{(r-r_{\mathrm{in}})}^2{tan}^2\varphi },\hfill \end{array}$$(3)

where A is the normalization of the truncated Gaussian and β is the standard deviation of the dust distribution in the plane of the disk. For P(z), we defined the disk height to be 3σ of the normal distribution.

Along the wind, we distribute N_w clouds along the hyperbolic cone in a similar manner to the disk. To give the wind a thickness, we randomly offset a_w for each clump in the wind by a normal distribution with a standard deviation of a_width, which can be a free parameter but is fixed in this work, and centred on a_w. The spatial probability distribution for clumps in the wind structure is therefore:

$$\begin{array}{cc}& P(x,y),P\left(z\right)=0,\mathrm{for}r<r_{\mathrm{in}}\hfill \\ \hfill & P(x,y)\propto r^{-{\alpha }_{\mathrm{w}}}\mathrm{for}{r}_{\mathrm{in}}<r<r_{\mathrm{out}},\hfill \\ \hfill & P\left(a_{\mathrm{w}}\right)\propto exp-\frac{a_{\mathrm{w}}^2}{2a_{\mathrm{width}}^2},\hfill \end{array}$$(4)

where α_w is a power law for the distribution of clumps in the wind in the (x, y) plane, and r_out is the outer radius defined as four times larger than the image field of view (FOV) of the longest wavelength being modelled. r_out is arbitrarily large and defined to save computation time by excluding points that fall too far away from the centre to contribute to the total AGN IR flux. When computing z from Eq. (2) for a clump, each point has an equal probability of being positive or negative to create a biconical outflow. An example schematic can be seen in Fig. 2.

Fig. 2. Schematic of the disk described in Eq. (1) and wind defined in Eq. (2) with clumps distributed with the probabilities described in Eqs. (3) and (4). We provide the example β and αw used to make the wind schematic.

We then incline the system towards the viewer (about the x-axis) by inc and rotate in the image plane by ang (about the y-axis); ang is equivalent to the on-sky from north through east position angle (PA) of the polar axis that is the axis perpendicular to the disk. In an aligned system, this is the direction of the jet.

3.3. Defining a clump

In the previous section, we describe how we distribute clumps in the geometry. Here we define what a clump is. We define all clumps as identical spheres of radius r_c, where the only unique properties of an individual clump is the position and temperature.

Each clump is given two temperatures; this is because AGNs are centrally powered. If a clump is optically thick, then the side facing the central engine is directly heated to T_p, while the other side is shielded and therefore at a lower temperature T_b. For the directly heated side of the clump, we set the temperature as a broken power law based on distance from the central engine. We set the power and break separately for each component. For the wind, we set:

$$\begin{array}{cc}\hfill T_{\mathrm{p}}& =T_0{r}_{\mathrm{in}}^{\alpha }{R}^{-\alpha },\mathrm{for}R<r_{\mathrm{in}}{\alpha }_{\mathrm{r}}\hfill \\ \hfill & =T_0{r}_{\mathrm{in}}^{\alpha +{\alpha }_{\mathrm{off}}}{\alpha }_{\mathrm{r}}^{{\alpha }_{\mathrm{off}}}{R}^{-(\alpha +{\alpha }_{\mathrm{off}})},\mathrm{for}R>r_{\mathrm{in}}{\alpha }_{\mathrm{r}}\hfill \end{array}$$(5)

where $R=\sqrt{x^2+y^2+z^2}$, α is the power before the break, α_off is the additional power after the break, and α_r is the break radius in units of r_in. In the disk, we set:

$$\begin{array}{cc}\hfill T_{\mathrm{p}}& =T_0{r}_{\mathrm{in}}^{\gamma }{R}^{-\gamma },\mathrm{for}R<r_{\mathrm{in}}{\gamma }_{\mathrm{r}}\hfill \\ \hfill & =T_0{r}_{\mathrm{in}}^{\gamma +{\gamma }_{\mathrm{off}}}{\gamma }_{\mathrm{r}}^{{\gamma }_{\mathrm{off}}}{R}^{-(\gamma +{\gamma }_{\mathrm{off}})},\mathrm{for}R>r_{\mathrm{in}}{\gamma }_{\mathrm{r}}\hfill \end{array}$$(6)

where γ is the power before the break, γ_off is the additional power after the break, and γ_r is the break radius in units of r_in. In this work, we do not fit the SED of NGC 1068 which makes the result insensitive to the absolute temperature; only the relative temperature is required. Therefore, we set the power law normalization of T₀ to 1800 K at r_in.

The second temperature, T_b, is also set separately for the disk and the wind. We define this temperature as a variable constant for each component. For the disk, we define the variable as T_bd and for the wind we define it as T_bw.

3.4. Optical depth of a clump

We assume a clump is uniform in density and we can define the optical depth, τ, along the LOS for a clump as a function of (x,z). The equation for the optical depth is:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{c}}(r_{\mathrm{c}}^{\prime },\lambda )=\left(2{\tau }_{c0}\left(\lambda \right)\sqrt{r_{\mathrm{c}}^2-{r_{\mathrm{c}}^{\prime }}^2}\right)\mathrm{for}{r}_{\mathrm{c}}^{\prime }\le r_{\mathrm{c}},\end{array}$$(7)

where τ_c0(λ) is the wavelength dependent attenuation coefficient of the clump, $r_{\mathrm{c}}^{\prime }=\sqrt{{(x-x_{\mathrm{c}})}^2+{(z-z_{\mathrm{c}})}^c}$, and (x_c, z_c) is the (x, z) location of the clump.

3.5. Defining emission

Once we have a clump temperature and optical depth, we can define the emission. When defining emission for a clump, we evaluate per pixel of the image. A pixel has a coordinate (i,j), where i and j are x and y scaled by the pixel size in mas. We assume that a clump emits thermally as two black bodies such as $\left(1-\exp \left(-{\tau }_{\text{c}}\left(r_{\text{c}}^{\prime }\left(i,j\right),\lambda \right)\right)\right)B_{\nu }\left(\nu ,T_{\text{p}}\right)$ and $\left(1-\exp \left(-{\tau }_{\text{c}}\left(r_{\text{c}}^{\prime }\left(i,j\right),\lambda \right)\right)\right)B_{\nu }\left(\nu ,T_{\text{b}}\right)$ where B_ν is Planck’s law. For T_p, we calculate the total flux from this temperature for a clump as $F_{\text{p}}\left(\lambda \right)\propto \left(1-\exp \left(-{\tau }_{\text{c}}\left(r_{\text{c}}^{\prime },\lambda \right)\right)\right)B_{\nu }\left(\nu ,T_{\text{p}}\right)$. However, we defined this temperature component as the side that faces the central engine and at some viewing angles the clump will self-obscure if optically thick. We approximate this phase angle effect (sometimes referred to as moon phase effect) by applying a factor to the flux of a clump based on its position relative to the centre of the model and the viewer:

$$\begin{array}{cc}\hfill {\theta }_l& =arctan\frac{z_{\mathrm{c}}}{-y_{\mathrm{c}}},\hfill \\ \hfill {\theta }_m& =arctan\frac{-y_{\mathrm{c}}}{x_{\mathrm{c}}}.\hfill \end{array}$$(8)

Therefore, the flux factor can be defined as S_moon = 0.5(cos θ_l sin θ_m+1). This is the same treatment that is used in some BLR cloud models (Rakshit et al. 2015). The effect scales with optical thickness, so we apply a factor of s_c(1−exp(−1.5r_cσ_c(λ))) to S_moon, where s_c is a free variable scaling, which scales the influence of the phase angle effect (P_a) with optical thickness. So F_p can be fully defined as:

$$\begin{array}{cc}& S_{\mathrm{moon}}=0.5(cos{\theta }_{l}sin{\theta }_m+1),\hfill \\ \hfill & P_{\mathrm{a}}\left(\lambda \right)=s_{\mathrm{c}}{S}_{\mathrm{moon}}(1-exp(-{\tau }_{\mathrm{c}}(1.5r_{\mathrm{c}},\lambda ))),\hfill \\ \hfill & F_{\mathrm{p}}\left(\lambda \right)=P_{\mathrm{a}}\left(\lambda \right)(1-exp(-{\tau }_{\mathrm{c}}(i,j)))B_{\nu }(\nu ,T_{\mathrm{p}}).\hfill \end{array}$$(9)

For the flux derived from the base temperature, we do not need to apply this phase angle effect. Therefore, the total flux contribution from the cooler side is:

$$\begin{array}{c}\hfill F_{\mathrm{cool}}\left(\lambda \right)=(1-exp(-{\tau }_{\mathrm{c}}(i,j)))B_{\nu }(\nu ,T_{\mathrm{b}}).\end{array}$$(10)

We note that the flux is per Steradian. The phase angle effect accounts for the differences in emitting area for the two flux components so they can be additively combined before accounting for surface area. We account for this area when creating the image.

3.6. Applying obscuration

When creating an image, we apply two forms of obscuration. The first is from the clumps which we have already defined. A clump from the viewers perspective attenuates the flux behind it by $1-\exp \left(-{\tau }_{\text{c}}\left(r_{\text{c}}^{\prime },\lambda \right)\right)$. The second form of obscuration is from a smooth dust structure (the blue solid lines in Fig. 2).

The smooth dust structure we define as a smooth medium which is only present in the disk (see Eq. (1)) and only acts in obscuration. This obscuration can be applied analytically. For simplicity, we define the obscuration in the plane of the disk such that (x_r,y_r,z_r) are the same as the (x,y,z) in Fig. 1 and $r_{\mathrm{r}}=\sqrt{x_{\mathrm{r}}^2+y_{\mathrm{r}}^2}$, namely, before the geometry was rotated by inc and ang. We chose the optical depth of the smooth structure to decrease in the (x_r,y_r) plane and be uniform in z_r such that:

$$\begin{array}{c}\hfill \begin{array}{cc}\tau (x_{\mathrm{r}},y_{\mathrm{r}},z_{\mathrm{r}},\lambda )\hfill & \propto {\tau }_0\left(\lambda \right)exp-\frac{x_{\mathrm{r}}^2}{2{\beta }^2}exp-\frac{y_{\mathrm{r}}^2}{2{\beta }^2},\hfill \\ [4pt]\hfill & \hfill & \mathrm{for}{r}_{\mathrm{r}}>r_{\mathrm{in}}\hfill \\ [1.5pt]\hfill & \hfill & \mathrm{and}(r_{\mathrm{r}}-r_{\mathrm{in}})tan-\varphi <z_{\mathrm{r}}<(r_{\mathrm{r}}-r_{\mathrm{in}})tan\varphi \hfill \\ [4pt]\tau (x_{\mathrm{r}},y_{\mathrm{r}},z_{\mathrm{r}},\lambda )\hfill & =0,\mathrm{elsewhere},\hfill \end{array}\end{array}$$(11)

to mimic the (x_r,y_r) clump distribution in the same structure where τ₀(λ) is the attenuation coefficient at r_r = 0 not accounting for r_in. If we rotate the smooth obscuring geometry by applying the rotation matrix:

$$\begin{array}{c}\hfill \left[\begin{array}{c}y\hfill \\ z\hfill \end{array}\right]=\left[\begin{array}{cc}cos\left(inc\right)\hfill & -sin\left(inc\right)\hfill \\ sin\left(inc\right)\hfill & cos\left(inc\right)\hfill \end{array}\right]\left[\begin{array}{c}{y}_{\mathrm{r}}\hfill \\ z_{\mathrm{r}}\hfill \end{array}\right]\end{array}$$(12)

the LOS from an emitter to the viewer is along y only allowing for a 1D integration. If the LOS from an emitter enters the disk at (x₁,y₁,z₁) and exits at (x₁,y₂,z₁) then the obscuration from the smooth structure can be defined as:

$$\begin{array}{cc}\hfill \tau (\mathit{inc},\lambda )& ={\tau }_0\left(\lambda \right)\sqrt{1+{tan}^2\left(\mathit{inc}\right)}exp(-\frac{x_1^2}{2{\beta }^2})\hfill \\ \hfill & \times {\int }_{y_1}^{y_2}exp(-\frac{y^2}{2{\beta }^2})\mathrm{d}y.\hfill \end{array}$$(13)

This can be redefined in terms of the error function:

$$\begin{array}{cc}\hfill \tau & ={\tau }_0\left(\lambda \right)\sqrt{1+{tan}^2\left(\mathit{inc}\right)}exp(-\frac{x_1^2}{2{\beta }^2})\hfill \\ \hfill & \times {\left[\sqrt{\frac{\pi }{2}}\beta \mathrm{erf}\left(\frac{y}{\sqrt{2}\beta }\right)\right]}_{y_1}^{y_2},\hfill \end{array}$$(14)

where we use the standard definition of the error function:

$$\begin{array}{c}\hfill \mathrm{erf}\left(t\right)=\frac{2}{\sqrt{\pi }}{\int }_0^{z}exp(-t^2)\mathrm{d}t.\end{array}$$(15)

In the completely face-on case where inc = π/2, τ is ill-defined by the given parameterization. Instead, the limit as inc → π/2 should become the uniform dust case. We do not expect to encounter this Blazar-like case so we do not include it.

In reality, we do not make τ₀(λ) or τ_c0(λ) a model parameter, instead we define βτ₀(λ) = n_dσ(λ) and τ_c0(λ) = n_cσ(λ). n_d and n_c are the model parameters and are in units of g cm⁻² and g cm⁻² mas⁻¹, respectively. σ(λ) is the optical cross-section per mass in cm² g⁻¹ of the dusty material used. We use the same dust mix for the clumps and the smooth structure. n_c and n_d are density scaling factors.

For the dust we apply in this work, we chose a mix of amorphous carbon and olivine silicates (Min et al. 2005) as was done in Gámez Rosas et al. (2022). There, we allowed for the relative amount of carbon to olivines to be a free parameter, p_g, in the same way as in their SED fitting. Finally, σ(λ) results for the two dust species can be seen in Fig. 3.

Fig. 3. Optical cross-section per gram of the dust materials used in this work as well as a 20% amorphous carbon and 80% olivine mix (pg = 0.2).

3.7. Model evaluation specifics

We have fully defined our model geometry and described the characteristics of its emission and absorption. However, there are some important assumptions we make when comparing the model to data. Here, we discuss how we evaluate the model.

3.7.1. PSF effects

There is one final flux scaling factor that we apply when determining fluxes that is due to the telescope PSF. The scaling factor is to account for the flux drop from the telescope PSF and is approximated by a 2D Gaussian on the plane of the sky with a FWHM of λ/D centered on the model centre. This reduces the flux of a point by the distance from the centre. This approximation does not account for the effects from fibre injection in GRAVITY and the pinhole and pupil stops of MATISSE but it is sufficient for our image sizes. It is noteworthy that we included this factor to ensure any bright points far from the centre are treated properly, but the effect is usually negligible due to the PSF of a UT. For a telescope diameter of 8 m, the FWHM at 2 μm is ∼50 mas, which is 1.7 times the FOV at the same wavelength.

3.7.2. Making an image

Before comparing our model to data, we must create an image per wavelength. To make our image, we take our defined model, with the analytical smooth obscuration applied, and determine the clump emission and attenuation for all clump occupied pixels in the image. We assume that the obscuration at the centre of a pixel is a good approximation of the average value in the pixel when evaluating the clumps emission and obscuration. For pixels that are only partially filled, we reduce the flux by the fraction of the pixel that is unfilled. To account for changes in pixel size between wavelengths, we normalize the clump flux such that:

$$\begin{array}{c}\hfill {\int }_{-r_{\mathrm{c}}}^{r_{\mathrm{c}}}{F}_{\mathrm{c},\lambda }\mathrm{d}{r}_{\mathrm{c}}^{\prime }=F_{\mathrm{p},\lambda }+F_{\mathrm{cool},\lambda },\end{array}$$(16)

in the case where τ → ∞. So for any given pixel we get:

$$\begin{array}{c}\hfill F_{i,j,k}=F_{i,j,k-1}exp(-{\tau }_c\left(r_{cijk}^{\prime }\right))+f_{\mathit{ijk}},\end{array}$$(17)

for pixel ij and clump-in-pixel k, where k ranges from 1 to the number of clumps in the pixel. Furthermore, F₀ = 0, f_ijk is the flux of clump k for pixel ij, and $r_{cijk}^{\prime }$ is the in-clump radius $r_{\text{c}}^{\prime }$ for clump k at pixel ij. Should a clump fall outside an image, the flux is not included within the image but it is accounted for.

3.7.3. Over-resolved structure

When observing extended objects with interferometry, flux can exist outside the largest angular scale the interferometer probes but still within the single dish photometric aperture. We define this as over-resolved structure. In AGNs, this can be from the dust itself or unrelated emission from the host galaxy. Instead of images, IR interferometry provides observable quantities that can be used to reconstruct the target source. The two we model are the squared visibility (V²) and closure phase (Φ_c). We compare the model V² and Φ_c, calculated from the model image using the python software package GALARIO, to the observed V² and Φ_c, respectively. Any over-resolved flux that falls outside the image does not contribute to Φ_c but it does contribute to V² as a PA independent scaling factor. V² is defined as the squared flux observed by the interferometer divided by the squared flux observed by a single telescope. If we assume all the flux observed by the telescope falls within the image when it does not, the single telescope model flux will be lower than it should be and, in turn, V² will be higher.

Additional scaling to V² can arise from instrumental effects causing flux losses in the interferometer. Therefore, when modelling a single wavelength or a small band it is common to apply a scaling factor to the model V² such that V² = (V₀V_model)² (e.g., Gámez Rosas et al. 2022; GRAVITY Collaboration 2020a; Leftley et al. 2018). This achromatic factor has been shown to be a good approximation for a single wavelength. However, it will fail to capture any chromatic over-resolved structure. Because our model is chromatic with multiple instruments, applying an appropriate scale factor is not trivial.

To simplify the problem, we make the assumption that the instrumental losses are insignificant or do not change significantly between GRAVITY and MATISSE; this is a relatively safe assumption because NGC 1068 is bright for an AGN and the data for both instruments has been very well reduced by the respective teams. Furthermore, the greatest source of losses for AGN IR interferometry, particularly at shortest wavelengths, is thought to be the MACAO AO system (Leftley et al. 2021). The GRAVITY data was observed under excellent conditions (optical coherence time of 7–13 ms GRAVITY Collaboration 2020b) and the MATISSE data was filtered to remove bad AO performance (Gámez Rosas et al. 2022). Therefore, the losses from AO should be minimal. For fainter AGNs this may not be appropriate and will need to be readjusted. Therefore, we just account for sources of over-resolved emission. We define our over-resolved structure with the following three components: over-resolved model flux, a hot component, and an achromatic component.

3.7.4. Over-resolved model flux

Some of this over-resolved flux is inherently accounted for in our model. We sum the flux of all clumps that fall outside an image which we use to scale the model V². For computational efficiency, we only applied clump obscuration for clumps inside the image at that wavelength. We assume that clumps that fall outside the image are too sparsely distributed to significantly overlap each other, therefore, we can sum just the emission from these clumps after only smooth obscuration. If we sum the flux inside and outside the image, we receive the total flux of the model AGN. This is the direct contribution of the model dust to the over-resolved structure, the following components are not necessarily related to our dust structure, but still contribute flux in the IR.

3.7.5. A hot component

We found that we do still require an additional factor to explain the shortest and longest wavelengths simultaneously. We apply a hot component to the over-resolved flux. Unlike the other two over-resolved components, we only apply this one when performing multi-band fitting.

If we assume the simplest case that the structure responsible for the shortest baseline (∼40 m or ∼18 Mλ) V² at 2.05 μm is the same for the equivalent baseline at 3.29 μm then we find V² doubles between the two wavelengths (V² increases from 0.05 to 0.1 between 2.05 μm and 3.29 μm, see Fig. 4). We fit a simplistic model of $V\left(\lambda \right)=V_{\mathrm{b}}\frac{(1-\mathrm{frac}){\mathrm{BB}}_{\mathrm{under}}(\lambda ,T_1)}{\mathrm{frac}{\mathrm{BB}}_{\mathrm{over}}(\lambda ,T_2)+(1-\mathrm{frac}){\mathrm{BB}}_{\mathrm{under}}(\lambda ,T_1)}$, where BB is a black body of temperature T normalized to 1 at 2 μm for an under- and over-resolved component; frac is the relative flux of the under-resolved component at 2 μm; and V_b is an achromatic visibility, to visibilities between 16 Mλ and 20 Mλ from λ = 2 μm to λ = 5 μm. We find that the change can be explained by T₂ ∼ 20 000 K, T₁ ∼ 440 K, frac ∼ 0.3, and V_b ∼ 0.26 as seen in Fig. 4 as the red line. The fit tells us that the over-resolved component can be explained by a hot black-body of 20 000 K responsible for 30% of the flux at λ = 2 μm. This does not rule out that the component could be colder with some change in visibility originating from a more complex structure than the unchanging single temperature geometry source assumed in the quick fit.

Fig. 4. Visibility of NGC 1068 between 16 Mλ and 20 Mλ compared to wavelength in purple. Overplotted is a best fit model of an over- and under-resolved black-body model in red.

As such, we introduce an over-resolved hot black body of 20 000 K to our model when fitting multiple bands. We scale the flux from this black body based on the AGN flux such that F_BBh(λ) = BB_hotF_AGN(λ_min)B_λ (λ,20 000 K)/ B_λ(λ_min,20 000 K). λ_min is the shortest wavelength modelled. We fix the temperature because we find it to be otherwise poorly constrained.

There are multiple possibilities for the physical origin of this component. Under the assumption that the component is the Raleigh-Jeans tail of a hot black body, while we cannot constrain the exact temperature, we can conclude that it is too hot to be thermal emission from dust and not a missing component of our AGN thermally emitting dust model; namely, it behaves as a power-law against wavelength with a power of −2 in our covered wavelength range. Previous work using Keck by Weinberger et al. (1999) finds that 50% of the K-band nuclear (sub 0″.1) flux comes from an extended region not attributed to direct heating by the AGN but instead scattered nuclear emission. Weigelt et al. (2004) followed this with K-band speckle interferometry which found an elongated Gaussian structure with a lower limit FWHM of 26 × 58 mas responsible for 54 ± 5% of the nuclear K-band flux with the remaining from a compact nuclear source. The extended component in their work is proposed to be a mix of scattered light from dust in an outflow and thermal emission. Such extended scattered light from the accretion disk would be sufficient to explain the over-resolved emission at 2 μm. Alternatively, star formation could reach sufficient temperatures of ∼10 000 K and Lopez-Rodriguez et al. (2015) estimate that 45 ± 16% of the polarized K-band flux in the central 0″.5 arises from stellar contamination.

Under the assumption that the component’s emission is non-thermal, we find that while a power law with wavelength is just as good a description of the data, the power is completely degenerate with the temperature of the unresolved component with an upper limit on the power of 0. This means the non-thermal radio emission, with a power of λ^0.3 is also a possible contributor (Wittkowski et al. 1998). In reality, the over-resolved structure is likely to be from a combination of sources.

In future work, the inclusion of the SED when fitting or further investigation into the central region of NGC 1068 with high angular resolution imaging may provide better constraints on the temperature and physical origin of this component. We note that, because we do not model absolute fluxes, the exact spectral shape of the component beyond being more influential in K- than L-Band will have minimal impact on our model result.

3.7.6. An achromatic component

We additionally introduce an achromatic term, V₀. We find we require some scaling at longer wavelengths that is not covered by the hot black body. We apply this scaling as a wavelength independent multiplicative factor to the visibility at all wavelengths. This could be replaced by a second cooler over-resolved black body in future versions of the model. If it is indeed physical (and not instrumental), the origin of this factor could be from an extended larger smooth structure, such as the northern component 3 revealed with MIDI data by López-Gonzaga et al. (2014). The component from their work is suggested to be a cool 300 K optically thin thermally emitting component but it is also mentioned that it could be a warmer (∼700 − 800 K) not heated by the central engine or instead part of the radio emission. The last possibility is not preferred because their component 3 is not well-aligned with radio emission. Additionally, Rouan et al. (2004) find 475 K dust at 450 mas from the centre from high angular resolution observations with NACO. The spectral shape of this component is therefore unclear and the visibility is sufficiently approximated by the simpler achromatic term. With the inclusion of AT data or the SED, it may be possible to constrain this component further. This will not have significant impact in this work and can be investigated in more detail in future work.

3.8. Defining terms

When performing the modeling, our choice of parameters differs in some cases from our definitions used to describe the model above. Here, we summarize the parameters we used when modeling and how they relate to the model described above.

In summary, the model has 27 free parameters, 11 of which are geometric parameters and 13 of which relate to the point source fluxes; we also have 1 additional fixed geometric parameter. The parameters, with descriptions, can be found in Table 1. The remaining two parameters are: lnf_v and lnf_c; here, f is the fractional amount for which the variance is underestimated by the likelihood function if the errors were assumed correct (Foreman-Mackey et al. 2013) and is defined separately for closure phase and V² as f_c and f_v, respectively. We define the likelihood for one observable as

$$\begin{array}{cc}\hfill 2ln\left(L\right)& =-\sum _m^N[\frac{{(y\left(x_m\right)-\mathrm{model}(x_m,\alpha ))}^2}{{\sigma }_m^2+f^2\mathrm{model}{(x_m,\alpha )}^2}+\hfill \\ \hfill & ln\left(2\pi ({\sigma }_m^2+f^2\mathrm{model}{(x_m,\alpha )}^2)\right)].\hfill \end{array}$$(18)

4. Method

We use the model described in Sect. 3 with the PYTHON Markov chain Monte Carlo (MCMC) fitting code EMCEE (Foreman-Mackey et al. 2013). We make use of the GALARIO (Tazzari et al. 2018) and ASTROPY (Astropy Collaboration 2013) libraries. While the distribution of points is different, the model fitting process is the same as we used in the point source model of Leftley et al. (2021) and Gámez Rosas et al. (2022). However, we do not penalize the addition of points because, unlike in the previous works, additional points do not require additional parameters.

Our model is polychromatic and can be fit to the entire wavelength range of GRAVITY and MATISSE simultaneously. However, the wider the wavelength range, the more likely we are to encounter complex chromatic structures not covered by our simple model. Therefore, we run the model in two different modes, an individual band fit and a full polychromatic fit. We can then compare the results from each to infer the presence of other structures if such structures are present.

Because we need to calculate an image cube of models when fitting to compare with the observables at each wavelength, it is necessary to bin the observational data in the wavelength direction when fitting due to computational resource limitations. However, too large a binning removes the chromatic information contained within a single band. Therefore, we decided on the following binning in each band: from 3.2 μm to 3.9 μm in 4 bins for the L band, 4.6 μm to 4.8 μm in two bins for the M band, and 8.40 μm to 9.0 μm in 1 bin, 9.40 μm to 11 μm in 3 bins, and 11.4 μm to 12.7 μm in two bins for the N band. The K band was unbinned due to the rapid change in structure in the band coupled with the low spectral resolution from the fringe tracker data. We do not use the shortest wavelength bin due to the known contamination by the metrology laser (GRAVITY Collaboration 2020b) so we have four wavelength points in the K band. We used the same binning in both the band-by-band and polychromatic modeling.

4.1. Band-by-band

At the shortest wavelengths, we are sensitive to the relative position of individual points; namely, we resolve the non-uniform distribution of dust in the AGN. Due to the random nature of the point distribution this means we have a dependence on the random seed. We expect this to be particularly influential in the edge-on cases because the inner radius ring-like distribution dominates in the face-on case which becomes smooth at smaller N_p than the more disperse wind distribution. To evaluate the effect per band, we perform the fit multiple times for different seeds in the band-by-band fit. From this we can determine the specific best fit, the median best fit, and the uncertainty in the model and the parameters introduced by randomness. We also perform the fit twice, once imposing that the central region be obscured (nc < ϕ) and once imposing an unobscured central region labelled Sy2 case and Sy1 case, respectively. When we do not impose this, we find that inclination is linearly correlated with wavelength which is potentially unphysical.

We defined two sets of boundaries (given in Table 1) and fit the model at the wavelength ranges given in Sect. 4. We randomly selected a list of unique random seeds and run an MCMC fit for each band for each seed. The initial distribution of walkers is distributed uniformly over the initial range given in Table 1, which are not the same as the boundaries given for the fit. The starting ranges were selected to cover the ranges to which the model is most sensitive, less computationally expensive, or most physically reasonable from hydrodynamic simulations (Williamson et al. 2019, 2020).

We selected the best result per wavelength from all samples of all seeds through Eq. (18). This becomes the starting point of a subsequent run of the same seeds with walkers distributed normally around the result with a standard deviation of 5%. Once fitted, when calculating the median fit model observables (V² and Φ_c) and their uncertainties, we randomly sampled the pool of walkers 1000 times with replacement and determine the observables for each choice. We then take the 16th, 50th, and 84th percentiles of the 1000 observables to find the error and median. The error on the parameters themselves were calculated for each seed from the 16th and 84th percentiles, with the given value being the 50th percentile not including burn-in samples. When giving average values of parameters all seeds together, we combined the sampler pools from all seeds before deriving the median and errors.

4.2. Polychromatic fit

We proceeded to fit all provided wavelengths simultaneously. We used the parameters found in the band-by-band modeling as prior information to set the initial walker distribution. Furthermore, we only fit for 1 seed. For the majority of the parameters, we used the results from the M band because this lays in the centre of our coverage. However, we selected r_c, s_c, and p_g from the K band as they are most sensitive to the shortest wavelengths. Furthermore, we selected n_d, n_c, ϕ, and β from the N band due to the silicate feature and larger FoV. While this was the procedure we decided upon, upon fitting, we found the increase in data coverage made the final model largely insensitive to the choice of initial walker position. Finally, similarly to the band-by-band fit, we performed a second longer modeling run from the best result. The results and statistics we present here come from the second run.

5. Results

In the initial modelling of the band-by-band fit, we found the Sy2 case to be a better fit that the Sy1 case by likelihood. Because the Sy2 case was preferred, we performed the final fit on this case only. As described in Sect. 4.1, we initiated the final fit for each seed from the best parameters of the best seed for each band. The initial modelling demonstrates interesting results beyond being an starting point for the final fit. This includes chromatic changes in geometry and comparisons between Sy1 and Sy2 cases. Therefore, we present and discuss these results in Appendix A.1.

5.1. Band-by-band fit

The final best fits for each seed compared to the data can be seen in Fig. C.1 and the images are shown in Fig. C.2. The final best-fitting parameters are plotted in Fig. C.3 as Xs.

It can be seen in Figs. C.2 and C.1 that the different seeds produce visually and interferometrically similar models. The similarity enforces the idea that the seed does not significantly impact the result. It can also be seen in Table C.1 that the uncertainty for each parameter in each band is far greater than the variations between seed. However, there is some variation between seeds so the model is not completely insensitive to it.

We find that the model provides a good description of the data for each band. The closure phase is reproduced with an error weighted mean residual of < 8° in every band. For V², we have an error weighted mean residual of 1.5 × 10⁻³ (7%) in the K band, 1.5 × 10⁻³ (23%) in the L band, and 4.5 × 10⁻³ (16%) in the M band. However, the N band is the hardest to reproduce, possibly due to the silicate feature, with an error weighted mean residual of 2.9 × 10⁻³ (37%) which could suggest that our description of absorption is incomplete. For example, we do not include the possibility of a foreground component that is reported in Gámez Rosas et al. (2022). Alternatively, or in addition, it could be caused by the base temperatures. Because the N band had a much more tightly defined base temperature for both the disk and the wind, the results suggest that the N band is far more sensitive to them. This is expected because a 300 K would peak in the N band. We chose a flat base temperature with no radial dependence; however, it is possible that this is an overly simplistic assumption and a radial dependency is required to prevent the N band becoming too over resolved when the silicate feature obscures the central region. Further modeling would be needed to see if this is due to the dust composition, geometry, temperature laws, or a combination of these properties.

5.2. Polychromatic fit

The polychromatic fit was initiated from the best starting parameters selected from the band-by-band fit (see Sect. 4.2). We do not enforce edge-on or face-on because only one inc can be found for all bands in a polychromatic fit removing the possibility of the type changing with band. The resulting parameters of the fit can be found in Table C.1 and Fig. C.3. We find that the polychromatic model can aptly explain the majority of the observable data at every wavelength simultaneously with an edge-on inclination. The V² (see Fig. C.4 for best result and Fig. C.5 for median result with errors) is well reproduced in the K and the L band with only minor discrepancy in the 100 m peak which is within 2σ of the median result. The same is true for the M band except for a significant discrepancy at the shortest baseline. The shortest baseline V² in the M band is lower than the next shortest baseline V² at the same PA, this is not seen at any other covered wavelength. Furthermore, in Extended Data Fig. 4 of Gámez Rosas et al. (2022), the V² of every baseline increases or is constant with wavelength except for the shortest baseline. Therefore, we suggest that the shortest baseline has been suppressed from correlated flux loss. This could explain why the shortest baseline has a much higher scatter in the L band than the other baselines. Finally, the model well reproduces the N band at 8.7 μm and 12.3 μm but it is too flat in the silicate feature; namely, the model is too extended when the central region is heavily obscured. This reinforces the idea that the current obscuration or base temperature description is not sufficient to completely explain the N band. To check if a different dust composition is responsible, we additionally tried using ISM dust. However, we find that the quality of the fit in the N band is worse when using ISM dust.

The overall closure phase (see Fig. C.6 for best result and Fig. C.7 for median result with errors) is well reproduced by the best fit model. The closure phase is more sensitive to the small detail in the model, namely, it is more sensitive the distribution of clumps. Therefore, the errors on the median fit are larger than for the V². In the K band, the larger triangles are so sensitive to the absolute position of the clumps that the closure phase is close to uniformly random with different clump distribution for the same geometric parameters. A different seed or a small change in a parameter triggers a different random draw of clump locations which returns a closure phase between ±180° with near uniform probability for the largest triangle. Hence, the median of these triangles are ∼0° ±100°. However, the triangles with a perimeter of 200 m or less are well reproduced. We also find that the image in the K band closely resembles the image reconstruction of GRAVITY Collaboration (2020b) (see Figs. 5 and 6). We determine that, in this case, the K-band geometry is difficult to interpret alone because of the inherent sensitivity to individual clump locations.

Fig. 5. Image of the polychromatic model per wavelength. The white dot is the centre of the model, that is, the location of the SMBH in an aligned system. All images are color-scaled by a power of 0.6 to highlight faint structure.

Figure 6 shows the polychromatic fit evaluated at the wavelengths of the image reconstructions of Gámez Rosas et al. (2022), GRAVITY Collaboration (2020b). The model at every wavelength can be seen to be visually similar to the image reconstruction. While the quality of the reproduction of the observables is the indicator of the goodness of the model, the high level of visual coincidence between the model and image reconstructions shows that the model can reproduce the observed brightness distributions at every wavelength simultaneously. No additional distinct, complex, or unique structures are required to explain the visual differences between bands.

Fig. 6. Image of the polychromatic model unconvolved (top row) and convolved with elongated beam (middle row) at the wavelengths of the image reconstructions of Gámez Rosas et al. (2022), GRAVITY Collaboration (2020b) and the image reconstructions from these works (bottom row). All images are color-scaled by a power of 0.6 to highlight faint structure.

One of the questions we aim to tackle is if the assumed alignment of the wavebands in Gámez Rosas et al. (2022) is accurate in our model scenario. Their aperture photometry, on which their thermal map of the dusty structures was based, was done with pivotal assumptions about the inter-band image alignments. When we compare the model centre of our best-fitting model to the brightest pixel over all wavelengths, we find that the photocentre does move (see Fig. 7). Relative to the model centre, we find that the brightest spot moves up to 1 mas over most wavelengths which is within their aperture. The exception is the silicate feature which moves up to ∼5 mas; however, we do not claim this position with certainty because the fit within the feature is not reliable. Therefore, the assumed alignment between bands in their work is consistent with our model.

Fig. 7. Location of the brightest spot from the polychromatic model with wavelength, (0,0) is the model centre.

6. Physical interpretation of results and comparisons with the literature

Earlier in this paper, we present the modelling results and shown that our model can provide a good reproduction of the observed data. Below, we further discuss the physical implications of the model and describe the comparison to different structures previously observed within the AGN.

6.1. Molecular gas comparison

One of the most direct comparisons which can be made is to the pc scale molecular gas. The dusty structure we observe is mostly gas by mass and such gas can be observed on the pc scales with the Atacama Large Millimeter Array (ALMA). García-Burillo et al. (2019) observe CO(3–2), CO(2–1), and HCO⁺(4–3) in NGC 1068 with an angular resolution of ∼50 − 100 mas. This overlaps the VLTI covered scales in the N band and is well within the spatial extent of our model. The extent of the line emitting structure, with a major-axis FWHM of ∼100 − 250 mas depending on the emission line, is also well within our model outer limit (r_out ≃ 888 mas for the polychromatic model) and the within the extent of our best fitting dust disk (FWHM of ∼300 mas) and can therefore be compared. They find that each molecular gas emission line is elongated with an average PA of 113° and, assuming the gas is in a disk, find an upper limit on the inclination of 30° −40°. The result of the molecular emission line analysis is in complete agreement with the disk component of the polychromatic model with an inclination of 16° ±1° and a PA of ${113}_{2^{°}}^2$ (PA of ${23}_{2^{°}}^2+{90}^{°}$). These gas lines predominantly tracing the disk component and not the wind is consistent with expectations (see Hönig 2019) as the gas lines trace denser regions and the molecular gas is mostly destroyed by radiation in the wind.

Gallimore et al. (2016) observe CO(6–5) with ALMA and also find a disk at 112° consistent with the other ALMA observations and our work. Furthermore, they do succeed in observing this gas line in an outflowing structure. This structure is on mas spatial scales with a PA of 33°. This is slightly offset by 10° with our model but it is in the same direction. An 10° offset is also on the same order as the combined uncertainty of our model and the estimated uncertainties in the CO(6–5) outflow direction from the beam size.

6.2. ALMA observed dust

ALMA has also observed dust continuum in NGC 1068 (García-Burillo et al. 2016). The observation had an angular resolution of 70 × 50 mas at 432 μm meaning the spatial scale overlaps with our N-band observations. The dominant dust temperature probed by ALMA should be cooler and, therefore, disk dominated at these scales. They find the dust has a PA of 142° ±23° which by itself is consistent with our disk. Furthermore, direct comparison between our model at longer wavelengths and their Figs. 1b and c contain additional possible aligning features. The southwest extension on larger scales in their Fig. 1b is in the same direction as our wind. Their Fig. 1c, which shows a slightly higher resolution image, shows an extension much closer to our model disk with a Northern spur that lies along one edge of the outflowing wind. We conclude that our result is completely consistent with the cooler ALMA observed dust.

6.3. Radio and H₂O maser emission

Both GRAVITY Collaboration (2020b) and Gámez Rosas et al. (2022) have made comparisons to the H₂O maser emission. If we assume that the 5 GHz position of S1 in Gallimore & Impellizzeri (2023) aligns with the center of our polychromatic model (i.e., the “central engine” position), we find that the maser disk (their R4+G1 maser group) traces the rim of the obscuring disk and the other groups, which they attribute to spiral arms, trace the outer wall of the dusty outflow. The blueshifted groups coincide with the other side of the outflow to the redshifted. We plot the maser emission against the unconvolved polychromatic model at 2 μm in Fig. C.8. This is an intriguing coincidence, but if the groups coincidental with the wind are associated, then we might expect the outflowing cone inclined toward the observer to coincide with the blueshifted emission whereas it is instead redshifted. We also find that with this alignment, the 22 GHz continuum emission that matches the brightness distribution in a similar manner to Gámez Rosas et al. (2022). Aligning by the dynamic center of the maser disk instead of the 5 GHz peak produces similar alignments. Compared to GRAVITY Collaboration (2020b), Gámez Rosas et al. (2022), and Gallimore & Impellizzeri (2023), our alignment is close to model 4 of GRAVITY Collaboration (2020b) and that of Gámez Rosas et al. (2022). It also agrees with the scenario presented in Sect. 5 of Gallimore & Impellizzeri (2023), where the maser disk (R4+G1) aligns with the obscuring disk.

6.4. Narrow line region

In the unified scheme of AGNs, the narrow line region (NLR) outflows along the polar direction, similar to the model dusty wind. Whether the NLR is coupled to the dust outflow or not, in the simplest case of an aligned system they should both outflow in the same direction. However, it should be noted that the NLR is currently only studied on much larger scales than we probe in this work making comparisons less precise.

Using a biconical outflow model, Das et al. (2006) find a NLR PA of 30°, inclination of 5°, and an outer opening angle of the wind of 80°. However, they probe ∼4″ angular scales which is much larger than our FoV. Poncelet et al. (2008) find that on smaller ∼1″ scales that the NLR is consistent with a bicone of PA = 10°, an inclination towards the observer of ∼10°, and an opening angle of ∼82° for the cone. Given the differences between the works, the difference between the spatial scales probed, and the inherent complexity of NGC 1068 we conclude that our wind structure, with a $\text{PA}={23}_{2^{°}}^2$; an inclination of 16° ±1°; and an opening angle of 90° ±8°, is consistent with the NLR.

6.5. Extended hot dust and scattered light

In Sect. 3.7.3, we discussed the extended ∼100 mas structure observed by Weigelt et al. (2004) in relation to the hot over-resolved component seen in our model and by GRAVITY Collaboration (2020b). However, Weigelt et al. (2004) produced K- and H-band image reconstructions with sub-hundred mas resolution. Our main K-band emitting structure would be unresolved in their images; however, should the hot over-resolved component indeed be scattered light then the scattering medium would be the outflowing dust and the geometry should be related to our wind structure. The extended structure of their images, particularly in the H band, show striking resemblance to our wind structure. Each branch of the X-shape in their Fig. 3 H-band image are closely aligned with the edges of our wind. Their K-band image shows a narrower, less X-like, elongation which aligns with our wind. Additionally, the wind and disk well align with the polarized light H- and K-band structures seen in Gratadour et al. (2015). This coincidental alignment could indicate that the extended K- and H-band structure, and therefore our over-resolved hot component, is indeed from scattered accretion disk light. However, a more in-depth study of the hundred mas scales is needed to make firm conclusions about the origin of this emission.

6.6. Comparison summary

From these comparisons and our model, we can place the VLTI observed structure in the bigger picture. We conclude that the data of NGC 1068 is consistent with the AGN containing a dusty obscuring molecular gas disk that obscures the dust sublimation region and the hotter part of the disk. From this disk, dust is ejected into a wind coincidental with the outflowing NLR and dominates the K-, L-, M-, and N-band emission in the central pcs but between 50%−70% of the K-band flux on hundred mas scales is not from AGN heated thermally emitting dust. The wind dust may then act as a scattering surface which scatters the accretion disk emission contributing to at least part of the missing K-band emission. Finally, the less obscured cooler dust further out in the disk then starts to become more significant in emission in the sub-mm regime seen by ALMA.

7. Comparison of models, improvements, and future prospects

In this paper, we present the results of the polychromatic modeling of NGC 1068. There are some notable features which are worth further consideration, which we discuss here.

7.1. Comparisons between models

Initially, we expected that the band-by-band modeling would be a better reproduction of the data than the polychromatic modeling because it would allow more complex chromatic changes between bands than our model can describe to be approximated by changes in geometry. When we compare the band-by-band results to the polychromatic results, we find that they reproduce the data comparably per wavelength. In the N band, the polychromatic model actually provides a better visibility reproduction and a more visually similar brightness distribution to the image reconstructions (see Fig. 6) which we attribute to the wider wavelength range providing better constraints on the inner geometry. Because the one polychromatic model uses fewer free parameters than three band-by-band models to explain the same data with a similar quality, we can conclude that the polychromatic model is a better general description. This is an encouraging sign that our model is not missing any significant components. Perhaps more importantly, it validates the registration of emission in each band to common components.

When we compare the recovered parameters between the two modeling methods we find good agreement in most cases. We use the band-by-band as prior information for the polychromatic fit and, therefore, we would expect the results to remain similar unless a wider wavelength range provides information not available within a single band. The only parameters we find contain a significant discrepancy are γ_off and f_w, the break radius of the disk temperature power law and the fraction of dust clumps in the wind component, respectively. We adjusted γ_off in the final model to determine the effect on the observables. We find that the parameter is poorly constrained in the edge-on case and the offset is in reality insignificant. We suggest that this is because the disk is mostly self obscured and, therefore, the inner break location is not seen. We find that the discrepancy in f_w is significant. Using N_p to convert f_w into an absolute number of points in the wind, we still find a discrepancy. This means the polychromatic model contains less dust in the wind. We attribute this difference to the additional over-resolved component. The polychromatic model prefers a more central distribution of points in the wind (higher α_w) which results in fewer points outside the image and fewer points needed to construct similar central regions. The fewer points outside the image leads to less over-resolved flux in the model which is compensated by the additional over-resolved parameter. This is additionally evidenced by the fact that V₀ does not increase despite the addition of BB_hot. In the band-by-band fits, the achromatic over-resolved component V₀ was compensating for both the achromatic and chromatic over-resolved structure.

The necessity of V₀ in the polychromatic fit is also interesting. This suggests that there is an over-resolved structure responsible for 32 ± 6% of the flux at every wavelength. However, we see in the band-by-band fit that V₀ seemingly peaks in the L band. It is possible that instead of an achromatic component, we have cool dust component of a few hundred K causing more over-resolved flux in the M and N bands than the L band and that the hot over-resolved component is brighter but is being partially compensated by V₀. However, this is difficult to conclude at this time. In future work, the addition of the total flux when modeling would provide an additional constraint to such a component. Due to the interactions between galactic material, the jets, and the winds, such a cool component is conceivable (e.g., AGN-heated 475 K dust at 450 mas from the center seen with NACO: Rouan et al. 2004).

7.2. Silicate absorption feature

The polychromatic modeling with the disk+wind scenario produced an answer that can reproduce the V² and Φ_c with a error weighted mean residual of 0.0035 (20%) and 11°, respectively. The largest deviation between the model and the data is within the silicate feature near 9.7 μm. Future iterations of the model should allow for different dust compositions to test if this is caused by composition. In this work, we used the olivine and amorphous carbon dust species used in Gámez Rosas et al. (2022). In future work, using the DUSTEM python package (Compiègne et al. 2011) would allow for greater flexibility in dust absorption mixes.

7.3. Future prospects

In this work, we use our model to describe one AGN; however, we developed the model with the intent to keep it general rather than specific to NGC 1068. Primarily, this was to test if NGC 1068 could be explained by a simple generic AGN geometry but it also allows the model to have other applications. We give some prospects towards what could be done in the future below.

7.3.1. The over-resolved components

Future works will be devoted to investigating the physical nature of the over-resolved components. This can be done in part by including the SED and correlated fluxes, which will provide additional constraints on the absolute flux of the components. The SED is particularly necessary to distinguish if the hot-component is indeed hot or instead non-thermal and if it is thermal what temperature it is. Furthermore, high angular resolution imaging in the IR will be studied which may resolved these structures. They could also potentially be resolved through the shorter baselines provided by VLTI AT observations, which we did not include in this work due to complex mismatched FoV effects, or through in-progress studies with shorter baseline LBTI (Large Binocular Telescope Interferometer Hinz et al. 2002) observations. Additionally, approved aperture masked data from 8 m class telescopes in the IR will provide access to the scales of the over-resolved structures (priv. comm.).

These components are interesting to investigate for NGC 1068. However, it is unclear whether they are required to explain VLTI observed AGNs in general. AGNs are extended objects in the IR and often have over-resolved structure but if the IR is dominantly centrally heated dust emission then the over-resolved structure would inherently be included in the model. If the additional components needed for NGC 1068 are a common feature related to the AGN dust (e.g., by scattering), it could be included in the model as a physical feature of the dust clumps themselves rather than a scaling factor. Our model will need to be applied to other AGNs to determine if this is the case. Furthermore, this structure is over-resolved in the relatively local NGC 1068 but, in more distant AGNs, such structures may become resolved.

7.3.2. Radiative transfer simulations

A line of investigation that we do not fully cover in this work is the physicality of the model. Our model is a “toy” model in the sense that we base our geometry and assumptions on previous work and physics but we do not perform detailed physical simulations such as hydrodynamics or full RT; namely, we calculate dust emission and absorption in a physical manner but we set dust temperatures, densities, and distribution parametrically. In a future work, we will collect our toy model parameters and translate them into an RT simulation. The goal of such a work would be to see if our model geometry and dust densities could reproduce our derived distribution of temperatures and brightness in RT. If so, we may better constrain the dust properties in a more “physically realistic” manner. Such a work would allow us to place better constraints on the realism of the properties we derive and test if our model could be a less computationally expensive way to obtain a geometry and dust composition, which can be followed up with a more expensive RT simulation once it has been constrained.

7.3.3. Larger AGN sample

In future works, we will utilize the model on other MATISSE-GRAVITY observations of Sy1 and Sy2 AGNs with less coverage than NGC 1068. Many AGNs have been observed with MATISSE and GRAVITY that do not have sufficient uv coverage to perform image reconstructions. Our model has proven sufficient to explain the observables of the NGC 1068 “test case” and it produces models visually similar to the image reconstructions. Therefore, we can now apply the model to AGN for which image reconstruction is not possible.

NGC 1068 is very well observed and we utilize that through a simple but still relatively complex model when comparing the number of free parameters to models used on less well covered AGNs (e.g., Leftley et al. 2021; GRAVITY Collaboration 2020a; López-Gonzaga et al. 2016; Burtscher et al. 2013). 25 free non-error parameters, equivalent to a fixed point source and four free elongated 2D Gaussians, is not excessive for NGC 1068 as demonstrated in Gámez Rosas et al. (2022). For AGNs observed with only one or two snapshots with MATISSE or GRAVITY alone, some parameters may need to be fixed to prevent over fitting and degeneracy. However, a lower required uv coverage is one advantage of modeling over image reconstruction and an additional advantage of a chromatic model is the significant increase in observational data in a fit. This leads to better constraints for the same uv coverage, as demonstrated by the improvement between the band-by-band and polychromatic fits. This should allow for all parameters to be free for smaller numbers of MATISSE snapshots than achromatic models. For a large survey, the complexity constraints from the data coverage will need to be investigated per object.

7.3.4. Improvements

While RT is a next step in modeling of NGC 1068, our model can be used on other objects. For this, there are improvements that can be made in the toy model itself. For example, it is clear that the model could benefit from the inclusion of more flexible dust compositions. It has been designed such that different dust species can be included easily by providing their optical cross-sections per gram. As such we can incorporate existing dust libraries, for instance, DUSTEM (Compiègne et al. 2011), to give users greater flexibility.

Furthermore, because dust is described by probability distributions, new geometries can be included as needed and added arbitrarily. We plan to streamline this process for ease of use with the goal of including the model as an option within OIMODELER².

7.3.5. Generalizations beyond AGNs

We focus on AGNs in this work, but our model can be applied to other dusty structures. Many dusty objects in astronomy can be described by a disk or a disk+wind. As such, the model can be applied to objects such as some B[e] stars, young stellar objects, and proto-planetary disks. The list can be further extended once we allow for user-defined geometry.

8. Conclusions

We present our successful attempt to reproduce the IR interferometric observables of NGC 1068 from GRAVITY and MATISSE simultaneously with a polychromatic model capable of covering the full K-band to N-band range. This represents a breakthrough in very broad band modeling of dusty structures at very high angular resolution.

We find that a simple clumpy disk and hyperbolic wind composed of thermally emitting non-ISM dust can reproduce the squared visibility (V²) and closure phase (Φ_c) with an error weighted mean residual of 0.0035 (20%) and 11°, respectively, with only one additional ad hoc structure. The resulting PA, inclination, and wind opening angle are consistent with those derived from observations of molecular gas (García-Burillo et al. 2019), ALMA observed dust (García-Burillo et al. 2016), polarized IR light (Gratadour et al. 2015), IR speckle interferometry (Weigelt et al. 2004), radio and maser emission (Gallimore & Impellizzeri 2023), and the NLR (Poncelet et al. 2008; Das et al. 2006) which strongly supports that the IR emission can be explained by the disk+wind iteration of the unification scheme. The additional structure is an over-resolved hot or non-thermal component (≫2000 K, approximated as a 20 000 K black body), for instance, a central sub-arcsecond star forming region or scattered light from the central accretion disk. Furthermore, our obscuration, from both dust clumps and a smooth dust disk, reproduces the data well but over-obscures within the N-band silicate feature; in future work we should investigate if a change in geometry or composition could improve this. One possibility could be using a puffed up inner rim in the disk or changing the dust composition.

When we model the K and L bands individually, we find it difficult to distinguish between two cases for the recovered system: an inclined ring (Sy1 like) or an edge-on wind with an obscured central region (Sy2-like). When fitting all wavelengths simultaneously, we find the model strongly favours a Sy2-like dust distribution where the central region is obscured, this supports the conclusions of Gámez Rosas et al. (2022) and demonstrates the additional constraining power of chromatic information when modeling IR interferometric observations. The model also yields an alignment between K, L, M, and N bands, allowing us to link the structures in the image reconstructions of GRAVITY Collaboration (2020b) and Gámez Rosas et al. (2022). Our model interprets the bright regions of each band as the wind whilst the sublimation region remains fully obscured at all observed wavelengths, this supports the unified theory of AGN presented by Antonucci (1993) and the disk+wind interpretation presented in Hönig (2019). We conclude that the disk+wind unified model of AGN is a good description of the inner structure of NGC 1068 in the IR, which can be translated into a full radiative transfer model to further investigate the dust properties.

In this work, we have demonstrated the efficacy of simple chromatic modeling as a powerful tool for IR interferometric data. In future works, we will extend our geometry for NGC 1068 to a full RT simulation and apply our versatile modeling approach to a larger sample of AGNs. This will further test our approach and construct a sample of AGNs for which we have a good description of their dust structure. However, our technique is not limited to AGNs. In the future, we aim to provide our model as a general tool for interferometric studies of dusty structure, such as B[e] stars, young stellar objects, and proto-planetary disks.

Acknowledgments

We would like to thank the referee for the kind comments and suggestions which helped improve this work. This work was supported by French government through the National Research Agency (ANR) with funding grant ANR AGN_MELBa (ANR-21-CE31-0011) and the Université Côte d’Azur JEDI investment in the Future project managed by the under the reference number ANR-15-IDEX-01. 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, GALARIO, and MATPLOTLIB.

References

Appendix A: Band-by-band results – Additional information

A.1. Initial band-by-band fit

Here we present the initial results of the band-by-band modeling for the Sy2 case in Sect. A.1.1 and Sy1 case in Sect. A.1.2. The initial fits contain some interesting and useful results beyond providing a starting point for the second run and we present them here. It should be noted that we present these results because their trends and scatter inform us about structure we might otherwise overlook, the dependency of the model on the seed, and how strongly defined each initial parameter is for the final fit. These results should not be used for further analysis; instead, the results in Sect. 5.1 should be used. The initial run is performed to determine a reliable starting position for the final run as a method to fully sample the parameter space in a more time efficient manner. It is not run for a sufficient length of time to determine reliable errors from a uniform start for our case.

A.1.1. Sy2 case

The initial run, without a well-defined starting position, produced results that were very consistent between seeds in the L, M, and N bands suggesting that the underlying geometry can be found in these bands for any seed. The K band produced a wider array of results, which suggests that it could be more sensitive to the individual clump positions and, therefore, the random seed, than the longer wavelengths. The best model from the initial fit for each seed and band can be seen in Fig. A.1. The parameters for the initial fits for each seed can be found in Fig. C.3. We compare the data from each best fit model to the data in Fig. C.9.

Fig. A.1. Best model for each band (row) and each seed (column) for the initial Sy2 band-by-band fit as determined by the maximum likelihood. The best seed for each band is marked with a red x. The K, L, M, N bands are evaluated at 2.1 μm,3.5 μm,4.7 μm, and 8.5 μm, respectively. All images are normalized and given a 0.6 power color scaling to match that of Gámez Rosas et al. (2022) for easy comparison.

Most parameters were consistent between bands within the scatter of the different seeds; however, there are some exceptions. Of particular interest are a_w, b, ang and inc. inc, the inclination angle, remains consistent within the inter-seed scatter between the all bands at an average value of ∼15° which is close to the NLR geometry inferred inclination of 10° Poncelet et al. (2008); ang, the PA of the cone or the system/polar axis, is consistent within inter-seed scatter in the L, M, and N bands at 24° ±5° which agrees with the system axis from Asmus et al. (2016) of 28° ±20° measured from the NLR, radio, masers, and mid-IR imaging. Furthermore, ang in the K band is highly scattered but consistently higher. Although ang is close between the L, M, and N bands, there is a tentative positive relation with wavelength that could be indicative of a chromatic anisotropy in the brightness distribution of the wind. For example, anisotropic illumination of the wind by the accretion disk causes one side to be hotter than the other which would be more apparent at shorter wavelengths and smaller spatial scales. Our model does not allow for anisotropic illumination and, therefore, this effect would be manifested as a rotation. Then, Parameter b tentatively increases with wavelength whereas a_w decreases more strongly. Also, b and a_w are both parameters that define the hyperbolic wind shape and a hyperbolic line of $\frac{z^2}{a_{\mathrm{w}}^2}-\frac{x^2}{b^2}=1$ approaches $z=\frac{a_{\mathrm{w}}}{b}x$ at a large x. If we define the opening angle of the wind as the angle of the approached line from the polar axis, we can define it as arctan(b/a_w). Then, the wavelength dependence of b and a_w instead suggests that the opening angle of the wind increases with wavelength from an opening angle in the K band of ∼40° to ∼50° in the L band. We show this in Fig. A.2.

Fig. A.2. Opening angle of the initial band-by-band fit in the Sy2 case (Sect. A.1.1) against wavelength.

Finally, we find that V₀ is consistently lower in the K band, which is expected from our initial assumption we would require an over-resolved hot component that is not included in the band-by-band fits. It is consistent with the over-resolved fraction of 0.5 found in GRAVITY Collaboration (2020b) when performing the image reconstruction in the K band.

The presence of the silicate feature makes the N-band observables difficult to reproduce without a wider wavelength coverage. We suspect that this is because much of the N band is heavily obscured and therefore much of the central geometry is not available; conversely, the wider FoV makes more of the extended geometry available. This results in the different seeds providing similar ang and geometry for the winds but the disk temperature distribution is unconstrained.

Compared to the other bands, n_d is high and better constrained in the N band. This is expected because of the presence of the silicate feature which is a strong constraint on the obscuration. However, n_c is no better constrained or higher in the N band than the other bands. This suggests that the smooth obscuration is relatively more impactful at longer wavelengths, possibly due to the lower resolution causing reduced distinction between smooth and clumpy structure. When translating the recovered n_d and n_c into τ along an equatorial LOS towards the model’s center, we find n_c is the dominant contribution to the disk obscuration due to the number of clumps in the LOS. However, n_d is still significant for some LOS. Particularly if the LOS is through a region with fewer clumps.

The initial L-, M-, and N-band images produced by the models (see Fig. A.1, the N band is evaluated at 8.5 μm) are all quite consistent between seeds and, more importantly, with the image reconstructions of Gámez Rosas et al. (2022). The resemblance between the two methods already confirms that the brightness distributions found in Gámez Rosas et al. (2022) can be reproduced by a disk+wind geometry.

At first glance, the K band is seemingly less comparable to the image reconstruction of GRAVITY Collaboration (2020b); although, the image produce by seed S4 is visually similar with a ring-like structure and surrounding clumps. We find that the K-band V² are dominated by a north-east to south-west binary-like structure with closure phase information introduced by a significantly clumpy structure. This is what is seen in the image of GRAVITY Collaboration (2020b), there is a bright complex north-east component, the ring-like structure, and a dim southern clump. We suggest that because the true brightness distribution is seemingly clumpy and that for a clumpy structure, due to finer angular resolution, the relative positions of the clumps would more important in the K band than in longer wavelengths, it is unlikely that our model, particularly in the initial fit, would find a perfect solution for a particular seed. This is supported by the higher variety of images and parameters in the K band indicating a higher sensitivity to the seed. Therefore, the more likely solution for the model to discover is dominated by the easy to reproduce binary-like structure which is a good description of the V² and the Φ_c at smaller triangle perimeters (< 220m). As such, we often find a bright dense cluster of clumps and a faint cluster, a solution that can be easily reproduced with any seed. This bright-dim spot structure can also be reproduced with a ring similar to the one in Vermot et al. (2021).

A.1.2. Sy1 case

The initial run for the Sy1 case starts from the same parameters as the Sy2 case, except that ϕ must be less than inc. We note that this initial run will not be likely to reproduce case of Vermot et al. (2021) because it would require an ang of 180° offset from the initial position. This Sy1 start is specifically to test the model 3 and model 4 of GRAVITY Collaboration (2020b). The best fits for each seed compared to the data can be seen in Fig. C.10 and the images are shown in Fig. C.11.

The Sy1 case, with the initial parameters given, converged to Sy2-like results. The model preferred inclinations that were as close to obscured as allowed or positioned such that the wind acted as the obscurer for the central region. The best fits produced were the same order of magnitude as the Sy2 case when comparing likelihood; however, they were consistently of lower likelihood.

A.1.3. Rotated Sy1 case

To test the ring model of Vermot et al. (2021), we perform the same initial modeling as the Sy1 case but offset the start of ang by 180°. It is important to note that Vermot et al. (2021) explains the ring as a distinct structure, such as a tidally disrupted clump, separate from the unification scheme. Our model can produce a ring that produces a similar brightness distribution but this ring will be the sublimation region, not a distinct structure.

For this case, we recover the ring structure from Vermot et al. (2021) in the K band and in the L band. In the M Band a ring is still recovered but closer to edge-on and in the N band, we recover the Sy2 case. We conclude that the ring fit is a valid solution in the K and L bands individually (see Fig. C.12). However, when we perform a polychromatic fit of all bands begun from the L-band ring fit, we find that it is not valid for when simultaneously explaining all wavelengths with our geometry. The polychromatic ring fit does not reproduce the observables as well as the Sy2 model, particularly in the N band, as can be seen in Figs. C.13 and C.14.

Appendix B: Technical check – Dust mass

In our model, we calculate the dust obscuration from clumps and smooth structure. This means we can translate this into a dust density and mass. While this will be inaccurate without the SED, we can check that we do not use unphysical amounts of dust for our best fitting model. In the band-by-band fits, n_c varies between 3 − 4; however, within a band the value is significantly more constrained. This is most likely due to the dust composition being poorly constrained by a single band and therefore n_c becomes semi-degenerate with p_g. In Fig. 3, we show the cross-section per gram of the two materials in our dust mix, it can be seen that outside the silicate absorption feature the two materials have differences of up to two orders of magnitude. For example, in the K band it is possible to achieve the same absorption with pure amorphous carbon and pure olivine dust if there is ∼45 times more dust by mass, equivalent to a change in n_c of 1.65. Ergo, for better constraints on a clumps mass we use the result of $n_{\mathrm{c}}=-4_{0.2}^{0.2}$ log₁₀(g cm⁻² mas⁻¹) and $r_{\mathrm{c}}=1.0_{0.1}^{0.1}$ mas from the polychromatic model. Using a distance to NGC 1068 of 14.4 Mpc from the NASA/IPAC Extragalactic Database³ (NED), we can use a conversion from mas to cm of 2.15×10¹⁷ cm mas⁻¹. This leads to the volume of a clump being 42 × 10⁵¹ cm³ which gives a mass per clump of ∼0.01 M_⊙. This is equivalent to a total dust mass in the disk of ${170}_{60}^{100}$ M_⊙, including the uncertainties from N_p and f_w. García-Burillo et al. (2019) measure a gas mass of the molecular torus of ∼1.5 × 10⁵ M_⊙ which would make the disk ∼0.1% dust. For comparison, this is below ISM values and, therefore, our model is not requiring an unphysical amount of dust.

Appendix C: Figures and tables

Fig. C.1. V2 and Φc from the best models in Sect. 5.1 at one wavelength per band for each seed compared to the observed data.

Fig. C.2. Best model for each band (row) and each seed (column) for the final band-by-band fit as determined by the maximum likelihood. The best seed for each band is marked with a red x at the top left corner. The K, L, M, and N bands are evaluated at 2.1 μm, 3.5 μm, 4.7 μm, and 8.5 μm, respectively. All images are normalized and given a 0.6 power color-scaling to match that of Gámez Rosas et al. (2022) for easy comparison.

Fig. C.3. Parameters of the best model for each band and each seed (represented by colors as given in legend) for the initial band-by-band fit (colored circles) and final band-by-band fit (colored crosses) as determined by the maximum likelihood. The seed used as the starting position for the final fit in each band is given a black border. Overplotted as lines are the best and median fit for the polychromatic model. The shaded region is the 1σ errors on the median fit parameters.

Fig. C.4. V2 of the polychromatic model. The colors represent the PA on sky of the telescope pair used for that point. The color scheme is continuous and cyclical about π, due to the fact that a V2 measurement at PA =x and x + π are equivalent. The line is the best polychromatic model by maximum likelihood evaluated at different PA in $\frac{\pi }{6}$ intervals.

Fig. C.5. V2 of the polychromatic model. The colors represent the PA on sky of the telescope pair used for that point. The color scheme is continuous and cyclical about π, due to the fact that a V2 measurement at PA =x and x + π are equivalent. The line is the median fit of the polychromatic model with 1 σ errors as the shaded region evaluated at different PA in $\frac{\pi }{6}$ intervals.

Fig. C.6. Closure phase of the polychromatic model per wavelength compared to observations. The closure phases for the best fitting model by maximum likelihood were calculated for each of the four telescope triplets, at all accessible uv positions for the UTs where NGC 1068 has a minimum altitude of 35°. The average vector PA is the average of the PA for each of the three baselines in the triplet, ensuring the telescopes are in ascending numerical order. For instance, the average PA for U1-U2-U3 is the is average of the PA of the baselines U1-U2, U2-U3, and U1-U3 (not U3-U1). For visualization purposes, we averaged the different BCD positions for the MATISSE data and reordered the GRAVITY closure phase telescope triplets to match the order of the captions. The changes in closure phase sign from reordering were taken into account.

Fig. C.7. Closure phase of the median polychromatic model per wavelength with 1 σ errors compared to observations. The closure phases were calculated for each of the four telescope triplets at all accessible uv positions for the UTs where NGC 1068 has a minimum altitude of 35°. The average vector PA is the average of the PA for each of the three baselines in the triplet ensuring the telescopes are in ascending numerical order. For instance, the average PA for U1-U2-U3 is the is average of the PA of the baselines U1-U2, U2-U3, and U1-U3 (not U3-U1). For visualization purposes, we have averaged the different BCD positions for the MATISSE data and reordered the GRAVITY closure phase telescope triplets to match the order of the captions. The changes in closure phase sign from reordering were taken into account.

Fig. C.8. Comparison of the H2O maser emission from Gallimore & Impellizzeri (2023) to the unconvolved polychromatic model at 2 μm. The model has a log color scaling to highlight faint structure. The model center is aligned with the 5 GHz position of S1 from their work. The velocity is their recession velocity.

Fig. C.9. V2 and Φc from the best models from Sect. A.1.1 at one wavelength per band for each seed compared to the observed data.

Fig. C.10. V2 and Φc from the best models in Sect. A.1.2 at one wavelength per band for each seed compared to the observed data.

Fig. C.11. Best model for each band (row) and each seed (column) for the initial Sy1 band-by-band fit as determined by the maximum likelihood. The best seed for each band is marked with a red x at the center. The K, L, M, and N bands are evaluated at 2.1 μm,3.5 μm,4.7 μm, and 8.5 μm, respectively. All images are normalized and given a 0.6 power color-scaling to match that of Gámez Rosas et al. (2022) for easy comparison.

Fig. C.12. Best model for the initial Sy1 ring-case band-by-band (top row) and polychromatic (bottom row) fit as determined by the maximum likelihood. All images are normalized and given a 0.6 power colour scaling to match that of Gámez Rosas et al. (2022) for easy comparison.

Fig. C.13. Same as Fig.C.5 for the ring case of Sect.A.1.3. The line is the median fit of the polychromatic model with 1 σ errors as the shaded region evaluated at different PA.

Fig. C.14. Same as Fig.C.7 for the ring case of Sect.A.1.3. The closure phases were calculated for each of the four telescope triplets at all accessible uv positions for the UTs where NGC 1068 has a minimum altitude of 35°. The average vector PA is the average of the PA for each of the three baselines in the triplet ensuring the telescopes are in ascending numerical order. For instance, the average PA for U1-U2-U3 is the is average of the PA of the baselines U1-U2, U2-U3, and U1-U3 (not U3-U1). For visualization purposes, we have averaged the different BCD positions for the MATISSE data and reordered the GRAVITY closure phase telescope triplets to match the order of the captions. The changes in closure phase sign from reordering were taken into account.

All Tables

All Figures

	Fig. 1. Schematic of the disk described in Eq. (1) and wind defined in Eq. (2). We provide the example aw and b used to make the wind schematic.
In the text

	Fig. 2. Schematic of the disk described in Eq. (1) and wind defined in Eq. (2) with clumps distributed with the probabilities described in Eqs. (3) and (4). We provide the example β and αw used to make the wind schematic.
In the text

	Fig. 3. Optical cross-section per gram of the dust materials used in this work as well as a 20% amorphous carbon and 80% olivine mix (pg = 0.2).
In the text

	Fig. 4. Visibility of NGC 1068 between 16 Mλ and 20 Mλ compared to wavelength in purple. Overplotted is a best fit model of an over- and under-resolved black-body model in red.
In the text

	Fig. 5. Image of the polychromatic model per wavelength. The white dot is the centre of the model, that is, the location of the SMBH in an aligned system. All images are color-scaled by a power of 0.6 to highlight faint structure.
In the text

	Fig. 6. Image of the polychromatic model unconvolved (top row) and convolved with elongated beam (middle row) at the wavelengths of the image reconstructions of Gámez Rosas et al. (2022), GRAVITY Collaboration (2020b) and the image reconstructions from these works (bottom row). All images are color-scaled by a power of 0.6 to highlight faint structure.
In the text

	Fig. 7. Location of the brightest spot from the polychromatic model with wavelength, (0,0) is the model centre.
In the text

Fig. A.1. Best model for each band (row) and each seed (column) for the initial Sy2 band-by-band fit as determined by the maximum likelihood. The best seed for each band is marked with a red x. The K, L, M, N bands are evaluated at 2.1 μm,3.5 μm,4.7 μm, and 8.5 μm, respectively. All images are normalized and given a 0.6 power color scaling to match that of Gámez Rosas et al. (2022) for easy comparison.

In the text

	Fig. A.2. Opening angle of the initial band-by-band fit in the Sy2 case (Sect. A.1.1) against wavelength.
In the text

	Fig. C.1. V2 and Φc from the best models in Sect. 5.1 at one wavelength per band for each seed compared to the observed data.
In the text

Fig. C.2. Best model for each band (row) and each seed (column) for the final band-by-band fit as determined by the maximum likelihood. The best seed for each band is marked with a red x at the top left corner. The K, L, M, and N bands are evaluated at 2.1 μm, 3.5 μm, 4.7 μm, and 8.5 μm, respectively. All images are normalized and given a 0.6 power color-scaling to match that of Gámez Rosas et al. (2022) for easy comparison.

In the text

Fig. C.3. Parameters of the best model for each band and each seed (represented by colors as given in legend) for the initial band-by-band fit (colored circles) and final band-by-band fit (colored crosses) as determined by the maximum likelihood. The seed used as the starting position for the final fit in each band is given a black border. Overplotted as lines are the best and median fit for the polychromatic model. The shaded region is the 1σ errors on the median fit parameters.

In the text

	Fig. C.4. V2 of the polychromatic model. The colors represent the PA on sky of the telescope pair used for that point. The color scheme is continuous and cyclical about π, due to the fact that a V2 measurement at PA =x and x + π are equivalent. The line is the best polychromatic model by maximum likelihood evaluated at different PA in $\frac{\pi }{6}$ intervals.
In the text

	Fig. C.5. V2 of the polychromatic model. The colors represent the PA on sky of the telescope pair used for that point. The color scheme is continuous and cyclical about π, due to the fact that a V2 measurement at PA =x and x + π are equivalent. The line is the median fit of the polychromatic model with 1 σ errors as the shaded region evaluated at different PA in $\frac{\pi }{6}$ intervals.
In the text

Fig. C.6. Closure phase of the polychromatic model per wavelength compared to observations. The closure phases for the best fitting model by maximum likelihood were calculated for each of the four telescope triplets, at all accessible uv positions for the UTs where NGC 1068 has a minimum altitude of 35°. The average vector PA is the average of the PA for each of the three baselines in the triplet, ensuring the telescopes are in ascending numerical order. For instance, the average PA for U1-U2-U3 is the is average of the PA of the baselines U1-U2, U2-U3, and U1-U3 (not U3-U1). For visualization purposes, we averaged the different BCD positions for the MATISSE data and reordered the GRAVITY closure phase telescope triplets to match the order of the captions. The changes in closure phase sign from reordering were taken into account.

In the text

Fig. C.7. Closure phase of the median polychromatic model per wavelength with 1 σ errors compared to observations. The closure phases were calculated for each of the four telescope triplets at all accessible uv positions for the UTs where NGC 1068 has a minimum altitude of 35°. The average vector PA is the average of the PA for each of the three baselines in the triplet ensuring the telescopes are in ascending numerical order. For instance, the average PA for U1-U2-U3 is the is average of the PA of the baselines U1-U2, U2-U3, and U1-U3 (not U3-U1). For visualization purposes, we have averaged the different BCD positions for the MATISSE data and reordered the GRAVITY closure phase telescope triplets to match the order of the captions. The changes in closure phase sign from reordering were taken into account.

In the text

	Fig. C.8. Comparison of the H2O maser emission from Gallimore & Impellizzeri (2023) to the unconvolved polychromatic model at 2 μm. The model has a log color scaling to highlight faint structure. The model center is aligned with the 5 GHz position of S1 from their work. The velocity is their recession velocity.
In the text

	Fig. C.9. V2 and Φc from the best models from Sect. A.1.1 at one wavelength per band for each seed compared to the observed data.
In the text

	Fig. C.10. V2 and Φc from the best models in Sect. A.1.2 at one wavelength per band for each seed compared to the observed data.
In the text

Fig. C.11. Best model for each band (row) and each seed (column) for the initial Sy1 band-by-band fit as determined by the maximum likelihood. The best seed for each band is marked with a red x at the center. The K, L, M, and N bands are evaluated at 2.1 μm,3.5 μm,4.7 μm, and 8.5 μm, respectively. All images are normalized and given a 0.6 power color-scaling to match that of Gámez Rosas et al. (2022) for easy comparison.

In the text

	Fig. C.12. Best model for the initial Sy1 ring-case band-by-band (top row) and polychromatic (bottom row) fit as determined by the maximum likelihood. All images are normalized and given a 0.6 power colour scaling to match that of Gámez Rosas et al. (2022) for easy comparison.
In the text

	Fig. C.13. Same as Fig.C.5 for the ring case of Sect.A.1.3. The line is the median fit of the polychromatic model with 1 σ errors as the shaded region evaluated at different PA.
In the text

Fig. C.14. Same as Fig.C.7 for the ring case of Sect.A.1.3. The closure phases were calculated for each of the four telescope triplets at all accessible uv positions for the UTs where NGC 1068 has a minimum altitude of 35°. The average vector PA is the average of the PA for each of the three baselines in the triplet ensuring the telescopes are in ascending numerical order. For instance, the average PA for U1-U2-U3 is the is average of the PA of the baselines U1-U2, U2-U3, and U1-U3 (not U3-U1). For visualization purposes, we have averaged the different BCD positions for the MATISSE data and reordered the GRAVITY closure phase telescope triplets to match the order of the captions. The changes in closure phase sign from reordering were taken into account.

In the text