GW Ori: circumtriple rings and planets

ABSTRACT

GW Ori is a hierarchical triple star system with a misaligned circumtriple protoplanetary disc. Recent Atacama Large Millimeter/submillimeter Array observations have identified three dust rings with a prominent gap at |$100\, \rm au$| and misalignments between each of the rings. A break in the gas disc may be driven by the torque from either the triple star system or a planet that is massive enough to carve a gap in the disc. Once the disc is broken, the rings nodally precess on different time-scales and become misaligned. We investigate the origins of the dust rings by means of N-body integrations and 3D hydrodynamic simulations. We find that for observationally motivated parameters of protoplanetary discs, the disc does not break due to the torque from the star system. We suggest that the presence of a massive planet (or planets) in the disc separates the inner and outer discs. We conclude that the disc breaking in GW Ori is likely caused by undetected planets – the first planet(s) in a circumtriple orbit.

accretion, accretion discs, hydrodynamics, planet–disc interactions, stars: individual: GW Ori

1 INTRODUCTION

Observations show that the majority of stars form in relatively dense regions within stellar clusters, naturally leading to multistar systems. It is estimated that more than |$40$|–|$50{{\ \rm per\ cent}}$| of stars are in a binary pair, while about |$20{{\ \rm per\ cent}}$| are observed in triple or higher order systems (Duquennoy & Mayor 1991; Raghavan et al. 2010; Tokovinin 2014a,b). To date, there have been planets found in 32 triple star systems (Busetti, Beust & Harley 2018; Fragione, Loeb & Ginsburg 2019). However, as yet no planet in a circumtriple orbit has been discovered.

GW Ori is one such hierarchical triple star system, with a spectroscopic binary (A and B) at a separation of about |$1\, \rm au$| and a tertiary stellar companion (C) at a projected distance of about |$8\, \rm au$| (Mathieu, Adams & Latham 1991; Berger et al. 2011; Czekala et al. 2017; Kraus et al. 2020). The system is at a distance of |$408\pm 10\, \rm pc$| (Gaia Collaboration 2021) and has an age of |$1.0\pm 0.1\, \rm Myr$| (Calvet et al. 2004). The derived triple star parameters from Czekala et al. (2017) are shown in Table 2. Stars A and B have masses of |$2.8\, $| and |$1.7\, \rm M_{\odot }$|⁠, respectively, while star C has a mass of |$1.2\, \rm M_{\odot }$| and the eccentricity of the outer companion is estimated to be 0.22 (Czekala et al. 2017). The system has a circumtriple protoplanetary disc in which the dust is observed to extend to about |$400\, \rm au$|⁠, while the gas extends farther to about |$1300\, \rm au$| (Fang et al. 2017).

Bi et al. (2020) presented Atacama Large Millimeter/submillimeter Array (ALMA) observations of the dust continuum and molecular gas emission of the circumtriple disc of GW Ori. They identified three dust rings in the circumtriple disc at radii of about 46, 188, and |$338\, \rm au$| from the triple star. The structure of the rings is outlined in Table 1. The radial width of dust ring 3 is the widest ring ever found in a protoplanetary disc. The dust masses of the rings 1, 2, and 3 were estimated to be 74 ± 8, 168 ± 19, and |$245\pm 28\, \rm M_{\oplus }$|⁠, respectively (Bi et al. 2020). Misalignments are present between each individual dust ring and triple star orbital plane from visibility modelling of dust continuum and CO kinematics. The inclination and the longitude of ascending node of the AB–C binary orbit were found to be 150 ± 7^○and 28 ± 9^○, respectively (Berger et al. 2011). Assuming that the entire disc has the same clockwise on-sky projected orbital direction, Bi et al. (2020) found that rings 1, 2, and 3 were misaligned by 11 ± 6^○, 35 ± 5^○, and 40 ± 5^○, respectively, relative to the orbital plane. Therefore, their results suggested that ring 1 and the AB–C binary plane are close to being coplanar, while rings 2 and 3 are misaligned but close to coplanar with each other.

Table 1.

The location of the dust rings in the GW Ori circumtriple disc. The values were obtained by Markov chain Monte Carlo (MCMC) fitting from Bi et al. (2020). The second column gives the inner edge of the ring. The third column denotes the centre location of the ring, and the fourth column represents the outer edge of the ring.

Ring	Inner radius	Centre radius	Outer radius
	(au)	(au)	(au)
1	|$36.9804^{+0.134}_{-0.152}$|	|$46.5516^{+ 0.05}_{-0.048}$|	|$56.1228^{+0.134}_{-0.152}$|
2	|$153.9865^{+0.2534}_{-0.8105}$|	|$188.2827^{+0.1065}_{-0.1037}$|	|$222.5790^{+0.2534}_{-0.2567}$|
3	|$270.5542^{+0.7514}_{-0.8105}$|	|$337.2438^{+0.3735}_{-0.3602}$|	|$403.9334^{+0.7514}_{-0.8105}$|

A potential explanation for the large misalignment between the inner and middle rings is the ‘disc breaking’ phenomenon, where the torque from the misaligned stars overcomes the viscous stresses and pressure holding the disc together and separates the disc into distinct planes (e.g. Nixon, King & Pringle 2011; Nixon et al. 2012; Facchini, Lodato & Price 2013; Nixon, King & Price 2013). However, Bi et al. (2020) ran smoothed particle hydrodynamics (SPH) simulations modelling the disc with an inner, misaligned binary (that approximates the triple star evolution) and found that the disc does not break only due to the triple star torque. They modelled the disc with an aspect ratio H/r = 0.05. For this case, the disc did undergo strong warping at the beginning of the simulation but later relaxed to a steady state resembling a flat disc. They conclude that the gaps in the dust cannot solely be caused by the gravitational torque of the system and thus there must be another cause for the gaps.

Recently, Kraus et al. (2020) presented higher resolution ALMA dust observations of GW Ori. They updated the triple star parameters that are shown in Table 2. They found that |$M_{\rm A} = 2.47 \pm 0.33\, \rm M_{\odot }$|⁠, |$M_{\rm B} = 1.43 \pm 0.18\, \rm M_{\odot }$|⁠, and |$M_{\rm C} = 1.36 \pm 0.28 \, \rm M_{\odot }$| and the eccentricity of the outermost companion was increased to e_AB–C = 0.379. They found three distinct coherent dust rings similar to Bi et al. (2020). The eccentricity of ring 1 is found to be ∼0.2. Moreover, Kraus et al. (2020) ran SPH simulations to investigate the origin of the substructures located in the circumtriple disc. They found that the torque from the triple system could effectively break the disc, finding a relative misalignment between rings 1 and 2 that was consistent with the observations.

Table 2.

The orbital elements of the GW Ori triple system from Czekala et al. (2017) and Kraus et al. (2020).

Parameter	Czekala et al. (2017)		Kraus et al. (2020)
	Orbit A–B	Orbit (AB)–C	Orbit A–B	Orbit (AB)–C
Period P (d)	241.50 ± 0.05	4246 ± 66	241.619 ± 0.05	4216.8 ± 4.6
Semimajor axis a (au)	1.25 ± 0.05	9.19 ± 0.32	1.20 ± 0.04	8.89 ± 0.04
Eccentricity e	0.13 ± 0.02	0.22 ± 0.09	0.069 ± 0.009	0.379 ± 0.003
Inclination i (○)	157 ± 1	150 ± 7	156 ± 1	149.6 ± 0.7
Longitude of the periastron ωa (○)	17 ± 7	130 ± 21	1 ± 7	105 ± 1
Longitude of the ascending node Ωa (○)	263 ± 13	282 ± 9	258.2 ± 1.3	230.9 ± 1.1
Total mass Mtot (⁠|$\rm M_{\odot }$|⁠)	|$4.48^{+0.42}_{-0.36}$|	|$5.63^{+0.58}_{-0.43}$|	3.90 ± 0.40	5.26 ± 0.22
Star A mass MA (⁠|$\rm M_{\odot }$|⁠)	|$2.80^{+0.36}_{-0.31}$|		2.47 ± 0.33
Star B mass MB (⁠|$\rm M_{\odot }$|⁠)	|$1.68^{+0.21}_{-0.18}$|		1.43 ± 0.18
Star C mass MC (⁠|$\rm M_{\odot }$|⁠)	|$1.15^{+0.40}_{-0.23}$|		1.36 ± 0.28

