# A Coplanar Circumbinary Protoplanetary Disk in the TWA 3 Triple M Dwarf System

Ian Czekala, Álvaro Ribas, Nicolás Cuello, Eugene Chiang, Enrique Macías, Gaspard Duchêne, Sean M. Andrews, Catherine C. Espaillat

► To cite this version:

Ian Czekala, Álvaro Ribas, Nicolás Cuello, Eugene Chiang, Enrique Macías, et al. A Coplanar Circumbinary Protoplanetary Disk in the TWA 3 Triple M Dwarf System. The Astrophysical Journal, 2021, 912, 10.3847/1538-4357/abebe3. insu-03705311

HAL Id: insu-03705311

https://insu.hal.science/insu-03705311v1

Submitted on 2 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A Coplanar Circumbinary Protoplanetary Disk in the TWA 3 Triple M Dwarf System

Ian Czekala¹,²,³,⁴,⁵,<fnref n="12" /> ID, Álvaro Ribas⁶, Nicolás Cuello⁷ ID, Eugene Chiang⁵,⁸ ID, Enrique Macías⁶,⁹ ID, Gaspard Duchêne⁵,⁷ ID, Sean M. Andrews¹⁰ ID, and Catherine C. Espaillat¹¹ ID

¹ Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA; iczekala@berkeley.edu

² Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

³ Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

⁴ Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA

⁵ Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA

⁶ European Southern Observatory (ESO), Alonso de Córdova 3107, Vitacura, Casilla 19001, Santiago de Chile, Chile

⁷ Univ. Grenoble Alpes, CNRS, IPAG (UMR 5274), F-38000 Grenoble, France

⁸ Department of Earth and Planetary Science, University of California at Berkeley, Berkeley, CA 94720-4767, USA

⁹ Joint ALMA Observatory, Alonso de Córdova 3107, Vitacura, Santiago, Chile

¹⁰ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

¹¹ Department of Astronomy, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA

Received 2021 February 21; revised 2021 March 1; accepted 2021 March 2; published 2021 April 29

## Abstract

We present sensitive ALMA observations of TWA 3, a nearby, young (~10 Myr) hierarchical system composed of three pre-main-sequence M3–M4.5 stars. For the first time, we detected ¹²CO and ¹³CO J = 2–1 emissions from the circumbinary protoplanetary disk around TWA 3A. We jointly fit the protoplanetary disk velocity field, stellar astrometric positions, and stellar radial velocities to infer the architecture of the system. The Aa and Ab stars (0.29 ± 0.01 M<sub>☉</sub> and 0.24 ± 0.01 M<sub>☉</sub>, respectively) comprising the tight (P = 35 days) eccentric (e = 0.63 ± 0.01) spectroscopic binary are coplanar with their circumbinary disk (misalignment <6° with 68% confidence), similar to other short-period binary systems. From models of the spectral energy distribution, we found the inner radius of the circumbinary disk (r<sub>inner</sub> = 0.50–0.75 au) to be consistent with theoretical predictions of dynamical truncation r<sub>cav</sub>/a<sub>inner</sub> ≈ 3. The outer orbit of the tertiary star B (0.40 ± 0.28 M<sub>☉</sub>, a ~ 65 ± 18 au, e = 0.3 ± 0.2) is not as well constrained as the inner orbit; however, orbits coplanar with the A system are still preferred (misalignment < 20°). To better understand the influence of the B orbit on the TWA 3A circumbinary disk, we performed SPH simulations of the system and found that the outer edge of the gas disk (r<sub>outer</sub> = 8.5 ± 0.2 au) is most consistent with truncation from a coplanar, circular, or moderately eccentric orbit, supporting the preference from the joint orbital fit.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Binary stars (154); Trinary stars (1714); Stellar dynamics (1596); Pre-main-sequence stars (1290)*

## 1. Introduction

The distribution of pre-main-sequence multiple system architectures informs our understanding of the mechanisms that govern star and planet formation. Recently, Czekala et al. (2019) found that the degree of alignment between the disk and its host binary (the mutual inclination, $\theta$) is a strong function of the orbital period. Circumbinary disks around short-period binaries (P < 40 days) are preferentially coplanar, while disks around longer period binaries exhibit a wide range of mutual inclinations, including polar configurations ($\theta \approx 90°$). These trends might be manifestations of the same physical mechanisms that produce close binaries. Because it is difficult to directly fragment stars on scales <5 au (Larson 1969; Bate et al. 2002; Offner et al. 2016), it is believed that tight binaries are instead produced from wider binaries that have hardened through star-disk interactions (Offner et al. 2010; Bate 2019). What mechanisms mediate this evolution, whether migration affects mutual inclination, and whether initial mutual inclination affects migration efficiency are all open questions.

Typically, a nearly coplanar disk surrounding a low eccentricity binary will precess around the binary angular momentum vector as dissipative forces damp the angular momentum vectors of the disk and the binary into alignment (Foucart & Lai 2013). But if

the binary is sufficiently eccentric, then the disk can precess around the eccentricity vector<fnref n="13" /> (Aly et al. 2015; Martin & Lubow 2017; Zanazzi & Lai 2018; Cuello & Giuppone 2019) and access polar mutual inclinations. Indeed, highly misaligned disks are preferentially found around highly eccentric (e > 0.7) binaries (Czekala et al. 2019; Kennedy et al. 2019). There are, however, several disks around eccentric short-period binaries that are coplanar, such as AK Sco, DQ Tau, and UZ Tau E (Czekala et al. 2015, 2016, 2019). Our understanding of how these period and eccentricity trends interrelate is limited by the sample size of circumbinary disk systems with well-measured architectures.

The pre-main-sequence system TWA 3 represents an opportunity to expand this sample to aid in the interpretation of binary formation and evolution mechanisms. TWA 3 consists of three young (10 ± 3 Myr; Bell et al. 2015), pre-main sequence M3–M4.5 stars in a hierarchical configuration; their spectral types correspond to ~0.3 M<sub>☉</sub> (Herczeg & Hillenbrand 2014; Tofflemire et al. 2019). The inner Aa–Ab binary has an orbital period of P = 34.8785 ± 0.0009 days, an eccentricity of e = 0.628 ± 0.006, and spectral types of M4 and M4.5, respectively (Kellogg et al. 2017). The Gaia DR2 parallax is $\sigma = 27.31 \pm 0.12$ mas (including a 0.02 mas

<footnote n="12">¹² NASA Hubble Fellowship Program Sagan Fellow.</footnote>

<footnote n="13">¹³ The eccentricity vector is drawn from binary apoapse to periapse.</footnote>
systematic term, Lindegren et al. 2018), corresponding to a
distance of 36.62 ± 0.16 pc (Bailer-Jones et al. 2018; Gaia
Collaboration et al. 2018). Tofflemire et al. (2019) noted that the
A and B components suffer from significant excess astrometric
noise—possibly due to photometric variability—but the parallax
distances for each source are consistent with each other and the
A–B orbit arcs in Kellogg et al. (2017). Time-series photometry
and emission-line spectroscopy revealed that accretion from the
circumbinary disk to the inner binary is phased with periastron,
and that material is preferentially accreted onto the primary star,
TWA 3Aa (Tofflemire et al. 2017, 2019). The gradual movement
of the outer triple companion (projected separation 1.′′55 or
57 au; Tokovinin et al. 2015) over a ~20 yr baseline suggests an
orbital period of ~200–800 yr (Kellogg et al. 2017).

Andrews et al. (2010) used the Submillimeter Array (SMA) to localize the submillimeter emission in the TWA 3 system to the A binary, measuring a flux density of 75 mJy at 340 GHz. As demonstrated by fits to the deprojected and azimuthally averaged baselines, the circumbinary disk itself was only marginally resolved (1.′′11 × 0.′′74 beam) but found to have a radius of ~20 au. Andrews et al. (2010) did not detect ¹²CO J = 3–2 emission to an upper limit of 0.6 Jy beam⁻¹ integrated over a 0.7 km s⁻¹ channel. Based upon the fit of an elliptical Gaussian to the visibilities, Andrews et al. (2010) derived a disk inclination of iₗdisk = 36° ± 10° (relative to the sky plane) and disk orientation of Ωₗdisk = 169° ± 15° (the position angle of the ascending node measured east of north).<fnref n="14" />

More recent orbital solutions suggested that the inner binary
orbit, circumbinary disk, and outer tertiary orbit may be
misaligned. Kellogg et al. (2017) combined an astrometric
observation of the inner binary (Anthonioz et al. 2015) with
their double-lined radial velocity solution to constrain the
position angle of the ascending node Ω<sub>inner</sub> [∈ 93°, 123°] and
inclination i<sub>inner</sub> [∈ 32°, 63°] or i<sub>inner</sub> [∈ 118°, 149°]. These
orbital parameters, together with the disk parameters reported
in Andrews et al. (2010), suggested that the planes of the
spectroscopic binary and the circumbinary disk were mis-
aligned by at least θ ~ 30° (Kellogg et al. 2017).

We acquired ALMA observations of the TWA 3A circum-
binary disk to better understand its size and orientation relative
to the stellar orbits. In Section 2 we describe the ALMA
observations and data reduction. In Section 3 we dynamically
model the gas rotation curve of the TWA 3A circumbinary
disk, fit the spectral energy distribution (SED) of TWA 3A, and
perform a joint stellar orbital fit to the radial velocities and
astrometric measurements of the Aa, Ab, and B stars. In
Section 4 we describe our smoothed-particle hydrodynamics
(SPH) simulations of the TWA 3 system and discuss how they
support our interpretation of the circumbinary disk as nearly
coplanar and dynamically truncated, both internally and
externally, by binary companions. We also briefly review
similar analog protoplanetary and exoplanetary systems in our
discussion of mutual inclination and disk truncation. We
conclude in Section 5.

