<figure><img src="image_1.png" /></figure>

HAL
open science

# Direct imaging of active galactic nucleus outflows and their origin with the 23 m Large Binocular Telescope

J W Isbell, S. Ertel, J.-U. Pott, G. Weigelt, M. Stalevski, J. Leftley, W. Jaffe, R G Petrov, N. Moszczynski, P. Vermot, et al.

► To cite this version:

J W Isbell, S. Ertel, J.-U. Pott, G. Weigelt, M. Stalevski, et al. Direct imaging of active galactic nucleus outflows and their origin with the 23 m Large Binocular Telescope. Nature Astronomy, 2025, 9 (3), pp.417-427. 10.1038/s41550-024-02461-y. hal-05123198

HAL Id: hal-05123198

https://hal.science/hal-05123198v1

Submitted on 20 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Direct imaging of AGN outflows and their origin with the 23 m Large Binocular Telescope

J. W. ISBELL,¹,² S. ERTEL,¹,³ J.-U. POTT,² G. WEIGELT,⁴ M. STALEVSKI,⁵,⁶ J. LEFTLEY,⁷ W. JAFFE,⁸ R. G. PETROV,⁷

N. MOSZCZYNSKI,⁹ P. VERMOT,⁹ P. HINZ,¹⁰ L. BURTSCHER,¹¹ V. GÁMEZ ROSAS,⁸ A. BECKER,³ J. CARLSON,³

V. FARAMAZ-GORKA,¹ W. F. HOFFMANN,¹ J. LEISENRING,¹ J. POWER,³ AND K. WAGNER¹

¹Department of Astronomy and Steward Observatory, The University of Arizona, 933 North Cherry Ave, Tucson, AZ, 85721, USA

²Max-Planck-Institut für Astronomie, Königstuhl 17, Heidelberg 69117, Germany

³Large Binocular Telescope Observatory, The University of Arizona, 933 North Cherry Ave, Tucson, AZ, 85721, USA

⁴Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

⁵Astronomical Observatory, Volgina 7, 11060 Belgrade, Serbia

⁶Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281-S9, Gent, 9000, Belgium

⁷Laboratoire Lagrange, Université Côte d'Azur, CNRS, Boulevard de l'Observatoire, CS 34229, 06304

Nice Cedex 4, France

⁸Leiden Observatory, Leiden University, Niels Bohrweg 2, NL-2333 CA Leiden, The Netherlands

⁹LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris Cité, 5 place Jules Janssen, 92190

Meudon, France

¹⁰University of California, Santa Cruz, 1156 High St, Santa Cruz, CA 95064, USA

¹¹Umweltinstitut München, Goethestraße 20, 80336 München, Germany

ABSTRACT

Active galactic nuclei (AGNs) are a key component of galaxy evolution due to feedback on the host from its supermassive black hole. The morphology of warm, in- and outflowing dusty material can reveal the nature of the onset of feedback, AGN feeding, and the unified model of AGN. Here we use the Large Binocular Telescope Interferometer (LBTI) to image the dense, obscuring disk and extended dusty outflow region of NGC 1068. In Fizeau imaging mode the LBTI synthesizes the equivalent resolution of a 22.8 m telescope. The 8.7 µm Fizeau images of NGC 1068 have an effective resolution of 47 × 90 mas (3.3 × 6.2 pc) in a 5" field of view after performing PSF deconvolution techniques described here. This is the only extragalactic source to be Fizeau imaged using the LBTI, and the images bridge the scales measured with the Very Large Telescope Interferometer (VLTI; 0.5-5 pc) and those of single telescopes such as JWST and Keck (> 15 pc). The images detect and spatially resolve the low surface brightness mid-infrared (MIR) features in the AGN disk/wind region that are over-resolved by the VLTI. The images show strong correlation between MIR dust emission and near-infrared (NIR) emission of highly excited atomic lines observed by SINFONI. Such LBTI imaging is a precursor to infrared imaging using the upcoming generation of extremely large telescopes, with angular resolutions up to 6x better than JWST, the largest space telescope in orbit.

Main

The study of AGN at high spatial resolution aims to test the Unified Model of AGN (Antonucci 1993), understand the AGN lifecycle, and resolve the onset of AGN feedback – crucial information for studies of AGN at high redshift and within cosmological simulations. Dusty outflows are of particular interest because they dominate the energy budget of AGN feedback, especially in highly accreting AGN (e.g., quasars; Fabian (2012)). The winds are thought to be driven by radiation pressure from high-energy photons in the vicinity of the SMBH (e.g., Wada 2012; Wada et al.
2016; Williamson et al. 2019; Leftley et al. 2019; Williamson et al. 2020). In simulations, anisotropic radiation from the accretion disk (Netzer 1987) drives the dust primarily in the polar direction, but the outflows are further impacted by and impact inflowing material in a so-called fountain flow (Wada 2012). Polar extended dust has been commonly observed on parsec scales in a number of AGN with VLTI/MIDI (Burtscher et al. 2013; López-Gonzaga et al. 2016) and studied more extensively in two nearby AGN with VLTI/MATISSE (Isbell et al. 2022, 2023; Gámez Rosas et al. 2022). However, how the AGN and its winds are fueled and how the polar ejected material interacts with the host galaxy remain open questions.

NGC 1068 is an important laboratory for the understanding of AGN feeding and feedback. Located only 14.4 Mpc away, it is the prototypical Seyfert 2 AGN (Antonucci 1993); it led to the development of the Unified AGN model. Its proximity and brightness have allowed for extremely high spatial resolution imaging of its nucleus in the near-infrared (NIR; GRAVITY Collaboration et al. 2020), MIR (Gámez Rosas et al. 2022), sub-mm (García-Burillo et al. 2016; Impellizzeri et al. 2019; García-Burillo et al. 2019), and radio (Gallimore et al. 2004; Mutie et al. 2024). The sub-mm observations probe cold dust and dense molecular gas on similar scales as the LBTI. Interpretation of the connections between dusty, molecular material, hot ionized outflows, and the AGN itself is complicated by interactions with the “bent” radio jet (May & Steiner 2017), which leaves the core oriented approximately North-South, but is seen at large scales to be oriented northeast-southwest (e.g., Mutie et al. 2024). Additionally, high-resolution NIR/MIR imaging at the VLTI with the 8.1 m Unit Telescopes resolves out structures larger than ≈ 4 pc at 8.7 µm, so the link between the small-scale and galactic structures is tenuous. Measurements on scales intermediate to the existing 8.1 m single-dish observations and the 34 m VLTI shortest unit telescope (UT) baseline are necessary to understand how material is delivered to the accretion disk and to the nascent outflows.

The LBTI (Hinz et al. 2016; Ertel et al. 2020) co-phases and interferometrically combines the beams from the two 8.4 m mirrors of the Large Binocular Telescope (LBT), separated 14.4 m. In Fizeau imaging, this results in a virtual telescope aperture with a maximum baseline of 22.8 m. While the Auxiliary Telescopes (ATs) at the VLTI offer similar baselines, their 1.8 m diameters provide limited sensitivity compared to the LBTI. The LBTI thus provides complementary data to the existing VLTI observations of NGC 1068. Open phase loop Fizeau interferometric imaging is the simplest interferometric mode of the LBTI. It has great potential for the high-fidelity, wide-field (~5 arcsec), intermediate-contrast imaging of faint near-/mid-IR targets or targets that are too extended for phase tracking (examples include Solar System bodies, AGNs, ultra-luminous infrared galaxies, and massive and evolved stars). The mode, however, has not been actively developed since its successful observations of Io (Leisenring et al. 2014; Conrad et al. 2015; de Kleer et al. 2017). Open loop Fizeau imaging at the LBTI is poised to bridge the gap between high-dynamic-range single-dish observations with e.g., JWST, and low-dynamic-range high-resolution interferometric images with e.g., VLTI and CHARA. Its 22.8 m resolution is comparable to the upcoming generation of 30 m class telescopes, and can serve as a key testing ground for science cases and techniques.

In this paper, we use open loop Fizeau imaging with the LBTI on an extragalactic target, providing spatial constraints on the structure and heating of the nuclear outflow and dust torus in NGC 1068. These observations reveal low surface brightness features at unprobed scales intermediate to single-dish telescopes and long baseline infrared interferometry.

# Observations and Data Reduction

On UTC 6 November 2022 we observed NGC 1068, 1 Tau, and Procyon using the NOMIC instrument (Hoffmann et al. 2014) on the LBTI in the Open-Loop Fizeau mode. We used the W08699-9.122 filter ($\lambda_c$ = 8.7 µm, range [8.13, 9.35] µm) and took 20 nodding cycles on the target and calibrators, obtaining 59° field rotation. Before and after the science target we observed 1 Tau and Procyon to serve as flux and point spread function (PSF) calibrators. We later selected images via “lucky fringing,” resulting in an exposure time of 200 × 0.0137 s = 2.74 s per cycle or 54.8 s total exposure on NGC 1068. The resulting corotated and stacked target images and PSF calibrator images are each shown in Fig. 1. See Methods §1 for details on frame selection, PSF calibration, and flux calibration.

