Resolving the FU Orionis System with ALMA: Interacting Twin Disks?

FU Orionis objects are low-mass pre-main sequence stars characterized by dramatic outbursts several magnitudes in brightness. These outbursts are linked to episodic accretion events in which stars gain a significant portion of their mass. The physical processes behind these accretion events are not yet well understood. The archetypal FU Ori system, FU Orionis, is composed of two young stars with detected gas and dust emission. The continuum emitting regions have not been resolved until now. Here, we present 1.3 mm observations of the FU Ori binary system using the Atacama Large Millimeter/submillimeter Array. The disks are resolved at 40 mas resolution. Radiative transfer modeling shows that the emission from FU Ori north (primary) is consistent with a dust disk with a characteristic radius of ∼11 au. The ratio between the major and minor axes shows that the inclination of the disk is ∼37°. FU Ori south is consistent with a dust disk of similar inclination and size. Assuming the binary orbit shares the same inclination angle as the disks, the deprojected distance between the north and south components is 0.″6, i.e., ∼250 au. Maps of 12CO emission show a complex kinematic environment with signature disk rotation at the location of the northern component, and also (to a lesser extent) for FU Ori south. The revised disk geometry allows us to update FU Ori accretion models, yielding a stellar mass and mass accretion rate of FU Ori north of 0.6 M⊙ and 3.8 × 10−5 M⊙ yr−1, respectively.


