# Resolving the Hot Dust Disk of ESO323-G77

James Leftley, Konrad Tristram, Sebastian Hönig, Daniel Asmus, Makoto

Kishimoto, Poshak Gandhi

► To cite this version:

James Leftley, Konrad Tristram, Sebastian Hönig, Daniel Asmus, Makoto Kishimoto, et al. Resolving the Hot Dust Disk of ESO323-G77. The Astrophysical Journal, 2021, 912 (2), pp.96. 10.3847/1538-4357/abee80. hal-03476857

HAL Id: hal-03476857

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

Submitted on 1 Aug 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Resolving the Hot Dust Disk of ESO323-G77

James H. Leftley<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, Konrad R. W. Tristram<fnref n="3" /> ID, Sebastian F. Hönig<fnref n="1" /> ID, Daniel Asmus<fnref n="1" />,<fnref n="4" /> ID, Makoto Kishimoto<fnref n="5" /> ID, and Poshak Gandhi<fnref n="1" /> ID

## Abstract

Infrared interferometry has led to a paradigm shift in our understanding of the dusty structure in the central parsecs of active galactic nuclei (AGNs). The dust is now thought to comprise a hot (~1000 K) equatorial disk, some of which is blown into a cooler (~300 K) polar dusty wind by radiation pressure. In this paper, we utilize the new near-IR interferometer GRAVITY on the Very Large Telescope Interferometer (VLTI) to study a Type 1.2 AGNs hosted in the nearby Seyfert galaxy ESO 323-G77. By modeling the squared visibility and closure phase, we find that the hot dust is equatorially extended, consistent with the idea of a disk, and shows signs of asymmetry in the same direction. Furthermore, the data is fully consistent with the hot dust size determined by K-band reverberation mapping as well as the predicted size from a CAT3D-WIND model created in previous work using the spectral energy distribution of ESO 323-G77 and observations in the mid-IR from VLTI/MID-infrared Interferometric instrument).

*Unified Astronomy Thesaurus concepts: Active galactic nuclei (16); Long baseline interferometry (932); Optical interferometry (1168); Astronomy data modeling (1859); Infrared astronomy (786); Near infrared astronomy (1093)*

## 1. Introduction

Infrared (IR) interferometry has played an integral part in the development of our understanding of active galactic nuclei (AGNs). Observations with the Mid-infrared Interferometer (MIDI; Leinert et al. 2003) at the Very Large Telescope Interferometer (VLTI) allowed the study of the central tens of parsecs in the mid-IR. At this scale, it was thought that a dusty toroidal structure, believed to be the source of obscuration required for the difference between Type 1 and Type 2 active galactic nucleus (AGN) in the unification scheme, would be found (e.g., Antonucci 1993). What was found instead was a polar dust structure (e.g., Tristram et al. 2007; Kishimoto et al. 2011b; Hönig et al. 2012, 2013; Burtscher et al. 2013; Tristram et al. 2014; Hönig et al. 2014; López-Gonzaga et al. 2016).

The discovery of the polar dust structure spurred the development of the disk-wind model (Emmering et al. 1992; Elitzur & Shlosman 2006; Hönig et al. 2013). In the model, the polar dust is an outflow that must have an origin where it receives material and energy. The source of the material in the disk-wind model is an equatorial hot dust disk that is close to the central engine of the AGN (Hönig & Kishimoto 2017). Such a disk would be on scales too small for MIDI to study in most AGN and is thought to be responsible for the unresolved emission observed with MIDI. The disk geometry was drawn from the warm dust disk seen with MIDI in the nearby AGN hosted by the Circinus Galaxy (Tristram et al. 2012, 2014).

The source of obscuration in the disk-wind scenario, necessary for AGN unification, is from the hot dust disk and the launch region of the dusty wind. The cool extended wind could also explain the apparent isotropic mid-IR appearance of AGN, irrespective of inferred line-of-sight obscuration (Gandhi et al. 2009; Asmus et al. 2015; Hönig & Kishimoto 2017). The hot dust should primarily emit in the near-IR, which makes the near-IR an

important regime to test the disk-wind model. So far, AGNs have been less studied with near-IR interferometry than with mid-IR interferometry. Weigelt et al. (2012) studied NGC 3783 using AMBER (Astronomical Multi-BEam combineR; Petrov et al. 2007); they reported a hot, 1400 K, thin dust ring fit. Kishimoto et al. (2009, 2011a) reported that near-IR interferometry was likely tracing the sublimation radius of the AGN using Keck observations, which agrees with the discovered thin ring/disk. The disk component can be explained in the disk-wind model as the hot dust disk, thought to be responsible for the bulk of the near-IR emission. However, the ring-like emission can also be explained in the classic torus model as the expected sublimation radius inside the torus. Further study of the near-IR dust distribution is needed to distinguish the two.

The new instrument GRAVITY at the VLTI presents us with the opportunity to constrain this hot dust disk. It operates in the near-IR, where such a disk would be brightest, and offers access to the spatial scales required. This had been attempted in a few objects already (GRAVITY Collaboration et al. 2020a, 2020b) and structure on the spatial scale of the hot dust disk has been detected. In this work we will attempt to constrain the hot (~1000 K) dust in the Type 1.2 nuclei of the local Seyfert galaxy ESO 323-G77. This object has shown a polar extension in the mid-IR as well as a particularly dominant unresolved component when studied with MIDI (Leftley et al. 2018). The unresolved component is thought to be the putative disk in this object, which makes it ideal to try and constrain the launch region of the dusty wind.

This paper is structured as follows. In Section 2, we detail the GRAVITY observations and the method of data reduction. Section 3 presents and discusses the results drawn directly from the observations as well as simplistic modeling of the squared visibility. Section 4 gives the models used in our attempt to

<footnote n="1">¹Department of Physics & Astronomy, University of Southampton, Southampton, SO17 1BJ, UK</footnote>

<footnote n="2">²Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Parc Valrose, Bat. H. Fizeau, 06108 Nice, France</footnote>

<footnote n="3">³European Southern Observatory, Alonso de Córdova 3107, Casilla 19001, Santiago, Chile</footnote>

<footnote n="4">⁴Gymnasium Schwarzenbek Europaschule, Buschkoppel 7, 21493 Schwarzenbek, Germany</footnote>

<footnote n="5">⁵Department of Astrophysics and Atmospheric Sciences, Kamigamo-Motoyama, Kita-ku, Kyoto, 603-8555, Japan</footnote>