We deconvolved the calibrator PSF from the science target PSF in order to recover the underlying NGC 1068 source flux distribution. While this had previously been done for Io (Leisenring et al. 2014) with the LMIRCam detector, NGC 1068 is much fainter and extended, and several changes were made to the AO system of the LBT (Pinna et al. 2016) in the years following the publication of those data. We therefore performed image deconvolution in three different ways (see Methods §2), finding good agreement between them. We show the results of deconvolution using the Richardson-
Lucy (R-L) method as in (Leisenring et al. 2014) in Fig. 2. Other deconvolutions are shown in Extended Data (ED). We estimate the resulting resolution using the southernmost and northernmost point sources in the image, measuring FWHM 47 × 90 mas after deconvolution. In this article, we focus on the R-L deconvolved images of NGC 1068 because of the robust legacy and fidelity of the method.

# The Dusty Wind

The circumnuclear region in NGC 1068 has been extensively studied in high resolution at numerous wavelengths. We compare the recovered 8.7 µm morphology at 47 mas (3.3 pc) resolution to existing radio maps (at 5 GHz with 105×51 mas resolution; Gallimore et al. (2004)), optical images ([OIII5007] emission with 52 mas resolution; Macchetto et al. (1994)), and NIR/MIR observations using single-dish deconvolution (~100 mas resolution; Bock et al. (2000)) and interferometry (Weigelt et al. 2004; Gámez Rosas et al. 2022, 74 mas and 3 mas resolution, respectively). We show this comparison in Fig. 3, specifically showing the 5 GHz radio emission, previous single-dish MIR imaging, [OIII] emission, and VLTI/MATISSE MIR reconstructed image with respect to our recovered MIR structures.

In the deconvolved images (Fig. 2 and ED Figs. 1, 2, & 3), we see significant extended flux along the north-south direction as well as a bright point sources to the northeast at 600 mas (45 pc) and to the south at 300 mas (20 pc, called the “SN” for southern *N*-band source). The 47 × 90 mas recovered resolution allows us to distinguish substructures near the AGN core. Here we see extensions of flux from the flux peak toward the north and northeast out to ≳ 10 pc (150 mas). We refer to this northeast extension hereafter as “NEX” and the northern extension as “NEP.” Between the flux peak and the SN there is a drop in emission we refer to as the “dark band.”

We find morphological agreement with both the radio emission and previous optical-infrared imaging. The authors in Gallimore et al. (2004) locate the SMBH in a structure called S1; similar to positions inferred from maser emission (Gallimore & Impellizzeri 2023), model fitting to MATISSE data(Leftley et al. 2024), and morphological correspondence(Bock et al. 2000; Gámez Rosas et al. 2022) we place the peak of the MIR emission 50 mas north of S1. We then find that the southernmost and northern 8.7 µm structures correspond roughly to peaks in the radio emission (called ‘S2’ and ‘C’ in Gallimore et al. (2004), respectively). However, we note in both cases there is an offset of ≈ 50 mas, with the MIR emission found exterior to the radio. The northernmost emission is near the ‘NE’ source of (Gallimore et al. 2004), but like in (Bock et al. 2000) we find a ~100 mas separation between the MIR and 5 GHz emission. The 7.9 µm emission of (Bock et al. 2000) includes a peak far to the NE that we do not recover. Assuming the flux of the northeastern component has remained the same since the (Bock et al. 2000) observations, we would expect a flux density of ≈ 0.68 mJy/px, which is above the 2σ flux level of the deconvolved image. Over the same period, the 5 GHz flux of cloud ‘C’ has been shown to decrease by ≈ 50% (Mutie et al. 2024); a similar flux change would put this source below our current detection limit. Indeed, when comparing our 144 mas aperture flux of cloud ‘C’ to the 200 mas aperture flux of (Bock et al. 2000), we measure approximately one sixth of the previous value (100 ± 6 mJy at 8.7 µm versus 661 mJy at 7.9 µm).

We see that the northern structures have counterparts in the [OIII] and NIR emission. Notably, our measured NEX feature traces the opening angle of the ionization cone (Macchetto et al. 1994). SINFONI observations of several atomic emission lines ([Fe II], [Si IV], Pa α, and HeI) and molecular hydrogen show a similar morphology to the arc of flux starting from ‘C’ and continuing from the tongue counterclockwise toward ‘NE’ (May & Steiner 2017). The infrared emission is thus coincident with highly ionized, relatively hot gas. The atomic lines are associated with a shock “bubble” being excavated by the jet(May & Steiner 2017; Vermot et al. 2023).

Infrared interferometry also shows morphological similarity to the features near the center of the Fizeau image. In the *K'-band* bispectrum speckle interferometry images of (Weigelt et al. 2004), significant flux extends from the AGN toward the north ≈ 400 mas. The *K'-band* emission corresponds in scale and shape to the northern extended *N*-band flux we image. The authors give a rough flux range of 30–230 mJy; in our northern aperture; we measure ~100 mJy. The *K'-band* flux range is too large to reasonably estimate a dust temperature at this location. The VLTI/MATISSE observations(Gámez Rosas et al. 2022) and modeling(Leftley et al. 2024) indicate the continuation of NEX down to sub-parsec scales. In the modeling work(Leftley et al. 2024), this is the southern edge of a dusty (outflow) cone. In the model-independent imaging(Gámez Rosas et al. 2022), the *N*-band emission also shows an ×-shape similar to the ALMA emission (Impellizzeri et al. 2019) and the edges of the [OIII] emission. The match in PA of the LBTI Fizeau
NEX, the VLTI/MATISSE “×”, and the ALMA 256 GHz map(Impellizzeri et al. 2019) strongly indicates that this is
the same structure measured at different scales.

Using apertures with a diameter of 10.1 pc, we measure the flux at two image locations: the flux of the southern component (SN) is $F_{\text{South}} = 39.6 \pm 2.4$ mJy, with the aperture centered 25.1 pc to the south of the center; at the same distance to the north we measure $F_{\text{North}} = 100.6 \pm 6.2$ mJy. Using Eq. 1 from (Barvainis 1987) and the relation $F(r) = \sigma AT(r)^4$, we calculate an expected dust flux in an aperture with area $A$ as a function of radius (see Methods §4). The observed values are far above the simple prediction, as is readily seen in Fig. 4. At distances of 20 and 45 pc the flux excesses are associated with 5 GHz features ‘C’ and ‘NE.’ Spectral index measurements of cloud C show steep spectra indicative of optically thin, non-thermal emission between 5 and 21 GHz (Mutie et al. 2024). We use the spectral index and 5 GHz integrated flux of ‘C’ to estimate a jet flux at 8.7 µm of 0.06 mJy. The non-thermal contribution is thus quite small. It is therefore likely that a mixture of non-thermal emission and shock heating is necessary to produce the MIR fluxes. The exact contributions from non-thermal radiation and shock-heated dust are beyond the scope of this work. The correspondence between the different observations indicates that the MIR emission simultaneously traces AGN heated dust, jet-cloud interactions in ‘C’ (and possibly ‘S2’), and highly ionized jet shocks in the ‘tongue.’ Our observations lack precise astrometry, but independent morphological similarity to the radio, MIR, and [OIII] observations constrains our positions to within a few pixels ($\le 3$ px= 54 mas, or one PSF width).

## Constraints on the dusty torus

The disk+wind model consists of a geometrically thin, optically thick equatorial disk and a conical structure perpendicular to it (usually with hyperbolic or parabolic walls) representing a radiation-driven wind. This morphology was recently shown to be a good representation of the VLTI/GRAVITY and VLTI/MATISSE data for NGC 1068 (Leftley et al. 2024), but the interferometric data resolved out structures larger than ~ 5 pc. We computed a disk+wind radiative transfer (RT) model out to 15 pc using the SKIRT code (Camps & Baes 2020). The parameters of the RT modeling are given in Methods §5. The model was based on the one presented in (Isbell et al. 2022) with clumpy graphite dust in a hyperboloid cone and graphite-silicate dust in the disk. Surface brightness profiles of the resulting simulated and observed flux distributions are shown in Fig. 4.

Many observed features on < 15 pc scales (e.g., NEX, the dark band, and SN) and the N-S asymmetry of the large-scale emission are reproduced by the disk+wind model (shown in ED Fig. 4): near the central flux peak we see an edge-brightened outflow cone like NEX; south of the flux peak is a dark absorption band due to the densest part of the disk; and the ~ 75° inclined disk out to ≥ 25 pc causes significant extinction toward the south. The observed morphology, however, is not purely a disk+wind. There is much more flux directly north of the AGN in the images than in the modeling – this follows the path of the radio jet and is likely emission from jet-cloud interactions. Based on deviations from a centrally heated flux profile, correspondence with radio jet and shock morphology, and 5-21 GHz spectral indices of the sources, we claim that a large fraction of the MIR emission we recover is actually related to the radio jet (either directly or through e.g., shock heating) rather than the radiation-pressure-driven wind. Through comparisons to RT models as well as to multi-wavelength observations, we are able to spatially distinguish two separate AGN-host interaction channels at 3.4 pc resolution, showing a disk+wind(+jet) morphology. Previous VLT imaging (Isbell et al. 2021) supported the disk+wind model but did not have sufficient resolution to distinguish jet-related features.