## 2. Data

We obtained deep Atacama Large Millimeter Array (ALMA)
observations of the TWA 3 system in 2018. We used a
correlator setup that assigned two 2 GHz wide spectral windows

to the dust continuum (centered on 220 and 232 GHz) and two
spectral windows at 122 kHz (0.16 km s⁻¹) resolution to target
the ¹²CO and ¹³CO J = 2 − 1 transitions. Two sets of observa-
tions (project code 2018.1.01545.S) were executed on 2018
October 16th and 27th (JD 2458408.1404 and 2458419.0483)
using 44 and 43 antennas of the main array, respectively.
The array was similarly configured for each observation, with
baselines ranging from 15 m to 2.4 km. Both observing sequences
used the J1107−4449 quasar as an amplitude and bandpass
calibrator, and used the J1126−3828 quasar as a phase calibrator.
Each execution spent 47.4 minutes on-source, for a total on-source
time of 1 hr 34.8 m. The mean precipitable water vapor for each
observation was 1.9 mm and 0.4 mm, respectively.

We began our data reduction with the pipeline-calibrated
measurement set provided by ALMA/NAASC staff. We used
the CASA 5.4 (McMullin et al. 2007) facility software and
followed common calibration and imaging procedures (e.g.,
the DSHARP reduction scripts, Andrews et al. 2018),<fnref n="15" /> the
pertinent details of which we now describe. To assess the
quality of each execution block we first reduced each
observation individually. We excised the channels with line
emission to create a continuum-only measurement set with a
total bandwidth of 4.6 GHz. We performed an initial round
of continuum imaging using the CASA tclean task with
robust = 0.5, an image size of 512 × 512 pixels, 0.′′015 pixel
size, deconvolver = `''multiscale''`, scales = [0,
15, 30, 45, 75] pixels, a threshold of 0.6 mJy, and an
elliptical mask with position angle 110°, semimajor axis of 0.′′45,
and semiminor axis of 0.′′36.

We fit the continuum emission with an elliptical Gaussian
using the `imfit` and `uvmodel` tasks. We found excellent
astrometric agreement between executions, with the emission
centroid located at ICRS 11:10:27.731–37.31.51.84, coincident
with the Gaia position of TWA 3A to within 0.′′05. We found
adequate agreement between the total continuum fluxes (31.9 mJy
and 36.1 mJy, respectively), only slightly more different than the
expected 10% amplitude calibration uncertainty.<fnref n="16" />

We proceeded to self-calibrate the combined measurement
set through a series of applications of tclean to the
continuum visibilities using threshold depths {0.6, 0.15, and
0.15} mJy beam⁻¹ interleaved with applications of the
gaincal and applycal CASA tasks using spectral-window
dependent solves with intervals {60 s, 30 s, and 18 s}. We then
cleaned to a final depth of 0.02 mJy beam⁻¹ and performed one
round of phase and amplitude gain solutions over the scan
length (8 minutes, solint = `inf`). We monitored the peak
flux, total flux, and rms of the images throughout the process
(Brogan et al. 2018), and found the peak continuum S/N
improved to 2130 from an initial value of 260.

We then used the applycal task to apply the self-
calibration solutions to channels containing the spectral line
observations<fnref n="17" /> (without propagating flags for failed solutions,
applymode = calonly). We estimated the continuum
from nearby line-free channels and subtracted it from the
spectral line observations using the uvcontsub CASA task.

<footnote n="14">¹⁴ The lack of a gas detection meant degenerate “flipped” solutions were also valid: $i_{\text{disk}} = 144^{\circ} \pm 10^{\circ}$ and $\Omega_{\text{disk}} = 349^{\circ} \pm 15^{\circ}$.</footnote>

<footnote n="15">¹⁵ https://bulk.cv.nrao.edu/almadata/lp/DSHARP/</footnote>

<footnote n="16">¹⁶ See the ALMA Technical Handbook (Remijan et al. 2019) and ALMA Memo 594 https://science.nrao.edu/facilities/alma/aboutALMA/Technology/ALMA_Memo_Series/alma594/memo594.pdf.</footnote>

<footnote n="17">¹⁷ Though part of a standard self-calibration workflow, we noticed that the S/N of the line channels did not measurably improve after applying the self-calibration solutions. We believe that this is because the fine resolution channels are thermal-noise dominated and not limited by residual phase errors.</footnote>
<table><caption>**Table 1**
Image Properties</caption><thead><tr><th></th><th>Beam Dimensions, P.A.</th><th>rms (mJy beam<sup>-1</sup>)</th></tr></thead><tbody><tr><td>226 GHz cont.</td><td>0"24 × 0"19, -64°</td><td>0.013</td></tr><tr><td><sup>12</sup>CO J = 2–1</td><td>0"24 × 0"19, -64°</td><td>1.1</td></tr><tr><td><sup>13</sup>CO J = 2–1</td><td>0"25 × 0"20, -63°</td><td>1.3</td></tr></tbody></table>

Note. The rms noise levels for the spectral line cubes correspond to the values per 0.8 km s⁻¹ channel. All images were synthesized with robust = 0.5.

The line channels were imaged using the tclean task with the auto-multithresh masking algorithm with cell = `0.015arcsec''`, gain = 0.1, deconvolver = `multiscale''`, scales = `[0, 10, 30, 100, 200, 300]` pixels, robust = 0.5, and deconvolved to a depth of threshold = `0.1mJy''`. The beam dimensions (FWHM) and image-plane rms are described in Table 1.

We summed the pixels within the continuum CLEAN mask to measure a total continuum flux of 36.7 mJy. The dust continuum emission is compact, with nearly all of the flux contained within the central beam. We fit elliptical Gaussians to the continuum emission with the uvmodelfit and imfit tasks and derived FWHM dimensions of (0"12 × 0"12) and (0"16 × 0"10), respectively. Since these dimensions are on the order of the beam size, to better quantify whether we resolved the dust continuum, we also deprojected and azimuthally averaged the continuum visibilities using the uvplot package (Tazzari 2017), as shown in Figure 2. We used deprojection values of $i_{disk} = 49°$ and $\Omega_{disk} = 116°$, which were derived from our modeling effort of the $^{12}$CO J = 2–1 line emission, since the more extended emission provided better constraints for these parameters than the dust continuum emission (see Section 3.1). The declining visibility profile demonstrates that the outer extent of the disk is resolved: if the deprojected continuum emission profile is represented by a Gaussian, then a Fourier-domain FWHM $\approx 1200 k\lambda$ implies that the image-plane Gaussian has FWHM $\approx 0''15$, or 5.5 au.

In preliminary line imaging, we identified $^{12}$CO emission from [-9, 10] km s⁻¹(LSRK). To save computational complexity, we used the task mtransform to average the channels to 0.8 km s⁻¹ width. The $^{12}$CO emission was strongly detected across this velocity range (peak channel S/N = 30, see Figure 3, top panel); the $^{13}$CO emission was detected at lower but still significant levels (peak channel S/N = 6, see Appendix). We used the immoments CASA task to sum all flux in each channel across the velocity dimension, producing the moment maps in the third and fourth panels of Figure 1. We used the imstat CASA task to sum all flux within the CLEAN mask across the velocity and spatial dimensions, yielding integrated fluxes of 526 mJy km s⁻¹ and 86 mJy km s⁻¹ for $^{12}$CO and $^{13}$CO, respectively.

We noticed a central cavity in the moment maps of the continuum-subtracted line emission, corresponding to the location of peak dust continuum. The depression is most apparent in the $^{12}$CO emission but is also visible in the $^{13}$CO emission (Figure 1, third and fourth panels). To investigate whether this may be a continuum-subtraction artifact, we produced channel maps and a moment map for the non-continuum-subtracted $^{12}$CO spectral channels (Figure 1, second panel). This moment map does not exhibit a central cavity, suggesting that the feature seen in $^{12}$CO and $^{13}$CO is indeed a continuum-subtraction artifact. In protoplanetary disks, such an

artifact can arise when gas emission on the near side of the disk is optically thick and absorbs continuum emission originating from the disk midplane. When the continuum emission (estimated from channels offset in velocity from the line emission) is subtracted, most if not all of the line flux is also (erroneously) subtracted (for a full description of the effect, see Weaver et al. 2018). For this artifact to be significant, the continuum emission needs to have a brightness temperature comparable to the line emission, suggesting that the continuum emission is also optically thick or nearly so.

No noticeable continuum or gas emission is present in the system beyond the disk surrounding TWA 3A. We used an aperture approximately three times the area of the beam to extract continuum photometry at the location of B (see Section 3.3; Mason et al. 2018) and did not detect anything ($37 \pm 27 \mu Jy$).

# 3. Analysis

## 3.1. Dynamical Gas Analysis

Following the approach in Czekala et al. (2015), we constructed a forward model of the continuum-subtracted $^{12}$CO J = 2–1 visibilities to derive constraints on the disk architecture and velocity field. Briefly, a 3D model of the disk density, temperature, and velocity is parametrically defined and then ray-traced using RADMC-3D (Dullemond 2012) to produce image cubes. These cubes are Fourier-transformed and sampled at the spatial frequencies corresponding to the baselines of the ALMA measurement set to compute the likelihood of the data set.<fnref n="18" /> The posterior of the model parameters, defined by the data likelihood and any additional prior probability distributions, is explored using Markov Chain Monte Carlo (MCMC).

We initially explored models with self-similar surface density profiles $\Sigma(r)$ of the form described in Czekala et al. (2015) and references therein (see also Sheehan et al. 2019). We found that the gradual exponential taper of this profile at large radii was not well-matched to the TWA 3A $^{12}$CO J = 2–1 emission, which appears to decay quickly at large radii (see Figure 1). Drawing inspiration from dust continuum studies that flexibly parameterized the surface brightness profile $I(r)$ using the Nuker function (e.g., Tripathi et al. 2017), we adapted this to a surface density profile as

