Sites of Planet Formation in Binary Systems. I. Evidence for Disk−Orbit Alignment in the Close Binary FO Tau

Benjamin M. Tofflemire; Lisa Prato; Adam L. Kraus; Dominique Segura-Cox; G. H. Schaefer; Rachel Akeson; Sean Andrews; Eric L. N. Jensen; Christopher M. Johns-Krull; J. J. Zanazzi; M. Simon

doi:10.3847/1538-3881/ad354d

1. Introduction

Binary stars can profoundly shape the dynamics of the circumstellar environment. In the epoch of planet formation, their orbital motion can alter the distribution and internal kinematics of protoplanetary disks. Stellar binarity is generally expected to impede the formation of planetary systems (Nelson 2000; Zsom et al. 2011; Müller & Kley 2012; Picogna & Marzari 2013; Jordan et al. 2021; Zagaria et al. 2023), but the detailed orbital properties and formation pathways of binary systems that foster planet formation remain largely unknown. Submillimeter observations at high angular resolution, possible with the Atacama Large Millimeter/submillimeter Array (ALMA) observatory, facilitate the investigation of the binary−disk interaction and its effect on planet formation.

The impact binarity has on the planet-forming reservoir is, to first order, dependent on the binary separation. Population studies of (sub)millimeter continuum emission (tracing roughly millimeter-size dust grains) in young binaries find that beyond separations of a few hundred au binary disks are indistinguishable from those orbiting single stars (Akeson et al. 2019). Interior to this value, circumstellar disks become less massive (Jensen et al. 1996; Harris et al. 2012; Barenfeld et al. 2019) and/or smaller (Cox et al. 2017; Zurlo et al. 2020, 2021) as the binary projected separation decreases. Large dust grains can be biased tracers of disk properties owing to processes like radial drift (Weidenschilling 1977; Ansdell et al. 2018), but the trend in disk radii follows the general predictions of binary−disk truncation theory (e.g., Artymowicz & Lubow 1994), in that systems with smaller separations have smaller disks. For completeness, we note that disk material can also reside in circumbinary orientations (e.g., Czekala et al. 2019), but these architectures become less common at larger binary separations (e.g., Akeson & Jensen 2014; Akeson et al. 2019).

Submillimeter observations combining dust and gas provide a more complete description of the binary−disk interaction. For a sample of Taurus binaries, Manara et al. (2019) find that dust disks are smaller and have sharper cutoffs in their radial profiles than single stars. Including gas observations for this sample, Rota et al. (2022) find gas disk radii that are consistent with dynamical truncation theory, assuming modest orbital eccentricities (since the actual eccentricities are not uniformly known). Additionally, they find the average ratio of the gas radius to the dust radius is larger for binaries (∼4) than single stars (∼3; Ansdell et al. 2018; Kurtovic et al. 2021; Sanchis et al. 2021). Models of radial drift in truncated disks have been shown to reproduce this difference (Zagaria et al. 2021a, 2021b) without invoking any additional dynamical effects (e.g., dynamical stirring; Quintana et al. 2007).

Beyond the characterization of disk-bearing binary systems, the absence of disk material is also relevant. The occurrence rate of disks in binary systems is known to decrease with decreasing binary separation (Cieza et al. 2009; Kraus et al. 2012; Cheetham et al. 2015). It is unclear whether this trend is driven by the rapid exhaustion of a truncated disk, a dependence of disk survival on specific orbital parameters (e.g., high eccentricity, mutual disk−orbit inclination), or a product of the binary formation mechanism (e.g., disk fragmentation or core fragmentation, or either with subsequent migration; see Offner et al. 2023).

Inferring the consequence these observational results carry for planet formation is challenging. Empirically, the occurrence rate of transiting planets in field-age binary systems declines as the binary separation decreases (Wang et al. 2014; Kraus et al. 2016; Moe & Kratter 2021). The trend closely follows that of the binary−disk occurrence rate, signaling an intuitive link between disk retention and planet formation. Yet, for the significant minority of binary systems with separations less than 50 au that retain long-lived disks and form planets, it is unclear what properties these systems have that allow them to do so. Emerging results from the binary transiting-planet population suggest a preference for mutual alignment of the planet and binary orbital planes (Behmard et al. 2022; Christian et al. 2022; Dupuy et al. 2022; Lester et al. 2023).

A critical obstacle toward developing a predictive theory of planet formation in binary systems is linking protoplanetary disk properties we can now measure with ALMA to the binary orbital and detailed stellar parameters. The disk studies above have had little more than the projected binary separation to interpret their observations. This limitation is rooted in the difficulty of determining orbital solutions for binaries with separations of a few tens of au. At the distances of the closest star-forming regions (∼140 au), decades of high angular resolution imaging are required to measure full orbits or sufficient arcs. While challenging, these are the prime targets where the binary−disk interaction is pronounced and where binary separations are large enough to allow disk sizes, in principle, that can be resolved with current observatories.

To connect the properties of protoplanetary disks to those of the binary orbit, we have observed 1.3 mm continuum and the ¹²CO J = 2–1 transition using the ALMA observatory for a sample of eight disk-hosting binary systems with known orbital solutions. Observing in extended array configurations, we seek to characterize the binary−disk interaction by resolving the location and distribution of circumstellar material. In this study, we present our analysis of FO Tau, a ∼150 mas separation binary (∼20 au). FO Tau has the largest projected separation in our sample, making it a convenient target to demonstrate our methodology. A comprehensive analysis of the full sample will follow. We present the first full orbital solution for the system, combining astrometric and radial velocity (RV) measurements, and forward-model the ALMA visibilities to determine the obliquity of the individual circumstellar disks with the binary orbital plane. This work provides one of the first analyses linking stellar and protoplanetary dynamics, marking a critical step toward understanding planet formation in the binary environment.

2. Observations and Data Reduction

2.1. ALMA

2.1.1. Observations

Our ALMA observations of FO Tau took place in Cycle 7 (Project Code 2019.1.01739.S) using Band 6 receivers in dual-polarization mode. Three spectral windows were placed to sample continuum at 231.6, 234.0, and 245.9 GHz, each with a bandwidth of 1.875 GHz sampled at 31.25 MHz resolution. A fourth spectral window was centered on the ¹²CO J = 2–1 transition (230.538 GHz) at the source's heliocentric velocity (∼16 km s⁻¹; Kounkel et al. 2019). This CO window covered 938 MHz with a two-channel velocity resolution of 488 kHz (0.64 km s⁻¹). This setting was chosen to yield a velocity resolution required to resolve a Keplerian rotation profile while maintaining a broad bandwidth for continuum sensitivity.

To achieve high angular resolution while maintaining sensitivity to extended emission, observations were made in compact and extended array configurations. The compact configuration had baselines between 14 m and 3.6 km, corresponding to an angular resolution and maximum recoverable scale of 0 farcs 1 and 2 farcs 1, respectively, in Band 6 (roughly the C6 configuration). The extended configuration baselines spanned 40 m to 11.6 km, reaching an angular resolution of 0 farcs 034 and maximum recoverable scale of 0 farcs 83 (roughly the C8 configuration). Observations in the compact and extended configurations took place on 2021 July 18 11:17 and 2021 August 22 09:15 with 6- and 50-minute on-source times, respectively, following the standard calibrator observing sequence. The angular scales these observations correspond to a spatial resolution of ∼5 au and a maximum recoverable extent of emission of ∼280 au at FO Tau's distance (d = 135 ± 4 pc). This sensitivity range covers the physical scales of interest for our analysis, specifically probing both circumstellar and circumbinary emission.

2.1.2. Calibration and Imaging

Extended and compact configurations were calibrated separately by the standard ALMA pipeline (CASA v6.1.1.15; McMullin et al. 2007) and provided as a fully calibrated measurement set by the North American ALMA Science Center. The calibration includes data flagging, phase corrections from water vapor radiometer measurements, bandpass and amplitude calibration, and a temporal gain correction. From the pipeline-calibrated data we followed the post-processing procedures developed for the DSHARP program, which focus on combining array configurations and self-calibration. A full description of this approach with accompanying scripts can be found in Andrews et al. (2018) and on the program's data release web page.¹² We provide an outline of our procedure below using CASA v6.5.1.23.

First, we created a continuum data set for each measurement set (compact, extended), removing channels within ±20 km s⁻¹ of the ¹²CO J = 2–1 transition at the FO Tau systemic velocity. We imaged each data set using tclean (Briggs weighting, robust = 0) to test the astrometric alignment and flux scaling. A shift was made to align the phase center of each measurement set on the FO Tau A continuum source. The flux scaling between the data sets was consistent within the uncertainty (flux ratio of 0.999 ± 0.006).

Next, we began an iterative phase-only self-calibration process on the compact configuration. We imaged the compact configuration with a broad cleaning mask that encompassed both stellar components, which created a source model we used to calibrate the visibilities against. An annulus outside this mask was used to measure the signal-to-noise ratio (S/N). We repeated this process, reducing the solution interval until one of the following conditions was met: the S/N does not improve from the previous iteration, the number of flagged solutions for any time interval exceeds 20%, or we reach the native integration step (6 s). Solution intervals of the total duration and of 100 s improved the S/N from 52 to 84 before plateauing at 60 s. One iteration of amplitude-only self-calibration over the full integration was then performed, which improved the S/N to 85.

We then combined the extended and self-calibrated compact measurement sets and repeated this process. Phase-only self-calibration was performed for the full integration, 1500 s, 1000 s, and 600 s, bringing the S/N from to 60 to 150 before leveling out. An iteration of amplitude-only self-calibration did not significantly improve the data quality, leaving the last iteration of phase-only self-calibration as our fiducial continuum data set. Final continuum imaging was performed with interactive masking using Briggs weighting (robust = 0) with a 3σ cleaning threshold. This resulted in a synthesized continuum beam of 45 × 22 mas with a position angle of 22°. The rms deviation of the final image is 0.03 mJy beam⁻¹. This image and its analysis are presented in Section 4.1.

From the initial pipeline-calibrated measurement sets, we also created ¹²CO J = 2–1 data cubes for the compact and extended configurations, including channels within ±20 km s⁻¹ from FO Tau's systemic velocity. The same self-calibration steps above were applied, followed by continuum subtraction. Here we adopted the two-channel spacing of 0.635 km s⁻¹. The resultant combined cube was imaged with interactive masking using Briggs robust = 1 weighting (which is larger than the continuum imaging owing to a lower average S/N in the CO cube). The synthesized CO beam is 79 × 36 mas with a position angle of 25°. The CO channel maps have rms values of 0.9 mJy beam⁻¹ and are cleaned to a depth of 3σ. These maps are presented in Sections 4.2 and 4.3.

2.2. Keck—NIRSPEC behind AO