The dusty wind, ionization cone, and radio jet are all presumed to be symmetric above and below the accretion disk. Based on correspondence with the radio data, we place the SMBH at the location marked with a green ‘×’ in Fig. 5, near the bottom of the brightest flux region. If we then assume that similar structure should extend to the north and south, it is immediately apparent that some absorbing medium is present – likely the disk-like component of the dusty torus. We estimate the optical depth of this dusty structure by relating the flux at an equal distance (accounting for inclination) above and below the SMBH (see Methods §6). Using the same apertures as above, we find that the southern emitting region is being obscured by an object with an optical depth of $\tau_{8.7} = 0.41 \pm 0.04$. This obscuring structure is assumed to be the disk-like component of the disk+wind model (e.g., Stalevski et al. 2019; Leftley et al. 2024). The optical depth of such a structure is expected to increase approaching the SMBH. We have measured the average value at 25.1 pc within a 10.1 pc diameter aperture. Assuming a linear change in $\tau$, we obtain
an optical depth gradient of $-0.05 \text{ pc}^{-1} \le \Delta\tau_{8.7\mu\text{m}} \le -0.02 \text{ pc}^{-1}$ when comparing to VLTI/MATISSE < 1 pc scale measurements (Gámez Rosas et al. 2022).

The dark band between the flux peak in our images and the southernmost emission is inferred to be due to an obscuring dust disk centered on 'S1.' VLTI/MATISSE observations revealed a similarly oriented obscuring dust disk (Gámez Rosas et al. 2022; Leftley et al. 2024, plotted as yellow ellipses in Fig. 3, panel d; reproduced from (Gámez Rosas et al. 2022)). A molecular torus observed with CO(2-1), CO(3-2), and HCO$^+$(4-3) (García-Burillo et al. 2019) is oriented approximately perpendicular to the radio jet, similar to the orientation of the obscuring "dark band" we and VLTI/MATISSE observe. We estimate that the $75^\circ$ inclined obscuring disk has a diameter $\ge 70$ pc in order to obscure SN. This is more than twice as large as the 30 pc CO disk measured by ALMA (García-Burillo et al. 2019); the discrepancy may indicate that CO traces the densest, inner parts of the disk, but a more extended region still contributes to obscuration.

Finally, using the assumed geometry and inclination of the circumnuclear structures, we estimate a lower limit on the dust density within the obscuring "dark band" to be $\rho_{\text{dust}} \ge 9 \times 10^{-27} \text{ g cm}^{-3}$. Compared to average ISM gas density, this corresponds to 60% of the standard ISM dust-gas-ratio of $\approx 1\%$ (Bohlin et al. 1978). A study of 36 X-ray observed AGN yielded dust-to-gas ratios in the range [1%-100%] the standard ISM value (Esparza-Arredondo et al. 2021). Our measurement falls within this range, but we note that the local density can be much higher if the dust is located within a thin disk as models predict. Future RT modeling is necessary to test the implications of the optical depth gradient on the disk's density profile.

The methods used herein can be applied to other AGN in order to understand the applicability of the disk+wind(+jet) as a unified model of AGN. Based on the measured fringe SNR of these observations, $324.26 \pm 76.17$ for the 20 cycles of NGC 1068 ($10.2 \pm 3.6$ Jy) and $307.22 \pm 28.17$ for the 7 good cycles of 1 Tau (11.2 Jy), we estimate that we could achieve SNR$\ge 10$ for sources $\ge 0.3$ Jy with roughly one-minute total integration. This is certainly optimistic due to decreased AO performance on faint sources, but repeated observations, more cycles, or fringe tracking on a nearby star could overcome this. We cross-matched the AGN catalog of (Véron-Cetty & Véron 2010) with the ALLWISE catalog (Cutri et al. 2021) to search for sources. Given the $0.3$ Jy N-band flux cutoff we use $W3$ as a proxy to identify 24 sources which are bright enough in N and have $V < 15$ for AO tracking. This sample of northern AGN includes roughly equal numbers of Seyfert 1 and Seyfert 2 systems, and would allow for a population study comparable in scope to the VLTI/MIDI sample (Burtscher et al. 2013; López-Gonzaga et al. 2016) which caused the paradigm shift to disk+wind models.

## Conclusions

In this paper we present extragalactic results from open loop Fizeau imaging with LBTI/NOMIC. We observed NGC 1068 with a $22.8 \times 8.4$ m synthetic aperture at $8.7 \mu\text{m}$; after PSF deconvolution we obtained a measured $47 \times 90$ mas ($3.3 \times 6.2$ pc) resolution. Through the resulting images we resolved substructures of the circumnuclear dust, tested the current paradigm of disk+wind AGN models, and spatially distinguished the dusty wind from jet-host interactions. More specifically, in this work we:

1. Developed the methods to deconvolve the PSF of complex, large-scale emission in LBTI/NOMIC Fizeau images,

2. Used the correspondence between the Fizeau images and existing observations of the 5 GHz, [OIII], NIR lines, and MIR dust emission to identify different emission mechanisms and to create a self-consistent schematic of the MIR dust emission from subparsec to $\sim 100$ pc scales,

3. Identified the NEX-a morphological feature which extends NE from the subparsec-scale MATISSE structure and traces the southern edge of the ionization cone- which we claim is the resolved onset of AGN feedback via dusty winds and in agreement with the current disk+wind unified model paradigm,

4. Measured the radial flux distribution of the $8.7 \mu\text{m}$ emission, identifying clear excesses beyond the radiative dust heating of the classical AGN torus model (by the accretion disk) at 20 and 45 pc from the AGN, and postulating that this excess emission is due to dust-jet interactions and shocks based on correspondence with the 5 GHz emission and NIR emission lines,
5. Measured the optical depth of the inferred obscuring disk of the disk+wind model to be $\tau_{8.7\mu m} = 0.41 \pm 0.04$ between 20 and 30 pc radial distance from the AGN. From this we infer a lower limit on the dust density of the obscuring material to be $\rho_{dust} \ge 9 \times 10^{-27} \text{ g cm}^{-3}$.

These results show the interplay between jet-mode and quasar-mode feedback, as we find signatures of dust-jet interaction and dusty, radiation-pressure driven winds in the ionization cone. The disentanglement of jet emission, shocks, and dusty winds is important for the interpretation of observations with telescopes such as JWST in which the processes are blended due to lower resolution. We additionally constrain the obscuring disk in the disk+wind flavor of the Unified Model of AGN. Future observations with this mode utilizing the L-band and/or other narrow N-band filters will allow us to measure spatially resolved dust temperatures and further constrain emission processes. Finally, we have presented unprecedented ELT-scale images of an extragalactic source, developing the methods and tools necessary to make this LBTI observing mode usable for general science cases and tests for upcoming facilities.
Main Text Figures

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Empirical PSF estimates. Panel a: Median PSF of 1 Tau throughout the night. Panel b: Median PSF of 1 Tau stacked and rotated to match the NGC 1068 field rotation. This serves as an estimate of the empirical PSF of the corotated and stacked science image. Panel c: Stacked and corotated Fizeau image of NGC 1068 combining the 20 nodding cycles. In all panels, North is up and East is left.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** LBTI Fizeau image of NGC 1068 after deconvolution. Panel a) Richardson-Lucy deconvolved LBTI Fizeau image of NGC 1068 at 8.7 µm using Method 1. Panel b) Contours of the same image with feature labels. Contours start at 90% of the peak flux density and decrease by a multiplicative factor of 2 down to a factor 128 below the maximum contour. This image was deconvolved from the stacked-corotated cycles using the Richardson-Lucy method. The resulting resolution FWHM is 46.8 × 90 mas (3.3 × 6.3 pc), shown as an ellipse in the left panel. Red and green labels indicate features from the literature (Gallimore et al. (2004) and Bock et al. (2000), respectively) and blue labels indicate features from the literature (Gallimore et al. (2004) and Bock et al. (2000), respectively) and green labels indicate features from the literature (Gallimore et al. (2004) and Bock et al. (2000), respectively).</figcaption></figure>
indicate new features from this work: the northern extended peak (NEP), the northeast extension (NEX), the dark band, and the southern N-band source (SN). The results of other deconvolution approaches are given in Methods.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of deconvolved LBTI Fizeau image to observations at other wavelengths. Deconvolved LBTI Fizeau image of NGC 1068 at 8.7 µm compared to observations of the 5 GHz radio emission (panel a; Gallimore et al. (2004)), the single-dish 7.9 µm emission (panel b; Bock et al. (2000)), and the ionization cone traced by [OIII] (panel c; Macchetto et al. (1994)). Panel d shows the MIR emission at sub-parsec scales as observed with VLTI/MATISSE (L-band in blue, M-band in green, and N-band in red; figure reproduced from (Gámez Rosas et al. 2022)). The “spur” to the northeast of center found in the LBTI image is seen at the subparsec scale in the VLTI/MATISSE image. The yellow ellipses represent the orientation of the obscuring dust disk. The red contours are the same as in Fig. 2.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

