Diagnosing 0.1–10 au Scale Morphology of the FU Ori Disk Using ALMA and VLTI/GRAVITY

We report new Atacama Large Millimeter/submillimeter Array Band 3 (86–100 GHz; ∼80 mas angular resolution) and Band 4 (146–160 GHz; ∼50 mas angular resolution) observations of the dust continuum emission toward the archetypal and ongoing accretion burst young stellar object FU Ori, which simultaneously covered its companion, FU Ori S. In addition, we present near-infrared (2–2.45 μm) observations of FU Ori taken with the General Relativity Analysis via VLT InTerferometrY (GRAVITY; ∼1 mas angular resolution) instrument on the Very Large Telescope Interferometer (VLTI). We find that the emission in both FU Ori and FU Ori S at (sub)millimeter and near-infrared bands is dominated by structures inward of ∼10 au radii. We detected closure phases close to zero from FU Ori with VLTI/GRAVITY, which indicate the source is approximately centrally symmetric and therefore is likely viewed nearly face-on. Our simple model to fit the GRAVITY data shows that the inner 0.4 au radii of the FU Ori disk has a triangular spectral shape at 2–2.45 μm, which is consistent with the H2O and CO absorption features in a 10−4 M⊙ yr−1, viscously heated accretion disk. At larger (∼0.4–10 au) radii, our analysis shows that viscous heating may also explain the observed (sub)millimeter and centimeter spectral energy distribution when we assume a constant, ∼10−4 M⊙ yr−1 mass inflow rate in this region. This explains how the inner 0.4 au disk is replenished with mass at a modest rate, such that it neither depletes nor accumulates significant masses over its short dynamic timescale. Finally, we tentatively detect evidence of vertical dust settling in the inner 10 au of the FU Ori disk, but confirmation requires more complete spectral sampling in the centimeter bands.


Introduction
Understanding the physical mechanisms of protostellar accretion is fundamentally important in studies of star formation. Optical and near-infrared surveys have shown that young stellar objects (YSOs) are 10-100 times underluminous with respect to the expected luminosity from steady accretion (Kenyon & Hartmann 1995;Evans et al. 2009), which indicates that YSOs may accrete episodically (Dunham & Vorobyov 2012). If episodic accretion is a widespread phenomenon during the YSO phases, it should manifest observationally. This is consistent with the discoveries of two types of YSOs in outburst: the FU Orionis (FUor) and the EX Lupi (EXor) objects, which are characterized by a rapid large increase in their optical and infrared (OIR) brightness (for reviews, see Hartmann & Kenyon 1996;Herbig 2007;Audard et al. 2014).
FUors have outburst durations of decades to centuries (Hartmann & Kenyon 1996). During the outburst state, their optical brightness can increase by ∼4 mag or more. Models suggest that the accretion rates of these YSOs vary from 10 −7 M e yr −1 in the low (T Tauri) accretion state to 10 −4 M e yr −1 in the high (FUors) accretion state (Hartmann & Kenyon 1996). While accretion processes in quiescent T Tauri stars are generally understood as magnetospheric streams from the inner disk (e.g., Koenigl 1991;Calvet et al. 2000), how the gas and dust reservoirs immediately around the FUors are different (or altered) compared with quiescent T Tauri stars is not yet well understood. This limits our understanding of the outburst triggering mechanisms and the consequences of them.
To shed light on this issue, we have performed high angular resolution observations toward the archetypal FU Orionis object, FU Ori, using the Atacama Large Millimeter/submillimeter Array (ALMA) and the General Relativity Analysis via VLT InTerferometrY (GRAVITY) instrument of the Very Large Telescope Interferometer (VLTI). Throughout this paper, we assume the distance of FU Ori to be d∼416 pc, according to the parallax measurement published in the second data release of the Gaia space telescope (Gaia Collaboration et al. 2018). According to the prior-assisted parallax distances measurements of Bailer-Jones et al. (2018), we quote a nominal±2% distance uncertainty, which will not qualitatively affect our analysis. The observations are introduced in Section 2, and the results are presented in Section 3. By jointly analyzing these new observations with the previous (sub)millimeter observations of the ALMA, the Submillimeter Array (SMA), 15 and the NRAO 16 Karl G. Jansky Very Large Array (JVLA) and the infrared spectra taken with the Spitzer and Herschel 17 space telescopes, our interpretation and the further discussion about the physical implications are provided in Section 4. Our conclusion is given in Section 5. We refer to Berger et al. (2012) for a review of the convention and terminology for the optical and infrared interferometry technique.

Observations
We introduce the archival Spitzer and Herschel spectra and the VLTI/GRAVITY observations in Sections 2.1 and 2.2. We provide details of our ALMA observations in Section 2.3.

