11institutetext: Konkoly Observatory, Research Centre for Astronomy and Earth
Sciences, Eötvös Loránd Research Network (ELKH), Konkoly-Thege Miklós út
15-17, 1121 Budapest, Hungary 11email<EMAIL_ADDRESS>22institutetext:
ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter sétány 1/A,
1117 Budapest, Hungary 33institutetext: Leiden Observatory, Leiden University,
P.O. Box 9513, 2300 RA, Leiden, The Netherlands 44institutetext: Max-Planck-
Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany
55institutetext: Laboratoire Lagrange, Université Côte d’Azur, Observatoire de
la Côte d’Azur, CNRS, Boulevard de l’Observatoire, CS 34229, 06304, Nice Cedex
4, France 66institutetext: Max-Planck-Institut für Radioastronomie, Auf dem
Hügel 69, 53121 Bonn, Germany 77institutetext: Scottish Universities Physics
Alliance (SUPA), School of Physics and Astronomy, University of St Andrews,
North Haugh, St Andrews, KY16 9SS, UK 88institutetext: Department of Physics
and Astronomy, University of Nevada, Las Vegas, 4505 S. Maryland Parkway, Las
Vegas, NV 89154, USA 99institutetext: 11F of AS/NTU Astronomy-Mathematics
Building, No.1, Sect. 4, Roosevelt Rd, Taipei 10617, Taiwan, R.O.C.
1010institutetext: Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
1111institutetext: Anton Pannekoek Institute for Astronomy, University of
Amsterdam, Science-Park 904, 1098 XH, Amsterdam, The Netherlands
1212institutetext: AIM, CEA, CNRS, Université Paris-Saclay, Université Paris
Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France 1313institutetext:
Institute for Mathematics, Astrophysics and Particle Physics, Radboud
University, PO Box 9010, MC 62, 6500 GL, Nijmegen, The Netherlands
1414institutetext: SRON Netherlands Institute for Space Research, Niels
Bohrweg 4, 2333 CA, Leiden, The Netherlands 1515institutetext: Institut für
Theoretische Physik und Astrophysik, Christian-Albrechts-Universität zu Kiel,
Leibnizstraße 15, 24118, Kiel, Germany

# The disk of FU Orionis viewed with MATISSE/VLTI††thanks: Based on
observations collected at the European Southern Observatory under ESO programs
0104.C-0782(B), 0104.C-0016(D), 0106.C-0501(D), and 0106.C-0501(F).

First interferometric observations in $L$ and $M$ bands.

F. Lykou 11 P. Ábrahám 1122 L. Chen 11 J. Varga 3311 Á. Kóspál 112244 A.
Matter 55 M. Siwak 11 Zs. M. Szabó 667711 Z. Zhu 88 H. B. Liu 99 B. Lopez 55
F. Allouche 55 J.-C. Augereau 1010 P. Berio 55 P. Cruzalèbes 55 C. Dominik
1111 Th. Henning 44 K.-H. Hofmann 66 M. Hogerheijde 331111 W. J. Jaffe 33 E.
Kokoulina 55 S. Lagarde 55 A. Meilland 55 F. Millour 55 E. Pantin 1212 R.
Petrov 55 S. Robbe-Dubois 55 D. Schertl 66 M. Scheuck 44 R. van Boekel 44 L.
B. F. M. Waters 13131414 G. Weigelt 66 S. Wolf 1515

(Received Dec 01, 2021; accepted Apr 20, 2022)

###### Abstract

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 configurations. We also obtained contemporaneous photometry in the
optical ($UBVRIr^{\prime}i^{\prime}$; SAAO and Konkoly Observatory) and near-
infrared ($JHK_{s}$; NOT). Our results were compared with radiative transfer
simulations modeled by radmc-3d.

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 $\sim 1.3\pm 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 $\sim 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-3d provide adequate fits to the
spectral energy distribution from the optical to the submillimeter and to the
interferometric observables when opting for an accretion rate $\dot{M}\sim
2\times 10^{-5}$ $M_{\sun}$ yr-1 and assuming $M_{*}=0.6$ $M_{\sun}$. Most
importantly, the hot inner accretion disk’s outer radius can be fixed 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\times 10^{-4}$ $M_{\sun}$, and for a typical
gas-to-dust ratio of 100, the total mass in the disk is approximately 0.02
$M_{\sun}$. We did not find 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 ($\leq 1$ mas).

###### Key Words.:

Stars: individual: FU Ori – protoplanetary disks – circumstellar matter –
infrared: stars – techniques: interferometric – radiative transfer