**Figure 4.** Surface brightness profiles of the observations and RT modeling. Distances have been de-projected using an inclination of 75° from face on. Error bars on the observed surface brightness come from uncertainty in flux calibration, typically ≲ 10%. We estimate an upper limit to the thermal dust flux as a function of radius (Methods §4). The dashed line shows predicted upper limits, and the dotted line shows upper limits obscured by $\tau_{8.7} = 0.41$. We also plot the measured flux values from a disk+wind radiation transfer model out to 15 pc. Corresponding features from Fig. 2 are labeled.

a

<figure><img src="image_7.png" /></figure>

**Figure 5.** Schematic summary of the sub-parsec to ≳ 100 pc scale processes measured with the MIR emission Panel a: A cartoon decomposition of the emission processes we identify. Panel b: The image from Fig. 1 with labels identifying features observed in different wavelengths. The direction of the radio jet (including the bend at cloud ‘C’) is shown in red. The edges of the ionization cone are shown in cyan. The location of the SMBH is marked with a green ‘×’ based on the ‘S1’ position in (Gallimore et al. 2004). Finally, an obscuring dust disk is indicated with white ellipses; this disk extends far enough to cover the assumed southern emission features with an optical depth of $\tau_{8.7\mu m} = 0.41 \pm 0.04$.
# Methods

## 1. OBSERVATIONS AND CALIBRATION

On 6 November 2023, we observed NGC 1068, 1 Tau, and Procyon using the NOMIC instrument on the LBTI in the Open-Loop Fizeau mode. We used the W08699-9_122 filter ($\lambda_c$ = 8.7 µm, range [8.13, 9.35] µm) and took a number of nodding cycles on the target and calibrators consisting of 4000 exposures of 0.014 seconds (2000 exposures on target, 2000 exposures on sky). Between 05:34:42 UTC and 08:48:58 observed 20 nodding cycles toward NGC 1068. Before and after the science target we observed 1 Tau and Procyon to serve as flux and point spread function (PSF) calibrators.

In order to identify good and bad chop-nod cycles, we inspect the peak flux of each cycle throughout the night. Bad cycles are expected to be those with significant thermal background variations or thin clouds. We plot the peak flux over time of both NGC 1068 and 1 Tau in Methods Supplementary Fig. 1. From this we see that 1 Tau produced the same peak flux values at the beginning and end of the night, indicating a stable thermal background. We can, however, identify two cycles with low peak fluxes. Upon inspection, we find these were affected by thin cirrus clouds causing a variable background and bad adaptive optics performance. We discard these when computing statistics for the 1 Tau PSF. For both 1 Tau and NGC 1068, the cycle-to-cycle flux variations are found to be similar to the peak flux variations within a cycle. We conclude that neither the mean PSF nor the flux calibration changes throughout the night. This allows us to work with stacked images.

The ideal PSF of the LBTI is the interference pattern of two 8.4 m circular apertures separated by 14 m. The resulting 22.8 m maximum baseline gives 78.7 mas resolution (5.5 pc) at 8.7 µm. The interference fringes exhibit a narrow central fringe peak, and a set of side lobes which rapidly decrease in amplitude as one moves further from the center. We refer to the interference fringe pattern as a Fizeau PSF. Over the span of approximately 3 hours on target, there was 59° field rotation. Field rotation causes the 22.8 m baseline of the LBTI to probe the target at different angles; after 180° field rotation this approximates a single 22.8 m telescope. Simulated field rotation on a real Fizeau PSF exhibits an approximately circular, 22.8 m PSF after ≈ 60° (Isbell et al. 2024). Robust estimates of both the individual cycle PSF and the PSF after field rotation are necessary for proper image deconvolution to recover the true flux distribution of the science target. The resulting corotated and stacked target and PSF calibrator are each shown in Fig. 1.

Fringe tracking was done in the so-called “open loop” mode, which means that not all frames are necessarily at the zero-order (or white-light) fringe and may jump due to atmospheric shifts to higher order fringes that are not useful for this study. However, the source was also imaged on the LMIRCam (Skrutskie et al. 2010) detector in the L band (2.8-4.1 µm) and the Fizeau PSF was spectrally dispersed. The interferometric phase on LMIRCam (which in Fizeau imaging has the same optical path delay as NOMIC) was set to the white-light fringe and phase shifts were corrected manually by keeping the fringes straight. While the NOMIC fringes moved due to atmospheric piston during and between exposures, this ensured that the movement was generally centered around the white-light fringe.

The Fizeau PSF may also move in the image plane or the two individual telescope PSFs may become slightly misaligned due to residual jitter from the AO system, though this effect is generally small compared to the 8.4 m N-band PSF and the precision at which the two PSFs can be overlapped (~10 mas differential tip-tilt jitter Ertel et al. (2018, 2020)). We account for both of these effects using “lucky Fizeau fringing” – i.e., we discard exposures with misaligned fringes or poorly overlapped telescope beams. We first compute the 2D cross-correlation function between the theoretical PSF and each frame of our observed PSF calibration stars. In each nodding cycle, we keep the frames which have the 10% largest cross-correlation values; in other words the 10% of frames most similar to the theoretical PSF (a similar approach was used in Leftley et al. (2021)). The selection of 10% is somewhat arbitrary, but larger values tend to decrease the peak SNR of the PSF image by “blurring” the PSF, and smaller values have too few frames to build up a robust empirical PSF estimate. Within each nodding cycle, the top 10% of frames are stacked, with each frame having been shifted to the pixel location with the highest 2D cross-correlation value, typically ≤ 3 pixels (pixel scale is 17.9 mas/px).

Two cycles of 1 Tau (cycles 5 and 6) showed significant deviations from the theoretical PSF, with broad, blurred artifacts. Upon closer inspection, these cycles suffered from light cirrus passing through. We therefore discard them as outliers. The remaining target and calibrator stacked PSFs show stability throughout the night. Finally, we compute the mean PSF for each calibration star. The mean PSF for 1 Tau is used in the following as our empirical PSF estimate because its flux (11.3 Jy; Cruzalèbes et al. (2019)) is more similar than Procyon's (74 Jy; Cruzalèbes et al. (2019))
to that of NGC 1068 (10.2 ± 3.6 Jy; Asmus et al. (2014)), and so noise estimates are similar. The empirical PSF is shown in Fig. 1, illustrating the characteristic triple-peak fringe pattern of the LBTI.

Good frames of NGC 1068 were selected in a similar way as for the PSF calibrators. Within each nodding cycle, we compute the 2D cross-correlation function between NGC 1068 and a nominal PSF (in this case, the measured PSF from the calibrator star). The frames with the 90th-percentile largest cross-correlation values are kept. As before, the target frames are re-centered to the pixel location with the peak 2D cross-correlation value and then stacked.

We obtain 20 stacked images of NGC 1068, one for each nodding cycle throughout the night. The resulting stacked image from each nodding cycle is shown in Methods. 10% of frames corresponds to 200 frames per cycle; this results in exposure time of 200 × 0.0137 s = 2.74 s per cycle or 54.8 s total exposure on NGC 1068. Before stacking each cycle, we corotated each sub-frame such that North is oriented along the y-axis in the frame. To create an estimate of the final, rotated PSF, we corotated the PSF calibrator frames by the same amount as the target sub-frames within each cycle and stacked.

We use the Mid-infrared stellar Diameters and Fluxes compilation Catalogue (MDC; Cruzalèbes et al. (2019)) to obtain the N-band fluxes of the calibrators. As is typical, we assume the unresolved stellar flux is distributed completely within the Fizeau PSF of the calibrator. This gives us a scaling between counts and Jy in each pixel that we apply to the science target. The results are consistent between Procyon (74 Jy) and 1 Tau (11.3 Jy), but we use the values from 1 Tau because it is closer to the nuclear flux of NGC 1068 (10.2 ± 3.6 Jy; Asmus et al. (2014)) and thus has a similar signal-to-noise ratio.

As LBTI Fizeau imaging is rather new to the community, we summarize the main differences from traditional, single aperture observations. First, the Fizeau PSF is characterized by a “triple peaked” fringe pattern, yielding 22.8 m resolution in one direction and 8.4 m resolution in the orthogonal direction. We show that even small amounts of field rotation allows one to obtain resolution better than 8.4 m in all directions, though a circular PSF would require 180° rotation (Isbell et al. 2024). Second, for non-ideal field rotation (and generally for the non-circular PSF of any segmented primary mirror), the effects of the resulting elongated PSF can be removed using image deconvolution. While image deconvolution is rather commonplace for observations with e.g., Keck(Bock et al. 2000) or JWST(Leist et al. 2024), we explore various methods below for the specific LBTI configuration. Third, we employed “lucky fringing,” keeping only the best 10% of data. This extremely conservative approach was feasible due to the brightness of the source, but recently published LBTI results have shown that ≥ 90% of frames can be kept in some cases for Fizeau imaging (Isbell et al. 2024), greatly increasing observation efficiency. For sources with K < 4.5 mag, fringe tracked imaging can be performed, removing the need for frame selection. A scheduled upgrade to the fringe tracking detector will increase this limit to K < 10 mag (Conrad et al. 2023). Finally, because the Fizeau aperture consists of two 8.4 m telescopes, the collecting area equivalent to that of a 11.9 m telescope (111 m²). Upcoming facilities, such as the ELT and GMT, will have larger effective collecting areas (978 m² and 368 m², respectively) and thus higher sensitivity than the LBTI.