<footnote>Received 2021 January 29; revised 2021 March 11; accepted 2021 March 12; published 2021 May 10</footnote>
<table><caption>**Table 1**
The GRAVITY Observations Used in This Paper and Their Total Observation Length</caption><thead><tr><th>Date</th><th>Time</th><th>Obs. Length (s)</th><th>V Group</th><th>C Group</th></tr></thead><tbody><tr><td>2019-04-23</td><td>03:04:24</td><td>300</td><td>...</td><td>...</td></tr><tr><td>2019-04-23</td><td>03:09:59</td><td>300</td><td>...</td><td>...</td></tr><tr><td>2019-04-23</td><td>03:25:07</td><td>300</td><td>...</td><td>...</td></tr><tr><td>2019-05-13</td><td>03:51:06</td><td>300</td><td>...</td><td>...</td></tr><tr><td>2019-05-13</td><td>04:00:33</td><td>300</td><td>...</td><td>...</td></tr><tr><td>2019-05-13</td><td>04:11:48</td><td>300</td><td>...</td><td>...</td></tr><tr><td>2019-06-20</td><td>23:38:04</td><td>160</td><td>1</td><td>1</td></tr><tr><td>2019-06-20</td><td>23:41:16</td><td>160</td><td>1</td><td>1</td></tr><tr><td>2019-06-20</td><td>23:47:47</td><td>160</td><td>2</td><td>2</td></tr><tr><td>2019-06-20</td><td>23:54:06</td><td>160</td><td>2</td><td>2</td></tr><tr><td>2019-06-21</td><td>00:00:40</td><td>160</td><td>3</td><td>3</td></tr><tr><td>2019-06-21</td><td>00:03:48</td><td>160</td><td>3</td><td>3</td></tr><tr><td>2020-03-15</td><td>07:31:13</td><td>160</td><td>4</td><td>4</td></tr><tr><td>2020-03-15</td><td>07:34:28</td><td>160</td><td>4</td><td>*</td></tr><tr><td>2020-03-15</td><td>07:41:07</td><td>160</td><td>4</td><td>4</td></tr><tr><td>2020-03-15</td><td>07:44:19</td><td>160</td><td>5</td><td>4</td></tr><tr><td>2020-03-15</td><td>07:47:31</td><td>160</td><td>5</td><td>5</td></tr><tr><td>2020-03-15</td><td>07:54:13</td><td>160</td><td>5</td><td>5</td></tr><tr><td>2020-03-15</td><td>07:57:28</td><td>160</td><td>6</td><td>5</td></tr><tr><td>2020-03-15</td><td>08:01:25</td><td>160</td><td>6</td><td>6</td></tr><tr><td>2020-03-15</td><td>08:04:37</td><td>160</td><td>*</td><td>*</td></tr><tr><td>2020-03-15</td><td>08:14:04</td><td>160</td><td>7</td><td>6</td></tr><tr><td>2020-03-15</td><td>08:17:31</td><td>160</td><td>7</td><td>6</td></tr></tbody></table>

**Note.** If a group is given, objects of the same group were combined when used. V Group is the grouping used when determining the visibility squared, C Group is the grouping used when determining the closure phase.\* This data set was not included in the analysis.

recreate the squared visibility and closure phase. Section 5 presents the results of the modeling and discusses the reliability and interpretation of these results.

## 2. Observations and Reduction

We obtained data from the GRAVITY instrument (GRAVITY Collaboration et al. 2017). In total, there was 8 hr of scheduled GRAVITY time split into four observations. We used a nonstandard observation strategy to compensate for any biases that were introduced from the adaptive optics (AO) systems on the VLTI. The bias from AO was discovered when observing faint objects with AMBER. It was found that the poor performance of the AO systems, on the UTs at the VLT, when observing R-band faint targets, such as AGNs, led to visibility losses (see also Burtscher et al. 2012 for the effect of poor AO correction in the mid-IR). In principle, this loss can be corrected with a faint red calibrator observed temporally and spatially close to the science target. The faint calibrator is a star that has similar R-band properties to the science target. The star would then be calibrated with the normal bright calibrator and any visibility loss present could be quantified corrected in the science.

The observations are listed in Table 1. Each observation block was structured RCAL-SCI-RCAL-CAL, where RCAL is the red calibrator. The SCI component of the 2019 April and 2019 May block is split into three 300 s observations, and the 2019 June block was instead split into six 160 s observations. Like the June block, the 2020 March block was due to be split into six parts; however, due to poor atmospheric conditions on the night, some observations were repeated and 11 total

<figure><img src="image_3.png" /><figcaption>**Figure 1.** The GRAVITY uv plane for ESO 323-G77.</figcaption></figure>

observations were obtained. Due to the poor conditions, observation 2020-03-15 08:04:37 was unusable and discarded. The UT3-UT2 and UT3-UT1 visibility squared of 2020-03-15 07:34:28 were discarded for the same reason as well as the closure phase for this observation. When performing the data reduction, all observations were reduced and calibrated individually. For the 2019 June night, observations were then paired together in time, and averaged into three observations. For the 2020 March night, the 10 sets of visibility squared, in time order, were collected into groups of three-three-two-two and averaged to produce four observations; the nine sets of closure phases were collected into three groups of three and averaged. The groupings are given in Table 1. For every observation, the CAL is HD 120271, the first RCAL is UCAC4 237-064201, and the second RCAL is NOMAD1 0480-0335664. The achieved uv plane for the observations of ESO 323-G77 can be found in Figure 1.

The nucleus of ESO 323-G77 is dim by the standards of the GRAVITY fringe tracker and there is no off-axis bright object with which to perform dual-field observations. When attempted in single field mode, 50% of the light was not enough to successfully track fringes. However, the fringe tracker data can be used for more than just tracking. Because we are interested in the continuum more than the line emission, we did not require high spectral resolution. The fringe tracker itself provides low spectral resolution ($R \sim 20$) data, albeit for short integration times. Therefore, we did not need the science channel at all. Instead, 100% of the target light can be sent to the fringe tracker using the dual-field off-axis mode. The target for the fringe tracker was set to be the nucleus of ESO 323-G77 and the science target for the science channel was set to be an arbitrary piece of the sky far enough away. Detector Integration Time (DIT) = 10 s for the science combiner was chosen for the 160 s exposures in order to allow for a full cycle of optical path difference offsets with the fringe tracker to de-bias the data (Lacour et al. 2019). DIT = 30 s for the science combiner was chosen for the 300 s exposures. For the data reduction and the analysis, the science channel was ignored.

### 2.1. GRAVITY Data Reduction

We make use of the ESOREX pipeline through the python tools provided by the Gravity Consortium to perform the initial
reduction. We calculate the final visibility and phase from the reduced products manually.

To reduce the data, we first use the run_gravi_reduce script from the python tools, which utilizes ESOREX with the gravity vis recipe. The tool produces intermediary products (given by setting the parameter -p2vmreduced-file = TRUE), which are a set of files that contain the reduced data before post processing as well as some extra relevant information such as the group delay. Most importantly for our reduction, the intermediary products contain the complex visibility for each frame, a rejection flag for each frame, and a geometric flux for each frame. We compute the visibility and the closure phase separately from this point.

2.1.1. Visibility Determination

To determine the visibility and visibility squared we initially used the default pipeline frame rejection without any smooth- ing over frames. This uses a signal-to-noise ratio cutoff of 3. We then compared the visibility of all calibrators to the Strehl ratio, calculated from the acquisition camera images by the pipeline using the -reduce-acq-cam option in the gravityvis ESOREX recipe. The Strehl ratio is our best proxy for the performance of the AO systems. Also considered was the header provided coherence time, the FWHM of a Gaussian fitted to the object as seen by acquisition camera, the core width and power index of a Moffat function fitted to the object, and the variance in central position of the source as determined by the fitted Moffat function, i.e., the “jitter,” in the acquisition camera over the course of the observation. All these measures correlated with loss in visibility and visibility squared from AO performance; the strongest correlation was with the Strehl ratio and the FWHM.