$$ \Sigma(r) = \Sigma_c \left(\frac{r}{r_c}\right)^{-\gamma} \left[1 + \left(\frac{r}{r_c}\right)^{\alpha}\right]^{(\gamma-\beta)/\alpha} \quad (1) $$

We followed Tripathi et al. (2017) and chose to sample in $\log_{10} \alpha$, imposing priors on shape parameters $\log_{10} \alpha$, $\beta$, and $\gamma$. We restricted $\log_{10} \alpha$ and $\beta$ with uniform probability to the ranges of $0 \le \log_{10} \alpha \le 2$ and $2 \le \beta \le 10$, respectively, and imposed a tapered prior on $\gamma$ of the form

$$ p(\gamma) = \frac{1}{1 + e^{-5(\gamma+3)}} - \frac{1}{1 + e^{-15(\gamma-2)}} \quad (2) $$

As discussed in Tripathi et al. (2017), these priors are practically motivated to allow a broad range of surface density profiles, including those with interior cavities and sharp outer

<footnote n="18">¹⁸ Using the DiskJockey.jl package (Czekala et al. 2015), https://github.com/iancze/DiskJockey.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** From left to right: the 226 GHz dust continuum, the channels centered on the ¹²CO emission (including the dust continuum) integrated over the full range of the ¹²CO emission, the continuum-subtracted ¹²CO emission, and the continuum-subtracted ¹³CO emission. The FWHM beam is shown in the lower left of each panel. Full channel maps for ¹²CO and ¹³CO are in Section 3.1 and the Appendix, respectively.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Continuum visibilities deprojected from the disk inclination. That the flux drops with increasing uv distance indicates the continuum emission is spatially resolved. The Gaussian profile with FWHM ≈ 1200 kλ implies an image-plane morphology with FWHM ≈ 0.15, or 5.5 au.</figcaption></figure>

edges, while restricting the parameter space sufficiently to avoid pathological surface density profiles ill-suited to protoplanetary disks. Because the Nuker profile increases the dimensionality of the posterior by adding two new parameters, we kept computational demands tractable by zeroing out the phase-center offsets $\delta_{\alpha}$, $\delta_{\delta}$ typically used with the standard model.

We sampled the posterior distribution using the DiskJockey package (Czekala et al. 2015) and assessed convergence both visually and by applying the Gelman-Rubin convergence diagnostic to independent chain ensembles. The resulting marginal posteriors on the disk model parameters are listed in Table 2. Most marginal distributions are well-described by Gaussians. The posteriors on the temperature profile exponent $q$ and Nuker shape parameters $\log_{10}\alpha$ and $\beta$ ran up against the range of their prior bounds, so 68% upper or lower confidence intervals are quoted for these parameters

<table><caption>**Table 2**
Inferred Disk Parameters</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>M<sub>A</sub> (M<sub>☉</sub>)</td><td>0.534 ± 0.010</td></tr><tr><td>r<sub>c</sub> (au)</td><td>6.8 ± 0.2</td></tr><tr><td>T<sub>10</sub> (K)</td><td>38 ± 1</td></tr><tr><td>γ</td><td>-3.5 ± 0.5</td></tr><tr><td>q</td><td>≤ 0.05</td></tr><tr><td>log<sub>10</sub> α</td><td>≥ 1.6</td></tr><tr><td>β</td><td>≥ 9.8</td></tr><tr><td>log<sub>10</sub> M<sub>disk</sub> log<sub>10</sub>(M<sub>☉</sub>)</td><td>-6.37 ± 0.09</td></tr><tr><td>ξ (km s<sup>-1</sup>)</td><td>0.60 ± 0.03</td></tr><tr><td>i<sub>disk</sub> (°)</td><td>48.8 ± 0.7<sup>a</sup></td></tr><tr><td>Ω<sub>disk</sub> (°)</td><td>116.5 ± 0.4</td></tr><tr><td>v<sub>r</sub> (km s<sup>-1</sup>)</td><td>1.22 ± 0.02<sup>b</sup></td></tr></tbody></table>

**Notes.** Using ¹²CO. The 1D marginal posteriors are well-described by a Gaussian, so we report symmetric error bars here (statistical uncertainties only). These parameters were inferred using a distance of $d = 36.62$ pc.

<sup>a</sup> Ambiguity with $i_{\text{disk}} = 131.2 \pm 0.7$.

<sup>b</sup> LSRK reference frame.

instead. Figure 3 shows a realization of the model and residual visibilities drawn from the posterior distribution, imaged in the same way as the data. The residual channel maps are broadly consistent with residual thermal noise, demonstrating that the synthesized model is a good fit for the data. We define the outer edge of the disk using the radius that contains 95% of the mass: $r_{\text{outer}} = 8.5 \pm 0.2$ au. That the 226 GHz continuum emission (Figure 1) is more compact than the extent of the gas emission is consistent with the expectation that radial drift has moved large dust grains (mm or cm sized) inward (Andrews 2020).

The inferred surface density profile constrains the peak density to approximately 7 au in radius; the $q \approx 0$ temperature exponent indicates that the gas in the ring is nearly constant temperature across its narrow radial extent. However, the gas-depleted central cavity implied by this surface density profile is not real, but instead reflects the continuum-subtraction artifact that removed ¹²CO $J=2-1$ emission from central radii. As shown in Figure 3, the flexible Nuker profile provided an excellent fit to the continuum-subtracted ¹²CO $J=2-1$ emission and allowed us to achieve our primary objectives of
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Continuum-subtracted ¹²CO J = 2 – 1 data, model, and residual channel maps. Model and residual visibilities were imaged the same way as the data. The velocity scale is labeled in the LSRK frame.</figcaption></figure>

inferring the disk velocity field and measuring the outer extent
of the circumbinary disk. A joint dust and gas model
simultaneously fit to the dust continuum and ¹²CO J = 2 – 1
emission could in principle recover a more accurate gas surface
density profile. However, the behavior of the surface density
profile at small radii is effectively a nuisance component to our
dynamical analysis and does not justify a more sophisticated
model, especially since its significantly expanded computa-
tional requirements would curtail our ability to thoroughly

explore the posterior distributions of key parameters (*M*<sub>A</sub>, *i*<sub>disk</sub>,
and Ω<sub>disk</sub>) via MCMC.

We constrained the disk position angle to be Ω<sub>disk</sub> = 116.5 ± 0.4, which is significantly different from the value found by Andrews et al. (2010; Ω<sub>disk</sub> = 144° ± 10°). We constrained the disk inclination to be either i<sub>disk</sub> = 48.8 ± 0.7 or i<sub>disk</sub> = 131.2 ± 0.7, which is also in disagreement with the values found by Andrews et al. (2010; i<sub>disk</sub> = 36° ± 10° or i<sub>disk</sub> = 144° ± 10°). Andrews et al. (2010) derived the disk
<figure><img src="image_6.png" /><figcaption>**Figure 4.** The SED of TWA 3A: photometric observations are shown in red, and the Spitzer/IRS spectrum is shown in orange. The Q-band observation (black empty circle) from Jayawardhana et al. (1999) was not used during the fitting. The best-fit model is shown in blue, along with the 100 highest- likelihood SED models (gray).</figcaption></figure>

inclination and position angle by fitting an elliptical Gaussian to marginally resolved sub-mm continuum observations, so it is only mildly surprising that their simplistic model deviates from the new values we derived using a more realistic dynamical gas model fit to higher quality observations. We constrained the central stellar mass to be $M_A = M_{Aa} + M_{Ab} = 0.534 \pm 0.010 M_\odot$ and the systemic velocity of the TWA 3A circum- binary disk to be $1.22 \pm 0.02 \text{ km s}^{-1}$ in the LSRK frame ($10.48 \pm 0.02 \text{ km s}^{-1}$ in the BARY frame).<fnref n="19" /> This is fully consistent with the radial velocity of the TWA 3A barycenter $v_{LSRK} = 0.91 \pm 0.40 \text{ km s}^{-1}$ ($\gamma_A = 10.17 \pm 0.40 \text{ km s}^{-1}$ on the CfA system; Kellogg et al. 2017).

## 3.2. SED Modeling and (sub)mm Spectral Index

In light of our new 226 GHz ALMA flux density measure- ment, we updated the SED of the TWA 3A system to learn about the properties of the circumbinary disk and its interior cavity. We sourced photometric fluxes from the SED compilation in Kellogg et al. (2017). We also incorporated the mid-IR spectrum from the IRS spectrograph on board the Spitzer Space Telescope (Andrews et al. 2010) with the contribution from the B component subtracted. The SED (shown in Figure 4) continuously decreases redward of ~20 µm, suggesting that the circumbinary disk is truncated at larger radii, highly settled, or both. The observed spectral slope between the SMA 880 µm and the 1.3 mm ALMA data is

$$ \alpha = \frac{d \log F_{\nu}}{d \log \nu} = 1.7 \pm 0.3, \qquad (3) $$

where the uncertainty is calculated by adopting a 10% calibration uncertainty for both points. That the spectral index is $\alpha \lesssim 2$ suggests the dust continuum emission is optically thick at these wavelengths, consistent with the continuum-subtrac- tion artifact described in Section 2. If the maximum grain size in the disk is between 500 µm and 1 cm, self-scattering from high albedo dust grains can reduce the (sub)mm-wave emission from an optically thick region and produce $\alpha < 2$ (Zhu et al. 2019).

