GW Ori: Interactions Between a Triple-star System and its Circumtriple Disk in Action

GW Ori is a hierarchical triple system which has a rare circumtriple disk. We present Atacama Large Millimeter/submillimeter Array (ALMA) observations of 1.3 mm dust continuum and 12CO J=2-1 molecular gas emission of the disk. For the first time, we identify three dust rings in the disk at ~46, 188, and 338 AU, with estimated dust mass of ~70-250 Earth masses, respectively. To our knowledge, the outer ring in GW Ori is the largest dust ring ever found in protoplanetary disks. We use visibility modelling of dust continuum to show that the disk has misaligned parts and the innermost dust ring is eccentric. The disk misalignment is also suggested by the CO kinematics modelling. We interpret these substructures as evidence of ongoing dynamical interactions between the triple stars and the circumtriple disk.


INTRODUCTION
GW Ori is a hierarchical triple system (Berger et al. 2011) at a distance of 402 ± 10 parsec (Gaia Collaboration et al. 2018). Two of the stars (GW Ori AB) jiaqing.bi@gmail.com rbdong@uvic.ca compose a spectroscopic binary with a separation of ∼ 1 AU (Mathieu et al. 1991). A tertiary component (GW Ori C) was detected by near-infrared interferometry at a projected distance of ∼ 8 AU (Berger et al. 2011). The stellar masses have been constrained to be ∼ 2.7, 1.7, and 0.9 M , respectively (Czekala et al. 2017). The system harbours a rare circumtriple disk, with dust extending to ∼ 400 AU, and gas extending to ∼ 1300 AU (Fang et al. 2017). Spectral energy distribution (SED) modelling indicates a gap in the disk at 25−55 AU (Fang et al. 2014).
Here we present high resolution ALMA observations in the disk around GW Ori at 1.3 mm dust continuum emission and 12 CO J = 2 − 1 emission, where we find new substructures of the disk that indicate ongoing diskstar interactions. We arrange the paper as follows: In Section 2, we describe the setups of the ALMA observations and data reduction. In Section 3, we present the imaged results of dust continuum and 12 CO J = 2 − 1 observations. In Section 4, we present results of dust continuum visibility modelling. In Section 5, we discuss the possible origins of the observed substructures. In Section 6, we summarize our findings and raise some open questions.

OBSERVATION AND DATA REDUCTION
The observations were taken on December 10, 2017 (ID: 2017.1.00286.S). The disk was observed in Band 6 (1.3 mm) by 46 antennas, with baseline lengths ranging from 15 to 3321 meters. The total on source integration time was 1.6 hours. There were two 1.875-GHz-wide basebands centered at 217 and 233 GHz for continuum emission, and three basebands with 117 MHz bandwidths and 112 kHz resolution, centered at 230.518, 219.541 and 220.380 GHz to cover the 12 CO, 13 CO and C 18 O J = 2 − 1 lines.
The data were calibrated by the pipeline calibration script provided by ALMA. We used the Common Astronomy Software Applications package (casa; version 5.1. 1-5;McMullin et al. 2007) to process the data. We adopted casa task clean to image the continuum map ( Figure 1a), with the uniform weighting scheme and a 0".098 circular restoring beam. We performed phase self-calibration onto the image with a solution interval of 20 seconds. This resulted in an RMS noise level of ∼ 40 µJy beam -1 and an enhanced peak signal-to-noise ratio (SNR) of ∼ 157, compared with ∼ 63 µJy beam -1 and ∼ 90 before self-calibration, respectively. The integrated flux density of the disk (195 ± 20 mJy) is consistent with the result from previous ALMA observations (202 ± 20 mJy; Czekala et al. 2017).
The 12 CO J = 2 − 1 line data were obtained (after subtracting the continuum on the self-calibrated data) in the briggs weighting scheme with robust = 0.5, and velocity resolution 0.5 km s -1 . The resulting line cube has a beam size of 0".122 × 0".159 at the position angle -32.3 • . The noise level is ∼ 1.4 mJy beam -1 per channel and the peak signal-to-noise ratio is ∼ 83. Line emission was detected between 7.0 and 20.0 km s -1 with a central velocity of 13.5 km s -1 . The integrated flux is 60.6 Jy km s -1 , assuming a 2" radius. The intensity-weighted velocity map (a.k.a., the first-moment map) was constructed by calculating the intensity-weighted velocity with a threshold of three times the noise level. The averaged uncertainty of the twisted pattern in the firstmoment map (i.e., the inset in Figure 1b) is ∼ 0.2 km s -1 , derived from error propagation theory, assuming the uncertainty of velocity due to the channel resolution is 0.25 km s -1 . The observations of the CO isotopologues C 18 O and 13 CO J = 2 − 1 emission will be presented in future work. Figure 1a shows the continuum map with spatial resolution of 0".098 (∼ 39 AU). We identify three dust rings with different apparent shapes in the disk at ∼ 46, 188, and 338 AU (hereafter the inner, middle, and outer ring). The location of the inner ring coincides with the predicted cavity size from SED modelling (Fang et al. 2017). The continuum flux densities of the inner, middle, and outer ring are 42 ± 4, 95 ± 10 and 58 ± 6 mJy, respectively. To our knowledge, the outer ring is the largest ever found in protoplanetary disks.