The major focus of this work is to investigate the dynamical origin of the dust rings that are present in the GW Ori inclined circumtriple disc. In Section 2, we begin by considering the behaviour of test particle orbits around the triple star versus a binary. In Section 3, we describe our numerical set-up for the hydrodynamical simulations, and the results are provided in Section 4. We discuss our findings in the context of the observations of GW Ori in Section 5 and summarize in Section 6.

2 FOUR-BODY DYNAMICS

In this section, we simulate the GW Ori system using the whfast integrator that is a second-order symplectic Wisdom Holman integrator with 11th-order symplectic correctors in the N-body simulation package, rebound (Rein & Tamayo 2015). We explore the effect of modelling the close spectroscopic binary (stars A and B) as a single star, making the GW Ori system a binary (stars AB and C). We consider both a triple system (stars A, B, and C) and a binary system (stars AB and C). We solve the gravitational equations for the three (two) bodies in the frame of centre of mass of the three- (two-) body system. The parameters used for both the triple system and the binary system are given in Table 2.

Test particle orbits around binary and higher order systems can show qualitatively similar behaviour to rings in a gas disc (e.g. Doolin & Blundell 2011; Martin et al. 2014; Chen et al. 2019; Smallwood et al. 2019). Thus, in order to investigate the effect of simplifying the three-star configuration to a binary system, we simulate the evolution of test particles around each stellar geometry in Fig. 1. The stars’ orbits are similar regardless of using three stars or a binary, and additionally, the orbits are stable over a long time-scale. We select four different semimajor axes of the particles of 47, 188, 337, and |$100\, \rm au$|⁠. The first three separations denote the centre radius of the three observed dust rings, and the last one represents the centre of the observed inner gap. The tilt of the test particle, i_tp, is measured from the z-axis and the phase angle of the test particle, ϕ_tp, is measured using the angular momentum unit vector of the test particle in the y and x directions. For the binary case, the test particles nodally precess about the angular momentum vector of the binary and show small tilt oscillations as a result of the eccentricity of the AB–C binary (e.g. Smallwood et al. 2019). Test particles around the triple star system also undergo small tilt oscillations that are driven by the precession of the triple stars.

Figure 1.

Test particle evolution around the GW Ori triple star system (left-hand panel) and around the approximated binary star system (right-hand panel). We show the semimajor axis a_tp, inclination i_tp, and longitude of the ascending node ϕ_tp of the test particles as a function of time. The stellar parameter set used is from Kraus et al. (2020). We show four different test particles that begin with an initial semimajor axis of 47 (blue), 188 (green), 337 (red), and |$100\, \rm au$| (yellow). The first three separations denote the centre radius of the three observed dust rings, and the last one represents the centre of the observed inner gap.

The time-scale for these tilt oscillations increases with test particle separation. The particle dynamics are not significantly affected by the presence of the AB binary (the left-hand panel includes the triple star system and the right-hand panel includes the binary approximation). The amplitude of the tilt oscillations is slightly greater around the triple star, by only a few degrees, than around the binary star approximation. The frequency of the tilt oscillations occurs on a faster time-scale around a binary than around the triple star. However, in both cases, the amplitude of the tilt oscillations is small.

The nodal precession time-scale for the test particles is very similar around the triple star compared to the binary star approximation. Since the nodal phase angle is changing on a much faster time-scale than the tilt, we expect that in a disc simulation the disc will remain relatively unwarped (constant inclination with radius) but it may become twisted (variation of nodal phase angle with radius). Disc breaking as a result of the binary is therefore likely driven by the precession. Since the precession is very similar in the binary and triple star cases, we expect a similar amount of warp/twist to the disc around two stars as around a triple star system. In Section 4, we first use hydrodynamical gas disc simulations around the binary and triple star systems to confirm this.

3 NUMERICAL METHODS

For our numerical simulations, we use the 3D SPH (e.g. Price 2012) code phantom (Price et al. 2018). We simulate an accretion disc around a binary and triple star system. phantom has been well tested and used to model misaligned accretion discs in binary systems (e.g. Nixon et al. 2013; Martin et al. 2014; Franchini, Martin & Lubow 2019a). A misaligned disc feels a gravitational torque exerted by a binary or triple star. This causes the disc to undergo differential precession that can lead to disc warping, ‘breaking’ (Facchini et al. 2013; Nixon et al. 2013) or ‘tearing’ (Nixon 2012). Dissipation within the misaligned disc causes the differential precession to generate a warp that will evolve in the diffusive regime or the bending wave regime, depending on the disc thickness and viscosity (Papaloizou & Pringle 1983; Papaloizou & Terquem 1995; Ogilvie 1999). In the bending wave regime, the disc aspect ratio is larger than the viscosity coefficient (H/r ≳ α), and the warp induced in the disc by the binary torque propagates as a pressure wave with speed c_s/2 (Papaloizou & Pringle 1983; Papaloizou & Terquem 1995). For the diffusive regime, the disc aspect ratio is less than the viscosity (H/r ≲ α) with a diffusion coefficient inversely proportional to the disc viscosity. Most of our simulations are in the bending wave regime. Although we do consider some simulations in the diffusive regime for demonstrative purposes, we note that protoplanetary discs are expected to be in the bending-wave regime (Hartmann et al. 1998). We discuss this further in Section 5.3.

In total, we conduct 11 simulations that are summarized in Table 3 and we present additional information of our simulation parameters in Table A1. First, we compare the disc structure around three stars versus two stars. Secondly, we compare the two different sets of stellar parameters from Czekala et al. (2017) and Kraus et al. (2020). We then compare various disc aspect ratios and Shakura & Sunyaev (1973) α values. Furthermore, we also run two simulations that include a planet that is initially coplanar to the disc. The simulations without a planet are simulated for |$3000\, \rm \mathrm{ \mathit{ P}}_{\rm orb}$| and simulations with a planet are simulated for |$230\, \rm \mathit{ P}_{\rm planet}$|⁠, where P_orb is the binary orbital period and P_planet is the planet orbital period. For a |$1\, \rm \mathit{ M}_{\rm J}$| planet at |$100\, \rm au$|⁠, 1P_planet ≈ 36P_orb, so therefore the planet simulations ran for |${\sim}8000\, \rm \mathit{ P}_{\rm orb}$|⁠. The simulations without a planet have reached a steady state within |$3000\, \rm \mathrm{ \mathit{ P}}_{orb}$|⁠. Kraus et al. (2020) ran their simulation for a much shorter time of |${\sim}860\, \rm \mathrm{ \mathit{ P}}_{\rm orb}$|⁠.

Table 3.

The set-up of the SPH simulations that lists the number of sink particles, inner circumbinary (circumtriple) disc radius r_in, outer circumbinary (circumtriple) disc radius r_out, initial circumbinary (circumtriple) disc tilt i₀, disc aspect ratio H/r, the input Shakura & Sunyaev (1973) α-viscosity parameter, surface density power-law index p, the mean shelled-averaged smoothing length per scale height 〈h〉/H, mass of the planet M_p, initial inclination of the planet |$i_{0,\rm p}$|⁠, and the specific set of binary parameters. The bold highlights the defining parameter for each simulation.

Simulation	# of Stars	rin	rout	i0	H/r	α	p	(〈h〉/H)mean	Mp	i0, p	Stellar parameters
		(au)	(au)	(○)					(MJ)	(○)
run0	|$\bf 3$|	40	200	38	0.05	0.01	0.5	0.30	–	–	Kraus et al. (2020)
run1	2	40	400	38	|$\bf 0.05$|	0.01	1.5	0.32	–	–	Czekala et al. (2017)
run2	2	40	400	40	0.05	0.01	1.5	0.32	–	–	Kraus et al. (2020)
run3	2	40	400	40	|$\bf 0.1$|	0.01	1.5	0.20	–	–	Czekala et al. (2017)
run4	2	40	400	40	|$\bf 0.05$|	0.01	1.5	0.32	1	40	Czekala et al. (2017)
run5	2	40	400	40	|$\bf 0.1$|	0.01	1.5	0.20	1	40	Czekala et al. (2017)
run6	2	|$\bf 40$|	200	38	0.05	0.01	0.5	0.30	–	–	Kraus et al. (2020)
run7	2	40	200	38	0.05	|$\bf 0.05$|	0.5	0.30	–	–	Kraus et al. (2020)
run8	2	40	200	38	0.05	|$\bf 0.1$|	0.5	0.30	–	–	Kraus et al. (2020)
run9	2	|$\bf 30$|	200	38	0.05	0.1	0.5	0.32	–	–	Kraus et al. (2020)
run10	2	|$\bf 20$|	200	38	0.05	0.01	0.5	0.34	–	–	Kraus et al. (2020)

