Resolving the Hot Dust Disk of ESO323-G77

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 recreate the squared visibility and closure phase. Section 5 presents the results of the modeling and discusses the reliability and interpretation of these results.

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 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.

Zoom In Zoom Out Reset image size

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 ∼ 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 smoothing 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 gravity_vis 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{ijk}}={\psi }_{{ij}}\cdot {\psi }_{{jk}}\cdot {\psi }_{{ik}}^{* },\end{eqnarray} \tag{ 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.

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{{ijk}}=\arg \left(\displaystyle \frac{1}{n}\displaystyle \sum _{N=0}^{n}\displaystyle \frac{1}{300}\displaystyle \sum _{f=300N}^{299+300N}{{\rm{\Psi }}}_{{ijkf}}\right).\end{eqnarray} \tag{ 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 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.

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.

Zoom In Zoom Out Reset image size

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 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.

Zoom In Zoom Out Reset image size

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:

$$\begin{eqnarray}&&{V}^{2}(u,v,\lambda )=(1-{s}_{f})\ast F{\left(u,v,\lambda ,{\rm{\Theta }}\right)}^{2},\end{eqnarray} \tag{ 3 }$$

where

$$\begin{eqnarray}&&F(u,v,\lambda ,{\rm{\Theta }})=\exp \left[-{\left(\displaystyle \frac{C{\rm{\Theta }}}{\lambda }\right)}^{2}({u}^{2}+{v}^{2})\right]\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&C=\displaystyle \frac{{\pi }^{2}}{1.296\times {10}^{9}\cdot \,\sqrt{\mathrm{ln}2}},\end{eqnarray} \tag{ 5 }$$

where Θ is the FWHM value of the Gaussian in milliarcseconds, λ is the wavelength, and s_f 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 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 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 pc mas⁻¹, 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σ of the reverberation mapped radius of ${89}_{-18}^{+11}$ days, which also assumes a thin ring model. The 2.15 μm bin was 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σ of the reverberation mapped size for the 2.34 μm bin. The remaining two wavelength bins are instead within 2σ 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 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σ of the current luminosity distance of 70 ± 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.

Zoom In Zoom Out Reset image size

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)

Wavelength	Ring Radius	FWHM
(μm)	(lt-day)	(mas)
Reverberation	${89}_{-18}^{+11}$	${0.37}_{-0.07}^{+0.05}$
2.06	${116}_{-13}^{+11}$	${0.48}_{-0.05}^{+0.05}$
2.15	${105}_{-12}^{+11}$	${0.44}_{-0.05}^{+0.05}$
2.26	${121}_{-6}^{+6}$	${0.5}_{-0.02}^{+0.02}$
2.34	${80}_{-11}^{+10}$	${0.33}_{-0.05}^{+0.04}$
Average	${100}_{-5}^{+5}$	${0.41}_{-0.02}^{+0.02}$

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

Download table as:  ASCIITypeset image

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 () and the position angle of the major axis (PA) for the other bins can be found in Table 3.

Zoom In Zoom Out Reset image size

Table 3. The Results of the Elongated Gaussian and Point-source Fit

Wavelength (μm)		PA (°)	Θy
2.06*	1	${13}_{-27}^{+31}$	${0.46}_{-0.11}^{+0.07}$
2.15	${0.72}_{-0.11}^{+0.1}$	${24}_{-22}^{+19}$	${0.42}_{-0.09}^{+0.06}$
2.26	${0.7}_{-0.04}^{+0.05}$	${51}_{-5}^{+3}$	${0.57}_{-0.02}^{+0.02}$
2.34	${0.72}_{-0.11}^{+0.11}$	${49}_{-23}^{+13}$	${0.36}_{-0.06}^{+0.04}$
2.15 − 2.34	${0.73}_{-0.05}^{+0.06}$	${46}_{-9}^{+5}$	${0.48}_{-0.02}^{+0.02}$

Note. is the minor to major axis ratio, PA is the position angle of the major axis, and Θ_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.

Download table as:  ASCIITypeset image

The polar axis of the AGN system, i.e., the axis perpendicular to the plane of the accretion disk, is 174° ± 2° or 155° ± 14° east of north as defined by polarization measurements (Schmid et al. 2003; Smith et al. 2004; Batcheldor 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° ± 2° or 65° ± 14°, 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σ of both directions. While the position angle of the 2.26 and 2.34 μm bins agree within 1σ of the MIDI equatorial direction, none of the results agree within 1σ uncertainty of the equatorial direction of 84° ± 2° 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 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.

Zoom In Zoom Out Reset image size

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 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

$$\begin{eqnarray}&&\mathrm{BIC}=\mathrm{ln}(n)\,k-2\,\mathrm{ln}(L),\end{eqnarray} \tag{ 6 }$$

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

$$\begin{eqnarray}&&\begin{array}{l}2\mathrm{ln}(L)=-\displaystyle \sum _{m}^{N}\left[\displaystyle \frac{{\left(y({x}_{m})-\mathrm{model}({x}_{m},\alpha )\right)}^{2}}{{\sigma }_{m}^{2}+{f}^{2}\mathrm{model}{\left({x}_{m},\alpha \right)}^{2}}\right.\\ \quad \left.+\ \mathrm{ln}\left(2\pi \left({\sigma }_{m}^{2}+{f}^{2}\mathrm{model}{\left({x}_{m},\alpha \right)}^{2}\right)\right)\right]\end{array}\end{eqnarray} \tag{ 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_f ). A point source in this case is not a delta function but a 2D Gaussian with am 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_f ) 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_p −1) free parameters, not including s_f , where N_p is the number of points. When calculating the BIC the number of free parameters does not include s_f 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 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 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 × 10⁵ 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 σ = 10⁻² 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 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.

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 $\tfrac{\lambda }{2{B}_{\max }}$, where λ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5. The Modeling Result for the Combination of the 2.15, 2.26, and 2.34 μm Wavelength Bins

Component	pf	ΔR.A.	Δdecl.
1*	0.15	0	0
2	${0.84}_{-0.001}^{+0.001}$	$-{0.4}_{-0.04}^{+0.03}$	$-{0.36}_{-0.06}^{+0.07}$
3	${0.01}_{-0.001}^{+0.001}$	${3.7}_{-0.07}^{+0.06}$	${1.88}_{-0.08}^{+0.07}$

Note. *This component was fixed.

Download table as:  ASCIITypeset image

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 pc mas⁻¹. 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 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.2. Gaussian Models

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).

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 parameters for the fixed Gaussian. The chosen setup gives the model 6N_G −3 free parameters when calculating the BIC. Both Gaussian models sampled 2–16 Gaussians.

Table 7. The Bounds for Each Free Parameter of the Elongated Symmetric Gaussian Model

Parameter	Lower	Upper
ΔR.A.	−10 mas	10 mas
Δdecl.	−10 mas	10 mas
Af	0	10
FWHMmajor	10−4 mas	5 mas
	0	1
PA	0°	180°

Note. A_f is the relative amplitude flux, FWHM_major is the FWHM of the major axis, and is the axis ratio.

Download table as:  ASCIITypeset image