By modeling this loss, we could remove the majority of the visibility loss from the science exposures but a significant scatter, approximately 0.2 in visibility squared on the best night, still rendered most of the visibility data useless. Therefore, we derived a separate method for accounting for this loss.

We utilized the individual frame complex visibility and selected frames based on our own criteria. We attempted to use the group delay to select frames, as done in GRAVITY Collaboration et al. (2020a); however, we found that selecting by geometric flux provided the best correction in our data. We selected the 3% brightest frames in geometric flux, not including any frames rejected by the gravity vis recipe. All non-selected frames were then flagged and the visibility and visibility squared were calculated using the GRAVITY python tool run_gravi_reduce_from_p2vmred with the ESOREX recipe gravity vis from p2vmred using the flag - use-existing-rejection = TRUE. This removed some, but not all, of the AO loss. To remove the remaining bias, we calibrated the science data with the red calibrators. We show the visibility squared data with and without this cut in Figure 2. Both the science and calibrators were reduced by the same method.

It is suspected that the reason the AO loss can be removed this way is that when the AO is performing, well more flux is injected into the fiber. Therefore, by selecting the highest geometric fluxes, we select the frames where the AO is performing best or where the atmosphere was least turbulent. This method is reminiscent of “lucky imaging” sometimes used

for ground-based photometry (Fried 1978). The disadvantage of this method is we discard most of the data.

Finally, it is possible that the visibility squared at zero baseline length is not 1. This scaling issue is possibly caused by coherence loss and effects all baselines of the same observation equally. For a single night of observations, where the instrument setup remains constant, this scaling should be negligibly variant between similar uv positions (GRAVITY Collaboration et al. 2020a) and can be accounted for when modeling. However, between nights the scaling could change. This scaling does not appear to be completely removed through calibration and therefore may be caused by an instrumental response to dimmer K-band objects. The calibrators, including the red calibrators, are brighter in K than the science. To account for the scale variation, we take the median visibility squared value of all science observations during a night and scale it to the night of highest median visibility squared at a similar position angle. This can only be done for observations with similar position angles and baseline lengths because any scaling would be degenerate with angular dependent structure. This equates to the median visibility squared of the March and June nights being scaled to the median visibility squared of the April and May nights, respectively. While we cannot account for difference in scaling at different position angles, by normalizing to the highest visibility squared nights the different scaling between the different position angles is minimized. The remaining average scaling is then included in the models when fitting.

2.1.2. Closure Phase Determination

Accurate determination of the closure phase is important to determine the asymmetric detail in the geometry of ESO 323- G77. While it is piston invariant, there are still some residual offsets that need correcting after the initial pipeline reduction.

To calculate the closure from a single triplet of cotemporal frames in the p2vmred file, we used a similar method to the ESOREX pipeline. We calculate the complex bispectrum using

$$
\Psi_{ijk} = \psi_{ij} \cdot \psi_{jk} \cdot \psi_{ik}^*, \quad (1)
$$

where i, j, and k are three unique telescopes and ψ is the complex visibility of a telescope pair. We proceed to follow the method set by the fringe tracker (Lacour et al. 2019). When observing, the fringe tracker data are averaged over 300 frames when producing the closure phase estimates. Therefore, we split the data into 300 frame segments. In each 300 frame segment, we remove any trio of frames that have any of their three frames flagged with the default pipeline flagging. We do not make any flux cuts for the closure phase determination. We then calculate the mean bispectrum of the segment. We calculate the final closure phase from the angle of the mean of the bispectrum of every segment in the observation (Equation (2)) and the standard error from their distribution.

$$
\Phi_{ijk} = \arg \left( \frac{1}{n} \sum_{N=0}^{n} \frac{1}{300} \sum_{f=300N}^{299+300N} \Psi_{ijkf} \right). \quad (2)
$$

We calculate this for every science and calibrator observation and find that all calibrators have nonzero closure phase. The bright calibrator has a very large phase of ≈7° on every night. We do not know what causes this and more bright calibrators
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The visibility squared observations, at 2.15, 2.26, and 2.34 µm, using the default pipeline values and the 3% flux cut. The baseline lengths for each wavelength are adjusted to 2.15 µm equivalent. Both are calibrated using the same calibrators.</figcaption></figure>

would need to be checked to see if this is instrumental or if HD 120271 is a bad calibrator. The closure phase of the red calibrators on the same night are nonzero and approximately equal; furthermore, the common phase is clearly present in the science. As multiple red calibrators show the same phase component we can conclude this is instrumental. We therefore subtract the mean closure phase of the red calibrators from the science target. The remaining closure phase was taken to be truly from the science object.

### 3. Observational Results and Discussion

After reduction and calibration of the visibility squared (shown in Figure 3 for 2.15 µm) and closure phase, we can already draw some conclusions. There is a clear drop in visibility squared of approximately 0.1 from the 60–130 m baseline at all covered wavelengths except for the 1.99 µm bin. The shortest baseline (UT3-UT2) has a lower visibility squared than the second shortest baseline in the May, June, and March observations at all observed wavelengths. The lower squared visibility could suggest that the shortest baseline suffers an extra component of loss due to an unknown reason that was not mirrored in the calibrator. The suppressed first baseline is not present in the April observation, which is at a similar position angle to the June observation, suggesting it is not a real structure. Similar losses are seen in other objects in GRAVITY Collaboration et al. (2020a). Therefore, we do not consider the UT3-UT2 baseline when fitting squared visibility.

The first wavelength bin, 1.99 µm, is heavily suppressed in visibility and visibility squared by ~0.5 after the 3% flux cut. It also has a large, ~10°, scattered closure phase. Furthermore, both the visibility and closure phase seem disordered with baseline length when compared to all other wavelength bins. Without the 3% flux cut, the visibility and visibility squared are often negative. The reason is most likely the contamination from back-scattered light from the GRAVITY metrology laser at 1.908 µm (see GRAVITY Collaboration et al. 2017). This contamination has an intensity equivalent to a ~10 mag source, determined from the sky observations, for the setup we used, i.e., its contribution is comparable to the flux from our science target. We therefore do not consider results from the shortest wavelength bin to be reliable.

The second bin (2.06 µm) also has squared visibilities that are ~0.1 lower than the latter three wavelength bins. By itself it

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The visibility squared observations at 2.15 µm are in red. In purple is the best-fit model to the observations with 1σ error and in orange is the reverberation mapped radius with 1σ error. We only show visibility squared from 0.5–1 for clarity.</figcaption></figure>

can be considered reliable because it shows a similar structure to the latter bins but scaled and the scale of the squared visibilities is separately considered in fitting. As a precaution, we do not use the 2.06 µm bin either when performing multiwavelength fitting so as not to introduce any erroneous structure. Over the last three wavelength bins (2.15, 2.26, and 2.34 µm), the visibility squared has a standard deviation of at most 0.05 at similar baseline lengths. Therefore, the latter three bins can be used together in multiwavelength fitting.

