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
DDraft version March 3, 2021
Typeset using L A TEX twocolumn style in AASTeX63
A coplanar circumbinary protoplanetary disk in the TWA 3 triple M dwarf system
Ian Czekala,
1, 2, 3, 4, 5, ∗ ´Alvaro Ribas, Nicol´as Cuello, Eugene Chiang,
5, 8
Enrique Mac´ıas,
9, 6
Gaspard Duchˆene,
5, 7
Sean M. Andrews, and Catherine C. Espaillat Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802,USA 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´ordova 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´ordova 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
ABSTRACTWe present sensitive ALMA observations of TWA 3, a nearby, young ( ∼
10 Myr) hierarchical systemcomposed of three pre-main sequence M3–M4.5 stars. For the first time, we detected CO and CO J =2-1 emission from the circumbinary protoplanetary disk around TWA 3A. We jointly fit theprotoplanetary disk velocity field, stellar astrometric positions, and stellar radial velocities to infer thearchitecture of the system. The Aa and Ab stars (0 . ± . M (cid:12) and 0 . ± . M (cid:12) , respectively)comprising the tight ( P = 35 days) eccentric ( e = 0 . ± .
01) spectroscopic binary are coplanarwith their circumbinary disk (misalignment < ◦ with 68% confidence), similar to other short-periodbinary systems. From models of the spectral energy distribution, we found the inner radius of thecircumbinary disk ( r inner = 0 . − .
75 au) to be consistent with theoretical predictions of dynamicaltruncation r cav /a inner ≈
3. The outer orbit of the tertiary star B (0 . ± . M (cid:12) , a ∼ ±
18 au, e = 0 . ± .
2) is not as well constrained as the inner orbit, however, orbits coplanar with the A systemare still preferred (misalignment < ◦ ). To better understand the influence of the B orbit on theTWA 3A circumbinary disk, we performed SPH simulations of the system and found that the outeredge of the gas disk ( r outer = 8 . ± . Keywords: protoplanetary disks – stars: pre-main sequence – orbits – classical T Tauri stars – Trinarystars – M dwarf stars INTRODUCTIONThe distribution of pre-main sequence multiple systemarchitectures informs our understanding of the mecha-nisms that govern star and planet formation. Recently,Czekala et al. (2019) found that the degree of alignmentbetween the disk and its host binary (the mutual inclina-tion, θ ) is a strong function of orbital period. Circumbi-nary disks around short-period binaries ( P <
40 days)
Corresponding author: Ian [email protected] ∗ NASA Hubble Fellowship Program Sagan Fellow are preferentially coplanar, while disks around longer pe-riod binaries exhibit a wide range of mutual inclinations,including polar configurations ( θ ≈ ◦ ). These trendsmight be manifestations of the same physical mecha-nisms that produce close binaries. Because it is diffi-cult to directly fragment stars on scales < a r X i v : . [ a s t r o - ph . E P ] M a r Typically, a nearly coplanar disk surrounding a loweccentricity binary will precess around the binary an-gular momentum vector as dissipative forces damp theangular momentum vectors of the disk and the binaryinto alignment (Foucart & Lai 2013). If the binary issufficiently eccentric, however, then the disk can pre-cess around the eccentricity vector (Aly et al. 2015;Martin & Lubow 2017; Zanazzi & Lai 2018; Cuello &Giuppone 2019) and access polar mutual inclinations.Indeed, highly misaligned disks are preferentially foundaround highly eccentric ( e > .
7) binaries (Kennedyet al. 2019; Czekala et al. 2019). There are, however,several disks around eccentric short period binaries thatare coplanar, such as AK Sco, DQ Tau, and UZ Tau E(Czekala et al. 2015, 2016, 2019). Our understanding ofhow these period and eccentricity trends interrelate islimited by the sample size of circumbinary disk systemswith well-measured architectures.The pre-main sequence system TWA 3 represents anopportunity to expand this sample to aid in the in-terpretation of binary formation and evolution mecha-nisms. TWA 3 consists of three young (10 ± ∼ . M (cid:12) (Herczeg & Hillenbrand 2014; Tofflemireet al. 2019). The inner Aa–Ab binary has an orbitalperiod of P = 34 . ± . e = 0 . ± . (cid:36) = 27 . ± .
12 mas (including a 0 .
02 mas system-atic term, Lindegren et al. 2018) corresponding to a dis-tance of 36 . ± .
16 pc (Gaia Collaboration et al. 2018;Bailer-Jones et al. 2018). Tofflemire et al. (2019) notedthat the A and B components suffer from significantexcess astrometric noise—possibly due to photometricvariability—but the parallax distances for each sourceare consistent with each other and the A–B orbit arcsin Kellogg et al. (2017). Time-series photometry andemission line spectroscopy revealed that accretion fromthe circumbinary disk to the inner binary is phased withperiastron, and that material is preferentially accretedonto the primary star, TWA 3Aa (Tofflemire et al. 2017,2019). The gradual movement of the outer triple com-panion (projected separation 1 . (cid:48)(cid:48)
55 or 57 au; Tokovininet al. 2015) over a ∼
20 yr baseline suggests an orbitalperiod of ∼ −
800 yr (Kellogg et al. 2017).Andrews et al. (2010) used the Submillimeter Array(SMA) to localize the submillimeter emission in theTWA 3 system to the A binary, measuring a flux densityof 75 mJy at 340 GHz. As demonstrated by fits to thedeprojected and azimuthally averaged baselines, the cir-cumbinary disk itself was only marginally resolved (1 . (cid:48)(cid:48) × . (cid:48)(cid:48)
74 beam), but found to have a radius of ∼
20 au.Andrews et al. (2010) did not detect CO J = 3 − The eccentricity vector is drawn from binary apoapse to periapse. emission to an upper limit of 0.6 Jy beam − integratedover a 0 . − channel. Based upon the fit of an ellip-tical Gaussian to the visibilities, Andrews et al. (2010)derived a disk inclination of i disk = 36 ◦ ± ◦ (relative tothe sky plane) and disk orientation of Ω disk = 169 ◦ ± ◦ (the position angle of the ascending node measured eastof north). More recent orbital solutions suggested that the in-ner binary orbit, circumbinary disk, and outer ter-tiary orbit may be misaligned. Kellogg et al. (2017)combined an astrometric observation of the inner bi-nary (Anthonioz et al. 2015) with their double-lined ra-dial velocity solution to constrain the position angle ofthe ascending node Ω inner ∈ [93 ◦ , ◦ ] and inclination i inner ∈ [32 ◦ , ◦ ] or i inner ∈ [118 ◦ , ◦ ]. These orbitalparameters, together with the disk parameters reportedin Andrews et al. (2010), suggested that the planes ofthe spectroscopic binary and the circumbinary disk weremisaligned by at least θ ∼ ◦ (Kellogg et al. 2017).We acquired ALMA observations of the TWA 3A cir-cumbinary disk to better understand its size and orien-tation relative the stellar orbits. In § § § § DATAWe obtained deep Atacama Large Millimeter Array(ALMA) observations of the TWA 3 system in 2018.We used a correlator setup that assigned two 2 GHzwide spectral windows to the dust continuum (centeredon 220 GHz and 232 GHz) and two spectral windows at122 kHz (0.16 km s − ) resolution to target the CO and CO J = 2 − The lack of a gas detection meant degenerate “flipped” solutionswere also valid: i disk = 144 ◦ ± ◦ and Ω disk = 349 ◦ ± ◦ . Table 1.
Image propertiesbeam dimensions, P.A. RMS [mJy beam − ]226 GHz cont. 0 . (cid:48)(cid:48) × . (cid:48)(cid:48)
19, -64 ◦ CO J =2 − . (cid:48)(cid:48) × . (cid:48)(cid:48)
19, -64 ◦ CO J =2 − . (cid:48)(cid:48) × . (cid:48)(cid:48)
20, -63 ◦ Note —The RMS noise levels for the spectral line cubes correspondto the values per 0.8 km s − channel. All images were synthesizedwith robust=0.5 . spent 47.4 minutes on-source, for a total on-source timeof 1 hr 34.8 m. The mean precipitable water vapor foreach observation was 1.9mm and 0.4mm, respectively.We began our data reduction with thepipeline-calibrated measurement set provided byALMA/NAASC staff. We used the CASA 5.4 (Mc-Mullin et al. 2007) facility software and followed com-mon calibration and imaging procedures (e.g., theDSHARP reduction scripts, Andrews et al. 2018) , thepertinent details of which we now describe. To assessthe quality of each execution block we first reduced eachobservation individually. We excised the channels withline-emission to create a continuum-only measurementset with a total bandwidth of 4.6 GHz. We performedan initial round of continuum imaging using the CASA tclean task with robust=0.5 , an image size of 512 × . (cid:48)(cid:48)
015 pixel size, deconvolver="multiscale" , scales=[0, 15, 30, 45, 75] pixels, a threshold of0.6 mJy, and an elliptical mask with position angle110 ◦ , semi-major axis of 0 . (cid:48)(cid:48)
45 and semi-minor axis of0 . (cid:48)(cid:48) imfit and uvmodelfit tasks. Wefound excellent astrometric agreement between exe-cutions, with the emission centroid located at ICRS11:10:27.731 -37.31.51.84, coincident with the Gaia po-sition of TWA 3A to within 0 . (cid:48)(cid:48)
05. We found adequateagreement between the total continuum fluxes (31.9 mJyand 36.1 mJy, respectively), only slightly more differ-ent than the expected 10% amplitude calibration uncer-tainty. We proceeded to self-calibrate the combined measure-ment set through a series of applications of tclean tothe continuum visibilities using threshold depths { } mJy beam − interleaved with applica-tions of the gaincal and applycal CASA tasks using https://bulk.cv.nrao.edu/almadata/lp/DSHARP/ 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 spectral-window dependent solves with intervals {
60 s,30 s, and 18 s } . We then cleaned to a final depth of0 .
02 mJy beam − and performed one round of phase andamplitude gain solutions over the scan length (8 min-utes, solint=‘inf’ ). We monitored the peak flux, to-tal flux, and RMS of the images throughout the process(Brogan et al. 2018) and found the peak continuum S/Nimproved to 2130 from an initial value of 260.We then used the applycal task to apply the self-calibration solutions to channels containing the spectralline observations (without propagating flags for failedsolutions, applymode=calonly ). We estimated the con-tinuum from nearby line-free channels and subtracted itfrom the spectral line observations using the uvcontsub CASA task. The line channels were imaged usingthe tclean task with the auto-multithresh mask-ing algorithm with cell="0.015arcsec" , gain=0.1 , deconvolver="multiscale" , scales=[0, 10, 30,100, 200, 300] pixels, robust=0.5 , and deconvolvedto a depth of threshold="0.1mJy" . The beam dimen-sions (FWHM) and image-plane RMS are described inTable 1.We summed the pixels within the continuum CLEANmask to measure a total continuum flux of 36.7 mJy.The dust continuum emission is compact, with nearlyall of the flux contained within the central beam. Wefit elliptical Gaussians to the continuum emission withthe uvmodelfit and imfit tasks and derived FWHMdimensions of (0 . (cid:48)(cid:48) × . (cid:48)(cid:48)
12) and (0 . (cid:48)(cid:48) × . (cid:48)(cid:48) uvplot pack-age (Tazzari 2017), shown in Figure 2. We used de-projection values of i disk = 49 ◦ and Ω disk = 116 . ◦ ,which were derived from our modeling effort of the CO J = 2 − § ≈ λ implies that the image planeGaussian has FWHM ≈ . (cid:48)(cid:48)
15, or 5.5 au.In preliminary line imaging, we identified CO emis-sion from [-9, 10] km s − (LSRK). To save computationalcomplexity, we used the task mstransform to averagethe channels to 0.8 km s − width. The CO emissionwas strongly detected across this velocity range (peakchannel S/N= 30, see Figure 3, top panel); the CO Though part of a standard self-calibration workflow, we noticedthat the S/N of the line channels did not measurably improve af-ter applying the self-calibration solutions. We believe that this isbecause the fine resolution channels are thermal-noise dominatedand not limited by residual phase errors. − . . . α cos δ [ ] − . . . ∆ δ [ ] dust dust + CO . . . . CO .
00 0 .
05 0 . CO − .
02 0 .
00 0 . − Jy beam − · km s − Jy beam − · km s − Jy beam − · km s − 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 § r e a l[ m J y ] λ ] − i m ag [ m J y ] Figure 2.
The continuum visibilities deprojected from the diskinclination. That the flux drops with increasing uv distance in-dicates that the continuum emission is spatially resolved. TheGaussian profile with FWHM ≈ λ implies an image planemorphology with FWHM ≈ . (cid:48)(cid:48)
15, or 5.5 au. emission was detected at lower but still significant levels(peak channel S/N= 6, see Appendix A). We used the immoments
CASA task to sum all flux in each channelacross the velocity dimension, producing the momentmaps in the third and fourth panels of Figure 1. Weused the imstat
CASA task to sum all flux within theCLEAN mask across the velocity and and spatial dimen-sions, yielding integrated fluxes of 526 mJy km s − and86 mJy km s − for CO and CO, respectively. We noticed a central cavity in the moment maps ofthe continuum-subtracted line emission, correspondingto the location of peak dust continuum. The depressionis most apparent in the CO emission but is also visiblein the CO emission (Figure 1, third and fourth pan-els). To investigate whether this may be a continuumsubtraction artefact, we produced channel maps and amoment map for the non-continuum subtracted COspectral channels (Figure 1, second panel). This mo-ment map does not exhibit a central cavity, suggestingthat the feature seen in CO and CO is indeed a con-tinuum subtraction artefact. In protoplanetary disks,such an artefact can arise when gas emission on the near-side of the disk is optically thick and absorbs continuumemission originating from the disk midplane. When thecontinuum emission (estimated from channels offset invelocity from the line emission) is subtracted, most ifnot all of the line flux is also (erroneously) subtracted(for a full description of the effect, see Weaver et al.2018). For this artefact to be significant, the continuumemission needs to have a brightness temperature compa-rable to the line emission, suggesting that the continuumemission is also optically thick, or nearly so.No noticeable continuum or gas emission is present inthe system beyond the disk surrounding TWA 3A. Weused an aperture approximately three times the areaof the beam to extract continuum photometry at thelocation of B (see § ± µ Jy). ANALYSIS3.1.
Dynamical gas analysis
Following the approach in Czekala et al. (2015),we constructed a forward model of the continuum-subtracted CO J = 2 − RADMC-3D (Dullemond 2012) to produceimage cubes. These cubes are Fourier-transformed andsampled at the spatial frequencies corresponding to thebaselines of the ALMA measurement set to compute thelikelihood of the dataset. The posterior of the modelparameters, defined by the data likelihood and any ad-ditional prior probability distributions, is explored usingMarkov Chain Monte Carlo (MCMC).We initially explored models with self-similar surfacedensity profiles Σ( r ) of the form described in Czekalaet al. (2015) and references therein (see also Sheehanet al. 2019). We found that the gradual exponential ta-per of this profile at large radii was not well-matched tothe TWA 3A CO J = 2 − I ( r ) using theNuker function (e.g., Tripathi et al. 2017), we adaptedthis to a surface density profile asΣ( r ) = Σ c (cid:18) rr c (cid:19) − γ (cid:20) (cid:18) rr c (cid:19) α (cid:21) ( γ − β ) /α . (1)We followed Tripathi et al. (2017) and chose to sample inlog α and imposed priors on shape parameters log α , β , and γ . We restricted log α and β with uniformprobability to the ranges of 0 ≤ log α ≤ ≤ β ≤
10, respectively, and imposed a tapered prior on γ of the form p ( γ ) = 11 + e − γ +3) −
11 + e − γ − . (2)As discussed in Tripathi et al. (2017), these priors arepractically motivated to allow a broad range of surfacedensity profiles, including those with interior cavitiesand sharp outer edges, while restricting the parameterspace sufficiently to avoid pathological surface densityprofiles ill-suited to protoplanetary disks. Because theNuker profile increases the dimensionality of the poste-rior by adding two new parameters, we kept computa-tional demands tractable by zeroing out the phase-centeroffsets δ α , δ δ typically used with the standard model.We sampled the posterior distribution using the DiskJockey package (Czekala et al. 2015) and assessedconvergence both visually and by applying the Gelman-Rubin convergence diagnostic to independent chain en-sembles. The resulting marginal posteriors on the diskmodel parameters are listed in Table 2. Most marginaldistributions are well-described by Gaussians. The pos-teriors on the temperature profile exponent q and Nukershape parameters log α and β ran up against the rangeof their prior bounds, so 68% upper or lower confidence Using the
DiskJockey.jl package (Czekala et al. 2015), https://github.com/iancze/DiskJockey
Table 2.
Inferred Disk ParametersParameter Value M A [ M (cid:12) ] 0 . ± . r c [au] 6 . ± . T [K] 38 ± γ − . ± . q ≤ . α ≥ . β ≥ . M disk log [ M (cid:12) ] − . ± . ξ [km s − ] 0 . ± . i disk [ ◦ ] 48 . ± . a Ω disk [ ◦ ] 116 . ± . v r [km s − ] 1 . ± . b a Ambiguity with i disk = 131 . ± . ◦ b LSRK reference frame.
Note —Using CO. The 1D marginal pos-teriors are well-described by a Gaussian,so we report symmetric error bars here(statistical uncertainties only). These pa-rameters were inferred using a distance of d = 36 .
62 pc. intervals are quoted for these parameters instead. Fig-ure 3 shows a realization of the model and residual visi-bilities drawn from the posterior distribution, imaged inthe same way as the data. The residual channel maps arebroadly consistent with residual thermal noise, demon-strating that the synthesized model is a good fit to thedata. We define the outer edge of the disk using theradius that contains 95% of the mass: r outer = 8 . ± . q ≈ CO J = 2 − CO J = 2 − -8.4 km s − -8.4 km s − -8.4 km s − -7.6-7.6-7.6 -6.8-6.8-6.8 -6.0-6.0-6.0 -5.2-5.2-5.2 -4.4-4.4-4.4 -3.6-3.6-3.6 -2.8-2.8-2.8-2.0-2.0-2.0 -1.2-1.2-1.2 -0.4-0.4-0.4 0.40.40.4 1.21.21.2 2.02.02.0 2.82.82.8 3.63.63.64.44.4 − . . . α cos δ [ ] − . . . ∆ δ [ ] − − −
10 0 10 20 30 − − −
10 0 10 20 30 × RMSmJy beam − d a t a m o d e l r e s i du a l Figure 3.
Continuum-subtracted CO J = 2 − and gas model simultaneously fit to the dust continuumand CO J = 2 − M A , i disk , and Ω disk ) via MCMC.We constrained the disk position angle to be Ω disk =116 . ◦ ± . ◦ , which is significantly different from thevalue found by Andrews et al. (2010) (Ω disk = 144 ◦ ± ◦ ). We constrained the disk inclination to be either i disk = 48 . ◦ ± . ◦ or i disk = 131 . ◦ ± . ◦ , which isalso in disagreement with the values found by Andrewset al. (2010) ( i disk = 36 ◦ ± ◦ or i disk = 144 ◦ ± ◦ ). An-drews et al. (2010) derived the disk inclination and posi-tion angle by fitting an elliptical Gaussian to marginallyresolved sub-mm continuum observations, so it is onlymildly surprising that their simplistic model deviatesfrom the new values we derived using a more realistic dy-namical gas model fit to higher quality observations. Weconstrained the central stellar mass to be M A = M Aa + M Ab = 0 . ± . M (cid:12) and the systemic velocity ofthe TWA 3A circumbinary disk to be 1 . ± .
02 km s − in the LSRK frame (10 . ± .
02 km s − in the BARYframe). This is fully consistent with the radial velocityof the TWA 3A barycenter v LSRK = 0 . ± .
40 km s − ( γ A = 10 . ± .
40 km s − on the CfA system, Kellogget al. 2017).3.2. SED modeling and (sub)mm spectral index
In light of our new 226 GHz ALMA flux density mea-surement, we updated the spectral energy distribution(SED) of the TWA 3A system to learn about the prop-erties of the circumbinary disk and its interior cavity.We sourced photometric fluxes from the SED compila-tion in Kellogg et al. (2017). We also incorporated themid-IR spectrum from the IRS spectrograph onboardthe
Spitzer Space Telescope (Andrews et al. 2010) withthe contribution from the B component subtracted. TheSED (shown in Figure 4) continuously decreases red-ward of ∼ µ m, suggesting that the circumbinary diskis truncated at larger radii, highly settled, or both. Theobserved spectral slope between the SMA 880 µ m andthe 1.3 mm ALMA data is: α = d log F ν d log ν = 1 . ± . , (3)where the uncertainty is calculated by adopting a 10 %calibration uncertainty for both points. That the spec-tral index is α (cid:46) In the direction of TWA 3, v LSRK = v BARY − .
26 km s − . λ [ µ m]10 -13 -12 -11 -10 -9 ν F ν [ e r g c m − s − ] Figure 4.
The SED of TWA 3A: Photometric observations areshown in red, and the
Spitzer /IRS spectrum 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 shownin blue along with the 100 highest likelihood SED models (grey). the continuum subtraction artifact described in §
2. Ifthe maximum grain size in the disk is between 500 µ mand 1 cm, self-scattering from high albedo dust grainscan reduce the (sub)mm-wave emission from an opti-cally thick region and produce α < λ ∼ µ m) in theSED, we followed Andrews et al. (2010) and constructeda simple disk model with two zones: an inner cavitywith a constant surface density from r in to r cav , and adisk with surface density profile Σ( r ) ∝ /r from r cav to r out . We fixed r in to 0.2 au, following Kellogg et al.(2017)’s constraint on the inner binary semi-major axisof 0 .
19 au. We fixed r out = 8 . i disk = 49 ◦ fol-lowing our analysis of the CO J = 2 − M dust , cavity radius r cav , cavity depletion factor δ , flar-ing parameter β , and scale height at 10 au H . Thesurface density within the cavity ( r in < r < r cav ) wasset to a constant surface density Σ cav = δ Σ disk ( r cav ).The scale height as a function of radius is defined as H ( r ) = H (cid:16) r
10 au (cid:17) β . (4)We followed the approach in Andrews et al. (2010) andcombined the Aa and Ab components into a single stel-lar photosphere with an effective temperature of 3350 Kand radius of 1 R (cid:12) (updated to the new Gaia distance),equivalent to 0.11 L (cid:12) . We adopted the dust compositionvalues used in Pinte et al. (2016): the grain size distri-bution followed d n/ d a ∝ a − . , where a is the grainsize. The distribution ranged from a min = 0 . µ m to a max =1 cm. We computed model SEDs using the MC-FOST radiative transfer code (Pinte et al. 2006) assum-ing 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] × − M (cid:12) ; r cav : [0.2,0.5, 0.75, 1, 1.5, 2] au; δ : 1/[1, 10, 100, 1000]; β : [1.025,1.05, 1.075]; and H : [0.2, 0.3, 0.4, 0.5, 0.6] au.For each model, we calculated a χ figure of merit us-ing the observed photometry and IRS spectrum. Sincewe were only interested in disk properties, we onlyused photometric points λ > µ m in the followingSED fit. We also excluded the Q-band (18.2 µ m) ob-servation from Jayawardhana et al. (1999) from the fitsince it is an outlier compared to the other photometricpoints. A set of consistently calibrated photometric andspectroscopic covariance matrices does not exist for theTWA 3A spectroscopic dataset, and so we were unableto use per-datapoint flux uncertainties in our construc-tion of a χ fit metric. Instead, we explored the con-sistency of the model grid with the SED data by thefollowing procedure.First, because adjacent pixels in spectroscopic fluxesare frequently correlated due to residual calibration er-rors, we subsampled the spectrum and only fit everythird point. Then, we explored the consistency of themodel grid with the SED data by assigning relative un-certainties of 5%, 10%, and 20% for each photometricand spectroscopic datapoint and calculated the χ met-ric. We found that r cav = 0 .
20 au was excluded at highsignificance ( >
99% probability) for all choices of un-certainty reweighting factors. The 20% reweighting forboth photometric and spectroscopic datasets yielded thelowest reduced χ ν of 1.16. Based on the models withhigh figures of merit (see Figure 4) we conclude that thedisk around TWA 3A has a dust mass of 1 − × − M (cid:12) and a disk cavity r cav of ≈ Stellar orbits
From the literature, we collected a diverse orbitaldataset including radial velocity and astrometric mea-surements of all three stars in the TWA 3 hierarchicaltriple (Reipurth & Zinnecker 1993; Webb et al. 1999;Weintraub et al. 2000; Brandeker et al. 2003; Correiaet al. 2006; Janson et al. 2014; Tokovinin et al. 2015;Anthonioz et al. 2015; Kellogg et al. 2017; Knapp & Nan-son 2018; Mason et al. 2018). Our goal was to extendthe comprehensive analysis of Kellogg et al. (2017) byincorporating the new disk-based dynamical constrainton M A (see § Gaia parallax.We modelled these diverse datasets using the exoplanet software package (Foreman-Mackey et al.2020). Briefly, exoplanet is designed to unify routinesneeded for orbital parameter inference within the
PyMC3 (Salvatier et al. 2016) framework. Posterior gradients Following Kellogg et al. (2017), we assigned a date of 1992.0216to the observation by Reipurth & Zinnecker (1993). are provided through the
Theano framework (TheanoDevelopment Team 2016), enabling the usage of pow-erful MCMC samplers like Hamiltonian Monte Carlo(HMC; Hoffman & Gelman 2011) to efficiently explorehigh dimensional spaces. To fit the TWA 3 datasets,we extended exoplanet to include functionality for as-trometric orbits; these routines have been available inthe main exoplanet package as of v0.2.0. We con-structed the hierarchical model by nesting a Keplerianorbit for Aa–Ab (“inner”) inside of a wider orbit forA–B (“outer”). In the following analysis, we adoptedorbital conventions where the argument of periastron ω is reported as the value of the “primary” star and Ω de-scribes the position angle of the ascending node, whichis the node where the secondary is receding from the ob-server. For the inner orbit, ω Aa refers to the argumentof periastron of TWA 3Aa and Ω inner refers to the po-sition angle of the TWA 3Aa–Ab ascending node. Forthe outer orbit, ω A refers to the argument of periastronof TWA 3A (under the assumption that the Aa and Abstars can be treated as a single star, A) and Ω outer refersto the position angle of the TWA 3A–B ascending node.Following standard radial velocity analysis, we in-cluded “jitter” and offset terms for each of the instru-ments. In keeping with Kellogg et al. (2017), we de-rived orbital parameters using the CfA radial veloc-ity reference scale. Because there may be a small butunknown systematic radial velocity offset between theCfA and ALMA velocity scales, we did not use the sys-temic velocity of the TWA 3A circumbinary disk (at theepoch of the ALMA measurement) in the joint model.We applied uniform priors on the following quantities:log P inner , log P outer , cos i inner , cos i outer , and log jitterterms. We applied broad Gaussian priors on the sam-pled stellar masses of M Ab of 0 . ± . M (cid:12) and M B of 0 . ± . M (cid:12) , loosely corresponding to the spectraltypes, and truncated to positive values only. We alsoapplied broad Gaussian priors of 0 . ± . − on theinstrument offset terms.Table 3 lists a full description of the inferred orbitalparameters, where posterior means and standard devi-ations are provided for parameters whose posteriors areapproximately Gaussian. Table 3 covers two scenarios:1) the “primary” solution where i inner > ◦ and 2) the“alternate” solution where i inner < ◦ . The alternatesolution yields retrograde orbits between the inner andouter orbits. We show the phase-folded spectroscopicbinary orbit (identical for both scenarios) in Figure 5,which is in good agreement with that found by Kel-logg et al. (2017). The joint hierarchical fit deliveredmore precise posteriors for many outer orbital parame-ters, though some of the degeneracies noted in Kellogget al. (2017) still remain.Three outer orbit parameters have bimodal posteriordistributions: ω A , Ω outer , and T , outer . The posteriorsfor these outer orbit parameters are identical betweenthe primary and alternate solutions, so the single cor- Table 3.
Orbital parametersParameter primary solution alternate solutionSampled P inner [days] 34 . ± .
001 34 . ± . a inner [mas] 4 . ± .
04 4 . ± . M Ab [ M (cid:12) ] 0 . ± .
01 0 . ± . e inner . ± .
01 0 . ± . i inner [ ◦ ] 131 . ± . . ± . ω Aa a [ ◦ ] 81 ± ± inner b [ ◦ ] 104 ± ± T , inner [JD - 2,450,000] 2704 . ± .
07 2704 . ± . P outer [yrs] 548 ±
244 555 ± M B [ M (cid:12) ] 0 . ± .
28 0 . ± . e outer . ± . . ± . i outer [ ◦ ] 139 ±
13 139 ± ω A a [ ◦ ] · · · c · · · c Ω outer b [ ◦ ] · · · c · · · c T , outer [JD - 2,450,000] · · · c · · · c (cid:36) [ (cid:48)(cid:48) ] 27 . ± .
12 27 . ± . σ ρ [ (cid:48)(cid:48) ] 0 . ± .
004 0 . ± . σ θ [ ◦ ] 0 . ± .
008 0 . ± . σ CfA [km s − ] 3 . ± . . ± . σ Keck [km s − ] 0 . ± . . ± . σ FEROS [km s − ] 3 . ± . . ± . σ du Pont [km s − ] 2 . ± . . ± . − ] − . ± . − . ± . − ] 1 . ± . . ± . − ] − . ± . − . ± . M Aa [ M (cid:12) ] 0 . ± .
01 0 . ± . M A [ M (cid:12) ] 0 . ± .
01 0 . ± . a inner [au] 0 . ± .
001 0 . ± . a outer [au] 63 ±
18 64 ± r p, outer [au] 45 ±
21 45 ± a The argument of periastron of the primary. ω secondary = ω primary + π . b The ascending node is identified as the point where the secondarybody crosses the sky plane receding from the observer. c Posterior is non-Gaussian and not accurately represented by a sum-mary statistic; see Figure 6 v A a [ k m s − ] CfAKeckFEROSdu Pont − O − C − v A b [ k m s − ] . . . . . . − O − C Figure 5.
The phase-folded inner binary orbit and radial veloc-ity residuals. In blue are ten realizations of the inner orbit withthe velocity trend from B removed. The small scatter demon-strates that the inner orbital parameters are tightly constrainedby the data. ner plot in Figure 6 is valid for both scenarios. Rep-resentative outer orbits drawn from the posterior dis-tribution are shown in Figure 7. Although the formaluncertainties and inferred jitter values of the Keck v B measurements are large (0 .
59 km s − and 0 .
82 km s − ,respectively), the actual scatter of the four measuredvalues is substantially smaller. The first two measure-ments (separated by 48 days in 2002/2003) differ by only0 .
05 km s − . The second two measurements (separatedby one year in 2009/2010) differ by only 0 .
15 km s − .While circumstantial, this does raise the possibility thatthe uncertainties on the v B measurements are overes-timated (potentially driven by the scatter in v Aa and v Ab ) and that the increase in v B over the Keck measure-ment baseline ( ≈ . − over 8 years) may be sig-nificant. If true, a monotonically increasing v B clearlyfavors a single posterior mode, highlighted in blue inFigures 6 & 7 (the mode corresponding to decreas-ing v B is highlighted in orange). The MCMC samplesand PyMC3 models corresponding to both scenarios areavailable online. We used these orbital posteriors and the inferred val-ues of i disk and Ω disk to calculate the mutual inclina- https://zenodo.org/record/4568830 − Ω o u t e r [ ◦ ] − − ω A [ ◦ ] T , o u t e r [ y r ] − Ω outer [ ◦ ] T , outer [yr] v B ↑ v B ↓ Figure 6.
Corner plot of the sampled outer orbital parametersthat have bimodal posterior distributions; contours are 1 σ and2 σ within each mode. This figure is identical for both the “pri-mary” and “alternate” solutions. Samples are color-coded basedon 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 con-tours) was sampled simultaneously, the samples have been bifur-cated for plotting purposes only. tions θ between the disk and the stellar orbits (e.g.,Equation 1, Czekala et al. 2019) under different com-binations of the degenerate orientations. As discussedin Czekala et al. (2019), the sensible imposition ofspherically isotropic priors (i.e., the uniform priors oncos i disk and cos i inner ) results in effective mutual incli-nation priors of p ( θ ) ∝ sin( θ ). These isotropic priorshave the consequence of strongly disfavoring coplanararchitectures—simply because of the small phase spacevolume. To quantify the constraining power of the dataonly, we also report the posteriors re-weighted such thatthe effective prior is flat (i.e., similar to the marginallikelihood p (data | θ )). DISCUSSION4.1.
Mutual inclinations Technically there are yet two more degenerate scenarios where i disk < ◦ but i inner > ◦ , or vice-versa. Since the inferredinclinations of the inner binary and circumbinary disk are alreadyso similar (modulo the degeneracy), we do not consider thesescenarios. As established in § θ disk − inner (cid:46) ◦ under a flat mutual inclination prior( θ disk − inner (cid:46) ◦ under sin θ disk − inner prior). Czekalaet al. (2019) found that circumbinary planets, debrisdisks, and protoplanetary disks around short-period bi-naries ( P (cid:46)
40 days) all have low mutual inclinations.So, the low mutual inclination for the TWA 3A is con-sistent with expectations given its 35-day orbital period.That circumbinary disk mutual inclination trendswith binary period is likely a byproduct of the formationmechanism for tight binary stars, which requires forma-tion at larger distances and migration to present-dayconfigurations (Bate et al. 2002). The close binary frac-tions of T Tauri stars and field stars are similar (Kounkelet al. 2019), implying that this migration occurs quickly,before the class II T Tauri phase. It is unlikely thattertiary interactions (e.g., Fabrycky & Tremaine 2007)are responsible for the majority of tight binaries (Moe &Kratter 2018); rather, migration via a circumbinary diskappears to be the dominant pathway (Tokovinin & Moe2020). TWA 3A represents both the lowest mass binary( M A = 0 . ± . M (cid:12) ) hosting a circumbinary disk andthe longest period binary ( P inner = 34 . ± .
001 days)before the population of mutual inclinations transitionsfrom entirely coplanar systems to a broad dispersion ofmutual inclinations (see Figure 14, Czekala et al. 2019).Though the range of binary periods over which this tran-sition occurs is not yet well defined, recent observationshave shed light on the dispersion of mutual inclinationsat slightly longer binary periods of several months. Therecent measurements of WW Cha by VLTI/GRAVITYhave demonstrated that coplanar, truly low mutual in-clination ( θ < ◦ ) circumbinary disks can and do still Table 4.
Inferred mutual inclinations i disk v B p ( θ ) θ disk − inner θ inner − outer θ disk − outer [ ◦ ] [ ◦ ] [ ◦ ] [ ◦ ] < ↑ sin( θ ) < +17 − +18 − < ↑ flat < +17 − +18 − < ↓ sin( θ ) < +10 − +7 − < ↓ flat < +10 − +7 − > ↑ sin( θ ) <
13 22 +14 − +19 − > ↑ flat < < < > ↓ sin( θ ) <
13 76 +14 − +13 − > ↓ flat < +15 − +14 − Note —One σ asymmetric error bars are reported for all uni-modal distributions. Sixty eight percent confidence upperlimits are reported for one-sided distributions. − − α cos δ [ ] − − ∆ δ [ ] . . . . . v BARY [km s − ] 1 . . . ρ [ ] θ [ ◦ ] v B [ k m s − ] Keck
Figure 7.
Thirty representative outer orbits drawn from the highlighted posterior modes in Figure 6. left : the sky plane, centered onTWA 3A (Aa and Ab are represented by the black star), with the the velocity field of the surrounding circumbinary disk. right : theastrometric 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. exist around longer period ( P = 207 days), eccentric( e = 0 .
45) binaries (GRAVITY Collaboration et al.2021). The measurement of WW Cha’s orbital prop-erties is important because it enhances the contrast ofthe stark transition from coplanar systems to a broaddistribution of mutual inclinations. As a rule, all plan-ets and disks orbiting binaries with P (cid:46)
40 days arecoplanar, yet at binary periods of ∼ −
10 months therealready exists both a coplanar system (WW Cha) anda polar-oriented system (HD 98800B, P = 315 days;Kennedy et al. 2019). Discovering new circumbinarydisks with host binary periods near the “transition re-gion” ( P = 30 −
300 days) and measuring their mutualinclinations will help map out the regimes where migra-tion may deliver a coplanar system.Substantial degeneracies remain in the relationshipbetween 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 < ◦ or i > ◦ and 2) which mode of the( ω A , Ω outer , T , outer ) posterior is correct. The inferredmutual inclinations corresponding to the four permu-tations of these degeneracies are delineated in Ta-ble 4, under two different prior assumptions (sin( θ ) or flat). These separate into a coplanar configu-ration ( θ inner − outer ≈ ◦ ), two orthogonal configura-tions ( θ inner − outer ≈ ◦ ) and a retrograde configuration( θ inner − outer > ◦ ). The field population of triples withprojected outer separations <
50 au exhibits a high de-gree of alignment between inner and outer orbital planes(average mutual inclination < ◦ ; Tokovinin 2017).Considering this backdrop, we suspect that the TWA 3Aa–Ab and A–B orbital planes are nearly coplanar aswell. We provide further hydrodynamical evidence forthis scenario in the next subsection.4.2. Disk truncation
The time-dependent gravitational potential of a bi-nary star will influence the radial extent of a protoplan-etary disk—clearing an interior cavity in a circumbinaryconfiguration and truncating the outer edge in a cir-cumstellar configuration. TWA 3 is noteworthy becauseboth types of disk truncation are present in the samesystem.Several analytical and numerical works have derivedthe radius (usually conveyed as a ratio relative to the bi-nary semimajor axis, r/a ) to which an interior, coplanar,and eccentric binary like Aa–Ab ( e inner = 0 . ± . Since the value of θ disk − inner is low, values of θ inner − outer and θ disk − outer are similar in all cases. We will refer to θ inner − outer inwhat follows but the points apply equally to θ disk − outer as well. r/a ∼ −
3, though they also noted that nonaxisym-metric waves at the rim of the disk make it difficult todefine the edge location uniquely. Miranda & Lai (2015)studied circumbinary disk truncation across a range ofmutual inclinations, finding that the truncation radiuswas smaller for more misaligned systems. For coplanarsystems, their results agreed with Artymowicz & Lubow(1994). Using a suite of numerical simulations, Mirandaet al. (2017) found that the truncation radius of copla-nar circumbinary disks is r/a ∼ . − . e ≈ . r/a ∼ . r/a ∼ −
6. Morerecently, Hirsh et al. (2020) used a smoothed-particlehydrodynamics (SPH) setup to derive 50% thresholdsfor an e ≈ . r inner /a inner ≈ . − .
5. Hirshet al. (2020) noted that Thun et al. (2017)’s discrep-ancy with previous results could be attributable to theirchoice of inner polar grid boundary. For the TWA 3Acircumbinary disk, the best-fit interior cavity radius of r cav /a inner ≈ § r cav /a inner ≈ . e inner = 0 .
63 binary, Thunet al. (2017) found the inner edge would have e disk ≈ . e disk ≈ . − .
35. Mu˜noz & Lithwick (2020) usedlinear theory of perturbed, pressure-supported disks tosolve for the eccentricity profile and showed that the ec-centric modes are concentrated to within r/a < r/a ∼
10 (corresponding to1.7 au for TWA 3A). In simulations, the disk eccentricityis e disk < .
05 after r/a = 10 (Thun et al. 2017; Ragusaet al. 2020). Unfortunately, the scale of the inner rim ofthe circumbinary disk is below the resolution that canbe meaningfully probed by the current ALMA obser-vations. If the outer disk ( r (cid:38) e (cid:46) . r outer /a outer ∼ . − . r outer /a outer ∼ . − .
25 for e = 0 . − .
4. Coplanar configurations are also the most effec-tive at truncating the disk: as the mutual inclination θ inner − outer increases, r outer also increases (about 20%larger for θ inner − outer = 90 ◦ and about 40% larger for θ inner − outer = 135 ◦ ; see Figure 4 of Miranda & Lai2015). Using the 95% enclosed mass limit of 8 . ± . a outer = 63 ±
18 au means that r outer /a outer ≈ . − . Phan-tom code (Price et al. 2018a). The SPH method is wellsuited for misaligned disk simulations given that thereis no preferred geometry and angular momentum is con-served to the accuracy of the time-stepping scheme (seee.g. Price 2012). We used 10 gas particles to modelthe circumbinary disc and set the initial inner and outerradii of the disk to r in = 2 au and r out = 20 au, respec-tively. The surface density initially followed a power-lawprofile (Σ ∝ R − ) and the temperature profile followedthe power law profile T = 34 K ( r/
10 au) − . , as in Priceet al. (2018b). The disc total mass was set to 0 . M (cid:12) ,which allowed us to neglect the disk self-gravity. Fur-thermore, we assumed that the disc is locally isothermal,where the sound speed follows a power-law c s ∝ r − / with H/r = 0 .
05 at r = 10 au. Finally, we adopted amean Shakura-Sunyaev disc viscosity α SS ≈ × − .We included all three stars in the simulation(TWA 3Aa, Ab, and B, each represented by a sink par-ticle), which interact with the gas via gravity and accre-tion (Bate et al. 1995). The accretion radii of Aa andAb were set equal to 0 .
01 au, while the accretion radiusof B was set equal to 10% of the Hill radius. These val-ues ensure that the inner and outer regions of the diskare properly modeled while keeping computational costsreasonable (see e.g. Price et al. 2018b, Cuello et al. 2019and M´enard et al. 2020).The masses and orbits of stars Aa and Ab were ini-tialized to the best-fit values listed in Table 3. The cir-cumbinary disc around Aa and Ab was initially in thesame plane as the Aa–Ab binary orbit. The outer com-panion B was set on a circular ( e outer = 0) orbit withsemi-major axis a outer and mutual inclination with re-spect to the binary orbital plane of θ inner − outer . In prin-ciple, the stellar orbits are allowed to change as mass isaccreted, but since an insignificant amount of materialwas accreted throughout the simulation, the change instellar orbits was negligible. We ran six simulations cor-responding to the coplanar, orthogonal, and retrograde3 i0 o -a45au20 au i90 o -a45au20 au i130 o -a45au20 aui0 o -a80au20 au i90 o -a80au20 au i130 o -a80au20 au Figure 8.
Gas surface density in SPH simulations of the TWA 3A binary, circumbinary disk, and exterior companion TWA 3B, over arange of θ inner − outer mutual inclinations ( i ) and semi-major axes a outer ( a ). The red dots at the center of the circumbinary disk representthe positions of the stars Aa and Ab (on top of each other) after 5 orbits of the outer companion B (not shown in frame). The images areconvolved with an 8 au Gaussian beam (white circle in the top left) for better comparison with the observations shown in Fig. 1. The orbitof B which most closely reproduces the outer radius of the circumbinary disk corresponds to the coplanar one with the smaller semi-majoraxis ( θ inner − outer = 0 ◦ , a outer = 45 au; i0 ◦ -a45 ). Each representation is shown from the perspective of an Earth observer. orbits of θ inner − outer = [0 ◦ , ◦ , ◦ ] and two differentvalues of a outer = [45 au ,
80 au]. We used M B = 0 . M (cid:12) and e outer = 0 for all simulations. A visualization ofthe results (convolved with an 8 au beam) is providedin Figure 8, showing that the outer edge of the diskis sensitive to choices of θ inner − outer and a outer . Thesmallest disks were produced in the coplanar simula-tions ( θ inner − outer = 0 ◦ ). The simulation parameterscorresponding to the i disk > ◦ and v B ↑ scenario( θ inner − outer = 0 ◦ , a = 45 au, e = 0: i0 ◦ -a45au ) deliveran outer disk edge (10 au) that most closely matches theALMA observations. This suggests that the true semi-major axis of the A–B binary lies closer to the lowerrange of its estimate ( a outer = 63 ±
18 au), the trueeccentricity of the A–B binary is significantly non-zero( e outer = 0 . ± . Comparison to other systems
The HD 98800 multiple system, coincidentally in thesame TW Hydra association as TWA 3 (and thus,similarly aged at ∼
10 Myr), is an apt analog sys-tem to TWA 3. Because the HD 98800B circumbi-nary disk is in a hierarchical multiple system, it alsoexperiences both interior and exterior disk truncationforces (from the Ba–Bb binary and from the wider Acompanion, respectively). Kennedy et al. (2019) spa-tially resolved the circumbinary disk with ALMA andconvincingly demonstrated that it is in a circumpo-lar ( θ disk − inner ≈ ◦ ) configuration. Franchini et al.(2019) demonstrated that the relatively small cavity size( r cav /a inner ∼ . − .
5) is a consequence of the reducedtorques from an orthogonally-oriented binary (see alsoMiranda & Lai 2015).The exterior companion HD 98800A (which is itselfa spectroscopic binary Aa–Ab, but here is effectivelytreated as a single star) orbits with a outer = 54 au,4 e = 0 .
52 and θ outer ≈ ◦ . The HD 98800B diskouter edge is (cid:46) r outer /a outer ∼ . . M (cid:12) A9Ve; Dominik et al. 2003) is surrounded by a diskwhose 1.4 mm continuum emission extends to ≈
40 auand whose CO emission extends to ≈
100 au (Wagneret al. 2018; van der Plas et al. 2019). Scattered light ob-servations of the disk revealed spiral arms (Wagner et al.2015) and narrow-lane shadows (Benisty et al. 2017).HD 100453B is an external 0 . M (cid:12) companion locatedat ≈
100 au projected distance (Chen et al. 2006; Collinset al. 2009). Definitive conclusions about the influenceof B on the disk are made difficult by the uncertaintyin its orbit; however, recent analysis by Gonzalez et al.(2020) supports a scenario where the disk and binaryplane are substantially misaligned (60 ◦ ), since coplanarorbits consistent with the astrometric data would be oth-erwise inconsistent with the disk morphology, includingthe spirals and observed velocity field. The spiral fea-tures and narrow lane shadows seen in scattered lightalso suggest a complicated inner disk structure inducedby an undetected, substellar companion interior to thedisk (van der Plas et al. 2019; Rosotti et al. 2020; Nealonet al. 2020). The misalignments in this potentially triplesystem can be explained as being driven by the outer B,which drives the outer disk, substellar companion, andinner disk to precess and occasionally undergo Kozai-Lidov oscillations (Nealon et al. 2020).The LTT 1445ABC triple system, which hosts a tran-siting exoplanet (Winters et al. 2019), also bears men-tioning in the context of TWA 3. LTT 1445ABC consistsof three mid to late M dwarfs in a hierarchical configu-ration: B–C forms the inner binary and A is an outertertiary and the most massive star in the system. Theplanet transits A; the entire stellar system is co-planar.The TWA 3 system is something of a pre-main sequencecounterpoint to LTT 1445ABC, in particular the config-uration of its circumbinary disk contrasts with the factthat in LTT 1445A the planet transits the single star.Though if the TWA 3Aa-Ab stars were considered to-gether, TWA 3A would also be the primary star in thesystem. Kennedy et al. (2019) were unable to break the i disk < ◦ or i disk > ◦ degeneracy, but the computed value of θ outer is sim-ilar for both cases. CONCLUSIONSOur main conclusions from this study of the TWA 3system are as follows. • We detected CO J = 2 − CO J = 2 − • We forward modeled the CO J = 2 − disk =116 . ± . i disk = 48 . ± .
7) and infer the stel-lar mass enclosed by the disk M A = M Aa + M Ab =0 . ± . • We combined the disk dynamical constraints withextant radial velocity and astrometric measure-ments of TWA 3Aa, Ab, and B to infer individualstellar masses 0 . ± . M (cid:12) , 0 . ± . M (cid:12) , and0 . ± . M (cid:12) , respectively. • We drew constraints on the orbital architectureof the system, and inferred that the plane of theinner Aa–Ab binary and its circumbinary diskare coplanar (misalignment < ◦ with 68% confi-dence). There are several degenerate solutions forthe mutual inclination between the orbital planesof the inner (Aa–Ab) and outer (A–B) stellar or-bits, however, SPH simulations lend support to thecoplanar solution. • We found the inner and outer radii of the circumbi-nary disk (0 . − .
75 au and 8 . ± . Spitzer spectrum of TWA 3.IC was supported by NASA through the NASA Hub-ble Fellowship grant HST-HF2-51405.001-A awarded bythe Space Telescope Science Institute, which is oper-ated by the Association of Universities for Researchin Astronomy, Inc., for NASA, under contract NAS5-26555. This project has received funding from theEuropean Union’s Horizon 2020 research and inno-vation programme under the Marie Sk(cid:32)lodowska-Curiegrant agreement No 210021. EC acknowledges NASAgrants 80NSSC19K0506 and NNX15AD95G/NEXSS.The National Radio Astronomy Observatory is a fa-cility of the National Science Foundation operated un-der cooperative agreement by Associated Universities,Inc. This paper makes use of the following ALMAdata: ADS/JAO.ALMA
Software:
CASA (v4.4; McMullin et al. 2007),DiskJockey (Czekala et al. 2015; Czekala et al.2019), RADMC-3D (Dullemond 2012), emcee (Foreman-Mackey et al. 2013), Astropy (Astropy Collaborationet al. 2013), PyMC3 (Salvatier et al. 2016), Theano(Theano Development Team 2016), MCFOST (Pinteet al. 2006),
Phantom (Price et al. 2018a), splash (Price 2007), uvplot (Tazzari 2017) .6 REFERENCES
Aly, H., Dehnen, W., Nixon, C., & King, A. 2015, MNRAS,449, 65, doi: 10.1093/mnras/stv128Andrews, S. M. 2020, ARA&A, 58, 483,doi: 10.1146/annurev-astro-031220-010302Andrews, S. M., Czekala, I., Wilner, D. J., et al. 2010, ApJ,710, 462, doi: 10.1088/0004-637X/710/1/462Andrews, S. M., Huang, J., P´erez, L. M., et al. 2018, ApJ,869, L41, doi: 10.3847/2041-8213/aaf741Anthonioz, F., M´enard, F., Pinte, C., et al. 2015, A&A,574, A41, doi: 10.1051/0004-6361/201424520Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651,doi: 10.1086/173679Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,et al. 2013, A&A, 558, A33,doi: 10.1051/0004-6361/201322068Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M.,Mantelet, G., & Andrae, R. 2018, AJ, 156, 58,doi: 10.3847/1538-3881/aacb21Bate, M. R. 2019, MNRAS, 484, 2341,doi: 10.1093/mnras/stz103Bate, M. R., Bonnell, I. A., & Bromm, V. 2002, MNRAS,336, 705, doi: 10.1046/j.1365-8711.2002.05775.xBate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS,277, 362, doi: 10.1093/mnras/277.2.362Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015,MNRAS, 454, 593, doi: 10.1093/mnras/stv1981Benisty, M., Stolker, T., Pohl, A., et al. 2017, A&A, 597,A42, doi: 10.1051/0004-6361/201629798Brandeker, A., Jayawardhana, R., & Najita, J. 2003, AJ,126, 2009, doi: 10.1086/378057Brogan, C. L., Hunter, T. R., & Fomalont, E. B. 2018,arXiv e-prints, arXiv:1805.05266.https://arxiv.org/abs/1805.05266Chen, X. P., Henning, T., van Boekel, R., & Grady, C. A.2006, A&A, 445, 331, doi: 10.1051/0004-6361:20054122Collins, K. A., Grady, C. A., Hamaguchi, K., et al. 2009,ApJ, 697, 557, doi: 10.1088/0004-637X/697/1/557Correia, S., Zinnecker, H., Ratzka, T., & Sterzik, M. F.2006, A&A, 459, 909, doi: 10.1051/0004-6361:20065545Cuello, N., & Giuppone, C. A. 2019, A&A, 628, A119,doi: 10.1051/0004-6361/201833976Cuello, N., Dipierro, G., Mentiplay, D., et al. 2019,MNRAS, 483, 4114, doi: 10.1093/mnras/sty3325Czekala, I. 2021, MCMC samples for TWA 3 orbitalanalysis, Zenodo, doi: 10.5281/zenodo.4568830Czekala, I., Andrews, S. M., Jensen, E. L. N., et al. 2015,ApJ, 806, 154, doi: 10.1088/0004-637X/806/2/154Czekala, I., Andrews, S. M., Torres, G., et al. 2016, ApJ,818, 156, doi: 10.3847/0004-637X/818/2/156 Czekala, I., Chiang, E., Andrews, S. M., et al. 2019, ApJ,883, 22, doi: 10.3847/1538-4357/ab287bCzekala, I., Jensen, E., & Huang, J. 2019,iancze/DiskJockey: Upgrades to Julia v1.0,doi: 10.5281/zenodo.3235028Dominik, C., Dullemond, C. P., Waters, L. B. F. M., &Walch, S. 2003, A&A, 398, 607,doi: 10.1051/0004-6361:20021629Dullemond, C. P. 2012, RADMC-3D: A multi-purposeradiative transfer tool, Astrophysics Source CodeLibrary. http://ascl.net/1202.015Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298,doi: 10.1086/521702Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman,J. 2013, PASP, 125, 306, doi: 10.1086/670067Foreman-Mackey, D., Luger, R., Czekala, I., et al. 2020,exoplanet-dev/exoplanet v0.3.2,doi: 10.5281/zenodo.1998447Foucart, F., & Lai, D. 2013, ApJ, 764, 106,doi: 10.1088/0004-637X/764/1/106Franchini, A., Lubow, S. H., & Martin, R. G. 2019, ApJL,880, L18, doi: 10.3847/2041-8213/ab2fd8Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al.2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051Gonzalez, J.-F., van der Plas, G., Pinte, C., et al. 2020,MNRAS, 499, 3837, doi: 10.1093/mnras/staa2938GRAVITY Collaboration, Eupen, F., Labadie, L., et al.2021, arXiv e-prints, arXiv:2102.00122.https://arxiv.org/abs/2102.00122Herczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, 97,doi: 10.1088/0004-637X/786/2/97Hirsh, K., Price, D. J., Gonzalez, J.-F., Ubeira-Gabellini,M. G., & Ragusa, E. 2020, MNRAS, 498, 2936,doi: 10.1093/mnras/staa2536Hoffman, M. D., & Gelman, A. 2011, arXiv e-prints,arXiv:1111.4246. https://arxiv.org/abs/1111.4246Janson, M., Bergfors, C., Brandner, W., et al. 2014, ApJS,214, 17, doi: 10.1088/0067-0049/214/2/17Jayawardhana, R., Hartmann, L., Fazio, G., et al. 1999,ApJL, 521, L129, doi: 10.1086/312200Kellogg, K., Prato, L., Torres, G., et al. 2017, ApJ, 844,168, doi: 10.3847/1538-4357/aa7c60Kennedy, G. M., Matr`a, L., Facchini, S., et al. 2019, NatureAstronomy, 3, 230, doi: 10.1038/s41550-018-0667-xKnapp, W., & Nanson, J. 2018, Journal of Double StarObservations, 14, 503Kounkel, M., Covey, K., Moe, M., et al. 2019, AJ, 157, 196,doi: 10.3847/1538-3881/ab13b1 Kuruwita, R. L., & Federrath, C. 2019, MNRAS, 486, 3647,doi: 10.1093/mnras/stz1053Larson, R. B. 1969, MNRAS, 145, 271,doi: 10.1093/mnras/145.3.271Lindegren, L., Hern´andez, J., Bombrun, A., et al. 2018,A&A, 616, A2, doi: 10.1051/0004-6361/201832727Martin, R. G., & Lubow, S. H. 2017, ApJL, 835, L28,doi: 10.3847/2041-8213/835/2/L28Mason, B. D., Hartkopf, W. I., Miles, K. N., et al. 2018,AJ, 155, 215, doi: 10.3847/1538-3881/aab9b8Mason, B. D., Wycoff, G. L., Hartkopf, W. I., Douglass,G. G., & Worley, C. E. 2001, AJ, 122, 3466,doi: 10.1086/323920McJunkin, M., France, K., Schneider, P. C., et al. 2014,ApJ, 780, 150, doi: 10.1088/0004-637X/780/2/150McMullin, J. P., Waters, B., Schiebel, D., Young, W., &Golap, K. 2007, in Astronomical Society of the PacificConference Series, Vol. 376, Astronomical Data AnalysisSoftware and Systems XVI, ed. R. A. Shaw, F. Hill, &D. J. Bell, 127M´enard, F., Cuello, N., Ginski, C., et al. 2020, A&A, 639,L1, doi: 10.1051/0004-6361/202038356Miranda, R., & Lai, D. 2015, MNRAS, 452, 2396,doi: 10.1093/mnras/stv1450Miranda, R., Mu˜noz, D. J., & Lai, D. 2017, MNRAS, 466,1170, doi: 10.1093/mnras/stw3189Moe, M., & Kratter, K. M. 2018, ApJ, 854, 44,doi: 10.3847/1538-4357/aaa6d2Mu˜noz, D. J., & Lithwick, Y. 2020, ApJ, 905, 106,doi: 10.3847/1538-4357/abc74cNealon, R., Cuello, N., Gonzalez, J.-F., et al. 2020,MNRAS, 499, 3857, doi: 10.1093/mnras/staa2721Offner, S. S. R., Dunham, M. M., Lee, K. I., Arce, H. G., &Fielding, D. B. 2016, ApJL, 827, L11,doi: 10.3847/2041-8205/827/1/L11Offner, S. S. R., Kratter, K. M., Matzner, C. D., Krumholz,M. R., & Klein, R. I. 2010, ApJ, 725, 1485,doi: 10.1088/0004-637X/725/2/1485Pinte, C., Dent, W. R. F., M´enard, F., et al. 2016, ApJ,816, 25, doi: 10.3847/0004-637X/816/1/25Pinte, C., M´enard, F., Duchˆene, G., & Bastien, P. 2006,A&A, 459, 797, doi: 10.1051/0004-6361:20053275Price, D. J. 2007, PASA, 24, 159, doi: 10.1071/AS07022—. 2012, Journal of Computational Physics, 231, 759,doi: 10.1016/j.jcp.2010.12.011Price, D. J., Wurster, J., Tricco, T. S., et al. 2018a, PASA,35, e031, doi: 10.1017/pasa.2018.25Price, D. J., Cuello, N., Pinte, C., et al. 2018b, MNRAS,477, 1270, doi: 10.1093/mnras/sty647 Ragusa, E., Alexander, R., Calcino, J., Hirsh, K., & Price,D. J. 2020, MNRAS, 499, 3362,doi: 10.1093/mnras/staa2954Reipurth, B., & Zinnecker, H. 1993, A&A, 278, 81Remijan, A., Biggs, A., Cortes, P., et al. 2019, ALMA Cycle7 Technical Handbook, doi: 10.5281/zenodo.4511522Ribas, ´A., Mac´ıas, E., Espaillat, C. C., & Duchˆene, G.2018, ApJ, 865, 77, doi: 10.3847/1538-4357/aad81bRosotti, G. P., Benisty, M., Juh´asz, A., et al. 2020,MNRAS, 491, 1335, doi: 10.1093/mnras/stz3090Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJComputer Science, 2, e55Sheehan, P. D., Wu, Y.-L., Eisner, J. A., & Tobin, J. J.2019, ApJ, 874, 136, doi: 10.3847/1538-4357/ab09f9Tazzari, M. 2017, mtazzari/uvplot: v0.1.1,doi: 10.5281/zenodo.1003113Theano Development Team. 2016, arXiv e-prints,abs/1605.02688. http://arxiv.org/abs/1605.02688Thun, D., Kley, W., & Picogna, G. 2017, A&A, 604, A102,doi: 10.1051/0004-6361/201730666Tofflemire, B. M., Mathieu, R. D., Herczeg, G. J., Akeson,R. L., & Ciardi, D. R. 2017, ApJ, 842, L12,doi: 10.3847/2041-8213/aa75cbTofflemire, B. M., Mathieu, R. D., & Johns-Krull, C. M.2019, AJ, 158, 245, doi: 10.3847/1538-3881/ab4f7dTokovinin, A. 2017, ApJ, 844, 103,doi: 10.3847/1538-4357/aa7746Tokovinin, A., Mason, B. D., Hartkopf, W. I., Mendez,R. A., & Horch, E. P. 2015, AJ, 150, 50,doi: 10.1088/0004-6256/150/2/50Tokovinin, A., & Moe, M. 2020, MNRAS, 491, 5158,doi: 10.1093/mnras/stz3299Tripathi, A., Andrews, S. M., Birnstiel, T., & Wilner, D. J.2017, ApJ, 845, 44, doi: 10.3847/1538-4357/aa7c62van der Plas, G., M´enard, F., Gonzalez, J. F., et al. 2019,A&A, 624, A33, doi: 10.1051/0004-6361/201834134Wagner, K., Apai, D., Kasper, M., & Robberto, M. 2015,ApJL, 813, L2, doi: 10.1088/2041-8205/813/1/L2Wagner, K., Dong, R., Sheehan, P., et al. 2018, ApJ, 854,130, doi: 10.3847/1538-4357/aaa767Weaver, E., Isella, A., & Boehler, Y. 2018, ApJ, 853, 113,doi: 10.3847/1538-4357/aaa481Webb, R. A., Zuckerman, B., Platais, I., et al. 1999, ApJL,512, L63, doi: 10.1086/311856Weintraub, D. A., Saumon, D., Kastner, J. H., & Forveille,T. 2000, ApJ, 530, 867, doi: 10.1086/308402Winters, J. G., Medina, A. A., Irwin, J. M., et al. 2019, AJ,158, 152, doi: 10.3847/1538-3881/ab364dZanazzi, J. J., & Lai, D. 2018, MNRAS, 473, 603,doi: 10.1093/mnras/stx2375 Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJL, 877,L18, doi: 10.3847/2041-8213/ab1f8c A. CO J = 2 − -8.4 km s − -7.6 -6.8 -6.0 -5.2 -4.4 -3.6 -2.8-2.0 -1.2 -0.4 0.4 1.2 2.0 2.8 3.6 − . . . α cos δ [ ] − . . . ∆ δ [ ] − . − . − . . . . . − − − × RMSmJy beam − Figure 9.
Continuum-subtracted CO J = 2 −−