## 2. IMAGE DECONVOLUTION

As can be clearly seen in Fig. 1, the stacked and corotated science cycles exhibit significant artifacts due to the Fizeau PSF. While extended flux features are still discernible, small-scale features near the center of the image are dominated by these artifacts. Using empirical measurements of the Fizeau PSF in individual cycles and for the stacked image, we can deconvolve the known PSF features to reveal the underlying source flux. This was previously done in (Leisenring et al. 2014) for observations of Io, but in the intervening years there have been changes to the adaptive optics performance of the LBTI(Pinna et al. 2016) and we used the NOMIC detector instead of LMIRCam. We therefore attempt Fizeau PSF deconvolution in three different ways: the first two use Richardson-Lucy deconvolution as in (Leisenring et al. 2014) and the third was developed for this work based on the Högbom CLEAN algorithm (Högbom 1974).

### 2.1. Richardson-Lucy Deconvolution

We proceed two different ways using Richardson-Lucy deconvolution (Richardson 1972; Lucy 1974):

1. Stack all corotated frames and deconvolve with a stacked, corotated PSF estimate, or

2. Deconvolve the image from each cycle with the PSF calibrator and then stack the resulting images.
The first method has the benefit of doing as few convolutions as possible (each introduces noise) at the expense of requiring that the estimate of the final PSF is robust. The second method uses a robust, single-cycle PSF but requires 20 (lower SNR) deconvolutions instead of one. In both cases, we use Richardson-Lucy deconvolution (from scikit-image_restore) because of its robustness to noise and efficiency. We find that both sequences give similar results, which is strong indication of the validity and robustness of the target flux distribution we recover. In practice, stacking before deconvolution (Method 1) increased the SNR of faint features, as noise is added destructively. This led to more stable results when changing the parameters n_iter and filteridar (the number of iterations and the clipping value for avoiding division by zero, respectively).

We estimate the uncertainty of each pixel in the resulting images using a bootstrap method; we perturb the empirical PSF by Gaussian noise and re-performed the deconvolution numerous time. For the images presented in Main we used Method 1 with n_iter=32 and filterid=0.005; we show the results in Fig. 2. The total flux of the image was not conserved during deconvolution, and therefore the deconvolved image is re-normalized so that the total flux of the resulting image equals the total flux of the corotated+stacked image (specifically, it is scaled by a factor of 0.7). The resulting PSF is 46.8 × 90 mas FWHM.

As expected from previous observations, we see significant extended flux along the north-south direction as well as a bright point source to the northeast at 600 mas (45 pc). However, the 78 mas PSF from the 22.8 m baseline allows us to resolve inner substructure near the AGN core. Here we see a “spur” of flux pointing toward the northeast out to ~10 pc. In the deconvolved-then-stacked (Method 2) image, we see evidence of an additional, smaller spur to the northwest. These spurs are found near the expected AGN location and are oriented in a similar direction as extended N-band dust measured on smaller scales with MATISSE (Gámez Rosas et al. 2022). In Method 2, we find that the bad channel (visible in Supplementary Fig. 2 in the upper third of the images) causes more significant artifacts than in Method 1. In ED Fig. 1, we see a negative contour at the location of the bad channel when using Method 2. It seems that the stacking of images in Method 1 washes out the effects of this channel. Due to the stability of the PSF throughout the night which allows for robust image stacking and the strong artifacts present in Method 2, we use Method 1 for analysis throughout this article. We note, however, that results only change by a few percent when using the other deconvolution method.

## 2.2. CLEAN-inspired method

The Högbom CLEAN algorithm is an aperture synthesis algorithm, developed to deconvolve a known “beam” from a measured radio emission map affected by prominent and extended sidelobes (Högbom 1974). Similarly, the known ideal Fizeau PSF produced by the LBTI has prominent and extended sidelobes. Due to the high fidelity of CLEANed radio images, we developed an experimental implementation for LBTI Fizeau images which iteratively removes the known, extended PSF.

We use the following notation: $Im_0$ is the initial, “dirty” image to deconvolve, beam is the known PSF estimate, gain is scaling applied to the PSF to remove flux from the dirty image in each iteration, CUTOFF is the flux level to which the image is cleaned, MAXITER is the maximum number of iterations to compute, and peak is the maximum flux in the current image, $I_i$. The algorithm then proceeds as follows:

START

WHILE peak > CUTOFF AND n_iter < MAXITER

1. Find peak flux in image, at coordinates $(x_p, y_p)$

2. Create $Im_{i+1}$

* Create a zero-image with same dimensions as $Im_0$

* Place a point source with flux = gain × peak at $(x_p, y_p)$ and convolve with beam

* Subtract this from $Im_i$ to get $Im_{i+1}$

3. Put a point source at $(x_p, y_p)$ in the CLEAN image

4. $n_{iter} = n_{iter} + 1$

ENDWHILE

END
We applied the above method with gain = 10⁻⁴ and MAXITER = 10⁶ on each of the 20 exposure cycles. We also set CUTOFF = Mean(Im₀), i.e., stop when the iteration peak is just above the background level. We finally corotate and stack the images as above. The resulting, deconvolved image is shown in ED Fig. 2. The results of this method show more low-surface brightness features than the Richardson-Lucy deconvolution in Fig. 2; in particular, the flux extending to the south of the photocenter and the flux within the northernmost arc are more prominent in the CLEAN-inspired deconvolution. Importantly, the “spur” extending to the NE of the photocenter appears in this deconvolution as well, with a similar opening angle and length. The fact that this feature appears in this independent deconvolution method gives yet more credence to its fidelity.

### 2.3. Neural Network Deconvolution

The method presented here is a new algorithm developed to perform image deconvolution without using an experimental PSF. Instead, we use a neural network trained on a large dataset of simulated LBTI observations. This method is derived from the interferometric image reconstruction algorithm proposed in (Sanchez-Bermudez et al. 2022). The results are presented in ED Fig. 3.

The algorithm steps are as follows:

1. Generate a large set of mock astrophysical images.

2. Simulate random LBTI PSFs, convolve the images with these PSFs, and introduce noise.

3. Train a neural network to restore the original images using only the simulated data.

4. Apply the trained neural network to the actual astronomical observation.

We elaborate on these steps in the subsequent paragraphs.

## I. Generation of mock astrophysical images

The objective is to rapidly produce a diverse set of images containing varied structures with different spatial scales. We create each image by summing several component images, each representing structures at a characteristic spatial scale, $\lambda$, randomly selected between 1 pixel and the full field of view (FoV). The number of component images, N, varies randomly from 1 to 20. For each component we:

1. Compute a 2D discrete exponential power spectrum P with characteristic scale $\lambda$.

2. Assign a random phase $\phi$ to each pixel.

3. Calculate the real part of the Fourier Transform of $P \times e^{i\phi}$, and mask the negative values.

The N images are combined using random weights, and the total flux of the resultant image is normalized to unity. This process is repeated to generate 2²⁰ unique images of dimensions ($n_{pix} \times n_{pix}$), with $n_{pix} = 100$.

## II. Simulation of Observations

Each generated image is then used to simulate an observation cube with dimension ($n_{frames}, n_{pix}, n_{pix}$), with $n_{frames} = 32$. The decision to use data cubes with $n_{frames}$ rather than a single frame is driven by two factors: first, this approach allows the neural network to prioritize data from the highest quality frames, thereby optimizing image selection; second, it enables differentiation of PSF variations—which change from frame to frame—from static astrophysical elements.

Instantaneous PSFs are computed as the Fourier Transform of $Pup \times \exp(i \times \phi_{res})$, with $Pup$ the telescope pupil and $\phi_{res}$ the residual atmospheric phase post-adaptive optics correction. The pupil is modeled as two 8.4 m disks, each with a central 0.9 m obstruction, and a 14.4 m center-to-center separation. $\phi_{res}$ is calculated from a sum of Zernike polynomials with coefficients randomly selected following a power-law distribution $P(n) \propto n^{-\alpha}$, with n the radial order of the polynomial.

A number $n_{\tau}$ of these instantaneous PSFs, drawn uniformly between 1 and 100, are averaged to create a semi-long exposure PSF. Each image is then convolved with $n_{frames}$ such PSFs, creating an observation data cube.

Two types of normal noise are added to the datacube, one proportional to the square root of the surface brightness of the images (photon noise of the source), and one uniform across the FoV (photon noise of the subtracted sky). The overall noise level is randomly drawn for each observation to produce a SNR between 1 and 50 for individual frames.
### III. Neural Network and training

After generating the training dataset, consisting of 2²⁰ random images and the corresponding observations, we train a neural network to reconstruct the original images from these convolved, simulated datasets.

We use a simplified version of the U-NET architecture, which is fully convolutional and includes skip layers. A schematic is shown in Supplementary Fig. 4.