The closure phase in each reliable wavelength bin shows structure. The phase for the 2.15, 2.26, and 2.34 µm wavelength bins is plotted in Figure 4 in purple. We find no phase for the larger triangles, as measured by the sum of the three baselines of the telescope triplet. The phase gradually decreases to ~-1.5° for smaller triangles. A small change in phase could be interpreted as a small shift of photocenter at larger spatial scales. The component responsible for the phase would have to be resolved out for the higher spatial resolution triplets.

#### 3.1. Comparison to Reverberation Mapping

As an initial test for structure, we make the assumption that the visibility squared in a single wavelength bin can be explained by a simple 1D Gaussian model. Additionally included in the model is a scale factor; this factor scales all the visibility squared measurements so that the zero baseline value may be different from one due to either visibility loss or an over-resolved structure. The model is as follows:

$$V^2(u, v, \lambda) = (1 - s_f) * F(u, v, \lambda, \Theta)^2, \quad (3)$$

where

$$F(u, v, \lambda, \Theta) = \exp \left[ - \left( \frac{C\Theta}{\lambda} \right)^2 (u^2 + v^2) \right] \quad (4)$$

and

$$C = \frac{\pi^2}{1.296 \times 10^9 \cdot \sqrt{\ln 2}}, \quad (5)$$

where Θ is the FWHM value of the Gaussian in milliarcseconds, λ is the wavelength, and s<sub>f</sub> is the scale factor. The constant C comes from the conversion of degrees to milliarcseconds and the conversion from σ to FWHM with an extra factor of π originating from the Fourier transform. For fitting we use EMCEE and the same method employed in the
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The closure phase for the 2.15, 2.26, and 2.34 µm wavelength bins are plotted in purple against the baseline length in Mλ. The red depicts the best-fitting model sampled at the same uv locations.</figcaption></figure>

geometric modeling of Leftley et al. (2019). We model three variables: $Θ$, $s_f$, and $f$, where $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). We do not employ the point-source fraction from their modeling since visibility squared shows no clear evidence of a constant. However, a point-source component is predicted due to contributions from the accretion disk; we do not see any sign of a turnover due to the point source because the extended component is not sufficiently resolved. To compare to the reverberation mapped radius from Boulderstone et al. (2021), we utilize the method used by the GRAVITY Collaboration et al. (2020a), in which a fitted Gaussian FWHM is converted to a ring radius assuming a fixed point-source contribution. While it is possible to directly fit a ring model, a Gaussian model allows us to keep our models consistent and comparable when searching for more complex geometry in latter modeling. Similar to the method used by the GRAVITY Collaboration et al. (2020a), when converting from Gaussian FWHM to ring size, we do include the expected point-source contribution. We use a point-source luminosity fraction of 0.15 instead of the 0.2 they reported. A point-source fraction of 0.15 is chosen here because the point-source fraction of 0.2 is a general value determined by Kishimoto et al. (2007) but a specific value for ESO 323-G77 was later determined by M. Kishimoto et al. (2021, in preparation).

From the fitting of the 2.15 µm bin, we find an FWHM of the Gaussian of $0.44_{-0.05}^{+0.05}$ mas (Figure 3). Using an angular size for ESO 323-G77 of $0.311 \text{ pc mas}^{-1}$, to match the distance used for the MIDI modeling in Leftley et al. (2018), this translates to a ring + point-source radius of $105_{-12}^{+11}$ lt-day. The result is within $1\sigma$ of the reverberation mapped radius of $89_{-18}^{+11}$ days, which also assumes a thin ring model. The 2.15 µm bin was

<table><caption>**Table 2**
The Size of the Hot Dust in ESO 323-G77 from a Thin Ring Model Fit to the Visibility Squared and the Hot Dust Size from Reverberation Mapping (Boulderstone et al. 2021)</caption><thead><tr><th>Wavelength (µm)</th><th>Ring Radius (lt-day)</th><th>FWHM (mas)</th></tr></thead><tbody><tr><td>Reverberation</td><td>89<sup>+11</sup><sub>-18</sub></td><td>0.37<sup>+0.05</sup><sub>-0.07</sub></td></tr><tr><td>2.06</td><td>116<sup>+11</sup><sub>-13</sub></td><td>0.48<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td>2.15</td><td>105<sup>+11</sup><sub>-12</sub></td><td>0.44<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td>2.26</td><td>121<sup>+6</sup><sub>-6</sub></td><td>0.5<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>2.34</td><td>80<sup>+10</sup><sub>-11</sub></td><td>0.33<sup>+0.04</sup><sub>-0.05</sub></td></tr><tr><td>Average</td><td>100<sup>+5</sup><sub>-5</sub></td><td>0.41<sup>+0.02</sup><sub>-0.02</sub></td></tr></tbody></table>

**Note.** The average is the mean of the visibility sizes.

chosen because the central wavelength best matches the reverberation mapping wavelength.

As a further test, we make the assumption that the regions probed by each wavelength bin in the GRAVITY data are the same within the angular resolution of GRAVITY. We fit all the wavelength bins, except the 1.99 µm bin, in the same manner as the 2.15 µm bin; the resulting ring radii are given in Table 2 and shown in Figure 5. In addition to the 2.15 µm bin, we find that the radial size agrees within $1\sigma$ of the reverberation mapped size for the 2.34 µm bin. The remaining two wavelength bins are instead within $2\sigma$ of the reverberation mapped size. Hence, there is a good agreement between both methods; although, the visibility size is slightly larger for three of four bins (see Figure 5). A larger interferometric size, as compared to the reverberation mapped size, is expected because interferometry gives the average size of the brightness distribution, while reverberation mapping is more sensitive to the inner edge of the hot dust disk, which is necessarily smaller (Hönig et al. 2014). If the 2.15 µm bin is assumed to be the best candidate for matching the reverberation mapped radii and the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The results of the circular Gaussian fits presented in Table 2. $r_{\text{ring}}$ on the left of the plot is the radius of the thin ring calculated from the fitted Gaussian FWHM on the right. The points in purple are the fits to the interferometric data, the orange area is the size measured from reverberation mapping.</figcaption></figure>

disk of ESO 323-G77 is well described by a thin ring, then a distance to ESO 323-G77 can be calculated by matching the two results. Comparing the sizes results in a distance of $54_{-11}^{+11}$ Mpc, which is within $2\sigma$ of the current luminosity distance of $70 \pm 5$ Mpc although the errors are large. A larger angular radii gives a smaller distance so the 2.06 and 2.26 µm agree less with current cosmology and the best agreement with current cosmology comes from the 2.34 µm bin. The calculation is done naively because we do not account for the difference in the region probed by the two techniques. Simultaneous modeling of the light curves used for reverberation mapping, and the interferometric data must be performed to determine the true source geometry and a more accurate distance (Hönig et al. 2014). Furthermore, we ignore the differences between the different distance measures in cosmology and the true geometry of the source which are negligible at such low redshift.

## 3.2. Angular Visibility Structure

