The disk of FU Orionis viewed with MATISSE/VLTI (cid:63) First interferometric observations in L and M bands.

Aims. We studied the accretion disk of the archetypal eruptive young star FU Orionis with the use of mid-infrared interferometry, which enabled us to resolve the innermost regions of the disk down to a spatial resolution of 3 milliarcseconds (mas) in the L band, that is, within 1 au of the protostar. Methods. We used the interferometric instrument MATISSE / VLTI to obtain observations of FU Ori’s disk in the L , M , and N bands with multiple baseline conﬁgurations. We also obtained contemporaneous photometry in the optical ( UBVRIr (cid:48) i (cid:48) ; SAAO and Konkoly Observatory) and near-infrared ( JHK s ; NOT). Our results were compared with radiative transfer simulations modeled by radmc -3 d . Results. The disk of FU Orionis is marginally resolved with MATISSE, suggesting that the region emitting in the thermal infrared is rather compact. An upper limit of ∼ 1 . 3 ± 0 . 1 mas (in L ) can be given for the diameter of the disk region probed in the L band, corresponding to 0.5 au at the adopted Gaia EDR3 distance. This represents the hot, gaseous region of the accretion disk. The N -band data indicate that the dusty passive disk is silicate-rich. Only the innermost region of said dusty disk is found to emit strongly in the N band, and it is resolved at an angular size of ∼ 5 mas, which translates to a diameter of about 2 au. The observations therefore place stringent constraints for the outer radius of the inner accretion disk. Dust radiative transfer simulations with radmc -3 d provide adequate ﬁts to the spectral energy distribution from the optical to the submillimeter and to the interferometric observables when opting for an accretion rate ˙ M ∼ 2 × 10 − 5 M (cid:12) yr − 1 and assuming M ∗ = 0 . 6 M (cid:12) . Most importantly, the hot inner accretion disk’s outer radius can be ﬁxed at 0.3 au. The outer radius of the dusty disk is placed at 100 au, based on constraints from scattered-light images in the literature. The dust mass contained in the disk is 2 . 4 × 10 − 4 M (cid:12) , and for a typical gas-to-dust ratio of 100, the total mass in the disk is approximately 0.02 M (cid:12) . We did not ﬁnd any evidence for a nearby companion in the current interferometric data, and we tentatively explored the case of disk misalignment. For the latter, our modeling results suggest that the disk orientation is similar to that found in previous imaging studies by ALMA. Should there be an asymmetry in the very compact, inner accretion disk, this might be resolved at even smaller spatial scales ( ≤ 1 mas).