Dust Continuum Emission
The three rings harbor an enormous amount of solid material. We estimate the dust (solid) mass M dust of the rings with the equation provided in Hildebrand (1983) where F ν is the continuum surface brightness at a submm frequency ν, d is the distance from the observer to the source, B ν (T dust ) is the Planck function at the dust temperature T dust , and κ ν is the dust opacity. The dust temperature is estimated using a fitting function provided by Dong et al. (2018a) T dust = 30 where L is the total stellar luminosity, and r is the location of the ring. The stellar luminosity modified by the distance provided by GAIA DR2 is 49.3 ± 7.4 L (Calvet et al. 2004;Gaia Collaboration et al. 2018). We assume a dust grain opacity of 10 cm 2 g -1 at 1000 GHz with a power-law index of 1 (Beckwith et al. 1990). We estimate the dust masses of the rings to be 74 ± 8, 168 ± 19, and 245 ± 28 M ⊕ , respectively, with the uncertainties incorporating the uncertainties in the surface brightness of the rings, source distance, stellar luminosity, and radial location of the rings.
3.2. 12 CO J=2-1 Emission Figure 1b shows the first-moment map of 12 CO J = 2−1 emission (with the zeroth-moment map provided in Appendix A). For regular Keplerian rotating disks, we expect a well-defined butterfly-like pattern in the firstmoment map. However, we find a twisted pattern inside ∼ 0".2, which may result from a misalignment between the inner and outer parts of the disk (i.e., having different inclinations and orientations; Rosenfeld et al. 2014), as has been found in the disks around, e.g., HD 142527 Marino et al. 2015) and HD 143006 Benisty et al. 2018).

MODELLING OF DUST AND GAS EMISSION
The different apparent shapes of the rings could result from a few scenarios, such as coplanar rings with different eccentricities, circular rings with different inclinations, or rings with both different eccentricities and inclinations. Here we present evidence for disk misalignment and disk eccentricity found in modelling the dust and gas emission.