We plot the uv plane for 2.15 µm with the points colored by visibility squared in Figure 6. From visual inspection, the plot suggests a dependence of visibility on position angle. As a preliminary test for position angle dependency, we attempt to fit the visibility squared data of the 2.06, 2.15, 2.26, and 2.34 µm bins individually with a centro-symmetric elongated Gaussian and point-source model. The point source is fixed to 0.15 to match the previous modeling and the method follows Leftley et al. (2018). Although the point source is fixed, the introduction of the scale factor means that the number of free parameters is the same as in their work. Fitting the latter four wavelength bins individually shows that a nonradially symmetric Gaussian is preferred for all but the 2.06 µm bin. The distribution of the fitted parameters for each of the fits can be found in Appendix A. The 2.06 µm bin prefers a radially symmetric Gaussian. The minor to major axis ratio of the elongated Gaussian ($\epsilon$) and the position angle of the major axis (PA) for the other bins can be found in Table 3.

The polar axis of the AGN system, i.e., the axis perpendicular to the plane of the accretion disk, is $174^{\circ} \pm 2^{\circ}$ or $155^{\circ} \pm 14^{\circ}$ east of north as defined by polarization measurements (Schmid et al. 2003; Smith et al. 2004; Batchelor et al. 2011) and the mid-IR extension from MIDI (Leftley et al. 2018), respectively. Therefore, the dust disk should have a position angle of $84^{\circ} \pm 2^{\circ}$ or $65^{\circ} \pm 14^{\circ}$, assuming there is no misalignment between the dust disk and the respective responsible medium. When an object is equatorially scattered,

the accretion disk is the scattering medium giving a very good descriptor of the accretion disk’s equatorial position angle (Smith et al. 2004). However, ESO 323-G77 is polar scattered, and the scattering medium in this case is thought to be dust clouds along the polar axis (Schmid et al. 2003). In the polar scattered case, the determined axis may not be a good descriptor of the precise polar axis. The MIDI extension traces the warm (~300 K) dust, the location of which is model dependent. In previous work, it has been found that the warm dust is generally polar extended, although often somewhat misaligned with the polarization determined polar axis (López-Gonzaga et al. 2016).

The results of the fitting show that the latter two wavelength bins prefer an extension that is closer to equatorial than polar, for both polar axis measurements, while the 2.15 µm bin is closer to polar by polarization and equatorial by MIDI; although, the 2.15 µm bin has relatively large uncertainty and is within $3\sigma$ of both directions. While the position angle of the 2.26 and 2.34 µm bins agree within $1\sigma$ of the MIDI equatorial direction, none of the results agree within $1\sigma$ uncertainty of the equatorial direction of $84^{\circ} \pm 2^{\circ}$ or the corresponding polar direction.

In an attempt to reduce the uncertainties, we perform a multiwavelength fit. This fit is performed using the 2.15, 2.26, and 2.34 µm bins under the assumption that the different wavelength bins probe the same dust with different angular resolution. The result of this fit is an elongated Gaussian with

<figure><img src="image_8.png" /><figcaption>**Figure 6.** The uv plane at 2.15 µm, the points are colored by their squared visibility.</figcaption></figure>

<table><caption>**Table 3**
The Results of the Elongated Gaussian and Point-source Fit</caption><thead><tr><th>Wavelength (µm)</th><th>ε</th><th>PA (°)</th><th>Θ<sub>y</sub></th></tr></thead><tbody><tr><td>2.06*</td><td>1</td><td>13<sup>+31</sup><sub>-27</sub></td><td>0.46<sup>+0.07</sup><sub>-0.11</sub></td></tr><tr><td>2.15</td><td>0.72<sup>+0.1</sup><sub>-0.11</sub></td><td>24<sup>+19</sup><sub>-22</sub></td><td>0.42<sup>+0.06</sup><sub>-0.09</sub></td></tr><tr><td>2.26</td><td>0.7<sup>+0.05</sup><sub>-0.04</sub></td><td>51<sup>+5</sup><sub>-5</sub></td><td>0.57<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>2.34</td><td>0.72<sup>+0.11</sup><sub>-0.11</sub></td><td>49<sup>+13</sup><sub>-23</sub></td><td>0.36<sup>+0.04</sup><sub>-0.06</sub></td></tr><tr><td>2.15 – 2.34</td><td>0.73<sup>+0.06</sup><sub>-0.05</sub></td><td>46<sup>+5</sup><sub>-9</sub></td><td>0.48<sup>+0.02</sup><sub>-0.02</sub></td></tr></tbody></table>

**Note.** $\epsilon$ is the minor to major axis ratio, PA is the position angle of the major axis, and $\Theta_y$ is the FWHM of the major axis*. While the object is not extended within errors, a position angle is provided for the cases where it is fitted with an extension.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** The overlapped elongated Gaussian fits to the squared visibility at 2.15, 2.26, and 2.34 µm compared to the mid-IR extension direction and the polar axis from polarization. The 2.15, 2.26, and 2.34 µm wavelength bins are plotted in blue, green, and red, respectively.</figcaption></figure>

an axis ratio of $0.73_{-0.05}^{+0.06}$ and a position angle of $46_{-9}^{+5}$. The result is inconsistent with both the polarization determined polar axis and equatorial axis (see Figure 7); however, it is closer to equatorial than polar. The position angle is consistent with the equatorial direction determined with MIDI.

## 4. Modeling

Simple models fit to the squared visibility gave us an overview of the dominant source structure; we will now proceed to model the visibility and closure phase simultaneously. Fitting visibility and phase will allow us to extract more information from the observations. The modeling in this chapter was performed in `python` and makes use of the packages: `ASTROPY` (Astropy Collaboration et al. 2013), `GALARIO` (Tazzari et al. 2018), and `EMCEE`. `GALARIO` allows us to perform Fourier transforms on images and derive complex visibility for a given set of uv points with little computational expense. Less computational expense makes expensive fitting methods such as Markov Chain Monte Carlo (MCMC) compatible with highly flexible models. For the MCMC method, we use `EMCEE` with flat priors for all given variables.

We attempt to explain the GRAVITY observations in the three longest wavelength bins simultaneously. Any attempt to fit the wavelength bins separately failed to constrain any reliable geometry. In a single wavelength bin, the data are too sparse for the relatively complex models, as compared to the ones used in Section 3. Therefore, we make the assumption that

<table><caption>**Table 4**
The Bounds for Each Free Parameter of the Point-source Model</caption><thead><tr><th>Parameter</th><th>Lower</th><th>Upper</th></tr></thead><tbody><tr><td>ΔR.A.</td><td>-10 mas</td><td>10 mas</td></tr><tr><td>Δdecl.</td><td>-10 mas</td><td>10 mas</td></tr><tr><td>p<sub>f</sub></td><td>0</td><td>0.85</td></tr></tbody></table>