Table A1.

Additional information for the SPH simulations that lists the disc aspect ratio at the disc inner and outer edges H/r, the minimum and maximum values of the Shakura & Sunyaev (1973) α-viscosity parameter, and the minimum and maximum values of the shelled-averaged smoothing length per scale height 〈h〉/H.

Simulation	αmin	αmax	(〈h〉/H)min	(〈h〉/H)max
run0	0.006	0.021	0.21	0.71
run1	0.008	0.013	0.27	0.42
run2	0.008	0.013	0.27	0.42
run3	0.008	0.013	0.17	0.27
run4	0.008	0.013	0.27	0.42
run5	0.008	0.013	0.17	0.27
run6	0.006	0.021	0.21	0.71
run7	0.007	0.017	0.21	0.49
run8	0.035	0.083	0.21	0.49
run9	0.071	0.165	0.21	0.49
run10	0.007	0.0184	0.21	0.58

We additionally generate synthetic CO maps for two of our simulations for comparison with the ALMA ¹²CO J = 2–1 first moment map from Bi et al. (2020). For this, we use the Monte Carlo radiative transfer code mcfost (Pinte et al. 2006, 2009). mcfost is particularly well suited for use with a particle-based numerical method because it uses a Voronoi mesh (rather than a cylindrical or spherical grid) to generate a grid from the particles. Because the mesh follows the particle distribution it does not require any interpolation.

3.1 Triple star parameters

We set-up one hydrodynamical simulation with three stars that are modelled as sink particles. The triple star parameters that we adopt are from Kraus et al. (2020), which is summarized in Table 2. The second column in Table 3 details the stellar set-up (either a binary or triple star system) used in each simulation. The inner binary and tertiary star begin at apastron. The components of the tight A-B binary have an accretion radius of |$R_{\rm acc}=0.5\, \rm au$|⁠, while the tertiary companion has |$R_{\rm acc}=2.3\, \rm au$|⁠. Particles within the hard accretion radius are considered accreted and their mass and angular momentum are added to the respected star.

3.2 Binary parameters

In Section 2, we showed that we can model the GW Ori hierarchical triple system as an AB–C binary. To model the binary we use the two sets of binary orbital parameters measured by Czekala et al. (2017) and Kraus et al. (2020). Table 2 compares the binary parameters from each study. The last column in Table 3 details which simulation uses which set of binary parameters. The binary begins at apastron and the accretion radius of each binary component is |$R_{\rm acc}=4\, \rm au$|⁠, regardless of which binary parameters are adopted. Particles within this radius are accreted and their mass and angular momentum are added to the star.

3.3 Disc set-up

Each simulation consists of 10⁶ equal mass Lagrangian particles initially distributed from the inner disc radius, r_in, to the outer disc radius, r_out. We consider two values of the inner radius, |$20\, \rm au$| and |$40\, \rm au$|⁠. The latter is farther out than the radius where the tidal torque truncates the disc (Artymowicz & Lubow 1994) – although we note that for a misaligned disc the tidal torque produced by the binary is much weaker, allowing the disc to survive closer to the binary (e.g. Lubow, Martin & Nixon 2015; Miranda & Lai 2015; Nixon & Lubow 2015; Lubow & Martin 2018). The observed outer radius of the gas disc is |${\sim}1300\, \rm au$| (Bi et al. 2020), which means that the majority of the angular momentum is in the outer regions of the disc. In our simulations, we truncate the outer radius to be |$r_{\rm out}=200\, \rm au$| or |$400\, \rm au$| in order to speed up computational time and increase the resolution. We note that this truncated outer radius preserves the angular momentum balance (i.e. the outer disc still holds most of the angular momentum).

The total disc mass is set to |$0.1 \, \rm M_{\odot }$| assuming a dust to gas ratio of 0.01. The value of the disc mass comes from the observations of the dust mass (e.g. Bi et al. 2020). Kraus et al. (2020) estimated a lower disc mass because they did not recover the total flux due to missing short baselines in their ALMA observations. We ignore the effect of self-gravity since it has no effect on the nodal precession rate of flat circumbinary discs and the inferred disc mass is not large enough for self-gravity to be important.

The surface density profile is initially a power-law distribution given by

$$\begin{eqnarray*} \Sigma (R) = \Sigma _0 \bigg (\frac{r}{r_{\rm in}} \bigg)^{-p}, \end{eqnarray*}$$

(1)

where Σ₀ is the density normalization and p is the power-law index. Note that the density normalization is set from the total disc mass above (the simulated total mass is similar to the amount of mass from the three dust rings inferred from Bi et al. 2020, assuming a gas-to-dust ratio of 100). We use a locally isothermal disc with a disc thickness that is scaled with radius as

$$\begin{eqnarray*} H = \frac{c_{\rm s}}{\Omega } \propto r^{3/2-q}, \end{eqnarray*}$$

(2)

where |$\Omega = \sqrt{GM/r^3}$| and c_s is the sound speed. We choose q = 0.5 to ensure that |$H/r = \rm const$| over the radial extent of the disc. We consider two values of the disc aspect ratio |$H/r = 0.05,\, 0.1$| at r = r_in. We take the Shakura & Sunyaev (1973) α to be either 0.01, 0.05, or 0.1. We use the α prescription detailed in Lodato & Price (2010) given as

$$\begin{eqnarray*} \alpha \approx \frac{\alpha _{\rm AV}}{10}\frac{\langle h \rangle }{H}, \end{eqnarray*}$$

(3)

where α_AV is the artificial viscosity and 〈h〉 is the mean smoothing length on particles in a cylindrical ring at a given radius (Price et al. 2018). The viscosity and mean smoothing length are not constant over the disc because we set the disc aspect to be constant. Table A1 shows the minimum and maximum values for the α viscosity parameter and the shelled-averaged smoothing length per scale height 〈h〉/H.

3.4 Adding a planet

We also consider two simulations with a planet that is inclined to the binary orbit but coplanar to the initial disc. The disc has the same surface density profile from equation (1) but we implement a pre-carved gap in the disc in order to prevent excessive accretion of material on to the planet (see e.g. Lubow & Martin 2016; Martin et al. 2016). The inner and outer boundaries of the gap are taken to be |$56\, \rm au$| and |$153\, \rm au$|⁠, which are taken from observations (e.g. Bi et al. 2020; Kraus et al. 2020). The initial semimajor axis of the planet is set to be roughly at the centre of the gap between the inner and middle rings at |${\sim}100\, \rm au$|⁠. We consider a planet mass of |$M_{\rm p} = 1\, \rm M_{J}$|⁠. This mass is sufficient enough to open a gap in the gas (Lin & Papaloizou 1986; Marsh & Mahoney 1992; Nelson et al. 2000). The dynamics are qualitatively the same regardless of planet mass. We keep the viscosity constant across these two simulations (α ≈ 0.01), however, the gap opening is dependent on the viscosity meaning that at lower viscosity planets are able to open gaps easier (e.g. Duffell & MacFadyen 2013). The planet has an initially circular orbit with an accretion radius of |$0.25\, \rm r_{\rm H} = 3.82\, \rm au$| (e.g. Nealon et al. 2018), where r_H is the Hill radius. The pre-carved gap can be seen in the first and third panels in Fig. 9 for the two different disc aspect ratios, 0.1 and 0.05. The initial surface density profile of the pre-carved gap for H/r = 0.05 is given by the black line in the upper panel of Fig. 10.

3.5 CO maps

The hydrodynamical simulations do not include dust grains. Instead we simply assume that the dust distribution contains small, well coupled grains that represents the distribution of the gas. We construct a dust population with a grain-size distribution dn/ds ∝ s^−m between s_min = 0.03 and |$s_{\rm max} = 1000\, \mu \mathrm{ m}$| with m = 3.5. We assume the gas-to-dust ratio value of 100 and calculate the total mass of dust from the gas mass in the simulation. The dust grain opacities are calculated assuming spherical and homogeneous grains (according to Mie theory) and are temperature independent, and we assume astrosilicates composition (Weingartner & Draine 2001). The stars are both assumed to be |$1\, \rm Myr$| old and the mass of the stars comes directly from the simulation. The stellar luminosities are calculated using isochrones from Siess, Dufour & Forestini (2000). We assume that the dust and gas are in thermal equilibrium given that the circumbinary disc is passively heated. We use 10⁸ photon packets on a Voronoi mesh built directly on the particle distribution. To calculate the moment maps, we additionally assume a uniform abundance ratio of |$^{12}\rm CO$|-to-|$\rm H_2$| of 10⁻⁴ and we use |$80\, \rm m\,s^{-1}$| resolution. Consistent with Bi et al. (2020), our synthetic maps are convolved with the ALMA CLEAN beam of 0.122 arcsec × 0.159 arcsec (the beam is shown in the lower left of the synthetic maps).