Introduction
FU Orionis is the archetype of the FUor class of young stellar objects (YSOs), which experience eruptive events initiated by increased accretion (order of ∼ 10 −4 M yr −1 ) of material from their circumstellar disks onto the stellar surface (Herbig 1977;Hartmann & Kenyon 1996). These eruptions subsequently heat the material near the protostar (i.e., the inner accretion disk) and result in a brief rise in brightness with a subsequent decay. For the case of FU Orionis, the eruption increased its brightness by ∆B ∼ 6 magnitudes within a single year (1936-37;Hoffleit 1939;Herbig 1966;Kenyon et al. 1988). According to Kenyon et al. (2000), the B-band brightness has been declining slowly ever since by approximately 0.015 mag per year. Not more than 20 other stars have been found to belong to the FUor class thus far (e.g., Connelley & Reipurth 2018). It is unclear whether FUors are an entirely unique class of YSOs or a common step in the evolutionary path for all YSOs (e.g., Hartmann 2009). FU Ori is located at a distance of 402.3 +3.0 −3.7 pc according to Gaia Early Data Release 3 (EDR3; Bailer-Jones et al. 2021). It is thought to be a young low-mass object based on theoretical estimates (e.g., 0.6 M ; Pérez et al. 2020;Zhu et al. 2020). The star is surrounded by a circumstellar disk (Hartmann & Kenyon 1985) and, on a larger scale, by a reflection nebulosity (Wachmann 1954;Herbig 1966). The properties of the circumstellar accretion disk around FU Orionis, such as its orientation, accretion rate, and composition, have been extensively studied in the past. The majority of earlier works (e.g., Calvet et al. 1991;Hartmann & Kenyon 1996;Henning et al. 1998;Kenyon et al. 2000;Zhu et al. 2007) provided parametric studies of the disk based on fitting the spectral energy distribution (SED) and spectroscopic observations with accretion disk models; this continued until the arrival of high-spatial-resolution observations, which allowed the disk to be directly detected.
Scattered-light images in the near-infrared reveal an arcshaped stream of material appearing at about 0.3 directly east of FU Ori, with additional gaps in the northern and western directions and a tentative suggestion of a jet ejection (Liu et al. 2016;Takami et al. 2018;Laws et al. 2020). Furthermore, Pérez et al. (2020) suggested that the gaseous component of FU Ori's disk is in fact in Keplerian rotation and possibly related to the arc-shaped stream seen in scattered light.
FU Orionis is a member of a wide binary system. The companion star, FU Ori S, was located approximately ∼0.5 away at a position angle (PA) of ∼ 161 • in 2002 (Wang et al. 2004). Its presence has been confirmed by many follow-up imaging studies (e.g., Reipurth & Aspin 2004), and it is believed to be a premain-sequence K star. Far-millimeter continuum emission suggests that the companion also hosts its own circumstellar disk, and that both disks are aligned in similar directions Liu et al. 2017Liu et al. , 2019Pérez et al. 2020).
The introduction of infrared interferometry since the late 1990s has allowed an even closer inspection of the disk, as multiwavelength studies have found the apparent size of the accretion disk's innermost hot region to be less than 2 au (Malbet et al. 1998(Malbet et al. , 2005Quanz et al. 2006;Monnier et al. 2009;Liu et al. 2019;Labdon et al. 2021). These earlier studies explored the case of a third companion near FU Orionis, although a concrete detection is yet to be made. More recently, Liu et al. (2019) and Labdon et al. (2021) measured near-zero closure phases, which indicate an inclined, centro-symmetric distribution and would exclude the presence of a nearby companion. Nevertheless, Labdon et al. (2021) proposed an upper limit of about 1.3% for the H-band flux ratio for a subsolar companion with a separation between 0.5 and 50 milliarcseconds (mas), or between 0.5 and 20 au for our adopted Gaia distance.
Nearly all interferometric studies find different geometric properties, that is, inclination and PA, for the FU Orionis disk. Since these measurements were inferred at different wavelengths, from the near-infrared to the radio regimes, it is worth questioning whether individual components of the circumstellar disk may be misaligned. The most precise measurements thus far were made via direct imaging of the disk in submillimeter wavelengths . Since the geometric properties are a crucial input in models, the results of disk simulations (such as inferring the stellar radius and mass, the disk's inner and outer radius, and the mass accretion rate) can often vary. For the mass accretion rate in particular, such variations are found to be on the order of ∼ 10 −4 M yr −1 . For instance, the boundary between the hot inner disk and the passive, cooler component is found at 0.76 ± 0.35 au with the outer disk radius fixed at 7.7 au by Labdon et al. (2021), assuming a distance of 416 pc; this may be in agreement with Liu et al. (2019), suggesting that submillimeter and centimeter emission originates within 10 au radii. On the other hand, Zhu et al. (2007) place the outer radius of the hot, inner disk at ∼ 0.5 au (at a distance of 500 pc), with a not well-confined truncation radius of ∼ 1 au between the inner and the dusty disk, as it is not clear which disk component dominates in emission in the 4 -8 µm region (see also Zhu et al. 2008). Recently, Liu et al. (2021) found that the disk's mid-plane is dominated by millimeter-sized grains and posited that a temporal change in the source's brightness may be the result of dynamical changes in the disk.
In this work we present the first observations of FU Orionis with the Multi AperTure mid-Infrared SpectroScopic Experiment (MATISSE) instrument of the Very Large Telescope Interferometer (VLTI) obtained in three bands, namely L, M (3 -5 µm), and N (8 -13 µm) (Lopez et al. 2014;Matter et al. 2016a,b;Petrov et al. 2020;Lopez et al. 2022). They were complemented by nearly simultaneous photometric observations in the optical and near-infrared from the South African Astronomical Observatory (SAAO), the Nordic Optical Telescope (NOT), and the Konkoly Observatory. Our findings are presented in Sects. 2 and 3. The interferometric results are compared with analytical and radiative transfer models for the inner and outer parts of the accretion disk, respectively, and can be found in Sect. 5, which is followed by a discussion on our findings in Sect. 6. Our conclusions are described in the final section.

MATISSE/VLTI long-baseline interferometry
Guaranteed Time Observations (GTO) were obtained at the European Southern Observatory with the 4-beam recombiner instrument MATISSE (Lopez et al. 2014;Matter et al. 2016a,b;Petrov et al. 2020;Lopez et al. 2022) on the VLTI over five epochs (Table B.1) with different baseline configurations with the 8.2 m unit telescopes (UTs) and the 1.8 m auxiliary telescopes (ATs). These configurations produce baseline coverage between 40 and 130 meters, with angular resolutions of ∼ 3 -9 mas in L and ∼ 8 -26 mas in N, respectively. Hereafter, the term "baseline" will always refer to the "projected" baselines (length and/or PA) on the sky.
Here, we present results from two epochs (i.e., epochs 4 and 5, Table B.1) for which we obtained simultaneous photometric observations (see the following sections). A brief description of the data quality of all epochs is given in Appendix B.
We opted for low-spectral resolution (R ∼ 30) in L and M and for high-spectral resolution (R ∼ 220) in the N band in epoch 4 with the UTs. These observations were obtained in the standard hybrid mode, which utilizes the simultaneous photometry (SiPhot) mode in the L band and the high sensitivity (High-Sens) mode in the N band. In the former, photometry is recorded simultaneously to the interferometric fringes in a separate channel after splitting the incoming beam, while in the latter mode the photometry is collected after measuring the interferometric fringes. A similar spectral setup was chosen for epoch 5 (ATs), that is, low-spectral resolution in LM and high in N. However, these data sets were obtained in the new GRAVITY for MA-TISSE (GRA4MAT) mode (Lopez et al. 2022), which is designated for faint targets. In this mode, the GRAVITY/VLTI instrument is used as a fringe-tracker, thereby providing improved sensitivities for the ATs and increasing the spectral coverage. This was essential for FU Ori since we found that, for such a faint target, observations in the standard hybrid mode and the ATs were very difficult (cf. Appendix B).
The observations benefited from excellent atmospheric conditions with an average seeing as low as 0.5 (cf.  K2II). The data sets were reduced with the standard MATISSE data reduction software (DRS) pipeline (versions 1.5.2 and 1.6) and designated MATISSE python tools 1 . Since FU Ori is quite faint in all bands, that is, below the minimum suggested from MATISSE specifications 2 , when reducing the N-band UT data we opted to use a larger spectral bin (16 pixels) as opposed to the standard method (7 pixels). This averaging minimizes the noise in the correlated fluxes and in the closure phases, in order to reveal any deviations from spherical symmetry. A detailed description of the data calibration for MATISSE has been presented in previous works (e.g., Varga et al. 2021); therefore, we refrain from repeating it here. The interferometric products are presented in Sect. 3.

Optical and infrared photometry
With the goal to construct an SED close in time to MATISSE interferometric observations, we collected optical and nearinfrared photometry in the period from January 2021 to February 2021. A description of these observations follows below.

Nordic Optical Telescope
Photometric observations were obtained with the NOT nearinfrared camera and spectrograph (NOTCam; Abbott et al. 2000) on 25 January and 18 February 2021 (program 62-410, PI: F. Lykou), with the high-resolution (HR) imaging mode in standard J, H, and K s filters. The detector array (1024 × 1024 pixels) offers a field-of-view of 80 ×80 in HR mode with a plate scale of 0.078 /pixel. To avoid saturation due to the brightness of FU Ori in the near-infrared, a 5 mm pupil mask was intro-duced that reduced the transmission down to 10%. The average seeing on both nights was ≤ 1.2 .
FU Orionis exposures were bracketed by two comparison stars (2MASS J05452885+0901452, 2MASS J05451389+0904443) to assist in photometric calibration. In the first epoch, all three stars were exposed for 4.5 sec per dither in a five-point dithering pattern (read-read-reset mode) in all three filters. Due to final low signal-to-noise of the faint calibrator on the first epoch (January), we repeated the observations in February using the same camera setup but with a nine-point dithering pattern and 3.6 sec per dither with the ramp-sampling readout mode.
A median sky frame was subtracted from each science frame, and a differential flat-field was used to correct for pixel-to-pixel variations. Photometry was computed in the Two Micron All Sky Survey (2MASS) system for each dithering position against both comparison stars in varying aperture radii between 10 and 50 pixels, and it stabilized at 32 pixels (∼ 2.5 ). FU Ori's brightness was found to be consistent between the two epochs at 6.90, 6.21, and 5.65 mag in J, H, and K s, respectively, with a typical uncertainty δm ∼ 0.03 mag. The NOTCam photometry is included in Table A.1.

South African Astronomical Observatory
Optical photometric measurements were obtained at SAAO simultaneously with the NOT observations (25 January 2021; P.I.: M. Siwak, program ID: Siwak-2020-05-40-inch-317). We utilized the Sutherland Highspeed Optical Cameras (SHOC; Coppejans et al. 2013) with Bessel U BVRI filters, mounted on the 1-meter Lesedi telescope. The instrument offers a 5.72 ×5.72 field-of-view with a 0.335 /pixel plate scale; however, 2x2 binning was used to match to the seeing conditions. The secondary standards GSC 00714-00203 and GSC 00715-00188 were used as calibrators (Siwak et al. 2018). The photometry is tabulated in Table A.1, where U BV are in the Bessel system, and R c I c in the Cousins system.

Konkoly Observatory
Photometry was obtained with the fully automated Astro Systeme Austria (ASA) AZ800 alt-azimuth, direct-drive, 80centimeter Ritchey-Chrtien (RC) telescope at the Piszkéstető Observatory, Hungary, on 13-18 February 2021. We used a 2048×2048 pixel FLI PL230 CCD camera with the E2V CCD230-42 detector. The optical setup with an effective focal distance of f = 5600 mm yields a pixel scale of 0.55 and a field-of-view of 18.8 ×18.8 . We obtained three images per night with Bessel B, V and Sloan r and i filters.
A standard data reduction method was followed by applying bias, dark, and flat-field corrections. Aperture photometry was obtained for both the science target (FU Ori) and about 30 nearby, comparison stars using an aperture radius of 10 pixels (5.5 ) and a sky annulus of 10-15 pixels (5.5 -8.25 ). Photometric calibration was done by fitting a color term using the magnitudes and colors of the comparison stars from the AAVSO Photometric All-Sky Survey (APASS) DR9 catalog (Henden et al. 2016). The resulting photometry is presented in Table A.1, where BV are in the Bessel system and r i in the Sloan system. The typical uncertainty of the photometry is 0.01 mag.  (Table B.1). NEOWISE photometry from 2020 (corrected for saturation) is shown for comparison (green). Bottom: N-band spectrum (black; UTs) against earlier spectra from Spitzer (green; Green et al. 2006), MIDI (red; Quanz et al. 2006), and SOFIA (blue; Green et al. 2016). The flux uncertainty of MATISSE is about 10% for the case of FU Ori (error bar in the lower panel). Spectral regions where atmospheric transmittance hinders terrestrial observations are cut off from the spectra.

Interstellar extinction
The neutral hydrogen column density, N H , in a region within 0.1 • of FU Ori's location is 1.99 × 10 21 atoms/cm 2 based on HI 4-PI Survey (HI4PI) data (HI4PI Collaboration et al. 2016). Adopting the calculation of Güver & Özel (2009) where N H = (2.21 ± 0.09) × 10 21 A V , we find a low interstellar extinction of A V = 0.90 ± 0.04 mag. This should account for foreground extinction within a region of 0.1 • ; however, FU Ori is located inside the dust cloud B35 and therefore the extinction is suspected to be higher.
Conflicting values have been found in the literature with regard to the measurement of the interstellar extinction. Previous estimates were based on fitting FU Ori's broadband photometry and/or spectra, and placed A V between 1.5 (Zhu et al. 2007) and 3.2 mag (Herbig 1966). Here, we adopt the value derived from independent methods, and in particular from the 3D Dust Maps of Green et al. (2019) for the adopted Gaia distance of 402.3 +3.0 −3.7 pc (Bailer-Jones et al. 2021), that is, A V = 1.7 ± 0.1 mag.

Total flux
We draw caution on the absolute flux calibration of MATISSE in all three bands. On the one hand, the performance is inherent to the source's brightness itself (cf. Sect. 2.1), on the other hand it is dependent on atmospheric conditions. The latter can be assessed through the variability of the transfer function. It should be also optimized when using chopping in the L band to remove the thermal background, as this is affected more by the limited terrestrial atmospheric transmission between 3 and 5 µm. Figure 1 (top panel) shows the total flux spectra in LM from the UTs (blue) and from the ATs (red). A description of the data quality is given in Appendix B. Overall, the flux levels agree quite well despite being obtained in different modes and baseline configurations, while they are directly comparable to Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE) photometry obtained almost a year earlier (corrected for saturation; green points).
The N-band total flux spectrum is shown in the bottom panel of Fig. 1 (black line). We find that the flux uncertainty of MA-TISSE with the UTs for epoch 4 is about 10% for this faint target. This is mostly due to calibration errors and it can become worse beyond 11µm. Despite the higher sensitivity offered with the UTs on MATISSE, the S/N decreases beyond 11 µm for FU Ori (cf. Appendix B); therefore, all interferometric products will be noisier at the tail end of the N band. No photometry was obtained in epoch 5.
When compared with earlier spectra from satellite, airborne, and terrestrial telescopes -Spitzer, (green line; Green et al. 2006), the Stratospheric Observatory for Infrared Astronomy (SOFIA; blue line; Green et al. 2016), and the MID-infrared Interferometric (MIDI) instrument (red line; Quanz et al. 2006), respectively -we note that the total flux appears to have dimmed by at least 0.5 Jy since 2004. However, due to the flux calibration uncertainties from all instruments, this is not clear. Further discussion on this dimming can be found in Sect. 6.1.

L and M bands
The squared visibilities, V 2 , and correlated fluxes, F corr , in L and M bands are shown in the left and middle columns in Figs. 2 and 3. The top panels are the data from the UTs, and the bottom panels are the data from the ATs. The influence of the terrestrial atmosphere is evident in the abrupt changes in the L-band V 2 and F corr below 3.2 µm and beyond 3.9 µm, that is, at the edges of this atmospheric transmission window (see also Appendix B for a further description on data quality). Despite the different observational setups between the two epochs, we find that the results are quite similar in both bands.
In both the L and M bands, the squared visibilities approach unity at the shorter baselines (i.e., B ≤ 90 m or else spatial scales > 4 mas) and thus indicate that FU Ori's inner disk region is almost unresolved at those spatial scales. However, that region is marginally resolved at the longest baselines (spatial scales < 4 mas) where 0.65 ≤ V 2 ≤ 0.80. This suggests that the emitting region in L and M (i.e., the accretion disk) is rather compact and confined within an area with a global size ≤ 1.6 au; therefore, its radius should be smaller than 0.8 au. The correlated fluxes show a similar behavior to the squared visibilities with minute variations between baselines. Overall, the correlated spectra are flat and do not show any of the typical spectral features of FUors (e.g., Brα, CO), but this should be expected in this case where low-spectral-resolution setups were used.

N band
The N-band data from epoch 4 are shown in Fig. 4. FU Ori is known to have the 10 µm silicate feature in emission (e.g., Green et al. 2006;Quanz et al. 2006) and that is indeed seen in the MATISSE correlated spectra and the visibilities. Although the sensitivity of the instrument was increased with the introduction of GRA4MAT in epoch 5, the N-band brightness of FU Orionis is still below the GRA4MAT sensitivity limit; therefore, the observations suffered from very low S/N and are not presented here.
The visibilities also indicate that a part of the dusty disk is resolved by MATISSE more at the longest baselines (i.e., 90 ≤ B ≤ 120 m range). This suggests that the innermost region of the dusty disk is located inside an area of 9 mas in size (at 10 µm), otherwise within a radius of 2 au from the protostar. However, that area accounts for less than 50% of the 10 µm emission since the average visibility is ≤ 0.5. The majority of the 10 µm silicate dust emission arises from a larger region of the disk, which can be probed by MATISSE with the shortest UT baselines (B ∼ 40 m). Its size ought to be smaller than roughly 25 mas at 10 µm, corresponding to an area smaller than 10 au. The physical size of the disk is expected to be much larger, perhaps a few hundred astronomical units, as hinted by nearinfrared scattered-light images. However, the extended disk is not seen here by MATISSE because either (a) it could only be probed at shorter baseline configurations (e.g., B ≤ 20 m), or (b) the extended disk is not a strong emitter in the mid-infrared. Overall, the N-band visibilities measured by MATISSE are consistent with those measured by MIDI (Quanz et al. 2006).

Closure phases
The MATISSE closure phases in L, M, and N are shown in the rightmost columns in Figs. 2, 3, and 4. In the L and M bands, the closure phases are consistent between epochs and observational setups, and remain constant within each band with an average value of 1.5 • . One exception are the M-band data of epoch 5 (ATs; GRA4MAT), which are noisy overall. Nevertheless, this average value is below the expected performance of MATISSE of 5 • (Petrov et al. 2020). Since FU Ori's disk is marginally resolved in two baselines, we refrain from interpreting this closure phase signal as anything but a quasi-resolved point source.
The N-band closure phase signal is also nonzero with some small variations within the band, and it approaches 5 • at about 10 µm at the largest baseline triangle, which includes the longest UT baseline (green line, top row in Fig. 4). This nonzero closure phase signal could point at an asymmetric brightness distribution. As we show further on, a disk inclined in our line of sight can reproduce such a signal (Sect. 6.2).

Geometric sizes
To estimate the geometric properties of the accretion disk from the MATISSE observations, simple centro-symmetric models were fitted to the interferometric data to obtain the angular size and orientation of the detected structure. The Gaussian full width at half maximum (FWHM) for the L-band data gave a "marginally resolved" angular size of 1.3 ± 0.1 mas at 3.5 µm, which suggests an apparent size of about 0.5 au for the adopted Gaia distance (Fig. 5).
For the the N band, we opted to model the disk with 2D elliptical Gaussian functions following the method of Varga et al. (2021). The best-fit to the data is provided with a Gaussian FWHM of 5.3 ± 1.4 mas at 10.5 µm, a PA of 15 ± 25 • for the minor axis, and an inclination, i, of 55 ± 15 • .
As mentioned in Sect. 3.2, FU Ori's disk is marginally resolved at the longest baselines in the L band, while in the N band MATISSE was able to resolve a portion of the dusty, passive disk. Therefore, from the geometric sizes estimated above, we can place upper limits on the apparent outer radii of the hot, accretion disk of ≤ 0.3 au in L and 1 au in N.
To assist the reader, throughout the text we accept the conventional orientations for the minor axis PA and disk inclination angle. That is, the PA is measured east-of-north starting at zero degrees north, while i = 90 • in our line of sight describes an edge-on structure, and a pole-on structure is found at i = 0 • . As such, wherever necessary for the context of this work, literature values were adapted to fit these conventions.

Analytical and radiative transfer modeling
To further constrain other disk parameters, such as accretion rate, flaring, density distribution, grain sizes, and chemical composition, we opted to apply analytical and radiative transfer models. A description of this process follows below.
In this analysis, we explore two different disk orientations: one with the values derived by MATISSE observations as mentioned above (Sect.4), and the other following the geometric properties derived by the Atacama Large Millimeter/submillimeter Array (ALMA) observations , that is, a disk with PA= 43.6 • (±1.7) and an inclination i = 37.7 • (±0.8). We used both orientations to model both the inner (hot; >1000 K) region and the dusty (passive, cooler) region of the accretion disk.

A steady-state accretion disk model
The majority of the circumstellar emission at shorter wavelengths (λ ≤ 4µm) originates from the hot, inner accretion disk (i.e., within a radius of 1 au). In our first attempt to model this inner disk, we adopt a model of an optically thick but geometrically thin, viscous accretion disk, where the accretion rate remains constant irrespective of radial distance.
We follow an approach similar to Zhu et al. (2008), although here the synthetic SED is calculated by integrating black-body emission in concentric annuli between the disk's inner radius R in , and outer radius R out . Therefore, we focus on reproducing the broadband emission, unlike Zhu et al. (2008) that could simulate spectroscopic features. The disk temperature profile T (r) of this analytical model follows Hartmann & Kenyon (1996): Article number, page 5 of 17 A&A proofs: manuscript no. 42788corr_ver2 where r is the distance from the star, R * is the stellar radius, M is the stellar mass,Ṁ is the accretion rate, and G, σ are the wellknown gravitational and Stefan-Boltzmann constants.
We constrain the model SED based on the contemporary photometry up to 3 µm (Sect. 2.2, Table A.1) corrected for interstellar reddening (Sect. 2.3). The initial parameters for the steady-state disk model, that is, the radial extent of the hot inner accretion disk (R in , R out ) and the accretion rate (MṀ), were taken from Zhu et al. (2007) and Zhu et al. (2008). However, since the SED, the reddening, and the distance are now fixed to contemporary values, we found that the aforementioned model parameters had to be adjusted. As such, we opted for a range of R out at 0.2, 0.3, 0.6, 0.8, and 1.0 au, and adjusted the accretion rate (1 × 10 −5 ≤ MṀ[M 2 /yr] ≤ 3 × 10 −5 ) to obtain an adequate fit to the SED. The best results were achieved for outer radii at 0.3 and 0.6 au.
Although this simple, analytical approach provided a good description of the hot, inner disk of FU Ori, and could adequately fit the SED up to 2.5 µm, such a model cannot simulate the entire accretion disk, and especially its dusty component. We expect that the dusty passive disk would dominate at the mid-infrared and at longer wavelengths, while having a minor contribution in the near-infrared.
We therefore opted to model both the inner and the dusty (passive) accretion disk components with the Monte Carlo radiative transfer simulation tool radmc-3d (Dullemond et al. 2012), which can treat different geometries and dust compositions. The radiative transfer model is described in the next section.

Radiative transfer model
We modeled FU Ori's disk by employing the Monte Carlo radiative transfer simulation tool radmc-3d (Dullemond et al. 2012), which can treat different geometries and dust compositions. For each model setup, including the user-specified spatial distribution of circumstellar material, and the knowledge of the optical properties of the material, radmc-3d could solve the radiative equilibrium problem and find out the radiative equilibrium temperature distribution. Multiwavelength brightness distribution could then be computed, using the information of density and temperature distribution and opacity properties. We adopt a model setup similar to Zhu et al. (2008) and Pérez et al. (2020), including two components, that is, an inner accretion disk and an outer passive disk. The input parameters for the inner accretion disk were initially based on the best-fit results provided by the analytical model (see Sect. 5.1), but we allow for minor adjusting.
The spatial distribution of dust and the accretion heating source is described in the following. For each disk component, we assume a power-law density distribution. The disk density profile is described as where z is the height above the disk mid-plane, H is the scale height, and Σ is the surface density profile. The last two param- where q is the disk flaring parameter and h ref is the ratio of the scale height at the reference radius R ref , while where Σ ref is the surface density at the disk's inner radius and p is a power-law exponent. We set R ref equal to the radius where the inner and outer disk are separated.
To account for the accretion heating in the inner disk, we put the heating source in the mid-plane with power per area of

Dust composition
In the following we explain our choice of the dust composition of the passive disk. The shape of FU Ori's silicate feature indicates dust grain growth and deficiency of submicron grains. Quanz et al. (2006) fitted the silicate feature using a mixture composed mainly of amorphous silicate with olivine and pyroxene stoichiometry, with grain sizes from 0.1 to 6 µm. However, their study does not exclude the existence of larger grains (>10 µm). Based on this, we adopt a mixture of two dust components, which we named "small silicate" and "larger." The small silicate component, in turn, is the same mixture as described in the Table 4 of Quanz et al. (2006). The larger component is a mixture of carbonaceous material and "astronomical silicate" (Draine & Lee 1984) with a mass ratio of C:Si=1:2, and a powerlaw grain size distribution f (a) ∼ a −3.5 (Mathis et al. 1977), from a min = 10 µm to a max = 1000 µm. The mass ratio was set to this value because recent studies generally indicate that there is more silicate than carbon in protoplanetary disks (e.g., Habart et al. 2021). The small silicate and larger dust were then mixed with a mass ratio of f small : (1f small ), where f small is a free parameter and represents the percentage of small grains in the mixture. It is constrained by the strength of the silicate feature. The dust opacities were derived using the OpacityTool software (Toon & Ackerman 1981;Min et al. 2005;Woitke et al. 2016), which is based on the distribution of hollow spheres (DHS) theory (Min et al. 2005). In the DHS theory, for a ensemble of dust grains with irregular shape, the optical effects of grain shape are represented by a single shape parameter f max . OpacityTool can compute the optical properties for a given dust mixture, if the dust grain size and shape distribution, as well as the complex refractive index of the dust species as a function of wavelength are known. The complex refractive index data were collected from the literature (Dorschner et al. 1995 For the inner accretion disk, the opacity (τ) is dominated by gas free-free absorption, which is temperature-dependent, and therefore the radiative transfer process cannot be accurately simulated with radmc-3d. However, this does not significantly affect the optical and near-infrared continuum, as long as the gas density is high enough to ensure τ 1. As such, and for simplicity, we represent the inner disk with an artificial gray material with κ abs = 1 cm 2 g −1 and κ sca = 0 cm 2 g −1 at any wavelength.

Disk geometry
As described above (Sect. 5.1), the inner disk was modeled as a geometrically thin one. This presumes that its scale height is quite small compared with the dusty component. We first attempted to model the inner component as a "thick slab" geometry in radmc-3d, that is, with scale height h 0 and flaring index q 1 = 0. However, this could not reproduce the optical to near-infrared SED well, and we opted for a slightly flared inner disk.
Scattered-light images of FU Ori suggested a disk size of ∼80 au (Laws et al. 2020). However, at radio wavelengths the disk appears slightly more compact. Pérez et al. (2020) modeled the 1.3 mm continuum image using a disk with outer radius of 100 au, but with a characteristic radius of just 11.3 au, outside of which the surface density drops quickly. Considering this information, we also adopt in our model an outer radius of 100 au, but use a steep power-law for the surface density Σ ∝ r −2 , so that the disk is optically thick at 1.3 mm only in the inner ∼10 au region. We also tested models with smaller (60 au) and larger (150 au) outer radii, but these failed to reproduce the mid-infrared to submillimeter emission in the SED.
Under the framework described above, we adjusted the model parameters to fit the observations, including the contemporary SED, the correlated fluxes, and the closure phases in the MATISSE observations. As mentioned in the previous section, we found that analytical models with inner disk outer radii at 0.3 and 0.6 au, provided an adequate fit to the SED for the two disk orientations. However, when these were used as startup parameters for the radmc-3d model, models at 0.6 au produced more flux below 5µm, and models at 0.3 au provided better fits once the accretion rate was slightly reduced (order of 1 × 10 −5 M 2 yr −1 ). Table 1 lists the parameters for the two disk orientations. In general, we have not optimized our parametric search (e.g., with a χ 2 minimization or Markov chain Monte Carlo process) but selected parameters based on visual inspections of the data versus model comparisons. Therefore, we do not provide errors (or a range of values) for each parameter; we only present our initial parameter-space search (cf., Sect. 5.1). Certain parameters, such as the disk's orientation and the inner component's surface density, are fixed (boldface in Table 1), while some of the outer component's parameters are inherited from the inner one (e.g., R in,2 ≡ R out,1 ).
We illustrate in Fig. 6 the disk structure (i.e., disk density and temperature profiles) of the radmc-3d model with the ALMA orientation, noting that the model with the MATISSE orientation has a similar structure. We also provide a rudimentary sketch of FU Ori's disk in Fig. 7. Here, we indicate the regions of the disk that have been detected by MATISSE and by other highangular-resolution instruments. These are marked based either on direct images of the disk -for example from ALMA and the Gemini Planet Imager (GPI) -or on geometric sizes -from the Michigan InfraRed Combiner-eXeter (MIRC-X) of the Center for High Angular Resolution Astronomy (CHARA), GRAV-ITY/VLTI, and MATISSE/VLTI. The synthetic SEDs of the models described in Table 1 are shown in Fig. 8. The synthetic correlated fluxes at 3.5µm and at 10µm are compared with observations in Figs. 9 and 10, respectively, while the N-band synthetic closure phases are shown in Fig. 11. Furthermore, we show the synthetic visibilities and correlated fluxes in the entire N band in Figs. C.1 and C.2. Table 1. radmc-3d model parameters for two disk orientations. They correspond to the models shown in Figs. 6,[8][9][10][11]and C.

The current SED of FU Orionis
FU Orionis has been steadily fading since it reached peak brightness in 1937. Kenyon et al. (2000) estimated a decline rate of 0.015 mag yr −1 in the B band. This would suggest a drop in magnitude of 1.3 mag since 1937 (when the peak brightness was ∼9.8 mag and assuming m pg is equivalent to B) placing the current brightness at 11.1 mag in B.
Following this, we opted to complement our interferometric observations with contemporaneous photometric observations (Table A.1; Sect. 2.2), which enables us to model the current properties of FU Ori's accretion disk. Our photometric measurements are in agreement with the expected decline as above, taking into account that FU Ori itself is a short-term variable source (∆V ∼ 0.035 mag within 1 day; Kenyon et al. 2000).
Whilst comparing the SED from our observations to literature studies, we note a shift, a decline in brightness, in the nearinfrared as well. For example, the system has faded by approximately 0.5 mag in the K band within the last 20 years. Since the photometry indirectly influences the parametric search in the analytical and radiative transfer simulations (Fig. 8, upper panel;Sect. 5), it is expected that our model results will differ from previous works (e.g., Zhu et al. 2007;Pérez et al. 2020;Labdon et al. 2021). We would therefore caution the reader in comparing any new disk simulations to photometric measurements from the past literature. A further analysis of the past evolution will be presented in a forthcoming paper (Lykou et al., in preparation).

Radiative transfer model results
Our approach to simulating FU Ori's disk with different models and the subsequent "best-fit" parameters 3 was presented in Sect. 5. An initial parametric search was performed with a steady-state accretion disk model for the inner hot disk. This provided an adequate fit to the optical and near-infrared broadband photometry, but the models under-estimate the flux beyond 3 µm. We later used these parameters from the analytical approach, as a starting model for the radmc-3d simulations, where the (featureless) inner, hot accretion disk acts like an illuminating source for the dusty disk. These radiative transfer simulations required further adjustment of initial parameters, such as the accretion disk's outer radius and mass accretion rate, to achieve an adequate fit for the entire SED and interferometric data. In the following, we focus on the radmc-3d model results.
Despite the uncertainty in the absolute photometric calibration of the MATISSE N-band spectrum, we find that both models can adequately fit the MATISSE L-, M-, and N-band spectra within the uncertainties (Fig. 8). This includes the flat-topped silicate feature that is evident in all FU Ori spectra (Fig. 1).
A first major result of the radiative transfer analysis is that the inner disk's outer radius is found to be smaller (0.3 au) than literature values (e.g., Zhu et al. 2007 : 0.5 -1 au;Labdon et al. 2021 : 0.74 ± 0.35 au). This appears to agree with the MATISSE results where we find that the L-band emitting region ought to be smaller than 0.5 au in diameter. Liu et al. (2021) hypothesized that the hot (inner) disk region may be cooling down, which may be supported by a temporal decline in the SED, and could lead to a shrinkage of that region. When opting for larger (outer) radii, we found that the model flux was much higher in the 2 -6 µm and 10 -30 µm regions irrespective of the disk orientation. The "best-fit" models shown here (Fig. 8) overestimate the flux by In this work, we presume that the inner disk is optically thick and devoid of dust, although it could be surrounded by a halo of gaseous material. The dusty disk is flared, with the bulk of large grains settling in the midplane. We roughly indicate the extent of the regions (geometric sizes and/or direct images) that have been detected by MATISSE and by other high-angular-resolution instruments, such as MIRC-X and the Precision Integrated-Optics Near-infrared Imaging ExpeRiment (PIONIER) (H), GRAVITY (K), ALMA (1.3mm continuum), and the GPI and High-Contrast Coronographic Imager for Adaptive Optics (HiCIAO) (scattered light; J, H). The physical scales are given here in astronomical units, and they are converted to angular sizes, with the Gaia EDR3 distance to FU Ori adopted. Other than that, this drawing is not to scale. approximately 10% in the 2 -6 µm region. However, the fits fall within the uncertainties of the MATISSE photometry. Figure 9 shows a comparison of the synthetic correlated fluxes at 3.5 µm from the two disk orientations (ALMA as crosses, and MATISSE as squares) with the MATISSE data from the UTs and the ATs. At first glance, there is no clear distinction between the ALMA and MATISSE geometries. On the other hand, both disk models overestimate the correlated fluxes. We attribute this to the overestimation of the synthetic total flux in the 2 -6 µm region as discussed above.
By fixing the outer radius at 0.3 au for both disk orientations, we first had to adjust the accretion rate 4 for each model to 1.4 × 10 −5 and 2.7 × 10 −5 M yr −1 for the ALMA and MATISSE orientations, respectively (Table 1). The remaining parameters were kept similar for both models, except for f small , the percentage of small silicate grains in the dust mixture (Sect. 5.2.1). We find that by reducing f small from 10% (as used in the ALMA model) to 5%, the MATISSE model can also provide an adequate fit to the SED and visibilities.
A distinction between the two disk orientations can be seen in the N-band data. It appears that the ALMA orientation provides a better fit to both the correlated fluxes at 10 µm (Fig.10) and the closure phases (Fig. 11). This suggests that a pole-on geometry is more favorable. This reinforces our argument that the nonzero closure phases (Sect. 3.3) can be interpreted by a flared disk inclined to our line of sight.  Table A.1) and the MATISSE spectrum (red; epoch 4) against the analytical (green) and radmc-3d (black) models with the two disk orientations: ALMA (dashed lines) and MATISSE (solid lines). The bottom two panels are enlargements of the L-and N-band regions. The shaded regions represent the flux uncertainties of the MATISSE spectra, while the 9.4 -9.9µm spectral region is heavily affected by the terrestrial atmosphere.