Visibility Modelling of the Dust Continuum Emission
We fit the dust continuum map assuming that there are three dust rings in the disk with Gaussian radial profiles of surface brightness where F i is the surface brightness as a function of the distance to the center r, with i = 1, 2, 3 denoting parameters for the inner, middle, and outer ring, respectively. F 0,i is the peak surface brightness, r i is the radius of the ring (i.e., where the ring has the highest surface brightness), and σ i is the standard deviation. Initially, we assume all three rings are intrinsically circular when viewed face-on, and their different apparent shapes entirely originate from different inclinations. For each ring, we assume an independent set of peak surface brightness, center location, radius, width, inclination, and position angle as the model parameters. We call this combination of assumptions Model 1.
After projecting the rings according to their position angles and inclinations, we calculate the synthetic visibility of the models using galario (Tazzari et al. 2018), and launch MCMC parameter surveys to derive posterior distribution of model parameters using emcee (Foreman-Mackey et al. 2013). In the MCMC parameter surveys, the likelihood function L is defined as where V obs is the visibility data from ALMA observations, V mod is the synthetic model visibility, N is the total number of visibility data points in V obs , and m j is the weight of each visibility data point in V obs . The prior function is set to guarantee the surface brightness, ring radius, and ring width do not go below zero, the position angle does not go beyond (-90, 90) degrees, and the inclination does not go beyond (0, 90) degrees. For each model, there are 144 chains spread in the hyperspace of parameters. Each chain has 15000 iterations including 10000 burn-in iterations. The results of the parameter surveys are listed in Table 1. The fitting result of Model 1, listed in Table 1a, suggests that the three rings have statistically different centers. Particularly, the center of the inner ring differs from the centers of the outer two by ∼ 20% of the inner ringâĂŹs radius. This non-concentricity indicates nonzero intrinsic eccentricities in the rings, particularly the inner ring (see Section 5.1).
We explore the non-zero intrinsic eccentricity in the inner ring with two models. In both models, the outer two rings are intrinsically circular and concentric. Their center coincides with one of the two foci of the inner ring. In Model 2, that center is set free, while in Model 3 it is assumed to coincide with the stellar position provided by GAIA DR2 (Gaia Collaboration et al. 2018). In those two models, we introduce two more parameters for the intrinsic eccentricity and apoapsis angle of the inner ring. The position angle only indicates the direction to the ascending node on the axis along which the ring is inclined. The fitting results are listed in Table 1b and 1c, and the following calculations are based on the result of Model 3. Figure 1c and 1e show the model image and the residual map of Model 3, respectively. The residual map is produced by subtracting model from data in the visibility plane, and then imaging the results in the same way used for the observations. We interpret the residuals as additional substructures on top of the ideal model (e.g., a warp within the ring; Huang et al. 2020).
All three models yield roughly consistent inclinations and position angles of each ring. However, we cannot determine the mutual inclinations between them (i.e., the angles between their angular momentum vectors) from dust emission modelling alone, due to the unknown direction of orbital motion.

Kinematics Modelling of the
Following the prescription and parameter values used to fit low resolution CO isotopologue data of GW Ori (Fang et al. 2017), we set up a gas surface density model using a power-law profile with an exponential tail and the aspect ratio h/r parametrized as where Σ c and (h/r) c are corresponding values at the characteristic scaling radius r c . The disk mass is taken as 0.12 M , corresponding to Σ c = 3 g cm −2 for r c = 320 AU, with γ = 1.0, (h/r) c = 0.18 and ψ = 0.1. The dust surface density profile is set by assuming a gas-to-dust ratio of 100, and decreasing the dust surface density by a factor of 1000 inside the derived gap radii: inside 37 AU, from 56 to 153 AU, and from 221 to 269 AU. The 12 CO channel maps are then computed and ray-traced by the physical-chemical modeling code dali (Bruderer 2013), which simultaneously solves the heating-cooling balance of the gas and chemistry to determine the gas temperature, molecular abundances and molecular excitation for a given density structure. Similar to Walsh et al. (2017), we model the misaligned disk with an inner cavity and three annuli each with its own inclination and position angle 1 , as listed in Table 2. The channel map is run through the ALMA simulator using the settings of the ALMA observations. The resulting first-moment map is shown in Figure 1d. In Figure 1f we show the simulated first-moment map for another model as a comparison, in which the disk is coplanar with an inclination of 37.9 • and a position angle of -5 • .
The models show that the 12 CO J = 2 − 1 firstmoment map in the ALMA observation cannot be reproduced by a coplanar disk. Instead, the misaligned disk model described in Table 2 matches the observed first-moment map better, indicating the presence of misalignment in the GW Ori disk.

DISCUSSIONS
Several disks have been observed to have non-zero eccentricity and/or misalignment (e.g., MWC 758, Dong et al. 2018b;HD 142527, Casassus et al. 2015;Marino et al. 2015;and HD 143006, Pérez et al. 2018;Benisty et al. 2018). Unlike most of them, in which the origin is uncertain, the GW Ori system provides strong and 1 dali is unable to vary inclination as a function of radius. The final channel map is constructed by concatenating three channel maps, each for one component, ray-traced at its inclination and position angle and cut out at the specified radius range listed in Table 2.
direct link between substructures and star-disk gravitational interactions. Therefore it offers a unique laboratory to probe three-dimensional star-disk interactions.
In this section, we discuss the possible origins of the observed substructures due to star-disk interactions.