4 RESULTS

Here, we show the results of the hydrodynamical simulations. First, we compare the disc evolution around a triple star system versus a binary system. Secondly, we simulate a disc around a binary (rather than the triple, see Section 2) and compare the two sets of system parameters from Czekala et al. (2017) and Kraus et al. (2020). Thirdly, we examine in detail the evolution of the tilt, longitude of the ascending node, and surface density profile in a disc without a planet. Fourthly, we introduce a giant planet at 100 au and again consider the disc evolution. Finally, we compare our results to the recent work by Kraus et al. (2020).

When we analyse the SPH simulations, we separate the disc into 300 radial bins that span from the inner-most bound particle to the initial outer disc radius. Within each bin, we calculate the azimuthally averaged surface density, longitude of ascending node, tilt, twist, and eccentricity. The tilt, i, is defined as the angle between the initial angular momentum vector of the binary (the z-axis) and the angular momentum vector of the disc. The twist, ϕ, is measured relative to the x-axis (the initial binary eccentricity vector).

4.1 Three-star hydrodynamical simulations

In this section, we further compare the results of our N-body calculations of modelling a triple star system as a binary star system by using hydrodynamical simulations. We use the stellar and disc parameters from Kraus et al. (2020) from table 2, however, we model our circumtriple disc at a higher resolution (10⁶ Lagrangian particles). Fig. 2 shows the evolution of the disc tilt and longitude of the ascending node as a function of time for a disc viscosity α = 0.01 (bending-wave regime). The disc tilt follows a similar evolution regardless of the stellar prescription. The outer parts of the discs (dotted lines) are closer to each other in tilt than inner parts of the discs (solid lines). However, the tilt changes are only a few degrees.

Figure 2.

Evolution of the inclination i, and longitude of the ascending node ϕ, both as a function of time in units of P_C, where P_C is the orbital period of star C. The disc viscosity is α = 0.01. The black lines represent a circumbinary disc (run6 from Table 3), while the red lines represent a circumtriple disc (run0). The disc is evaluated at two radii, |$45\, \rm au$| (solid) and |$180\, \rm au$| (dotted). The evolution of a circumbinary disc around GW Ori is similar to a circumtriple disc.

We do not see any evidence for the disc breaking due to the tidal torque of the triple star, or the binary star. The precession of the circumtriple disc is similar to that of the circumbinary disc (lower panel in Fig. 2). Fig. 3 compares the surface density, tilt, and twist profiles for the circumbinary and circumtriple discs at three different times. The surface density profiles at all times for both discs are smooth, confirming that neither disc is broken. The warp in the circumtriple disc is consistent with the warp profile in the circumbinary disc. Lastly, the precession rate of the circumtriple disc as a function of radius is slightly slower at later times than the circumbinary disc. This suggests that the circumtriple disc is more stable against breaking than the circumbinary disc for the parameters of these simulations.

Figure 3.

Evolution of the disc surface density Σ (top panel), tilt i (middle panel), and twist ϕ (bottom panel), as a function of radius for a disc viscosity α = 0.01. The black lines represent a circumbinary disc (run6 from Table 3), while the red lines represent a circumtriple disc (run0). We evaluate the disc at three different times, |$t = 0\, P_{\rm C}$| (solid), |$10\, P_{\rm C}$| (dashed), and |$1000\, P_{\rm C}$| (dotted), where P_C is the orbital period of star C. The tilt and twist are only shown in the range from |$25\, {\rm au} \lt r \lt 200\, \rm au$|⁠. There is less mass in the inner parts of the circumtriple disc than the circumbinary disc and the warp profiles are similar for the two discs.

With the above comparison, we are confident that we can model the GW Ori circumtriple disc as a circumbinary disc. With this assumption, we can explore a larger span of the parameter space while increasing computational efficiency. The remaining simulations discussed in this work will be modelling a circumbinary disc.

4.2 Binary parameters from Czekala et al. (2017) versus Kraus et al. (2020)

Recently, Kraus et al. (2020) presented additional high resolution observations of the GW Ori circumtriple system. Before their updated stellar parameters became public, we used the stellar parameters presented in Czekala et al. (2017) to approximate GW Ori as a binary system. Here, we test the difference in the disc evolution between using the binary parameters from Czekala et al. (2017) compared to Kraus et al. (2020).

Fig. 4 shows the disc surface density profile (upper panels) and disc tilt (lower panel) as a function of disc radius for runs 1 and 2 from Table 3. The disc surface density profiles are similar irrespective of the binary parameters used. There is |${\lesssim}20{{\ \rm per\ cent}}$| difference in the tilt evolution and the shape of the tilt profile is similar in both models. Based on this, we conclude that the disc will evolve in a similar fashion when using either set of binary parameters. Thus, unless otherwise indicated, we use the binary parameters given by Czekala et al. (2017).

Figure 4.

The disc surface density profile (Σ, upper panel) and disc tilt (i, lower panel) as a function of disc radius for runs 1 and 2 from Table 3. The black lines correspond to a simulation with the binary parameters of Czekala et al. (2017) and the red lines represent a simulation with the binary parameters of Kraus et al. (2020). The line style denotes the time at which the disc measurements are taken, with the solid, dashed, and dotted corresponding to times t = 0, 1000, and |$3000\, \rm \mathrm{ \mathit{ P}}_{orb}$|⁠, respectively. The surface density profile and tilt evolution show similar structures independent of the binary parameters used.

4.3 Effect of disc parameters

Here, we test the dynamical effects the binary has on two different disc aspect ratios, H/r = 0.1 and 0.05. The upper of these values represents the thicker aspect ratio that has been inferred for GW Ori (Bi et al. 2020) while the lower is a more typical value expected for protoplanetary discs (e.g. D’Alessio et al. 1998).

The left-hand panel of Fig. 5 shows the evolution of the disc inclination and longitude of the ascending node for a disc aspect ratio H/r = 0.1 (run3 from Table 3). We probe the disc at three different radii, 47, 188, and |$337\, \rm au$| corresponding to the centres of the three observed dust rings (Table 1). The tilts of the rings show little warping and a slow alignment towards the binary orbital plane during the simulation. From the evolution of the longitude of the ascending node, the disc shows a slow precession rate. Furthermore, the three measured radii are precessing at roughly the same rate, which is further evidence that no disc breaking has occurred.

Figure 5.

Evolution of the inclination, i, and longitude of the ascending node, ϕ, both as a function of time for two different disc aspect ratios. Left-hand panel: H/r = 0.1 (run3 from Table 3). Right-hand panel: H/r = 0.05 (run1). The disc is evaluated at three radii, |$47\, \rm au$| (black), |$188\, \rm au$| (blue), and |$337\, \rm au$| (red). A thinner disc shows strong warps when compared to a thicker disc.

The right-hand panel of Fig. 5 shows the evolution of the GW Ori disc with a smaller disc aspect ratio H/r = 0.05 (run1 from Table 3). Unlike the thicker disc, the thinner disc is more prone to warping (and breaking since the disc communication time-scale is longer). The disc tilt (upper, right-hand panel) decreases substantially across the entire disc as time increases. The lower sub-panel shows the longitude of the ascending node as a function of time. The three measured radii are precessing in a similar fashion that suggest the disc is not broken but strongly warped. The decrease in the disc tilt in time is caused by the disc aligning to the orbital plane of the binary. Depending on the disc misalignment and binary eccentricity, due to dissipation the disc will evolve to one of two possible alignments. For sufficiently small initial inclination the disc precesses about the binary angular momentum vector and moves towards a coplanar alignment with the binary orbital plane (Papaloizou & Terquem 1995; Lubow & Ogilvie 2000; Nixon et al. 2011; Facchini et al. 2013; Foucart & Lai 2014; Smallwood et al. 2019). The alignment time-scale is very sensitive to the disc aspect ratio, H/r. For small H/r, the alignment time-scale is shorter versus being longer for a more thicker disc (e.g. Lubow & Martin 2018; Smallwood et al. 2020). Therefore, this is why the thicker disc in the left-hand panel of Fig. 5 is not aligning within our time domain, while the thinner disc in the right-hand panel is aligning to the binary orbital plane.