The training is done on 50 epochs, each one including the entire dataset in batches of 16 observations. We use the Adam optimizer with a mean absolute loss function and an initial learning rate of 5×⁻⁴. This learning rate is halved each time there is no improvement in the loss after five consecutive epochs. The final mean absolute error, indicating the difference between the network’s prediction and the actual original images, is 0.09 for both the validation and the training dataset, demonstrating that the training did not result in overfitting.

### IV. Application to the actual observations

The observations consist of 20 cycles, corresponding to 20 different orientations of the source in the plane of the sky, each of these cycles corresponding to 2000 individual short-exposure frames.

We split each of these 2000 frame sets into 62 data cubes each with 32 frames, which were fed to the trained neural network, hence producing 62 deconvolved images for each orientation. These 62 deconvolutions are then averaged to produce one image per orientation. The 20 resulting images are then de-rotated, and the final image is computed as the median of these images. It is shown in ED Fig. 3. The low surface-brightness features are better recovered than in Fig. 2, and resemble those of the CLEAN approach in ED Fig. 2. The central emission is more compact than the CLEAN deconvolution results, but shows a similar NW-SE extension. The NEX feature, however, is not as readily apparent in these results, though the flux envelope near the center traces roughly the opening angle of the ionization cone.

#### 2.4. Estimating Image Significance

The standard deviation of the cycle PSFs is used as an uncertainty estimate, allowing us to bootstrap errors of the deconvolution we perform on the target frames. In Supplementary Fig. 5 we show all 1 Tau cycles’ PSF images and flux profiles. We compute the per pixel mean and standard deviation of all cycles (excluding cycles 5 and 6 due to thin cirrus). In 1000 bootstrap iterations, we perturb the empirical PSF (Fig. 1) with random noise set by the standard deviation at each pixel. We then stack and corotate the perturbed PSF to match the rotations of the 20 target observing cycles. From the 1000 iterations, we can compute the uncertainty per pixel for the recovered deconvolved images below; from this we conclude that all primary structures in the image (the central component with a “spur”, the southern point source, the northern extended structure, and the northernmost point source) are robustly recovered at $\gtrsim 10\sigma_{bootstrap,px}$. As an example, we display the signal-to-noise ratio (SNR) map for deconvolution Method 1 with 8.7 µm image contours overlaid in Supplementary Fig. 6; this shows that all primary features of the image are found at SNR > 10.

### 3. (NON)THERMAL DUST EMISSION ALONG THE RADIO JET

With apertures having a diameter of 8 pixels (0.144 arcsec, 10.1 pc), we measure the flux at several image locations. We measure the flux of the southern component (SN) to be $F_{South} = 39.6 \pm 2.4$ mJy, with the aperture centered 25.1 pc to the south of the center. At the same distance to the north (in cloud ‘C’), we measure the flux to be $F_{North} = 100.6 \pm 6.2$ mJy.

Utilizing Eq. 1 from (Barvainis 1987), we can estimate dust temperature as a function of radius:

$$T_{gr}(r) = 1650 \left( \frac{L_{acc}}{r^2} \frac{\text{pc}^2}{10^{10} L_{\odot}} \right)^{1/5.6} e^{-\tau_{uv}/5.6} \text{K}, \quad (1)$$

where $L_{acc}$ is the luminosity of the accretion disk in $L_{\odot}$, r is the distance from the accretion disk in parsec, and $\tau_{uv}$ is the optical depth to the ultraviolet continuum. The luminosity of grains of this temperature is then $L_{\nu,\Omega}(r) = \sigma/\nu\pi \times T_{gr}(r)^4$ W m⁻² Hz⁻¹ sr⁻¹. The predicted thermal dust flux in an aperture is $F_{8.7\mu m,therm} = L_{\nu,\Omega}(r)A$, where A is the aperture area. We show in Fig. 4 that there are significant flux excesses 20 and 45 pc from the core which are associated with 5 GHz features ‘S2,’ ‘C,’ and ‘NE.’
As stated in the Main Text, the contribution from non-thermal sources within cloud C is very small in the MIR. The authors in (Villar-Martín et al. 2001) present several possible dust-heating mechanisms within shocks for NGC 1068, but there has been some disagreement in the literature about the prevalence of shocks along the radio jet in NGC 1068. The authors of (Kraemer et al. 1998) claim that shocks are necessary to heat the dense clouds at their “position 2,” (cloud C in Gallimore et al. (2004)) but (Martins et al. 2010) claim photoionization actually could be sufficient with different assumptions of cloud densities. More recently, (Kakkad et al. 2018) show high cloud column densities, indicating that shocks are indeed required.

The MIR emission of NGC 1068 in the central 50 pc seemingly requires contributions from several distinct thermal and non-thermal processes. Future LBTI observations of this and other AGN, either with a number of L-, M-, or N-band filters to get dust color temperatures or with narrow filters targeting emission lines (e.g., [S IV] 10.5 µm and [Ne II] 12.8 µm), will be able to distinguish the emission or heating mechanisms.

## 4. RADIATION TRANSFER MODELING

The models which most closely match interferometric data of nearby Seyferts (e.g., (Stalevski et al. 2019; Leftley et al. 2024)) consist of a compact, dusty disk and a hollow hyperbolic cone extending in the polar direction. Here we modify the best-fitting disk+wind model of (Stalevski et al. 2019) to use the luminosity of NGC 1068 (10$^{12}$ L$_{⊙}$) and the opening angle and inclination fit by (Leftley et al. 2024). We then increase the box size of the pre-existing small-scale circumnuclear model to measure the flux at 15 pc scale. Except when otherwise listed, the parameters of the disk+wind model are the same as used in (Isbell et al. 2022). We use $i = 75^{\circ}$, a hyperbolic cone opening angle $\theta_{open} = 45^{\circ}$, an outer extent of the hyperbolic cone of $r_{hyp,out} = 24$ pc, a disk optical depth of $\tau_{disk,9.7} = 15$ (the same as in (Isbell et al. 2022)), and a disk outer radial extent of $r_{disk,out} = 10$ pc. We also set the instrument field of view to 30 pc and the number of pixels in each direction to 600 pixels. The resulting model is shown in ED Fig. 4, and we also extract radial flux profiles to compare to our measured emission profiles.

The main goal of the RT modeling was to determine whether the NEX in our images could realistically be centrally-heated dust in a windy outflow. Morphologically, there is indeed a similar structure due to edge-brightening of the wind. To test whether the flux was plausible we placed an aperture (3.8 pc radius) at approximately the same location in both the Fizeau image and the RT model. The aperture was placed on the southern edge of the ionization cone in both cases, 7 pc from the nucleus. The RT model aperture contains ~100 mJy of flux with a peak flux density of 1.2 mJy/px. The observation aperture contains $49 \pm 3$ mJy with a peak flux density of 1.5 mJy/px. While a detailed comparison of the resulting fluxes is not meaningful, the disk+wind RT model shows that it is plausible that the emission in the NEX is AGN heated dust at the edge of the outflow cone. A detailed RT model parameter study, comparing to the combined VLTI/GRAVITY, VLTI/MATISSE, and LBTI data, will be the subject of future work.

## 5. OPTICAL DEPTH

The dusty wind, ionization cone, and radio jet are all presumed to be symmetric above and below the accretion disk. Based on correspondence with the radio data, we place the SMBH at the location marked with a green ‘×’ in Fig. 5, near the bottom of the brightest flux region. If we then assume that similar structure should extend to the north and south, it is immediately apparent that some absorbing medium is present – likely the disk-like component of the dusty torus.

The optical depth is given by

$$ \tau = -\ln \frac{F_{\text{South}}}{F_{\text{North}}} = lN\sigma = l\kappa\rho, \qquad (2) $$

where $F_{\text{South}}$ and $F_{\text{North}}$ are the fluxes equidistant south and north of the SMBH, $N$ is the number density of the obscuring material, $\sigma$ is the material’s cross-section, $\kappa$ is the mass extinction coefficient of the material, $l$ is the path length of the absorbing material, and $\rho$ is the density of the material. Assuming a specific dust composition, we are left with two free parameters representing the density of the dust and the path length (i.e., the thickness of the obscuring structure). Taking the measurements from the same apertures as above (Methods §4), we find that the southern emitting region would need to be obscured by an object with an optical depth of $\tau_{8.7} = 0.41 \pm 0.04$.

The optical depth of such a structure should increase as one gets closer to the SMBH. We measure the average value at 25.1 pc within a 10.1 pc diameter aperture. Previous measurements (Gámez Rosas et al. 2022) at $\lesssim 1$ pc give $\tau_{8.7\mu m} = 0.94 - 1.51$ (their measurement is at a slightly longer wavelength, so we scaled their reported optical depth by a factor 0.74 to compare at 8.7 µm, based on the extinction curve from Schartmann et al. (2005)). Assuming a linear change in $\tau$, we obtain an optical depth gradient of $-0.05 \text{ pc}^{-1} \le \Delta\tau_{8.7\mu m} \le -0.02 \text{ pc}^{-1}$.
Since the optical depth $\tau$ is related to dust density, composition, and path length through the obscuring material, we explore what a gradient might imply. The dust composition is likely approximately constant, so we focus on dust density and path length. Current disk+wind RT models should naively lead to an *increased* optical depth at larger radii because the constant-density disk becomes thicker with an opening angle of 5°, leading to a longer path length. We observe, however, the opposite trend. The RT modeling in this work assumes a constant dust density in the disk (based on comparisons to Circinus on much smaller scales Stalevski et al. (2019); Isbell et al. (2023)). Our new measurements imply the necessity of a decreasing density gradient counteracting the increase in path length. An exploration of RT parameters would be greatly aided by more constraints on dust temperatures and emission mechanisms, so this will be explored in future work after LBTI observations with more filters have been obtained.