Disk mass
Our models suggest a disk mass (dust) of 2.4 × 10 −4 M . This is in line with previous studies , which placed the dust mass at 2 × 10 −4 M for optically thick dust in the submillimeter regime. On the other hand, our result is a factor of 3 higher than the dust mass predicted by Pérez et al. (2020), and an order of magnitude less than what Liu et al. (2019) derived from radio observations (1.8 × 10 −3 M ).
For a typical gas-to-dust ratio of 100, similar to the interstellar medium value, we obtain a lower limit for the total disk mass (gas and dust) of ∼ 0.02 M . FU Ori's disk is therefore still one of the most compact and least massive in its class Kóspál et al. 2021).
The derived total disk mass is slightly higher than the expected amount of material accreted onto FU Ori since its eruption, but it is still a small amount. That is, if one assumes that the accretion rate remained constant at ∼ 1 × 10 −5 M yr −1 (at the order of the accretion rate derived by our models) for the first 85 years since the eruption, then the total accreted mass would have been 8.5×10 −4 M . The remaining disk mass, if it is indeed 0.02 M , would be accreted within 2000 years at this rate.

Dust mineralogy and radial variations
Here we explore any potential radial variations in the silicate feature within the disk itself. The radmc-3d models predict a small radial dependence of the silicate feature with respect to spatial sampling. That is, correlated spectra from the shortest baselines, which sample larger portions of the disk, indicate stronger silicate emission as opposed to the longest baselines, where the emission diminishes. This is of course expected, since the silicate emission is stronger when integrating over larger areas of the silicate-rich disk. On the other hand, the longest MATISSE baseline is probing a region of the simulated disk that is silicatefree in this case, that is, the inner accretion disk, and hence no silicate emission arises from that area. This hypothesis will be examined in the following section.
We compare the MATISSE/VLTI N-band correlated spectra with the Spitzer spectrum from 2004 (Green et al. 2006) in Fig. 12. For the sake of clarity, we normalize all spectra following the method of van Boekel et al. (2003). That is, we obtain a linear fit between 8.3 and 12.9 µm, which is defined as the "continuum," which we then subtract from the spectrum. As last step, we normalize the subtracted spectrum with the mean of the continuum. The final product is scaled to be larger than unity for clarity. Overall, we do not notice any significant differences in the shape of the silicate feature per baseline. In fact, the silicate feature appears to be flat-topped, and thus resembles emission arising from larger-sized grains, such as α ≥ 2 µm (e.g., Fig.2 of Bouwman et al. 2001).
Unlike other types of eruptive star disks (e.g., EXors), FUors do not seem to show any crystalline silicate features. FU Orionis itself was found to be devoid of crystallines by Quanz et al. (2006), and it would appear that the MATISSE observations corroborate the MIDI results. Malbet et al. (2005) argued that FU Ori might be a triple system, since modeling of their interferometric data (H and K bands) required an additional component from a point source (i.e., spot). The spot was estimated to be located at approximately 10 au (i.e., 25 mas) from the protostar and at a PA of about 130 • . Liu et al. (2019) reported a small deviation within the CO band in their GRAVITY/VLTI data, although the closure phase signal within errors is close to zero. Their geometric model fit to the data required an additional off-centered source, which they argue could be produced either by the presence of a nearby companion or by an inclined disk. Liu et al. (2019) and Pérez et al. (2020) have shown that the CO gas in the disk is in Keplerian rotation, and therefore one could argue that the GRAVITY/VLTI observations probed the inner part of that gaseous component of the disk. Furthermore, A&A proofs: manuscript no. 42788corr_ver2  Labdon et al. (2021) argue that their closure phase measurements are null within error bars and, as such, more consistent with a centro-symmetric distribution that is inclined to our line of sight.