FO Tau was observed on UT 2010 December 12 (UTC 09:32:18.66) at the Keck II 10 m telescope using NIRSPEC (McLean et al. 2000) behind the adaptive optics (AO) system. The spectrograph slit position angle was set to 230° in order to align both components of the close binary on the slit simultaneously; the 2-pixel slit yielded a resolution of ∼30,000. The AO system reimages the instrumental dimensions onto the detector reduced by a factor of ∼11; thus, the post-AO slit width × length used was 0 farcs 027 × 2 farcs 26. At the time of observations, the air mass was 1.01 and the natural seeing was 0 farcs 5, facilitating robust correction from the AO system. The AO frame rate for FO Tau was 149 Hz, and wave front sensor counts were around 265.

We observed FO Tau in the H band with the N5 filter, with 1.555 μm near the middle of the central order (49). This spectral region is rich in atomic and molecular lines suitable for the measurement of fundamental atmospheric and other properties in low-mass stars spanning effective temperatures of 3000–5000 K. Keck's cold, high-altitude, and low-humidity site on Maunakea ensures that no significant telluric absorption is present in this near-infrared (NIR) spectral range; thus, the data shown in Figure 1 have not been divided by an A0 star.

**Figure 1.** Spatially resolved spectra of FO Tau A (top) and B (bottom). The continuum-normalized spectra are presented in colored lines, and the best-fit model is overlaid in black.
Download figure:
Standard image High-resolution image

Internal comparison lamp lines were used to obtain a second-order dispersion solution; median-filtered stacks of 10 internal flat lamp and dark frame exposures yielded a master flat and dark. A continuum S/N of ∼130 was achieved with six 300 s exposures of the target binary, taken in nodded pairs along the slit to facilitate background subtraction. Differenced nod pairs were divided by the final dark-subtracted flat. We used the REDSPEC reduction package (Kim et al. 2015) for the array rectification and spectral extraction.

2.3. High Angular Resolution Imaging

We obtained three new AO observations of FO Tau using the NIRC2 camera at the Keck Observatory (Wizinowich et al. 2000) on UT 2016 October 20, 2019 January 20, and 2022 October 19. On each night, we obtained sets of 12 images dithered across the detector in the H_cont and K_cont filters. Each image consisted of 1.0 s exposures with 10 coadded frames. The images were flat-fielded using dark-subtracted, median-combined dome flats. Pairs of dithered images were subtracted to remove the sky background. We observed a point-spread function (PSF) reference star using the same AO frame rate either before or after the observations of FO Tau. On the first two nights we used DN Tau as a PSF reference, while DH Tau was used on the third night. Our three observations were take at air masses of 1.07, 1.12, and 1.02, respectively. In each case the AO system easily separated two sources, producing a clear airy ring for each component.

We fit a binary model to the images of FO Tau using the PSF fitting technique described by Schaefer et al. (2014). The procedure uses the PSF star to construct a binary model and searches through a grid of separations and flux ratios to solve for the best-fitting binary separation, position angle, and flux ratio. On UT 2016 October 20, the shape of the AO-corrected images varied over time because of an oscillation problem with the Keck II telescope (C. Alvarez 2024, private communication). To account for the changing PSF shape between individual images, we created an effective PSF for each image using the two components in FO Tau following the approach described by Schaefer et al. (2018), instead of using the observed PSF reference star to model the binary. We corrected the binary positions using the geometric distortion solutions published by Service et al. (2016), applied a plate scale of 9.971 ± 0.004 mas pixel⁻¹, and subtracted 0 fdg 262 ± 0 fdg 020 from the measured position angles. The positions are averaged over the measurements from individual images, and uncertainties are computed from the standard deviation.

We also obtained AO observations of FO Tau with NIRC2 in the J, H_cont, and K_cont filters on UT 2009 November 21. We obtained two images in a single dither position, each consisting of 40 coadded exposures of 0.5 s each. The seeing was better than average (∼0 farcs 4), and the system was observed at moderate air mass (∼1.55), yielding diffraction-limited PSFs even at J band. The images were processed and analyzed with PSF fitting photometry following the methods described by Kraus et al. (2016). To briefly recap, the relative astrometry and photometry of the FO Tau binary were fit within each image by testing the χ² goodness of fit for a library of potential empirical templates (encompassing the 1000 images of single stars that were taken closest to that date) using an initial estimate of the photometry/astrometry, and then the best empirical template was used to optimize the fit of the photometry and astrometry, and we iterated between the two stages until the same empirical template produced the same best-fit values. The values for the two images were then averaged, with the nominal two-value rms being adopted as the uncertainty because it matches the typical rms (∼0.005 mag) for contrast measurements of other bright, similar-contrast binaries separated by ∼3 λ/D.

Table 1 presents the astrometric measurements from the literature and this work, which are included in our fit to the FO Tau binary orbit. The flux ratios of FO Tau B relative to A measured from our Keck AO observations are given in Table 2. The AO observations are available for further inspection on the Keck Archive.¹³

Table 1. High Angular Resolution Orbit Monitoring

Date	ρ (arcsec)	P.A. (deg)	Provenance	References
1991.7139	0.165 ± 0.005	180.0 ± 4.0	Calar Alto 3.5 m-Speckle	(1)
1991.7933	0.161 ± 0.001	181.7 ± 0.6	Palomar Hale 5 m-Speckle	(2)
1993.7591	0.183 ± 0.005	182.0 ± 0.9	Calar Alto 3.5 m-Speckle	(3)
1994.7995	0.153 ± 0.002	190.6 ± 0.4	Palomar Hale 5 m-Speckle	(2)
1994.9500	0.159 ± 0.004	189.7 ± 0.9	Calar Alto 3.5 m-Speckle	(3)
1994.9637	0.154 ± 0.002	191.2 ± 0.4	Palomar Hale 5 m-Speckle	(2)
1995.7687	0.156 ± 0.004	193.1 ± 1.3	Calar Alto 3.5 m-Speckle	(3)
1996.7379	0.152 ± 0.004	194.3 ± 0.9	Calar Alto 3.5 m-Speckle	(3)
1996.9103	0.143 ± 0.004	200.0 ± 1.1	Calar Alto 3.5 m-Speckle	(3)
1996.9268	0.145 ± 0.007	193.7 ± 1.0	IRTF-Speckle	(4)
1996.9295	0.160 ± 0.001	188.4 ± 5.9	IRTF-Speckle	(4)
1997.1759	0.153 ± 0.003	194.7 ± 0.4	HST-WFPC2	(4)
1997.8741	0.150 ± 0.005	198.9 ± 0.8	Calar Alto 3.5 m-Speckle	(3)
1999.79	0.146 ± 0.001	203.8 ± 4.1	HST-STIS	(5)
2001.107	0.145 ± 0.004	207.6 ± 0.6	Calar Alto 3.5 m-Speckle	(6)
2001.838	0.149 ± 0.004	206.7 ± 1.3	Calar Alto 3.5 m-Speckle	(6)
2009.888	0.138 ± 0.002	234.9 ± 0.6	Keck-NIRC2 AO	(7)
2011.7780	0.1372 ± 0.0002	241.1 ± 0.1	Keck-NIRC2 AO	(8)
2013.0729	0.134 ± 0.002	245.9 ± 0.8	Keck-NIRC2 AO	(8)
2015.7151	0.124 ± 0.012	253.2 ± 3.7	ALMA	(9)
2016.8024	0.1375 ± 0.0018	258.1 ± 0.7	Keck-NIRC2 AO	This Study
2019.0523	0.1369 ± 0.0022	266.0 ± 0.9	Keck-NIRC2 AO	This Study
2021.5437	0.1370 ± 0.0005	274.7 ± 0.2	ALMA	This Study
2022.7982	0.1359 ± 0.0010	278.3 ± 0.2	Keck-NIRC2 AO	This Study

References. (1) Leinert et al. 1993; (2) Ghez et al. 1995; (3) Woitas et al. 2001; (4) White & Ghez 2001; (5) Hartigan & Kenyon 2003; (6) Tamazian et al. 2002; (7) Keck Archive D N2.20091121.2475 (PI: Hillenbrand); (8) Schaefer et al. 2014; (9) Akeson et al. 2019.

Download table as: ASCII Typeset image

Table 2. Flux Ratios for FO Tau from Keck Adaptive Optics Observations

Date		Flux Ratio (f_B/f_A)
	J	H_cont	K_cont
2009.888	0.932 ± 0.005	0.885 ± 0.003	0.821 ± 0.003
2016.8024	⋯	0.880 ± 0.005	0.810 ± 0.004
2019.0523	⋯	0.862 ± 0.054	0.821 ± 0.010
2022.7982	⋯	0.891 ± 0.006	0.793 ± 0.006

Download table as: ASCII Typeset image

2.4. Time-series Photometry

FO Tau was observed by the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) in Sectors 19, 43, and 44 with a 2-minute cadence. The photometry was reduced with the SPOC pipeline (Jenkins 2015; Jenkins et al. 2016). Our analysis makes use of the simple aperture photometry (SAP) and the presearch data conditioning SAP (PDCSAP; Smith et al. 2012; Stumpe et al. 2012, 2014) light curves. All of the TESS data used in this paper can be found in MAST:10.17909/zpka-v291.

2.5. Literature Photometry and Astrometry

Our analysis makes use of photometric measurements from the literature and astrometry from Gaia DR3. The relevant measurements are included in Table 3.

Table 3. Known and Derived Properties of FO Tau

Parameter	FO Tau A	FO Tau B	References
Identifiers

Gaia DR3	163183644576299264		Gaia DR3
2MASS	04144928+2812305		2MASS
TIC	56627416		Stassun et al. (2018)

Astrometry, Distance, and Photometry

α R.A. (J2000)	04:14:49.3		Gaia DR3
δ Decl. (J2000)	+28:12:30.5		Gaia DR3
μ_α (mas yr⁻¹)	6.8 ± 0.3		Gaia DR3
μ_δ (mas yr⁻¹)	−23.9 ± 0.2		Gaia DR3
ϖ (mas)	7.3 ± 0.2		Gaia DR3
RUWE	11.499		Gaia DR3
Distance (pc)	${135}_{-3}^{+4}$		Bailer-Jones et al. (2021)
F_phot at 6100 Å (10⁻¹⁴ erg s⁻¹ cm⁻² Å⁻¹)	0.66 ± 0.11	0.62 ± 0.11	Hartigan & Kenyon (2003)
F_phot at 7510 Å (10⁻¹⁴ erg s⁻¹ cm⁻² Å⁻¹)	6.9 ± 1.2		Herczeg & Hillenbrand (2014)
F_phot in H-band (10⁻¹⁴ erg s⁻¹ cm⁻² Å⁻¹)	2.3 ± 0.3	2.2 ± 0.3	This work
ΔJ (mag)	0.077 ± 0.014		This work

Stellar Properties