Fig. 6 shows the disc evolution for our simulation with the thinner aspect ratio of H/r = 0.05 (run1 from Table 3). The upper panels show the initial conditions where the disc is tilted by 40^○. The lower panels show the disc evolution at a time |$t = 1000\, \rm \mathrm{ \mathit{ P}}_{\rm orb}$|⁠. The warp in the disc can be clearly seen in the middle panel and the continuous nature of the disc indicates that there is no break present in the disc.

Figure 6.

Disc evolution for a circumbinary disc with i₀ = 40^○ with a disc aspect ratio H/r = 0.05 (run1 from Table 3). Upper panels: initial set-up for the GW Ori disc around an eccentric binary with separation of |$9.2\, \rm au$|⁠. The bottom panels: the disc at a time of |$t = 1000\, \rm \mathit{ P}_{orb}$|⁠. The colour bar denotes the gas density. The left-hand panels show the view looking down on to the binary orbital plane, the x–y plane. The middle panels show the x–z plane and the right-hand panels show the y–z plane. The binary torque causes the disc to become strongly warped but does not break.

To investigate the warping in further detail, in Fig. 7, we show the surface density (top panel), and eccentricity (bottom panel) as function of disc radius. As the disc evolves over time, the inner regions viscously spread inwards and the outer portions outwards. The surface density profile at all times is smooth confirming that the disc is not broken. The bottom panel shows the disc eccentricity as a function of disc radius. The disc initially starts circular, however, there is eccentricity growth that occurs as the warp propagates outwards.

Figure 7.

Surface density Σ (top panel) and eccentricity e (bottom panel) as a function of radius at different times for the H/r = 0.05 simulation (run1 from Table 3) shown in Fig. 4. The black, blue, red, and yellow curves correspond to |$t = 0, 100, 1000, 3000\, \rm \mathit{ P}_{orb}$|⁠, respectively. The disc maintains a smooth surface density profile that is indicative of no disc breaking.

We have demonstrated that the observed parameters in GW Ori do not lead to disc breaking. This motivates us to consider an alternative mechanism to provide a break in the disc, separating the disc into the distinct misaligned planes that are observed.

4.4 Effect of a giant planet

In this section, we carry out SPH simulations with a planet that is initially coplanar with respect to the disc. Both the disc and the planet begin misaligned to the binary orbital plane. We again consider two disc aspect ratios, H/r = 0.1 (run5 from Table 3) and 0.05 (run4), where α = 0.01 in both simulations (such that both discs are in the wave-like regime). Due to the pre-carved gap (described in Section 3.4), initially the inner disc has a mass of |$0.0133\, \rm M_{\odot }$| and the outer disc has |$0.0867\, \rm M_{\odot }$|⁠.

Fig. 8 shows the evolution of the tilt and longitude of the ascending node for disc aspect ratios H/r = 0.1 (left-hand panel) and H/r = 0.05 (right-hand panel). We also include the evolution of the tilt and longitude of the ascending node for the planet. For the thicker disc simulation (left-hand panel), the whole disc remains at an inclination similar to the initial conditions. There is a small deviation in the twist of the disc in the inner parts because the simulation begins with a pre-carved gap that allows the inner disc to freely precess initially. The initial planet evolution is dominated by the binary as it begins to undergo a tilt oscillation due to the binary eccentricity (e.g. Smallwood et al. 2019). Once the planet tilt evolves out of the plane of the disc, the torque from the planet is no longer strong enough to maintain the gap. The large disc aspect ratio means that viscous spreading of the disc occurs on a fast time-scale. Therefore, the gap then quickly fills with gas and the disc evolves rigidly. The disappearance of the gap can be clearly seen in Fig. 9, where the two left-most panels show the initial disc set-up (with H/r = 0.1 and the pre-carved gap centred on |$100\, \rm au$|⁠) and the disc evolution at a time |$t = 230\, \rm \mathrm{ \mathit{ P}}_{planet}$|⁠. The gap has dissipated due to the faster disc spreading coupled with the misalignment of the planet to the disc. The later evolution of the planet is dominated by the interaction of the planet with the disc. Secular planet–disc tilt oscillations occur in discs in binaries (Lubow & Martin 2016; Martin et al. 2016). In circumbinary discs around circular binaries, the planet tends to be closer to the binary orbital plane as a result of planet–disc interactions (Pierens & Nelson 2018).

Figure 8.

Evolution of the inclination, i, and longitude of the ascending node, ϕ, both as a function of time for two different disc aspect ratios. Left-hand panel: H/r = 0.1 (run5 from Table 3). Right-hand panel: H/r = 0.05 (run4). The disc is evaluated at three radii, |$47\, \rm au$| (black), |$188\, \rm au$| (blue), and |$337\, \rm au$| (red). The planet is given by the yellow lines.

Figure 9.

Disc evolution for a circumbinary disc with i₀ = 40^○ along with a circumbinary planet with |$i_{0,\rm p} = 40^{\circ }$|⁠. We show the results for two different aspect ratios H/r = 0.1 (two left-most panels, run5 from Table 3) and 0.05 (two right-most panels, run4). The first and third panel beginning from the left shows the evolution at a time of |$t = 1\, \rm \mathit{ P}_{planet}$| with the pre-carved gap for H/r = 0.1 and H/r = 0.05. The second and fourth panels from the left show the disc evolution at a time of |$t = 230\, \rm \mathit{ P}_{planet}$| for the two different disc aspect ratios. The colour bar denotes the gas density. We show the view in the y–z plane. A planet is able to carve a gap within a thin disc.

The right-hand panel of Fig. 8 shows the evolution for a thinner disc with H/r = 0.05. The lower disc aspect ratio means it is easier for the planet to create and hold open a gap. Similar to the thick disc case, the planet becomes misaligned to the plane of the disc. The evolution of the inner ring, centred at |$47\, \rm au$| (black lines), is dominated by tilt oscillations that are primarily driven by the outer disc. The planet and the inner disc ring undergo tilt oscillations that are driven by the outer disc. Both the inner ring and planet undergo similar evolution but on different time-scales. The evolution of both is dominated by the massive outer disc part since their inclinations are lower than the outer disc (e.g. Pierens & Nelson 2018). The middle (⁠|$188\, \rm au$|⁠) and outer (⁠|$337\, \rm au$|⁠) portions of the disc show evolutionary changes (decrease in inclination), due to binary-disc alignment. Had we not initially truncated the outer radius we would expect their inclinations to remain more constant. Since this simulation has a lower disc aspect ratio than the simulation described in the previous paragraph, the viscous spreading of the disc into the planet gap is slower that allows the planet to break the disc. The two right-most panels in Fig. 9 show the initial disc set-up (with H/r = 0.05 and the pre-carved gap centred at |$100\, \rm au$|⁠) and the disc evolution at a time |$t = 230\, \rm \mathit{ P}_{planet}$|⁠. The formation of a |$1\, \rm \mathit{ M}_{\rm J}$| planet in the disc is able to maintain a long lived strongly warped disc structure when h/r < 0.05.

We further examine the evolution of the disc with H/r = 0.05 (run4) at times |$t = 0\, , 100\, , 180\, , 230\, \rm \mathit{ P}_{planet}$| in Fig 10. There is initially a pre-carved gap in the surface density profile, shown by the trough in the curve centred on |$100\, \rm au$|⁠. Since we start with an initial gap, the disc is essentially broken and this break propagates outwards to about |$150\, \rm au$| at a time |$t = 100\, \rm \mathrm{ \mathit{ P}}_{\rm planet}$|⁠. At |$t = 180\, \rm \mathit{ P}_{\rm planet}$| the break has propagated to the outer regions of the disc and slowly dissipates. After |$t = 180\, \rm \mathit{ P}_{\rm planet}$|⁠, the initial break has been fully dissipated, however, the |$1\, \rm M_{\rm J}$| planet is close to coplanar to the disc and it begins to open a new gap which is shown by the dip in the surface density profile at |$t = 230\, \rm \mathrm{ \mathit{ P}}_{\rm planet}$|⁠. During the simulation, the planet migrates inward to |${\sim}75\, \rm au$| by a time of |$t = 230\, \rm \mathrm{ \mathit{ P}}_{\rm planet}$|⁠. The centre of the break is located at |${\sim}75\, \rm au$|⁠. The peaks in the eccentricity profile correspond to the break location. At |$t = 230\, \rm \mathit{ P}_{\rm planet}$| the inner regions of the disc have a larger eccentricity than the outer parts of the disc. From observations, the inner ring is more eccentric than the middle and outer rings, with an estimated eccentricity of ∼0.2 (Bi et al. 2020; Kraus et al. 2020). Moreover, Fig. 11 shows the planet mass as a function of time for an initially |$1\, \rm \mathit{ M}_{\rm jup}$| planet. The planet mass increases significantly after |$t = 200\, \rm \mathit{ P}_{planet}$|⁠, which corresponds to the time the planet has realigned with the disc.