Spitzer and Herschel Spectra
The Herschel/PACS and SPIRE spectra were taken from the COPS-DIGIT-FOOSH (CDF) archive, a high-level data product provided to the Herschel Science Archive (see Green et al. 2016b for details). Because the source size at these wavelengths is comparable to the SPIRE beam size, and because of the lack of background subtraction, the spectra of the two modules of SPIRE instruments (SLW and SSW) are often mismatched. To resolve this discrepancy, Green et al. (2016b) apply the Semi-Extended Correction Tool in HIPE (Ott 2010;Wu et al. 2013) to calibrate the SPIRE spectra by modeling the source size. The best-fit source size, 23 5 for FU Ori, is then convolved with the beam profile of SPIRE, which is a function of wavelength (Makiwa et al. 2013). Therefore, the resulting SPIRE spectrum represents the emission from different apertures at given wavelengths, which correspond to the convolved sizes of the beam and the source size (Yang et al. 2018). For example, the aperture sizes are 29 9, 34 5, and 43 6 at 250 μm, 350 μm, and 500 μm, respectively. We refer to Section 2.2 of Green et al. (2016b) for a complete description of the data reduction. The observations were carried out at both medium and high spectral resolution, although only the medium resolution data achieved sufficient signal-to-noise ratio (S/N) for analysis. The medium spectral resolution setting covered the whole near-infrared K band with a spectral resolution of ∼500.

VLTI/GRAVITY Observations
The telescopes chosen for these observations were the medium Auxiliary Telescope (AT) configurations. This configuration includes the stations K0-G2-D0-J3, which led to baselines ranging from ∼40 to ∼100 m ( Figure 1).
The calibrator star observed concurrently was HD 38494, which is a K2 star of unknown luminosity class. Its photometric angular diameter was estimated to be θ UD =0.71±0.06 mas according to the Jean-Marie Mariotti Center Stellar Diameters Catalogue (JSDC; Bourges et al. 2017). It is nearly unresolved for our observations. The visibility of the fringes is expected to range from 0.971 to 0.996. The uncertainty of the photometric angular diameter of HD 38494 leads to a bias in the reduced data of at most 0.005 in visibility. We reduced the data using the GRAVITY pipeline (Lapeyrere et al. 2014), version 1.0.11. We note that the direct observables from VLTI/GRAVITY are normalized to the total flux. The reduced spectrum does not clearly present emission lines (Figure 1).
The Band 3 observations were carried out on 2017 November 8. The four spectral windows were centered on the sky frequencies of 86. 000, 87.863, 98.196, and 100.001 GHz. We observed the quasar J0510+1800 for absolute flux and passband calibrations, and J0547+1223 for complex gain calibrations. The Band 4 observations were carried out on 2017 November 7. The four spectral windows were centered on the sky frequencies of 146. 001, 147.863, 158.196, and 160.001 GHz. We observed the quasar J0510+1800 for absolute flux and passband calibrations, and J0536+0944 for complex gain calibrations. We manually calibrated and phase self-calibrated these data using the CASA software package (McMullin et al. 2007), version 5.4.0. When performing absolute flux scaling, we assumed that J0510+1800 has a 2.0 Jy absolute flux and a −0.30 spectral index at the reference frequency 93.015 GHz, and has a 1.6 Jy absolute flux and a −0.4 spectral index at the reference frequency 153.016 GHz. These assumptions were based on interpolating the calibrator grid survey measurements. We produced the Briggs Robust=0 weighted continuum images from line-free spectral channels using the CASA task clean. For each of the two observed bands, we created images for each of the four spectral windows separately using the multifrequency synthesis (MFS) method, setting the parameter nterm=1. The four spectral windows in each band achieved comparable rms noise levels and angular resolutions. At Band 3, the spectral window centered at 100 GHz achieved a q q maj min =0 082×0 075 (P.A.=−79°) synthesized beam and a 69 μJy beam −1 rms noise level; at Band 4, the spectral window centered at 160 GHz achieved a q q maj min =0 047×0 043 (P.A.=57°) synthesized beam and a 72 μJy beam −1 rms noise level. In each band, the synthesized beam sizes at other spectral windows are inversely proportional to their central frequencies. Figure 2 shows the reduced VLTI/GRAVITY data on FU Ori. The closure phases (CPs) are smaller than±2°.5. In addition, the overall scatter of the CP is less that 1°. This indicates that, on the spatial scales resolved by our VLTI/ GRAVITY observations, FU Ori appears approximately centrosymmetric. Overall, the squared visibilities have a fairly high level, ranging from 0.7 to 0.9. The variations of the squared visibilities with wavelength are similar for all baselines, irrespective of baseline lengths and orientations. The squared visibilities are approximately constant from wavelength 2.0 to 2.2 μm, and then drop by about 0.1 from 2.2 to 2.45 μm.