T_eff (K)^a	3475 ± 50(±100)	3450 ± 50(±100)	This work
R_⋆ (R_⊙)^b	${1.48}_{-0.09}^{+0.10}\,(\pm 0.19)$	1.42 ± 0.09(±0.18)	This work
L_⋆ (L_⊙)^b	0.25 ± 0.02(±0.05)	0.24 ± 0.02(±0.05)	This work
M_⋆ (M_⊙)	${0.35}_{-0.05}^{+0.06}$	0.34 ± 0.05	This work
log(g)	3.9 ± 0.1	3.8 ± 0.1	This work
B (kG)	1.8 ± 0.1	1.8 ± 0.1	This work
v sin i (km s⁻¹)	14.0 ± 1.0	13.5 ± 0.8	This work
Veiling at 15600 Å	0.35 ± 0.05	0.27 ± 0.05	This work
$\dot{M}$ _⋆ (M_⊙ yr⁻¹)	2.6 × 10⁻⁸	1.7 × 10⁻⁸	Hartigan & Kenyon (2003)
Spectral type (NIR)	M2.5	M2.5	This work
Spectral type (optical)	M3.5	M3.5	Hartigan & Kenyon (2003)

Binary Orbital Properties

T₀ (JD)	${2465000}_{-3000}^{+2000}$		This work
P (yr)	120 ± 20		This work
a (au)	${22}_{-1}^{+2}$		This work
e	${0.21}_{-0.03}^{+0.04}$		This work
i_orbit (deg)	${33}_{-5}^{+4}$		This work
Ω (deg)	${303}_{-10}^{+7}$		This work
ω_A (deg)	${30}_{-20}^{+30}$		This work
Minimum Separation (a(1 − e); au)	18 ± + 1		This work

Disk Properties

1.3 mm F (mJy)	2.96 ± 0.07	2.69 ± 0.07	This work
1.3 mm Peak I (mJy beam⁻¹)	1.31 ± 0.02	1.28 ± 0.02	This work
${\mathrm{log}}_{10}\,({M}_{\mathrm{dust}}/{M}_{\odot })$	≳ − 5.26	≳–5.29	This work
i_disk,dust (deg)	27.3 ± 0.5	26 ± 1	This work
PA_disk,dust (deg)	121 ± 1	121 ± 2	This work
R_{eff,95% dust} (au)	3.7 ± 0.1	3.6 ± 0.5	This work
i_disk,CO (deg)	${27}_{-6}^{+7}$	${40}_{-20}^{+30}$	This work
PA_disk,CO (deg)	${120}_{-20}^{+30}$	${120}_{-50}^{+100}$	This work
R_c,CO (au)	10 ± 2	${8}_{-3}^{+4}$	This work
Ω_disk (deg)	${300}_{-20}^{+30}$	${300}_{-50}^{+100}$	This work
v_sys (km s⁻¹ LSRK)	${8.2}_{-0.5}^{+0.6}$	6 ± 2	This work

Notes.

^aSpectroscopically determined; additional uncertainties from systematic and astrophysical effects are included in the parenthetical (see Section 3.1.1).^bDetermined from SED fitting in Section 3.1.2.

Download table as: ASCII Typeset image

Close binaries (ρ < 1'') are known to complicate the Gaia astrometric pipeline, realized most notably through the renormalized unit weight error (RUWE; Lindegren et al. 2018) parameter. RUWE values above ∼1.2 have been shown to indicate the presence of unresolved companions in main-sequence stars (Rizzuto et al. 2018; Belokurov et al. 2020; Bryson et al. 2020). In disk-bearing stars, this threshold is closer to 2.5 (Fitton et al. 2022). FO Tau's RUWE is 11.499, in line with its known binarity.

The degree to which this elevated uncertainty introduces additional random or systematic errors is unclear. Most relevant for the current analysis is the parallax and the distance derived from it. The Bailer-Jones et al. (2021) geometric distance is ${135}_{-3}^{+4}$ pc, while clustering analyses of Taurus members place FO Tau's subgroup center at ∼130 pc with a radial dispersion of a few parsecs (Esplin & Luhman 2019; Krolikowski et al. 2021). Given that these values are in fair agreement, we adopt the geometric distance in our analysis.

3. Stellar and Binary Orbital Parameters

In this section we describe our measurements and modeling of the stellar and orbital parameters of the FO Tau binary. A summary of our results is provided in Table 3.

3.1. Stellar Characterization

3.1.1. Modeling High-resolution NIR Spectra

The H band is an information-rich spectral region for cool stars, containing diagnostics for the typical stellar parameters (T_eff, v sin i, log(g)), as well as the magnetic field strength (e.g., Han et al. 2023). Figure 1 presents a telluric-free region of the FO Tau spectra. The two spectra are very similar. The most readily apparent difference is the line depths, which are shallower for FO Tau A owing to a high contribution from veiling.

To model these spectra, we develop a grid of H-band spectra using the Spectroscopy Made Easy spectral synthesis code of Valenti & Piskunov (1996) along with the NextGen atmosphere models of Allard & Hauschildt (1995). Our grid covers the range of T_eff, v sin i, magnetic field strength, surface gravity, and veiling expected for young, low-mass stars. The spectra are synthesized using laboratory atomic transition data from the Vienna Atomic Line Database (Piskunov et al. 1995; Ryabchikova et al. 2015) following the procedure of Johns-Krull et al. (1999) and Johns-Krull (2007). The synthetic spectra are then calibrated against the solar spectrum (Livingston & Wallace 1991) and the spectrum of 61 Cyg B (Wallace & Hinkle 1996). Following an iterative procedure, oscillator strengths and van der Waals broadening constants are adjusted in order to obtain the optimum match of the synthesized spectra to the observed standard spectra. Using initial values already spectroscopically estimated, we run our spectra through this grid to obtain estimates for the complete set of desired properties for the young binary component stars in our sample. The best-fitting models are shown as black lines in Figure 1. The best-fit values are compiled in Table 3.

The primary errors quoted in Table 3 are internal fitting uncertainties and do not include systematic or additional astrophysical sources of error. These sources of uncertainty are particularly relevant for the T_eff determination in this work. Systematic T_eff uncertainties can arise from the choice of the model grid. One clear example comes from López-Valdivia et al. (2021), where including the magnetic field as a fit parameter increases the measured T_eff by ∼40 K compared to a nonmagnetic model (for stars with similar magnetic field strengths to FO Tau). A wavelength dependence in the T_eff is also observed, attributed to the chromatic dependence of the spot contribution (e.g., Gully-Santiago et al. 2017). Astrophysically, young stars also have high spot covering fractions (e.g., Fang et al. 2016; Gully-Santiago et al. 2017; Cao & Pinsonneault 2022) that vary with the stellar rotation. In a sample of about a dozen young stars in Taurus, S.-Y. Tang et al. (2024, in preparation; Tang et al. 2023) find changes in T_eff on the order of ∼100 K, based on the rotational phase from H-band line equivalent width ratios. Given these effects, a T_eff uncertainty of 100 K is a more conservative value, which we include as the parenthetical in Table 3.

3.1.2. Modeling the Spectral Energy Distribution

In this section we constrain the stellar masses, specifically the mass ratio, q, and total mass, M_total, to place informed priors on the binary orbit fit in Section 3.3. With only approximately one-quarter of the orbit traced with astrometry, the current observational baseline does not precisely constrain the stellar masses on its own. To establish these priors, we fit model spectra to measurements of the stellar photospheric flux to infer stellar effective temperatures and radii. With these, we derive mass priors from the H-R diagram using stellar evolution models.

Both FO Tau components hosts disks (NIR excess), are actively accreting (variable optical/UV excess), and have significant extinction, which complicates a direct fit of the broadband spectral energy distribution. With the current data set, we can make these corrections for these effects to compute an H-band photospheric flux for each component. Beginning with the unresolved Two Micron All Sky Survey (2MASS) H-band apparent magnitude (Skrutskie et al. 2006), we correct for extinction by computing the H-band extinction (A_H) adopting A_H/A_V = 0.19(± 0.03) (Fitzpatrick 1999 extinction law, updated in the NIR by Indebetouw et al. 2005) and the average visual extinction measured by Hartigan et al. (1995) and Herczeg & Hillenbrand (2014) from flux-calibrated, low-resolution optical spectra (A_V = 1.98 ± 0.20). The average (error-weighted) flux ratio from the spatially resolved NIRC2 H-band photometry (Table 2) is then used to split the corresponding flux between the two components, where we have propagated the flux ratio rms as the uncertainty. Finally, the veiling measurements from our spatially resolved, H-band spectra (Section 3.1.1; Table 3) are used to remove the contribution of each disk. A 30% error on the veiling is propagated following the average level of H-band veiling variability found in Sousa et al. (2023). The resulting photospheric H-band flux is provided in Table 3.

To supplement these values, we also adopt literature measurements of photospheric fluxes that have been corrected for extinction and accretion, namely, spatially resolved fluxes at 6100 Å from Hubble Space Telescope–Space Telescope Imaging Spectrograph (HST-STIS) spectra (Hartigan & Kenyon 2003) and a spatially unresolved flux at 7510 Å from ground-based spectra (Herczeg & Hillenbrand 2014). The latter was modeled as a single star, but that is incidental for our use here given the similarity of the primary and secondary (Section 3.1.1). Both studies measure consistent V-band extinctions, 1.9 and 2.05, respectively, but given the separate analyses, we inflate the uncertainties on these measurements, corresponding to ±0.2 mag of extinction using the Cardelli et al. (1989) law, which are added in quadrature to the quoted errors. These values and the inflated errors are presented in Table 3.

Additionally, we include a spatially resolved J-band contrast from Keck-NIRC2 (Section 2.3). This observation is not corrected for accretion or disk veiling, but we include it in our analysis because the J band falls at the minimum of each process's contribution. An additional uncertainty based on the H-band contrast variability is added in quadrature. This measurements and its inflated uncertainty are presented in Table 3.

We fit these measurements (two sets of resolved flux values, an unresolved flux, and a magnitude contrast) with synthetic BT-Settl (CFIST) atmospheric models (Allard et al. 2013), interpolating over effective temperature and log(g) values between 3.0 and 4.5. The model has six parameters: a T_eff and radius for each component, a shared log(g), and a distance. Gaussian priors are placed on the effective temperatures, informed by the H-band spectral fit including the inflated uncertainty (Section 3.1.1, Table 3), and on the distance (Table 3). The fit is made within a Markov Chain Monte Carlo (MCMC) framework using emcee (Foreman-Mackey et al. 2013) with 28 walkers. Fit convergence is assessed on the fly by measuring the chain autocorrelation. The fit ends once the chain autocorrelation changes by less than 5% or the number of steps exceeds the autocorrelation time by a factor of 50. The first five autocorrelation times are discarded as burn-in.