Introduction
Protostars build a significant fraction of their final mass through an accretion disk fed by a surrounding envelope. As this envelope dissipates, accretion rates tend to drop below 10 −7 or even 10 −8 M e yr −1 . While most young stellar objects have luminosities that are significantly lower than expected from steady protostellar disk accretion (the so-called luminosity problem; Kenyon et al. 1990;Evans et al. 2009), some exhibit episodes of high activity, increasing their optical brightness by several orders of magnitude on timescales as short as one year (see Audard et al. 2014 for a review). Accretion at this stage can reach up to 10 −4 M e yr −1 . A possible solution to the socalled luminosity problem is that most protostars also go through an "FUor" phase (named after the class prototype FU Orionis) early in their lifetime. This would imply that a significant fraction of the mass of a young protostar is accreted in an episodic way, with prolonged periods of low accretion Vorobyov & Basu 2010;Bae et al. 2014).
An Early Science campaign to characterize FUor sources with the Atacama Large Millimeter/submillimeter Array (ALMA) shows that they are compact in radio continuum . Their disk emission has a smaller characteristic radii for a given disk mass than TTauri disks and they have disk sizes and masses that are more similar to those of Class I than to Class II sources (e.g., Sheehan & Eisner 2017;Tobin et al. 2018). FUors appear surrounded by prominent outflows of molecular gas with hourglass-like morphologies (e.g., V2775 Ori, Zurlo et al. 2017). These outflows can sometimes display wide opening angles (e.g., HBC 494, Ruíz-Rodríguez et al. 2017a, 2017b, and even show notable misalignments with respect to the orientation of the dust disk (e.g., V1647 Ori, Principe et al. 2018).
Episodic accretion and its implications for star and planet formation are not well understood. Several physical processes have been proposed to explain such dramatic accretion events. The most favored mechanisms include disk fragmentation and the subsequent inward migration of the fragments (Vorobyov & Basu 2010), gravitational instability, and magneto-rotational instabilities (Armitage et al. 2001;Zhu et al. , 2010Martin & Lubow 2011;Bae et al. 2014), among others. These models predict the presence of distinct features (such as clumps and spiral arms) that should be detectable with current observational capabilities, although these have not been found on high resolution images of FUor disks (e.g., V883 Ori, Cieza et al. 2016).
Another possibility is that, since these young disks dwell in crowded star-forming regions, their outbursts may be connected to, or driven by, binary interactions (Bonnell & Bastien 1992;Reipurth & Aspin 2004), or star-disk flybies (e.g., Pfalzner 2008;Cuello et al. 2019) as stars are common and disks have large cross-sections, and possibly even disk-disk interactions (e.g., Muñoz et al. 2015). A few FUor-like protostars have now been identified as having binary components (e.g., L1551 IRS5, RNO 1B/C, AR 6A/B; Pueyo et al. 2012 and references therein), feeding the discussion on whether the rapid accretion could be caused by close-companion interactions. Perturbers such as planets have also been proposed as an explanation for the quick rise in the outburst lightcurves of FU Ori and V1057 Cyg (Clarke et al. 2005).
ALMA Early Science observations of the FUOri binary system show that each binary component is surrounded by compact continuum emission (unresolved at ∼0 6 resolution) and widespread molecular emission (Hales et al. 2015). Using the Karl G. Jansky Very Large Array (VLA) in its most extended configuration, Liu et al. (2017) separate the two binary disks, but the individual disks remain unresolved at 29GHz with a 0 07 resolving beam. Combining these observations with lower (8-10 GHz) and higher frequency data, Liu et al. (2017) show that the 8-346GHz spectral energy distributions (SEDs) of FU Ori north (or simply "FU Ori", the primary) and FU Ori south (secondary) cannot be reproduced by constant spectral indices, suggesting a multi-component circumstellar environment (bearing in mind that these sources are time variable). Liu et al. (2017) also show that the observations can be explained by a circumstellar environment consisting of free-free emission from ionized gas, and thermal emission from two dust components, a compact optically thick disk and an optically thin region. On the other hand, Zhu et al. (2007) show that the optical-infrared (IR) spectrum of FUOri requires a hot disk with a mass accretion of 2×10 −4 M e yr −1 .
In this work, we present 40mas resolution ALMA observations of the FUOri system which resolve the continuum emitting regions around each stellar component, providing the first direct measure of their sizes at millimeter wavelengths (Section 3). An analysis of the local kinematics around the FUOri components using molecular line observations is presented in Section 4. Implications are discussed in Section 5, while Section 6 presents our concluding remarks. Throughout this work we use a distance to FUOri of 416±9pc (Gaia Collaboration et al. 2018, DR2), consistent with the distance to the Orion Nebula Cluster from a parallax of non-thermal radio emission (Menten et al. 2007).

Observations and Data Reduction
FUOri was observed during ALMA Cycle 4 on 2017 September 6 using ALMA's 12 m array in extended (C40-8) configuration, with 42 antennas with baselines ranging from 41m to 7.5 km. This configuration yields an angular resolution of ∼0 05 (∼21 au) and a maximum recoverable scale of ∼1 05 (∼437 au) at 1.3 mm. The median precipitable water vapor column was 0.35 mm, with clear sky conditions and an average phase rms of 28°. 6 after water vapor correction.
The primary flux calibrator was the quasar J0423−0120, while J0510+1800 was used as a bandpass calibrator. The quasar J0532+0732, located 3°. 5 from the target was observed as a phase calibrator, by alternating with the science target every 54 seconds to calibrate the time-dependent variation of the complex gains. The total time spent on source was 39.8 minutes. A secondary phase calibrator (J0551+0829), located 4°. 7 from the phase calibrator, was observed regularly as a check source to assess the quality of the phase transfer. The ALMA correlator was configured in frequency division mode. Two spectral windows with 1.875GHz bandwidth were set up for detecting the dust continuum, centered at 232.005 and 218.505GHz, respectively. The 12 CO(2-1), 13 CO(2-1) and C 18 O(2-1) transitions of carbon monoxide were targeted by configuring three spectral windows at rest frequencies of 230.538GHz, 220.399GHz and 219.560GHz, respectively. The spectral resolution for the line observations was 122.070kHz (∼0.2 km s −1 ). All data were calibrated by the ALMA staff using the ALMA Science Pipeline version r39732 in the CASA package (McMullin et al. 2007).

ALMA Continuum Imaging
Imaging of the 1.3 mm continuum emission was performed using the CLEAN task in CASA with uniform weighting, resulting in a synthesized beam size of 60×42mas and position angle (PA) of 68°.8 (Figure 1). Self-calibration of the data was performed to improve coherence. A single iteration of phase-only (i.e., no amplitude) calibration resulted in an improvement in the resulting signal-to-noise ratio (S/N) by a factor of 1.7 (30 μJy rms) (compared to the S/N of a CLEAN image before selfcalibration). Another iteration of self-calibration was found to only marginally improve the image quality.
As visibilities are spatially integrated quantities, one of the components needs to be subtracted to inspect the visibilities of the other. This is achieved by subtracting the CLEAN model of one of the components from the visibility data. There are no indications of substructures in either disk, or of emission in between the stellar components in continuum. The observed visibilities of each component (right panel in Figure 1) show the decreasing profile characteristic of extended sources (e.g., Wilson et al. 2009).
By imaging each continuum spectral window separately (218 and 232 GHz), we derive in-band spectral indices (using the full bandwidth of each spectral window, i.e., 1.875 GHz) of 2.1 for both FUOri and FUOri south. This is consistent with the spectral indexes derived in recent SED fitting (Liu et al. 2019). The ALMA data confirm both disks have almost identical spectral indices, consistent with the emission being optically thick.

ALMA 12 CO Imaging
To construct channel maps of CO emission, a CLEAN procedure was performed on the continuum-subtracted data. We used a Briggs weighting scheme with a robust parameter of 1.0. This weighting yields the best results in terms of achieving good signal-to-noise without compromising on resolution. Channel maps were produced with a spectral resolution equivalent to 1 km s −1 . These broad channels are needed to pinpoint fast Keplerian kinematics as opposed to slow outflows. Each channel map has an rms noise of 1.2 mJy beam −1 , for a CLEAN beam of 0 1×0 09. The 1σ noise level is 2 mJy beam −1 if systematics (large scale fringes in the central channels) are also included. As shown in Hales et al. (2015), the kinematics is complex and it is heavily influenced by large scale cloud emission and absorption, and also possibly a slow outflow. A full analysis of the global kinematics, i.e., how outflowing and cloud material connect to the binary FUOri kinematics, will be presented in a future publication (S. Hales et al. 2020, in preparation). Here we focus on probing the kinematics of the gas in the vicinity of the FUOri stars.

Radiative Transfer (RT)
The RT modeling procedure described in Cieza et al. (2018) and Hales et al. (2018) is used to derive disk structural parameters. The model consists of a passively heated disk characterized by a power-law surface density profile with slope γ and a characteristic radius R c , defined by: The scale height as a function of radius is given by: where h c is the scale height at the characteristic radius R c , and Ψ defines the degree of flaring in the disk. Here we use to enable comparison with previous modeling of FUor sources (e.g., Cieza et al. 2018). Equation (1) can be integrated to calculate the disk mass as This simple disk can therefore be described by five free parameters M d , γ, R c , H 100 , and Ψ (flaring). The flux emerging from the parametric disk model is computed using the RT code RADMC-3D (Dullemond et al. 2012). Since FUOri is a variable embedded object, the information about its stellar spectrum is not well known. Here we assume an effective stellar temperature of 10 4 K and a stellar radius of 5R e to account for the stellar photosphere and the accretion luminosity. The computational grid extends from the radius at which the dust temperature is higher than the sublimation temperature of 1200 K, up to an outer radius of 100 au. Since the model is axisymmetric, the computation is only done in the radius and colatitude. The latter extends from 0 (pole) to π/2 (midplane) radians. Radial and colatitude domains are sampled with a grid of 256 by 64 cells, respectively.
We adopt a distribution of dust grains with a power law in size a, given by n(a) ∝ a −3.5 , extending from 0.1μm to 3 mm. For the dust optical properties, we use a mix of amorphous carbon grains and astrosilicate grains (see Cieza et al. 2018 for details). The absorption opacity at 1.3 mm is thus κ abs = 2.2 cm 2 g −1 . The temperature of the dust particles is calculated using RADMC3Dʼs mctherm module. Since we aim to explore thousands of different models, the modified random walk option is enabled to speed up calculation over optically thick regions.
The model parameters {M dust , γ, R c , H 100 , Ψ} were constrained using a Bayesian approach. In addition to the disk structure parameters, a centroid shift (δx, dy) is also optimized. The inclination angle i and PA of the model are fixed to the values obtained from imfit (see Section 2.1). This is a compromise we think necessary as the emission subtends only two to three resolution beams. As will be shown in Section 4, this disk orientation is consistent with signatures of Keplerian rotation seen in 12 CO channel maps. A 1.3 mm image is produced via ray tracing, using RADMC3Dʼs second-order integration of the RT equation. The model visibilities are interpolated to the same uv points as the observations via a fast Fourier transform (we use the same algorithm used and described in Cieza et al. 2018 andMarino et al. 2017). These model visibilities are compared against the observed visibilities of each component in FUOri.
The posterior distributions of each parameter are sampled with the EMCEE Markov Chain Monte Carlo (MCMC) algorithm (Foreman-Mackey et al. 2013). This allows us to determine the set of parameters which maximizes the likelihood function (via a χ 2 comparison computed on the visibility plane). Since the disk model neglects viscous heating, the derived temperatures serve only as a crude approximation. Nevertheless, it allows for an estimate of the size and bulk mass of the dust disk to be obtained.

Modeling Results
The results of the MCMC RT modeling, after running 2000 iterations (∼10 times the autocorrelation time) with 240 walkers, are shown in Figure 2. The posterior distributions of the disk structural parameters are single peaked and relatively narrow, for both components. The 1.3 mm observations of FUOri north and FUOri south can be described by disk profiles, with characteristic radii of approximately 11 au for both disks. The slope of the surface density distribution is ∼0.5, also similar for both disks. The northern component requires a slightly hotter disk with a scale height H 100 ≈2.8, compared to the southern dust emission which requires a scale height of ∼2.3 au at 100 au. See the corner plots in Figure 2 for parameter uncertainties.
The masses of the dust disks derived from our best RT models are relatively small, with 22±2M ⊕ for FU Ori north and 8.8±1.4M ⊕ for FU Ori south. Our RT calculations show that the disks around the northern and southern components become optically thick at  r 11 au and r6 au, respectively, therefore the inferred total dust masses are likely underestimated. This can be explained by the temperatures in our RT calculation, which reach higher values than the 20K assumed in previous works. If the bulk of the mass comes from a region at ∼50-80K, the mass estimates agree within uncertainties. Figure 3 shows the 12 CO channel maps spanning 11 km s −1 in velocity, centered near FUOri north's systemic speed of 11.4 km s −1 . 13 The closest channel to the systemic velocity (12 km s −1 ) suffers from foreground absorption. The location of each disk is indicated with white (dashed) crosses, with major axes parallel to the disk PA inferred from the continuum data.

Evidence of Disk Rotation around Each Component
The local kinematics around the location of FU Ori north show the signature of a rotating disk. The emission changes from blueshifted (velocities 11 km s −1 , in a southeasterly direction) to redshifted (13 km s −1 , toward the northwest). The red-and blueshifted emission loci are symmetrical with respect to the peak in dust continuum, along a disk PA of 135 • east of north, i.e., the disk orientation derived from the continuum emission.
Evidence of spatially resolved disk rotation around FUOri south is also present in the channel maps, although less clear than the emission patterns in the vicinity of the northern component. There is a switch from blue-to redshifted emission along the southern continuum disk PA. The switch happens at a systemic speed (for the southern disk) between the 9 and 10 km s −1 channels. The southern disk emission is much fainter than its northern counterpart in channels with velocities >11 km s −1 . This could be due to intracloud absorption around the southern component. The larger scale kinematics are complex and heavily influenced by cloud emission and absorption, and also by a slow outflow (Hales et al. 2019, in preparation).
CO channels show a prominent elongated emission feature between 9 and 14 km s −1 , immediately east of the primary. At 9 km s −1 , the extended emission is connected (via the lower level flux countour) to the southern component's kinematics. Yet at 10 and 11 km s −1 the emission appears connected to the northern component's blueshifted disk kinematics. Interestingly, the brightest feature seen in the reflected light also appears as an arc-like elongation to the east of FU Ori north (Takami et al. 2018). As shown in Figure 3, the extended arc in CO emission could indeed be a counterpart to the bright arc in the scattered light. This would imply that the arc-like structure could be moving at (projected) speeds of ∼2-3 km s −1 , with respect to the north component's systemic speed.
Although using wide spectral channels allows us to pinpoint disk kinematics, this resolution is not ideal to determine the kinematics of the bright blue arc. To study the nature of this feature, whether it is an outflowing or inflowing structure, requires deep gas observations in optically thinner tracers.
The presence of a blueshifted large scale outflow to the northeast of FU Ori, as well as a redshifted counterpart to the southwest, will be presented in a future work (Hales et al. 2019, in preparation). Assuming that the outflow is launched from the FUOri north disk implies that the near side of the disk is to the southwest. The evidence for rotation seen in the 12 CO channel maps allows us to infer a counterclockwise disk rotation in the plane of the sky. Figure 4 shows a cartoon representation of the inferred geometries of the FUOri binary. The inclinations and PAs correspond to those calculated from the continuum image. The kinematics of each disk, inferred from the 12 CO channel maps, are shown as a red/blue gradient. In the case of FUOri south, the disk's near/far side and the sense of rotation are unknown.

FU Ori Disk Sizes and Masses
Modeling of the continuum emission yields disk sizes with a characteristic radius of ∼11 au for both components. These small disk sizes are comparable to those of protostellar disks, which have typical radii of less than 10 au at 8 mm  FUOri disk is ∼0.4, closer to T Tauri disks (Andrews et al. 2010) than other FUor sources observed at 1.3 mm ). The dust masses and disk sizes inferred from the RT modeling are small compared to the continuum emission from other FUor and EX Lupi objects (EXor) sources. Cieza et al. (2018) fit the continuum emission of eight FUor and EXor targets (their Figure 6), and found that FUor are more massive than EXor. Interestingly, our resolved observations of FUOri north and south suggest that the FUor prototype system has dust masses comparable to the EXor sources in the sample in Cieza et al. (2018).