##  1 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
$\sim 10^{-4}$ $M_{\sun}$ 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
$\Delta B\sim 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_{\sun}$;
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 arc-shaped 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 $\sim$0.5″ away at a position angle (PA) of $\sim
161\degr$ 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 pre-main-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 (Hales et al., 2015; Liu et al., 2017, 2019;
Pé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, 2005; Quanz 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 (Pérez
et al., 2020). 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 $\sim 10^{-4}$ $M_{\sun}$ yr-1. For instance, the
boundary between the hot inner disk and the passive, cooler component is found
at $0.76\pm 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 $\sim 0.5$ au (at a distance of 500
pc), with a not well-confined truncation radius of $\sim 1$ au between the
inner and the dusty disk, as it is not clear which disk component dominates in
emission in the 4 - 8 $\rm\mu$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 $\rm\mu m$), and $N$ (8 - 13 $\rm\mu 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.

##  2 Observations

###  2.1 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 4) 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 $\sim$ 3 - 9 mas in $L$ and $\sim$ 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 4) 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\sim 30$) in $L$ and $M$ and for high-
spectral resolution ($R\sim 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 (HighSens) 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
MATISSE (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. Table 4). The atmospheric transfer function
was estimated with Mid-infrared stellar Diameters and Fluxes compilation
Catalogue (MDFC; Cruzalèbes et al., 2019) calibrators HD48433 (K0III), HD28413
(K4III), HD47886 (M1III), HD37160 (G9.5III), and HD31767 (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
tools111https://gitlab.oca.eu/MATISSE/tools/-/wikis/home. Since FU Ori is
quite faint in all bands, that is, below the minimum suggested from MATISSE
specifications222https://www.eso.org/sci/facilities/paranal/instruments/matisse/inst.html,
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.

###  2.2 Optical and infrared photometry

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

####  2.2.1 Nordic Optical Telescope

Photometric observations were obtained with the NOT near-infrared 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 $Ks$ filters. The detector array ($1024\times 1024$
pixels) offers a field-of-view of 80″$\times$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 introduced that reduced the transmission
down to 10%. The average seeing on both nights was $\leq 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 ($\sim 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 $Ks$, respectively, with a typical uncertainty
$\delta m\sim 0.03$ mag. The NOTCam photometry is included in Table 3.

####  2.2.2 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 $UBVRI$ filters, mounted on
the 1-meter Lesedi telescope. The instrument offers a 5.72′$\times$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 3, where $UBV$ are in the Bessel system, and
$R_{c}I_{c}$ in the Cousins system.

####  2.2.3 Konkoly Observatory

Photometry was obtained with the fully automated Astro Systeme Austria (ASA)
AZ800 alt-azimuth, direct-drive, 80-centimeter Ritchey–Chrt́ien (RC) telescope
at the Piszkéstető Observatory, Hungary, on 13–18 February 2021. We used a
2048$\times$2048 pixel FLI PL230 CCD camera with the E2V CCD230-42 detector.
The optical setup with an effective focal distance of $f=5600\,{\rm mm}$
yields a pixel scale of 0.55″ and a field-of-view of 18.8″$\times$18.8″. We
obtained three images per night with Bessel $B$, $V$ and Sloan $r^{\prime}$
and $i^{\prime}$ 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 3, where $BV$ are in the Bessel system and
$r^{\prime}i^{\prime}$ in the Sloan system. The typical uncertainty of the
photometry is 0.01 mag.

Figure 1: MATISSE total flux spectra. Top: $L$ and $M$ band from the UTs
(blue) and the ATs (red; GRA4MAT) from epochs 4 and 5, respectively (Table 4).
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.

###  2.3 Interstellar extinction

The neutral hydrogen column density, $N_{H}$, in a region within 0.1° of FU
Ori’s location is $1.99\times 10^{21}$ atoms/cm2 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\pm 0.09)\times 10^{21}\,A_{V}$, we find
a low interstellar extinction of $A_{V}=0.90\pm 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\pm
0.1$ mag.

##  3 Interferometric data

###  3.1 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 $\rm\mu$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 MATISSE 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$\rm\mu$m. Despite the higher sensitivity
offered with the UTs on MATISSE, the S/N decreases beyond 11 $\rm\mu$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.

###  3.2 Visibilities and correlated fluxes

####  3.2.1 L and M bands

The squared visibilities, $V^{2}$, and correlated fluxes, $F_{\rm 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_{\rm corr}$ below 3.2 $\rm\mu$m
and beyond 3.9 $\rm\mu$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\leq 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\leq V^{2}\leq 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 $\leq 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$\alpha$, CO), but this
should be expected in this case where low-spectral-resolution setups were
used.

####  3.2.2 N band

The $N$-band data from epoch 4 are shown in Fig. 4. FU Ori is known to have
the 10 $\rm\mu$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\leq B\leq 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 $\rm\mu$m), otherwise within a radius of 2 au from the
protostar. However, that area accounts for less than 50% of the 10 $\rm\mu$m
emission since the average visibility is $\leq 0.5$. The majority of the 10
$\rm\mu$m silicate dust emission arises from a larger region of the disk,
which can be probed by MATISSE with the shortest UT baselines ($B\sim 40$ m).
Its size ought to be smaller than roughly 25 mas at 10 $\rm\mu$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 near-infrared 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\leq 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).

Figure 2: MATISSE squared visibilities, correlated fluxes, and closure phases
in the $L$ band. The first row shows data from epoch 4 (UTs) and the bottom
row data from epoch 5 (ATs). The visibilities and correlated fluxes are color-
coded with respect to the baselines, as shown in the legends, while closure
phases have a different color scheme (per baseline triangle). The apparent
flaring in the data below 3.2$\rm\mu$m and beyond 3.9$\rm\mu$m is a known
instrumental artifact at the edges of the atmospheric transmission windows.

Figure 3: As in Fig. 2 but for the $M$ band.
Figure 4: As in Fig. 2 but for $N$-band data of epoch 4.

###  3.3 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 $\rm\mu$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).

##  4 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\pm 0.1$ mas at
3.5 $\rm\mu$m, which suggests an apparent size of about 0.5 au for the adopted
Gaia distance (Fig. 5).

Figure 5: Gaussian (1D) model fit compared with the epoch 4 and 5 correlated
fluxes at 3.5$\rm\mu$m.

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\pm 1.4$ mas at 10.5 $\rm\mu$m, a
PA of $15\pm 25$° for the minor axis, and an inclination, $i$, of $55\pm 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 $\leq 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.

##  5 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 (Pérez et al., 2020), that
is, a disk with PA$=43.6^{\circ}(\pm 1.7)$ and an inclination
$i=37.7^{\circ}(\pm 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.

###  5.1 A steady-state accretion disk model

The majority of the circumstellar emission at shorter wavelengths
($\lambda\leq 4\rm\mu 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_{\rm in}$, and outer radius $R_{\rm
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):

$T(r)=\left[\frac{3GM\dot{M}}{8\pi
R^{3}\sigma}\left(1-\sqrt{\frac{R_{*}}{r}}\right)\right]^{1/4},$ (1)

where $r$ is the distance from the star, $R_{*}$ is the stellar radius, $M$ is
the stellar mass, $\dot{M}$ is the accretion rate, and $G,\sigma$ are the
well-known gravitational and Stefan-Boltzmann constants.

We constrain the model SED based on the contemporary photometry up to 3
$\rm\mu$m (Sect. 2.2, Table 3) 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_{\rm in},\,R_{\rm out}$) and
the accretion rate ($M\dot{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_{\rm out}$
at 0.2, 0.3, 0.6, 0.8, and 1.0 au, and adjusted the accretion rate ($1\times
10^{-5}\leq M\dot{M}[M_{\odot}^{2}/yr]\leq 3\times 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
$\rm\mu$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.

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

$\rho=\frac{\varSigma(r,\phi)}{H\sqrt{2\pi}}\exp{\left(-\frac{z^{2}}{2H^{2}}\right)}\hskip
5.0pt,$ (2)

where $z$ is the height above the disk mid-plane, $H$ is the scale height, and
$\varSigma$ is the surface density profile. The last two parameters can be
described by power-law relations as

$\frac{H(r)}{r}=h_{\rm ref}\left(\frac{r}{R_{\rm ref}}\right)^{q},$ (3)

where $q$ is the disk flaring parameter and $h_{\rm ref}$ is the ratio of the
scale height at the reference radius $R_{\rm ref}$, while

$\varSigma(r)=\varSigma_{\rm ref}\left(\frac{r}{R_{\rm ref}}\right)^{-p},$ (4)

where $\Sigma_{\rm ref}$ is the surface density at the disk’s inner radius and
$p$ is a power-law exponent. We set $R_{\rm 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

$P=\frac{3GM\dot{M}}{4{\pi}r^{3}}\left[1-\left(\frac{R_{*}}{r}\right)^{1/2}\right],\,\mathrm{where}\,\,R_{\mathrm{in,1}}<r<R_{\mathrm{out,1}}.$
(5)

####  5.2.1 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 $\mu$m. However, their
study does not exclude the existence of larger grains
(${>}10\leavevmode\nobreak\ \mu$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 power-law grain size distribution $f(a)\sim a^{-3.5}$ (Mathis
et al., 1977), from $a_{\rm min}=10$ $\rm\mu$m to $a_{\rm max}=1000$
$\rm\mu$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_{\mathrm{small}}$:(1-$f_{\mathrm{small}}$),
where $f_{\mathrm{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_{\mathrm{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, amorphous silicate; Draine & Lee 1984, Laor & Draine
1993, Weingartner & Draine 2001, astronomical silicate; Jager et al. 1998,
carbonaceous material). We adopt a grain shape parameter of $f_{\rm max}=0.8$
following Woitke et al. (2016).

For the inner accretion disk, the opacity ($\tau$) 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 $\tau\gg 1$. As such, and for
simplicity, we represent the inner disk with an artificial gray material with
$\kappa_{\rm abs}=1$ cm2g-1 and $\kappa_{\rm sca}=0$ cm2g-1 at any wavelength.

####  5.2.2 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\neq 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
${\sim}80\leavevmode\nobreak\ $au (Laws et al., 2020). However, at radio
wavelengths the disk appears slightly more compact. Pérez et al. (2020)
modeled the $1.3\leavevmode\nobreak\ $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 $\Sigma\propto r^{-2}$, so that the
disk is optically thick at $1.3$ mm only in the inner ${\sim}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 start-up parameters for the radmc-3d model,
models at 0.6 au produced more flux below 5$\rm\mu$m, and models at 0.3 au
provided better fits once the accretion rate was slightly reduced (order of
$1\times 10^{-5}\,M_{\sun}^{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 $\chi^{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_{\rm in,2}\equiv R_{\rm 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 high-angular-
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), GRAVITY/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$\rm\mu$m and at 10$\rm\mu$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.
13 and 14.

Table 1: radmc-3d model parameters for two disk orientations. They correspond
to the models shown in Figs. 6, 8-11, and 13-14. Fixed parameters are marked
in boldface font.

Parameter | ALMA | MATISSE
---|---|---
PA (deg) | 43.6 | 15
$\bm{i}$ (deg) | 37.7 | 55
Inner disk (Component 1)
$M\dot{M}$ ($M_{\sun}^{2}/\rm{yr}$) | $0.85\times 10^{-5}$ | $1.64\times 10^{-5}$
$R_{\rm in,1}\,(R_{\sun})$ | 1.98 | 2.33
$R_{\rm out,1}$ (au) | 0.30 | 0.30
$\bm{\Sigma}_{\rm in,1}$ ($\,\rm{g/cm^{2}}$) | 10.0 | 10.0
$\bm{p}_{1}$ | 0 | 0
$h_{\rm ref,1}$ | 0.06 | 0.06
$q_{1}$ | 0.13 | 0.13
Outer disk (Component 2)
$R_{\rm in,2}$ (au) | $\equiv R_{\rm out,1}$ | $\equiv R_{\rm out,1}$
$R_{\rm out,2}$ (au) | 100 | 100
$\varSigma_{\rm in}$ ($\rm{g/cm^{2}}$) | 940 | 904
$p_{2}$ | $-2.0$ | $-2.0$
$h_{\rm ref,2}$ | $\equiv h_{\rm ref,1}$ | $\equiv h_{\rm ref,1}$
$q_{2}$ | 0.18 | 0.18
$f_{\mathrm{small}}$ | 0.10 | 0.05
333The dust mass derived from these models is $2.4\times 10^{-4}$ $M_{\sun}$
and $2.3\times 10^{-4}$ $M_{\sun}$ for the ALMA and MATISSE orientations,
respectively.

Figure 6: Disk structure in the best-fit radmc-3d radiative transfer model
with ALMA orientation. Top: Disk density profile. The colorbar indicates units
of g cm-3. Bottom: Disk temperature profile. The colorbar is in units of
kelvin. The vertical dashed lines mark the edge between the inner accretion
disk and outer passive disk (0.3 au). In the temperature distribution, contour
lines for temperatures of $3000$ K, $1500\leavevmode\nobreak\ $K, and $300$ K
are plotted.  Figure 7: Rudimentary sketch of the FU Orionis accretion disk,
indicating the inner, hot region (dark gray area; component 1) and the dusty,
passive disk (light gray area; component 2), which at 0.3 au are separated. 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 mid-plane. 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. Figure 8: Spectral energy distribution comparing the contemporary
photometry (blue circles; Table 3) 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$\rm\mu$m spectral region is heavily affected by the terrestrial
atmosphere.

##  6 Discussion

###  6.1 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 $\sim$9.8 mag and assuming $m_{\rm 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 3; 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 ($\Delta V\sim
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 near-infrared 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).

Figure 9: Synthetic correlated fluxes at 3.5 $\rm\mu$m from the best-fit
radmc-3d models of Table 1 (ALMA: crosses; MATISSE: squares) compared with
MATISSE data, color-coded per respective baseline: UTs (top panel) and ATs
(bottom panel).  Figure 10: As in Fig. 9 but for 10 $\rm\mu$m and only for the
UT (epoch 4) data. Figure 11: Synthetic $N$-band closure phases from the two
radmc-3d disk models (ALMA model as a solid black line, and MATISSE model as a
dotted line) compared with the MATISSE data (red; epoch 4).

###  6.2 Radiative transfer model results

Our approach to simulating FU Ori’s disk with different models and the
subsequent “best-fit” parameters444These models are not unique, given the
large number of free parameters, but we note them here as “best-fitting.”
Future observations may provide better constraints for these models. 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 $\rm\mu$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\pm 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 $\rm\mu$m and 10 - 30 $\rm\mu$m regions
irrespective of the disk orientation. The “best-fit” models shown here (Fig.
8) overestimate the flux by approximately 10% in the 2 - 6 $\rm\mu$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
$\rm\mu$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 $\rm\mu$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 rate555Adopting a stellar mass of 0.6 $M_{\sun}$
(Pérez et al., 2020).  for each model to $1.4\times 10^{-5}$ and $2.7\times
10^{-5}$ $M_{\sun}$ yr-1 for the ALMA and MATISSE orientations, respectively
(Table 1). The remaining parameters were kept similar for both models, except
for $f_{\rm small}$, the percentage of small silicate grains in the dust
mixture (Sect. 5.2.1). We find that by reducing $f_{\rm 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 $\rm\mu$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.

####  6.2.1 Disk mass

Our models suggest a disk mass (dust) of $2.4\times 10^{-4}$ $M_{\sun}$. This
is in line with previous studies (Hales et al., 2015), which placed the dust
mass at $2\times 10^{-4}$ $M_{\sun}$ 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\times
10^{-3}$ $M_{\sun}$).

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 $\sim
0.02$ $M_{\sun}$. FU Ori’s disk is therefore still one of the most compact and
least massive in its class (Liu et al., 2018; 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
$\sim 1\times 10^{-5}$ $M_{\sun}$ 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\times 10^{-4}$ $M_{\sun}$. The
remaining disk mass, if it is indeed 0.02 $M_{\sun}$, would be accreted within
2000 years at this rate.

###  6.3 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 silicate-free 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 $\rm\mu$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 $\alpha\geq
2\,\rm\mu$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.

Figure 12: Normalized correlated spectra from MATISSE (color-coded per
baseline) against Spitzer (in black; Green et al., 2006).

###  6.4 On a potential third companion

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

Similarly, we cannot confirm the presence of a third companion666The only
confirmed companion, FU Ori South, is outside the field-of-view of the UT
array. 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 modulation cycle. We opted for a companion similar to the
separation (25 mas) and flux ratio ($\leq 3\%$) of Malbet et al. (2005) at
multiple PAs, but find that if such a companion existed, its visibility
amplitude ($\leq 0.05$) would be undetectable with our current MATISSE
$N$-band observations considering the large measurements in error ($\delta
V\sim 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 $\sim 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%.

###  6.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
($\geq 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 near-infrared, mid-infrared and
radio interferometric observations (Table 2). For example, Malbet et al.
(2005) mention a best-fit model to their $H$\- and $K$-band data, with an
inclination of $\sim 55$°, as opposed to an earlier work (Malbet et al., 1998)
where a value of $\sim 30$° was measured in the $K$ band. Radio continuum
observations at 33 GHz (Karl G. Jansky Very Large Array (JVLA); Liu et al.,
2017) provide a highly uncertain disk PA, while the axis ratio suggests an
inclination of $34^{+18}_{-8}$ deg, which is similar to the ALMA results.

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 $\rm\mu$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$\sim 34$° and inclination of $\sim 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 (Pérez et al., 2020). 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 ($\theta\leq 1$ mas) can provide more insight
into that.

Table 2: Accretion disk PA (minor axis) and inclination derived from interferometric data. Reference | Band | PA (°) |  $i$ (°)
---|---|---|---
Malbet et al. (1998) | $K$ | – | $\sim 30$
Malbet et al. (2005) | $H,K$ | 47${}^{+7}_{-11}$ | 55${}^{+5}_{-7}$
Quanz et al. (2006) | $N$ | $19.1/3.4$ | $55.4\pm 2.4$
Monnier et al. (2009) | $N$ | $-25\pm 35$ | –
Liu et al. (2017) | 33 GHz | $7.9\pm 66$ | $34^{+18}_{-8}$
Pérez et al. (2020) | 1.3 mm | $43.6\pm 1.7$ | $37.7\pm 0.8$
Labdon et al. (2021) | $J,H,K$ | $34\pm 11$ | $32\pm 4$
this work (UT; 2021) | $N$ | $15\pm 25$ | $55\pm 15$
this work (AT; 2021) | $L$ | $67\pm 10$ | $48^{+7}_{-8}$

##  7 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 radiative 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 $N$-band correlated fluxes
suggest a disk orientation with a minor-axis PA of $15\pm 25$ degrees and an
inclination of $55\pm 15$ degrees, which are similar to literature values from
MIDI/VLTI (Quanz et al., 2006) but differ from more recent imaging studies by
ALMA (Pérez et al., 2020).

  * •

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\times 10^{-5}$
$M_{\sun}$ yr-1 (if we assume a stellar mass of 0.6 $M_{\sun}$), 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$\geq 2\,\rm\mu$m).

  * •

Our model predicts a dust mass of $2.4\times 10^{-4}\,M_{\sun}$ 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^{\rm
total}_{\rm disk}\sim 0.02$ $M_{\sun}$. 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 $\geq 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 $\leq$1 au at 3.5$\rm\mu 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).

Modeling parameters could be constrained further by using observations in
different array configurations than the ones presented here, should further
advances in instrumentation allow observations of such faint targets. High-
angular-resolution and high-spectral-resolution observations in the $K$ band
(e.g., GRAVITY/VLTI) may be able to disentangle the gaseous and continuum
emission at the center of the disk, and thus allow a combination of gas and
dust modeling from the near- to the mid-infrared.

###### Acknowledgements.

F.L. would like to thank the NOTCam team, and especially A. A. Djupvik, for
their assistance in the planning and the qualitative analysis of the
observations. F.L. and P.A. are funded from the Hungarian NKFIH OTKA project
no. K-132406. A.K., L.C., and M.S. are supported through ERC grant agreement
No 716155 (SACCRED). J.V. is supported by NOVA, the Netherlands Research
School for Astronomy. Our study is supported by the Lendület grant LP2012-31
of the Hungarian Academy of Sciences (MTA), and by the grant GINOP
2.3.2-15-2016-00033 of the National Research, Development and Innovation
Office (NKFIH), Hungary, funded by the European Union. Based on observations
collected at the European Southern Observatory under ESO program
0104.C-0782(B), 0104.C-0016(D), 0106.C-0501(D) and 0106.C-0501(F), and on
observations made with the Nordic Optical Telescope, owned in collaboration by
the University of Turku and Aarhus University, and operated jointly by Aarhus
University, the University of Turku and the University of Oslo, representing
Denmark, Finland and Norway, the University of Iceland and Stockholm
University at the Observatorio del Roque de los Muchachos, La Palma, Spain, of
the Instituto de Astrofisica de Canarias. MATISSE was designed, funded and
built in close collaboration with ESO, by a consortium composed of institutes
in France (J.-L. Lagrange Laboratory – INSU-CNRS – Côte d’Azur Observatory –
University of Côte d’Azur), Germany (MPIA, MPIfR and University of Kiel), the
Netherlands (NOVA and University of Leiden), and Austria (University of
Vienna). The Konkoly Observatory and Cologne University have also provided
some support in the manufacture of the instrument. This research has made use
of the services of the ESO Science Archive Facility.

## References

  * Abbott et al. (2000) Abbott, T. M., Aspin, C., Sorensen, A. N., et al. 2000, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4008, Optical and IR Telescope Instrumentation and Detectors, ed. M. Iye & A. F. Moorwood, 714–719
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147
  * Bi et al. (2020) Bi, J., van der Marel, N., Dong, R., et al. 2020, ApJ, 895, L18
  * Bouwman et al. (2001) Bouwman, J., Meeus, G., de Koter, A., et al. 2001, A&A, 375, 950
  * Calvet et al. (1991) Calvet, N., Hartmann, L., & Kenyon, S. J. 1991, ApJ, 383, 752
  * Connelley & Reipurth (2018) Connelley, M. S. & Reipurth, B. 2018, ApJ, 861, 145
  * Coppejans et al. (2013) Coppejans, R., Gulbis, A. A. S., Kotze, M. M., et al. 2013, PASP, 125, 976
  * Cruzalèbes et al. (2019) Cruzalèbes, P., Petrov, R. G., Robbe-Dubois, S., et al. 2019, MNRAS, 490, 3158
  * Dorschner et al. (1995) Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503
  * Draine & Lee (1984) Draine, B. T. & Lee, H. M. 1984, ApJ, 285, 89
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93
  * Green et al. (2013) Green, J. D., Evans, Neal J., I., Kóspál, Á., et al. 2013, ApJ, 772, 117
  * Green et al. (2006) Green, J. D., Hartmann, L., Calvet, N., et al. 2006, ApJ, 648, 1099
  * Green et al. (2016) Green, J. D., Jones, O. C., Keller, L. D., et al. 2016, ApJ, 832, 4
  * Güver & Özel (2009) Güver, T. & Özel, F. 2009, MNRAS, 400, 2050
  * Habart et al. (2021) Habart, E., Boutéraon, T., Brauer, R., et al. 2021, A&A, 649, A84
  * Hales et al. (2015) Hales, A. S., Corder, S. A., Dent, W. R. D., et al. 2015, ApJ, 812, 134
  * Hartmann (2009) Hartmann, L. 2009, Accretion Processes in Star Formation: Second Edition (Cambridge: Cambridge University Press)
  * Hartmann & Kenyon (1985) Hartmann, L. & Kenyon, S. J. 1985, ApJ, 299, 462
  * Hartmann & Kenyon (1996) Hartmann, L. & Kenyon, S. J. 1996, ARA&A, 34, 207
  * Henden et al. (2016) Henden, A. A., Templeton, M., Terrell, D., et al. 2016, VizieR Online Data Catalog, II/336
  * Henning et al. (1998) Henning, T., Burkert, A., Launhardt, R., Leinert, C., & Stecklum, B. 1998, A&A, 336, 565
  * Herbig (1966) Herbig, G. H. 1966, Vistas in Astronomy, 8, 109
  * Herbig (1977) Herbig, G. H. 1977, ApJ, 217, 693
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116
  * Hoffleit (1939) Hoffleit, D. 1939, Harvard College Observatory Bulletin, 911, 41
  * Jager et al. (1998) Jager, C., Mutschke, H., & Henning, T. 1998, A&A, 332, 291
  * Kenyon et al. (1988) Kenyon, S. J., Hartmann, L., & Hewett, R. 1988, ApJ, 325, 231
  * Kenyon et al. (2000) Kenyon, S. J., Kolotilov, E. A., Ibragimov, M. A., & Mattei, J. A. 2000, ApJ, 531, 1028
  * Kóspál et al. (2021) Kóspál, Á., Cruz-Sáenz de Miera, F., White, J. A., et al. 2021, ApJS, 256, 30
  * Kraus et al. (2020) Kraus, S., Kreplin, A., Young, A. K., et al. 2020, Science, 369, 1233
  * Labdon et al. (2021) Labdon, A., Kraus, S., Davies, C. L., et al. 2021, A&A, 646, A102
  * Laor & Draine (1993) Laor, A. & Draine, B. T. 1993, ApJ, 402, 441
  * Laws et al. (2020) Laws, A. S. E., Harries, T. J., Setterholm, B. R., et al. 2020, ApJ, 888, 7
  * Liu et al. (2018) Liu, H. B., Dunham, M. M., Pascucci, I., et al. 2018, A&A, 612, A54
  * Liu et al. (2019) Liu, H. B., Mérand, A., Green, J. D., et al. 2019, ApJ, 884, 97
  * Liu et al. (2016) Liu, H. B., Takami, M., Kudo, T., et al. 2016, Science Advances, 2, e1500875
  * Liu et al. (2021) Liu, H. B., Tsai, A.-L., Chen, W. P., et al. 2021, ApJ, 923, 270
  * Liu et al. (2017) Liu, H. B., Vorobyov, E. I., Dong, R., et al. 2017, A&A, 602, A19
  * Lopez et al. (2014) Lopez, B., Lagarde, S., Jaffe, W., et al. 2014, The Messenger, 157, 5
  * Lopez et al. (2022) Lopez, B., Lagarde, S., Petrov, R. G., et al. 2022, A&A, 659, A192
  * Malbet et al. (1998) Malbet, F., Berger, J.-P., Colavita, M., et al. 1998, ApJ, 507, L149
  * Malbet et al. (2005) Malbet, F., Lachaume, R., Berger, J.-P., et al. 2005, A&A, 437, 627
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
  * Matter et al. (2016a) Matter, A., Lagarde, S., Petrov, R. G., et al. 2016a, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9907, Optical and Infrared Interferometry and Imaging V, ed. F. Malbet, M. J. Creech-Eakman, & P. G. Tuthill, 990728
  * Matter et al. (2016b) Matter, A., Lopez, B., Antonelli, P., et al. 2016b, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9907, Optical and Infrared Interferometry and Imaging V, ed. F. Malbet, M. J. Creech-Eakman, & P. G. Tuthill, 99070A
  * Min et al. (2005) Min, M., Hovenier, J. W., & de Koter, A. 2005, A&A, 432, 909
  * Monnier et al. (2009) Monnier, J. D., Tuthill, P. G., Ireland, M., et al. 2009, ApJ, 700, 491
  * Pérez et al. (2020) Pérez, S., Hales, A., Liu, H. B., et al. 2020, ApJ, 889, 59
  * Petrov et al. (2020) Petrov, R. G., Allouche, F., Matter, A., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11446, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114460L
  * Quanz et al. (2006) Quanz, S. P., Henning, T., Bouwman, J., Ratzka, T., & Leinert, C. 2006, ApJ, 648, 472
  * Reipurth & Aspin (2004) Reipurth, B. & Aspin, C. 2004, ApJ, 608, L65
  * Siwak et al. (2018) Siwak, M., Winiarski, M., Ogłoza, W., et al. 2018, A&A, 618, A79
  * Takami et al. (2018) Takami, M., Fu, G., Liu, H. B., et al. 2018, ApJ, 864, 20
  * Toon & Ackerman (1981) Toon, O. B. & Ackerman, T. P. 1981, Appl. Opt., 20, 3657
  * van Boekel et al. (2003) van Boekel, R., Waters, L. B. F. M., Dominik, C., et al. 2003, A&A, 400, L21
  * Varga et al. (2021) Varga, J., Hogerheijde, M., van Boekel, R., et al. 2021, A&A, 647, A56
  * Wachmann (1954) Wachmann, A. 1954, ZAp, 35, 74
  * Wang et al. (2004) Wang, H., Apai, D., Henning, T., & Pascucci, I. 2004, ApJ, 601, L83
  * Weingartner & Draine (2001) Weingartner, J. C. & Draine, B. T. 2001, ApJ, 548, 296
  * Woitke et al. (2016) Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103
  * Zhu et al. (2007) Zhu, Z., Hartmann, L., Calvet, N., et al. 2007, ApJ, 669, 483
  * Zhu et al. (2008) Zhu, Z., Hartmann, L., Calvet, N., et al. 2008, ApJ, 684, 1281
  * Zhu et al. (2020) Zhu, Z., Jiang, Y.-F., & Stone, J. M. 2020, MNRAS, 495, 3494

##  Appendix A Contemporary SED

The contemporary SED of FU Orionis, used as a comparison to the analytical and
radiative transfer disk models, is listed in Table 3. Here, we present our own
photometric observations from SAAO, NOT and Konkoly Observatory (Sect. 2.2),
as well as archival mid- and far-infrared photometry obtained after 2010. The
archival submillimeter continuum photometry was obtained after 2008. For the
NEOWISE photometry in particular, we use data from the latest data release,
and correct it for saturation following the instrument’s
documentation777https://wise2.ipac.caltech.edu/docs/release/neowise/expsup/sec2_1civa.html.

Table 3: FU Ori photometry, not corrected for extinction.

Band |  $F_{\nu}$ (Jy) |  $\sigma$ (Jy) | Source | Epoch | Ref.
---|---|---|---|---|---
$U$ | 0.0305 | 0.0011 | SAAO | 2021 | t.w.
$B$ | 0.1359 | 0.0014 | SAAO | … | t.w.
$V$ | 0.4061 | 0.0037 | SAAO | … | t.w.
$R_{c}$ | 0.6924 | 0.0082 | SAAO | … | t.w.
$I_{c}$ | 1.2423 | 0.0125 | SAAO | … | t.w.
$B$ | 0.1292 | 0.0014 | RC80 | … | t.w.
$V$ | 0.4338 | 0.0051 | RC80 | … | t.w.
$r^{\prime}$ | 0.6841 | 0.0083 | RC80 | … | t.w.
$i^{\prime}$ | 1.1341 | 0.0172 | RC80 | … | t.w.
$G_{\rm BP}$ | 0.3322 | 0.0011 |  Gaia DR2 | 2015 | (1)
$G$ | 0.5811 | 0.0008 |  Gaia DR2 | 2015 | (1)
$G_{\rm RP}$ | 1.1387 | 0.0047 |  Gaia DR2 | 2015 | (1)
$J$ | 2.9491 | 0.0751 | NOT | 2021 | t.w.
$H$ | 3.5083 | 0.0988 | NOT | … | t.w.
$Ks$ | 4.0031 | 0.0803 | NOT | … | t.w.
$W1$ | 3.7126 | 0.2542 | NEOWISE | 2020 | (2)
$W2$ | 3.3039 | 0.2431 | NEOWISE | … | (2)
$\rm\lambda\,(\mu m)$ |  $F_{\nu}$ (Jy) |  $\sigma$ (Jy) | Source | Epoch | Ref.
5.6 | 4.24 | 2.06 | SOFIA | 2016 | (3)
7.7 | 3.64 | 1.91 | SOFIA | … | (3)
11.1 | 4.21 | 0.00 | SOFIA | … | (3)
25.3 | 7.68 | 2.79 | SOFIA | … | (3)
31.5 | 8.74 | 2.96 | SOFIA | … | (3)
34.8 | 8.54 | 2.94 | SOFIA | … | (3)
70 | 6.5 | 0.7 | Herschel | 2012 | (4)
100 | 5.9 | 0.6 | Herschel | … | (4)
160 | 4.9 | 0.5 | Herschel | … | (4)
250 | 2.8 | 0.5 | Herschel | 2011 | (4)
350 | 2.5 | 0.6 | Herschel | … | (4)
500 | 0.9 | 0.3 | Herschel | … | (4)
853 | 0.0501 | $3e-04$ | ALMA | 2012 | (5)
1100 | 0.0350 | 0.0070 | SMA | 2008 | (6)
1300 | 0.0143 | 0.0015 | ALMA | 2017 | (7)

“t.w.” = this work; (1) Gaia Collaboration et al. (2018); (2) corrected for
saturation; (3) Green et al. (2016); (4) Green et al. (2013); (5) Hales et al.
(2015); (6) Liu et al. (2017) ; (7) Pérez et al. (2020)

##  Appendix B Observations log

The observation log for all five epochs is presented in Table 4, where DIT is
the integration time, and $\tau_{0}$ is the atmospheric coherence time. A
short description of the data quality for all five epochs is given below.

##### Epoch 1:

The detection limits for the AT configurations in low-spectral-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, $\leq 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 $\leq 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$\rm\mu$m.
This may be the result of incorrect flat-fielding, manifesting as correlated
noise features, because this spectral mode utilizes a larger portion of the
detector that has not been characterized well yet888At the time this work was
submitted to the journal, delays due to the COVID-19 pandemic had hindered the
optimization of the instrument..

##### 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 $L$-band medium-spectral resolution data are beyond the scope of this work
and will be analyzed in a future publication.

Table 4: Observing log. Epoch | Date | Band | $R$ | DIT | Configuration | Seeing | $\tau_{0}$ | Calibrator
---|---|---|---|---|---|---|---|---
|  |  | ($\lambda/\Delta\lambda$) | (ms) | array | (″) | (ms) |
1 | 2019-12-22T04:51:39 | LM | low | 111 | K0-G2-D0-J3 | 0.46 | 7.83 | HD49161
| ” | N | low | 20 | K0-G2-D0-J3 | 0.53 | 8.56 | HD49161
2 | 2020-03-15T00:23:24 | LM | low | 111 | U1-U2-U3-U4 | 1.04 | 3.14 | HD48433
| ” | N | low | 20 | U1-U2-U3-U4 | 1.10 | 2.90 | HD48433
3 | 2021-01-07T04:40:21 | LM | low | 111 | U1-U2-U3-U4 | 0.71 | 5.23 | –
| ” | N | high | 75 | U1-U2-U3-U4 | 0.71 | 5.23 | –
4 | 2021-01-08T03:33:01 | LM | low | 111 | U1-U2-U3-U4 | 0.63 | 7.87 | HD28413, HD47886
| ” | N | high | 75 | U1-U2-U3-U4 | 0.55 | 10.12 | HD28413, HD47886
5* | 2021-01-16T02:36:10 | LM | med | 3000 | A0-G1-J2-J3 | 0.84 | 7.31 | HD31767
| 2021-01-16T02:41:21 | N | high | 75 | A0-G1-J2-J3 | 0.81 | 7.32 | HD31767
5 | 2021-01-16T02:47:06 | LM | low | 111 | A0-G1-J2-J3 | 0.70 | 7.30 | HD31767
| ” | LM | low | 111 | A0-G1-J2-J3 | 0.66 | 8.42 | HD31767

##  Appendix C Extended model comparison

Further comparisons of the two radmc-3d models against the MATISSE $N$-band
data are provided in Figs. 13 and 14, which show the synthetic visibilities
and correlated fluxes with respect to wavelength, respectively. Both models
can fit the visibilities within observational uncertainties, but the ALMA-
orientation model provides better fits for the correlated fluxes.

Figure 13: Synthetic visibilities (squared) for the two radmc-3d models at two
disk orientations (ALMA, solid black; MATISSE, dashed) against the MATISSE
$N$-band data.  Figure 14: As in Fig. 13, but the correlated fluxes.