Figure 2 presents the result of our best-fit model compared to the photospheric flux measurements. The fit favors lower effective temperatures than those measured in the H band, highlighting, perhaps, some of the model and/or wavelength dependence of the measurement described in Section 3.1.1. We find effective temperatures of ${3370}_{-60}^{+50}$ K and ${3380}_{-60}^{+50}$ K. The uncertainties represent the 68% posterior confidence interval, which only includes formal uncertainties. The fit returns radii of ${1.48}_{-0.09}^{+0.10}$ R_⊙ and 1.42 ± 0.09 R_⊙ for the primary and secondary, respectively. With these posteriors we derive bolometric luminosities of 0.25 ± 0.02 L_⊙ and 0.23 ± 0.02 L_⊙ for the primary and secondary, respectively, using the Stefan–Boltzmann law. From this analysis we adopt the stellar radii and luminosities and present them in Table 3, but we adopt the T_eff value derived in Section 3.1.1.

The uncertainties of this fit are likely underestimated. The spectroscopically informed T_eff prior constrains the posterior, which would favor lower temperatures in its absence. The sources of systematic uncertainties described in Section 3.1.1 are also at play here. A conservative ±100 K uncertainty on the T_eff is more appropriate, placing this value in fair agreement with the spectroscopically determined value. This value is placed in parentheses in Table 3. We estimate the associated systematic uncertainty on the radius and luminosity values by smoothing the T_eff–R and T_eff–L posteriors with a Gaussian kernel. We set the kernel size to achieve a ±100 K T_eff uncertainty and adopt the 68% confidence intervals along the R and L axes in order to capture the covariance between the parameters. These uncertainties are also included in the parentheticals in Table 3.

We map our SED fit T_eff and L distributions in the H-R diagram to stellar mass, interpolating isochrones from the Dartmouth Stellar Evolution Program (DSEP; Dotter et al. 2008) using the scipy griddata function. They correspond to a mass ratio of 1.0 ± 0.1 and a total mass of 0.58 ± 0.05 M_⊙, taking the median and standard deviation of the distributions. We confirm that the output component ages are coeval within observed ranges ( ${\sigma }_{{\rm{\Delta }}| {\mathrm{log}}_{10}\;\tau | }\lt 0.16$ dex; Kraus & Hillenbrand 2009) and find an age of 1.0 ± 0.3 Myr, which is in good agreement with the color–magnitude isochronal age of FO Tau's Taurus kinematic subgroup ( ${1.34}_{-0.19}^{+0.18}$ Myr; Krolikowski et al. 2021).

To estimate model-dependent bias in these measurements, we compare the mass ratio and total-mass values to those derived from different model suites, and we also assess their agreement with dynamical mass measurements from the literature. For the mass ratio, we find consistent values across seven other stellar evolution models (MESA Isochrones and Stellar Tracks (MIST), Choi et al. 2016; Dotter 2016; BHAC 2015, Baraffe et al. 2015; DSEP-magnetic, Feiden & Chaboyer 2012, 2013; Feiden 2016; SPOTS (f_spot = 0, 0.17, 0.34, 0.51), Somers et al. 2020). As such, we adopt the value above as the mass ratio prior in our orbit fit.

For the total mass, we find that values are highly model dependent. The standard stellar evolution models are largely self-consistent (DSEP, MIST, BHAC, SPOTS f_spot = 0), while the models that include prescriptions for the effect of magnetic fields (DSEP-magnetic, SPOTS f_spot > 0) predict masses up to 65% higher. Empirically, model agreement with dynamical masses is mass dependent at young ages. For solar-type stars (FGK spectral types), the DSEP-magnetic models typically perform better than standard models (David et al. 2019; Simon et al. 2019; Braun et al. 2021). At lower masses, mostly as the outcome of small number statistics, it is unclear which model predictions to favor. Although the results are on the border of statistical significance, the studies cited above suggest that the magnetic models overestimate the masses of objects <0.4 M_⊙ based on their H-R diagram locations. Rizzuto et al. (2016, 2020), for instance, find that even the standard DSEP models overpredict dynamical measurements of binary total masses by 25%–85% below 1 M_⊙. This does not suggest that the standard models provide a more physically realistic description of low-mass stars (which are certainly magnetic), only that there appears to be no single set of models that performs best at all masses and ages. In the case of young, low-mass stars, the shortcomings of standard models may simply conspire to provide more accurate mass predictions from the H-R diagram than magnetic models.

Given the variation between models and evidence for systematic offsets with measured dynamical masses, we adopt the total-mass mean from the DSEP model above and a 1σ width of 20% (i.e., ${ \mathcal N }\,(0.58,0.12)$ ) as our prior in the orbit fit. With this prior, we provide support for a scenario in which models with magnetic prescriptions are more accurate and one in which some or all models overpredict the masses of young, low-mass stars.

3.2. Relative Stellar Radial Velocity

With angularly resolved NIRSPEC spectra, we can compute the relative RV of the FO Tau A and B components, which provides a valuable constraint when fitting the orbital solution. We measure the relative RV by computing the spectral-line broadening function (BF; Rucinski 1992) between the primary and secondary spectra. An in-depth description of the BF, as implemented here, can be found in Tofflemire et al. (2019). In short, the BF is the linear inversion between a target and template spectrum and is typically used with an observed spectrum and a narrow-line synthetic template. In that case, the BF represents a reconstruction of the average RV profile of photospheric absorption lines, providing the means to measure the star's RV and v sin i. In the present case, where the primary and secondary spectra are very similar, the BF of these two spectra returns a Gaussian profile centered at the RV offset between the two spectra with a width of the spectral resolution.

We compute the BF for four ∼100 Å regions using the saphires Python package (Tofflemire et al. 2019). The regions are devoid of telluric contamination and together span 15450–16280 Å. A Gaussian profile is fit to each region, and we take the mean and standard deviation of these measurements as our relative RV and its internal RV uncertainty: RV_B − RV_A = − 1.8 ± 0.2 km s⁻¹. The negative value indicates that FO Tau B is moving toward Earth compared to FO Tau A.

3.3. Binary Orbital Solution

We fit the FO Tau orbital solution using the orvara Python package (Brandt et al. 2021). The fit uses a parallel-tempering MCMC framework (Foreman-Mackey et al. 2013; Vousden et al. 2016) with 10 temperatures and 100 walkers, each taking 5 × 10⁵ steps. The orbit model includes nine parameters: the primary mass (M_A), secondary mass (M_B), semimajor axis (a), eccentricity (e), primary star's argument of periastron (ω_A), inclination (i), position angle of the ascending node (Ω), the reference longitude at 2010.0 (λ_ref), and the system parallax (ϖ). The eccentricity and argument of periastron are fit as $\sqrt{e}$ sin ω_A and $\sqrt{e}$ cos ω_A, respectively. For reference, the Ω convention used here is the angle measured east from north where the secondary crosses the plane of the sky toward the observer. Gaussian priors are placed on the parallax (Gaia EDR3) and on the mass ratio (q) and total-mass (M_total) hyperparameters, following the modeling in Section 3.1.2. The chains for each parameter are saved every 50th step, and chain convergence is assessed following the method outlined in Section 3.1.2.

We include the relative RV measurement from Section 3.2 in our fit but inflate its uncertainty to include astrophysical RV jitter. Both stellar components are actively accreting and likely spotted, resulting in surface heterogeneities that drive RV scatter. With only one RV epoch, we cannot measure a jitter value for either star, nor can we fit for it within the orbital solution. As such, we adopt an empirical NIR RV jitter value from the Crockett et al. (2012) RV survey of T Tauri stars. The median NIR jitter of their nine-star sample is ∼400 m s⁻¹, which we add in quadrature to our relative RV error.

We perform an initial orbit fit using the relative astrometry data in Table 1 and the NIRSPEC relative RV. This fit provides posteriors for the stellar masses, which inform our modeling of the ALMA CO visibilities in Section 4.3. The disk model, in turn, fits a center-of-mass velocity for each component, providing an additional relative RV epoch. Our final fit includes the ALMA relative RV, which breaks degeneracies in the radial component of the orbital motion. Figure 3 presents the on-sky orbital solution. Our relative astrometry is shown with blue circles, with the highest-likelihood orbit in black. Dashed and dotted lines mark the line of nodes and argument of periastron, respectively. The Ω symbol marks the ascending node. Colored lines represent random draws from the parameter posteriors, color-coded by their eccentricity. FO Tau B is currently on the far side of its orbit (i.e., farther from Earth than FO Tau A). The derived eccentricity posterior, following the same color scheme, is provided in Figure 4. Figure 5 presents our full measurement set as a function of time compared to the orbit fit. Select parameters of the final fit are included in Table 3. A full table of the fit parameters and their priors (Table 4) and a corner plot (Figure 12) are provided in Appendix A.

**Figure 3.** On-sky projection of the FO Tau orbit. Observations are shown with blue circles. The orbit with the highest likelihood is shown in black. The dotted line signifies periastron passage. The dashed line is the line of nodes, with the Ω symbol signifying the ascending node. Random draws from the posterior distribution are shown in the background, color-coded by the orbital eccentricity. The secondary is farthest from Earth in the southwest portion of the orbit and closest in the northeast.
Download figure:
Standard image High-resolution image

**Figure 4.** FO Tau orbital eccentricity posterior. The probability distribution function is represented with a Gaussian kernel density estimate. Line color matches the color bars in Figures 3 and 5.
Download figure:
Standard image High-resolution image

**Figure 5.** The FO Tau orbit fit. As a function of time, we present the projected separation, position angle, and RV difference from top to bottom, each with its associated residuals. Astrometry measurements are shown with blue circles. RV difference measurements are shown with red circles. The highest-likelihood orbital solution is shown in black. Random posterior draws are shown in the background, color-coded by the orbital eccentricity.
Download figure:
Standard image High-resolution image

We note that near-term improvement in the orbital solution would benefit most from additional relative RV epochs, where the model is less constrained. Spatially resolved, high-resolution spectra from AO facilities are required to make these measurements.

3.4. Photometric Variability

Figure 6 presents spatially unresolved FO Tau light curves from three TESS sectors. The variability is dominated by quasiperiodic brightening events that last multiple days. Robinson et al. (2022) define the system as a "burster," based on the Sector 19 light curve using the Q (periodicity) and M (symmetry) variability metrics developed by Cody et al. (2014). The burster population is characterized by high accretion rates (large UV continuum excesses) and warm inner disks (NIR excess). Variability in the TESS bandpass for these objects does not directly trace the accretion rate as measured at shorter wavelengths (Robinson et al. 2022) and is likely the contribution of accretion- and inner-disk-related processes. Their periodogram analysis produced a peak at ∼4.2 days, in agreement with analysis of TReS light curves in Xiao et al. (2012). We analyze the TESS light curves here in the search for periodic signals that could be plausibly linked to the stellar rotation signal.