Update on the FU Ori Accretion Disk Model
With the assumption that the inclination of FUOri north is 55 • , previous SED fitting to its optical-IR spectrum suggests that FUOri north harbors a disk accreting at a rate of  »´-M 2.4 10 4 M e yr −1 around a 0.3 M e star (Zhu et al. 2007). The inner radius is 5 R e in this accretion disk model.
Our observations presented here suggest an inclination angle of ∼35 • . This new inclination allows us to update the disk and star parameters for FUOri north. First, to produce the same line width that explains the optical spectrum of FU Ori ), the stellar mass has to increase to 0.6M e . With the smaller inclination angle, the disk luminosity is only 70% of the previous estimate, which leads to 84% of the previous disk inner radius and 59% of the previous  M M value. Thus, the disk inner radius is 3.5 R e and  M is 3.8×10 −5 M e yr −1 with the new central star mass.
Given the very low mass of the disks and the high accretion rate inferred above, the accretion event must last for a short time compared to the disk lifetime. The high accretion rate suggests that the mechanism behind the outbursts in luminosity happens in an episodic way. There is the possibility that the disk mass is replenished efficiently by inflow of material from the cloud or by cloudlet capture (e.g., Dullemond et al. 2019). Our observations do not rule out the presence of inflowing material and we hope to further explore this scenario in follow-up observations connecting Figure 3. Complex kinematic environment and evidence of disk rotation seen in the 12 CO channel maps of the FUOri system. The orientation (PA and inclination) inferred from the continuum maps are indicated by the white dashed lines for each disk. The northern component shows signatures of rotation as blue-and redshifted 12 CO emission loci along the disk PA. The red and blue points show the stellar locations of FU Ori north and south, respectively. The velocities, in km s −1 , are given in the lower right corners. Each panel shows 12 CO emission via color-filled contour maps for 13,17,21,25,29,33,37, and 41 mJy beam −1 . The thicker unfilled white contour shows the 5σ level, i.e., 10 mJy beam −1 (σ is 2 mJy beam −1 ). The beam size is shown as a gray ellipse in the lower left corner. The systemic velocity (∼11.7 km s −1 ) channel at 12 km s −1 suffers heavy foreground absorption and displays no emission above 5σ. Relative R.A. and decl. are shown with respect to the location of FUOri north (peak of the continuum). The background image shows the High-Contrast Coronographic Imager for Adaptive Optics (HiCIAO) polarized scattered light map published in Takami et al. (2018), see also Liu et al. (2016).
the large scale outflows with higher sensitivity and resolution local kinematics. In the following subsection we discuss the possibility of interaction scenarios.