<table><caption>**Table 5**
The Modeling Result for the Combination of the 2.15, 2.26, and 2.34 µm Wavelength Bins</caption><thead><tr><th>Component</th><th>p<sub>f</sub></th><th>ΔR.A.</th><th>Δdecl.</th></tr></thead><tbody><tr><td>1<sup>*</sup></td><td>0.15</td><td>0</td><td>0</td></tr><tr><td>2</td><td>0.84<sup>+0.001</sup><sub>-0.001</sub></td><td>-0.4<sup>+0.03</sup><sub>-0.04</sub></td><td>-0.36<sup>+0.07</sup><sub>-0.06</sub></td></tr><tr><td>3</td><td>0.01<sup>+0.001</sup><sub>-0.001</sub></td><td>3.7<sup>+0.06</sup><sub>-0.07</sub></td><td>1.88<sup>+0.07</sup><sub>-0.08</sub></td></tr></tbody></table>

Note. *This component was fixed.

the narrow range of wavelengths covered by the GRAVITY fringe tracker probes the same hot dust component at different baseline lengths, essentially allowing us to fill more of the uv plane. For simplicity, the models in this section are fitted in wavelength space.

We attempt to recreate the squared visibility and closure phase simultaneously using point sources with a method analogous to clean-like image reconstruction. To compare different numbers of point sources, the Bayesian information criteria (BIC) is used. The BIC is defined as

$$BIC = \ln(n) k - 2 \ln(L), \quad (6)$$

where *n* is the number of utilized data points, *k* is the number of parameters of the model, and *L* is the likelihood

$$2 \ln(L) = - \sum_{m}^{N} \left[ \frac{(y(x_m) - \text{model}(x_m, \alpha))^2}{\sigma_m^2 + f^2 \text{model}(x_m, \alpha)^2} + \ln(2\pi(\sigma_m^2 + f^2 \text{model}(x_m, \alpha)^2)) \right] \quad (7)$$

for a set of data of length *N* with values *y*, positions *x* and a model with parameters *α*. The variance is denoted by *σ*, and *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).

Each point source has a ΔR.A., Δdecl., and point-source fraction (*p*<sub>f</sub>). A point source in this case is not a delta function but a 2D Gaussian with an FWHM of a pixel, which prevents an image with finite resolution from producing a pixelated probability space. After the squared visibility is calculated, a scale factor (*s*<sub>f</sub>) is introduced in the same manner as in Equation (3). The sum of all point-source fractions is defined to be 1. A single point source was fixed to the central pixel with a fixed fraction of 0.15, the value used for the accretion disk for the reverberation mapping comparison. All position parameters are made relative to the central point source. The point-source model contained *k* = 3(*N*<sub>p</sub>−1) free parameters, not including *s*<sub>f</sub>, where *N*<sub>p</sub> is the number of points. When calculating the BIC the number of free parameters does not include *s*<sub>f</sub> because we are interested in the relative BIC, not the absolute, and it is a nonphysical component that is present in all models and only affects the visibility squared. The bounds for the free parameters are given in Table 4. When creating an image from the point-source model, the model is shifted so the central image pixel is the location of the mean of all point-source
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The visibility squared for the 2.15, 2.26, and 2.34 µm wavelength bins are plotted in purple against the triangle perimeter in Mλ. The red depicts the best-fitting model sampled at the same uv locations.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** The three point-source best-fit model to the GRAVITY observations. This image is not convolved with the beam size. The color bar represents normalized flux in log scale. The beam is given in white.</figcaption></figure>

positions weighted by point-source fraction to reduce boundary effects when performing the Fourier transform due to the finite image size. The bounds set for each parameter are as follows:

The model sampled from 2–32 point sources, and the initial walker positions were selected from a uniformly random distribution over the available parameter space. To prevent point-source stacking, the bound is set so that two point sources may not occupy the same pixel. When using the BIC, two point sources in the same pixel will always be a worse fit than one point source with the summed point-source fractions of both due to the use of additional unnecessary parameters. Therefore, it is inefficient to consider this solution. The pixel size is 0.0697 mas with an image size of 512 × 512 pixels. The

<figure><img src="image_12.png" /><figcaption>**Figure 10.** The 2.2 µm prediction of the visibility squared made by the CAT3D-WIND model created in Leftley et al. (2018) compared to the data. The model has a sublimation radius of 0.065 pc, which was used in their work.</figcaption></figure>

squared visibility and closure phase are then calculated from the image and compared to the data.

The MCMC analysis of each $N_p$ is run for 4000 iterations with 200 walkers. For the last 2000 iterations, or $4 \times 10^5$ samples, we calculate the BIC. All BICs from each $N_p$ are compared and the best combination of values are recorded. The lowest BIC parameter set are then used as the starting point for a final MCMC analysis by using 200 walkers scattered by a normal distribution of $\sigma = 10^{-2}$ centered on the lowest BIC parameter set for 4000 iterations. The last 2000 iterations are used to calculate the uncertainties and position of each parameter. We use the described method because the relatively sparse data coverage, as compared to the GRAVITY imaged AGN NGC 1068 (GRAVITY Collaboration et al. 2020b), requires an extremely large number of iterations to find the true result from a random starting position. There are also many “local minima” in the probability space preventing us from
<figure><img src="image_13.png" /><figcaption>**Figure 11.** The 2.2 µm CAT3D-WIND model created in Leftley et al. (2018) compared to the central point-source separation of the point-source model (white). The model has a sublimation radius of 0.065 pc, which was used in their work.</figcaption></figure>

finding a good starting position with a minimizer. Furthermore, it is too computationally expensive to calculate the BIC for every parameter combination. Therefore, the given method provided the best trade-off between accuracy and computational expense for the data.

Utilizing the same fitting method, we also employ a second set of more complex models. The second set of models are composed of Gaussians with variable FWHMs instead of the fixed FWHM point sources. Both elongated and radially symmetric Gaussians were attempted. These models are described in more detail in Section B in the Appendix. Both Gaussian models were compared to the point-source model using the BIC and provide no improvement with the current available data.

# 5. Modeling Results and Discussion

We find that the best-fitting model, out of the three different models presented, is the point-source model using three point-source components (Table 5). The point-source model provides a good description of the closure phase (Figure 4) and the squared visibility (Figure 8). The model, unconvolved with the beam size, is shown in Figure 9. We find that two of the sources, one of which is the fixed component, are separated by 0.5 mas and responsible for 99% of the total flux. While the separation is less than the Rayleigh criterion for a resolution of $2B_{max}/\lambda$, where $\lambda$ is the wavelength and $B_{max}$ is the maximum baseline length, this resolution denotes the separation at which the correlated flux of two point sources becomes zero. In reality, the resolution depends on the signal to noise of the observations and partially resolved structure can be constrained through modeling on scales less than the beam size. From the Gaussian modeling of the visibility squared, we determine that the physical size at the 0.5 mas scale can be robustly constrained. Therefore, we conclude the measured separation is reliable. The central pair lie along a position angle of 50°, which agrees more with an equatorial extension than a polar extension. This would suggest that there is structure on the 0.5 mas scale in an equatorial direction. The direction agrees with that found with the elongated Gaussian fits and the

separation agrees with the size from reverberation mapping. In previous work it has been shown that an elongated Gaussian, for the uv plane of ESO323-G77, can be constrained in any direction with the fitting method utilized (Leftley et al. 2018). Therefore, we conclude that the direction is reliable and the structure could correspond to the accretion disk and the dust disk or to either side of the dust disk's inner rim.