We use the Lomb–Scargle periodogram (LSP; Scargle 1982) and autocorrelation function (ACF) to search for periodic behavior in the full three-sector light curve (19, 43, and 44) and in multiple smaller regions, specifically the individual sectors, the union of back-to-back Sectors 43 and 44, and a subset of Sectors 43 and 44 that excludes the two large events at the end of Sector 43. We repeat the analysis for both the SAP and PDCSAP reductions given that additional systematics removed in the latter have been shown to affect rotation period measurements for stars with periods between 6 and 12 days (e.g., returning P_rot/2; Magaudda et al. 2022; Newton et al. 2022).

We do not find a single period across the temporal regions analyzed that has a consistently dominant LSP peak with positive autocorrelation, regardless of the reduction. The ∼4.2-day peak, for instance, is only present in the Sector 19 PDCSAP light curve. The strongest LSP peak is found in Sector 43 at 5.2 days, which is largely driven by the two brightening events in the latter half of the sector (Figure 6, middle panel). Periodogram peaks near ∼5.8 days are present in multiple time ranges (Sectors 19, 44, 43+44, and the 43+44 subset) but are not consistently the primary peak. Analyzing the full light curve, we find the strongest LSP peak at 5.9 days and a single ACF peak at 5.0 days.

The present analysis is suggestive of some astrophysical periodicity between 4 and 6 days; however, we cannot confidently assign it an origin. Assuming that some component of the periodicity is stellar rotation, for the stars' v sin i and radius values, this period range would correspond to inclinations between ∼50° and 75°. Although the formal uncertainty on this derivation is large (e.g., Masuda & Winn 2020), it would disfavor alignment with the binary orbit. Conversely, if we assume that the stars inherit the binary orbital inclination, the rotation period would be ∼2.9 days, which corresponds to negative ACF values (anticorrelation). Given the large amplitude variability and lack of spatial information, it seems most likely that the true stellar rotation periods are unrecoverable from the TESS light curves. As such, we are not able to comment on the alignment of the stellar angular momenta with that of the binary orbit or the protoplanetary disks.

4. Disk Properties

4.1. 1.3 mm Continuum

The 1.3 mm continuum image in Figure 7 displays the clear detection of circumstellar dust associated with each stellar component. The dust emission is very compact, and the observation is likely just on the boundary of resolving the dust disks. We fit the continuum map with a two-dimensional Gaussian profile for each component using the CASA imfit function. The fit provides the peak and integrated fluxes, which are presented in Table 3. The centers of each component fit (α_A = 04:14:49.30089, δ_A = +28:12:29.9844; α_B = 04:14:49.29056, δ_B = +28:12:29.9957) are used to compute a separation and position angle for the fit to the binary orbit (Table 1; Section 3.3) and to fix the disk centers for our CO visibility modeling in Section 4.3. Depending on the initial guess provided to the function, the profiles are either narrowly resolved or consistent with point sources, indicating that we cannot robustly measure disk radii or their projections (e.g., inclination, position angle) in the image plane.

To measure continuum disk properties below the image-plane resolution, we forward-model the continuum visibilities with an exponentially tapered power-law intensity profile,

$\begin{eqnarray}&&I\,(r)={I}_{0}{\,\left(\displaystyle \frac{r}{{R}_{c}}\right)}^{-{\gamma }_{1}\,}{e}^{-{\left(r/{R}_{c}\,\right)}^{{\gamma }_{2}}}.\end{eqnarray} \tag{ 1 }$

Here γ₁ is the power-law index and γ₂ is the exponential-taper index at the transition radius, R_c. This model describes disks with sharply decreasing outer profiles that are seen in wider binary systems (Manara et al. 2019) and allows us to make comparisons with studies using the same model (e.g., Long et al. 2019; Manara et al. 2019). For each disk, this radial profile is modulated by an inclination, position angle, and positional offset resulting in an on-sky model image. Complex visibilities at the observed u-v baselines are computed using the galario package (Tazzari et al. 2018) and finally fit to our observations using emcee. Our fit employs 80 walkers, each taking 7000 steps. The computational expense of this method does not lend itself to the chain autocorrelation convergence approach used in Section 3.1.2, but inspection of the parameter chains shows clear convergence. The last 3500 steps are used for parameter estimation.

The 68% confidence interval of the parameter posteriors only includes systematic uncertainties and produces underestimated uncertainties. To obtain more robust errors, we fit each of the four spectral windows separately. Three are continuum spectral windows with 1.875 GHz of bandwidth, and the fourth is the CO-centered higher-resolution spectral window with 0.9 GHz of bandwidth. Our adopted value is computed as the average of these four fits, weighted by their frequency bandwidth. The uncertainty is the bandwidth-weighted standard deviation.

For brevity, the detailed results of the fits are provided in Appendix B. Table 5 contains the individual fits and the adopted values. Figure 13 presents the visibility and image-plane data−model comparisons and displays the intensity profile for one of the spectral window fits. In short, the dust radii enclosed by 95% of the total flux (R_eff,95%) are 3.7 ± 0.1 pc and 3.6 ± 0.5 au, the inclinations are 27 fdg 3 ± 0 fdg 5 and 26° ± 1°, and positions angles are 121° ± 1° and 121° ± 2° for FO Tau A and B, respectively. These results are included in Table 3.

The integrated 1.3 mm fluxes are consistent with those from Akeson et al. (2019) that were observed with lower angular resolution. Under the assumption that the dust emission is optically thin, we can compute a dust mass with the following equation:

$\begin{eqnarray}&&{M}_{\mathrm{dust}}=\displaystyle \frac{{F}_{\nu }\,{d}^{2}}{{\kappa }_{\nu }\,{B}_{\nu }\,({T}_{\mathrm{dust}})},\end{eqnarray} \tag{ 2 }$

where F_ν is the 230 GHz flux density, κ_ν is the dust opacity (2.3 cm² g⁻¹ at 230 GHz; Andrews et al. 2013) and B_ν(T_dust) is the Planck function evaluated at 230 GHz for the dust temperature, T_dust. The dust temperature is computed assuming ${T}_{\mathrm{dust}}\sim 25{\,({L}_{\star }/{L}_{\odot })}^{1/4}$ K (Andrews et al. 2013 assume optically thin dust). Including an additional 5% uncertainty on the ALMA absolute flux calibration, we compute dust masses of 5.4(± 0.4) × 10⁻⁶ M_⊙ and 5.1(± 0.4) × 10⁻⁶ M_⊙ for the primary and secondary, respectively. These values are taken to be lower limits on the total dust mass considering how compact the disks are and that the inner regions of protoplanetary disks have been shown to have nonnegligible dust optical depths (e.g., Birnstiel et al. 2018; Huang et al. 2018). We do not attempt to compute a total disk mass for this reason.

We also note that at our detection limits we do not detect any extended dust emission that might be indicative of a circumbinary disk. This is consistent with previous observations of FO Tau that were sensitive to larger spatial scales (Akeson et al. 2019).

4.2. ¹²CO J = 2–1