Figure 10.

Beginning from the top we show the surface density Σ, disc tilt i, longitude of the ascending node ϕ, and eccentricity e as a function of radius at different times. The black, blue, red, and yellow curves correspond to |$t = 0, 100, 180, 230\, \rm \mathit{ P}_{planet}$|⁠, respectively. The initial conditions for the circumbinary disc are for run4 from Table 3, which has a disc aspect ratio H/r = 0.05 and pre-carved gap. A planet is able to maintain a gap within the disc seen by the dips in the surface density profile.

Figure 11.

The planet mass, M_p, as a function of time in planet orbital periods, P_planet from run5 from Table 3. The planet mass increases drastically after |$t \gt 200\, \rm \mathit{ P}_{planet}$|⁠.

4.5 CO kinematics

As we have recovered the disc structure inferred by Bi et al. (2020), we now consider the synthetic CO moment 1. Fig. 12 shows the comparison between the ¹²CO J = 1–2 moment 1 maps of our simulations. The upper panel reproduces the observation from Bi et al. (2020) for convenience. For a regular Keplerian disc, the green regions should represent a well-defined ‘butterfly-like’ pattern. However, the observations show a twisted pattern inside ∼0.2 arcsec, which can be accounted for by a warp in the disc where there is a misalignment between the inner and outer portions. The twist is outlined in the insert. This twisted pattern has also been seen in the discs around HD 142527 (Casassus et al. 2015; Marino, Perez & Casassus 2015), HD 143006 (Benisty et al. 2018; Pérez et al. 2018), HD 97048 (van der Plas et al. 2017), AB Aurigae (Poblete et al. 2020), IRS 48 (Calcino et al. 2019), and MWC 758 (Calcino et al. 2020). The lower-left panel shows the CO velocity map for the simulation of only a circumbinary disc (run3 from Table 3) and the lower-right panel in Fig. 12 shows the CO velocity map for our simulation that includes a planet (run4).

Figure 12.

The comparison between the ¹²CO J = 1–2 moment 1 maps of the observations from Bi et al. (2020) (upper panel) and our two hydrodynamical simulations (lower-left panel, run3, and lower-right panel, run4). The observation and synthetic images are performed with a beam of 0.122 arcsec × 0.159 arcsec with a position angle of −32.3^○ (bottom left corner). The dot–dashed line highlights the shape of the twist. The velocities in the upper and lower panel are shown with respect to the rest velocity of 13.5 km s⁻¹. We plot the binary (white stars) and planet (white dot) in the synthetic images. There are no localized artefacts surrounding the planet. We have copied the twist line in the observation panel and displayed it on our synthetic images.

The twist in the inner regions of the synthetic image with a planet is more consistent with the observation.

In both synthetic images, the inner portions of the map, ∼0.2 arcsec, shows a twist in inner regions which is consistent with the misaligned inner binary. While the broad features in both maps from the simulations are roughly in agreement with the features identified in the observations, only the simulation with the planet additionally has a matching disc structure.

4.6 Comparison with Kraus et al. (2020)

Up to this point we have shown, from our SPH simulations modelling a circumtriple and circumbinary disc, that the disc does not break due to the binary or triple star system. This motivated our previous section, where we invoked a planet to break the disc at the required radius. However, Kraus et al. (2020) conducted SPH simulations and found that the torque from the triple system can effectively break the disc. Here, we explore the differences between our model and the model from Kraus et al. (2020).

Kraus et al. (2020) also conduct their simulation using gas particles in a Lagrangian SPH code (Bate, Bonnell & Price 1995; Price 2007). They use fewer particles over a smaller radial extent, with 8 × 10⁵ between |$20\, \rm au$| and |$200\, \rm au$| where we use 1 × 10⁶ particles between |$40\, \rm au$| and |$400\, \rm au$|⁠. Their disc surface density profile is shallower than ours with Σ ∝ r^−1/2. Both of our models use a fixed disc aspect ratio of H/r = 0.05. One difference of their set-up compared to this work is that their disc is initially set-up orbiting a single mass of |$5.26\, \rm M_{\odot }$|⁠, which is the mass of the total triple star system. The disc is then simulated with this single gravitational mass so that any transient features due to the initial conditions are dissipated. Once the disc reaches a steady state, the central mass is replaced with the three stars and the disc is reoriented on the centre of mass of the system and inclined by 38^○ relative to the plane of the (AB) – C orbit. Along with this, any material within |$40\, \rm au$| is removed. They also neglect disc self-gravity and the gravitational effect from the disc on to the stars is ignored. The results of their simulation showed an inner ring break off from the outer disc and precess independently. Moreover, the eccentricity measured in the inner ring of their simulation is ∼0.15, the eccentricity measured from observations of the innermost ring of ∼0.2.

Since we have provided adequate evidence that the binary set-up and parameters do not cause the discrepancy in the disc evolution between the two models, we now explore the effects of varying the disc parameters. The simulations discussed here are runs 6–10 from Table 3. Here, we repeat the simulation set-up from Kraus et al. (2020) but we use a binary system rather than a triple system (we tested the triple system back in Section 4.1) and do not clear particles within |$40\, \rm au$|⁠. We explore how sensitive the α parameter is to induce disc breaking. Furthermore, Kraus et al. (2020) initially clear any particles within |$40\, \rm au$|⁠, however, we also investigate how tuned the disc breaking is to the initial inner disc radius. Our effort in exploring how these parameters contribute to disc breaking will resolve the discrepancies in the competing models of Bi et al. (2020) and Kraus et al. (2020).

4.6.1 Viscosity

To better visualize the breaking when varying α, we show the disc structure in Fig. 13. Beginning from the left-most panel, we show the initial condition, followed by the disc structure at a time |$t=500\, P_{\rm orb}$| for the simulations with α = 0.01 (run6 from Table 3), α = 0.05 (run7), and α = 0.1 (run8). The disc is initially misaligned by 38^○ and we view each disc in the x–y plane. The two discs that are in the diffusive regime (the two right-most panels) show the disc breaking. The disc that is in the bending wave regime, α = 0.01, has no signs of disc breaking. In the wave-like regime the communication throughout the disc is rapid and allows the disc to maintain a coherent disc-like structure while the diffusive disc cases break. The breaking criteria in hydrodynamical simulations are also dependent on resolution; however, the number of initial particles used here provides adequate resolution (comparing the values in our Table 3 with fig. 8 in Nealon, Price & Nixon 2015).

Figure 13.

Disc evolution for a circumbinary disc with varying the α–viscosity parameter. The binary components are shown by the red dots. Beginning from the left-most panel we show the initial conditions, α = 0.01 (run6 from Table 3), α = 0.05 (run7), and then α = 0.1 (run8). The disc evolution is shown at a time |$t = 500\, \rm P_{orb}$|⁠. The colour bar denotes the gas density. We show the view looking down on to the binary orbital plane, the x–y plane. For higher viscosity values, the disc is more prone to breaking.

Fig. 14 shows the surface density profile for discs with three different α–viscosity parameters. The lower viscosity value provides a smooth surface density profile which means that the disc has not broken. For |$\alpha = 0.05,\, 0.1$| (diffusive-type discs), there is a dip in the surface density profile which means that the disc is broken. For the lower α and H/r expected in GW Ori, Fig. 13 confirms that although warping is expected, breaking is not.

Figure 14.

The surface density profile taken from run6 (α = 0.01, black), run7 (α = 0.05, blue), and run8 (α = 0.1, red) at a time |$t = 500\, \rm P_{orb}$|⁠. The wave-like regime is shown by the solid line and the diffusive regime is denoted by the dotted lines. At lower viscosities more typical of protoplanetary discs, the smooth surface density profile shows no sign of breaking.