On a potential third companion
Similarly, we cannot confirm the presence of a third companion 5 from our data. At a first glance, none of the interferometric observables appear to show the sinusoidal signal expected from a binary companion. However, since this could be hidden inside the disk's signal, we attempted to test this. Essentially, faint companions at small separations would have produced a sinusoidal signal with a very low visibility amplitude and a very broad mod- 5 The only confirmed companion, FU Ori South, is outside the fieldof-view of the UT array.  Green et al. 2006). ulation cycle. We opted for a companion similar to the separation (25 mas) and flux ratio (≤ 3%) of Malbet et al. (2005) at multiple PAs, but find that if such a companion existed, its visibility amplitude (≤ 0.05) would be undetectable with our current MA-TISSE N-band observations considering the large measurements in error (δV ∼ 0.1). We obtain similar results for companions with similar flux ratios but at larger separations (e.g., 100 mas).
Coincidentally, MATISSE can detect companions down to a flux ratio of 2% in the L band with the ATs, which offer a larger field-of-view, but the targets were wide binaries (orbital separation ∼ 100 mas; see Lopez et al. 2022). FU Ori is marginally resolved in the L band, while the visibilities are overall flat without showing any sinusoidal modulation (Fig. 2). In a similar fashion to the approach presented above, we can exclude any companions in the separation range of 20 − 100 mas with L-band flux ratios larger than 5%.