Rearranging Eq. 2, we find the relation $\rho = \tau_{8.7 \ \mu m}/(l\kappa_{8.7 \ \mu m})$. Because NGC 1068 is inclined $\approx 75^\circ$ from face-on, we can compute an upper limit on the path length ($l$) corresponding to a lower limit on the dust density: $\rho_{lower} = \tau_{8.7 \ \mu m} \sin i/(d \ \kappa_{8.7 \ \mu m})$. For $i = 75^\circ$, $\kappa_{8.7 \ \mu m} = 12116.5 \ \text{m}^2 \ \text{kg}^{-1}$ from (Schartmann et al. 2005), $\tau_{8.7 \ \mu m} = 0.41 \pm 0.04$ as above, and $d \approx 25$ pc as measured from the image, we obtain $\rho_{dust} \ge 9 \times 10^{-27} \ \text{g cm}^{-3}$. The ISM has a mean gas density of 1 atom per cubic centimeter, or $\bar{\rho}_{ISM} = 1.6 \times 10^{-24} \ \text{g/cm}^{-3}$ (Draine 2011). This corresponds to a dust-to-gas ratio 25 pc from the AGN of $6 \times 10^{-3}$. This is 0.6 times the mean ISM dust-gas-ratio of $\approx 1\%$ (Bohlin et al. 1978). A study of 36 X-ray observed AGN yielded dust-to-gas ratios in the range [0.01-1] times the standard ISM value (Esparza-Arredondo et al. 2021); our measurement falls well within this range.

**Data availability:** Raw and processed data corresponding to the current study are available from the corresponding author upon request.

**Code availability:** The algorithms used for data processing are described in this work, but data processing Python scripts and data processing instructions are available within the public github repository located at https://github.com/jwisbell/lbti_fizeau.

**Acknowledgments:** We thank the anonymous referees for their feedback which greatly improved this work. VFG acknowledges funding from the National Aeronautics and Space Administration through the Exoplanet Research Program under Grant No. 80NSSC21K0394 (PI: S. Ertel). RGP, NM, JL, and PV been supported by the French Agence Nationale de la Recherche (ANR) through the grant "AGN MELBa" ANR-21-CE31-0011.

**Author contributions:** JWI, J-UP, GW, MS, JLef, WJ, RP, NM, and VGR contributed to scientific analysis, modeling, and interpretation. JWI and PV developed and tested deconvolution methods. JWI, SE, AB, JC, VFG, WFH, JLei, JP, and KW took part in observations setup, execution, and processing.

**Competing interests:** The authors declare that there are no competing interests.
Extended Data and Supplementary Information

<figure><img src="image_8.png" /><figcaption>**Extended Data Figure 1.** LBTI Fizeau image of NGC 1068 after deconvolution using Method 2. *Panel a*) Deconvolved LBTI Fizeau image of NGC 1068 at 8.7 µm *Panel b*) The [2,3,4,5,10,20,50]×σ contours of the same image. Method 2 used Richardson-Lucy on each observing cycle and then stacking. The results are overall similar to Method 1 (Fig. 2), but with slightly narrower features and a negative flux contour to the north due to a bad pixel channel (visible in Supplementary Figure 2).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Extended Data Figure 2.** Results of the CLEAN-inspired deconvolution algorithm. *Panel a*) the point source locations produced by the CLEANing iterations. *Panel b*) the CLEAN-inspired image convolved with the restoring beam. *Panel c*) the [8, 16, 32, 64]×σ contours of the restored CLEAN-image, with σ measured from the image background. The FWHM of the restoring beam is given in the bottom left corner. Resulting structures resemble those of Fig. 2, albeit with more prominent low-surface-brightness features.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Extended Data Figure 3.** Stacked and corotated results of the neural network deconvolution. *Left*) surface brightness of the recovered emission, normalized to the peak flux. *Right*) Contours of the surface brightness, beginning from 90% of the peak and decreasing by factors of 2 down to a factor 128 below the maximum contour. Low surface brightness features similar to those in Extended Data Fig. 2 are recovered, though the central emission appears to be slightly less extended.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Extended Data Figure 4.** Comparison between radiation transfer modeling results and the LBTI observation. *Left*) Radiation transfer of a disk+wind model. Contours are set at 4, 8 × 10⁻¹⁰ W m⁻²sr⁻¹ to emphasize the edge-brightening effect and to show that a similar structures to the NEX and southernmost MIR flux we observe are reproduced. *Right*) The same as Fig. 2. In both panels, the white lines represent the extraction aperture used to produce the surface brightness profile in Fig. 4.</figcaption></figure>
<figure><img src="image_12.png" /></figure>

**Supplementary Figure 1.** Transfer function throughout the night. The measured flux of both the calibrator and the target are shown to vary little throughout the night. Each point is the peak flux of the fringe pattern, and error-bars represent the standard deviation of the flux within each frame. They are comparable to the statistical flux variations shown over the night. The two points in red were affected by thin cirrus clouds and are discarded.
<figure><img src="image_13.png" /><figcaption>**Supplementary Figure 2.** Individual exposure cycles of NGC 1068 after frame selection. The selected frames are stacked within each cycle. The summed exposure time is approximately 2.75 s for each cycle, resulting in a fringe peak SNR of ≈ 320 in each stacked cycle image. The effects of field rotation are not corrected in these panels.</figcaption></figure>
<figure><img src="image_14.png" /></figure>

**Supplementary Figure 3.** Comparison of CLEAN-deconvolved LBTI Fizeau image to observations at other wavelengths. Deconvolved LBTI Fizeau image of NGC 1068 at 8.7 µm compared to observations of the 5 GHz radio emission (panel a; [18]), the single-dish 7.9 µm emission (panel b; [29]), and the ionization cone traced by [OIII] (panel c; [28]). Panel d shows the MIR emission at sub-parsec scales as observed with VLTI/MATISSE (L-band in blue, M-band in green, and N-band in red; figure reproduced from [13]). The “spur” to the northeast of center found in the LBTI image is seen at the subparsec scale in the VLTI/MATISSE image. The yellow ellipses represent the orientation of the obscuring dust disk. The red contours are the same as in Fig. 2.
<figure><img src="image_15.png" /><figcaption>**Supplementary Figure 4.** Schematic representation of the U-NET architecture used for deconvolution. Starting from an observation cube with dimensions ($n_{pix} \times n_{pix} \times n_{frames}$), the convolutional network concatenates layers both with and without MaxPool2D. The convolutions on different scales are then concatenated into a layer with dimensions ($n_{pix} \times n_{pix} \times 96$) before a final set of convolutions which produce the output image of dimension ($n_{pix} \times n_{pix}$).</figcaption></figure>
<figure><img src="image_16.png" /></figure>

**Supplementary Figure 5.** Individual PSF cycles of 1 Tau. Cycles 1 and 2 were obtained before the science target, and cycles 3-9 were obtained immediately after. Cycles 8 and 9 are presented in red because they are excluded from use in calibration due to cirrus clouds passing by during those exposures.
<figure><img src="image_17.png" /><figcaption>**Supplementary Figure 6.** Bootstrapped SNR estimate from uncertainties in the PSF. The black contours are the same 8.7 µm image contours as in Figs. 2 and 3 ([2,3,4,5,10,20,50]×σ). The high SNR region (in yellow) encloses all major features examined in this work.</figcaption></figure>
REFERENCES

Antonucci, R. 1993, Annu. Rev. Astron. Astrophys., 31,
473, doi: 10.1146/annurev.aa.31.090193.002353

Asmus, D., Hönig, S. F., Gandhi, P., Smette, A., & Duschl,
W. J. 2014, Mon. Not. R. Astron. Soc., 439, 1648,
doi: 10.1093/mnras/stu041

Barvainis, R. 1987, Astrophys. J., 320, 537,
doi: 10.1086/165571

Bock, J. J., Neugebauer, G., Matthews, K., et al. 2000,
Astron. J., 120, 2904, doi: 10.1086/316871

Bohlin, R. C., Savage, B. D., & Drake, J. F. 1978,
Astrophys. J., 224, 132, doi: 10.1086/156357

Burtscher, L., Meisenheimer, K., Tristram, K. R. W., et al.
2013, Astron. Astrophys., 558, A149,
doi: 10.1051/0004-6361/201321890

Camps, P., & Baes, M. 2020, Astronomy and Computing,
31, 100381, doi: 10.1016/j.ascom.2020.100381

Conrad, A., Stone, J., & Ertel, S. 2023, in Adaptive Optics
for Extremely Large Telescopes 7th Edition, ONERA,
Avignon, France, doi: 10.13009/AO4ELT7-2023-039

Conrad, A., de Kleer, K., Leisenring, J., et al. 2015, Astron.
J., 149, 175, doi: 10.1088/0004-6256/149/5/175