VLTI/GRAVITY Data
The fact that the differential visibility variations do not seem to depend on baseline lengths indicates that the intensity distributions may be approximated by a compact component (hereafter VLTI-compact) at the center and a more extended centrosymmetric component (hereafter VLTI-extended). 18 Including a structure that is nearly resolved out by all baselines leads to the observed <1.0 squared visibilities. In this case, the detected values of the squared visibilities depend mostly on the flux ratios of the VLTI-compact and the VLTI-extended components.
To give a qualitative sense, if we define VLTI-unresolved as having a visibility higher than 0.99 and VLTI-resolved as having a visibility amplitude less than 0.01, then for our longest, ∼100 m baseline, a uniform disk with 0.3 mas diameter is VLTI-unresolved. For our shortest, 30 m baseline, a two-dimensional Gaussian with ∼17 mas FWHM is VLTIresolved. The visibility amplitude of a 1 mas compact uniform disk ranges from 0.93 to 0.99 for the baselines ranging from 60 to 100 m. The visibility amplitude of an FWHM=8 mas Gaussian ranges from ∼0.0001 to ∼0.15. Figure 3 shows the ALMA 100 GHz (Band 3) and 160 GHz (Band 4) images. These ALMA observations detected FU Ori and FU Ori S (Reipurth & Aspin 2004;Wang et al. 2004) at high significances. However, they only marginally spatially resolved the structures. Their fluxes determined by fitting twodimensional Gaussians are summarized in Table 1. By quoting the previous JVLA observations at 29-37 GHz (Liu et al. 2017) and the ALMA observations at ∼225 GHz (Pérez et al. 2019) and at ∼346 GHz (Hales et al. 2015), the derived (sub) millimeter spectral indices (α) at various frequency ranges are summarized in Table 2. Figures 4 and 5 summarize the spectral energy distributions (SEDs) of these two protostars at wavelengths from 2 μm to 33 mm ((9-1.5)×10 5 GHz).

ALMA Data
Both FU Ori and FU Ori S show spectral index values lower than 2.0 19 at 29-37 GHz (see Liu et al. 2017); the averaged α are ∼2.5 over the frequency range of 29-100 GHz; the averaged α is approximately 2.0 from 100 to 232 GHz, and are higher than 2.0 at higher frequencies. We require multiple emission components with distinct physical properties to fit the complex submillimeter spectral slopes in the observed SEDs. Our detailed SED models for all data presented in Figures 4 and 5 are described in Section 4.

Discussion
In Section 4.1 we introduce a simple geometric model to interpret the VLTI/GRAVITY observations. In addition, we have generated simple radiative transfer models to interpret the SEDs of FU Ori and FU Ori S. In Section 4.2 we introduce how we produced the spectra for individual dust or free-free (i.e., from ionized gas) emission components in our radiative transfer model. In Section 4.3, we introduce how we integrate each of the emission components to the abstracted geometric models to reproduce the integrated SEDs, and how we optimized the model free parameters using the Markov Chain Monte Carlo (MCMC) method. We discuss the physical implications of our models in Section 4.4.

Interpreting VLTI/GRAVITY Data
The fact that the observed squared visibilities in the VLTI/ GRAVITY data vary with wavelength ( Figure 2, right panel) implies that the flux ratio of the VLTI-compact and the VLTIextended components has a wavelength dependence.
We can quantify this dependence by fitting the data. For simplicity, we assumed that the VLTI-compact and the VLTIextended components have constant sizes over the wavelength range covered by the VLTI/GRAVITY observations. In addition, we assumed that the VLTI-compact component is a uniform disk, while the VLTI-extended component is a twodimensional Gaussian of which the aspect ratio is ∼1. We then performed chi-squared fits to determine the sizes of the two 18 Note that the two infrared emission components resolved by the VLTI/ GRAVITY observations are both more compact than what were detected by JVLA and ALMA. Our terminology is to distinguish them from the spatially more extended (sub)millimeter and centimeter sources. 19 The spectral index α was measured assuming that the flux F ν around a reference frequency ν 0 can be expressed as F ν =F 0 (ν/ν 0 ) α . components, and to determine the flux ratios as a linear interpolation between eight equally spaced wavelengths between 2.0 and 2.45 μm (R∼40). We tried various combinations of sizes, allowing the two components to be unresolved, partially resolved, or fully resolved. Given that the observed closure phases are less than 2°, to avoid overfitting, the two components were concentric in most trials. However, in one of the trials, we also explored how much their centers can deviate.
Our best-fit geometric models are summarized in Table 3. The model in best agreement with the data is a uniform disk of ∼1 mas in diameter, and a two-dimensional Gaussian with FWHM∼8 mas (i.e., solid angle ∼1.7×10 −15 sr). The fit is further improved if the VLTI-extended component is slightly offset to the northeast by ∼1.2 mas ( Figure 6, right panel).
Assuming that the VLTI/GRAVITY detections arose predominantly from the circumstellar disk, this spatial offset can be interpreted either as a disk that is geometrically thick (e.g., flared) and is slightly inclined (e.g., Figure 5 of Zhu et al. 2008), or as a disk that includes substructures or companions (e.g., Malbet et al. 2005). Using the flux ratios from Table 4, we were able to reproduce the observed slight closure phase signal that increases with wavelength, and the wavelengthdependent variations of the squared visibilities ( Figure 2). By implementing an absolute flux scaling, the spectral shape of the VLTI-compact component in our best-fit model ( Figure 6) appears fully consistent with the viscous accretion disk model of Calvet et al. (1991), which assumed a ∼10 −4 M e yr −1 mass accretion rate. The triangular shape of the spectrum presented in Figure 6, following the framework of Calvet et al. (1991), is due to the absorption of the water band and the first-overtone . Visibility data for FU Ori taken with VLTI/GRAVITY. The left and right panels show the closure phase (in degree units) and the squared visibility (normalized to 1) as function of wavelength, respectively. The red curves are our best-fit model to the VLTI/GRAVITY data, which is composed of the marginally spatially resolved compact and extended sources (i.e., Model 4 in Table 3; see Figure 6 for more details of the model). For each closure phase, the names of the involved telescopes are labeled, which can be referenced from Figure 1. In each panel, the name, length, and orientation of the baseline are labeled. The vertical red line segments indicate the wavelengths of the Bracket γ transition of hydrogen (2.16612 μm) and the CO band heads (2.2935, 2.3227, 2.3535, 2.3829, 2.4142, 2.4461 μm).
vibration-rotation CO band against the bright continuum emission from the viscously heated midplane. Calvet et al. (1991) suggested that these absorption features are predominantly produced at radii of 0.1-0.3 au (i.e., 0.96±0.48 mas angular diameter assuming d∼416 pc) around the host protostar, which is consistent with the angular sizes in our model fits.