Figure 8 presents maps of ¹²CO J = 2–1 emission. The left panel shows the CO integrated intensity, where extended circumstellar gas disks are observed toward each component. An asinh image stretch is used to emphasize faint emission. Cloud absorption is apparent in the map, affecting the northwest (blueshifted) side of the primary disk and the southeast (redshifted) side and center of the secondary disk. The affected velocity channels correspond to ∼5 and 8 km s⁻¹ (LSRK radio; hatched region in the right panel's velocity color bar). The right panel provides the intensity-weighted velocity map using the bettermoments quadratic method (Teague & Foreman-Mackey 2018). Our observations resolve the blue- and redshifted orbital motion of each disk, but only marginally. The appearance of overlap of the two disks may be the result of the large beam size in comparison to their inherent size, or non-Keplerian flows that arise from tidal interactions. In both maps we employ a conservative Keplerian binary mask, based on the disk modeling in Section 4.3, to suppress background noise for clarity. No emission on the scale of the resolution element or larger is removed by the mask. An examination of CO emission at larger spatial scales did not reveal extended emission beyond the binary orbit that could be indicative of a circumbinary disk.

4.3. Modeling the CO Gas Disk

The spatial and spectral resolution of our CO data allows us to model the physical properties of the FO Tau protoplanetary disks. In the present work, we are primarily interested in the individual disk inclinations and their systemic velocities. At the same time, the sensitivity and angular resolution of our data limit our ability to constrain a full suite of disk parameters (or deviations from a single-star disk model). As such, we fit the continuum-subtracted CO visibilities with an isolated disk model for each component, adopting the stellar parameters above and fixing or marginalizing over the more granular disk parameters.

We adopt the parametric disk model outlined in Flaherty et al. (2015).¹⁴ Given our focus on the bulk disk properties, we provide a brief overview of the model here, but we refer interested readers to the citation above or Appendix C for more details. The model assumes a temperature and gas density structure in hydrostatic equilibrium (Dartois et al. 2003; Rosenfeld et al. 2013) and a Keplerian velocity field with modifications from gas pressure support. With the disk structure defined, line radiative transfer is computed along the line of sight assuming local thermodynamic equilibrium with Hanning smoothing applied to the resultant data cube. Each disk's structure is defined by 12 parameters, four of which are fit (but marginalized over), with the rest adopting fiducial values. The on-sky projection and radiative transfer of the disk are defined by nine parameters, six of which are defined by our ALMA observations or prior measurements, leaving three as free parameters: the disk inclination (i_disk), its systemic velocity (v_sys), and position angle (PA_disk). For reference, the definition of PA in this model is 180° less than the disk's longitude of the ascending node, Ω_disk. A full list of the parameters, their adopted or fit values, and any prior on their fit is included in Table 6 in Appendix C.

The fit of the disk model is made in the u-v plane. We realize an on-sky model of two disks centered at the locations of the continuum centroids and compute complex visibilities at the observed u-v baselines and velocity channels using the galario package. The model visibilities are fit to ALMA observations in an MCMC framework using emcee. The fit is made to velocity channels between −3.0 and 17 km s⁻¹, ignoring channels between 5 and 8 km s⁻¹ that have foreground cloud absorption. In total, 27 channels are used for the fit, each with a velocity width of 0.635 km s⁻¹. The fit is made with 80 walkers, each taking 7000 steps. The expense of the model computation does not allow us to assess convergence following the autocorrelation scheme in Section 3.3. Instead, we assess convergence by visually inspecting the parameter chains and confirming that parameter values are constant in 100 step increments moving backward from the last step. We conservatively use only the last 2000 steps to compute parameter values and 68% confidence intervals. A subset of the fit parameters is provided in Table 3, while the full suite is presented in Table 6 in Appendix C, along with a corner plot.

Channel maps of our CO observations, the best-fit model, and fit residuals can be found in Figure 9. White circles indicate the location of the primary and secondary continuum peaks, and the four channels with white diagonal lines indicate the presence of cloud absorption. As seen in the residuals, the model is able to describe the bulk of the CO emission; however, there is extended emission, particularly in the 10.51 km s⁻¹ channel that is not described by a Keplerian disk model.

In the present fit, we marginalize over the disk radii given the angular resolution of our observations. We place a uniform prior between 1 and 14 au in order to prevent the disks from overlapping along the line of sight. The median and 68% confidence interval for the primary and secondary disk radii are 10 ± 2 and ${8}_{-3}^{+4}$ au, respectively. These radii exceed the predicted dynamical truncation radii for FO Tau's orbital parameters (see Section 5.1.2). In the event that there exist extended CO emission features (tidal or otherwise) that are not captured by our model, the disk radii and other parameters may be affected by our simple disk model attempting to fit extended emission.

To explore this behavior further, we perform two additional MCMC fits that limit the upper bound of the radius priors to 8 and 5.5 au. The disk parameters show some variation, but all have agreement within the 68% confidence interval to the fiducial fit. The radius-constrained fits are unable to describe the emission at larger radii from the continuum peaks, signaling that the large disk radii are not solely the result of the angular resolution of the observation.

5. Discussion

5.1. The Binary−Disk Interaction

In the following subsections we discuss our measurements of the FO Tau orbit and protoplanetary disks in the context of theory describing the binary−disk interaction.

5.1.1. Disk−Orbit Alignment

The primary feature of the binary−disk interaction that the present study allows us to comment on is the alignment between the binary orbit and the individual circumstellar disks. We find evidence for alignment among the binary and protoplanetary orbital planes within the uncertainty of our measurements. The left panel of Figure 10 presents the inclination measurements for the three components. The orbit and CO disk inclinations are presented as the fit posteriors. The dust disk inclinations are presented as vertical bands with widths that correspond to 3σ. Our measurement for the secondary disk's CO inclination has the largest uncertainty, but the peak in the probability distribution aligns with that of the other, more precise measurements.

**Figure 10.** Left: comparison of the projected inclination posteriors of the binary orbit and the individual protoplanetary disks. Measurements from both the continuum (dust) and CO visibility fits are included for each disk. Middle and right: disk−orbit obliquity posteriors for the FO Tau A and FO Tau B CO disks, respectively. In each panel, the dark and faint lines represent the aligned and misaligned solutions, respectively (as labeled). A schematic displaying these scenarios in presented in Figure 11. The gray curve in the right panel corresponds to the obliquity distribution of a randomly oriented disk.
Download figure:
Standard image High-resolution image

Agreement in the projected inclinations, however, does not necessarily correspond to alignment. A more direct measure of binary and disk orbital planes is their obliquity (Θ), the angle between the binary and disk angular momentum vectors. Functionally,

$\begin{eqnarray}\begin{array}{rcl}\cos \ {\rm{\Theta }} & = & \cos \ {i}_{\mathrm{disk}}\ \cos \ {i}_{\mathrm{orbit}}\\ & & +\,\sin \ {i}_{\mathrm{disk}}\ \sin \ {i}_{\mathrm{orbit}}\ \cos \,({{\rm{\Omega }}}_{\mathrm{disk}}-{{\rm{\Omega }}}_{\mathrm{orbit}}).\end{array}\end{eqnarray} \tag{ 3 }$

The three-dimensional constraints on the binary orbit (direction of orbital motion, inclination, and longitude of the ascending node) allow for a direct measure of the binary orbit's angular momentum vector. The disks, however, have a degeneracy in the direction of rotation. Counterclockwise rotation corresponds to an angular momentum vector pointing toward the observer. Clockwise rotation corresponds to an angular momentum vector pointing away the observer. Framed in another way, our observation does not reveal the near and far sides of the disks. This degeneracy can be lifted with future observations of disk gas at high angular resolution that resolve the front and back of the disk atmosphere (e.g., Rosenfeld et al. 2013), or with angularly resolved, scattered-light imaging that is sensitive to forward-scattered light (e.g., Avenhaus et al. 2018).

Until then, two solutions exist: an aligned solution in which the disk rotates counterclockwise with the binary orbit, and a misaligned solution where the disk rotates clockwise, retrograde with the binary orbit. We compute the disk−orbit obliquity for both solutions and present them for FO Tau A and B in the middle and right panels of Figure 10, respectively. The right panel also includes the distributions of obliquities for random disk orientations in gray, for reference. (For clarity, we only include the obliquity distributions for the CO disks.) A schematic aid to visualize these two scenarios is presented in Figure 11. Although we cannot strictly rule out the pathological misaligned scenario, it seems unlikely given the consistency between projected inclinations and the longitudes of the ascending node for the orbit and disks.

**Figure 11.** Schematics of the binary−disk orbital planes for the two possible disk orientations allowed by our observation. The binary orbit is shown in black. The primary and secondary disks are displayed in blue and orange, respectively. Arrows indicate the angular momentum vectors. The dotted line marks periastron passage. The left column is the aligned solution in which the disks are rotating counterclockwise with the binary orbit. The right column is the misaligned solution where the disks rotate clockwise, retrograde with the binary orbit. The top row shows the on-sky projection where the projected disk inclinations are the same, despite having different directions of rotation. The bottom row presents an alternate projection where the difference between the solutions is more readily apparent.
Download figure:
Standard image High-resolution image

Assuming the aligned solution, the immediate question becomes whether FO Tau formed in this way, or whether this quality is the result of subsequent dynamical evolution. If the former, it may be signaling a formation pathway that favors disk retention in close binaries. The fraction of young, close binaries with disks (or mature binaries with planets) might then reflect the frequency of this formation pathway. Coplanarity between the binary and disk planes, a small binary separation, a low eccentricity, and a stellar mass ratio near unity are all expectations from the disk fragmentation paradigm (e.g., Bate & Bonnell 1997; Ochi et al. 2005; Young et al. 2015). At the same time, stellar pairs that form at wider separations, in a core fragmentation scenario, can rapidly migrate and achieve binary−disk alignment via dynamical interactions (Zhao & Li 2013; Lee et al. 2019; Guszejnov et al. 2023). Additionally, disk fragmentation is expected to require massive disks (e.g., Kratter et al. 2010; Zhu et al. 2012), and for the relatively low mass in the FO Tau system, most Class 0/I protostellar analogs are compact and stable to gravitational collapse (Tobin et al. 2020).

There exist many mechanisms in the literature that have been suggested to bring binary−disk systems into or out of alignment (e.g., envelope accretion, ejections, inclination oscillations, damping; see Offner et al. 2023 and references therein). One mechanism we explore, given the small semimajor axis of the system, is the damping of mutual disk−orbit inclinations through viscous warped disk torques (Bate 2000; Lubow & Ogilvie 2000). Following Zanazzi & Lai (2018), for FO Tau's nominal parameters, the disk damping rate takes the form

$\begin{eqnarray}\begin{array}{rcl}{\gamma }_{b} & = & 4.5\times {10}^{-6}\,\left(\displaystyle \frac{\alpha }{{10}^{-2}}\right)\,{\left(\displaystyle \frac{0.1}{{h}_{\mathrm{out}}}\right)}^{2}{\,\left(\displaystyle \frac{{M}_{2}}{0.34\ {M}_{\odot }}\right)}^{2}\\ & & \times \,{\left(\displaystyle \frac{0.35\ {M}_{\odot }}{{M}_{1}}\right)}^{3/2}{\,\left(\displaystyle \frac{22\ \mathrm{au}}{a}\right)}^{6}{\,\left(\displaystyle \frac{{r}_{\mathrm{out}}}{7\ \mathrm{au}}\right)}^{9/2}\,{\mathrm{yr}}^{-1},\end{array}\end{eqnarray} \tag{ 4 }$

where α is the viscosity parameter, h_out is the aspect ratio at the outer disk edge, M₁ and M₂ are the primary and secondary stellar masses, respectively, a is the semimajor axis, and r_out is the outer disk edge. (Here we have assumed the dimensionless viscous coefficient, ${{ \mathcal V }}_{b}\approx 1;$ see Table 1 of Zanazzi & Lai 2018 for justification.) For the representative values above, the alignment timescale $({\gamma }_{b}^{-1})$ is less than the age of the system (≲1 Myr) but would quickly exceed it for smaller viscosity values or larger disk aspect ratios. The latter may be particularly relevant in light of numerical simulations by Picogna & Marzari (2013) that find that disks are dynamically heated by companions, especially in the vertical direction.

The present analysis does not distinguish between formation and subsequent evolution as the source of FO Tau's disk−orbit alignment. A larger sample, particularly at larger separations (where the damping timescales are longer but orbital parameters are admittedly more uncertain), may help to address this question. FO Tau has the largest separation in our full ALMA sample, so our current data set is unlikely to address this ambiguity. It will, however, address whether disk−orbit alignment is common in protoplanetary disk-hosting binaries with small separations.

5.1.2. Disk Truncation

Theory predicts that the binary orbit will drive resonances that truncate the outer edges of circumstellar disks (Artymowicz & Lubow 1994; Lubow et al. 2015; Miranda & Lai 2015). With some dependence on the binary orbital parameters, the disk parameters, and the mutual inclination between the two, the disk truncation radius is on the order of ∼30% of the orbital semimajor axis. Manara et al. (2019) present an equation for the disk truncation radius following the work of Artymowicz & Lubow (1994), which has the form

$\begin{eqnarray}&&{R}_{\mathrm{trun}}=\displaystyle \frac{0.49{{aq}}^{-2/3}}{0.6{q}^{-2/3}+\mathrm{ln}\,(1+{q}^{-1/3})}\left(\,,{{be}}^{c}+0.88{\mu }^{0.01}\right).\end{eqnarray} \tag{ 5 }$

Here a is the semimajor axis, q is the stellar mass ratio (M_B/M_A), e is the orbital eccentricity, and μ is the secondary-to-total-mass ratio (M_B/(M_A + M_B)). The remaining variables, b and c, depend on μ and the disk Reynolds number. Using the range of b and c coefficient values compiled by Manara et al. (2019, their Appendix C.1), we compute disk truncation radii of 5–6 au for FO Tau's orbital parameters.

The dust disk sizes we measure (95% effective radius) are ∼4 au. Radii below the dynamical prediction are consistent with the expectation that dust should be more compact than gas resulting from radial drift (e.g., Ansdell et al. 2018). In the case of close binaries, radial drift may be operating at higher rates given the impact of disk truncation (Zagaria et al. 2021a; Rota et al. 2022).

As for the gas disk, our fit to the CO visibilities only loosely constrains the individual disk radii (Section 4.3) but favors values that are larger than those predicted above (R_c ∼ 8–10 au). While the observations do resolve the CO emission, they do not have the angular resolution to distinguish whether the fit radii are indeed larger than truncation models predict, or whether our simple model is not capturing more complex emission. The latter seems likely given the extended emission described in the next subsection, but determining the degree to which the radii may be inflated will require higher angular resolution observations.

The average ratio of the gas to dust radius (R_CO/R_dust) in wide binaries has been shown to be larger than that for single stars, when using the 95% effective radius: 4.2 for binaries (Rota et al. 2022) and 2.8 for single stars (Sanchis et al. 2021). We do not measure the 95% effective radius for the CO disks given the complicated on-sky projection, but from the critical radius of our CO visibility fit, we compute ratios of gas to dust disk radius of 2.6 ± 0.5 and 2 ± 1 for FO Tau A and B, respectively. These values are below the population averages but are not outliers in the distribution of either binary or single stars.

Lastly, we note that mutual inclinations between the binary and disk planes can lead to larger disk truncation radii (Lubow et al. 2015; Miranda & Lai 2015). However, the effect is relatively small; a 90° mutual inclination would correspond to a ∼30% increase in the disk radius. Given that the best-fit disk radii are still larger than this prediction and our measurements are consistent with alignment, we do not find that the large gas disk radii measured necessarily support this scenario.

5.1.3. Extended CO Emission

In the residual image of the CO visibilities fit (Figure 9, bottom panel set), we find an arc of extended emission to the east of the primary disk (LSRK = 10.51 km s⁻¹). The entire arc is more than 3σ above the background rms noise, with a peak >4σ. The feature is also visible as a spur in the CO intensity and first moment maps (Figure 8).

The arc-like feature is reminiscent of the extended emission that has been seen in gas maps (Rodriguez et al. 2018; Zapata et al. 2020; Cuello et al. 2023) and scattered-light imaging (Weber et al. 2023; Zhang et al. 2023) of wide binary systems at similar ages to FO Tau. These observations can offer a lens into the broader, low-density environment and trace structures linked to formation, which are not probed by the larger dust grains with (sub)millimeter continuum. Some of the systems referenced above may result from flybys (Cuello et al. 2023), rather than more stable binary systems like FO Tau. Still, tidal features like extended spiral arms and bridges connecting circumstellar disks are also expected to arise from the binary−disk interaction (e.g., Nelson 2000; Zsom et al. 2011; Müller & Kley 2012; Picogna & Marzari 2013; Jordan et al. 2021). The present feature, if resolved with future observations, may provide a test of the disk parameters (e.g., viscosity) that determine the size and locations of extended features.

5.1.4. Future Submillimeter Observations

Despite the challenges in observing compact binaries, FO Tau remains one of the best candidates to study the binary−disk interaction. Now with a full orbital solution, it is one of the widest binaries with precisely known orbital parameters and two bright circumstellar disks. Additionally, its relatively low orbital eccentricity, compared to the current sample of young binary orbits (Prato 2023), should allow for larger disk truncation radii compared to other young systems. Future observations at higher angular resolution that target less optically thick gas tracers will fully illuminate the discussion above.

6. Summary and Conclusions

In this study we analyze the orbital, stellar, and protoplanetary disk properties of the young pre-main-sequence binary FO Tau. We have combined decades of high angular resolution imaging, angularly resolved high-resolution NIR spectra, and (sub)millimeter interferometry to develop a comprehensive view of the system's dynamical state. With this we can make some of the first tests of the binary−disk interaction and its effect on planet formation. The main conclusions of our work are as follows:

1.
FO Tau is a young binary system with a semimajor axis of 22 au. The stars are nearly equal in mass, with dynamical masses of 0.35 and 0.34 M_⊙ for the primary and secondary, respectively. The orbital eccentricity (0.21) is low compared to other young binaries at similar separations.
2.
ALMA observations in Band 6 detect continuum and ¹²CO J = 2–1 line emission, resolving the contribution from each binary component. The disks are compact and likely truncated by the binary orbit. Circumstellar dust disks are not resolved in the image plane (∼5 au resolution), and the CO gas disks are only marginally resolved. Foreground cloud absorption contaminates the CO emission from parts of both disks. The integrated 1.3 mm fluxes for each component are nearly equal at 2.96 ± 0.07 mJy and 2.69 ± 0.07 mJy for the primary and secondary, respectively.
3.
Measuring inclinations of the gas and dust disks with fits made in the u-v plane, we find evidence for alignment between both protoplanetary disks and the binary orbit. Our gas disk fit does not provide strong constraints on the disk radii but suggests evidence for either larger radii than predicted by truncation models or a more complex emission structure that is not captured by our model.
4.
We do not find evidence for a circumbinary disk or extended emission detected beyond the binary orbit in dust or CO above an rms of 0.02 and 0.9 mJy beam⁻¹, respectively.

FO Tau is the first young binary system where the properties of the circumstellar protoplanetary disks can be compared to precisely known orbital parameters. The synthesis of these results above points to the FO Tau system as a relatively placid environment (mutually aligned, low eccentricity). Most young binaries at this separation do not host disks, so it is tempting to tie this property to the retention of a substantial reservoir of circumstellar material and, perhaps, a common quality of binaries that form and retain planets. Indeed, there is growing evidence for preferential alignment among circumstellar disks in wider binaries (Jensen et al. 2020, with unknown orbital parameters) and low mutual inclinations in the planet−binary orbital planes of field-age systems (Behmard et al. 2022; Christian et al. 2022; Dupuy et al. 2022; Lester et al. 2023). Although only suggestive now, coupling ALMA interferometry with the growing number of young, disk-bearing binaries with known orbital parameters will facilitate the same detailed dynamical study presented here for a representative population.

Acknowledgments

We would like to thank the referee for a constructive and helpful report. B.M.T. would like to thank Sam Factor, Kendall Sullivan, Neal Evans, and Stella Offner for helpful conversations. B.M.T. is supported by the Heising-Simons Foundation's 51 Pegasi b Postdoctoral Fellowship in Planetary Astronomy.

G.H.S. acknowledges support from NASA Keck PI Data Awards administered by the NASA Exoplanet Science Institute (PI Schaefer; 2016B-N046N2, 2019B-N166). Time at the Keck Observatory was also granted through the NOIRLab (PropID: 2022B-970020; PI: G. Schaefer) supported by the NSF Mid-Scale Innovations Program. L.P. was supported in part by NSF awards AST-1313399 and AST-2109179. Contributions to the development of the synthetic spectral grid were made by Thomas Allen, Ian Avilez, Kyle Lindstrom, Cody Huls, and Shih-Yun Tang, all formerly at Lowell Observatory.

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.1.01739.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), 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.

Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.