Motivated by the near-infrared dip ($\lambda \sim 10$ µm) in the SED, we followed Andrews et al. (2010) and constructed a simple disk model with two zones: an inner cavity with a constant surface density from $r_{in}$ to $r_{cav}$, and a disk with surface density profile $\Sigma(r) \propto 1/r$ from $r_{cav}$ to $r_{out}$. We fixed $r_{in}$ to 0.2 au, following Kellogg et al.’s (2017) constraint on the inner binary semimajor axis of 0.19 au. We fixed $r_{out} = 8.5$ au and $i_{disk} = 49^\circ$, following our analysis of the $^{12}$CO J = 2–1 line emission. Our model had five free parameters: total disk dust mass $M_{dust}$, cavity radius $r_{cav}$, cavity depletion factor $\delta$, flaring parameter $\beta$, and scale height at 10 au $H_{10}$. The surface density within the cavity ($r_{in} < r < r_{cav}$) was set to a constant surface density $\Sigma_{cav} = \delta \Sigma_{disk}(r_{cav})$. The scale height as a function of radius is defined as

$$ H(r) = H_{10} \left( \frac{r}{10 \text{ au}} \right)^{\beta}. \qquad (4) $$

We followed the approach in Andrews et al. (2010), and combined the Aa and Ab components into a single stellar photosphere with an effective temperature of 3350 K and radius of $1 R_\odot$ (updated to the new Gaia distance), equivalent to $0.11 L_\odot$. We adopted the dust composition values used in Pinte et al. (2016): the grain size distribution followed $dn/da \propto a^{-3.5}$, where a is the grain size. The distribution ranged from $a_{min} = 0.01$ µm to $a_{max} = 1$ cm. We computed model SEDs using the MCFOST radiative transfer code (Pinte et al. 2006), assuming no interstellar extinction (McJunkin et al. 2014). Our grid of models spanned parameter ranges $M_{dust}$: [0.5, 1, 2, 4, 8, 16, 32, 64, 128] $\times 10^{-6} M_\odot$; $r_{cav}$: [0.2, 0.5, 0.75, 1, 1.5, 2] au; $\delta$: $1/[1, 10, 100, 1000]$; $\beta$: [1.025, 1.05, 1.075]; and $H_{10}$: [0.2, 0.3, 0.4, 0.5, 0.6] au.

For each model, we calculated a $\chi^2$ figure of merit using the observed photometry and IRS spectrum. Since we were only interested in disk properties, we only used photometric points $\lambda > 1$ µm in the following SED fit. We also excluded the Q- band (18.2 µm) observation from Jayawardhana et al. (1999) from the fit, since it is an outlier compared to the other photometric points. A set of consistently calibrated photometric and spectroscopic covariance matrices does not exist for the TWA 3A spectroscopic data set, and so we were unable to use per-data point flux uncertainties in our construction of a $\chi^2$ fit metric. Instead, we explored the consistency of the model grid with the SED data by the following procedure.

First, because adjacent pixels in spectroscopic fluxes are frequently correlated due to residual calibration errors, we subsampled the spectrum and only fit every third point. Then, we explored the consistency of the model grid with the SED data by assigning relative uncertainties of 5%, 10%, and 20% for each photometric and spectroscopic data point, and calculated the $\chi^2$ metric. We found that $r_{cav} = 0.20$ au was excluded at high significance (>99% probability) for all choices of uncertainty reweighting factors. The 20% reweight- ing for both photometric and spectroscopic data sets yielded the lowest reduced $\chi^2$ of 1.16. Based on the models with high figures of merit (see Figure 4), we conclude that the disk around TWA 3A has a dust mass of $1-8 \times 10^{-6} M_\odot$ and a disk cavity $r_{cav}$ of ≈0.5–0.75 au. The other model parameters were not well constrained.

## 3.3. Stellar Orbits

From the literature, we collected a diverse orbital data set, including radial velocity and astrometric measurements of all

<footnote n="19">¹⁹ In the direction of TWA 3, $v_{LSRK} = v_{BARY} - 9.26 \text{ km s}^{-1}$.</footnote>
three stars in the TWA 3 hierarchical triple (Reipurth & Zinnecker 1993; Webb et al. 1999; Weintraub et al. 2000; Brandeker et al. 2003; Correia et al. 2006; Janson et al. 2014; Anthonioz et al. 2015; Tokovinin et al. 2015; Kellogg et al. 2017; Knapp & Nanson 2018; Mason et al. 2018).<fnref n="20" /> Our goal was to extend the comprehensive analysis of Kellogg et al. (2017) by incorporating the new disk-based dynamical constraint on $M_A$ (see Section 3.1) into a joint hierarchical triple fit with the extant radial velocity and astrometric data sets and Gaia parallax.

We modeled these diverse data sets using the exoplanet software package (Foreman-Mackey et al. 2020). Briefly, exoplanet is designed to unify the necessary routines for orbital parameter inference within the PyMC3 (Salvatier et al. 2016) framework. Posterior gradients are provided through the Theano framework (Theano Development Team 2016), enabling the usage of powerful MCMC samplers like Hamiltonian Monte Carlo (HMC; Hoffman & Gelman 2011) to efficiently explore high dimensional spaces. To fit the TWA 3 data sets, we extended exoplanet to include functionality for astrometric orbits; these routines have been available in the main exoplanet package as of v0.2.0. We constructed the hierarchical model by nesting a Keplerian orbit for Aa–Ab (“inner”) inside of a wider orbit for A–B (“outer”). In the following analysis, we adopted orbital conventions where the argument of periastron $\omega$ is reported as the value of the “primary” star and $\Omega$ describes the position angle of the ascending node, which is the node where the secondary is receding from the observer. For the inner orbit, $\omega_{Aa}$ refers to the argument of periastron of TWA 3Aa and $\Omega_{\text{inner}}$ refers to the position angle of the TWA 3Aa–Ab ascending node. For the outer orbit, $\omega_A$ refers to the argument of periastron of TWA 3A (under the assumption that the Aa and Ab stars can be treated as a single star, A) and $\Omega_{\text{outer}}$ refers to the position angle of the TWA 3A–B ascending node.

Following standard radial velocity analysis, we included “jitter” and offset terms for each of the instruments. In keeping with Kellogg et al. (2017), we derived orbital parameters using the CfA radial velocity reference scale. Because there may be a small but unknown systematic radial velocity offset between the CfA and ALMA velocity scales, we did not use the systemic velocity of the TWA 3A circumbinary disk (at the epoch of the ALMA measurement) in the joint model. We applied uniform priors on the following quantities: $\log P_{\text{inner}}$, $\log P_{\text{outer}}$, $\cos i_{\text{inner}}$, $\cos i_{\text{outer}}$, and log jitter terms. We applied broad Gaussian priors on the sampled stellar masses of $M_{\text{Ab}}$ of $0.29 \pm 0.50 M_{\odot}$ and $M_{\text{B}}$ of $0.3 \pm 0.5 M_{\odot}$, loosely corresponding to the spectral types, and truncated to positive values only. We also applied broad Gaussian priors of $0.0 \pm 5.0 \text{ km s}^{-1}$ on the instrument offset terms.

Table 3 lists a full description of the inferred orbital parameters, where posterior means and standard deviations are provided for parameters whose posteriors are approximately Gaussian. Table 3 covers two scenarios: (1) the “primary” solution where $i_{\text{inner}} > 90^{\circ}$ and (2) the “alternate” solution where $i_{\text{inner}} < 90^{\circ}$. The alternate solution yields retrograde orbits between the inner and outer orbits. We show the phase-folded spectroscopic binary orbit (identical for both scenarios) in Figure 5, which is in good agreement with that found by Kellogg et al. (2017). The joint hierarchical fit delivered more precise posteriors for many outer orbital parameters, though