Individual Emission Components
To evaluate the emission properties of dust, we utilized the DSHARP dust optical constants published in Birnstiel et al. (2018). For simplicity, we assumed a constant 170 K water ice sublimation temperature (Pollack et al. 1994). Therefore, we adopted the default DSHARP optical constants for dust emission sources that are cooler than 170 K, and adopted the ice-free optical constants for those that are warmer. Given that the physical conditions of the observed sources (in particular, FU Ori) may be out of equilibrium in various ways, it is not possible for us to evaluate the detailed form of the grain size distribution function from first principles. Therefore, when evaluating the size-averaged dust absorption (k n abs ) and effective scattering (k n sca,eff ) opacities, we simply assumed the typical power-law grain size distribution with a power-law index q=3.5, the minimum grain size a min =10 −4 mm, and the maximum grain size a max . Before considering mutual obscuration, the SEDs of individual dust emission components ( n F dust ) were evaluated based on the analytic radiative transfer solutions published in Birnstiel et al. (2018). Motivated by the small (or negligible) angular offsets of the unresolved and resolved components in the VLTI/GRAVITY models (Table 3; see discussion in Section 4.1), we considered all dust slabs to be approximately face-on. We note that introducing inclinations of the dust slabs will not change the conclusion from our radiative transfer models qualitatively. However, this would increase the total number of free parameters. Figure 7 shows examples of the SEDs produced for the dust slabs with a dust column density of 50 g cm −2 , temperature of 100 K, solid angle of 1 arcsec 2 , and maximum grain sizes of 0.002 mm (top panel), 0.2 mm (middle panel), and 2 mm  (bottom panel). In the low-frequency, low optical depth regime, the SEDs deviate from the blackbody emission model (i.e., Planck function) as dust grains cannot emit/absorb efficiently at wavelengths that are much longer than a max . In addition, for a max =0.2 mm or 2 mm, the SED deviates from a blackbody curve at higher frequencies. As frequency increases, the spectral indices fall below a blackbody curve, and then become steeper; thus, the flux in this frequency regime is below that of a blackbody curve. The dust slabs are optically thick in this frequency regime, and the effects of dust (self-) scattering are not necessarily negligible. We attribute this deviation from blackbody emission in the high-frequency regime to the frequency variations of albedo, which was addressed in detail in Liu (2019) and Zhu et al. (2019). For example, in Figure 7, the SED of the a max =2 mm dust slab shows a rather flat spectral index at ∼20-50 GHz, which is because the albedo increases with frequency; the spectral index is steepened at ∼50-1000 GHz because the albedo decreases with frequency.
Following Mezger & Henderson (1967) and Keto (2003), we approximated the optical depth of the free-free emission components t n ff by where EM is the emission measure defined as EM= ò n dℓ e 2 , with n e being the electron number density, and ℓ is the linear size scale of the free-free emission component along the line of sight. Fluxes of individual free-free emission components were evaluated based on where Ω ff is the solid angle of the free-free emission component, and B ν (T) is the Planck blackbody function.