The third component is 4.2 mas away from the fixed component and is responsible for 1% of the flux. The faint component is essential for explaining the closure phase. It causes the small shift in photocenter required at lower spatial resolutions. It resides to the northeast of the group, which is, temptingly, ~90° from the mid-IR extension found in Leftley et al. (2018). However, it is on a larger scale than expected of the dusty disk and ESO 323-G77 is thought to be inclined at 60° from face-on, when modeled with CAT3D-WIND (Hönig & Kishimoto 2017; Leftley et al. 2018), which makes interpretation more difficult. Therefore, it would be misleading to claim that the faint component is in any way related to the dusty disk at this stage. Alternatively, the blob may not be true geometry but instead it could be indicative of an asymmetry or substructure in the equatorial direction. In a three point-source model, the distance to the faint point is well constrained; however, this does not rule out a more complex asymmetry on a different spatial scale. For example, a nonuniform distribution of dust clumps in the disk or a warm dust blob further from the sublimation radius. This would require further modeling to check.

## 5.1. Comparison of Results to Model Predictions

In terms of spatial distance, the central point sources are separated by ~0.17 pc using $0.311 \text{ pc mas}^{-1}$. The separation distance is remarkably close to the expected K-band size of the hot dust from the CAT3D-WIND model (Hönig & Kishimoto 2017) given in Leftley et al. (2018) of 0.15 pc. The size is predominantly defined by the drop of the squared visibility. The central separation is essentially what was measured when we calculated the ring size from the visibility squared in Section 3.1. With the limited uv coverage the two point sources would appear similar to a ring.

The CAT3D-WIND model that was created in Leftley et al. (2018) can be used to create what would be seen by a near-IR interferometer in the same manner as was performed for the mid-IR. The CAT3D-WIND prediction is shown in Figure 10. The CAT3D-WIND model is consistent with the data although an improvement to the fit can be found by scaling the model size from the predicted sublimation radius of 0.065–0.04 pc. The similarity between the predicted radius of 0.065 pc and half the point separation in the recovered point-source model supports the idea each point is either side of the hot dust disk’s inner rim (see Figure 11) but does not exclude the accretion disk and hot dust suggestion. While the point-source separation is larger than the sublimation radius of the CAT3D-WIND model, it is expected to be larger because we are using a point-source model to explain extended structure and will therefore recover the average brightness size not the innermost rim separation. The central group does also contribute to the closure phase, however, the dominant component is the interplay between the group and the distant faint source.
## 5.2. Gaussian Models

The multiple radially symmetric Gaussians model does reproduce the data, however, it is not a better fit than the point-source model with the same number of free parameters or fewer. It also does not provide a significantly better fit, as described by the BIC, using a greater number of free parameters. Therefore, the BIC is larger for the radially symmetric Gaussian model when compared to the best point-source model in all cases. In the elongated Gaussian model, the axis ratio for every component is unconstrained. Therefore, the point-source model is the best descriptor based on the currently available data.

We conclude that the central group of points represents structure on the scale of the putative sublimation radius and hot dust size predicted by the CAT3D-WIND model in Leftley et al. (2018). The central extension is 73° away from the MIDI determined polar axis and 54° away from the polarization axis, which suggests the extension is equatorial. This agrees with the elongation found from Section 3.2. If the far point, at ~1.3 pc from the fixed component, is real structure, it could be a hot cloud blown out from the sublimation radius, a large cooler clump, or a more distant hot spot heated by some other means. In the case of the cooler clump, the faint component should be prominent in the *L* or *M* band with MATISSE. It may also be indicative of more complex substructure in the disk not captured by the simple models. The data are currently too noisy and sparse to determine if the classical torus model or the dusty wind model is a better explanation of the observations. However, our results support the idea that the *K*-band hot dust is located in an equatorial dust structure consistent with the sublimation radius.

# 6. Summary

We present a study of ESO 323-G77 in the near-IR using the VLTI instrument GRAVITY. We constrain the geometry of the hot (~1000 K) dust through the use of simple squared visibility models and more complex squared visibility and closure phase models. Our main conclusions are:

1. The average radial size of the hot dust emission, as inferred by the visibility squared, is consistent with the size of the sublimation radius as inferred from *K*-band reverberation mapping.

2. The squared visibility shows angular structure at 2.15, 2.26, and 2.34 µm. The 2.26 and 2.34 µm extensions are more consistent with the equatorial direction than polar when compared to both the polar axis inferred by polarization and the MIDI extension from Leftley et al. (2018). The extension at 2.15 µm is consistent with the equatorial direction as inferred by the MIDI extension but not polarization.

3. The squared visibility and closure phase is best explained by three point sources, when using the BIC test, which are consistent with the interferometric size predicted by the CAT3D-WIND model from Leftley et al. (2018); Hönig & Kishimoto (2017), the sublimation radius size from reverberation mapping, and the expected direction of the disk.

We would like to thank the referee for the kind comments and suggestions that helped improve this work.

J.H.L., S.F.H., and D.J.W. acknowledge support from the Horizon 2020 ERC Starting Grant DUST-IN-THE-WIND (ERC-2015-StG-677117). J.H.L. acknowledges the support of the French government through the UCA JEDI investment in the Future project managed by the National Research Agency (ANR) under reference number ANR-15-IDEX-01. M.K. acknowledges support from JSPS under grant 16H05731. P.G. acknowledges support from STFC and a UGC/UKIERI Thematic Partnership. D.A. acknowledges funding through the European Union's Horizon 2020 and Innovation program under the Marie Sklodowska-Curie grant agreement number 793499 (DUSTDEVILS).

Based on European Southern Observatory (ESO) observing program 0103.B-0096(A).

*Software: Astropy (Astropy Collaboration et al. 2013), Corner (Foreman-Mackey 2016), Emcee (Foreman-Mackey et al. 2013), Galario (Tazzari et al. 2018), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020).*

# Appendix A
## Elongated Gaussian Modeling Statistics Plots

In Section 3.2 we fit the visibility squared with an elongated Gaussian model. Here we provide the corner plots of those fits in Figures 12–16.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Corner plot (Foreman-Mackey 2016) of the probability density function (PDF) of each parameter used in the elongated Gaussian fit in Section 3.2 for 2.06 µm.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Corner plot (Foreman-Mackey 2016) of the PDF of each parameter used in the elongated Gaussian fit in Section 3.2 for 2.15 µm.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Corner plot (Foreman-Mackey 2016) of the PDF of each parameter used in the elongated Gaussian fit in Section 3.2 for 2.26 µm.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 15. Corner plot (Foreman-Mackey 2016) of the PDF of each parameter used in the elongated Gaussian fit in Section 3.2 for 2.34 µm.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Corner plot (Foreman-Mackey 2016) of the PDF of each parameter used in the elongated Gaussian fit in Section 3.2 for 2.15–2.34 µm.</figcaption></figure>

## Appendix B
### Alternative Models