Investigating disk misalignments
Multiwavelength observations of protoplanetary disks have shown that under the influence of nearby companions (stellar or planetary), the inner regions of said disks can often become misaligned or warped. A striking example of disk misalignment is GW Ori (Bi et al. 2020;Kraus et al. 2020). Here, we investigate the possibility of disk misalignment for FU Orionis.
Near-infrared scattered-light images (Liu et al. 2016;Takami et al. 2018;Laws et al. 2020) show a spiral-arc feature eastward of FU Ori, while they suggest cavities in the north and northeastern directions. These large-scale (≥ 0.5 ) features could suggest an interaction with a companion, perhaps the known companion FU Ori S, or even signify infalling material from the extended circumstellar environment onto the accretion disk.
Observations by ALMA (Liu et al. 2019;Pérez et al. 2020) suggest a different orientation from any of the previous nearinfrared, mid-infrared and radio interferometric observations (Table 2). For example, Malbet et al. (2005) mention a best-fit Previously, Quanz et al. (2006) estimated from their MIDI/VLTI data an inclination similar to that of Malbet et al. (2005); however, the disk's PA was quite different. Puzzling results can also be found in previous works, such as the reconstructed 10.7 µm image of Monnier et al. (2009) and the model image of Malbet et al. (2005), which indicate that the disk's minor axis is oriented toward the northwest. Labdon et al. (2021) reanalyzed archival H-and K-band interferometric observations and presented new J-band data from CHARA/MIRC-X. Their temperature gradient model suggests a disk with a PA∼ 34 • and inclination of ∼ 32 • (cf. Table 2), which within errors agrees with Pérez et al. (2020).
We measured FU Ori disk's inclination and PA by fitting geometrical models to the N-band data. We attempted the same process for the GRA4MAT L-band data; however, since the disk is marginally resolved in this band, we focus more on the geometric properties derived from the N-band data. Within the uncertainties, the results agree very well with literature values (Table 2).
With the exception of the imaging results by Pérez et al. (2020) in the submillimeter regime and by Liu et al. (2017) in the radio regimes, all estimates for the disk's orientation were based on geometric model fits to interferometric data. In this work, we present radiative transfer models that used both the orientation derived by geometric fits to the MATISSE data and the orientation by ALMA . The latter appears to be more favorable, suggesting that the dusty disk is more pole-on in our line of sight. We cannot strongly conclude any disk misalignment for this system at present, especially for the hot accretion disk. We expect that even higher-angular-resolution instruments (θ ≤ 1 mas) can provide more insight into that.