This research has made use of the Keck Observatory Archive (KOA), which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science Institute (NExScI), under contract with the National Aeronautics and Space Administration.

The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources that have contributed to the research results reported within this paper.¹⁵

Figures in this manuscript were created using color-impaired-friendly schemes from ColorBrewer 2.0.¹⁶

The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

We would like to acknowledge the Alabama-Coushatta, Caddo, Carrizo/Comecrudo, Coahuiltecan, Comanche, Kickapoo, Lipan Apache, Tonkawa, and Ysleta Del Sur Pueblo, and all of the American Indian and Indigenous Peoples and communities who have been or have become a part of the lands and territories of Texas.

Facilities: ALMA - Atacama Large Millimeter Array, Texas Advanced Computing Center (TACC) - , Keck:II (NIRC2, NIRSPEC).

Software: astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), emcee (Foreman-Mackey et al. 2013), saphires (Tofflemire et al. 2019), scipy (Jones et al. 2001; Virtanen et al. 2020), RotBroadInt (Carvalho & Johns-Krull 2023).

Appendix A: Binary Orbital Parameter Posteriors

This appendix provides a summary of the parameters in the FO Tau orbit modeling and the full results of our fit. Table 4 presents the parameters that are fit, the hyperparameters used to place additional priors on the fit, and the orbital parameters derived from the fit. A corner plot of the fit parameters is included in Figure 12.

**Figure 12.** Corner plot for parameters in the FO Tau orbit fit.
Download figure:
Standard image High-resolution image

Table 4. FO Tau Orbital Parameters

Parameter	Prior	Value	Description
Fit Parameters

M_A (M_⊙)	Uniform	${0.346}_{-0.049}^{+0.055}$	Primary mass
M_B (M_⊙)	Uniform	${0.340}_{-0.049}^{+0.054}$	Secondary mass
a (au)	1/a (log flat)	${22.0}_{-1.3}^{+1.9}$	Semimajor axis
$\sqrt{e}$ sin ω_A	Uniform	${0.19}_{-0.17}^{+0.15}$
$\sqrt{e}$ cos ω_A	Uniform	${0.390}_{-0.082}^{+0.077}$
i_orbit (deg)		${33.0}_{-4.8}^{+4.2}$	Inclination
Ω (deg)	Uniform	${303.4}_{-10}^{+6.7}$	Position angle of ascending nodes
λ_ref (deg)	Uniform	${311.0}_{-8.7}^{+10}$	Reference longitude at UTC 2010.0
ϖ (mas)	${ \mathcal N }\,(7.3,0.2)$	${7.42}_{-0.20}^{+0.20}$	Parallax

Hyperparameters

M_total = M_A + M_B (M_⊙)	${ \mathcal N }\,(0.57,0.11)$	${0.687}_{-0.093}^{+0.10}$	Total mass
q = M_B/M_A	${ \mathcal N }\,(1.00,0.05)$	${0.99}_{-0.11}^{+0.11}$	Mass ratio
sin i_orbit	sini i, i ∈ [0, 180]	${0.544}_{-0.072}^{+0.060}$

Derived Parameters

P (yr)		${124}_{-15}^{+23}$	Period
T₀ (JD)		${2464669}_{-2645}^{+1677}$	Time of periastron passage
e		${0.213}_{-0.032}^{+0.040}$	Eccentricity
ω_A		${33}_{-22}^{+25}$	Argument of periastron for primary orbit

Download table as: ASCII Typeset image

Appendix B: Continuum Disk Modeling

The results of our continuum visibility fitting are presented here. As described in Section 4.1, the four spectral windows are fit separately. The adopted value is the average of the fit values, weighted by their respective bandwidths. Table 5 presents the best-fit model parameters for each spectral window, as well as the adopted values. Figure 13 presents diagnostic plots for our fitting procedure, specifically for continuum SPW 3. It includes a panel or the observed and model visibilities, maps of the data, model and residuals, and the intensity profiles of the best-fit models.

**Figure 13.** Continuum visibility fitting. Results are shown for SPW 3. Left panel: data and model visibilities as a function of *u-v* distance. Middle panels: cleaned images of the data, model, and residuals. Contours in the residual image are set at −3σ and 3σ in dashed and solid lines, respectively. Right panel: best-fit radial intensity profile of the primary (blue) and secondary (orange).
Download figure:
Standard image High-resolution image