Disk Eccentricity
The A-B binary and the C component can be dynamically viewed as an AB-C binary. The eccentricity of the circumbinary disk may increase through resonant interactions with the binary (Papaloizou et al. 2001). In the case of no binary-disk misalignment, the binaryâĂŹs perturbing gravitational potential on the midplane of the disk is given by Lubow (1991). The coupling of this perturbing potential with the imposed eccentricity of the disk excites density waves at the 1:3 outer eccentric Lindblad Resonance, which lead to angular momentum removal in the inner parts of the disk. As no energy is removed along with the angular momentum in this process, the disk orbit cannot remain circular (Papaloizou et al. 2001). In the case of GW Ori, the inner dust ring is the most susceptible to this effect, which could explain why its center in Model 1 is more deviated from those of the other two rings.

Binary-disk Misalignment
Our dust and gas observations alone cannot break the degeneracy in the mutual inclination between different parts in the disk due to the unknown on-sky projected orbital direction of the disk. Previous studies indicate that the on-sky projected gas motion is likely to be clockwise (Czekala et al. 2017), same as the orbital motion of GW Ori C given by astrometric observations (Berger et al. 2011). Given the inclination and longitude of ascending node of the AB-C binary orbit being 150 ± 7 and 282 ± 9 degrees (Czekala et al. 2017), we assume that the entire disk has the same clockwise on-sky projected orbital direction. Following Fekel (1981), we find out that the binary-disk misalignments at 46 AU (the inner ring), 100 AU (a gap), 188 AU (the middle ring), and 338 AU (outer ring) are 11 ± 6, ∼28 2 , 35 ± 5, and 40 ± 5 degrees, respectively. A schematic diagram of our disk model is displayed in Figure 2. Therefore, the inner ring and the AB-C binary plane are close to being coplanar, and there is a monotonic trend of binary-disk misalignment from ∼ 10 • at ∼ 50 AU to ∼ 40 • at ∼ 340 AU, consistent with the expected outcome of the disk misalignment (see Section 5.3).
Several mechanisms could produce an initial binarydisk misalignment, such as turbulence in the starforming gas clouds (Bate 2012), binary formation in the gas cloud whose physical axes are misaligned to the rotation axis (Bonnell & Bastien 1992), and accretion of cloud materials with misaligned angular momentum with respect to the binary after the binary formation (Bate 2018).

Misalignment Within the Disk
A test particle orbiting a binary on a misaligned orbit undergoes nodal precession due to gravitational perturbations from the binary (Nixon et al. 2011;Facchini et al. 2018). For a protoplanetary disk in the bending-wave regime (i.e., where the aspect ratio is higher than the α-prescription of viscosity; Shakura & Sunyaev 1973), disk parts at different radii shall undergo global precession like a rigid body with possibly a small warp (Smallwood et al. 2019). Therefore the timescale of radial communication of disk materials (i.e., for pressure-induced bending waves propagating at half of the sound speed) and the timescale of global precession determine whether the disk can develop a misalignment inside.
Assuming an inner radius at 32 AU (3 − 4 times the AB-C binary semi-major axis; Czekala et al. 2017;Kraus et al. 2020) and an outer radius at 1300 AU (size of the gas disk, Fang et al. 2017), the global precession timescale of the entire disk is ∼ 0.83 Myr, and the radial communication time-scale is ∼ 0.06 Myr (see Appendix D.1 and D.2 for detailed calculations). The radial communication is able to prevent the disk from breaking or developing a significant warp, and we would not expect the observed large deviations in the inclination and position angle between the inner and middle ring. Therefore, we propose that the gap between the inner and middle ring is deep enough to break the disk into two parts (hereafter the inner and outer disk), undergoing nodal precession independently, due to another mechanism.
Due to the viscous dissipation, the precessing disk is torqued toward either polar alignment (i.e., the binarydisk misalignment becomes 90 • ; Martin & Lubow 2017; Zanazzi & Lai 2018), or coplanar alignment/counteralignment. The minimum critical initial binary-disk misalignment for which a disc moves towards polar alignment is ∼ 63 • in the limit of zero disk mass. Since a higher disk mass will lead to a larger critical angle (Martin & Lubow 2019), the GW Ori disk is most likely moving towards coplanar alignment. As we propose that the disk breaks into two parts undergoing global precession independently, they are also aligning to the binary independently on different time-scales.
Assuming the radial communication is blocked at 60 AU, we estimate the alignment time-scales to be ∼ 1 Myr for the inner disk and above 100 Myr for the outer disk (see Appendix D.3). This is consistent with the observed significantly smaller inclination of the inner ring with respect to the binary than those of the outer rings. The latter are likely inherited from birth and have not evolved much given the system age.
If the radial communication is also blocked at the gap between the middle and outer ring (e.g., ∼ 250 AU), the nodal precession time-scales for the middle and outer ring would be ∼ 0.6 Myr and ∼ 120 Myr, respectively, and the two rings are likely to develop significantly different position angles. Thus we propose that the gap between the middle and outer ring does not cut off the radial communication, and the two rings precess roughly as a rigid body with only a small warp between them.