Conclusions
We have presented new insights into the accretion disk of FU Ori in the mid-infrared with the use of the interferometric instrument MATISSE/VLTI. Our results are complemented by radia-tive transfer simulations that attempt to constrain the properties of the disk.
In summary, we find that: -The accretion disk is very compact in the L band as it is marginally resolved at <2 mas in size, or 0.8 au at the adopted Gaia distance. Therefore, the hot accretion disk's radius ought to be smaller than 0.4 au. Similarly, MATISSE was able to detect just the innermost, and possibly warmer, part of the passive dusty disk, with an approximate size of 5 mas in the N band, which translates to about 2 au. -Geometric-model fits (2D Gaussian) to the MATISSE Nband correlated fluxes suggest a disk orientation with a minor-axis PA of 15 ± 25 degrees and an inclination of 55 ± 15 degrees, which are similar to literature values from MIDI/VLTI (Quanz et al. 2006) but differ from more recent imaging studies by ALMA ). -Two radiative transfer models were explored for the two different disk orientations (MATISSE and ALMA). Both can provide relatively good fits to the SED and L-band data; however, a distinction can be seen when the models are compared with the N-band interferometric data (correlated fluxes and closure phases). It appears that an orientation that is more pole-on (i.e., ALMA) is more favorable. Since this discrepancy could allude to potential disk misalignment, we opt to reexplore this in future observations at even higher angular resolutions (e.g., GRAVITY/VLTI). -Our model fits suggest an average accretion rate of about 2 × 10 −5 M yr −1 (if we assume a stellar mass of 0.6 M ), which is somewhat lower than literature values. -There are no signatures of crystalline silicate features in either the total flux or the correlated N-band spectra, corroborating earlier results from MIDI/VLTI. The silicate feature itself is "flat-topped," suggesting it emanates from large-sized grains (size≥ 2 µm). -Our model predicts a dust mass of 2.4×10 −4 M for the outer passive disk. A lower limit for the total disk mass (gas and dust) can be obtained by assuming a typical gas-to-dust ratio of 100, that is, M total disk ∼ 0.02 M . However, we do not obtain any meaningful constraint on the mass of the inner gaseous disk other than it ought to be massive enough to be optically thick.
-The current MATISSE N-band observations cannot constrain the presence of an, as yet un-detected, tertiary companion for FU Ori. In the L-band data, the absence of any sinusoidal modulation excludes any companions with a separation of 20 − 100 mas and a flux ratio ≥ 5% in L.
FU Orionis, as the archetype of its class, is still an active laboratory for the study of post-eruptive accretion events, since it continues to fade in the optical and near-infrared and is not expected to return to its pre-eruption phase within the next 25 years.
We conclude that the MATISSE observations indicate that the hot, inner accretion disk is rather compact, with a diameter ≤1 au at 3.5µm. Since earlier studies suggested a region of 1-2 au, this could indicate that the hot emitting region has been shrinking as the system fades.
Due to its low brightness and its compactness, FU Orionis is not an ideal candidate for imaging in the mid-infrared, either with MATISSE or with aperture masking at larger spatial scales (e.g., VISIR/VLT). Imaging has also proven to be difficult with other near-infrared interferometers such as CHARA (Labdon et al. 2021).
Epoch 1: The detection limits for the AT configurations in lowspectral-resolution mode are 1, 5, and 10 Jy, respectively, for the L, M, and N bands. We found that FU Ori was fainter than expected, that is, ≤ 5 Jy in all three bands, so the first attempt to observe it with an AT configuration was not successful. Although the star was slightly brighter than the L-band limit, the data were of low signal-to-noise. The N-band data are of very poor quality. We therefore reject these observations. Epoch 2: Data were obtained in adverse atmospheric conditions with seeing as high as 1.5 and atmospheric coherence time ≤ 4 ms. This affected mostly the L and M observations as opposed to the N band; however, we also found problems with the absolute photometric performance of the UTs in N (total flux spectra). This data set was used only as qualitative comparison for epoch 4 N-band data.
Epoch 3: No matching calibrator was observed in this epoch. Therefore, that data set is not included in this work.
Epochs 4 and 5: Technical issues in epoch 4 hindered the observations. The designated calibrator (HD37160) was found to be brighter than expected and saturated; therefore, an additional calibrator (HD47886) was observed immediately after. However, another technical error occurred, and the connection to UT2 telescope was lost; therefore, only three of the six baselines could be calibrated for that run (HD47886). We therefore opted to recalibrate epoch 4 data with the first calibrator of that night's GTO program (HD28413), which were obtained at a similar airmass as FU Ori. In general, we find that the chopped interferometric data are of better quality in L and M, since temporal variations in atmospheric transmission in the mid-infrared could be corrected with chopping. We also find that the total spectra, the correlated fluxes, and the visibilities from this third calibrator agree within 10% with the results from HD47886, and moreover the total spectra are nearly identical to those obtained with a different mode in epoch 5. Therefore, here we show the data sets calibrated with HD28413. We also note that all N-band data obtained here in high-spectral-resolution mode suffer from poor S/N beyond 11µm. This may be the result of incorrect flatfielding, manifesting as correlated noise features, because this spectral mode utilizes a larger portion of the detector that has not been characterized well yet 7 .
Epoch 5*: In addition to the low-spectral-resolution data in LM for epoch 5, we obtained data in medium-spectral resolution in LM and tested the performance of the high-spectral-resolution mode in N. It proved unsuccessful since the brightness limit for that mode is found to be 26 Jy for the correlated fluxes. The Lband medium-spectral resolution data are beyond the scope of this work and will be analyzed in a future publication.