Table 5. FO Tau Continuum Disk Parameters

Data Set	F_ν	R_eff,68		R_eff,95		R_c		γ₁	γ₂	i	PA
	(mJy)	(arcsec)	(au)	(arcsec)	(au)	(arcsec)	(au)			(deg)	(deg)
FO Tau A

SPW 1	2.9	0.013	1.7	0.029	3.9	0.029	4.0	0.5	3.2	27.7	121
SPW 2	3.2	0.008	1.1	0.027	3.7	0.034	4.6	0.7	4.6	27.7	120
SPW 3	3.2	0.010	1.3	0.027	3.7	0.032	4.3	0.7	4.2	27.0	122
CO SPW	2.9	0.009	1.3	0.027	3.7	0.031	4.2	0.6	3.2	26.4	121

Adopted value	3.1 ± 0.2	0.010 ± 0.002	1.4 ± 0.3	0.028 ± 0.001	3.7 ± 0.1	0.032 ± 0.002	4.3 ± 0.3	0.6 ± 0.1	3.9 ± 0.7	27.3 ± 0.5	121 ± 1

FO Tau B

SPW 1	2.8	0.003	0.4	0.022	3.0	0.033	4.4	0.9	2.2	25	120
SPW 2	3.1	0.011	1.4	0.030	4.1	0.021	2.8	0.4	1.4	25	121
SPW 3	3.1	0.009	1.2	0.028	3.8	0.029	3.9	0.6	2.1	28	124
CO SPW	2.8	0.007	1.0	0.028	3.7	0.016	2.2	0.5	1.1	27	120

Adopted value	3.0 ± 0.2	0.007 ± 0.003	1.0 ± 0.5	0.027 ± 0.004	3.6 ± 0.5	0.026 ± 0.007	3.5 ± 0.9	0.6 ± 0.2	1.8 ± 0.5	26 ± 1	121 ± 2

Download table as: ASCII Typeset image

Appendix C: CO Disk Modeling and Parameter Posteriors

In this appendix we describe the structure of the disk model used in Section 4.3 and present the complete results of our fit. A summary of the model parameters is provided in Table 6, and a corner plot is presented in Figure 14.

**Figure 14.** Corner plot for parameters in the FO Tau disk fit. Posteriors and two-parameter correlation contours are represented with a Gaussian kernel density estimate. The results for FO Tau A and B (fit simultaneously) are displayed in blue and orange, respectively.
Download figure:
Standard image High-resolution image

Table 6. FO Tau CO Disk Parameters

Parameter	Prior	FO Tau A	FO Tau B	Description	Note
Disk Structure Parameters

q_T	${ \mathcal U }\,(-1,1)$	–0.1 ± 0.2	–0.1 ± 0.2	Temperature radial power-law index	Marginalized over
M_disk (M_⊙)	${ \mathcal U }\,({10}^{-4},0.05$ )	0.02 ± 0.02	0.02 ± 0.02	Disk mass	Marginalized over
R_c (au)	${ \mathcal U }\,(1,14)$	10 ± 2	${8}_{-3}^{+4}$	Critical radius	Marginalized over
T_atm,0 (K)	${ \mathcal U }\,(20,180)$	${80}_{-50}^{+60}$	${90}_{-50}^{+60}$	Atmosphere temperature normalization (r = 150 au)	Marginalized over
γ		1	1	Surface density radial power-law index	Fixed
R_in (au)		1	1	Inner boundary of disk density/temperature calculation	Fixed
R_out (au)		100	100	Outer boundary of disk density/temperature calculation	Fixed
X_CO		10⁻⁴	10⁻⁴	CO gas fraction relative to H₂	Fixed
Z_q,0		33.9	33.9	Vertical temperature normalization (r = 150 au)	Fixed
T_mid,0 (K)		17.5	17.5	Midplane temperature normalization (r = 150 au)	Fixed
M_⋆ (M_⊙)		0.36	0.36	Stellar mass	Fixed
v_turb (km s⁻¹)		0	0	Turbulent velocity	Fixed

Disk Observation Parameters

i_disk (deg)	${ \mathcal U }\,(0,90)$	${27}_{-6}^{+7}$	${40}_{-20}^{+30}$	Disk inclination
PA_disk (deg)	${ \mathcal U }\,(0,360)$	${120}_{-20}^{+30}$	${120}_{-50}^{+100}$	Disk position angle
v_sys (km s⁻¹)	${ \mathcal U }\,(2,12)$	${8.2}_{-0.5}^{+0.6}$	6 ± 2	Systemic velocity
d (pc)		135	135	Distance	Fixed
α_offset (arcsec)		0	−0.137	Spatial offset in R.A.	Fixed
δ_offset (arcsec)		0	0.011	Spatial offset in decl.	Fixed
${v}_{\min }$ (km s⁻¹)		−2.82	−2.82	Velocity of first channel	Fixed
v_step (km s⁻¹)		0.635	0.635	Channel width	Fixed
n_chan		32	32	Number of velocity channels	Fixed

Download table as: ASCII Typeset image

The disk model has a temperature structure following

$\begin{eqnarray}{T}_{\mathrm{gas}}\,(r,z)=\left\{\begin{array}{ll}{T}_{\mathrm{atm}}+({T}_{\mathrm{mid}}-{T}_{\mathrm{atm}})\,{\left(\cos \;\displaystyle \frac{\pi z}{2{Z}_{q}}\right)}^{2}\; & z\lt {Z}_{q}\\ {T}_{\mathrm{atm}} & z\geqslant {Z}_{q},\end{array}\right.\end{eqnarray} \tag{ C1 }$

where the atmospheric temperature (T_atm), midplane temperature (T_mid), and local scale height (Z_q) are a function of the disk radius with the following form:

$\begin{eqnarray}&&{T}_{\mathrm{atm}}\,(r)={T}_{\mathrm{atm},0}{\,\left(\displaystyle \frac{r}{150\ \mathrm{au}}\right)}^{{q}_{T}},\end{eqnarray} \tag{ C2 }$

$\begin{eqnarray}&&{T}_{\mathrm{mid}}\,(r)={T}_{\mathrm{mid},0}{\,\left(\displaystyle \frac{r}{150\ \mathrm{au}}\right)}^{{q}_{T}},\end{eqnarray} \tag{ C3 }$

$\begin{eqnarray}&&{Z}_{q}\,(r)={Z}_{q,0}{\,\left(\displaystyle \frac{r}{150\ \mathrm{au}}\right)}^{1.3}.\end{eqnarray} \tag{ C4 }$

We adopt a simplified power-law gas surface density model of the form

$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{gas}}\,(r\lt {R}_{c})=\displaystyle \frac{{M}_{\mathrm{gas}}\,(2-\gamma )}{2\pi \,({R}_{c}^{2-\gamma }-{R}_{\mathrm{in}}^{2-\gamma })}{r}^{-\gamma },\end{eqnarray} \tag{ C5 }$

where M_gas, R_c, and R_in are the disk gas mass, its outer radius, and its inner calculation boundary, respectively. The disk volume density is computed assuming hydrostatic equilibrium, which in turn sets the deviation from a pure Keplerian rotation profile:

$\begin{eqnarray}&&\displaystyle \frac{{v}_{\phi }^{2}}{r}=\displaystyle \frac{{{GM}}_{\star }\,r}{{\left({r}^{2}+{z}^{2}\right)}^{3/2}}+\displaystyle \frac{1}{{\rho }_{\mathrm{gas}}}\displaystyle \frac{\partial {P}_{\mathrm{gas}}}{\partial r}.\end{eqnarray} \tag{ C6 }$

Radiative transfer assumes local thermodynamic equilibrium with CO level populations set by the Boltzmann equation and a source function approximated by the Planck function. The emergent line shape is a Gaussian with a width of $\sqrt{2{k}_{{\rm{B}}}\,T\,(r,z)/{m}_{\mathrm{CO}}}$ . The modeling package includes optional parameters for disk turbulence, a parameterized disk wind, and dust radiation, which are not utilized here.

Sites of Planet Formation in Binary Systems. I. Evidence for Disk−Orbit Alignment in the Close Binary FO Tau

Article metrics

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. Observations and Data Reduction

2.1. ALMA

2.1.1. Observations

2.1.2. Calibration and Imaging

2.2. Keck—NIRSPEC behind AO

2.3. High Angular Resolution Imaging

2.4. Time-series Photometry

2.5. Literature Photometry and Astrometry

3. Stellar and Binary Orbital Parameters

3.1. Stellar Characterization

3.1.1. Modeling High-resolution NIR Spectra

3.1.2. Modeling the Spectral Energy Distribution

3.2. Relative Stellar Radial Velocity

3.3. Binary Orbital Solution

3.4. Photometric Variability

4. Disk Properties

4.1. 1.3 mm Continuum

4.2. ¹²CO J = 2–1

4.3. Modeling the CO Gas Disk

5. Discussion

5.1. The Binary−Disk Interaction

5.1.1. Disk−Orbit Alignment

5.1.2. Disk Truncation

5.1.3. Extended CO Emission

5.1.4. Future Submillimeter Observations

6. Summary and Conclusions

Acknowledgments

Appendix A: Binary Orbital Parameter Posteriors

Appendix B: Continuum Disk Modeling

Appendix C: CO Disk Modeling and Parameter Posteriors

Footnotes

Sites of Planet Formation in Binary Systems. I. Evidence for Disk−Orbit Alignment in the Close Binary FO Tau

Article metrics

Share this article

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. Observations and Data Reduction

2.1. ALMA

2.1.1. Observations

2.1.2. Calibration and Imaging

2.2. Keck—NIRSPEC behind AO

2.3. High Angular Resolution Imaging

2.4. Time-series Photometry

2.5. Literature Photometry and Astrometry

3. Stellar and Binary Orbital Parameters

3.1. Stellar Characterization

3.1.1. Modeling High-resolution NIR Spectra

3.1.2. Modeling the Spectral Energy Distribution

3.2. Relative Stellar Radial Velocity

3.3. Binary Orbital Solution

3.4. Photometric Variability

4. Disk Properties

4.1. 1.3 mm Continuum

4.2. 12CO J = 2–1

4.3. Modeling the CO Gas Disk

5. Discussion

5.1. The Binary−Disk Interaction

5.1.1. Disk−Orbit Alignment

5.1.2. Disk Truncation

5.1.3. Extended CO Emission

5.1.4. Future Submillimeter Observations

6. Summary and Conclusions

Acknowledgments

Appendix A: Binary Orbital Parameter Posteriors

Appendix B: Continuum Disk Modeling

Appendix C: CO Disk Modeling and Parameter Posteriors

Footnotes

4.2. ¹²CO J = 2–1