Hydrodynamic Simulations
The analytical results suggest that the radial communication is able to prevent the binary from breaking the disk (e.g. Nixon et al. 2013). As a result, we propose a break at ∼ 60 AU that is due to other mechanisms in order to explain the observed structures. We carry out a demonstrative smoothed particle hydrodynamic (SPH) simulation with the phantom code (Lodato & Price 2010;Price & Federrath 2010;Price et al. 2018) to test the non-breaking hypothesis in the non-linear regime. The results are shown in Figure 3.
We model the triple star system as the outer binary in order to speed up the simulation. The simulation consists of 10 6 equal mass Lagrangian SPH particles initially distributed from r in = 40 AU to r out = 400 AU. The initial truncation radius of the disk does not affect the simulation significantly, since the material moves inwards quickly due to the short local viscous timescale. A smaller initial outer truncation radius r out than what is observed is chosen in order to better resolve the disk. The binary begins at apastron with e b = 0.22 and a b = 9.2 AU (Czekala et al. 2017). The accretion radius of each binary component is 4 AU. Particles within this radius are accreted, and their mass and angular momentum are added to the star. We ignore the effect of self-gravity since it has no effect on the nodal precession rate of flat circumbinary disks.
The initial surface density profile is taken by where Σ 0 is the density normalization at r 0 = 40 AU, corresponding to a total disk mass of 0.1 M . We take a locally isothermal disk with a constant aspect ratio h/r = 0.05, where h is the scale-height. The Shakura & Sunyaev (1973) α parameter varies in the range 0.008 − 0.013 over the disk. The SPH artificial viscosity α AV = 0.31 mimics a disk with α ≈ α av 10 H h (8) Lodato & Price (2010), where H is the mean smoothing length on particles in a cylindrical ring at a given radius and we take β AV = 2. The disk is resolved with average smoothing length per scale height of 0.32. The evolution of surface density, binary-disk misalignment, and longitude of the ascending node at different radii suggest that the disk does not show any sign of breaking in 3000 binaryâĂŹs orbital periods (∼ 0.04 Myr), which is sufficiently long to tell if the disk would break or not since the radial communication time-scale in the simulation is ∼ 0.01 Myr. Instead, the disk presents a global warp. The warp is not taken into account in the analytic estimates. The small outer truncation radius in the simulation leads to a faster precession time-scale than that predicted by the analytic model, comparable to the radial communication timescale. The simulations suggests that unless the disk is very cool (i.e., low aspect ratio) and in the viscous regime, some other mechanism, e.g., a companion, is needed to break the disk at the gap between the inner and middle ring. This mechanism may also be responsible for producing the observed misalignment in the disk.
A disk with a lower aspect ratio or a higher α value, such that it falls into the viscous regime (h/r < α), may break due to the binaryâĂŹs torque (Nixon et al. 2013). However, observations have suggested lower α values than our adaption here (e.g., α 10 −3 , Flaherty et al. 2017). A lower viscosity leads to a larger binary truncation radius (Artymowicz & Lubow 1994), and a longer global precession time-scale (Equation D4). Therefore we expect even less warping (and no break) in the GW Ori disk than in our simulation.