Dynamically Perturbed Kinematics?
The disk rotation patterns around both components are skewed, in the sense that a unique PA is not sufficient to represent the emission loci over the full velocity range. In the FUOri binary system, this could suggest that the disks are being perturbed by mutual interactions. An internal perturbation, e.g., driven by self-gravity, is unlikely given the low mass of the disks. Given the proximity of both stars, we think that gravitational interactions are likely to be at play. Although the dust disks share similar geometries, the orbit of the perturber, FU Ori south, need not be in the same plane. The bright arc, seen in scattered light and CO in Figure 3, may be explained as a dynamical feature out of the plane of the FU Ori north disk if created by flyby or disk-disk encounter (see Figure 4 in Cuello et al. 2020).
The compact size of the dust disks can also be explained within a flyby scenario. An (inclined) prograde encounter, besides stripping the disk of dust material in the outer regions, leads to an enhancement of the density in the inner regions of the disk (see the red curves in Figure 11 in Cuello et al. 2019). This prograde encounter could also potentially explain the outburst event via enhanced accretion. However, in this case the perturbation only happens once.
An alternative scenario to a disk-disk interaction has been proposed by Dullemond et al. (2019). Here, the capture of a cloudlet or cloud fragment also leads to arc-shaped reflection nebulae. The capture of this cloud fragment also replenishes the disk allowing for a fresh supply of material to maintain the high accretion rate.