Abstracted Geometric Model and Integrated SEDs
The overall fluxes (F ν ) of FU Ori and FU Ori S were determined from the following formulation: where n F i is the flux of the dust or free-free emission component i, and t n i j , is the optical depth of the emission component j to obscure the emission component i. The abstracted geometric information is provided by t n i j , . We chose this approach instead of fitting analytical solutions of (gaseous) disks because dusty protoplanetary disks are commonly composed of substructures (e.g., rings, crescent, etc.). Our SED fitting procedure for the spatially unresolved target sources effectively decomposed them into substructures of certain projected areas, but without explicitly constraining the shapes. In this work we considered a simple implementation, such that t = We tried fitting the observed SEDs with the least number of emission components to minimize the total number of free parameters, for various configurations of τ i,j . The free parameters and the configurations of τ i,j were varied interactively, informed by the results of previous trials. Our interactive fits focused on matching the interferometric data. Nevertheless, we found that once a good fit for the interferometric data was achieved, the infrared spectrum predicted from the model is also very close to the Spitzer and Herschel observations.  Green et al. 2006Green et al. , 2013Green et al. , 2016a, and the fluxes of the au scales structures around FU Ori taken with the VLTI/GRAVITY (i.e., the "Extended" column of Table 4). Black lines show our model of the combined fluxes of these two sources. Lines with other colors are the fluxes of individual dust or free-free emission components in our model (see also Figure 5; see Table 5). We assumed that the envelope component was only detectable from Herschel and was resolved out by any of our interferometric observations. Model components that are labeled but cannot be found in the right panel are due to that their fluxes are below the plotted range.
After we obtained an approximate fit, we used MCMC to simultaneously optimize all free parameters (i.e., all the parameters in Table 5 except the column of overall dust masses). We assumed flat priors, which permitted each parameter to vary from half of its initial value to two times the initial value. To prevent the MCMC routine from sampling large unlikely portions of parameter space, we provided an additional constraint from the FU Ori S 9 GHz nondetection ( Figure 5). We forced the logarithmic likelihood to be negative infinity when the integrated flux of FU Ori S at 9 GHz is higher than three times the rms noise of the observations, a condition in the likelihood function to force rejecting such MCMC samples.
The MCMC fittings were initialized with 84 walkers with 1500 iterative steps each; in the end, the results from the first 500 steps were discarded. The Herschel and Spitzer data have very good S/Ns, such that their contribution to the likelihood largely outweighed the contribution from interferometric data. To avoid overfitting the Herschel and Spitzer data without achieving a good fit for the interferometric data, we needed to reduce the weight of Herschel and Spitzer data. This is implemented by artificially assigning the flux errors of the Herschel/SPIRE, Herschel/PACS, and Spitzer/IRS to be 1000, 10, and 1 times the detected fluxes. We have monitored how the likelihood evolved over the MCMC iterations to make sure that the contribution of the Herschel and Spitzer data are on the same order with the rest of the data. We note that during the steps of MCMC, some dust emission components may switch from being based on the default DSHARP optical constants to being based on the ice-free optical constants (i.e., the walkers "walked" from below to above the 170 K dust temperature). Because of this midroutine shift, it is very difficult to implement fitting methods other than MCMC.  Hales et al. 2015). Throughout this paper we assumed a nominal 10% error for the 346 GHz fluxes of FU Ori and FU Ori S since they were not clearly separated in the previous ALMA image due to the limited angular resolution. The inverted triangle shows the 3σ upper limit for FU Ori S at 9 GHz. Colored lines show fluxes of our model for each of these two resolved sources (see Table 5; for the labels see Figure 4). For both sources, blue lines show the free-free emission component; cyan lines show the dense and hot inner disks of a few au scales; red lines show the outer disks on few tens of au scales; the light green line shows a spatially compact dust component that is enclosing the hot inner disk of FU Ori and has a lower dust temperature than that of the hot inner disk.
We found at least four dust emission components are required to fit the JVLA, ALMA, and VLTI/GRAVITY data points for FU Ori (for more discussion see Section 4.4). Therefore, we also adopted a four emission component fit for FU Ori S.
In addition, we included an extended common-envelope component that is required to fit the far-infrared fluxes detected by the Herschel space telescope. The common-envelope component has an extended angular scale, such that it is filtered out by all interferometric observations presented in this work. We note that the envelope component must be included since the previous Herschel photometric imaging observations have spatially resolved complicated structures on subparsec scales, which connect to FU Ori and FU Ori S (Green et al. 2013). We used a simplified parametric model for the envelope ( Table 5); detailed modeling of the envelope is beyond the scope of our present study.
During the iterations, we found that we can obtain a reasonably good fit to the JVLA and ALMA data of FU Ori S by including only two dust emission components and a freefree emission component. We tentatively assign one additional, ∼140 K dust emission component to FU Ori S to better explain the Herschel or Spitzer spectra at (5-10)×10 3 GHz. Qualitatively, the fact that we need this extra component to explain the mid-far-infrared spectra indicates that the dust components in our models are not isothermal. We considered whether each dust component in our models should be allowed to have a small (e.g., 10%-20%) temperature range, which could yield better fits to the infrared spectra. However, our ability to measure any temperature variation is fundamentally limited by the wavelength-dependent aperture used to extract the Herschel/SPIRE spectra. This is because we applied semiextended source correction to align the two SPIRE modules (Wu et al. 2013;Green et al. 2016b). Residual artifacts from this process can bias our SED fits, although we have mitigated this by artificially lowering the weighting of the Herschel data. Nevertheless, we do not consider it to be meaningful to use a further detailed parameterization for dust temperature profiles to improve the fittings to our infrared spectra.
Parameters of our best-fit model are summarized in Table 5. The SEDs of the individual components after incorporating the effect of obscuration, and the integrated SEDs from all emission components, are presented in Figures 4 and 5.

Model Parameters and Their Physical Implications
In Sections 4.4.1 and 4.4.2 we discuss qualitatively the fitting parameters for FU Ori and FU Ori S. The overall geometric picture and the physical implications are discussed in Section 4.4.3.