CONCLUSIONS
We present ALMA 1.3 mm dust continuum observation and 12 CO J = 2 − 1 emission of the circumtriple disk around GW Ori. Our main conclusions are the following: 1. For the first time, we identify three dust rings in the GW Ori disk at ∼ 46, 188, and 338 AU, with their estimated dust mass being ∼ 74, 168, and 245 M ⊕ , respectively. The three dust rings have enough solids to make many cores of giant planets (∼ 10 M ⊕ ; Pollack et al. 1996).
2. We built three models under various assumptions to fit the dust continuum observations using MCMC fitting. Our results (Table 1) suggest that the inner ring has an eccentricity of ∼ 0.2, and the three rings have statistically different on-sky projected inclinations. The inner, middle, and outer ring are likely misaligned by ∼ 11, 35, and 40 degrees to the orbital plane of the GW Ori AB-C binary system, respectively.
3. A twisted pattern is identified in the first-moment map, suggesting the presence of a warp in the disk, consistent with what we have found in the dust continuum emission.
4. Using analytical analysis and hydrodynamic simulations, we find that the torque from the GW Ori triple stars alone cannot explain the observed large misalignment between the inner and middle dust rings. The disk would not break due to the torque, and a continuous disk is unlikely to show the observed large misalignment. Therefore this hints at some other mechanism that breaks the disk and prevents radial communication of bending waves between the inner and middle ring.
There are still open questions associated with the system. For example, are there any companions in the disk? Dust rings and gaps have been shown to be common in protoplanetary disks (Long et al. 2018;Andrews et al. 2018;Huang et al. 2018;van der Marel et al. 2019), and one of the most exciting hypotheses is that they are produced by embedded companions ranging from stellarmass all the way to super-Earths (Artymowicz & Lubow 1994;Dong et al. 2015;Zhang et al. 2018). Specifically, a companion may be opening the gap between the inner and middle ring and break the disk there. Companions at hundreds of AU from their host stars have been found before (e.g., HD 106906 b; Bailey et al. 2014). But how they form, i.e., forming in situ or at closer distances, or followed by scattering or migration to the outer regions, is unclear. If GW OriâĂŹs dust rings are in the process of forming companions, there will be circumtriple companions, which have not been found before (excluding quadruple systems; Busetti et al. 2018). The system will offer direct clues on the formation of distant companions.
We thank Sean Andrews, Myriam Benisty, John Carpenter, Ian Czekala, Sheng-Yuan Liu, Feng Long, Diego MuÃśoz, Henry Ngo, Laura PÃľrez, John Zanazzi, and Zhaohuan Zhu for discussions. We also thank the anonymous referee for constructive suggestions that largely improved the quality of the paper. J.B. thanks Belaid Moa for helps on the numerical implementation. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00286.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Numerical calculations are performed on the clusters provided by Com-puteCanada. This work is in part supported by JSPS KAKENHI Grant Numbers 19K03932, 18H05441 and      Note-The annuli are concentric with the center located at the stellar position provided by GAIA DR2 (ICRS R.A. = 5 h 29 m 08 s .390 and Dec. = 11 • 52'12".661). The longitude of the ascending node is measured East of North, and the inclination is defined in the range from 0 • to 90 • with 0 • meaning face-on. The disk model is composed of an empty inner cavity and three annuli, inside out. The inner annulus (from 32 to 48 AU) is for the inner dust ring. The outer annulus (from 153 to 1000 AU) is for the middle and outer dust rings. The middle annulus (from 48 to 153 AU) is for the gap in between.

APPENDIX
A. ALMA 12 CO J = 2 − 1 ZEROTH-MOMENT MAP Figure 4 shows the ALMA zeroth-moment map of 12 CO J = 2 − 1 emission. The spiral-like structure at ∼ 0.75 to the northwest is likely due to cloud contamination, as reported by previous studies (Czekala et al. 2017;Fang et al. 2017).