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

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 688 (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 self-calibration). Another iteration of self-calibration was found to only marginally improve the image quality.

Zoom In Zoom Out Reset image size

As shown in Figure 1, both components of the FU Ori binary are detected in continuum emission, separated by ∼05. FU Ori north has a total flux density of 14.30 ± 1.5 mJy, while the southern component has a flux density of 7.94 ± 0.8 mJy (errors include the 10% uncertainty in the absolute flux calibration). The peak fluxes are 7.7 ± 0.8 and 4.4 ± 0.5 mJy beam⁻¹ for FU Ori north and south, respectively. The emitting regions around each component are resolved at 60 × 42 mas angular resolution. The deconvolved source sizes (as FWHM) measured using the CASA task imfit are 50.4 ± 0.5 × 39.9 ± 0.6 (mas) at PA 1336 ± 17 and 47.7 ± 0.8 × 38.4 ± 0.9 (mas) at 1377 ± 37 for FU Ori north and FU Ori south, respectively. This corresponds to disk sizes of 7.5 au × 5.8 au and 7.2 au × 5.7 au, north and south, respectively. The deconvolved sizes translate to inclination angles of 377 ± 08 and 364 ± 12 for the north and south components, respectively. As will be shown in Section 4, the PA inferred from the continuum is in good agreement with the rotation seen in the kinematics of FU Ori north.

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.

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⁻¹. 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⁻¹, for a CLEAN beam of 01 × 009. The 1σ noise level is 2 mJy beam⁻¹ 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.

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:

$$\begin{eqnarray}&&{\rm{\Sigma }}(R)={{\rm{\Sigma }}}_{c}{\left(\displaystyle \frac{R}{{R}_{c}}\right)}^{-\gamma }\exp \left[-{\left(\displaystyle \frac{R}{{R}_{c}}\right)}^{2-\gamma }\right].\end{eqnarray} \tag{ 1 }$$

The scale height as a function of radius is given by:

$$\begin{eqnarray}&&h(R)={h}_{c}{\left(\displaystyle \frac{R}{{R}_{c}}\right)}^{{\rm{\Psi }}},\end{eqnarray} \tag{ 2 }$$

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 ${H}_{100}\equiv h(100\,\mathrm{au})$ 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

$$\begin{eqnarray}&&{M}_{d}=\displaystyle \frac{2\pi {R}_{c}^{2}{{\rm{\Sigma }}}_{c}}{2-\gamma }.\end{eqnarray} \tag{ 3 }$$

This simple disk can therefore be described by five free parameters M_d, γ, R_c, H₁₀₀, 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⁴ K and a stellar radius of 5 R_⊙ 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² g⁻¹. 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₁₀₀, Ψ} were constrained using a Bayesian approach. In addition to the disk structure parameters, a centroid shift (δx, $\delta y$) 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 ¹²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 and Marino 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 χ² 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.

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₁₀₀ ≈ 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.

Zoom In Zoom Out Reset image size

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\leqslant 11\,\mathrm{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 ¹²CO channel maps spanning 11 km s⁻¹ in velocity, centered near FU Ori north's systemic speed of 11.4 km s⁻¹.¹³ The closest channel to the systemic velocity (12 km s⁻¹) 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.

Zoom In Zoom Out Reset image size

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⁻¹, in a southeasterly direction) to redshifted (≥13 km s⁻¹, 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⁻¹ channels. The southern disk emission is much fainter than its northern counterpart in channels with velocities >11 km s⁻¹. 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⁻¹, immediately east of the primary. At 9 km s⁻¹, the extended emission is connected (via the lower level flux countour) to the southern component's kinematics. Yet at 10 and 11 km s⁻¹ 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⁻¹, 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 ¹²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 ¹²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.

Zoom In Zoom Out Reset image size

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 (Segura-Cox et al. 2018). The slope of the density distribution of each FU Ori disk is ∼0.4, closer to T Tauri disks (Andrews et al. 2010) than other FUor sources observed at 1.3 mm (Cieza et al. 2018). 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).

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 $\dot{M}\approx 2.4\times {10}^{-4}$ M_⊙ yr⁻¹ around a 0.3 M_⊙ star (Zhu et al. 2007). The inner radius is 5 R_⊙ 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 (Zhu et al. 2009), the stellar mass has to increase to 0.6 M_⊙. 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\,\dot{M}$ value. Thus, the disk inner radius is 3.5 R_⊙ and $\dot{M}$ is 3.8 × 10⁻⁵ M_⊙ yr⁻¹ 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 the large scale outflows with higher sensitivity and resolution local kinematics. In the following subsection we discuss the possibility of interaction scenarios.

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.