Cruzalèbes, P., Petrov, R. G., Robbe-Dubois, S., et al.
2019, Mon. Not. R. Astron. Soc., 490, 3158,
doi: 10.1093/mnras/stz2803

Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR
Online Data Catalog: AllWISE Data Release (Cutri+
2013), VizieR On-line Data Catalog: II/328. Originally
published in: IPAC/Caltech (2013)

de Kleer, K., Skrutskie, M., Leisenring, J., et al. 2017,
Nature (London), 545, 199, doi: 10.1038/nature22339

Draine, B. T. 2011, Physics of the Interstellar and
Intergalactic Medium (Princeton University Press)

Ertel, S., Absil, O., Defrère, D., Augereau, J. C., &
Mennesson, B. 2018, Experimental Astronomy, 46, 401,
doi: 10.1007/s10686-018-9600-7

Ertel, S., Hinz, P. M., Stone, J. M., et al. 2020, in Society
of Photo-Optical Instrumentation Engineers (SPIE)
Conference Series, Vol. 11446, Optical and Infrared
Interferometry and Imaging VII, ed. P. G. Tuthill,
A. Mérand, & S. Sallum, 1144607,
doi: 10.1117/12.2561849

Esparza-Arredondo, D., Gonzalez-Martín, O., Dultzin, D.,
et al. 2021, Astron. Astrophys., 651, A91,
doi: 10.1051/0004-6361/202040043

Fabian, A. C. 2012, Annu. Rev. Astron. Astrophys., 50,
455, doi: 10.1146/annurev-astro-081811-125521

Gallimore, J. F., Baum, S. A., & O'Dea, C. P. 2004,
Astrophys. J., 613, 794, doi: 10.1086/423167

Gallimore, J. F., & Impellizzeri, C. M. V. 2023, Astrophys.
J., 951, 109, doi: 10.3847/1538-4357/acd846

Gámez Rosas, V., Isbell, J. W., Jaffe, W., et al. 2022,
Nature, 602, 403, doi: 10.1038/s41586-021-04311-7

García-Burillo, S., Combes, F., Ramos Almeida, C., et al.
2016, Astrophys. J. Lett., 823, L12,
doi: 10.3847/2041-8205/823/1/L12

—. 2019, Astron. Astrophys., 632, A61,
doi: 10.1051/0004-6361/201936606

GRAVITY Collaboration, Pfuhl, O., Davies, R., et al. 2020,
Astron. Astrophys., 634, A1,
doi: 10.1051/0004-6361/201936255

Hinz, P. M., Defrère, D., Skemer, A., et al. 2016, in Society
of Photo-Optical Instrumentation Engineers (SPIE)
Conference Series, Vol. 9907, Optical and Infrared
Interferometry and Imaging V, ed. F. Malbet, M. J.
Creech-Eakman, & P. G. Tuthill, 990704,
doi: 10.1117/12.2233795

Hoffmann, W. F., Hinz, P. M., Defrère, D., et al. 2014, in
Society of Photo-Optical Instrumentation Engineers
(SPIE) Conference Series, Vol. 9147, Ground-based and
Airborne Instrumentation for Astronomy V, ed. S. K.
Ramsay, I. S. McLean, & H. Takami, 914710,
doi: 10.1117/12.2057252

Högbom, J. A. 1974, Astron. Astrophys. Suppl., 15, 417

Impellizzeri, C. M. V., Gallimore, J. F., Baum, S. A., et al.
2019, Astrophys. J. Lett., 884, L28,
doi: 10.3847/2041-8213/ab3c64

Isbell, J. W., Burtscher, L., Asmus, D., et al. 2021,
Astrophys. J., 910, 104, doi: 10.3847/1538-4357/abdfd3

Isbell, J. W., Meisenheimer, K., Pott, J.-U., et al. 2022,
A&A, 663, A35, doi: 10.1051/0004-6361/202243271

Isbell, J. W., Pott, J. U., Meisenheimer, K., et al. 2023,
Astron. Astrophys., 678, A136,
doi: 10.1051/0004-6361/202347307

Isbell, J. W., Ertel, S., Wagner, K., et al. 2024, in Optical
and Infrared Interferometry and Imaging IX, ed.
J. Kammerer, S. Sallum, & J. Sanchez-Bermudez, Vol.
13095, International Society for Optics and Photonics
(Proc. SPIE), 1309506, doi: 10.1117/12.3027270

Kakkad, D., Groves, B., Dopita, M., et al. 2018, Astron.
Astrophys., 618, A6, doi: 10.1051/0004-6361/201832790

Kraemer, S. B., Ruiz, J. R., & Crenshaw, D. M. 1998,
Astrophys. J., 508, 232, doi: 10.1086/306389

Leftley, J. H., Hönig, S. F., Asmus, D., et al. 2019,
Astrophys. J., 886, 55, doi: 10.3847/1538-4357/ab4a0b

Leftley, J. H., Tristram, K. R. W., Hönig, S. F., et al. 2021,
Astrophys. J., 912, 96, doi: 10.3847/1538-4357/abee80
Leftley, J. H., Petrov, R., Moszczynski, N., et al. 2024, Astron. Astrophys., 686, A204, doi: 10.1051/0004-6361/202348977

Leisenring, J. M., Hinz, P. M., Skrutskie, M., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9146, Optical and Infrared Interferometry IV, ed. J. K. Rajagopal, M. J. Creech-Eakman, & F. Malbet, 91462S, doi: 10.1117/12.2057290

Leist, M. T., Packham, C., Rosario, D. J. V., et al. 2024, Astron. J., 167, 96, doi: 10.3847/1538-3881/ad1886

López-Gonzaga, N., Burtscher, L., Tristram, K. R. W., Meisenheimer, K., & Schartmann, M. 2016, Astron. Astrophys., 591, A47, doi: 10.1051/0004-6361/201527590

Lucy, L. B. 1974, Astron. J., 79, 745, doi: 10.1086/111605

Macchetto, F., Capetti, A., Sparks, W. B., Axon, D. J., & Boksenberg, A. 1994, Astrophys. J. Lett., 435, L15, doi: 10.1086/187583

Martins, L. P., Rodríguez-Ardila, A., de Souza, R., & Gruenwald, R. 2010, Mon. Not. R. Astron. Soc., 406, 2168, doi: 10.1111/j.1365-2966.2010.17042.x

May, D., & Steiner, J. E. 2017, Mon. Not. R. Astron. Soc., 469, 994, doi: 10.1093/mnras/stx886

Mutie, I. M., Williams-Baldwin, D., Beswick, R. J., et al. 2024, Mon. Not. R. Astron. Soc., 527, 11756, doi: 10.1093/mnras/stad3864

Netzer, H. 1987, Mon. Not. R. Astron. Soc., 225, 55, doi: 10.1093/mnras/225.1.55

Pinna, E., Esposito, S., Hinz, P., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9909, Adaptive Optics Systems V, ed. E. Marchetti, L. M. Close, & J.-P. Véran, 99093V, doi: 10.1117/12.2234444

Richardson, W. H. 1972, Journal of the Optical Society of America (1917-1983), 62, 55

Sanchez-Bermudez, J., Mérand, A., Sallum, S., et al. 2022, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 12183, Optical and Infrared Interferometry and Imaging VIII, ed. A. Mérand, S. Sallum, & J. Sanchez-Bermudez, 121831G, doi: 10.1117/12.2629430

Schartmann, M., Meisenheimer, K., Camenzind, M., Wolf, S., & Henning, T. 2005, Astron. Astrophys., 437, 861, doi: 10.1051/0004-6361:20042363

Skrutskie, M. F., Jones, T., Hinz, P., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 77353H, doi: 10.1117/12.857724

Stalevski, M., Tristram, K. R. W., & Asmus, D. 2019, Mon. Not. R. Astron. Soc., 484, 3334, doi: 10.1093/mnras/stz220

Vermot, P., Barna, B., Ehlerová, S., et al. 2023, Astron. Astrophys., 678, A206, doi: 10.1051/0004-6361/202347227

Véron-Cetty, M. P., & Véron, P. 2010, Astron. Astrophys., 518, A10, doi: 10.1051/0004-6361/201014188

Villar-Martín, M., De Young, D., Alonso-Herrero, A., Allen, M., & Binette, L. 2001, Mon. Not. R. Astron. Soc., 328, 848, doi: 10.1046/j.1365-8711.2001.04916.x

Wada, K. 2012, Astrophys. J., 758, 66, doi: 10.1088/0004-637X/758/1/66

Wada, K., Schartmann, M., & Meijerink, R. 2016, Astrophys. J. Lett., 828, L19, doi: 10.3847/2041-8205/828/2/L19

Weigelt, G., Wittkowski, M., Balega, Y. Y., et al. 2004, Astron. Astrophys., 425, 77, doi: 10.1051/0004-6361:20040362

Williamson, D., Hönig, S., & Venanzi, M. 2019, Astrophys. J., 876, 137, doi: 10.3847/1538-4357/ab17d5
—. 2020, Astrophys. J., 897, 26, doi: 10.3847/1538-4357/ab989e