FU Ori Model
Qualitatively, the fact that the spectral index α of FU Ori is ∼2 at 86-232 GHz and is ∼3 at frequencies higher than 232 GHz ( Table 2) indicates that fluxes at intermediate frequency (e.g., ∼150 GHz) are a mix of one emission component with α>3 and the other emission source with α<2. The α>3 component (hereafter FUOri_dust3) becomes more prominent at higher frequencies, while the α<2 source becomes more prominent at lower frequencies; the observed spectral indices at the intermediate frequencies are weighted averages from these two components.
In order to fit the JVLA data at 29-33 GHz and the ALMA data at 86-160 GHz (Table 1; Figure 5), we need to realize the α<2 source by combining at least two dust components: a ∼400 K component with high dust column density and a max ∼2 mm (hereafter FUOri_dust1), obscured by a ∼130 K component with modest dust column density and a max ∼0.2 mm (hereafter FUOri_dust2). The dust temperature of FUOri_dust1 is consistent with the high dust brightness temperature observed at ∼33 GHz (Liu et al. 2017). The FUOri_dust1 component, whether or not it is mixed with some free-free emission, naturally explains the <2.0 spectral index at 29-37 GHz, and the ∼2.5 spectral index in between 29-100 GHz ( Table 2), due to the albedo effect introduced in Section 4.2 (see also Figure 7). Being obscured by FUOr-i_dust2 makes the spectral index of FUOri_dust1 much lower than 2.0 at ∼100-150 GHz. FUOri_dust2 has a max ∼0.2 mm because this a max value yields a high albedo at ∼200 GHz. In this case, FUOri_dust2, which is optically thick at ∼200 GHz and has a rather flat spectral distribution at this frequency, can scatter off the emission from FUOri_dust1 without contributing much of the emission. This is critical to fit the steeper spectral indices at 232-345 GHz including the optically thin dust component FUOri_dust3. If one or both of FUOri_dust1 and FUOri_dust2 contributes more emission at ∼200 GHz, it becomes impossible to reproduce the steep spectral index observed at 232-345 GHz.
The a max of FUOri_dust3 is not well constrained by the data presented in this paper. Consistent with previous reports of near-infrared scattered light (e.g., Liu et al. 2016;Takami et al. 2018), we presume that the a max of FUOri_dust3 is on the order of ∼2 μm. We cannot accurately determine the dust masses of these two components due to the uncertainties of the dust mass opacities. Note. Δx e , Δy e , and c r 2 are the horizontal and vertical offsets of the VLTI-extended component with respect to the phase referencing center, and the chi-square of the fittings, respectively.
Finally, by introducing another FUOri_dust4 component we can simultaneously fit the resolved VLTI-extended component ( Table 3) in the VLTI/GRAVITY data and the higherfrequency part of the Spitzer spectrum (Figure 4). The dust temperature of FUOri_dust4 (∼700 K) is higher than that of FUOri_dust1 (∼400 K), indicating that FUOri_dust4 is likely the closest component to the host protostar. The a max value of FUOri_dust4 is not constrained by the observations presented in this paper. The resolved VLTI-compact component ( Table 3) is hotter than the dust sublimation temperature, and therefore is not considered in our dust models. The thermal radiation from the VLTI-compact component may heat the VLTI-extended component (see Zhu et al. 2007). We note that an excellent fit to the 9 GHz observations with the free-free emission component was not necessary, because that particular measurement was impacted by poorly characterized delay errors, and is rather uncertain (see Liu et al. 2017;Pérez et al. 2019).
Emission at 9 GHz may also include a nonthermal emission contribution, which we do not have sufficient data to constrain.

FU Ori S Model
The spectral index of FU Ori S is ∼2 over a broad frequency range of 86-346 GHz ( Table 2). An optically thick dust component (FUOriS_dust1) with a max ∼0.2 mm can explain the slightly smaller than 2.0 spectral index at ∼90 GHz. Mixing FUOriS_dust1 with an optically thinner dust component (FUOriS_dust3) and a free-free emission component can better fit the observations at 29-37 GHz and at 346 GHz ( Figure 5). To reproduce the resolved SEDs for FU Ori S, there is no need of assuming mutual obscurations of the emission components since its spectral index at 232-345 GHz is not as steep as that of FU Ori (Table 5).