<table><caption>**Table 3**
Orbital Parameters</caption><thead><tr><th>Parameter<br>Sampled</th><th>Primary Solution</th><th>Alternate Solution</th></tr></thead><tbody><tr><td><i>P</i><sub>inner</sub> (days)</td><td>34.879 ± 0.001</td><td>34.879 ± 0.001</td></tr><tr><td><i>a</i><sub>inner</sub> (mas)</td><td>4.63 ± 0.04</td><td>4.63 ± 0.04</td></tr><tr><td><i>M</i><sub>Ab</sub> (<i>M</i><sub>☉</sub>)</td><td>0.24 ± 0.01</td><td>0.24 ± 0.01</td></tr><tr><td><i>e</i><sub>inner</sub></td><td>0.63 ± 0.01</td><td>0.63 ± 0.01</td></tr><tr><td><i>i</i><sub>inner</sub> (°)</td><td>131.5 ± 0.8</td><td>48.5 ± 0.8</td></tr><tr><td><i>ω</i><sub>Aa</sub><sup>a</sup> (°)</td><td>81 ± 1</td><td>81 ± 1</td></tr><tr><td><i>Ω</i><sub>inner</sub><sup>b</sup> (°)</td><td>104 ± 9</td><td>112 ± 9</td></tr><tr><td><i>T</i><sub>0,inner</sub> [JD-2,450,000]</td><td>2704.57 ± 0.07</td><td>2704.57 ± 0.07</td></tr><tr><td><i>P</i><sub>outer</sub> (yr)</td><td>548 ± 244</td><td>555 ± 249</td></tr><tr><td><i>M</i><sub>B</sub> (<i>M</i><sub>☉</sub>)</td><td>0.41 ± 0.28</td><td>0.40 ± 0.28</td></tr><tr><td><i>e</i><sub>outer</sub></td><td>0.3 ± 0.2</td><td>0.3 ± 0.2</td></tr><tr><td><i>i</i><sub>outer</sub> (°)</td><td>139 ± 13</td><td>139 ± 13</td></tr><tr><td><i>ω</i><sub>A</sub><sup>a</sup> (°)</td><td>...<sup>c</sup></td><td>...<sup>c</sup></td></tr><tr><td><i>Ω</i><sub>outer</sub><sup>b</sup> (°)</td><td>...<sup>c</sup></td><td>...<sup>c</sup></td></tr><tr><td><i>T</i><sub>0,outer</sub> [JD-2,450,000]</td><td>...<sup>c</sup></td><td>...<sup>c</sup></td></tr><tr><td><i>ω</i><sub>''</sub></td><td>27.31 ± 0.12</td><td>27.31 ± 0.12</td></tr><tr><td><i>σ</i><sub>ρ</sub><sup>''</sup></td><td>0.009 ± 0.004</td><td>0.009 ± 0.004</td></tr><tr><td><i>σ</i><sub>θ</sub> (°)</td><td>0.024 ± 0.008</td><td>0.024 ± 0.008</td></tr><tr><td><i>σ</i><sub>CFA</sub> (km s<sup>-1</sup>)</td><td>3.8 ± 0.3</td><td>3.8 ± 0.3</td></tr><tr><td><i>σ</i><sub>Keck</sub> (km s<sup>-1</sup>)</td><td>0.8 ± 0.1</td><td>0.8 ± 0.1</td></tr><tr><td><i>σ</i><sub>FEROS</sub> (km s<sup>-1</sup>)</td><td>3.4 ± 0.6</td><td>3.4 ± 0.6</td></tr><tr><td><i>σ</i><sub>du Pont</sub> (km s<sup>-1</sup>)</td><td>2.1 ± 0.4</td><td>2.1 ± 0.4</td></tr><tr><td>(Keck—CfA) (km s<sup>-1</sup>)</td><td>-1.3 ± 0.5</td><td>-1.3 ± 0.5</td></tr><tr><td>(FEROS—CfA) (km s<sup>-1</sup>)</td><td>1.3 ± 1.0</td><td>1.3 ± 1.0</td></tr><tr><td>(du Pont—CfA) (km s<sup>-1</sup>)</td><td>-0.2 ± 0.7</td><td>-0.2 ± 0.7</td></tr><tr><td>Derived</td><td></td><td></td></tr><tr><td><i>M</i><sub>Aa</sub> (<i>M</i><sub>☉</sub>)</td><td>0.29 ± 0.01</td><td>0.29 ± 0.01</td></tr><tr><td><i>M</i><sub>A</sub> (<i>M</i><sub>☉</sub>)</td><td>0.53 ± 0.01</td><td>0.53 ± 0.01</td></tr><tr><td><i>a</i><sub>inner</sub> (au)</td><td>0.170 ± 0.001</td><td>0.170 ± 0.001</td></tr><tr><td><i>a</i><sub>outer</sub> (au)</td><td>63 ± 18</td><td>64 ± 19</td></tr><tr><td><i>r</i><sub>p,outer</sub> (au)</td><td>45 ± 21</td><td>45 ± 21</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The argument of periastron of the primary. $\omega_{\text{secondary}} = \omega_{\text{primary}} + \pi$.

<sup>b</sup> The ascending node is identified as the point where the secondary body crosses the sky plane *receding* from the observer.

<sup>c</sup> Posterior is non-Gaussian and not accurately represented by a summary statistic; see Figure 6.

some of the degeneracies noted in Kellogg et al. (2017) still remain.

Three outer orbit parameters have bimodal posterior distributions: $\omega_A$, $\Omega_{\text{outer}}$, and $T_{0,\text{outer}}$. The posteriors for these outer orbit parameters are identical between the primary and alternate solutions, so the single corner plot in Figure 6 is valid for both scenarios. Representative outer orbits drawn from the posterior distribution are shown in Figure 7. Although the formal uncertainties and inferred jitter values of the Keck $v_B$ measurements are large ($0.59 \text{ km s}^{-1}$ and $0.82 \text{ km s}^{-1}$, respectively), the actual scatter of the four measured values is substantially smaller. The first two measurements (separated by 48 days in 2002/2003) differ by only $0.05 \text{ km s}^{-1}$. The second two measurements (separated by 1 yr in 2009/2010) differ by only $0.15 \text{ km s}^{-1}$. While circumstantial, this does raise the possibility that the uncertainties on the $v_B$ measurements are overestimated (potentially driven by the scatter in $v_{Aa}$ and $v_{Ab}$) and that the increase in $v_B$ over the Keck measurement baseline ($\approx 0.5 \text{ km s}^{-1}$ over 8 yr) may be significant. If true, a monotonically increasing $v_B$ clearly favors a single posterior mode, highlighted in blue in Figures 6 and 7 (the mode corresponding to decreasing $v_B$ is highlighted in orange). The

<footnote n="20">²⁰ Following Kellogg et al. (2017), we assigned a date of 1992.0216 to the observation by Reipurth & Zinnecker (1993).</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Phase-folded inner binary orbit and radial velocity residuals. In blue are 10 realizations of the inner orbit with the velocity trend from B removed. The small scatter demonstrates that the inner orbital parameters are tightly constrained by the data.</figcaption></figure>

MCMC samples and PyMC3 models corresponding to both scenarios are available online.<fnref n="21" />

We used these orbital posteriors and the inferred values of $i_{\text{disk}}$ and $\Omega_{\text{disk}}$ to calculate the mutual inclinations $\theta$ between the disk and the stellar orbits (e.g., Equation (1), Czekala et al. 2019) under different combinations of the degenerate orientations.<fnref n="22" /> As discussed in Czekala et al. (2019), the sensible imposition of spherically isotropic priors (i.e., the uniform priors on $\cos i_{\text{disk}}$ and $\cos i_{\text{inner}}$) results in effective mutual inclination priors of $p(\theta) \propto \sin(\theta)$. These isotropic priors have the consequence of strongly disfavoring coplanar architectures —simply because of the small phase space volume. To quantify the constraining power of the data only, we also report the posteriors re-weighted such that the effective prior is flat (i.e., similar to the marginal likelihood $p(\text{data}|\theta)$).

# 4. Discussion

## 4.1. Mutual Inclinations

As established in Section 3.3, the inner binary Aa–Ab orbit and circumbinary disk are very nearly coplanar with each other: across all scenarios listed in Table 4, $\theta_{\text{disk-inner}} \lesssim 6^\circ$ under a flat mutual inclination prior ($\theta_{\text{disk-inner}} \lesssim 13^\circ$ under $\sin\theta_{\text{disk-inner}}$ prior). Czekala et al. (2019) found that circumbinary planets, debris disks, and protoplanetary disks around short-period binaries ($P \lesssim 40$ days) all have low mutual inclinations. So, the

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Corner plot of the sampled outer orbital parameters that have bimodal posterior distributions; contours are $1\sigma$ and $2\sigma$ within each mode. This figure is identical for both the “primary” and “alternate” solutions. Samples are color-coded based on whether they deliver increasing (blue) or decreasing (orange) $v_B$ velocities over the Keck measurement baseline (see Figure 7). Note that the full posterior (the sum of the blue and orange contours) was sampled simultaneously; the samples have been bifurcated for plotting purposes only.</figcaption></figure>

low mutual inclination for the TWA 3A is consistent with expectations, given its 35-day orbital period.

That circumbinary disk mutual inclination trends with a binary period is likely a byproduct of the formation mechanism for tight binary stars, which requires formation at larger distances and migration to present-day configurations (Bate et al. 2002). The close binary fractions of T Tauri stars and field stars are similar (Kounkel et al. 2019), implying that this migration occurs quickly, before the class II T Tauri phase. It is unlikely that tertiary interactions (e.g., Fabrycky & Tremaine 2007) are responsible for the majority of tight binaries (Moe & Kratter 2018); rather, migration via a circumbinary disk appears to be the dominant pathway (Tokovinin & Moe 2020). TWA 3A represents both the lowest mass binary ($M_A = 0.53 \pm 0.01 M_\odot$) hosting a circumbinary disk and the longest period binary ($P_{\text{inner}} = 34.879 \pm 0.001$ days) before the population of mutual inclinations transitions from entirely coplanar systems to a broad dispersion of mutual inclinations (see Figure 14, Czekala et al. 2019). Though the range of binary periods over which this transition occurs is not yet well defined, recent observations have shed light on the dispersion of mutual inclinations at slightly longer binary periods of several months. The recent measurements of WW Cha by VLTI/GRAVITY have demonstrated that coplanar, truly low mutual inclination ($\theta < 8^\circ$) circumbinary disks can and still do exist around a longer period ($P = 207$ days), eccentric ($e = 0.45$) binaries (GRAVITY Collaboration et al. 2021). The measurement of WW Cha’s orbital properties is important because it enhances the contrast of the stark transition from coplanar systems to a broad distribution of mutual inclinations. As a rule, all planets and disks orbiting binaries with $P \lesssim 40$ days are coplanar, yet at binary periods of ~6–10 months, there already exists both a coplanar system (WW Cha)

<footnote n="21">²¹ https://zenodo.org/record/4568830#.YDvTb11Kida Czekala (2021) and https://github.com/iancze/TWA-3-orbit.</footnote>

<footnote n="22">²² Technically, there are yet two more degenerate scenarios where $i_{\text{disk}} < 90^\circ$ but $i_{\text{inner}} > 90^\circ$, or vice versa. Since the inferred inclinations of the inner binary and circumbinary disk are already so similar (modulo the degeneracy), we do not consider these scenarios.</footnote>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Thirty representative outer orbits drawn from the highlighted posterior modes in Figure 6. Left: the sky plane, centered on TWA 3A (Aa and Ab are represented by the black star), with the velocity field of the surrounding circumbinary disk. Right: the astrometric data and Kellogg et al. (2017) Keck measurements of $v_B$. Orbits in blue correspond to solutions that deliver increasing $v_B$ over the Keck observation baseline; orange orbits are decreasing.</figcaption></figure>

