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

Close binary systems present challenges to planet formation. As binary separations decrease, so too do the occurrence rates of protoplanetary disks in young systems and planets in mature systems. For systems that do retain disks, their disk masses and sizes are altered by the presence of the binary companion. Through the study of protoplanetary disks in binary systems with known orbital parameters, we seek to determine the properties that promote disk retention and, therefore, planet formation. In this work, we characterize the young binary-disk system, FO Tau. We determine the first full orbital solution for the system, finding masses of $0.35^{+0.06}_{-0.05}\ M_\odot$ and $0.34\pm0.05\ M_\odot$ for the stellar components, a semi-major axis of $22(^{+2}_{-1})$ AU, and an eccentricity of $0.21(^{+0.04}_{-0.03})$. With long-baseline ALMA interferometry, we detect 1.3mm continuum and $^{12}{\mathrm{CO}} \ (J=2-1)$ line emission toward each of the binary components; no circumbinary emission is detected. The protoplanetary disks are compact, consistent with being truncated by the binary orbit. The dust disks are unresolved in the image plane and the more extended gas disks are only marginally resolved. Fitting the continuum and CO visibilities, we determine the inclination of each disk, finding evidence for alignment of the disk and binary orbital planes. This study is the first of its kind linking the properties of circumstellar protoplanetary disks to a precisely known binary orbit. In the case of FO Tau, we find a dynamically placid environment (coplanar, low eccentricity), which may foster its potential for planet formation.


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 path-ways of binary systems that foster planet formation remain largely unknown.Sub-mm observations at high angular resolution, possible with the 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-)mm continuum emission (tracing ∼mm size dust grains) in young binaries find that beyond separations of a few 100 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. 2020Zurlo et al. , 2021) ) as the binary projected separation decreases.Large dust grains can be biased tracers of disk properties due to processes like radial drift (Weidenschilling 1977;Ansdell et al. 2018), but the trend in disk radii follows the general predictions of binarydisk 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).
Sub-mm 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 reproduce this difference (Zagaria et al. 2021a,b) 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 (Dupuy et al. 2022;Behmard et al. 2022;Christian 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 10s of AU.At the distances of the closest star-forming regions (∼ 140 AU), decades of highangular-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 12 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 milliarcsecond (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.

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 12 CO J = 2 − 1 transition (230.538GHz) 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 −1 ).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.1 ′′ and 2.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.034 ′′ and maximum recoverable scale of 0.83 ′′ (roughly the C8 configuration).Observations in the compact and extended configurations took place on 2021-07-18-11:17 and 2021-08-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.

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 fullycalibrated 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 the program's data release web page 1 .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 −1 of the 12 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 (SNR).We repeated this process, reducing the solution interval until one of the following conditions was met: the SNR does not improve from the previous iteration, the number of 1 https://almascience.eso.org/almadata/lp/DSHARP/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 SNR 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 SNR to 85.
We then combined the extended and self-calibrated compact measurement sets and repeated this process.Phaseonly self-calibration was performed for the full integration, 1500 s, 1000 s, and 600 s, bringing the SNR 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 three sigma cleaning threshold.This resulted in a synthesized continuum beam of 45 × 22 mas with a position angle of 22 • .The root-mean-square deviation (RMS) 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 12 CO J = 2 − 1 data cubes for the compact and extended configurations, including channels within ±20 km s −1 from FO Tau's systemic velocity.The same selfcalibration steps above were applied, followed by continuum subtraction.Here we adopted the 2-channel spacing of 0.635 km s −1 .The resultant combined cube was imaged with interactive masking using Briggs robust = 1 weighting (which is larger than the continuum imaging due to a lower average SNR in the CO cube).The synthesized CO beam is 79 × 36 mas with a position angle of 25 • .The CO channel maps have a RMS values of 0.9 mJy/beam and are cleaned to a depth of three sigma.These maps are presented in Section 4.2 and 4.3.

Keck -NIRSPEC behind AO
FO Tau was observed on UT 2010 December 12 (UTC 09:32:18.66)at the Keck II 10-meter telescope using NIR-SPEC (McLean et al. 2000) behind the adaptive optics (AO) system.The spectrograph slit position angle was set to 230 degrees 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 re-images the instrumental dimensions onto the detector reduced by a factor of ∼11, thus the post-AO slit width×length used was 0. ′′ 027 × 2. ′′ 26.At the time of observations, the airmass was 1.01 and the natural seeing was 0. ′′ 5, facilitating robust correction from the AO system.The AO frame rate for FO Tau was 149 Hz and wavefront 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 Mauna Kea 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.
Internal comparison lamp lines were used to obtain a second order dispersion solution; median filtered stacks of ten internal flat lamp and dark frame exposures yielded a master flat and dark.A continuum SNR 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.

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 Oct 20, 2019 Jan 20, and 2022 Oct 19.On each night, we obtained sets of 12 images dithered across the detector in the Hcont and Kcont filters.Each image consisted of 1.0 sec exposures with 10 co-added 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 airmasses 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 Oct 20, the shape of the AO corrected images varied over time because of an oscillation problem with the Keck II telescope (C.Alvarez; 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) and applied a plate scale of 9.971 ± 0.004 mas pixel −1 and subtracted 0. • 262 ± 0. • 020 from the measured position angles.The positions are av-eraged 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, Hcont, and Kcont filters on UT 2009 Nov 21.We obtained two images in a single dither position, each consisting of 40 co-added exposures of 0.5 sec each.The seeing was better than average (∼0.4 ′′ ) and the system was observed at moderate airmass (∼1.55), yielding diffractionlimited 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 was fit within each image by testing the χ2 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, 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 2 .

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 pre-search data conditioning simple aperture photometry (PDCSAP; Smith et al. 2012;Stumpe et al. 2012;Stumpe et al. 2014) light curves.All of the TESS data used in this paper can be found in MAST: 10.17909/zpka-v291.

Literature Photometry & Astrometry
Our analysis makes use of photometric measurements from the literature and astrometry from Gaia DR3.The relevant measurements are included in Table 3.   (Rizzuto et al. 2018;Belokurov et al. 2020;Bryson et al. 2020).In disk bearing stars, this thresh-old 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 +4 −3 pc, while clustering analyses of Taurus members places FO Tau's sub-group center at ∼130 pc with a radial dispersion of a few pcs (Esplin & Luhman 2019;Krolikowski et al. 2021).Given that these values are in fair agreement, we adopt the geometric distance in our analysis.

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.

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 , vsini, 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 , vsini, 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); 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 non-magnetic 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 prep, private communication;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.

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 ∼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 HR 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 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, lowresolution 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 corre-  3. ) sponding 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 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 magnitudes of extinction using the Cardelli et al. (1989) law, which are added in quadrature to the quoted er-rors.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 auto-correlation.The fit ends once the chain auto-correlation changes by less than 5% or the number of steps exceeds the auto-correlation time by a factor of 50.The first five auto-correlation 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 +50 −60 K and 3380 +50 −60 K. The uncertainties represent the 68% posterior confidence interval, which only includes formal uncertainties.The fit returns radii of 1.48 +0.10 −0.09R ⊙ and 1.42 ± 0.09 R ⊙ for the primary and secondary, respectively.With these posteriors we derive bolometric luminosities of 0.25 ± 0.02 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 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.1are 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 HR 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 (σ ∆|log10τ | < 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.18  −0.19 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 also assess their agreement with dynamical mass measurements from the literature.
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 (Simon et al. 2019;David 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 HR diagram locations.Rizzuto et al. (2016Rizzuto et al. ( , 2020)), for instance, find that even the standard DSEP models over-predict dynamical measurements of binary total masses by 25% to 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 short comings of standard models may simply conspire to provide more accurate mass predictions from the HR 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., 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 over-predict the masses of young, low mass stars.

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 vsini.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 to 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 −1 .The negative value indicates that FO Tau B is moving toward Earth compared to FO Tau A.

Binary Orbital Solution
We fit the FO Tau orbital solution using the orvara python package (Brandt et al. 2021).The fit uses a paralleltempering MCMC framework (Foreman-Mackey et al. 2013;Vousden et al. 2016) with 10 temperatures, 100 walkers, each taking 5×10 5 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 √ e sin ω A and √ e cos ω A .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 ) hyper-parameters, 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 −1 , 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 in 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., further 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), as well as a corner plot (Figure 12) are provided in Appendix A. 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.