Outbursting versus Quiescent Disks?
By assuming a geometrically thin, axisymmetric, Keplerian, hot inner disk around the center of FU Ori, Calvet et al. (1991) argued that the observed CO line widths at near-infrared bands are consistent with an inclination of ∼20°-60°. Based on analyzing the squared visibilities from near-and mid-infrared interferometric observations, Malbet et al. (2005), Zhu et al. (2008), and Quanz et al. (2006) suggested that the inclination of the disk is ∼50°. However, being an accretion outburst object, FU Ori may not be in equilibrium. The assumptions of geometrically thin, axisymmetric, and the Keplerian velocity fields all need to be tested by resolved observations. A great advantage of VLTI/GRAVITY over the previous generations of near-or mid-infrared interferometry is that we can anchor the hypothesis of a small inclination angle based on the resolved small closure phases (Figure 2; Section 4.1). Thus, on Figure 6. Results of model fits to the VLTI/GRAVITY data. The left panel shows spectral energy distributions (SEDs) of the unresolved (VLTI-compact) and resolved (VLTI-extended) components derived from our VLTI/GRAVITY model (Table 3; see also Table 4), and observed total SED quoted from Calvet et al. (1991) and Mould et al. (1978). The right panel plots Model 4 from Table 3, where the VLTI-compact and VLTI-extended components are shown in dark blue and red, respectively. The horizontal and vertical axes of the right panels are in units of milliarcseconds. the spatial scales of a few au, the morphology and the gas kinematics of the FU Ori disk may be more complicated than previously assumed, which can be further resolved by future observations with better uv coverage. If we assume an approximately face-on projection of FU Ori (and FU Ori S), then the abstracted geometry we introduced during the SED fits (Section 4.3) follows the picture in Figure 8. Overall, we interpret the observational data for FU Ori as follows: a >1000 K hot inner disk at 0.1-0.3 au radii (0.24-0.72 mas) that produces water and CO absorption features at near-infrared bands (unresolved by VLTI/GRAV-ITY); a ∼700 K, not very optically thick dust component with ∼3 au radius (∼7 mas; FUOri_dust4, resolved by VLTI/ GRAVITY); a very optically thick and a modestly optically thick dust component with up to ∼10 au radii (FUOri_dust1,2); an optically thin, cooler dust component on tens of au scales (FUOri_dust3), and some free-free emission. Assuming that the gas-to-dust mass ratio is ∼100, the mass surface density of the component FUOri_dust1 (Table 5) is reasonably consistent with the hydrodynamic simulations presented in Zhu et al. (2010) and Bae et al. (2014). However, the vertical thermal profile of FU Ori in its inner 10 au region appears opposite to the typical model of passive disks dominated by radiative heating (e.g., some T Tauri disks; Kama et al. 2009;Tapia & Lizano 2017). The role of viscous heating in dust thermal dynamics is presently uncertain as it is difficult to observationally constrain gas volume density and viscosity at the disk midplane.
FU Ori S can be interpreted as an optically thick dust component with ∼10 au radius (FUOriS_dust1) and an optically thin, cooler dust component on tens of au scales (FUOriS_dust3), potentially with contributions from free-free emission.
The qualitative difference between the inner ∼10 au region of the FU Ori and the FU Ori S disks, in particular the thermal profile, may be related to the thermal and magnetorotational instabilities triggered during the outburst of FU Ori. For a physical picture, we refer to Figures 1 and 2 of Zhu et al. (2009). The two-dimensional hydrodynamic simulations of Zhu et al. (2009) have demonstrated that during outburst, within 10 au of the protostar, viscous energy dissipation is sufficient to heat gas at the disk midplane to a considerably higher temperature than the gas at the disk surface. It is not yet very clear to us whether or not this can also explain the vertical dust temperature profile of the FU Ori disk.
For an order-of-magnitude estimate, we compare our fits of dust temperatures (Table 5) with the simplest analytic models of the dust temperature profiles (see Chiang & Goldreich 1997 and references therein) in Figure 9. We caution that many of the underlying assumptions of the simplest analytic models (e.g., axisymmetry, steady or stationary disk, etc.) contradict the observations of the FU Ori disk, which is likely asymmetric and may be undergoing instabilities over a broad spatial scale. Such a comparison can serve as a sanity check for whether or not a certain heating mechanism can potentially provide a sufficiently high heating rate to explain the observed dust radiation temperatures. However, the comparison is not yet sufficient for verifying or strictly falsifying a certain scenario.
To assess how the FU Ori disk can be heated due to viscous dissipation, we quoted the effective radiation temperature profile of a steady-state viscous disk ( ) T r viscous assuming that the disk is very dense and is optically thick such that dust and gas can be thermalized via inelastic collisions, and that there is no radiative heating (Pringle 1981): where G is the gravitational constant, M * =0.5 M e is the assumed host protostellar mass,Ṁ is the mass accretion rate that we assumed to be 10 −8 , 10 −6 , and 10 −4 M e yr −1 , σ is the Stephen-Boltzmann constant, and R * is the stellar radius that we assumed to be 2 R e . These profiles, which may be regarded as lower limits to the dust temperature in viscous disks, are presented as the blue lines in Figure 9.
To assess how the FU Ori S disk can be heated due to protostellar irradiation, we scaled the approximate solutions for the surface ( ( ) T r s ) and interior ( ( ) T r i ) dust temperature profiles of a radiative equilibrium disk (see Equations (11) and (4(a)) in Chiang & Goldreich 1997) according to the total protostellar luminosity. The upper and lower bounds of the yellow filled area are shown with respect to ( ) T r s and ( ) T r i in Figure 9. We assumed an effective stellar temperature T * =4000 K and stellar radius R * =2 R e , typical for T Tauri stars. We note that dust g cm −2 isothermal (100 K) dust slab of 1 arcsec 2 angular size, based on the analytic radiative transfer solution and the dust opacities published in Birnstiel et al. (2018). The solid gray line shows the case of blackbody emission. due to the Stefan-Boltzmann law the dust temperatures have a weak dependence on the protostellar luminosity.
In Figure 9, we also overplotted our fits of dust components (see Table 5). The inner and outer radii of these dust components were estimated to be 1% and 100% of their solid angle, assuming a circular geometry in a face-on projection. We found that it is plausible to interpret the observed radiation temperature of FUOri_dust1 based on ( ) T r viscous given the ∼10 −4 M e yr −1 accretion rate of FU Ori. If this is the case, a higher dust temperature at the disk midplane than at the surface can be expected, which explains why FUOri_dust1 has a higher temperature than FUOri_dust2 (Figure 8). Moreover, this explains how the 0.1-0.3 au scale hot inner disk with a 10 −4 M e yr −1 mass accretion rate (see Section 4.1) is being replenished by the up to ∼10 au scale gas reservoir at a modest rate, such that the hot inner disk neither becomes depleted nor accumulates mass over a short timescale. This may explain the relatively stable mid-infrared and (sub) millimeter fluxes in the previous monitoring observations (Green et al. 2016a;Liu et al. 2018).
The optically thinner components FUOri_dust3, FUOriS_dust3, FUOri_dust4, and FUOriS_dust4 are likely dominated by radiative heating. Radiative heating alone can reasonably explain the observed temperature distributions from FU Ori S.
The comparisons in Figure 9 are uncertain since the accretion rates of FU Ori and FU Ori S are not necessarily constant over all radii. In addition, FU Ori is unlikely to be in equilibrium, and it is not trivial to accurately estimate the disk scale-height and thus the radiative heating. Moreover, these comparisons have ignored other mechanical processes that can potentially be important in asymmetric or unstable systems (e.g., shocks, adiabatic compression, etc.; Sakai et al. 2014;Dong et al. 2016). More realistic considerations of dust and gas dynamics, grain growth, and dust heating/cooling would provide better comparison. We additionally hypothesize that, during the outburst, the inner 0.1-10 au disk may expand significantly in the vertical direction, may be partly thermally ionized, and some dust may be sublimated. The morphology of the 0.1-10 au disk may also become porous due to accretion and instabilities, allowing dust to be radiatively heated close to the disk midplane at a relatively large range of radii. Finally, why might we have detected millimeter-sized a max from FU Ori (i.e., from component FUOri_dust1) but not from FU Ori S? A tentative hypothesis is that at the quiescent stage, dust grains of millimeter or larger sizes may either be radially trapped in regions too small in a projected area to be detected by observations (e.g., Vorobyov et al. 2018;Okuzumi & Tazaki 2019), or areas that are fully obscured due to a combination of very high optical depth and the vertical dust settling. These mechanisms may be particularly efficient if the inner few au regions are effectively dead zones with a negligible ionization fraction during the quiescent stage. The instabilities during the outburst may help radially and vertically mix dust grains of various sizes, which make the millimetersized grains more easily detectable. That we find tentative evidence of vertical dust settling by comparing the a max values of FUOri_dust1 and FUOri_dust2 may also be because viscous heating is more efficient in heating the vertically settled grown dust from the midplane. This may be further tested by a systematic comparison of the (sub)millimeter and radio spectral indices of the inner disks of outbursting and the quiescent T Tauri sources.