<table><caption>**Table 4**
Inferred Mutual Inclinations</caption><thead><tr><th>$i_{\text{disk}}$ (°)</th><th>$v_B$</th><th>$p(\theta)$</th><th>$\theta_{\text{disk-inner}}$ (°)</th><th>$\theta_{\text{inner-outer}}$ (°)</th><th>$\theta_{\text{disk-outer}}$ (°)</th></tr></thead><tbody><tr><td rowspan="2">&lt;90</td><td>↑</td><td>sin($\theta$)</td><td>&lt;8</td><td>99<sup>+17</sup><sub>-16</sub></td><td>99<sup>+18</sup><sub>-17</sub></td></tr><tr><td>↑</td><td>flat</td><td>&lt;3</td><td>100<sup>+17</sup><sub>-16</sub></td><td>100<sup>+18</sup><sub>-18</sub></td></tr><tr><td rowspan="2">&lt;90</td><td>↓</td><td>sin($\theta$)</td><td>&lt;8</td><td>140<sup>+10</sup><sub>-10</sub></td><td>143<sup>+7</sup><sub>-10</sub></td></tr><tr><td>↓</td><td>flat</td><td>&lt;3</td><td>142<sup>+10</sup><sub>-11</sub></td><td>145<sup>+7</sup><sub>-10</sub></td></tr><tr><td rowspan="2">&gt;90</td><td>↑</td><td>sin($\theta$)</td><td>&lt;13</td><td>22<sup>+14</sup><sub>-10</sub></td><td>20<sup>+19</sup><sub>-12</sub></td></tr><tr><td>↑</td><td>flat</td><td>&lt;6</td><td>&lt;20</td><td>&lt;15</td></tr><tr><td rowspan="2">&gt;90</td><td>↓</td><td>sin($\theta$)</td><td>&lt;13</td><td>76<sup>+14</sup><sub>-17</sub></td><td>82<sup>+13</sup><sub>-16</sub></td></tr><tr><td>↓</td><td>flat</td><td>&lt;6</td><td>75<sup>+15</sup><sub>-17</sub></td><td>82<sup>+14</sup><sub>-17</sub></td></tr></tbody></table>

**Note.** One $\sigma$ asymmetric error bars are reported for all unimodal distributions. Confidence upper limits of 68% are reported for one-sided distributions.

and a polar-oriented system (HD 98800B, $P=315$ days; Kennedy et al. 2019). Discovering new circumbinary disks with host binary periods near the “transition region” ($P=30–300$ days) and measuring their mutual inclinations will help map out the regimes where migration may deliver a coplanar system.

Substantial degeneracies remain in the relationship between the inner and outer stellar orbits of TWA 3. These degeneracies stem from two unknowns: (1) whether the circumbinary disk and inner binary are $i < 90^{\circ}$ or $i > 90^{\circ}$, and (2) which mode of the ($\omega_A$, $\Omega_{\text{outer}}$, $T_{0,\text{outer}}$) posterior is correct. The inferred mutual inclinations corresponding to the four permutations of these degeneracies are delineated in Table 4, under two different prior assumptions (sin($\theta$) or flat).<fnref n="23" /> These separate into a coplanar configuration ($\theta_{\text{inner-outer}} \approx 0^{\circ}$), two orthogonal

configurations ($\theta_{\text{inner-outer}} \approx 90^{\circ}$), and a retrograde configuration ($\theta_{\text{inner-outer}} > 90^{\circ}$). The field population of triples with projected outer separations <50 au exhibits a high degree of alignment between inner and outer orbital planes (average mutual inclination <20°; Tokovinin 2017). Considering this backdrop, we suspect that the TWA 3Aa–Ab and A–B orbital planes are nearly coplanar as well. We provide further hydrodynamical evidence for this scenario in the next subsection.

## 4.2. Disk Truncation

The time-dependent gravitational potential of a binary star will influence the radial extent of a protoplanetary disk, clearing an interior cavity in a circumbinary configuration and truncating the outer edge in a circumstellar configuration. TWA 3 is noteworthy because both types of disk truncation are present in the same system.

Several analytical and numerical works have derived the radius (usually conveyed as a ratio relative to the binary semimajor axis, $r/a$) to which an interior, coplanar, and eccentric binary like Aa–Ab ($e_{\text{inner}} = 0.63 \pm 0.01$) is predicted to clear the inner edge of the circumbinary disk (Artymowicz & Lubow 1994; Miranda & Lai 2015; Miranda et al. 2017; Thun et al. 2017; Hirsh et al. 2020). Using an edge definition where the density falls to 50% of its peak value, Artymowicz & Lubow (1994) found that $r/a \sim 2–3$, though they also noted that nonaxisymmetric waves at the rim of the disk make it difficult to define the edge location uniquely. Miranda & Lai (2015) studied circumbinary disk truncation across a range of mutual inclinations, finding that the truncation radius was smaller for more misaligned systems. For coplanar systems, their results agreed with Artymowicz & Lubow (1994). Using a suite of numerical simulations, Miranda et al. (2017) found that the truncation radius of coplanar circumbinary disks is $r/a \sim 1.7–2.6$ (using a 10% of peak density definition), though

<footnote n="23">²³ Since the value of $\theta_{\text{disk-inner}}$ is low, values of $\theta_{\text{inner-outer}}$ and $\theta_{\text{disk-outer}}$ are similar in all cases. We will refer to $\theta_{\text{inner-outer}}$ in what follows, but the points apply equally to $\theta_{\text{disk-outer}}$ as well.</footnote>
there is ambiguity in both the sharpness of the inner edge and peak location. Thun et al. (2017) used a 2D grid-based setup to derive scalings for various inner edge thresholds. For an $e \approx 0.6$ binary, they found that the inner edge scales as $r/a \sim 3.5$ for a 10% of peak definition, but found that the 50% location scales as $r/a \sim 4-6$. More recently, Hirsh et al. (2020) used a smoothed-particle hydrodynamics (SPH) setup to derive 50% thresholds for an $e \approx 0.6$ binary, finding values more in line with the studies by Artymowicz & Lubow (1994) and Miranda & Lai (2015): $r_{\text{inner}}/a_{\text{inner}} \approx 3.0-3.5$. Hirsh et al. (2020) noted that Thun et al.'s (2017) discrepancy with previous results could be attributable to their choice of inner polar grid boundary. For the TWA 3A circumbinary disk, the best-fit interior cavity radius of $r_{\text{cav}}/a_{\text{inner}} \approx 3$ inferred from SED modeling (Section 3.2) is squarely near the median of the aforementioned theoretical predictions, though solutions with $r_{\text{cav}}/a_{\text{inner}} \approx 4.5$ are still consistent with the SED data.

The inner edge of a circumbinary disk is also expected to be eccentric, due to resonant interactions between the binary and the disk. For an $e_{\text{inner}} = 0.63$ binary, Thun et al. (2017) found the inner edge would have $e_{\text{disk}} \approx 0.4$; Ragusa et al. (2020; using SPH simulations of more extreme mass ratio binaries) also found similar results: $e_{\text{disk}} \approx 0.05-0.35$. Muñoz & Lithwick (2020) used the linear theory of perturbed, pressure-supported disks to solve for the eccentricity profile and showed that the eccentric modes are concentrated to within $r/a < 2$ and drop off exponentially after $r/a \sim 10$ (corresponding to 1.7 au for TWA 3A). In simulations, the disk eccentricity is $e_{\text{disk}} < 0.05$ after $r/a = 10$ (Thun et al. 2017; Ragusa et al. 2020). Unfortunately, the scale of the inner rim of the circumbinary disk is below the resolution that can be meaningfully probed by the current ALMA observations. If the outer disk ($r \gtrsim 4$ au) were eccentric, it would be readily observable as a strong flux asymmetry between the redshifted and blueshifted sides of the disk (Czekala et al. 2015). Such an effect is not seen in these observations.

For equal mass binaries, nearly circular ($e \lesssim 0.1$) coplanar orbits are expected to truncate the outer edge of a circumstellar disk at radii $r_{\text{outer}}/a_{\text{outer}} \sim 0.32-0.38$ (Artymowicz & Lubow 1994; Miranda & Lai 2015). More eccentric binaries are more effective at truncating the disk: $r_{\text{outer}}/a_{\text{outer}} \sim 0.20-0.25$ for $e = 0.2-0.4$. Coplanar configurations are also the most effective at truncating the disk: as the mutual inclination $\theta_{\text{inner-outer}}$ increases, $r_{\text{outer}}$ also increases (about 20% larger for $\theta_{\text{inner-outer}} = 90^\circ$ and about 40% larger for $\theta_{\text{inner-outer}} = 135^\circ$; see Figure 4 of Miranda & Lai 2015). Using the 95% enclosed mass limit of $8.5 \pm 0.2$ au from the CO modeling and our orbital constraint of $a_{\text{outer}} = 63 \pm 18$ au means that $r_{\text{outer}}/a_{\text{outer}} \approx 0.10-0.19$ for TWA 3.

To further investigate the role of external disk truncation in the TWA 3 system, we performed 3D SPH simulations of the interactions between the TWA 3A binary, circumbinary disk, and TWA 3 B using the PHANTOM code (Price et al. 2018b). The SPH method is well suited for misaligned disk simulations, given that there is no preferred geometry and angular momentum is conserved to the accuracy of the time-stepping scheme (see, e.g., Price 2012). We used $10^6$ gas particles to model the circumbinary disk and set the initial inner and outer radii of the disk to $r_{\text{in}} = 2$ au and $r_{\text{out}} = 20$ au, respectively. The surface density initially followed a power-law profile ($\Sigma \propto R^{-1}$) and the temperature profile followed the power-law profile $T = 34$ K($r/10$ au)$^{-0.5}$, as in Price et al. (2018a). The disk's total mass was set to $0.01 M_\odot$, which allowed us to