Photometric Variability
Figure 6 presents spatially unresolved FO Tau light curves from three TESS sectors.The variability is dominated by quasi-periodic 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 d, 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 auto-correlation 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-12 d (e.g., returning P rot /2; Newton et al. 2022;Magaudda et al. 2022).
We do not find a single period across the temporal regions analyzed that has a consistently dominant LSP peak with positive auto-correlation, regardless of the reduction.The ∼4.2 d peak, for instance, is only present in the Sector 19 PDC-SAP light curve.The strongest LSP peak is found in Sector 43 at 5.2 d, 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 d are present in multiple time ranges (Sector 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 d and a single ACF peak at 5.0 d.The present analysis is suggestive of some astrophysical periodicity between 4 and 6 d; however, we cannot confidently assign it an origin.Assuming some component of the periodicity is stellar rotation, for the stars' vsini and radius values, this period range would correspond to inclinations between ∼ 50 • to 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 the stars inherit the binary orbital inclination, the rotation period would be ∼ 2.9 d, which corresponds to negative ACF values (anti-correlation).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.

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 ob-servation 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 imageplane resolution, we forward model the continuum visibilities with an exponentially tapered power-law intensity profile, e −(r/Rc) γ 2 . (1) Here, γ 1 is power-law index and γ 2 is the exponential-taper index at the transition radius, R c .This model describes disks with sharply decreasing outer profiles which 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 uv 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 bandwidthweighted 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 imageplane 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.5AU, the inclinations are 27.3±0.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: where F ν is the 230 GHz flux density, κ ν is the dust opacity (2.3 cm 2 g −1 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 et al. 2013, assumes optically thin dust).Including an additional 5% uncertainty on the ALMA absolute flux calibration, we compute dust masses of 5.4(±0.4)× 10 −6 M ⊙ and 5.1(±0.4)×10−6 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 non-negligible dust optical depths (e.g., Huang et al. 2018;Birnstiel 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. 12CO J = 2-1 Figure 4.1 presents maps of 12 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 (blue shifted) side of the primary disk, and the southeast (red shifted) side and center of the secondary disk.The affected velocity channels correspond to ∼5 and 8 km s −1 (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 red-shifted 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.

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) 3 .Given our focus on the bulk disk properties, we provide a brief overview of the model here, but refer interested readers to the citation above or Appendix C for more detail.The model assumes a temperature and gasdensity 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 twelve 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 uv 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 uv baselines and velocity channels using the galario package.The model visibilities are fit to ALMA observations in a MCMC framework using emcee.The fit is made to velocity channels between −3.0 and 17 km s −1 , ignoring channels between 5 and 8 km s −1 that have foreground cloud absorption.In total, 27 channels are used for the fit, each with a velocity width of 0.635 km s −1 .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 auto-correlation 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 −1 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 +4 −3 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.

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.

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.
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: (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.Counter clockwise 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, scatteredlight 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 counter clockwise with the binary orbit, and a misaligned solution where the disk rotates clockwise, retrograde with the binary orbit.We compute the diskorbit 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 grey, 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.
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 where α is the viscosity parameter, h out is the aspect ratio at the outer disk edge, M 1 and M 2 are the primary and secondary stellar masses, respectively, a is the semi-major axis, and r out is the outer disk edge.(Here, we have assumed the dimensionless viscous coefficient, V b ≈ 1; see  an alternate projection where the difference between the solutions is more readily apparent.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 the light of numerical simulations by Picogna & Marzari (2013) that find disks are dynamically heated by companions, especially in the vertical direction.The present analysis does not distinguish between formation or subsequent evolution as the source of FO Tau's diskorbit alignment.A larger sample, particularly at larger sepa-rations (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 diskorbit alignment is common in protoplanetary disk-hosting binaries with small separations.

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 semi-major axis.Manara et al. (2019) present an equation for the disk truncation radius following the work of Artymowicz & Lubow (1994), which has the form, R trun = 0.49aq −2/3 0.6q −2/3 + ln(1 + q −1/3 ) be c + 0.88µ 0.01 .
(5) Here, a is the semi-major 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 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 sub-section, 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 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 gas to dust disk radius ratios 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 (Miranda & Lai 2015;Lubow et al. 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.

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 −1 ).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 4.1).
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 (Zhang et al. 2023;Weber 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-)mm 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.

Future Sub-mm Observations
Despite the challenges in observing compact binaries, FO Tau remains one of the best candidates to study the binarydisk 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.

SUMMARY & 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-angularresolution imaging, angularly-resolved high-resolution NIR spectra, and (sub-)mm 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 semi-major axis of 22 AU.The stars are near-equal mass, with dynamical masses of 0.35 M ⊙ 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 12 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.3mm fluxes for each component are near equal at 2.96 ± 0.07 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 uv 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 (Dupuy et al. 2022;Christian et al. 2022;Behmard 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.BMT would like to thank Sam Factor, Kendall Sullivan, Neal Evans, and Stella Offner for helpful conversations.BMT is supported by the Heising-Simons Foundation's 51 Pegasi b Postdoctoral Fellowship in Planetary Astronomy.
GHS 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 NOIR-Lab (PropID: 2022B-970020; PI: G. Schaefer) supported by the NSF Mid-Scale Innovations Program.LP 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 part-  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 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 2k B T (r, z)/m CO .The modeling package includes optional parameters for disk turbulence, a parameterized disk wind, and dust radiation, which are not utilized here.

Figure 1 .
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.

Figure 2 .
Figure 2. Two-component SED fit to FO Tau photospheric flux measurements.The primary and secondary best-fit BT-Settl CFIST spectra are shown in blue and orange, respectively, with a combined spectrum in black.Model spectra have been smoothed to a resolution of R ∼ 3000 for clarity.Red points present measurements of the photospheric flux as they are fit.Horizontal bars represent the FWHM of the filter or spectral region considered in the measurement.Leftmost are spatially resolved HST-STIS flux measurements at 6100 Å from Hartigan & Kenyon (2003), highlighted in the top-left subplot.The point at 7510 Å is the unresolved flux from Herczeg & Hillenbrand (2014).The lower subplot presents the model agreement with the J-band contrast determined in this work.The ∆J magnitude of the best-fit model is shown as the black point.Grey lines highlight the J-band spectral region.Rightmost are the H-band photospheric fluxes determined in this work.(Photometry values are provided in Table3.) vsys (km s −1 LSRK) 8.2 +0.6 −0.5 6 ± 2 This Work NOTE-(a) Spectroscopically determined; additional uncertainty from systematic and astrophysical effects are included in the parenthetical (see Section 3.1.1).(b) Determined from SED fitting in Section 3.1.2.

Figure 3 .
Figure 3. On-sky projection of the FO Tau orbit.Observations are shown in 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 furthest from Earth in the South-West portion of the orbit and closest in the North-East.

Figure 4 .
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.

Figure 5 .
Figure 5.The FO Tau orbit fit.As a function of time, we present the projected separation, position angle, and radial-velocity difference from top to bottom, each with its associated residuals.Astrometry measurements are show as blue circles.Radial-velocity difference measurements are shows as 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.

Figure 6 .
Figure 6.TESS PDCSAP light curves of FO Tau.The variability is dominated by stochastic bursts associated with accretion and inner disk related processes.

Figure 7 .
Figure 7.The FO Tau 1.3 mm continuum map.Circumstellar dust is detected, but not clearly resolved, toward each of the FO Tau binary components.The highest-likelihood orbital solution is overplotted in blue.The synthesized beam is presented in the bottom left corner and a 10 AU scale bar is shown in the bottom right.

Figure 8 .
Figure 8. Imaging of the FO Tau 12 CO J = 2 − 1 visibilities.Left: The CO integrated intensity.All emission shown is ≥3σ above the background noise.Absorption from intervening cloud material contaminates the north-west side of the primary disk and the south-east side and center of the secondary disk.The synthesized beam is presented in the bottom left and a 10 AU scale bar is provided in the bottom right.The highest-likelihood orbital solution is over-plotted in blue.White points mark the peak of continuum emission.Right: The CO intensity weighted velocity map using the bettermoments quadratic method (Teague & Foreman-Mackey 2018).Extended CO gas disks are observed toward each component.Velocity channels heavily affected by cloud absorption are shown in the grey hatched region of the velocity color bar.Black contours are the 1.3 mm continuum emission, beginning at 5σ.The CO and continuum synthesized beams are presented in the white and black ellipses, respectively.In both panels a Keplerian mask, based on our best-fit disk model, is employed to reduce background noise for clarity (see Section 4.2).

Figure 9 .
Figure 9. 12 CO J = 2 − 1 velocity channel maps (LRSK radio) of the observations (top panel set), the best-fit model (middle panel set), and the data−model residuals (bottom panel set).Each panel provides the channel's central velocity in the top left corner.White points signify peaks in the continuum emission.Channels contaminated by cloud absorption have a diagonal white line.

Figure 10 .
Figure10.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 Figure11.The grey curve in the right panel corresponds to the obliquity distribution of a randomly oriented disk.

Figure 11 .
Figure11.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 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 counter clockwise 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.
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.

Figure 13 .
Figure 13.Continuum visibility fitting.Results are shown for SPW 3. Left: Data and model visibilities as a function of uv-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: Best fit radial intensity profile of the primary (blue) and secondary (orange).

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.

Table 1 .
High-Angular Resolution Orbit Monitoring

Table 2 .
Flux ratios for FO Tau from Keck adaptive optics observations

Table 3 .
Known and Derived Properties of FO Tau