The Architecture of the GW Ori Young Triple Star System and Its Disk: Dynamical Masses, Mutual Inclinations, and Recurrent Eclipses
Ian Czekala, Sean M. Andrews, Guillermo Torres, Joseph E. Rodriguez, Eric L. N. Jensen, Keivan G. Stassun, David W. Latham, David J. Wilner, Michael A. Gully-Santiago, Konstantin N. Grankin, Michael B. Lund, Rudolf B. Kuhn, Daniel J. Stevens, Robert J. Siverd, David James, B. Scott Gaudi, Benjamin J. Shappee, Thomas W.-S. Holoien
DDraft version November 15, 2018
Typeset using L A TEX twocolumn style in AASTeX61
THE ARCHITECTURE OF THE GW ORI YOUNG TRIPLE STAR SYSTEM AND ITS DISK:DYNAMICAL MASSES, MUTUAL INCLINATIONS, AND RECURRENT ECLIPSES
Ian Czekala, ∗ Sean M. Andrews, Guillermo Torres, Joseph E. Rodriguez, Eric L. N. Jensen, Keivan G. Stassun,
4, 5
David W. Latham, David J. Wilner, Michael A. Gully-Santiago, Konstantin N. Grankin, Michael B. Lund, Rudolf B. Kuhn, Daniel J. Stevens, Robert J. Siverd, David James, † B. Scott Gaudi, Benjamin J. Shappee, ‡ and Thomas W.-S. Holoien § Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA 94305, USA Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA Department of Physics and Astronomy, Swarthmore College, 500 College Avenue, Swarthmore, PA 19081, USA Department of Physics and Astronomy, Vanderbilt University, 6301 Stevenson Center, Nashville, TN 37235, USA Department of Physics, Fisk University, Nashville, TN 37208, USA NASA Ames Research Center, Moffett Field, CA 94035, USA Crimean Astrophysical Observatory, pos. Nauchnyi, Crimea, 298409 Russia South African Astronomical Observatory, P.O. Box 9, Observatory 7935, South Africa Department of Astronomy, The Ohio State University, Columbus, OH 43210, USA Las Cumbres Observatory Global Telescope Network, 6740 Cortona Drive, Suite 102, Santa Barbara, CA 93117, USA The Observatories of the Carnegie Institution for Science, 813 Santa Barbara St., Pasadena, CA 91101 (Received October 9, 2017; Revised November 15, 2017; Accepted November 17, 2017)
Submitted to ApJABSTRACTWe present spatially and spectrally resolved Atacama Large Millimeter/submillimeter Array (ALMA) observationsof gas and dust orbiting the pre-main sequence hierarchical triple star system GW Ori. A forward-modeling of the CO and C O J =2–1 transitions permits a measurement of the total stellar mass in this system, 5 . ± . M (cid:12) ,and the circum-triple disk inclination, 137 . ± . ◦ . Optical spectra spanning a 35 year period were used to derive newradial velocities and, coupled with a spectroscopic disentangling technique, revealed that the A and B components ofGW Ori form a double-lined spectroscopic binary with a 241 . ± .
05 day period; a tertiary companion orbits thatinner pair with a 4218 ±
50 day period. Combining the results from the ALMA data and the optical spectra withthree epochs of astrometry in the literature, we constrain the individual stellar masses in the system ( M A ≈ . M (cid:12) , M B ≈ . M (cid:12) , M C ≈ . M (cid:12) ) and find strong evidence that at least one (and likely both) stellar orbital planes aremisaligned with the disk plane by as much as 45 ◦ . A V -band light curve spanning 30 years reveals several new ∼
30 dayeclipse events 0.1–0.7 mag in depth and a 0.2 mag sinusoidal oscillation that is clearly phased with the AB–C orbitalperiod. Taken together, these features suggest that the A–B pair may be partially obscured by material in the innerdisk as the pair approaches apoastron in the hierarchical orbit. Lastly, we conclude that stellar evolutionary modelsare consistent with our measurements of the masses and basic photospheric properties if the GW Ori system is ∼ Keywords: protoplanetary disks – stars: fundamental parameters – stars: pre-main sequence – stars:individual (GW Ori)
Corresponding author: Ian [email protected] a r X i v : . [ a s t r o - ph . E P ] N ov INTRODUCTIONPre-main sequence (pre-MS) stars in multiple systems—for which it is possible to precisely measure their fun-damental stellar properties through dynamical means—serve as touchstones for understanding the final stages ofstellar formation and the conditions under which plane-tary systems are assembled. While recent decades haveseen steady progress in understanding binary formationin general, lingering uncertainties still remain aboutthe characteristics of young spectroscopic binaries andhigher order systems (Duchˆene & Kraus 2013).GW Ori, a G-type star associated with the λ OrionisOB star-forming complex (Dolan 2000; Dolan & Math-ieu 2001, 2002), was one of the first T Tauri stars tobe revealed as a spectroscopic binary, with a period of240 days (Mathieu et al. 1991). Radial velocity (RV)monitoring hinted at the presence of a third body with aperiod of ∼
10 years; a tertiary was confirmed directly us-ing infrared interferometry (Berger et al. 2011). Circum-stellar material in the GW Ori system was first inferredfrom infrared excess emission (Mathieu et al. 1991); asubsequent detection of the dust continuum at submil-limeter wavelengths suggested the disk was especiallymassive and must be circumbinary ( M disk (cid:38) . M (cid:12) ;Mathieu et al. 1995).The disk material provided a natural explanationfor the quasi-periodic optical dimming of GW Ori over ∼
30 day durations: the suspicion was that a disk aroundthe secondary was eclipsing the primary, presuming anearly edge-on viewing angle (Shevchenko et al. 1992,1998). Fang et al. (2017) spatially resolved the diskmaterial with Submillimeter Array (SMA) observationsof the dust continuum and the line emission from COisotopologues, demonstrating its large radial extent andtherefore presumably circum-triple architecture. How-ever, they found the disk has an intermediate inclinationto the line of sight ( i disk ≈ ◦ ), in apparent conflictwith the eclipse model. Indeed, this adds to a collectionof indirect evidence for a more complicated geometry inthe inner disk, including mid-infrared fluxes that varyon ∼ year timescales (Fang et al. 2014), and CO rovi- ∗ KIPAC Postdoctoral Fellow † Event Horizon Telescope ‡ Hubble, Carnegie-Princeton Fellow § Carnegie Fellow brational emission lines with multi-component profiles(which requires a complicated geometry and/or temper-ature structure in the inner disk; Najita et al. 2003).Beyond resolving outstanding questions about its ar-chitecture, the GW Ori system presents an excellent op-portunity to obtain a precise dynamical mass measure-ment for an earlier type ( ∼ G8) star at a very youngage. Precise dynamical masses are crucial for calibrat-ing the photospheric predictions (e.g., T eff , L ) of stellarevolutionary models, and the region of the Hertzsprung-Russell (HR) diagram occupied by GW Ori is particu-larly sparsely populated with benchmark systems (Stas-sun et al. 2014). Especially interesting is the fact thatthree different dynamical mass measurement techniquescan be employed to study the GW Ori system: (1) RVmonitoring, which constrains the mass ratios of the stars(Mathieu et al. 1991; Fang et al. 2014); (2) astromet-ric monitoring, which provides the inclinations of theorbits and, when coupled with RV measurements, canreveal the individual component masses (Berger et al.2011); and (3) the disk-based dynamical mass technique(e.g., Simon et al. 2000, 2017; Rosenfeld et al. 2012;Czekala et al. 2015a, 2016), which measures the total stellar mass.In this paper, we combine information from each ofthese techniques to better understand the fundamentalproperties and underlying physical architecture of theGW Ori system. Section 2 describes new ALMA ob-servations of the GW Ori disk, an updated analysis of35 years’ worth of optical spectroscopic RV monitoring,and an extensive decades-long optical photometric cat-alog. Section 3 describes our tomographic reconstruc-tion of the disk velocity field, a new analysis of the RVdata, their combination with literature-based astromet-ric constraints (Berger et al. 2011), an assessment of thesystem parameters and geometry, and the connectionsto the observed photometric variations. Section 4 con-cludes by discussing the structure and orientation of thedisk with respect to the orbital architecture of the triplesystem, and considers the GW Ori system in the contextof other young multiple pre-MS systems. OBSERVATIONS AND DATA REDUCTION2.1.
Millimeter Interferometry
GW Ori was observed with the ALMA interferome-ter on 2015 May 14 (program ID 2012.1.00496.S), with ∆ α [ ] ∆ δ [ ] continuum CO CO C O 5 10 15 20
LSR velocity [km / s] F ν [ J y ] CO CO (x4)C O (x12)
Figure 1. ( left ) A 226 GHz continuum image. Contours start at 5 × the RMS noise level and increase by factors of 2. The synthesizedbeam geometry is shown in the lower left corner. ( middle, left to right ) Maps of the CO, CO, and C O velocity-integrated intensities(contours, starting at 10, 3, and 3 × the RMS noise levels, respectively, and increasing by factors of 2) overlaid on the intensity-weightedprojected velocities (color-scale). Note the prominent molecular cloud contamination in the CO map (see also Fig. 2). ( right ) Spatiallyintegrated spectra (inside the same
CLEAN mask, and smoothed with an 0.85 km s − Hanning kernel) for each CO line.
37 of the 12 m main array antennas configured to spanbaselines of 23–558 m. The double sideband Band 6 re-ceivers were employed in dual polarization mode, andthe ALMA correlator was set up to process data in 4spectral windows (SPWs). Two of these SPWs, centeredat 220.426 and 230.450 GHz to observe the CO and CO J =2–1 transitions (at rest frequencies of 220.399and 230.538 GHz, respectively), covered 234 MHz ofbandwidth in 3840 channels (a 61 kHz channel spac-ing). One other sampled 469 MHz around 219.763 GHzto observe the C O J =2–1 transition (at rest frequency219.560 GHz) with 3840 channels (a 122 kHz channelspacing). The last SPW sampled the continuum in a1.875 GHz range around 231.956 GHz using 128 coarsechannels (a 15.625 MHz channel spacing).The observations cycled between GW Ori and thequasar J0510+1800 with a 7 minute cadence. Thequasar J0423-0120 and Ganymede were observed asbandpass and flux calibration sources, respectively, atthe start of the execution block. The total on-sourceintegration time for GW Ori was 16 minutes. The ob-serving conditions were typical for Band 6 projects, witha precipitable water vapor level around 1.1 mm.The visibility data were calibrated with standard pro-cedures using the CASA software package (v4.4). Theraw, observed visibility phases were adjusted based onthe contemporaneous measurements of water vapor ra-diometers, flagged when applicable, and then the band-pass shape in each SPW was calibrated based on the ob-servations of J0423-0120. The absolute amplitude scalewas determined based on the observations of Ganymede.The complex gain behavior of the array and atmospherewas corrected based on the repeated observations ofJ0510+1800. The calibrated visibilities showed a strongcontinuum signal, suggesting that self-calibration couldsignificantly improve the data quality. An initial modelbased on a preliminary continuum image was used fortwo rounds of phase-only self-calibration (on 30 s, then 6 s intervals) and one additional round that includedthe amplitudes (on a 7 minute scan interval). Thisself-calibration reduced the RMS noise level in the con-tinuum by a factor of ∼
40. After applying the self-calibration tables to the entire dataset (channel by chan-nel), we parsed out data products for each individualemission tracer of interest. A set of continuum visibili-ties was constructed by spectrally averaging the line-freechannels in each SPW into ∼
125 MHz increments. Thespectral visibilities for the CO, CO, and C O lineswere continuum-subtracted and regridded into 170 ms − -wide channels in the LSRK restframe over a ∼
10 kms − range around the line centers.These fully reduced visibility sets were then imagedby Fourier inversion assuming a Briggs (robust=0.5)weighting scheme and deconvolution with the standard CLEAN algorithm. Some basic image properties for thesynthesized continuum image and spectral line imagecubes are listed in Table 1. The continuum and spectralline moment maps are shown together in Figure 1, alongwith a comparison of the integrated spectra. The chan-nel maps for individual lines are compiled in Figure 2.
Table 1.
ALMA Image Properties RMSbeam dimensions mJy beam −
226 GHz continuum 0 . (cid:48)(cid:48) × . (cid:48)(cid:48)
54, 126 ◦ CO J =2 − . (cid:48)(cid:48) × . (cid:48)(cid:48)
56, 126 ◦ CO J =2 − . (cid:48)(cid:48) × . (cid:48)(cid:48)
59, 126 ◦ O J =2 − . (cid:48)(cid:48) × . (cid:48)(cid:48)
58, 126 ◦ Note —The RMS noise levels recorded for the spectral linecubes correspond to the values per 170 m s − channel. C O − − . − . . . − −
20 0 20 40 − − α [ ] − − ∆ δ [ ] C O − − . − . . . . − −
10 0 10 20 − − α [ ] − − ∆ δ [ ] C O − − . − . − .
02 0 .
00 0 .
02 0 . − − − − α [ ] − − ∆ δ [ ] Jy/beam × RMSJy/beam × RMSJy/beam × RMS
Figure 2.
Channel maps of the CO, CO, and C O ( from top to bottom ) line emission from the GW Ori disk. Each channel representsthe emission in a 170 m s − -wide velocity bin. LSRK velocities are indicated in the upper left, and synthesized beam sizes in the lower leftof each panel. Scale bars are provided at the bottom right of each set of channel maps. The 226 GHz (1.3 mm) continuum map shows a bright(flux density = 202 ±
20 mJy), compact but marginallyresolved (deconvolved Gaussian FWHM ≈ . (cid:48)(cid:48)
9) sourcecentered on the GW Ori stellar system, with a peakintensity of 67 mJy beam − (S/N ≈ ±
60 mJy), but marginallydiscrepant with that of Fang et al. (2017, 320 ±
64 mJy).A crude estimate of the emission geometry (from a Gaus-sian fit to the visibilities) suggests an inclination of35–40 ◦ , with the major axis oriented ∼ ◦ E of N.The CO isotopologue channel maps reveal bright (in-tegrated intensities of 41.8, 5.7, and 0.8 Jy km s − for CO, CO, and C O, respectively) and extended(FWHM ∼ . (cid:48)(cid:48)
5) emission that is clearly in rotationaround the continuum centroid, spanning a projected ve-locity range of ± − from the line center. The lineemission is blueshifted to the south and redshifted to thenorth, consistent with the orientation estimated fromthe continuum emission. The peak intensities for eachline are ∼ − in the bright-est channels (peak S/N ≈ CO, CO, and C O, respectively. The CO channel mapsshow some clear evidence for structured contaminationfrom the surrounding molecular cloud, particularly as astreamer to the west at ∼ − and some diffuseclumps to the north around 13–14 km s − , confirmingthe “tail”-like feature seen by Fang et al. (2017). Thesestructures are much fainter, but still present, in COemission; they are not apparent in the C O maps.2.2.
Optical Spectroscopy
GW Ori was monitored spectroscopically at theHarvard-Smithsonian Center for Astrophysics for morethan 35 years, beginning in 1981 November. A total of203 usable spectra were gathered through 2009 Aprilusing three nearly identical echelle spectrographs (Dig-ital Speedometers, DS; now decommissioned) with aresolving power of R ≈ ,
000 mounted on three dif-ferent telescopes: the 1.5 m Tillinghast reflector at theFred L. Whipple Observatory (Mount Hopkins, AZ),the 4.5 m-equivalent Multiple Mirror Telescope (also onMount Hopkins) before conversion to a monolithic mir-ror, and occasionally on the 1.5 m Wyeth reflector atthe Oak Ridge Observatory (in the town of Harvard,MA). Each instrument was equipped with an intensifiedphoton-counting Reticon detector limiting the outputto a single echelle order 45 ˚A wide, which was cen-tered on the region of the Mg I b triplet at 5187 ˚A (seeLatham 1992). The signal-to-noise ratios of these ob-servations range from 14 to 59 per resolution element of8.5 km s − , with a median of 41. Wavelength calibra-tions were based on exposures of a thorium-argon lamptaken before and after each science exposure. Reduc-tions were performed with a dedicated pipeline, and thezero-point of the velocities was monitored regularly bymeans of exposures of the evening and morning twilightsky. The original analysis of Mathieu et al. (1991) useda subset of 45 of these spectra. A further 79 usablespectra of GW Ori were collected with the TillinghastReflector Echelle Spectrograph (TRES; F˝ur´esz 2008), a bench-mounted, fiber-fed echelle instrument attached tothe 1.5 m Tillinghast reflector. It provides a resolvingpower of R ≈ , − , with a median of74. Wavelength calibration was carried out as above,and reductions were performed as described by Buch-have et al. (2010). RV standard stars were observed eachnight to monitor the zero point and place it on the samesystem as the DS observations to within ∼ − .All of our spectra appear to be single-lined, with broadfeatures indicative of significant rotation. PreliminaryRV measurements were therefore made with standardone-dimensional cross-correlation techniques, as in theanalysis of Mathieu et al. (1991). However, severalpieces of evidence suggested it should be possible to de-tect the lines of the secondary in the 240 day binary. Inparticular, the large flux ratio of f B /f A = 0 . ± . H -band, when translated to the optical, would stillbe significant for any reasonable assumption of the effec-tive temperatures, making our non-detection of the sec-ondary somewhat surprising. Furthermore, those sameauthors proposed that the system is observed nearlyface-on, which would lead to strong line blending thatcould explain our lack of detection despite the sizablebrightness of the secondary.We embarked on a search for such a signature usingthe then-under-development PSOAP spectroscopic disen-tangling package (Czekala et al. 2017), with the assump-tion that it must be at or near the detection limit (e.g., q in (cid:46) .
2) since it had not been previously seen. Givena time-series of high resolution spectroscopic observa-tions covering the orbital phase of the binary or triplestar,
PSOAP simultaneously infers the intrinsic spectrumof each star along with the stellar orbit using Gaussianprocesses as a modeling basis. This provides a robustprobabilistic inference of both the orbits and spectrain a purely data-driven manner, which can further beused to measure fundamental properties with traditionalanalysis techniques. Preliminary results hinted at thedetection of the secondary, but for mass ratios muchlarger than expected ( q in > . double -linedbinary based on high resolution infrared spectroscopy( q in ∼ .
55, Prato et al. (2017), submitted ). Motivatedby that result, we renewed our efforts to search for thesecondary using
PSOAP and a targeted
TODCOR analysis,
Table 2.
Heliocentric RV measurements of GW Ori.HJD RV A σ A RV B σ B [2,400,000+] [km s − ] [km s − ] [km s − ] [km s − ]44919.0042 31.24 5.40 28.50 19.1345301.8865 25.10 5.18 25.29 18.3545336.7941 23.46 3.36 20.37 11.9245708.7038 33.05 5.83 20.37 20.6545709.6058 37.70 2.76 25.18 9.77 Note —Observations up to HJD 2,454,926.6573 were ob-tained with the DS, and the remainder with TRES. Thistable is available in its entirety in machine-readable form. which is a two-dimensional cross-correlation algorithmdesigned to minimize biases in the RVs due to line blend-ing.At present, one limitation of the
PSOAP framework(and Gaussian processes in general, to some extent) isthe extreme computational expense in performing largematrix calculations. This generally limits us to con-sidering fewer than 20 epochs of high resolution spec-tra at a time, which consequentially limits the complex-ity of the orbital model that can be used. Although itwas straightforward to extend the framework to utilizea hierarchical triple orbital model and three Gaussianprocess components, we found that we were unable toemploy enough spectroscopic epochs to sufficiently con-strain the more complex orbital model. Therefore, weexperimented using different subsets of the highest S/Ndata in the range of 5060 - 5290˚A to test our sensitivityto the presence of the secondary and tertiary spectralsignatures. In all of these tests, we clearly detected thefeatures of the secondary but found no obvious evidencefor spectroscopic signatures of the tertiary.Based on the guidance from these results, we re-examined our spectra with
TODCOR and succeeded in de-tecting the secondary via cross-correlation, as well. Asanticipated, the lines of the two stars are always heav-ily blended, which causes a strong degeneracy betweenthe adopted rotational line broadening for the templates(see below), the velocity amplitudes, the adopted tem-peratures, and the flux ratio. To measure RVs weadopted synthetic templates from the PHOENIX libraryof Husser et al. (2013), broadened to match the reso-lution of our spectra. For the TRES observations werestricted our analysis to the 100 ˚A order centered onthe Mg I b triplet, both for consistency with the analy- sis of the DS spectra, which cover only a 45 ˚A windowcentered on this region, and because experience showsthat it contains most of the information on the veloci-ties. The one-dimensional cross-correlations needed toconstruct the 2-D correlation function in TODCOR werecomputed using the IRAF task XCSAO (Kurtz & Mink1998). The template parameters were selected basedon an analysis of the stronger TRES spectra, as fol-lows. For the primary star we adopted a temperature of T eff = 5700 K proposed by Mathieu et al. (1991), alongwith log g = 3 . v sin i ) of the primary, the secondary tem-perature, and the secondary v sin i were then determinedby running extensive grids of 2-D cross-correlations overbroad ranges in each parameter in a manner similar tothat described by Torres et al. (2002), seeking the bestmatch between the templates and the observed spec-tra as measured by the peak cross-correlation coefficientaveraged over all exposures. For each combination oftemplate parameters we also determined the flux ratiothat maximizes the correlation.In this way we determined a best-fit secondary tem-perature of T eff = 4800 ±
200 K, and v sin i values forthe primary and secondary of 40 and 45 km s − , respec-tively, with estimated uncertainties of 5 km s − . Themeasured flux ratio in the Mg I b 5187 ˚A region is f B /f A = 0 . ± .
05. While in principle these tem-peratures and v sin i values are merely free parametersthat provide the best match to the observed spectra, inthe following we interpret them also as estimates of thephysical properties of the stars. The RVs we measuredfrom our DS and TRES spectra with these parametersare reported in Table 2, along with their uncertainties.Typical errors for the primary and secondary are 1.0 and2.7 km s − for TRES, and 2.5 and 8.7 km s − for the DSmeasurements, though individual errors can sometimesbe much larger. Despite the use of TODCOR , we reiteratethat the severe line blending at all phases of the innerorbit caused by a combination of rotational broaden-ing and small velocity amplitudes makes the RVs verysusceptible to errors in the template parameters (par-ticularly v sin i ) and in the adopted flux ratio, and asa result the orbital elements presented later may suf-fer from systematic errors not included in the statisti-cal uncertainties. Nevertheless, as a consistency check IRAF is distributed by the National Optical Astronomy Ob-servatories, which are operated by the Association of Universi-ties for Research in Astronomy, Inc., under cooperative agreementwith the National Science Foundation. we used PHOENIX spectra from Husser et al. (2013)for the primary and secondary stellar parameters givenabove to extrapolate our measured flux ratio at 5187 ˚Ato the near infrared, and obtained an H -band value of f B /f A = 0 . ± .
12. While less precise than the Bergeret al. (2011) measurement, the agreement is excellent.2.3.
Time-series Photometry
We have assembled a high cadence lightcurve ofGW Ori covering a ∼
30 year timespan by drawing fromseveral ongoing photometric surveys as well as archivalobservations. The details of the different surveys usedin compiling this catalog are described below.2.3.1.
Maidanak Observatory
Photoelectric
U BV R observations of GW Ori were ob-tained at Mount Maidanak Observatory in Uzbekistan.About 530
U BV R magnitudes were collected from 1987to 2003, although the number of U measurements is rel-atively small compared to the other photometric bands.All observations were performed with three telescopes(one 0.48 m and two 0.6 m reflectors) using identicalsingle-channel pulse-counting photometers with FEU-79photomultiplier tubes. The observations of GW Ori werecarried out as part of the ROTOR program, which wasdescribed by Shevchenko et al. (1993).The RMS uncertainty for a single measurement in theinstrumental system was 0 .
01 mag in
BV R and 0 .
02 magin U . Observations were carried out either differentiallyusing a nearby reference star or directly by estimatingthe nightly extinction. In the latter case, several refer-ence stars were observed every night to derive the extinc-tion coefficients in each filter. Selected standard starswere observed and used to calibrate instrumental mag-nitudes on the Cousins system. We then transformed themagnitudes to the Johnson U BV R system using the re-lationship: ( V − R ) C = − . . V − R ) J . Thesystematic uncertainty in this conversion is 0 .
01 mag.2.3.2.
KELT
The Kilodegree Extremely Little Telescope (KELT)project uses two telescopes to survey over 70% ofthe entire sky searching for transiting planets aroundbright stars (8 < V < ◦ × ◦ field-of-view(FOV), and a 23 (cid:48)(cid:48) pixel scale. Both telescopes use abroad R -band filter. KELT observes using a ParamountME German equatorial mount with a 180 ◦ meridian flip;therefore KELT observes in either an “east” or “west”orientation. The telescope optics are not perfectly ax-isymmetric, and so the point spread function (PSF) changes from one orientation to the other. Throughoutthe data reduction process, the east and west obser-vations are treated as though they were acquired fromseparate telescopes. For GW Ori specifically, the PSFasymmetry results in a 0.2 mag systematic offset be-tween the east and west orientations. See Siverd et al.(2012) and Kuhn et al. (2016) for a detailed descriptionof the KELT observing strategy and reduction process.GW Ori was located in KELT-South field 05 ( α = 06hr07m 48.0s, δ = +3 ◦ (cid:48) (cid:48)(cid:48) ) and was observed 2889times from UT 2010 February 28 until UT 2015 April09, with a median uncertainty of 0.005 mag.2.3.3. ASAS
Using two observing locations, in Las Campanas,Chile and Haleakala, Maui, the All-Sky Automated Sur-vey (ASAS) project was designed to observe the entiresky to a limiting optical magnitude of 14. The twoobservatory setups each contained a wide-field Minolta200/2.8 APO-G telephoto lenses with a 2K ×
2K ApogeeCCD and both observed simultaneously in B and V band. The telescope and camera set up correspond to a8 . ◦ × . ◦ V band from UT 2001 March 11 until UT 2009 November29, obtaining 480 observations with a median per-pointuncertainty of 0.036 mag.2.3.4. ASAS-SN
Focused on the discovery and characterization of su-pernovae, the All-Sky Automated Survey for Super-Novae (ASAS-SN; Shappee et al. 2014; Kochanek et al.2017) surveys the entire sky down to V ∼
17 mag ev-ery ∼ ×
2k thinned CCD (Brown et al. 2013). The tele-scopes have a 4 . ◦ × . ◦ . (cid:48)(cid:48) ANALYSIS AND RESULTSThese new datasets allow us to study the architec-ture of the GW Ori system in a comprehensive manner.First, we forward-model the ALMA CO and C Otransitions to reconstruct the disk velocity field, mea-sure its inclination, and determine the total stellar massof GW Ori. Next, we fit a hierarchical triple star modelto the RVs and archival astrometry to determine individ-ual stellar masses and orbital inclinations, and compare
Table 3.
Photometric measurements ofGW Ori.HJD m V σ V Telescope[2,400,000+] [mag] [mag]47031.4670 9.94 · · ·
Maidanak47032.4760 9.90 · · ·
Maidanak47034.4826 9.86 · · ·
Maidanak47035.4806 9.88 · · ·
Maidanak47036.4839 9.87 · · ·
Maidanak
Note —This table is available in its entiretyin machine-readable form. these properties to that of the disk. Last, we use theextensive lightcurve of GW Ori to identify new eclipseevents and oscillatory modes, and compare these to theorbital periods to inform theories of their physical origin.3.1.
A Reconstruction of the Disk Velocity Field
We use the spatially and spectrally resolved molecularline emission observed with ALMA to tomographicallyreconstruct the disk velocity field and make a dynami-cal measurement of the total stellar mass in the GW Orisystem. We follow the forward modeling procedures de-scribed by Czekala et al. (2015a, 2016) using the associ-ated open-source software package
DiskJockey . The basis of the parametric physical model adoptedin this approach is a radial surface density profile, Σ( r ),designed to mimic a simple theoretical description fora viscous accretion disk (Lynden-Bell & Pringle 1974;Hartmann et al. 1998). That profile decreases like 1 /r interior to a characteristic radius R c , and has an expo-nential taper e − r/R c at larger radii. The vertical dis-tribution of density is controlled by the temperaturestructure: we assume a vertically isothermal model witha radial profile T ( r ) = T ( r/
10 AU) − q . To convertthe total gas densities to the volume density of a givenspecies, we use abundance ratios that are representa-tive of the dense interstellar medium: [H / gas] = 0.8,[H/H ] = 2, [ CO/H] = 7 . × − , [ CO/ CO] =69, and [ CO/C O] = 557 (e.g., Henkel et al. 1994;Prantzos et al. 1996).The disk kinematics are assumed to be Keplerianand dominated by the total stellar mass M tot , with Available under an MIT license at https://github.com/iancze/DiskJockey . a velocity field that appropriately accounts for thetwo-dimensional distribution of the emitting layer (seeRosenfeld et al. 2013). The line-spread function is char-acterized with a width defined by the quadrature sumof thermal and non-thermal ( ξ ; presumably turbulent)contributions. For any physical structure specified bythese 6 parameters, { Σ c , R c , T , q , M tot , ξ } , we solvethe molecular rate equations (assuming LTE) and ray-trace the associated emission into a set of high resolu-tion channel maps using the radiative transfer package RADMC-3D (Dullemond 2012). That ray-tracing requiresthat we specify 5 additional geometric parameters: thedisk inclination to the line-of-sight ( i disk ), the positionangle of the disk rotation axis projected on the sky( ϕ ), the LSRK systemic radial velocity ( v r ), and a pairof positional offsets from the observed pointing ( µ α , µ δ ). We adopt a fixed distance to the GW Ori system, d = 388 pc (Kounkel et al. 2017), to make the problemmore computationally tractable; the effects of this as-sumption are discussed in Section 3.3. The GAIA DR1parallax to GW Ori is still rather uncertain: the meanestimate is slightly larger than our adopted distancebut its large uncertainty means that it is still consistentat the 1 σ level ( π = 2 . ± . χ likelihood function thatincorporates the nominal visibility weights. We assumeflat priors on all parameters except for i disk , where in-stead we adopt a simple geometric prior (the disk an-gular momentum vector is distributed uniformly on asphere; e.g., Czekala et al. 2016). The posterior distribu-tion of these parameters is explored using Markov ChainMonte Carlo (MCMC) simulations with the affine invari-ant ensemble sampler proposed by Goodman & Weare(2010), as implemented in the emcee code (Foreman-Mackey et al. 2013) and ported to the Julia program-ming language (included in
DiskJockey ).Compared to previous similar work, the modeling ofGW Ori is considerably more computationally expen-sive. This is primarily a consequence of the large phys-ical size of the disk, which makes the ray-tracing stepsubstantially more time-consuming. The inference foran individual spectral line takes ∼ CO lineis clearly contaminated by local cloud material, we re-strict our analysis to independent inferences of the modelparameters based on the CO and C O datasets. Forexpediency, we only model the data averaged to 25 chan- d a t a − − . − . . . . −
20 0 20 m o d e l − r e s i d u a l − − − α [ ] − − ∆ δ [ ] Jy/beam × RMS
Figure 3.
A comparison of the observed channel maps of the CO line emission ( top ) with a best-fit model ( middle ; constructed from asynthetic visibility set based on the inferred parameters listed in Table 4 and then imaged in the same way as the data) and the associatedresiduals ( bottom ; the imaged data − model residual visibilities). The annotation is the same as in Fig. 2. d a t a − − . − .
02 0 .
00 0 .
02 0 . − − m o d e l − r e s i d u a l − − − α [ ] − − ∆ δ [ ] Jy/beam × RMS
Figure 4.
A comparison of the observed channel maps of the C O line emission ( top ) with a best-fit model ( middle ; constructed from asynthetic visibility set based on the inferred parameters listed in Table 4 and then imaged in the same way as the data) and the associatedresiduals ( bottom ; the imaged data − model residual visibilities). The annotation is the same as in Fig. 2. nels of 0 . − width. Experiments modeling a sub- set of the channels at higher resolution (e.g., using everythird 0.17 km s − -wide channel) yielded similar results.0 050100150 z [ A U ] CO C O r [AU] Figure 5.
The maximum likelihood 2D temperature (top) and density (bottom) disk structures inferred using the CO (left) and C O(right) transitions. The temperature contours are in units of K. The density plots show the total gas density ( ρ gas ) and are in unitslog cm − . The color scales are normalized to the same limits for both transitions. . . . . . . . M tot [ M (cid:12) ]131132133134135136137138139140 i d i s k [ ◦ ] COC O Figure 6.
Posterior distributions for the model parameters fitto the CO and C O data independently, showing 1 and 2 σ contours. Dashed lines indicate constant values of M tot sin i disk . The parameter values inferred from each spectral linedataset are summarized together in Table 4. A compar-ison of the data and the best-fit models (and associatedresiduals) is presented in the form of channel maps inFigures 3 and 4 for CO and C O, respectively. Whileoverall the models successfully reproduce the observedemission, there are some interesting residuals, namely,an excess of emission in the center of the disk for thechannels between 13.1–14.3 km s − , seen in both COand C O. We will return to a discussion of a potentialorigin of those residuals in Section 3.4.Motivated by the presence of the aforementionedresiduals, we explored more sophisticated disk models,including a model with a vertical temperature gradientand CO depletion due to freeze-out and photodissocia-tion (after Rosenfeld et al. 2013), as well as a flexibletemperature model parameterized to mimic more so-phisticated (and computationally expensive) protoplan-etary disk models (Kamp & Dullemond 2004; Jonkheidet al. 2004). However, we found that neither of thesemodels resulted in a more satisfactory fit to the dataas measured by visual inspection and the Akaike infor-mation criterion (AIC; Akaike 1973). Encouragingly,however, they still yielded similar estimates of M tot as the standard model, which gives us confidence that1 Table 4.
Inferred Disk Model ParametersParameter CO C O M tot [ M (cid:12) ] 5 . ± .
06 5 . ± . r c [au] 237 ± ± T [K] 51 ± ± q . ± .
012 0 . ± . M disk log [ M (cid:12) ] − . ± . − . ± . ξ [km s − ] 0 . ± .
01 0 . ± . i d [deg] 137 . ± . . ± . a [deg] 90 . ± . . ± . v r b [km s − ] +13 . ± .
003 +13 . ± . µ α [ (cid:48)(cid:48) yr − ] − . ± . − . ± . µ δ [ (cid:48)(cid:48) yr − ] − . ± . − . ± . a For comparison with the stellar orbits, we note that the positionangle of the ascending node Ω disk is 90 ◦ offset from our PAconvention, i.e. Ω disk = PA + 90 ◦ ≈ . ◦ . b LSRK reference frame. In the barycentric reference frame, thedisk systemic velocity is +28 .
34 km s − . Note —The 1D marginal posteriors are well-described by aGaussian, so we report symmetric error bars here. the disk-based dynamical mass is sufficiently robustto choice of parameterization for the temperature anddensity structures.While the inferred physical structures inferred fromeach line are in mild disagreement, as might be expectedfor this simple parameterization, the spatial values of theinferred temperature and density in the disk are for themost part quite similar. To illustrate this, we plot the2D temperature and density profiles inferred from eachtransition in Figure 5. We attribute the small differencesin the structure parameters to the different layers of thedisk probed by the CO and C O transitions: the moreoptically thick CO probes the upper layers of the diskatmosphere while the more optically thin C O tran-sition is more sensitive to the colder, denser, layers ofthe disk midplane. In Figure 6, we plot the marginalizedposteriors for both transitions in the { M tot , i disk } -plane.Interestingly, the different transitions deliver differentinclinations (∆ i disk = 2 . ± . ◦ ), which we attributeto the previously mentioned model deficiencies and thefact that the CO and C O transitions probe differentlayers in the disk. With more computational power, itwould be worthwhile to explore a joint fit to both transi- tions to see if a single, more sophisticated, disk structurecould adequately fit both transitions simultaneously.Nevertheless, both transitions yield consistent con-straints on the total stellar mass, which is the most rel-evant parameter to our stated goals. The robustness ofthe dynamical mass technique is primarily because thekinematic morphology of the line emission (i.e., the dis-tribution of the emission in position–velocity space) isnot strongly dependent on the temperature and densitystructure of the disk, but is a rather strong function of M tot and i disk . When the disk is spatially resolved, thedependence of M tot on i disk is considerably diminished.We combine the inferred total masses from CO andC O, weighted by their uncertainties, to find M tot =5 . ± . M (cid:12) . The uncertainty in the distance toGW Ori (388 ± M tot = 5 . ± . M (cid:12) , which we report as thetotal mass estimate. Because the inferred disk inclina-tions are mutually inconsistent, we adopt a weightedaverage for the mean inclination and assume a largesystematic uncertainty, resulting in a final estimate of i disk = 137 . ± . ◦ . Our CO results are broadly con-sistent with that determined by Fang et al. (2017), whomeasure the disk inclination to be ∼ ◦ (modulo theabsolute inclination of the disk).3.2. An Updated Model of the Stellar Orbits
In this section, we present an orbital fit to the RVsdetermined in § v TRES , applied as a shift to all TRES RVsto place them on the DS reference frame; (2) ∆ v DS2 , to al-low for an offset between the primary and secondary DSvelocities, possibly caused by a mismatch between thetemplate parameters and those of the true stars; and (3)∆ v TRES2 , a similar primary/secondary offset for TRES.The residuals from our initial fit indicated that our for-mal velocity uncertainties are underestimated, and sothe uncertainties on each measurement are scaled toachieve a reduced χ ν = 1 for our final solution.The period of the inner orbit is consistent with thatof Mathieu et al. (1991) and Fang et al. (2014); how-ever, due to the SB2 nature of the system, most otherorbital parameters are significantly different. We finda larger semi-amplitude for the primary, K A = 8 . ± . . . V [ m a g ] MaidanakKELT-EastKELT-WestASASASAS-SN v r [ k m s − ] − [ k m s − ] date [JD] − [ k m s − ] Figure 7. top : photometric observations of GW Ori from 1987 until mid 2017. All photometric observations displayed here are in the V -band (ASAS, ASAS-SN, and Maidanak) or a broader filter (KELT) which has been shifted to align with V -band where the time seriesoverlap. bottom : primary (blue) and secondary (orange) radial velocities overlaid with several realizations of the most probable orbits,to show uncertainty in the orbit. Reticon velocities are shown with squares, TRES velocities are shown with circles, and the dotted linerepresents the center-of-mass velocity. Residuals for this orbit are shown in the panels below. .
14 km s − , a mass ratio of q ≡ M B /M A = 0 . ± .
02, and a statistically significant eccentricity e in =0 . ± .
02. The outer orbit has an orbital period of P out = 4218 ±
60 days (11.5 years), and a significant ec-centricity, e out = 0 . ± .
09. We find offset termsstatistically inconsistent with zero: a small but non-negligible offset between the DS and TRES zeropointsof 0 .
49 km s − and larger offsets for the secondary ve-locities of 8 .
77 km s − and 6 .
41 km s − , for the DSand TRES RVs, respectively. Given the large intrin-sic linewidth ( v sin i ≈
40 km s − ), these large offsetscan reasonably be ascribed to template mismatch. The systemic velocity inferred from the RV fit is nicely con-sistent with the systemic velocity of the circumtriple diskfrom the ALMA data. All parameters of the RV fit arelisted in the first column of Table 5. The full orbit asa function of time is shown in the second panel of Fig-ure 7. Graphical representations of our observations andthe inner and outer orbit models as a function of orbitalphase are shown in Figure 8 and Figure 9, respectively.Although there are only three epochs of published as-trometry in Berger et al. (2011), these points may stillhelp constrain the parameter space of possible orbits.Therefore, we explore a joint RV-astrometric analysis39 . . . V [ m a g ] − − v , i n − γ [ k m s − ] − [ k m s − ] − − v , i n − γ [ k m s − ] − . . . . . phase from inner periastron − − [ k m s − ] Figure 8. top : V -band light curve phased to the inner orbitalperiod. No discernible correlation is detected. The colors are thesame as in the top panel of Figure 7. bottom : RV measurementsof GW Ori and best-fit model for the inner orbit, after subtractingthe motion due to the outer orbit. built upon a model of the “three-dimensional orbit” fol-lowing Murray & Correia (2010), which adds new modelparameters including the semi-major axis, orbital incli-nation, and position angle of the ascending node for boththe inner and outer orbits. For a likelihood function, wecombine the χ RV likelihood and a new χ likelihoodfor the angular separation and position angle measure-ments of the B and C components relative to A. As withthe disk analysis, we also use a geometric prior on theorbital inclinations. For their last measurement epoch(2005), Berger et al. (2011) report an alternate position 9 . . . V [ m a g ] − v , o u t − γ [ k m s − ] − [ k m s − ] − − v , o u t − γ [ k m s − ] − . . . . . phase from outer periastron − − [ k m s − ] Figure 9. top : V -band light curve phased to the outer orbitalperiod, showing that the mean flux level oscillates by 0.2 magover the course of the outer orbit, and is lowest during apoastron(phase 0.5). The colors are the same as in the top panel of Figure 7. bottom : RV measurements of GW Ori and best-fit model for theouter orbit, after subtracting the motion due to the inner orbit. for the C component, and so we also perform a separatefit for this scenario.The jointly-constrained parameters are in the secondand third columns of Table 5. A graphical representa-tion of the orbit is shown in Figure 10. With the addi-tion of the astrometric dataset, we can measure the indi-vidual stellar masses and the inclinations of the orbitalplanes, which are also listed in Table 5. Depending onwhether the original or alternate position for C is used,we find the total stellar mass to be M tot = 5 . ± . M (cid:12) . ± . M (cid:12) , respectively. Both measurements areconsistent with the M tot independently measured withthe disk-based analysis ( M tot = 5 . ± . M (cid:12) ).To measure the degree of misalignment between theorbital planes and the circumtriple disk, we calculate theangle Φ between the angular momentum vectors of eachorbit according to Fekel (1981)cos Φ = cos i cos i + sin i sin i cos(Ω − Ω ) . (1)We find that the mutual inclination between the diskand the inner orbit is Φ in = 44 ± ◦ and the mu-tual inclination between the disk and the outer orbitis Φ out = 54 ± ◦ ; these values are similar if one uses the“alternate” C position (Φ in = 45 ± ◦ , Φ out = 50 ± ◦ ).Such a large misalignment is surprising given the naiveexpectation that the stellar orbits and disk would beroughly co-planar. Since these results only rest uponthree astrometric epochs, however, there is a possibilitythat the large inferred mutual inclinations may be theresult of unaccounted for systematic effects. In the nextsection, we use only the newly derived RVs and disk-based dynamical mass to formulate a more conservativeestimate of the mutual inclinations. We advocate con-tinued astrometric monitoring of the GW Ori system tofurther improve the three-dimensional orbit and defini-tively confirm the inclinations of the stellar orbits.3.3. Joint RV + Disk Constraints on IndividualComponent Masses
In this section, we combine the RV analysis with thedisk-based constraints on the total stellar mass to in-dependently infer the individual stellar masses of theGW Ori system without referencing the Berger et al.(2011) astrometry. We construct a joint likelihood func-tion with the following five parameters: M A , M B , M C , i in , and i out . The RV constraints are sufficiently cap-tured by the summary statistics M A sin i in , M B sin i in , M C sin i out / ( M tot /M (cid:12) ) / and the covariances betweenthem, which are well represented by a multivariateGaussian distribution. The disk constraint on the totalstellar mass M tot is well-represented by a Gaussian, aswell. We enforce flat priors on the stellar masses and ge-ometrical priors on the inclinations. We use the ensem-ble sampler MCMC (Goodman & Weare 2010; Foreman-Mackey et al. 2013) with 20 walkers to explore the poste-rior for 50,000 iterations, burn 25,000 iterations, and as-sess convergence by ensuring the Gelman-Rubin statistic(Gelman et al. 2014) is ˆ R < . Note that we do not use additional constraints on q in or otherderived parameters, as this would amount to double-counting theRV constraints. This analysis produces consistent but less precise con-straints on the stellar masses as the joint RV + astro-metric fits (see Table 6). Like the RV + astrometricanalysis, the disk + RV analysis also indicates that theinner stellar orbit must be significantly misaligned withthe disk, although the true misalignment is unknown:the difference between i in and i disk only provides a lowerlimit on the mutual inclination because the true mutualinclination must consider the position angles of the or-bits as well. To highlight these findings, we overplot ournewly derived constraints on the disk inclination andthe orbits in Figure 11. The measurements for the diskbased on the CO and C O data indicate the lowestinclinations (nearest to edge-on, i = 90 ◦ ). This is com-mensurate with the disk inclination measurements fromFang et al. (2017), and so we consider these results tobe robust. Interestingly, the constraints on i in differ be-tween the RV + astrometry and the RV + disk resultsat a significant level. We speculate that this differencemight be due to unknown systematics in the astrometryor RV datasets, or potentially an error in our assump-tion of the distance to GW Ori. While the astrometry+ RV analysis does not require a distance to the source,the disk-based analysis does require a distance in orderto break the M tot /d degeneracy. Although the exact de-gree of mutual inclination between the inner orbit andthe disk is unknown, we conclude that it is at least 10 ◦ and potentially as high as 45 ◦ .3.4. Variability: Eclipses and Periodic Behavior
The first extensive lightcurve of GW Ori was publishedby Shevchenko et al. (1992), based on several seasons ofphotoelectric photometry from Maidanak Observatory.Among other modulations typical of young stars, thosedata revealed two deep (∆ V = 0 . i in ≈ ◦ ).Maidanak Observatory continued their extensive pho-tometric monitoring of GW Ori until 1996 Shevchenkoet al. (1998). Those data identified two additionaleclipses, one in 1991, and one in 1992 (Figure 12, Gand H), which occurred exactly three orbital periods af-ter the 1990 event, but with a lower amplitude (∆ V =0 . V − R and B − V colors) during5 Table 5.
Orbital elements of GW Ori.Parameter RV RV + astrometry RV + astrometry † Inner orbit P [days] . . . . . . . . . . . . . . . . . . . . 241.49 ± ± ± K A [km s − ] . . . . . . . . . . . . . . . . 8.36 ± ± ± q . . . . . . . . . . . . . . . . . . . . . . . . . . 0.60 ± ± ± a [au] . . . . . . . . . . . . . . . . . . . . . . · · · ± ± e . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.13 ± ± ± i [deg] . . . . . . . . . . . . . . . . . . . . . · · · ± ± ω A [deg] . . . . . . . . . . . . . . . . . . . . 196 ± ± ± b [deg] . . . . . . . . . . . . . . . . . . . . · · · ±
13 264 ± T peri [HJD − ± ± ± γ [km s − ] . . . . . . . . . . . . . . . . . . +28.31 ± ± ± v TRES a [km s − ] . . . . . . . . 0.49 ± ± ± v Reticon a [km s − ] . . . . . . 8.77 ± ± ± v TRES a [km s − ] . . . . . . . 6.41 ± ± ± M A [ M (cid:12) ] . . . . . . . . . . . . . . . . . . · · · . +0 . − . . +0 . − . M B [ M (cid:12) ] . . . . . . . . . . . . . . . . . . · · · . +0 . − . . +0 . − . Outer orbit P [days] . . . . . . . . . . . . . . . . . . . . 4218 ±
60 4246 ±
66 4203 ± K AB [km s − ] . . . . . . . . . . . . . . . 2.47 ± ± ± a [au] . . . . . . . . . . . . . . . . . . . . . . · · · ± ± e . . . . . . . . . . . . . . . . . . . . . . . . . . . 0.22 ± ± ± i [deg] . . . . . . . . . . . . . . . . . . . . . · · · ± ± ω AB [deg] . . . . . . . . . . . . . . . . . . . 307 ±
18 310 ±
21 310 ± b [deg] . . . . . . . . . . . . . . . . . . . . · · · ± ± T peri [HJD − ±
565 53911 ±
260 53878 ± M C [ M (cid:12) ] . . . . . . . . . . . . . . . . . . · · · . +0 . − . . +0 . − . Derived propertiesInner time interval [cycles] . . 53.6 · · · · · ·
Outer time interval [cycles] . 3.1 · · · · · · M A sin i [ M (cid:12) ]. . . . . . . . . . . . . . . 0.30 ± · · · · · · M B sin i [ M (cid:12) ] . . . . . . . . . . . . . . . 0.18 ± · · · · · · M C sin i/ ( M tot /M (cid:12) ) / [ M (cid:12) ] 0.22 ± · · · · · · M tot [ M (cid:12) ] . . . . . . . . . . . . . . . . . · · · . ± . . ± . a We include parameters for a potential velocity offset between the primary and secondaryradial velocities for each instrument. In principle, this term should be consistent with 0; thenon-zero value likely indicates that there is some moderate template mismatch between thesecondary stellar spectrum and the synthetic spectrum used as a cross correlation template. b We follow the convention of the visual binary field and define the ascending node as thepoint where the secondary component (e.g., star B for the inner orbit, and star C for theouter orbit) crosses the plane of the sky moving away from the observer. † Fit using the “alternate” C position for the 2005 epoch of astrometry in Berger et al. (2011).6
Outer time interval [cycles] . 3.1 · · · · · · M A sin i [ M (cid:12) ]. . . . . . . . . . . . . . . 0.30 ± · · · · · · M B sin i [ M (cid:12) ] . . . . . . . . . . . . . . . 0.18 ± · · · · · · M C sin i/ ( M tot /M (cid:12) ) / [ M (cid:12) ] 0.22 ± · · · · · · M tot [ M (cid:12) ] . . . . . . . . . . . . . . . . . · · · . ± . . ± . a We include parameters for a potential velocity offset between the primary and secondaryradial velocities for each instrument. In principle, this term should be consistent with 0; thenon-zero value likely indicates that there is some moderate template mismatch between thesecondary stellar spectrum and the synthetic spectrum used as a cross correlation template. b We follow the convention of the visual binary field and define the ascending node as thepoint where the secondary component (e.g., star B for the inner orbit, and star C for theouter orbit) crosses the plane of the sky moving away from the observer. † Fit using the “alternate” C position for the 2005 epoch of astrometry in Berger et al. (2011).6 − − − α cos δ [mas] − − ∆ δ [ m a s ] − − α cos δ [mas] − − − ∆ δ [ m a s ] n e a r s i d e − −
10 0 10∆ α cos δ [mas] − − − − ∆ Z [ m a s ] − δ [mas] − − − ∆ Z [ m a s ] Figure 10.
Orbits from the joint RV-astrometric fit shaded according to their phase, where black represents periastron and color hueincreases with orbital phase. top left : orbits relative to the primary star A, in the plane of the sky, with the three epochs of astrometryfrom Berger et al. (2011). The light grey data point and outer orbit represent the fit to the “alternate” position for C. The following threeplots are relative to the center of mass of the system. A fictitious particle at 7 au on a circular orbit coplanar with the circumtriple diskis shown as grey dashed line. top right : the sky plane. For future discussion in § bottom left : looking down the North axis. bottom right : looking down the East axis. Positive Z points towards observer. Table 6.
Joint constraints on stellar massesand orbital inclinationsParameter RV + astrometry RV + disk M A [ M (cid:12) ] 2 . +0 . − . . +0 . − . M B [ M (cid:12) ] 1 . +0 . − . . +0 . − . M C [ M (cid:12) ] 1 . +0 . − . . +0 . − . i in [deg] 157 +1 − +1 − i out [deg] 150 +7 − +28 − Note —The RV + astrometry values are repli-cated from Table 5 for comparison purposes.We note that we are not able to infer theabsolute inclination of the stellar orbits di-rectly from the radial velocity data, so thereare in fact alternate solutions for the RV +disk results that yield i alt = 180 ◦ − i . Thesesolutions would be inconsistent with the as-trometric motion, however, so we opt to onlyreport the solutions with i ≥ ◦ .
130 140 150 160 i [ ◦ ]0 . . . . . . . p ( i | D ) [ d e g − ] prior COC O i in RV + disk i in RV + ast i out RV + ast
Figure 11.
The inclination posteriors on the disk inclination,inner stellar orbit, and outer stellar orbit, as determined fromvarious joint fits. Because i out is essentially unconstrained by theRV + disk analysis, it is not plotted for aesthetic reasons. Thegeometric prior on inclination (uniform orientation of orbits in 3Dspace) is shown as a thin grey dotted line. all of these eclipses; however, insufficient precision wasavailable to robustly constrain an associated extinctioncurve. Photoelectric observations by W. Herbst revealthree additional eclipses of similar depth between 1982and 1985 (events A - C; Shevchenko et al. 1998). Thosefeatures also appear to be separated by integer multi-ples of the inner orbital period, lending further support to the Algol-like variable hypothesis. With their longerphotometric time baseline, Shevchenko et al. (1998) alsonoted an overall decline of ∆ V = 0 . V -band: it was shifted to match the overlappingASAS-SN observations, meaning that it has a poten-tially problematic zeropoint uncertainty. With thesecaveats in mind, the clear rising and falling trends areseen within each individual dataset without the need8 Table 7.
V-band Photometric Eclipse CatalogLabel UT Mid Start End Duration Depth TelescopeJD JD days magA 1982 Oct · · · · · · · · · · · ·
HerbstB 1983 May · · · · · · · · · · · ·
HerbstC 1984 Nov · · · · · · · · · · · ·
HerbstD 1988 Sep 13 2447418 2447435 17 0.40 MaidanakE 1989 Oct 6 2447806 2447833 27 0.10 MaidanakF 1990 Sep 5 2448140 2448161 21 0.35 MaidanakG 1991 Oct 13 2448543 2448589 46 0.15 MaidanakH 1992 Sep 7 2448873 2448893 20 0.08 MaidanakI 2001 Oct 1 2452184 ≥ ≥
31 0.70 MaidanakJ 2007 Dec 15 2454450 2454508 58 0.10 ASASK 2009 Jan 30 2454862 2454912 50 0.15 ASASL 2012 Oct 15 2456216 2456259 43 0.11 KELTM 2014 Feb 16 2456705 2456744 39 0.20 KELTN 2014 Nov 5 2456967 2457080 30-130 0.25 KELT/ASAS-SNO 2015 Oct 14 2457310 2457405 95 0.10 ASAS-SN for any vertical shifts, suggesting that this modulationis likely real. The earlier photoelectric observations ofW. Herbst, stretching back to 1983, also clearly phaseup with the expected sinusoidal variation. This situatesthe long term dimming seen by Shevchenko et al. (1998)as part of an 11.5 yr period brightness oscillation.When considering the phase-folded lightcurve on theouter (AB–C) orbital period, it appears as if the deepeclipses preferentially occur between phases 0.4–0.8 fromperiastron. However, the three eclipses in the photoelec-tric observations of W. Herbst fall closer to 0.0 phase.Taken together, this suggests that the apparent cluster-ing of deep eclipses could be affected by the seasonalsampling in the dataset.The 11.5 yr variability reaches a minimum flux levelnear apoastron (phase = 0.5) of the outer orbit. Sincethe light from the A–B binary dominates the total opti-cal flux from the system, it must be one or both of thesestars that are either being partially occulted by circum-stellar material. Due to the fact that the disk is moreinclined than the stellar orbits, apoastron correspondsto the time when the A–B pair comes closest to beingscreened by the inner edge of the near side of the disk(see Figure 10, top right panel). From dynamical argu-ments, we expect the inner edge of the disk to be trun-cated out to 2–3 times the semi-major axis of the ter- tiary (Artymowicz & Lubow 1994), which correspondsto ∼ (cid:38)
10 au. Given the gradual dimming, it maybe more likely that the A–B pair is screened by tenuousmaterial residing inside this truncation radius, such asmicron-sized dust within the cleared region. There iscircumstantial support for this interpretation from thevariable infrared SED, which Fang et al. (2014) interpretas an indication of a variable reservoir of small grainsnear the A–B pairing which is cleared and replenisheddue to the actions of the tertiary.We find some CO and C O emission that exceedspredictions from the most probable standard disk modelat locations consistent with this near edge of the disk,but located at or near the systemic velocity (13.1–14.3 km s − LSRK). Those residuals could simply be anartifact from using an insufficiently complex disk model,but they may instead very well be probing the source ofthe eclipses, the long term dimming, or both. There isan outstanding question from this analysis as to whatthe disk looks like on the scales of the inner orbit. TheALMA observations do not have sufficient angular res-olution to probe the disk at the physical scales corre-9 −
25 0 2510 . . V [ m a g ] . . . . . . . . . . . . . . . . days from eclipse start . . . . . . . . . . Figure 12.
A gallery of the eclipse events noted in Table 7, labeled relative to the start of the eclipse. The colors are the sameas in the top panel of Figure 7. Note that the y-axis scale is significantly different from panel to panel, with some eclipses asdeep as 0.7 mag and others less than 0.1 mag. sponding to the tertiary orbit, so the distribution ofsolids within the GW Ori disk remains relatively uncon-strained. With its longest baselines, ALMA would havethe spatial resolution (0 . (cid:48)(cid:48)
02) to probe the disk on 8 auscales, more than sufficient to resolve a cleared regionconsistent with the tertiary orbit ( (cid:38)
18 au in diameter).Because the eclipse events are not synchronized withthe outer orbital period, it is not immediately clearwhether they share the same physical origin as thelonger-term brightness variations. The new dynamicalconstraints derived earlier indicate that the A–B orbitis not edge-on, and so we must consider alternatives tothe Algol mechanism. Any theory that seeks to explainthese inner eclipses must account for several pieces ofevidence, in addition to the updated orbital configura-tion. The eclipses span 10–50 days in duration, are ofvariable depth (between 0.08 and 0.70 magnitudes), ap-pear to be consistent with reddening by dust, are notstrictly periodic with the inner period, and seem to oc-cur at all phases of the outer orbital period. Moreover,the spectra (unwittingly) taken during times of eclipseshow no obvious changes in spectral features beyond thenormal variability described by Fang et al. (2014). Wespeculate that these quasi-periodic eclipses may be dueto an unstable circumbinary disk around A–B, or possi- bly the result of eclipses of A by accretion streams ontoeither the individual circumstellar disks of A or B.Finally, we searched for additional periodicty in thelightcurve beyond the inner and outer orbital periods.After excluding eclipses, we used the Lomb-Scargle (LS)periodicity search algorithm (Lomb 1976; Scargle 1982)(within the VARTOOLS analysis package; Hartman &Bakos 2016) to search for periodic modulations from1.1 to 100 days in the high-cadence KELT dataset. Themost significant period we recover is 2 . ± .
05 days,shown in Figure 13. We suggest that this correspondsto the rotation period of the primary, although we notethat Bouvier (1990) and Fang et al. (2014) derived al-ternate rotation periods of 3.3 days and 5.0–6.7 days,respectively. Future high-precision, high-cadence obser-vations of GW Ori will help to unambiguously identifythe rotation period of A. DISCUSSIONWith the newly derived component masses now estab-lished, we turn to discussion of the photospheric prop-erties of the stars and the age of the GW Ori system.With these in place, we discuss the system architecturein the context of other multiple systems.4.1.
Age and Photospheric Properties − . . . . . phase − . − . . . . r e l a t i v e m a g n i t u d e Figure 13.
The KELT photometric observations of GW Ori,with the three eclipses removed, phased to the 2.93 day pe-riod recovered from our LS analysis. Legend as in Figure 7
In order to place GW Ori A and GW Ori B on theHR diagram, we require updated measurements oftheir luminosities. To obtain those, we assembled anSED of GW Ori from the same sources listed in Fanget al. (2014), i.e., the
U BV R C I C photometry from Cal-vet et al. (2004) and the JHK s photometry from the2MASS survey (Skrutskie et al. 2006). We then man-ually adjusted the 2MASS J and H fluxes down by5% and 10%, respectively, to account for the approxi-mate contamination in those bands from star C basedon the flux ratios from Berger et al. (2011). We fittedthe SED using a two-component model based upon theNextGen atmospheres (Hauschildt et al. 1999) and theCardelli et al. (1989) reddening law with the followingconstraints: the distance to the source (388 ± . ± . T eff (5700 ±
200 K), andthe H -band flux ratio ( f B /f A = 0 . ± .
05; Bergeret al. 2011). That analysis yields an extinction of A V = 1 . ± . T eff = 4900 ±
200 K, lu-minosities of L A = 32 . ± . L (cid:12) and L B = 12 . ± . L (cid:12) ,radii of R A = 5 . ± . R (cid:12) and R B = 5 . ± . R (cid:12) ,and a V -band flux ratio of 0 . ± . . M (cid:12) and increase in incre-ments of 0 . M (cid:12) ; the isochrones are at ages of 0.1, 0.5,1.0, 2.0, and 3.0 Myr. The positions of A and B in thisdiagram are consistent with their measured dynamicalmasses at ages of 0.3–1.3 Myr. However, the evolution-ary tracks would imply that the A component is olderthan the B component by at least 0.3 Myr. This agediscrepancy is likely not real, but rather the result of in-naccuracies in the photospheric properties, evolutionarymodels, or both. The revised mass for GW Ori A is sig- nificantly lower than previous estimates in the literature(Berger et al. 2011; Fang et al. 2014), which, combinedwith the revised photospheric properties, means that theage of GW Ori is now slightly older than was previouslyimplied. Even with its lower mass and luminosity, how-ever, it is still true that GW Ori A is a Herbig Ae/Beprecursor (as was also noted by Fang et al. 2014), andwill evolve to a main sequence late-B/early-A star.In Figure 15, we utilize the three MIST models near-est the best-fit masses for each of the stellar componentsto compute the evolution of the effective temperaturesand the V - and H -band flux ratios, and compare thesequantities to the existing photospheric measurements.Beyond the primary and secondary effective tempera-tures and V -band flux ratio reported in this work, thereare two H -band flux ratio constraints from Berger et al.(2011), f B /f A = 0 . ± .
05 and f C /f A = 0 . ± . V -band flux contribution for the C component is verysmall ( (cid:46)
5% of the total flux), explaining why we wereunable to find optical spectroscopic signatures of thiscomponent even though it is nearly a solar mass star.As noted by Berger et al. (2011), it is possible thatone, two, or all three of the stars might show excess H -band emission due to the presence of a circumstellardisk and/or accretion signatures above the photosphericemission of these stars. While we find that the mea-sured photospheric properties are reasonably consistentwithout invoking such an excess, there is circumstantialevidence for the idea (beyond the eclipses). Najita et al.(2003) found blended CO fundamental emission peaks,which Bast et al. (2011) suggested might actually bethe blended profiles of emission originating from physi-cally distinct regions, e.g., individual circumstellar disksaround the A and B stars. However, it may also be pos-sible that the CO fundamental emission originates fromthe inner edge of the circumtriple disk, cleared by the ∼ v sin i A ) to infer its stellar obliquity, mod-ulo its absolute orientation. Solving for the inclinationand propagating the uncertainties in the respective pa-rameters yields an obliquity of i A = 23 ± ◦ , which isin remarkable agreement (modulo the absolute orien-tation) with the inclination of the inner stellar orbit1(157 ± ◦ ) as determined from the RV + astrometryanalysis (Sect. 3.3). 400050006000 T eff [ K ]10 L ∗ [ L (cid:12) ] . M (cid:12) . M (cid:12) . M (cid:12) . . Figure 14.
GW Ori A and B placed on the pre-main sequenceHR diagram, with evolutionary tracks from Choi et al. (2016).Mass tracks are in increments of 0 . M (cid:12) from 1 . M (cid:12) to 4 . M (cid:12) ,and isochrones label 0.1, 0.5, 1, 2 Myr ages. The GW Ori Triple System in Context
The many unique datasets presented in this paperhave enabled us to paint a detailed picture of theGW Ori system. It is young (0.3–1.3 Myr), containsa considerable amount of stellar mass ( M tot = 5 . ± . M (cid:12) ), and hosts a massive disk ( M disk ≈ . M (cid:12) ),which makes it an extremely interesting system to studyin the context of theories about star and planet forma-tion, migration, and stability.We estimated the total disk mass of GW Ori using theresults from our CO and C O modeling in § CO andH . We find somewhat larger disk masses when model-ing C O compared to CO (0 . M (cid:12) vs. 0 . M (cid:12) ,respectively), which is in conflict with the finding ofFang et al. (2017) that C O must be depleted rela-tive to CO. We attribute the differences in our diskmasses to insufficiently complex models of disk struc-ture and optical depth effects, and note that in generalestimating disk masses from CO is notoriously difficult(Yu et al. 2017), although in our case it is encouragingthat they are roughly consistent with estimates based onthe dust continuum emission (0 . M (cid:12) ; Fang et al. 2017).In the context of the large disk mass survey by Andrewset al. (2013), GW Ori’s disk mass is slightly larger thanthe mean predicted value for its stellar mass, althoughstill consistent with the large 1 σ envelope in this rela-tionship at high stellar masses. In light of this large disk 45005000550060006500 T e ff[ K ] M A = 2 . M (cid:12) M B = 1 . M (cid:12) M C = 0 . M (cid:12) . . . V r a t i o f B /f A f C /f A . . . . τ ∗ [Myr] . . . H r a t i o Figure 15.
Relative photospheric properties of the GW Ori stel-lar components as a function of age, using the MIST pre-main se-quence evolutionary models and assuming the stars are coeval.The measured effective temperatures for the primary and sec-ondary ( T eff = 5700 ±
200 K and T eff = 4900 ± V -band flux ratio ( f B /f A = 0 . ± . H -band fluxratios ( f B /f A = 0 . ± .
05 and f C /f A = 0 . ± .
01; Berger et al.2011) are shown here as dotted lines and are all roughly consistentwith the model predictions for ages of 0.3 - 1.3 Myr. mass, we investigate whether the disk is Toomre stabletoday. We use the more massive M disk values from theC O results to derive a lower bound on Toomre’s Qparameter Q ( r ) = c s Ω πG Σ = (cid:115) k B M tot π µm H G (cid:115) T ( r ) r Σ ( r ) . (2)For the range of disk parameters determined from ourCO fitting, the minimum value is Q ≈
100 at r ∼
300 au,which means that the disk is not currently undergoinga global gravitational instability ( Q ≈ ◦ . Thesefindings agree well with the low mutual inclinations2found for Kepler circumbinary planets (Winn & Fab-rycky 2015), and therefore have exciting implicationsfor a large circumbinary planet occurrence rate (Li et al.2016). Of course, severe selection effects are at work inboth samples, and so care must be taken when extrap-olating these results to the population at large. Withinthis context, GW Ori stands apart due to the fact thatits stellar orbit is not aligned with the disk. This maybe a consequence of the larger stellar masses involved,the longer orbital period(s), and/or the existence of aclose tertiary. Here we examine several other longer pe-riod systems that are characterized by their significantlynon-zero mutual inclinations, and show that the largemutual inclinations found in GW Ori are not as uniquewhen considering other longer period systems.KH 15D is an eccentric ( e = 0 .
6) binary system with aslightly longer period than the aforementioned systems(48 days), and hosts a circumbinary dust ring misalignedby 10–20 ◦ from the stellar orbit (Chiang & Murray-Clay2004; Capelo et al. 2012). The eccentric stellar orbit anddisk misalignment causes dramatic photometric eclipseevents as stars are screened by the edge of the dustdisk. The eclipses also come and go as the ring pre-cesses about the binary. The edge of the occulting diskmust be sharp, suggesting that the ring is confined to anarrow region by a planet at 4 au. The disk and stellarorbital misalignment may be driven by dynamical inter-actions between the eccentric binary and the disk (e.g.,Martin & Lubow 2017; Zanazzi & Lai 2018).Moving to still larger orbital separations, the fre-quency of disk-stellar orbital alignment becomes lessclear, mainly due to incomplete orbital coverage. Con-sider GG Tau A, which is a triple system with circum-stellar disks around each of three components, as wellas a larger circumtriple disk. The circumtriple disk iscomposed of a dense ring containing 80% of the massand an outer gas disk extending up to 800 au, and isof similar mass to that of the GW Ori disk (0 . M (cid:12) ;Guilloteau et al. 1999). The stellar architecture ofGG Tau A is rather different than GW Ori, however.The primary star Aa has a mass of 0 . M (cid:12) , and is sit-uated in an “outer” orbit with another binary, Ab1–Ab2, which together have a combined mass of less than0 . M (cid:12) (Dutrey et al. 2016). The orbital elements ofthe triple system still have some uncertainty. Nelson& Marzari (2016) make a dynamical argument that theouter orbit has a semi-major axis of 62 au and is likelycoplanar with the outer circumtriple ring; on the otherhand, Cazzoletti et al. (2017) argue from disk dynamicsthat the disk and binary planes are misaligned by 20-30 ◦ . Further astrometric observations are required todefinitively characterize this system. Recently, the transition disk system HD 142527 ( M ∗ =2 . M (cid:12) , M disk ∼ . M (cid:12) ) was discovered to have a Mdwarf companion orbiting inside its large disk cavity(Biller et al. 2012; Lacour et al. 2016). The presenceof this companion provides a possible explanation forwhy a smaller, inner disk in this system appears to behighly misaligned ( ∼ ◦ ; Avenhaus et al. 2014; Marinoet al. 2015), which may be driven by secular precessionresonance between the disk and the companion (Owen &Lai 2017). The extreme mass ratio between the M dwarfand the primary ( q = 0 . r ∼
500 au, M disk ≈ . M (cid:12) ; Eisner2012; Takakuwa et al. 2017) outside of two circumstel-lar disks. The binary is on a wide orbit ( a = 70 au), towhich the circumstellar disks are misaligned by up to25 ◦ . The relative proper motion of the binary over a15 yr baseline indicates that it contains 1 . M (cid:12) of stellarmass on a 246 yr orbit (Villa et al. 2017). The circum-stellar disks are probably aligned with the rotation ofthe outer circumbinary envelope, which may indicatethat the stars were formed by rotationally-driven frag-mentation, preserving this orientation (Lim et al. 2016).The most recent example of a misaligned circumbinarydisk is in the TWA 3A system (Kellogg et al. 2017),which hosts a circumbinary disk within a hierarchicaltriple system of stars of near equal mass (SpTs M3–M4). The “inner” binary Aa-Ab has a 35 day eccentric( e = 0 .
63) orbit and hosts a small disk extending 25 au insize (Andrews et al. 2010), while the “outer” orbit A-Btakes 200–800 years. Although the absolute inclinationsare not yet known, the parameter space is sufficientlyconstrained such that it is likely that all three planes(the inner orbit, outer orbit, and cirumbinary disk) aremisaligned by at least 30 ◦ . The Aa-Ab circumbinarydisk mutual inclination may be attributable to torquesfrom the distant B companion.Of all these T Tauri sources surveyed, GW Ori standsout in terms of stellar mass. Its architecture proved rel-atively difficult to probe via traditional detection tech-niques, requiring sustained, long-term radial velocitymonitoring over 35 years, as well as sophisticated careand attention to derive a radial velocity solution for theblended line profiles. Finally, resolved submm inter-3ferometric observations were necessary to measure theinclination of the circumtriple disk. Looking forward,sustained campaigns of astrometric monitoring will bemost helpful for definitively constraining the orbital in-clinations.Although there may still be significant dynamical evo-lution in the architecture of the GW Ori system before itreaches the main sequence (MS), it is also worth brieflyconsidering how its orbital parameters compare to thegeneral MS population of triple stars, even though thispopulation of older systems may have experienced sig-nificant dynamical evolution. Among late B star pri-maries (which GW Ori A will be on the MS), 13% ofsystems have multiplicity of three or higher (Eggleton &Tokovinin 2008); this fraction is roughly constant acrossspectral types B–G. In a detailed analysis of higher-ordermultiple systems, Tokovinin (1997) found that the ratio P long /P short was greater than 10 in almost all systems,presumably reflecting which orbits are stable. WhileGW Ori’s ratio (17) is smaller than most, it is not anoutlier among triple systems. Tokovinin (1997, 2017)find that the distribution of mutual inclinations betweenthe orbital planes in triple systems is inconsistent bothwith complete alignment of the inner and outer orbits(zero mutual inclination) as well as completely inde-pendent inclinations (randomly distributed). For tripleswith outer projected separation <
50 au the average mis-alignment is 20 ◦ , while orbits wider than 1000 au arenot preferentially aligned. The RV + astrometric fitsfor GW Ori suggest that the stellar orbits are consistentwith this picture (a mutual inclintion of 13 ± ◦ ). Thepopulation of misaligned triples may be the result of ac-cretion of gas with randomly aligned angular momentumat the epoch of star formation.Like most multiple stars, those in GW Ori likelyformed through turbulent fragmentation of the molec-ular cloud, possibly at larger separations than they arenow, and then hardened through decay via dynam-ical interactions, accretion, and the interaction of thecircumstellar disks (Offner et al. 2010; Bate 2012). Con-tinued study of the GW Ori system, including spatiallyresolving the innermost regions to discover circumstellardisks and their relative inclinations, will be valuable tofurther understanding its formation process. SUMMARY AND CONCLUSIONSIn this work we: • Used spatially and spectrally resolved ALMA ob-servations of the GW Ori circum-triple disk to de-rive a dynamical mass of M tot = 5 . ± . M (cid:12) .We find the disk is large, massive, and inclined at i disk = 137 . ± . ◦ . • Used 35 years of high resolution optical spectra toderive new radial velocities for both components inthe inner binary. We then fit a hierarchical tripleorbit and found a 241 day inner period (A-B) andan 11.5 year outer period (AB–C). When we com-bined the radial velocity constraints with the diskbased constraint on M tot , we found stellar massesof M A =2 . M (cid:12) , M B =1 . M (cid:12) , and M C =0 . M (cid:12) ,to a precision of ± . M (cid:12) . • Combined the radial velocity data with the astro-metric data from Berger et al. (2011) to performa joint RV-astrometric fit and found large mutualinclinations between the stellar orbits and the disk(Φ in = 44 ± ◦ , Φ out = 54 ± ◦ ). The stellarorbits may be mildly misaligned with each other(Φ in / out = 13 ± ◦ ). • Placed GW Ori A and B on the HR diagram, andfound that their stellar properties are broadly con-sistent with the predictions of pre-main sequencemodels for their measured masses at an age of ≈ • Compiled a lightcurve with a 30 year baseline andidentified several new eclipse events. We also iden-tified a 0.2 mag amplitude mode of variabilityphased with the outer orbital period, which sug-gests the A-B binary may be partially obscured bymicron-sized grains in the circumtriple disk cavityat outer apoastron. • Placed GW Ori in the context of other pre-mainsequence multiple systems. While short period ec-centric binary systems generally seem to have lowmutual inclinations with their circumbinary disks,there are a number of longer period systems thatexhibit significant mutual inclinations.Given its uniquely large stellar mass, massive cir-cumtriple disk, and puzzling eclipse behavior, GW Orishould remain a high priority target to study a uniqueclass of dynamical interactions in pre-main sequencemultiple systems.IC acknowledges support from the Smithsonian Insti-tution and Stanford KIPAC for the production of thismanuscript. IC and EJ acknowledge useful conversa-tions with Lisa Prato about GW Ori, and we greatlyappreciate her sharing her results in advance of pub-lication. IC acknowledges helpful conversations withMaxwell Moe and Jack Lissauer about GW Ori and re-lated systems; and would like to thank Eric Nielsen,4Bruce Macintosh, and Rebekah Dawson for helpful dis-cussions about astrometry and orbital dynamics. SAacknowledges the very helpful support provided by theNRAO Student Observing Support program related tothe early development of this project. GT acknowl-edges partial support for this work from NSF grantAST-1509375. This paper makes use of the followingALMA data: 2012.1.00496.S. ALMA is a partnershipof ESO (representing its member states), NSF (USA),and NINS (Japan), together with NRC (Canada) andNSC and ASIAA (Taiwan), in cooperation with the Re- public of Chile. The Joint ALMA Observatory is op-erated by ESO, AUI/NRAO, and NAOJ. This researchmade extensive use of the Julia programming language(Bezanson et al. 2017) and Astropy (Astropy Collabo-ration et al. 2013).
Software:
CASA (v4.4; McMullin et al. 2007), IRAF(Tody 1986, 1993), DiskJockey (Czekala et al. 2015a),RADMC-3D (Dullemond 2012), emcee (Foreman-Mackey et al. 2013), VARTOOLS (Hartman & Bakos2016), Astropy (Astropy Collaboration et al. 2013)REFERENCES
Akaike, H. 1973, Proceeding of the Second InternationalSymposium on Information Theory, B.N. Petrov and F.Caski, 267Andrews, S. M., Czekala, I., Wilner, D. J., et al. 2010, ApJ,710, 462Andrews, S. M., Rosenfeld, K. A., Kraus, A. L., & Wilner,D. J. 2013, ApJ, 771, 129Artymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,et al. 2013, A&A, 558, A33Avenhaus, H., Quanz, S. P., Schmid, H. M., et al. 2014,ApJ, 781, 87Bast, J. E., Brown, J. M., Herczeg, G. J., van Dishoeck,E. F., & Pontoppidan, K. M. 2011, A&A, 527, A119Bate, M. R. 2012, MNRAS, 419, 3115Berger, J.-P., Monnier, J. D., Millan-Gabet, R., et al. 2011,A&A, 529, L1Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B.2017, SIAM Review, 59, 65Biller, B., Lacour, S., Juh´asz, A., et al. 2012, ApJL, 753,L38Bouvier, J. 1990, AJ, 99, 946Brown, T. M., Baliber, N., Bianco, F. B., et al. 2013,PASP, 125, 1031Buchhave, L. A., Bakos, G. ´A., Hartman, J. D., et al. 2010,ApJ, 720, 1118Calvet, N., Muzerolle, J., Brice˜no, C., et al. 2004, AJ, 128,1294Capelo, H. L., Herbst, W., Leggett, S. K., Hamilton, C. M.,& Johnson, J. A. 2012, ApJL, 757, L18Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ,345, 245Cazzoletti, P., Ricci, L., Birnstiel, T., & Lodato, G. 2017,A&A, 599, A102Chiang, E. I., & Murray-Clay, R. A. 2004, ApJ, 607, 913Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 Czekala, I., Andrews, S. M., Jensen, E. L. N., et al. 2015a,ApJ, 806, 154Czekala, I., Andrews, S. M., Mandel, K. S., Hogg, D. W., &Green, G. M. 2015b, ApJ, 812, 128Czekala, I., Andrews, S. M., Torres, G., et al. 2016, ApJ,818, 156Czekala, I., Mandel, K. S., Andrews, S. M., et al. 2017,ApJ, 840, 49Dolan, C. J. 2000, PhD thesis, Department of Astronomy,University of Wisconsin–Madison, 475 North CharterStreet, Madison, WI 53706Dolan, C. J., & Mathieu, R. D. 2001, AJ, 121, 2124—. 2002, AJ, 123, 387Duchˆene, G., & Kraus, A. 2013, ARA&A, 51, 269Dullemond, C. P. 2012, RADMC-3D: A multi-purposeradiative transfer tool, Astrophysics Source CodeLibrary, ascl:1202.015
Dutrey, A., Di Folco, E., Beck, T., & Guilloteau, S. 2016,A&A Rv, 24, 5Eggleton, P. P., & Tokovinin, A. A. 2008, MNRAS, 389, 869Eisner, J. A. 2012, ApJ, 755, 23Fang, M., Sicilia-Aguilar, A., Roccatagliata, V., et al. 2014,A&A, 570, A118Fang, M., Sicilia-Aguilar, A., Wilner, D., et al. 2017, A&A,603, A132Fekel, Jr., F. C. 1981, ApJ, 246, 879Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman,J. 2013, PASP, 125, 306F˝ur´esz, G. 2008, Ph.D. Thesis, Univ. Szeged, HungaryGaia Collaboration, Brown, A. G. A., Vallenari, A., et al.2016, A&A, 595, A2Gelman, A., Carlin, J. B., Stern, H. S., et al. 2014, Bayesiandata analysis, Vol. 2 (CRC press Boca Raton, FL)Goodman, J., & Weare, J. 2010, Communications inApplied Mathematics and Computational Science, 5, 65Guilloteau, S., Dutrey, A., & Simon, M. 1999, A&A, 348,570 Hartman, J. D., & Bakos, G. ´A. 2016, Astronomy andComputing, 17, 1Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P.1998, ApJ, 495, 385Hauschildt, P. H., Allard, F., Ferguson, J., Baron, E., &Alexander, D. R. 1999, ApJ, 525, 871Henkel, C., Wilson, T. L., Langer, N., Chin, Y.-N., &Mauersberger, R. 1994, in Lecture Notes in Physics,Berlin Springer Verlag, Vol. 439, The Structure andContent of Molecular Clouds, ed. T. L. Wilson & K. J.Johnston, 72Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013,A&A, 553, A6Jensen, E. L. N., Dhital, S., Stassun, K. G., et al. 2007, AJ,134, 241Jonkheid, B., Faas, F. G. A., van Zadelhoff, G.-J., & vanDishoeck, E. F. 2004, A&A, 428, 511Kamp, I., & Dullemond, C. P. 2004, ApJ, 615, 991Kellogg, K., Prato, L., Torres, G., et al. 2017, ApJ, 844, 168Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017,PASP, 129, 104502Kounkel, M., Hartmann, L., Loinard, L., et al. 2017, ApJ,834, 142Kuhn, R. B., Rodriguez, J. E., Collins, K. A., et al. 2016,MNRAS, 459, 4281Kurtz, M. J., & Mink, D. J. 1998, PASP, 110, 934Lacour, S., Biller, B., Cheetham, A., et al. 2016, A&A, 590,A90Latham, D. W. 1992, in Astronomical Society of the PacificConference Series, Vol. 32, IAU Colloq. 135:Complementary Approaches to Double and Multiple StarResearch, ed. H. A. McAlister & W. I. Hartkopf, 110Li, G., Holman, M. J., & Tao, M. 2016, ApJ, 831, 96Lim, J., Yeung, P. K. H., Hanawa, T., et al. 2016, ApJ, 826,153Lomb, N. R. 1976, Ap&SS, 39, 447Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603Marino, S., Perez, S., & Casassus, S. 2015, ApJL, 798, L44Martin, R. G., & Lubow, S. H. 2017, ApJL, 835, L28Mathieu, R. D., Adams, F. C., Fuller, G. A., et al. 1995,AJ, 109, 2655Mathieu, R. D., Adams, F. C., & Latham, D. W. 1991, AJ,101, 2184McMullin, 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, 127Murray, C. D., & Correia, A. C. M. 2010, Keplerian Orbitsand Dynamics of Exoplanets, ed. S. Seager, 15 Najita, J., Carr, J. S., & Mathieu, R. D. 2003, ApJ, 589,931Nelson, A. F., & Marzari, F. 2016, ApJ, 827, 93Offner, S. S. R., Kratter, K. M., Matzner, C. D., Krumholz,M. R., & Klein, R. I. 2010, ApJ, 725, 1485Owen, J. E., & Lai, D. 2017, MNRAS, 469, 2834Prantzos, N., Aubert, O., & Audouze, J. 1996, A&A, 309,760Rosenfeld, K. A., Andrews, S. M., Hughes, A. M., Wilner,D. J., & Qi, C. 2013, ApJ, 774, 16Rosenfeld, K. A., Andrews, S. M., Wilner, D. J., &Stempels, H. C. 2012, ApJ, 759, 119Scargle, J. D. 1982, ApJ, 263, 835Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ,788, 48Shevchenko, V. S., Grankin, K. N., Ibragimov, M. A., &Melnikov, S. Y. 1992, Information Bulletin on VariableStars, 3746Shevchenko, V. S., Grankin, K. N., Ibragimov, M. A.,Mel’Nikov, S. Y., & Yakubov, S. D. 1993, Ap&SS, 202,121Shevchenko, V. S., Grankin, K. N., Mel’Nikov, S. Y., &Lamzin, S. A. 1998, Astronomy Letters, 24, 528Simon, M., Dutrey, A., & Guilloteau, S. 2000, ApJ, 545,1034Simon, M., Guilloteau, S., Di Folco, E., et al. 2017, ApJ,844, 158Siverd, R. J., Beatty, T. G., Pepper, J., et al. 2012, ApJ,761, 123Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ,131, 1163Stassun, K. G., Feiden, G. A., & Torres, G. 2014, NewAR,60, 1Takakuwa, S., Saigo, K., Matsumoto, T., et al. 2017, ApJ,837, 86Tody, D. 1986, in Proc. SPIE, Vol. 627, Instrumentation inastronomy VI, ed. D. L. Crawford, 733Tody, D. 1993, in Astronomical Society of the PacificConference Series, Vol. 52, Astronomical Data AnalysisSoftware and Systems II, ed. R. J. Hanisch, R. J. V.Brissenden, & J. Barnes, 173Tokovinin, A. 2017, ApJ, 844, 103Tokovinin, A. A. 1997, A&AS, 124Torres, G., Neuh¨auser, R., & Guenther, E. W. 2002, AJ,123, 1701Villa, A. M., Trinidad, M. A., de la Fuente, E., &Rodr´ıguez-Esnard, T. 2017, RMxAA, 53, 525Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409Yu, M., Evans, II, N. J., Dodson-Robinson, S. E., Willacy,K., & Turner, N. J. 2017, ApJ, 841, 396