In Section 4, we discussed the point-source model and mentioned the use of more complex Gaussian models. Here, we describe the Gaussian models. There are two Gaussian models we considered, a radially symmetric Gaussian model and an elongated Gaussian model. The radially symmetric Gaussian model has a central Gaussian with a fixed position and brightness, similar to the point-source model; however, its FWHM was not fixed. Each non-fixed point has four free parameters: FWHM, ΔR.A. position, Δdecl. position, and relative log amplitude flux (defined as the log of a Gaussian component’s amplitude relative to the fixed point). Unlike the point-source model, the fixed amplitude was not set to 0.15 because it is the integrated flux of the Gaussian that would need

to be set to this value. The integration adds computational expense. To minimize computational expense, we instead set the log amplitude of the fixed Gaussian to 10 and set the stipulation that the fixed point was the highest in amplitude. Each log amplitude was then set relative to the fixed point and the image was normalized to a summed flux of 1 in linear space. The bounds are given in Table 6. The radially symmetric Gaussian model has $4N_G-3$ free parameters, where $N_G$ is the number of Gaussians and not including $s_f$, when calculating the BIC.

The elongated Gaussian model was identical to the radially symmetric Gaussian model except that the FWHM parameter was replaced with the major axis FWHM and two extra free parameters were introduced which are the major–minor axis ratio and the position angle of the major axis. The bounds are in Table 7. The major axis FWHM and the axis ratio were free
<table><caption>**Table 6**
The Bounds for Each Free Parameter of the Radially Symmetric Gaussian Model</caption><thead><tr><td>Parameter</td><td>Lower</td><td>Upper</td></tr></thead><tbody><tr><td>ΔR.A.</td><td>-10 mas</td><td>10 mas</td></tr><tr><td>Δdecl.</td><td>-10 mas</td><td>10 mas</td></tr><tr><td>A<sub>f</sub></td><td>0</td><td>10</td></tr><tr><td>FWHM</td><td>10<sup>-4</sup> mas</td><td>5 mas</td></tr></tbody></table>

Note. A<sub>f</sub> is the relative amplitude flux.

<table><caption>**Table 7**
The Bounds for Each Free Parameter of the Elongated Symmetric Gaussian Model</caption><thead><tr><td>Parameter</td><td>Lower</td><td>Upper</td></tr></thead><tbody><tr><td>ΔR.A.</td><td>-10 mas</td><td>10 mas</td></tr><tr><td>Δdecl.</td><td>-10 mas</td><td>10 mas</td></tr><tr><td>A<sub>f</sub></td><td>0</td><td>10</td></tr><tr><td>FWHM<sub>major</sub></td><td>10<sup>-4</sup> mas</td><td>5 mas</td></tr><tr><td>ε</td><td>0</td><td>1</td></tr><tr><td>PA</td><td>0°</td><td>180°</td></tr></tbody></table>

Note. A<sub>f</sub> is the relative amplitude flux, FWHM<sub>major</sub> is the FWHM of the major axis, and ε is the axis ratio.

parameters for the fixed Gaussian. The chosen setup gives the model 6N<sub>G</sub>-3 free parameters when calculating the BIC. Both Gaussian models sampled 2–16 Gaussians.

## ORCID iDs

James H. Leftley https://orcid.org/0000-0001-6009-1803

Konrad R. W. Tristram https://orcid.org/0000-0001-8281-5059

Sebastian F. Hönig https://orcid.org/0000-0002-6353-1111

Daniel Asmus https://orcid.org/0000-0003-0220-2063

Makoto Kishimoto https://orcid.org/0000-0002-2216-3252

Poshak Gandhi https://orcid.org/0000-0003-3105-2615

## References

Antonucci, R. 1993, ARA&A, 31, 473

Asmus, D., Gandhi, P., Hönig, S. F., Smette, A., & Duschl, W. J. 2015, MNRAS, 454, 766

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Batchelder, D., Robinson, A., Axon, D. J., et al. 2011, ApJ, 738, 90

Boulderstone, B. 2021, PhD Thesis, Univ. Southampton

Burtscher, L., Meisenheimer, K., Tristram, K. R. W., et al. 2013, A&A, 558, A149

Burtscher, L., Tristram, K. R. W., Jaffe, W. J., & Meisenheimer, K. 2012, Proc. SPIE, 8445, 84451G

Elitzur, M., & Shlosman, I. 2006, ApJL, 648, L101

Emmering, R. T., Blandford, R. D., & Shlosman, I. 1992, ApJ, 385, 460

Foreman-Mackey, D. 2016, JOSS, 1, 1864

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fried, D. L. 1978, JOSA, 68, 1651

Gandhi, P., Horst, H., Smette, A., et al. 2009, A&A, 502, 457

GRAVITY Collaboration, Abuter, R., Accardo, M., et al. 2017, A&A, 602, A94

GRAVITY Collaboration, Dexter, J., Shangguan, J., et al. 2020a, A&A, 635, A92

GRAVITY Collaboration, Pfuhl, O., Davies, R., et al. 2020b, A&A, 634, A1

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hönig, S. F., & Kishimoto, M. 2017, ApJL, 838, L20

Hönig, S. F., Kishimoto, M., Antonucci, R., et al. 2012, ApJ, 755, 149

Hönig, S. F., Kishimoto, M., Tristram, K. R., et al. 2013, ApJ, 771, 87

Hönig, S. F., Watson, D., Kishimoto, M., & Hjorth, J. 2014, Natur, 515, 528

Hunter, J. D. 2007, CSE, 9, 90

Kishimoto, M., Hönig, S. F., Antonucci, R., et al. 2009, A&A, 507, L57

Kishimoto, M., Hönig, S. F., Antonucci, R., et al. 2011a, A&A, 527, A121

Kishimoto, M., Hönig, S. F., Antonucci, R., et al. 2011b, A&A, 536, A78

Kishimoto, M., Hönig, S. F., Beckert, T., & Weigelt, G. 2007, A&A, 476, 713

Lacour, S., Dembet, R., Abuter, R., et al. 2019, A&A, 624, A99

Leftley, J. H., Hönig, S. F., Asmus, D., et al. 2019, ApJ, 886, 55

Leftley, J. H., Tristram, K. R. W., Hönig, S. F., et al. 2018, ApJ, 862, 17

Leinert, C., Graser, U., Przygodda, F., et al. 2003, Ap&SS, 286, 73

López-Gonzaga, N., Burtscher, L., Tristram, K. R. W., Meisenheimer, K., & Schartmann, M. 2016, A&A, 591, A47

Petrov, R. G., Malbet, F., Weigelt, G., et al. 2007, A&A, 464, 1

Schmid, H. M., Appenzeller, I., & Burch, U. 2003, A&A, 404, 505

Smith, J. E., Robinson, A., Alexander, D. M., et al. 2004, MNRAS, 350, 140

Tazzari, M., Beaujean, F., & Testi, L. 2018, MNRAS, 476, 4527

Tristram, K. R. W., Burtscher, L., Jaffe, W., et al. 2014, A&A, 563, A82

Tristram, K. R. W., Meisenheimer, K., Jaffe, W., et al. 2007, A&A, 474, 837

Tristram, K. R. W., Schartmann, M., Burtscher, L., et al. 2012, JPhCS, 372, 012035

Weigelt, G., Hofmann, K.-H., Kishimoto, M., et al. 2012, A&A, 541, L9