Concluding Remarks
The new 1.3 mm ALMA data presented here resolve the continuum and gas emission around both FUOri components. We have derived sizes and orientations for each disk from continuum emission. The disk orientations suggest that FU Ori north and south share similar inclination angles. If we were to assume the disks are in a near-coplanar configuration, they are separated by a deprojected distance of 250 au. Their sizes are remarkably similar taking into account that the southern stellar component is at least twice as massive as the northern star (Beck & Aspin 2012).
The spatially resolved 12 CO kinematics allowed for disk rotation to be identified in the vicinity of each component. The emission revealing disk rotation also appears asymmetric and skewed, suggesting the disks are subject to interaction in the form of a flyby. This interaction could be due to a mutual diskdisk encounter and/or cloudlet capture within the cloud. Moreover, the slow channels (near systemic velocity) reveal CO emission which is spatially coincident with the bright reflection arc previously reported in scattered light observations. The elongated feature is, to some extent, connected to both the north and south components. To understand whether this emission corresponds to inflow or outflow to or from any of the FU Ori components will require deeper CO observations to sample the kinematics with finer spectral resolution and higher sensitivity. The gas observations also show that the systemic speeds of FU Ori north and south differ by at least 1 km s −1 , which could inform future dynamical models of the binary interaction. Several FUor type objects have known binary companions, therefore investigating the kinematics around these multiple systems can help connect binary interaction and flybys with accretion bursts triggers.
The disks in FU Ori north and south are compact and mostly optically thick, even at mm wavelengths (similarly to Class I, Sheehan & Eisner 2017;Segura-Cox et al. 2018). In order to determine the disk masses and dust properties within these central regions, observations at longer wavelengths with auscale resolutions will be required (i.e., at next generation VLA baselines, Murphy et al. 2018;White et al. 2018).
We thank the anonymous referee for a constructive report. We also thank Sebastián Marino for providing the tool to sample the model in visibilities, Ed Fomalont for useful discussions and suggestions during the analysis of the data, and Michihiro Takami   Science and Technology (MoST) of Taiwan (grant Nos. 108-2112-M-001-002-MY3 and 108-2923-M-001-006-MY3). N.C. acknowledges financial support provided by FONDECYT grant 3170680. We acknowledge support from the Millennium Science Initiative (Chile) through grant RC130007. This work used the Brelka cluster, financed by Fondequip project EQM140101, hosted at DAS/U. de Chile. This paper makes use of the following ALMA data: ADS/JAO.ALMA.2016.1.01228.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.