In the viscous regime, the dominant torques are the viscous torque and the precession torque. Therefore, for high α values the viscous torque dominates the precession torque, meaning that the disc is less likely to break (e.g. Nixon et al. 2013). For α = 0.05, the disc breaks more easily than when α = 0.1. The dip in the surface density is more prominent in the α = 0.05 case than when α = 0.1 because, in the diffusive regime, for high α values the disc becomes comparatively harder to break (Nixon et al. 2013).

In the wave-like regime, the breaking criteria depend upon the global precession rate compared to the communication time-scale (as described in Section 3.1). Neither of these time-scales depends directly upon α (see equations 3 and 8). However, as shown in Fig. 14, the inner edge of the disc is further out for lower α and for larger α values the inner edge can live closer to the binary. The global precession time-scale is sensitive to the value of the disc inner radius. The precession time-scale increases with r_in. Thus, in the wave-like regime, for lower α, r_in increases, the precession time-scale therefore increases and the disc is less likely to break. Since α in the GW Ori disc may be even lower than the 0.01 value we have considered, the disc may be even less susceptible to breaking that in the simulation shown here.

4.6.2 Initial location of the inner radius

Fig. 15 shows the surface density profiles of three discs with initial inner radii of |$40\, \rm au$| (black line, run6 from Table 3), |$30\, \rm au$| (purple line, run9), and |$20\, \rm au$| (yellow line, run10). The surface density profile of the larger inner radius simulation remains smooth which indicates that the disc is not broken. This disc structure is the same for when |$r_{\rm in} = 30\, \rm au$|⁠. Meanwhile, the disc with a smaller initial inner radius, |$20\, \rm au$|⁠, has a dip in the surface density profile. Fig. 16 shows the disc structure given two different initial radii (⁠|$40\, \rm au$| and |$20\, \rm au$|⁠). A clear break can be seen for when |$r_{\rm in} = 20\, \rm au$|⁠, while the disc structure with |$r_{\rm in} = 40\, \rm au$| remains smooth. Whether the disc breaks or not thus depends sensitively on the inner radius that the disc is initialized with.

Figure 15.

The surface density profile taken from run6 (⁠|$r_{\rm in} = 40\, \rm au$|⁠, black), run9 (⁠|$r_{\rm in} = 30\, \rm au$|⁠, purple), and run10 (⁠|$r_{\rm in} = 20\, \rm au$|⁠, yellow) at a time |$t = 1000\, \rm \mathit{ P}_{orb}$|⁠. Disc material close to the binary will cause the disc to break, showing a deep depression in the surface density profile.

Figure 16.

Disc evolution for a circumbinary disc with varying the initial inner disc radius. We show |$r_{\rm in} = 40\, \rm au$| (run6, left-hand panel) and |$r_{\rm in} = 20\, \rm au$| (run10, right-hand panel). The binary components are shown by the red dots. The disc evolution is shown at a time |$t = 1000\, \rm \mathit{ P}_{orb}$|⁠. The colour bar denotes the gas density. We show the view looking down on to the binary orbital plane, the x–y plane. Disc breaking occurs when there is too much material initially close to the binary.

5 DISCUSSION

5.1 Growth of an inclined planet

The results of the simulations described above show that if a planet forms in a misaligned disc and is massive enough to carve a gap, it can lead to an effectively broken disc. The inner disc precesses faster than the viscous spreading (shown in the third panel in Fig. 10) and so the disc parts remain misaligned, as seen in the second panel. This process may repeat itself – if the planet becomes misaligned again, the break will propagate outward with the warp and dissipate until the planet’s tilt oscillates back to a coplanar orientation with respect to the disc. The planet will then carve another gap, and so on. Each time the planet becomes aligned with the disc, the mass of the planet increases. From Fig. 11, the planet mass increases significantly after |$t = 200\, \rm \mathit{ P}_{planet}$|⁠, which corresponds to the time the planet has realigned with the disc. This implies that planets formed in a misaligned disc may become more massive than planets formed within a coplanar disc if they are able to carve multiple gaps in the disc due to their evolution where they become misaligned to the disc and later realigned.

The proposed inclined planet in GW Ori may be difficult to observe. Planets are large separations, regardless of their planetary radius, are more difficult to detect than giant planets at small separations. Giant planet detections are more common around A stars and that wide-orbit planets are more conducive around high-mass stars (Johnson et al. 2010; Reffert et al. 2015). Moreover, the occurrence rate of directly imaged giant planets is of the order of |$10{{\ \rm per\ cent}}$| (Galicher et al. 2016; Meshkat et al. 2017; Baron et al. 2019; Nielsen et al. 2019). This being said, the results displayed in the CO kinematics in Fig. 12 show no localized artefacts surrounding the giant planet, which suggest that detection would be challenging.

Furthermore, Kraus et al. (2020) presented SPHERE and GPI coronagraphic-polarimetric observations of GW Ori. They are best suited to reveal disc structures by exploiting the fact that direct starlight is not polarized but scattered light from the disc is. These observations are not ideal for searching for thermal emission from faint companions next to bright stars. High contrast imaging in total intensity employing various kinds of speckle suppression techniques (e.g. ADI, Marois et al. 2006) are needed for this purpose. In addition, searching for companions in discs faces the difficulty that disc signatures may compromise point source recoveries (Maire et al. 2017). So far, the only wildly accepted planet detections in discs are PDS 70b,c (Keppler et al. 2018; Haffert et al. 2019). The small stellocentric distance of the predicted gap opening planet in our model (∼0.25 arcsec) also poses a challenge. Additional challenges include the multiplicity of the central stars (potentially complicating coronagraph deployment) and their high luminosities (bolometric luminosity 50 L_⊙; Fang et al. 2014) that results in large flux ratio between the stars and planets.

5.2 Connecting the disc breaking and the dust structures

In this work, we exclusively simulated a gas disc. However, the observations by Bi et al. (2020) are of the dust in the gas disc. Dust particles undergo various degrees of coupling depending on their Stokes number (e.g. Birnstiel, Dullemond & Brauer 2010). Initially, well-coupled dust grains grow, over time, to higher Stokes number and gradually decouple from the gas disc. If significant decoupling occurs in a misaligned disc, the dust particle orbits may evolve independently of the gas disc due to differential precession and the dust structure will not maintain its coherent structure (e.g. Nesvold et al. 2016; Aly & Lodato 2020). As the disc rings in GW Ori are observed as coherent structures, the dust must be well coupled to the gas, justifying our use of gas-only simulations to infer the observed structures. We assume that there are two distinct rings in the GW Ori system, ring 1 and ring 23. Individual rings 2 and 3 have similar inclination and phase angle, which suggest that they are not broken but are instead mildly warped. An additional planet located at |$100\, \rm au$| is able to explain the misalignment between rings 1 and 23. Such a low-mass planet may be able to explain the warping in the outer ring 23, but we note there are also alternative scenarios that do not require planets that could explain the separation between the outer dust rings (e.g. Flock et al. 2015; Dullemond et al. 2018; Suriano et al. 2018, 2019; Riols & Lesur 2019; Tominaga, Takahashi & Inutsuka 2020).

5.3 Viscosity

In Section 5.4, we showed that our conclusions rely on a robust estimate for the viscosity in the disc (parametrized by α) and the aspect ratio (H/r). The viscosity is the fundamental process of how accretion discs evolve (Lynden-Bell & Pringle 1974; Pringle 1981). It determines the transport of mass and angular momentum within the disc, which in turn provides how much energy is released. Viscosity transports angular momentum outwards, allowing matter to spiral inwards in a disc. The disc viscosity parameter, α, can be estimated from observations. The simple estimate of α comes from comparing estimated disc masses, M_d, with estimated central accretion rates, |$\dot{M}_{\rm c}$|⁠, and from these deducing an accretion time-scale |$\tau _{\nu } \sim M_{\rm d}/\dot{M}_{\rm c}$| (e.g. Lodato et al. 2017; Martin et al. 2019). By measuring the sound speed c_s and the disc height H in the outer disc, the viscosity value can be determined (Hartmann et al. 1998).