Conclusion
We have analyzed unpublished archival data from the Guaranteed Time Observations of VLTI/GRAVITY at the near-infrared K band (2-2.45 μm) toward the archetypal accretion outburst YSO, FU Ori. In addition, we have performed high angular resolution ALMA observations at 86-100 GHz and 146-160 GHz bands, which simultaneously covered FU Ori and its companion, FU Ori S.
The observed small closure phases by VLTI/GRAVITY indicate that the FU Ori disk may be approximately face-on. In addition, by comparing with the squared visibilities resolved by previous generation near-and mid-infrared interferometry, we found that the inner few au region of FU Ori may not be simply an axisymmetric, Keplerian rotating thin disk. Instead, it may have a more complicated morphology, which may be related to the instabilities that occurred during the accretion outbursts. Combined analysis of all existing ALMA, SMA, and JVLA observations along with Spitzer and Herschel infrared spectra also points to an unconventional vertical dust thermal profile in the inner ∼10 au region of the FU Ori disks. This consistently suggests a complicated disk morphology in comparison to a quiescent T Tauri disk. The observed thermal profiles in the inner ∼10 au region may be explained by a viscously heated disk of which the mass inflow rate is ∼10 −4 M e yr −1 , which can explain how the 0.1-0.3 scales hot inner disk detected from infrared observations is being replenished.
This paper is based on data obtained from the ESO Science Archive Facility under request number AMERAND384481. This paper makes use of the following ALMA data: ADS/ JAO.ALMA #2011.0.00548.S, #2016.1.01228.S, and #2017.1.00388.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together  Figures 4 and 5. The shapes of individual components also do not have strict physical meanings since they were not very well spatially resolved by the observations presented in this paper. For FU Ori, a 1 mas angle corresponds to a spatial scale of 0.416 au. Figure 9. Comparison of the fits of dust components with the analytic models of dust temperature profiles due to viscous or radiative heating (see Section 4.4.3). The filled and hatched rectangles show the dust components in FU Ori and FU Ori S as listed in Table 5. The blue lines are the effective temperature profiles of the steady-state viscous disks (stellar mass M * =0.5 M e ) with accretion rates˙= --M 10 , 10 8 6 , and 10 −4 M e yr −1 , in the absence of radiative heating. The yellow filled region is bounded by the surface and interior temperature profiles of a radiative equilibrium disk illuminated by a protostar with a 2 R e radius and an effective 4000 K temperature (see Chiang & Goldreich 1997).