neglect its self-gravity. Furthermore, we assumed that the disk is locally isothermal, where the sound speed follows a power-law $c_s \propto r^{-1/4}$ with $H/r = 0.05$ at $r = 10$ au. Finally, we adopted a mean Shakura-Sunyaev disk viscosity $\alpha_{SS} \approx 5 \times 10^{-3}$.

We included all three stars in the simulation (TWA 3Aa, Ab, and B, each represented by a sink particle), which interact with the gas via gravity and accretion (Bate et al. 1995). The accretion radii of Aa and Ab were set equal to 0.01 au, while the accretion radius of B was set equal to 10% of the Hill radius. These values ensure that the inner and outer regions of the disk are properly modeled, while keeping computational costs reasonable (see, e.g., Price et al. 2018a; Cuello et al. 2019; and Ménard et al. 2020).

The masses and orbits of stars Aa and Ab were initialized to the best-fit values listed in Table 3. The circumbinary disk around Aa and Ab was initially in the same plane as the Aa–Ab binary orbit. The outer companion B was set on a circular ($e_{\text{outer}} = 0$) orbit with semimajor axis $a_{\text{outer}}$ and mutual inclination with respect to the binary orbital plane of $\theta_{\text{inner-outer}}$. In principle, the stellar orbits are allowed to change as mass is accreted, but since an insignificant amount of material was accreted throughout the simulation, the change in stellar orbits was negligible. We ran six simulations corresponding to the coplanar, orthogonal, and retrograde orbits of $\theta_{\text{inner-outer}} = [0^\circ, 90^\circ, 130^\circ]$ and two different values of $a_{\text{outer}} = [45$ au, $80$ au]. We used $M_B = 0.4 M_\odot$ and $e_{\text{outer}} = 0$ for all simulations. A visualization of the results (convolved with an 8 au beam) is provided in Figure 8, showing that the outer edge of the disk is sensitive to choices of $\theta_{\text{inner-outer}}$ and $a_{\text{outer}}$. The smallest disks were produced in the coplanar simulations ($\theta_{\text{inner-outer}} = 0^\circ$). The simulation parameters corresponding to the $i_{\text{disk}} > 90^\circ$ and $v_B \uparrow$ scenario ($\theta_{\text{inner-outer}} = 0^\circ, a = 45$ au, $e=0$: $i0^\circ-a45$au) deliver an outer disk edge (10 au) that most closely matches the ALMA observations. This suggests that the true semimajor axis of the A–B binary lies closer to the lower range of its estimate ($a_{\text{outer}} = 63 \pm 18$ au), the true eccentricity of the A–B binary is significantly non-zero ($e_{\text{outer}} = 0.3 \pm 0.2$), or both.

### 4.3. Comparison to Other Systems

The HD 98800 multiple system, coincidentally in the same TW Hydra association as TWA 3 (and thus similarly aged at ~10 Myr), is an apt analog system to TWA 3. Because the HD 98800B circumbinary disk is in a hierarchical multiple system, it also experiences both interior and exterior disk truncation forces (from the Ba–Bb binary and from the wider A companion, respectively). Kennedy et al. (2019) spatially resolved the circumbinary disk with ALMA and convincingly demonstrated that it is in a circumpolar ($\theta_{\text{disk-inner}} \approx 90^\circ$) configuration. Franchini et al. (2019) demonstrated that the relatively small cavity size ($r_{\text{cav}}/a_{\text{inner}} \sim 1.5-2.5$) is a consequence of the reduced torques from an orthogonally oriented binary (see also Miranda & Lai 2015).

The exterior companion HD 98800A (which is itself a spectroscopic binary Aa–Ab, but here is effectively treated as a single star) orbits with $a_{\text{outer}} = 54$ au, $e = 0.52$, and $\theta_{\text{outer}} \approx 65^\circ$.<fnref n="24" /> The HD 98800B disk outer edge is $\lesssim 7$ au, suggesting that the eccentric outer companion is much more effective at outer-disk truncation ($r_{\text{outer}}/a_{\text{outer}} \sim 0.13$), even though the disk and the outer companion are substantially misaligned. The disk may

<footnote n="24">²⁴ Kennedy et al. (2019) were unable to break the $i_{\text{disk}} < 90^\circ$ or $i_{\text{disk}} > 90^\circ$ degeneracy, but the computed value of $\theta_{\text{outer}}$ is similar for both cases.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Gas surface density in SPH simulations of the TWA 3A binary, circumbinary disk, and exterior companion TWA 3 B, over a range of $\theta_{\text{inner-outer}}$ mutual inclinations (i) and semimajor axes $a_{\text{outer}}$ (a). The red dots at the center of the circumbinary disk represent the positions of the stars Aa and Ab (on top of each other) after five orbits of the outer companion B (not shown in frame). The images are convolved with an 8 au Gaussian beam (the white circle in the top left) for better comparison with the observations shown in Figure 1. The orbit of B that most closely reproduces the outer radius of the circumbinary disk corresponds to the coplanar one with the smaller semimajor axis ($\theta_{\text{inner-outer}} = 0^{\circ}$, $a_{\text{outer}} = 45$ au; i0°-a45). Each representation is shown from the perspective of an Earth observer.</figcaption></figure>

have survived as long as it has due to the combined effect of the inner-binary stopping the accretion onto the central source (Kuruwita & Federrath 2019) and the outer companion (s) stopping the viscous spreading of the disk (Ribas et al. 2018).

The HD 100453 multiple system is also a useful reference point for exterior disk truncation with a misaligned companion. The primary star HD 100453A (1.7 $M_{\odot}$ A9Ve; Dominik et al. 2003) is surrounded by a disk whose 1.4 mm continuum emission extends to ≈40 au and whose CO emission extends to ≈100 au (Wagner et al. 2018; van der Plas et al. 2019). Scattered light observations of the disk revealed spiral arms (Wagner et al. 2015) and narrow-lane shadows (Benisty et al. 2017). HD 100453B is an external 0.2 $M_{\odot}$ companion located at ≈100 au projected distance (Chen et al. 2006; Collins et al. 2009). Definitive conclusions about the influence of B on the disk are made difficult by the uncertainty in its orbit; however, a recent analysis by Gonzalez et al. (2020) supports a scenario where the disk and binary plane are substantially misaligned (60°), since coplanar orbits consistent with the astrometric data would be otherwise inconsistent with the disk morphology, including the spirals and observed velocity field. The spiral features and narrow-lane shadows seen in scattered light also suggest a complicated inner disk structure induced by an undetected, substellar companion interior to the disk (van der Plas et al. 2019; Nealon et al. 2020; Rosotti et al. 2020). The misalignments in this potentially triple system can be explained as being driven by the outer B, which drives the outer disk,

substellar companion, and inner disk to precess and occasionally undergo Kozai-Lidov oscillations (Nealon et al. 2020).

The LTT 1445ABC triple system, which hosts a transiting exoplanet (Winters et al. 2019), also bears mentioning in the context of TWA 3. LTT 1445ABC consists of three mid- to late M dwarfs in a hierarchical configuration: B-C forms the inner binary, and A is an outer tertiary and the most massive star in the system. The planet transits A, and the entire stellar system is coplanar. The TWA 3 system is something of a pre-main sequence counterpoint to LTT 1445ABC; in particular, the configuration of its circumbinary disk contrasts with the fact that in LTT 1445A, the planet transits the single star. However, if the TWA 3Aa-Ab stars were considered together, TWA 3A would also be the primary star in the system.

## 5. Conclusions

Our main conclusions from this study of the TWA 3 system are as follows.

1. We detected $^{12}\text{CO}$ $J=2-1$ and $^{13}\text{CO}$ $J=2-1$ emission from the TWA 3A circumbinary disk for the first time.

2. We forward modeled the $^{12}\text{CO}$ $J=2-1$ visibilities to derive an updated disk orientation ($\Omega_{\text{disk}} = 116.5 \pm 0.4$ and $i_{\text{disk}} = 48.8 \pm 0.7$) and infer the stellar mass enclosed by the disk $M_{\text{A}} = M_{\text{Aa}} + M_{\text{Ab}} = 0.534 \pm 0.010$.

3. We combined the disk dynamical constraints with extant radial velocity and astrometric measurements of TWA 3Aa, 3Ab, and 3B to infer individual stellar masses
0.29 ± 0.01 $M_\odot$, 0.24 ± 0.01 $M_\odot$, and 0.40 ± 0.28 $M_\odot$, respectively.

4. We drew constraints on the orbital architecture of the system and inferred that the plane of the inner Aa–Ab binary and its circumbinary disk are coplanar (misalignment <6° with 68% confidence). There are several degenerate solutions for the mutual inclination between the orbital planes of the inner (Aa–Ab) and outer (A–B) stellar orbits; however, SPH simulations lend support to the coplanar solution.

5. We found the inner and outer radii of the circumbinary disk (0.50–0.75 au and 8.5 ± 0.2 au, respectively) to be consistent with theoretical predictions of dynamical truncation from coplanar orbits.

I.C. would like to thank Rob De Rosa and Eric Nielsen for help with questions about orbital inference, Brian Mason for assistance with the Washington Double Star Catalog (Mason et al. 2001), Ryan Loomis for assistance with self-calibration, and Christophe Pinte for assistance with MCFOST. I.C. and A.R. would like to thank Elise Furlan for providing the calibrated IRS/Spitzer spectrum of TWA 3. I.C. was supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51405.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. This project has received funding