In the past, the upper limit for the α value for protostellar discs was estimated to be α ∼ 0.01 on distance scales 10–|$100\, \rm au$| (Hartmann et al. 1998; Hartmann 2000; Trapman et al. 2020), but recent observations show that the upper limit is closer to α ∼ 0.001. Andrews et al. (2009) observed protoplanetary discs in Ophiuchus and found α ∼ 0.0005–0.08 for radius |$R = 10\, \rm au$|⁠. Hueso & Guillot (2005) found that 0.001 < α < 0.1 for DM Tau and 4 × 10⁻⁴ < α < 0.04 for GM Tau. More recently, Rafikov (2016) used a self-similar disc solution to calculate 0.0001 < α < 0.04 for resolved disc observations by ALMA. Ansdell et al. (2018) then refined these calculations by measuring the gas disc size and found 0.0003 < α < 0.09. Pinte et al. (2016) measured the dust scale height in HL Tau, and estimated the turbulent viscosity coefficient to be a few 10⁻⁴. The turbulence levels in discs using ALMA gas observations had been directly measured giving α ∼ 0.001 (e.g. Flaherty et al. 2015, 2017; Teague et al. 2018). Moreover, there is evidence that the characteristic ‘double gaps’ in HL Tau, TW Hya, and HD 169142 are produced by a low-mass planet, which requires a low-disc viscosity (e.g. Dong et al. 2017; Dong, Najita & Brittain 2018). From observations of protoplanetary discs, the disc is certainly expected to be in the bending-wave regime rather than the viscous regime. In the context of Sections 4.1 and 4.3 and Fig. 14, our results confirm that the gaps in the circumtriple disc around GW Ori are not produced by the binary (or triple star) torque.

5.4 Inner disc radius

We simulated three different initial disc radii, |$40,\, 30,\, \mathrm{ and}\, 20\, \rm au$|⁠, and showed that the binary torque is able to break the disc for |$r_{\rm in} = 20\, \rm au$| but not for |$r_{\rm in} = 40\, \rm au$| or |$r_{\rm in} = 30\, \rm au$|⁠. Material within our simulations is free to move inwards and so even though it begins at 30 or |$40\, \rm au$|⁠, it moves in closer to a radius that is determined by the balance of the tidal and viscous torques (e.g. Lubow et al. 2015; Miranda & Lai 2015; Franchini, Lubow & Martin 2019b). When the inner radius is smaller, there is much more material in the inner regions of the disc and the precession time-scale becomes smaller than the radial communication time-scale (Lubow & Martin 2018). The location of the inner radius is thus crucial to determining whether the disc breaks or not and should be considered carefully in future work. The simulation with initial |$r_{\rm in}=20\, \rm au$| places more material closer in than there would be in the quasi-steady state. This can only be achieved if the accretion of material is occurring in r < 30. In the case of GW Ori, observations suggest that the inner radius is located at ∼32 au. Bi et al. (2020) adopted an inner disc radius of |$32\, \rm au$|⁠, which is three to four times the AB–C binary semimajor axis, which is also supported by Czekala et al. (2017) and Kraus et al. (2020). With this larger inner radius, our simulations confirm there should be no breaking.

Formation scenarios suggest that misaligned discs are likely to be formed by misaligned gas falling on to existing binaries or triples (Bate 2018). The large outer radius of GW Ori (⁠|${\sim}1300\, \rm au$|⁠; Bi et al. 2020) and orientation of the outer disc are consistent with this formation mechanism. In this interpretation, as the gas and dust fall inwards, it crosses the region where warping and breaking can occur. As such, we do not predict that material should be found inside the warping radius with the observed orientation of the outer disc. Although not the focus of this study, this is highly relevant to the inner radius used in the initial conditions, which is used here and shared with Kraus et al. (2020).

5.5 Consideration of planets

A possible mechanism to produce gaps and misalignment in the GW Ori disc may be due to the presence of planets. High-mass planets exert a tidal torque that overpowers the local viscous torque and form a gap in the gas (Papaloizou & Lin 1984; Bryden et al. 1999). In Section 4.4, we showed that a giant planet that forms initially coplanar to the disc can produce a strong warp, if the disc is sufficiently thin. Therefore, the misalignment between the inner and middle rings in GW Ori may be caused by the presence of planets. However, if the disc aspect ratio is larger, the presence of the dust gaps in GW Ori must be produced by a low-mass planet that is well coupled to the gas disc. The low-mass planet will be well coupled to the gas forcing it to precess at the same rate as the disc. Dipierro et al. (2016) ran SPH simulations of a gas and dust disc with an embedded low-mass planet. The planet is effective in opening a gap in the dust but not in the gas. However, this would not explain the observed misalignment between the inner and middle rings since the gas disc would maintain a flat coherent structure.

If misaligned planets are present around the hierarchical triple star system, they would be difficult to detect. The torque produced by the binary affects the formation processes of planets embedded in the circumbinary gas disc compared to discs around single-star systems (Martin et al. 2014; Fu, Lubow & Martin 2015a,b, 2017). When a giant planet is formed within a misaligned disc, the torque from the binary prevents the planet from remaining coplanar to the binary orbital plane (Lubow & Martin 2016; Martin et al. 2016; Pierens & Nelson 2018). The probability of detecting inclined planets through the transit method is lower than coplanar planets. Follow-up observations have revealed that |${\sim}2.5{{\ \rm per\ cent}}$| of planets are in triple and multiple systems (Roell et al. 2012; Fragione et al. 2019). However, no planet in a circumtriple orbit has been detected. If a planet (or planets) is the cause of the dust gaps in the circumtriple disc around GW Ori, then they would be the first circumtriple planet(s).

6 CONCLUSIONS

We have examined the origin of the coherent dust structures around the GW Ori hierarchical triple system. Bi et al. (2020) first suggested that the break cannot be caused by the torque from the observed triple star system. More recently, Kraus et al. (2020) conducted SPH simulations and stated that the triple star torque was the cause for the break. In this work, we carried out extensive SPH simulations to explore the discrepancy between the two models.

First, we compared through 3D hydrodynamical simulations that a circumbinary disc evolves in a similar fashion as a circumtriple disc in the context of the GW Ori system. Secondly, we tested the differences in the binary parameters between the two models and found that the disc evolves in a similar fashion, independently of the binary parameters used. We then examined the disc viscosity and the inner radius of the disc since these two parameters heavily impact the criteria for disc breaking. We found that when α = 0.01, the disc is strongly warped but does not break. Lastly, we showed a small initial disc radius will cause the disc to break even in the bending-wave regime. However, this effect is due to the initial conditions of the simulation. Since α = 0.01 is an upper limit for protoplanetary discs (Hartmann et al. 1998) and the surface density profile is tapered within the inner regions, our results show that the break in the GW Ori circumtriple disc is not caused by the triple star system.

We present an alternative scenario to explain the origin of the dust rings in GW Ori, using a planet (or planets). We find that an initially massive planet can continuously open a gap within a thin disc as the planet’s tilt oscillates in and out of the disc plane. For a thicker disc, the viscous spreading is too fast for the planet to maintain the gap. However, a low-mass planet that is well coupled to the gas can still open a gap in the dust (e.g. Dipierro et al. 2016). In conclusion, we have shown that the break in the GW Ori circumtriple disc is not due to the torque imposed on to the disc by the stars. Therefore, the disc breaking must be caused by undetected planets, which would be the first planets in a circumtriple orbit.

ACKNOWLEDGEMENTS

We thank the anonymous referee for helpful suggestions that positively impacted the work. The authors thank Alison Young for providing details of her simulations. The three-body integrations in this work made use of the rebound code that can be downloaded freely from http://github.com/hannorein/rebound. Computer support was provided by UNLV’s National Supercomputing Center. We acknowledge the use of splash (Price 2007) for the rendering of Figs 6, 9, 13, and 16 as well as pymcfost for Fig. 12. RN acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 681601) and subsequent support from UKRI/EPSRC through a Stephen Hawking Fellowship (EP/T017287/1). RD acknowledges the financial support provided by the Natural Sciences and Engineering Research Council of Canada. CP acknowledges funding from the Australian Research Council via FT170100040 and DP180104235.

DATA AVAILABILITY

The data supporting the plots within this article are available on reasonable request to the corresponding author. A public version of the phantom and splash codes is available at https://github.com/danieljprice/phantom and http://users.monash.edu.au/~dprice/splash/download.html, respectively. mcfost is available on request.

APPENDIX A: SUPPLEMENTAL INFORMATION

In this appendix, we provide additional information on the set-up of the hydrodynamical simulations by showing the ranges of the α–viscosity parameter, and the shelled-averaged smoothing length per scale height in Table A1.