from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 210021. E.C. acknowledges NASA grants 80NSSC19K0506 and NNX15AD95G/NEXSS. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01545.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. This research has made use of NASA's Astrophysics Data System Bibliographic Services.

*Software:* CASA (v4.4; McMullin et al. 2007), DiskJockey (Czekala et al. 2015, 2019), RADMC-3D (Dullemond 2012), emcee (Foreman-Mackey et al. 2013), Astropy (Astropy Collaboration et al. 2013), PyMC3 (Salvatier et al. 2016), Theano (Theano Development Team 2016), MCFOST (Pinte et al. 2006), PHANTOM (Price et al. 2018b), SPLASH (Price 2007), uvplot (Tazzari 2017).

## Appendix

$^{13}$CO J = 2–1 Channel Maps

Channel maps for the $^{13}$CO J = 2–1 transition are provided in Figure 9.

<figure><img src="image_12.png" /><figcaption>Figure 9. Continuum-subtracted $^{13}$CO J = 2–1 data channel maps. The velocity scale is labeled in the LSRK frame.</figcaption></figure>
ORCID iDs

Ian Czekala https://orcid.org/0000-0002-1483-8811

Nicolás Cuello https://orcid.org/0000-0003-3713-8073

Eugene Chiang https://orcid.org/0000-0002-6246-2310

Enrique Macías https://orcid.org/0000-0003-1283-6262

Gaspard Duchêne https://orcid.org/0000-0002-5092-6464

Sean M. Andrews https://orcid.org/0000-0003-2253-2270

Catherine C. Espaillat https://orcid.org/0000-0001-9227-5949

References

Aly, H., Dehnen, W., Nixon, C., & King, A. 2015, MNRAS, 449, 65

Andrews, S. M. 2020, ARA&A, 58, 483

Andrews, S. M., Czekala, I., Wilner, D. J., et al. 2010, ApJ, 710, 462

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41

Anthonyoz, F., Ménard, F., Pinte, C., et al. 2015, A&A, 574, A41

Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58

Bate, M. R. 2019, MNRAS, 484, 2341

Bate, M. R., Bonnell, I. A., & Bromm, V. 2002, MNRAS, 336, 705

Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362

Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593

Benisty, M., Stolker, T., Pohl, A., et al. 2017, A&A, 597, A42

Brandeker, A., Jayawardhana, R., & Najita, J. 2003, AJ, 126, 2009

Brogan, C. L., Hunter, T. R., & Fomalont, E. B. 2018, arXiv:1805.05266

Chen, X. P., Henning, T., van Boekel, R., & Grady, C. A. 2006, A&A, 445, 331

Collins, K. A., Grady, C. A., Hamaguchi, K., et al. 2009, ApJ, 697, 557

Correia, S., Zinnecker, H., Ratzka, T., & Sterzik, M. F. 2006, A&A, 459, 909

Cuello, N., Dipierro, G., Mentiplay, D., et al. 2019, MNRAS, 483, 4114

Cuello, N., & Giuppone, C. A. 2019, A&A, 628, A119

Czekala, I. 2021, MCMC Samples for TWA 3 Orbital 957 Analysis, Zenodo, doi:10.5281/zenodo.4568830

Czekala, I., Andrews, S. M., Jensen, E. L. N., et al. 2015, ApJ, 806, 154

Czekala, I., Andrews, S. M., Torres, G., et al. 2016, ApJ, 818, 156

Czekala, I., Chiang, E., Andrews, S. M., et al. 2019, ApJ, 883, 22

Czekala, I., Jensen, E., & Huang, J. 2019, iancze/DiskJockey: Upgrades to Julia v1.5, Zenodo, doi:10.5281/zenodo.3235028

Dominik, C., Dullemond, C. P., Waters, L. B. F. M., & Walch, S. 2003, A&A, 398, 607

Dullemond, C. P. 2012, RADMC-3D: A Multi-purpose Radiative Transfer Tool v0.41, Astrophysics Source Code Library, ascl:1202.015

Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Foreman-Mackey, D., Luger, R., Czekala, I., et al. 2020, exoplanet-dev/exoplanet v0.3.2, Zenodo, doi:10.5281/zenodo.1998447

Foucart, F., & Lai, D. 2013, ApJ, 764, 106

Franchini, A., Lubow, S. H., & Martin, R. G. 2019, ApJL, 880, L18

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Gonzalez, J.-F., van der Plas, G., Pinte, C., et al. 2020, MNRAS, 499, 3837

GRAVITY Collaboration, Eupen, F., Labadie, L., et al. 2021, A&A, 648, A37

Herczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, 97

Hirsh, K., Price, D. J., Gonzalez, J.-F., Ubeira-Gabellini, M. G., & Ragusa, E. 2020, MNRAS, 498, 2936

Hoffman, M. D., & Gelman, A. 2011, arXiv:1111.4246

Janson, M., Bergfors, C., Brandner, W., et al. 2014, ApJS, 214, 17

Jayawardhana, R., Hartmann, L., Fazio, G., et al. 1999, ApJL, 521, L129

Kellogg, K., Prato, L., Torres, G., et al. 2017, ApJ, 844, 168

Kennedy, G. M., Matrà, L., Facchini, S., et al. 2019, NatAs, 3, 230

Knapp, W., & Nanson, J. 2018, JDSO, 14, 503

Kounkel, M., Covey, K., Moe, M., et al. 2019, AJ, 157, 196

Kuruwita, R. L., & Federrath, C. 2019, MNRAS, 486, 3647

Larson, R. B. 1969, MNRAS, 145, 271

Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2

Martin, R. G., & Lubow, S. H. 2017, ApJL, 835, L28

Mason, B. D., Hartkopf, W. I., Miles, K. N., et al. 2018, AJ, 155, 215

Mason, B. D., Wycoff, G. L., Hartkopf, W. I., Douglass, G. G., & Worley, C. E. 2001, AJ, 122, 3466

McJunkin, M., France, K., Schneider, P. C., et al. 2014, ApJ, 780, 150

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser., Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Ménard, F., Cuello, N., Ginski, C., et al. 2020, A&A, 639, L1

Miranda, R., & Lai, D. 2015, MNRAS, 452, 2396

Miranda, R., Muñoz, D. J., & Lai, D. 2017, MNRAS, 466, 1170

Moe, M., & Kratter, K. M. 2018, ApJ, 854, 44

Muñoz, D. J., & Lithwick, Y. 2020, ApJ, 905, 106

Nealon, R., Cuello, N., Gonzalez, J.-F., et al. 2020, MNRAS, 499, 3857

Offner, S. S. R., Dunham, M. M., Lee, K. I., Arce, H. G., & Fielding, D. B. 2016, ApJL, 827, L11

Offner, S. S. R., Kratter, K. M., Matzner, C. D., Krumholz, M. R., & Klein, R. I. 2010, ApJ, 725, 1485

Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, 459, 797

Price, D. J. 2007, PASA, 24, 159

Price, D. J. 2012, ICoPh, 231, 759

Price, D. J., Cuello, N., Pinte, C., et al. 2018a, MNRAS, 477, 1270

Price, D. J., Wurster, J., Tricco, T. S., et al. 2018b, PASA, 35, e031

Ragusa, E., Alexander, R., Calcino, J., Hirsh, K., & Price, D. J. 2020, MNRAS, 499, 3362

Reipurth, B., & Zinnecker, H. 1993, A&A, 278, 81

Remijan, A., Biggs, A., Cortes, P., et al. 2019, ALMA Cycle 7 Technical Handbook, Zenodo, doi:10.5281/zenodo.4511522

Ribas, Á, Macías, E., Espaillat, C. C., & Duchêne, G. 2018, ApJ, 865, 77

Rosotti, G. P., Benisty, M., Juhász, A., et al. 2020, MNRAS, 491, 1335

Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Comput. Sci., 2, e55

Sheehan, P. D., Wu, Y.-L., Eisner, J. A., & Tobin, J. J. 2019, ApJ, 874, 136

Tazzari, M. 2017, mtazzari/uvplot: v0.1.1, Zenodo, doi:10.5281/zenodo.1003113

Theano Development Team 2016, arXiv:1605.02688

Thun, D., Kley, W., & Picogna, G. 2017, A&A, 604, A102

Tofflemire, B. M., Mathieu, R. D., Herczeg, G. J., Akeson, R. L., & Ciardi, D. R. 2017, ApJL, 842, L12

Tofflemire, B. M., Mathieu, R. D., & Johns-Krull, C. M. 2019, AJ, 158, 245

Tokovinin, A. 2017, ApJ, 844, 103

Tokovinin, A., Mason, B. D., Hartkopf, W. I., Mendez, R. A., & Horch, E. P. 2015, AJ, 150, 50

Tokovinin, A., & Moe, M. 2020, MNRAS, 491, 5158

Tripathi, A., Andrews, S. M., Birnstiel, T., & Wilner, D. J. 2017, ApJ, 845, 44

van der Plas, G., Ménard, F., Gonzalez, J. F., et al. 2019, A&A, 624, A33

Wagner, K., Apai, D., Kasper, M., & Robberto, M. 2015, ApJL, 813, L2

Wagner, K., Dong, R., Sheehan, P., et al. 2018, ApJ, 854, 130

Weaver, E., Isella, A., & Boehler, Y. 2018, ApJ, 853, 113

Webb, R. A., Zuckerman, B., Platais, I., et al. 1999, ApJL, 512, L63

Weintraub, D. A., Saumon, D., Kastner, J. H., & Forveille, T. 2000, ApJ, 530, 867

Winters, J. G., Medina, A. A., Irwin, J. M., et al. 2019, AJ, 158, 152

Zanazzi, J. J., & Lai, D. 2018, MNRAS, 473, 603

Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJL, 877, L18