Inclined Massive Planets in a Protoplanetary Disc: Gap Opening, Disc Breaking, and Observational Signatures

TL;DR

This study uses 3D hydrodynamical simulations and radiative transfer calculations to explore how massive inclined planets can cause disc breaking, misalignment, and observable signatures like shadows and gaps, aiding in planet detection.

Contribution

It provides new analytical formulas and simulation results linking planet mass, gap properties, and disc breaking conditions, advancing understanding of planet-disc interactions.

Findings

Massive planets induce deep gaps and disc breaking.

Disc misalignment leads to observable shadows and kinematic signatures.

Deep gaseous gaps are associated with disc shadowing in observed systems.

Abstract

We carry out three-dimensional hydrodynamical simulations to study planet-disc interactions for inclined high mass planets, focusing on the disc's secular evolution induced by the planet. We find that, when the planet is massive enough and the induced gap is deep enough, the disc inside the planet's orbit breaks from the outer disc. The inner and outer discs precess around the system's total angular momentum vector independently at different precession rates, which causes significant disc misalignment. We derive the analytical formulae, which are also verified numerically, for: 1) the relationship between the planet mass and the depth/width of the induced gap, 2) the migration and inclination damping rates for massive inclined planets, and 3) the condition under which the inner and outer discs can break and undergo differential precession. Then, we carry out Monte-Carlo radiative transfer calculations for the simulated broken discs. Both disc shadowing in near-IR images and gas kinematics probed by molecular lines (e.g. from ALMA) can reveal the misaligned inner disc. The relationship between the rotation rate of the disc shadow and the precession rate of the inner disc is also provided. Using our disc breaking condition, we conclude that the disc shadowing due to misaligned discs should be accompanied by deep gaseous gaps (e.g. in Pre/Transitional discs). This scenario naturally explains both the disc shadowing and deep gaps in several systems (e.g. HD 100453, DoAr 44, AA Tau, HD 143006) and these systems should be the prime targets for searching young massive planets ($>M_J$) in discs.

Equations

$$T_{wave}\sim q^2{\left(\frac{H}{R}\right)}_{R_p}^{-3}{\Sigma }_0{\Omega }_p^2{R}_p^4,$$

$$T_{wave}\sim q^2{\left(\frac{H}{R}\right)}_{R_p}^{-3}{\Sigma }_{gap}{\Omega }_p^2{R}_p^4.$$

$$T_{\nu }\sim {\Sigma }_0\nu {\Omega }_p{R}_p^2,$$

$$\frac{{\Sigma }_{gap}}{{\Sigma }_0}\sim \alpha h^5{q}^{-2}.$$

$$\frac{{\Sigma }_{gap}}{{\Sigma }_0}=\frac{1}{1+0.04K},$$

$$K={\alpha }^{-1}{h}^{-5}{q}^2,$$

$$\frac{{\Delta }_{gap}}{R_p}=0.41K^{\prime 1/4},$$

$$K^{\prime }=K{h}^2={\alpha }^{-1}{h}^{-3}{q}^2.$$

$$\Phi =\frac{\int \frac{G{M}_p}{{\left((R_p\sin \varphi \sin i_p{)}^2+H^2\right)}^{1/2}}d\varphi }{\int d\varphi },$$

$$q\to q\frac{\int \frac{1}{\sqrt{(R_p\sin \varphi \sin (i_p)/H{)}^2+1}}d\varphi }{\int d\varphi }.$$

$$q\to q\frac{2}{\pi \sqrt{n+1}}\mathcal{K}(\frac{n}{n+1})$$

$$\frac{{\Sigma }_{gap}}{{\Sigma }_0}=\frac{1}{1+0.04K},$$

$$\frac{{\Delta }_{gap}}{R_p}=0.41K^{\prime 1/4},$$

$$K$$

$$K^{\prime }$$

$$-\frac{⟨{\partial }_t{R}_p⟩}{R_p}=\min \{\begin{array}{l}(2.7+1.1{\alpha }_s)\cdot t_{mig}^{-1}\\ 8.8\cdot h^2\cdot {\sin }^{-1}(i_p/2)\cdot {\sin }^{-1}(i_p)\cdot t_{mig}^{-1}\end{array}$$

$$t_{mig}={\Omega }_p^{-1}{h}^2{q}^{-1}{\left(\frac{{\Sigma }_0{R}_p^2}{M_{\ast }}\right)}^{-1}$$

$$-\frac{⟨{\partial }_t{i}_p⟩}{i_p}=\min \{\begin{array}{l}0.544\cdot t_{inc}^{-1}\\ 1.46\cdot h^4\cdot {\sin }^{-3}(i_p/2)\cdot i_p^{-1}\cdot t_{inc}^{-1}\end{array}$$

$$t_{inc}={\Omega }_p^{-1}{h}^4{q}^{-1}{\left(\frac{{\Sigma }_0{R}_p^2}{M_{\ast }}\right)}^{-1}$$

$$-\frac{⟨{\partial }_t{R}_p⟩}{R_p}=\min \{\begin{array}{l}(2.7+1.1{\alpha }_s)\cdot t_{mig}^{-1}\\ 8.8\cdot h^2\cdot {\sin }^{-1}(i_p/2)\cdot {\sin }^{-1}(i_p)\cdot t_{mig}^{-1}\end{array}$$

$$t_{mig}={\Omega }_p^{-1}{h}^2{q}^{-1}{\left(\frac{{\Sigma }_0{\left(\frac{{\Sigma }_{gap}}{{\Sigma }_0}\right)}^{qs}{R}_p^2}{M_{\ast }}\right)}^{-1}.$$

$$-\frac{⟨{\partial }_t{i}_p⟩}{i_p}=\min \{\begin{array}{l}0.544\cdot t_{inc}^{-1}\\ 1.46\cdot h^4\cdot {\sin }^{-3}(i_p/2)\cdot i_p^{-1}\cdot t_{inc}^{-1}\end{array}$$

$$t_{inc}={\Omega }_p^{-1}{h}^4{q}^{-1}{\left(\frac{{\Sigma }_0{\left(\frac{{\Sigma }_{gap}}{{\Sigma }_0}\right)}^{qs}{R}_p^2}{M_{\ast }}\right)}^{-1}.$$

$$\Sigma R^2\Omega \frac{\partial \mathbf{l}}{\partial t}=\frac{1}{R}\frac{\partial \mathbf{G}}{\partial R}+\mathbf{T}$$

$$\frac{\partial \mathbf{G}}{\partial t}-\left(\frac{{\Omega }^2-{\kappa }^2}{2\Omega }\right)\mathbf{l}\times \mathbf{G}+\alpha \Omega \mathbf{G}=\frac{P{R}^3\Omega }{4}\frac{\partial \mathbf{l}}{\partial R},$$

$$\int \Sigma R^2\Omega \frac{\partial \mathbf{l}}{\partial t}RdRd\varphi ={\int \mathbf{G}|}_{R_{in}}^{R_{out}}d\varphi +\int \mathbf{T}RdRd\varphi .$$

$$T(R)=-\frac{3}{2}\frac{G{M}_p\Sigma R^2}{R_p^3}\sin (i_p)\cos (i_p){\sin }^2\varphi$$

$${\omega }_{rb}=-\frac{3}{4}\frac{G{M}_p}{R_p^3}\cos (i_p)\frac{\int \Sigma R^{3}dR}{\int \Sigma R^3\Omega dR}.$$

$${\omega }_{rb}$$

$$\frac{{\Omega }_d}{{\Omega }_p}={\left(\frac{M_{\ast }}{M_{\ast }+M_p}\right)}^{1/2}{\left(\frac{R_p}{R_d}\right)}^{3/2},$$

$$\frac{{\omega }_{rb}}{{\Omega }_p}=\frac{3}{8}\cos (i_p)q{\left(\frac{1}{1+q}\right)}^{1/2}{\left(\frac{R_d}{R_p}\right)}^{3/2}.$$

$$\frac{{\omega }_{rb}}{{\Omega }_p}\sim \frac{3}{8}\cos (i_p)q,$$

Full text

††pagerange: Inclined Massive Planets in a Protoplanetary Disc: Gap Opening, Disc Breaking, and Observational Signatures–References††pubyear: 2018

1 Introduction

The planet-disc interaction theory has been developed over decades (Goldreich & Tremaine, 1979; Tanaka et al., 2002; Kley & Nelson, 2012; Baruteau et al., 2014). But most studies focus on planets that are coplanar with the disc. This is due to both the fact that planets in our Solar Systems are coplanar and the simplicity of the coplanar case where the vertically averaged two-dimensional equations can be derived and 2-D numerical simulations are normally sufficient (Müller et al., 2012; Fung & Chiang, 2016).

However, many discovered exoplanets have their orbital angular momentum vector misaligned with the stellar spin axis (see review by Winn & Fabrycky 2015). Especially for Hot Jupiters ($M_{p}>$0.3 $M_{J}$ and $P<10$ days), the misalignment spans the entire range from $0^{o}$ to $180^{o}$, revealed by measurements of the Rossiter-McLaughlin effect (Albrecht et al., 2012). The misalignment can be due to 1) processes through which the planet moves away from the disc midplane where it formed, such as planet-planet scattering (e.g. Chatterjee et al. 2008; Jurić & Tremaine 2008), Kozai-Lidov oscillations (e.g. Wu & Murray 2003; Petrovich 2015), and secular chaos (Wu & Lithwick, 2011), or 2) the primordial misalignment of the natal disc with respect to the stellar spin during the disc evolution, such as disc formation in a turbulent environment (Bate et al., 2010), magnetic torque from the star (Lai et al., 2011), or the gravitational torque from a misaligned companion star (Batygin, 2012; Martin et al., 2016). In either scenario, as long as the planet is misaligned with the disc, understanding how the misaligned planet can interact with the disc is crucial for studying the planet’s orbital evolution afterwards.

For an inclined low-mass planet which does not induce gaps in the disc, its interaction with the disc has been studied analytically in both the small inclination limit using the linear perturbation theory (Tanaka & Ward, 2004) and the large inclination limit with the dynamical friction theory (Rein, 2012). Numerically, Cresswell et al. (2007) and Bitsch & Kley (2011) have confirmed that low mass planets with small inclinations are undergoing exponential inclination decay, consistent with the linear theory. However, their measurements for moderately inclined planets are not consistent with the dynamical friction theory. This disagreement is likely due to the narrow range of the planet inclination explored by these simulations, as pointed out by Arzamasskiy et al. (2018). Arzamasskiy et al. (2018) have measured the planet’s migration and inclination damping rate for planets with inclination from 0o to 180o, and found a good agreement between the simulation results and the analytical theory.

While the inclined low-mass planets are relatively well studied, the situation is less clear for high mass planets which can induce gaps. Marzari & Nelson (2009) have found that both the eccentricity and inclination of giant planets damp very quickly, on the timescale of hundreds of orbits. On the other hand, Xiang-Gruess & Papaloizou (2013) and Bitsch et al. (2013) have found that the damping rates are reduced when a gap is opened in the disc. However, it is not quantitatively understood how migration and inclination damping rates depend on the gap depth and how the gap depth depends on the planet mass. Furthermore, there could be long-term secular interactions between the planet and the disc (Lubow & Ogilvie, 2001).

While most studies focus on the planet’s orbital evolution in discs, the planet can also affect the disc structure which has barely been studied. Xiang-Gruess & Papaloizou (2013) have found that the presence of a massive planet can cause a warp in the protoplanetary disc. However, the warp in the disc due to even a 6 $M_{J}$ planet only has $\lesssim$ 20o misalignment, which is not enough to explain some of the observational signatures of disc misalignment.

Observationally, some protoplanetary discs show dark spikes on their near-IR scattered light images. The dark spikes are best explained as the shadows cast by a misaligned inner disc blocking the stellar light (Marino et al., 2015; Stolker et al., 2016; Benisty et al., 2017; Debes et al., 2017; Long et al., 2017; Min et al., 2017; Casassus et al., 2018; Benisty et al., 2018). The inferred relative inclination between the inner and outer discs can be quite large: $\sim 70^{o}$ for HD 142527 (Marino et al., 2015), $\sim 72^{o}$ for HD 100453 (Benisty et al., 2017), $\sim 30^{o}$ for DoAr 44 (Casassus et al., 2018), and $\sim 30^{o}$ for HD 143006 (Benisty et al., 2018). Furthermore, some dark spikes vary with time, indicating the change of the inner disc on short timescales (Stolker et al., 2016; Debes et al., 2017). At much longer wavelengths, ALMA detect twisted flow patterns in the disc, directly suggesting a warped inner disc (Rosenfeld et al., 2014; Pineda et al., 2014; Casassus et al., 2015; Brinch et al., 2016; Loomis et al., 2017; Walsh et al., 2017). Finally, the optical/near-IR light curves of some YSOs show dimming events, indicating blocking the stellar light by a warped inner disc (Alencar et al., 2010; Bouvier et al., 2013; Lodato & Facchini, 2013; Cody et al., 2014; Facchini et al., 2016; Bodman et al., 2017; Schneider et al., 2018).

To explain these observations, Facchini et al. (2014); Juhász & Facchini (2017); Facchini et al. (2018); Price et al. (2018) propose that compact binaries can warp and break the circumbinary discs. These circumbinary disc simulations are imported into the radiative transfer code to generate the synthetic near-IR scattered light images and (sub-)mm molecular line channel maps.

In this paper, extending our previous study on planet-disc interactions for low mass planets on inclined orbits (Arzamasskiy et al., 2018), we use both the analytical theory and 3-D hydrodynamical simulations to study how a massive misaligned planet can interact with the disc. Different from previous simulations, we find that, when the planet mass is large enough and the induced gap is deep enough, the disc inside the planet’s orbit breaks from the outer disc. Then we use Monte-Carlo radiative transfer simulations to calculate the observational signatures of such misaligned discs. The analytical theory for gap depth/width, planet migration, and disc breaking is presented in §2. The numerical confirmation is presented in §3 and §4. The observational signatures are shown in §5. After a short discussion in §6, the paper is concluded in §7.

2 Theoretical Framework

The gap-opening process by planets has been studied extensively due to its importance in reducing the planet migration rate (Lin & Papaloizou, 1986) and explaining recent observations of discs with gaps and cavities (e.g. Espaillat et al. 2014). The quantitative relationships between the gap shape, migration rate, and the planet and disc properties have been worked out in great detail (Duffell et al., 2014; Fung et al., 2014; Kanagawa et al., 2015, 2016; Ginzburg & Sari, 2018; Dürmann & Kley, 2015). However, these studies only focus on planets that are coplanar with the disc. Numerical simulations by Xiang-Gruess & Papaloizou (2013) and Bitsch et al. (2013) have indicated that misaligned massive planets can behave very differently from coplanar planets. For example, a gap can be induced by a misaligned planet and the gap is shallower than the gap in the coplanar case. After a gap is induced, the planet’s inclination damping rate is reduced and a warp can develop throughout the disc.

In this work, we try to establish a theoretical framework for understanding the interaction between a misaligned high mass planet and the protoplanetary disc. We want to quantitatively answer following questions:

1. How are the depth and width of the gap determined by the planet mass, the planet inclination, and the disc properties ?

2. How does the gap opening process affect the planet migration rate and inclination damping rate?

3. *How does the gap opening process affect the disc evolution? How deep does the induced gap have to be for the inner disc to break from the outer disc and undergo differential precession? *

2.1 Gap Depth and Width

We first review the relationships between the planet mass and the gap profile (depth and width) for coplanar planets and then generalize them for inclined planets.

A coplanar planet excites density waves at Lindblad resonances. On each side of the planet, the total amount of angular momentum that is carried away by these density waves is derived by Goldreich & Tremaine (1979):

[TABLE]

where $\Omega_{p}$ is the disc orbital frequency at the planet position ($R_{p}$), $H$ is the disc scale height, $\Sigma_{0}$ is the unperturbed gas surface density at the planet position, and $q=M_{p}/M_{*}$ is the mass ratio between the companion and the central star. Later, we also use $h$ to represent $H/R$. Most of these waves are excited at radii with a distance of $H$ away from the planet, so-called torque cut-off (Goldreich & Tremaine, 1980). For a planet within a gap, the angular momentum of the excited waves can be calculated similarly but $\Sigma_{0}$ in Equation 1 is replaced with the surface density within the gap $\Sigma_{gap}$ (Fung et al., 2014):

[TABLE]

A more accurate derivation requires calculating the torque self-consistently considering the gap profile (Ginzburg & Sari, 2018). The induced density waves will propagate in the disc carrying a constant angular momentum flux. When the density waves steepen into shocks (Goodman & Rafikov, 2001), this angular momentum will be deposited to the disc, leading to gap opening.

In a viscous disc, a steady gap profile can be maintained when the angular momentum deposition into the disc balances the viscous stress (Fung et al., 2014):

[TABLE]

where $\nu=\alpha H^{2}\Omega$ is the kinematic viscosity. Equating equations 2 and 3, Fung et al. (2014) have shown that

[TABLE]

With a more detailed torque and angular momentum flux calculation, Kanagawa et al. (2015) has refined the gap-depth relationship:

[TABLE]

where

[TABLE]

which reproduces gaps from numerical simulations extremely well.

Regarding the gap width, Kanagawa et al. (2016) have found an empirical relationship

[TABLE]

where

[TABLE]

where $\Delta_{gap}$ is the radial extend of the gap at $\Sigma=\Sigma_{0}/2$.

For an inclined planet, we can simply extend these relationships using an averaged planet potential. Assuming the inclined planet’s orbital momentum vector is at an angle of $i_{p}$ from the disc’s angular momentum vector, we first calculate the projection of the planet’s position onto the disc. Then, we calculate the potential experienced by the disc element that is at a radial distance of $H$ away from this projected position. Finally, we average the derived potential over one planetary orbital period to get an effective planet potential onto the disc

[TABLE]

where $(R_{p}{\rm sin}\phi{\rm sin}i_{p})^{2}+H^{2})^{1/2}$ is the distance from the planet to the position in the disc that is one scale height away from the planet’s projection onto the disc, and $\phi$ is the angle between the vector to the planet and the line of nodes between the disc plane and the planet’s orbital plane. When $i_{p}$=0, this equation is $\Phi=GM_{p}/H$. Thus, we can generalize Equations 6 and 8 using the effective potential and replacing $q$ in these equations with

[TABLE]

The integration can be simplified as

[TABLE]

where $n$ is $(R_{p}{\rm sin}(i_{p})/H)^{2}$ and $\mathcal{K}(m)$ is the complete elliptic integral of the first kind ($m$=$k^{2}$ and $k$ is the elliptic modulus).

Thus, the gap depth and width for the inclined planet are

[TABLE]

with the new

[TABLE]

In the derivation above, we have ignored the effect of dynamical friction. Dynamical friction is important for planet migration (§2.2) but may not be important for gap opening (at least for moderately inclined planets) due to two reasons: 1) dynamical friction is a local effect which only occurs at the place where the planet moves through the disc, 2) the angular momentum exchange between the planet and the Mach cone left behind the planet is small compared with the one-side wave torque due to the Lindblad resonances. On the other hand, the detailed calculation on the gap opening due to dynamical friction still needs to be carried out in future and may be important for highly inclined planets.

2.2 Migration and Inclination Damping For Misaligned Planets

Arzamasskiy et al. (2018) found that the migration and the inclination damping rates for inclined low-mass planets can be described by the linear theory when the planet is mildly inclined (Tanaka et al., 2002; Tanaka & Ward, 2004) or by the dynamical friction theory when the planet is more inclined (Rein, 2012). The resulting migration rate can be described as the minimum between these two rates:

[TABLE]

where

[TABLE]

is the migration timescale, and $\alpha_{s}$ is the radial slope of the surface density profile ($\Sigma\propto R^{-\alpha_{s}}$).

The inclination damping rate can also be expressed as the minimum between these two limits:

[TABLE]

where

[TABLE]

is the inclination damping timescale. So the inclination damping timescale is shorter than the migration timescale by a factor of $h^{-2}$. Both the migration and inclination damping timescales are shorter for a more massive planet in a more massive and thinner disc.

On the other hand, for a massive planet which can induce a gap in the disc, the migration and inclination damping timescales should become longer since $\Sigma_{0}$ should be replaced by $\Sigma_{gap}$ considering that both wave launching and dynamical friction are proportional to the local disc surface density around the planet. Thus, we need to plug the new gap depth (Equations 12 and 13) into Equations 15 and 17 to get the migration and inclination damping rates. The migration rate is now

[TABLE]

with

[TABLE]

And the inclination damping rate is now

[TABLE]

with

[TABLE]

Here we add a free parameter $qs$ to represent the sensitivity of the planet migration and damping rates with respect to the gap depth. If the rates are proportional to $\Sigma_{gap}$, $qs$ is 1. However, since a misaligned planet can also interact with disc material outside the center of the gap, we expect $qs<1$. Later, by fitting results from numerical simulations, we empirically derive $qs\sim 2/3$. On the other hand, this fudge factor $qs$ is highly uncertain and may not apply to a large disc/planet parameter range.

2.3 Inner Disc Precession And The Disc Breaking Condition

As will be shown by numerical simulations in §4, when a deep gap is induced, the inner disc can lose connection with the outer disc so that the disc can break. The inner and outer discs will then precess at different rates driven by the misaligned planet, leading to a large disc misalignment.

To understand the disc breaking, we use the equations for small amplitude warping disturbances in a nearly inviscid disc (Lubow & Ogilvie, 2000) :

[TABLE]

where ${\bf l}(R,t)$ is the tilt vector at radius $R$. It is defined as a unit vector parallel to the angular momentum vector of each annulus in the disc. Equation 20 expresses the horizontal components of the angular momentum conservation. 2$\pi{\bf G(R,t)}$ is the horizontal internal torque in the disc, and ${\bf T(R,t)}$ is the horizontal torque density. Equation 21 determines the evolution of the internal torque G. The internal torque is affected by the shearing epicyclic motions, radial pressure gradients, and viscous decay described by the $\alpha$ viscosity parameter. We can integrate Equation 20 with $R$ to derive the precession rate of the inner disc:

[TABLE]

For a planet in the disc, its tidal torque consists of m=0 and m=2 components (Bate et al., 2002). Since the m=2 component only causes oscillation of the precession rate, we focus on the m=0 component. The m=0 component is from an axisymmetric potential, and can be written as

[TABLE]

where $i_{p}$ is the angle between the disc’s rotation axis and the planet’s angular momentum vector. $\phi$ is the azimuthal angle in the disc starting from the line of nodes between the disc plane and the planet’s orbital plane.

If the disc breaking occurs at $R_{p}$ so that the inner disc precesses on its own, the precession rate for the inner disc can be calculated using the torque from Equation 23. Assuming that the whole inner disc precesses as a rigid body at a rate of $\omega_{rb}$, $\partial{\bf l}/\partial t$ in Equation 22 is then sin$i\times\omega_{rb}$. The internal torque term (the $G$ term) in Equation 22 can also be dropped due to the disc breaking. Plugging Equation 23 into Equation 22, we can derive

[TABLE]

If the disc surface density follows $\Sigma\propto R^{-\alpha_{s}}$ from $R=0$ to $R_{d}$ with $R_{d}$ being the size of the inner disc, we have

[TABLE]

where $\Sigma_{d}$ and $\Omega_{d}$ are the surface density and orbital frequency at $R_{d}$. This equation has been derived by Terquem (1998), Bate et al. (2002), and Lai (2014). The rigid body approximation is a good approximation as long as the disc can communicate through wave (Larwood et al., 1996) or viscous torques (Wijers & Pringle, 1999). With $\alpha_{s}$=1 in our simulations and the relationship

[TABLE]

we can derive

[TABLE]

In our simulations with $q\ll 1$ where the gap inner edge is close to the planet position ($R_{d}\sim R_{p}$), we have

[TABLE]

and the precession timescale is $T_{pre}=1/\omega_{rb}$. When the inner disc is small and $R_{d}$ is significantly smaller than $R_{p}$, the precession rate can be a lot lower than the estimate from Equation 28, since $\omega_{rb}$ in Equation 27 has $(R_{d}/R_{p})^{3/2}$ dependence.

Now, let’s derive under what condition the disc breaking can occur at $R_{p}$. We follow Lubow & Ogilvie (2000) and adopt the complex representation for $\bf{l}$ and $\bf{G}$, which are $W=l_{x}+il_{y}$ and $G=G_{x}+iG_{y}$. Equations 20 and 21 become

[TABLE]

where $\Omega$ and $\kappa$ are defined as Equation (25) and (26) in Lubow & Ogilvie (2000), and $\Omega_{z}$ is

[TABLE]

where $b_{3/2}^{(1)}$ is the Laplace coefficient as in Lubow & Ogilvie (2000). We note that both $\frac{\Omega^{2}-\Omega_{z}^{2}}{2\Omega}$ and $\frac{\Omega^{2}-\kappa^{2}}{2\Omega}$ are on the order of $\omega_{rb}$ when the axisymmetric contribution of the companion potential has been taken into account.

When the disc breaking does not occur, the whole disc including both the inner and outer discs precesses as a whole and the resulting precession rate $\omega_{p}$ is much smaller than $\omega_{rb}$ in Equation 28. Thus, we drop the time dependent terms and integrate Equation 29 from [math] to $R$:

[TABLE]

If we plug this $G$ into Equation 30 and also drop the time dependent terms, we have

[TABLE]

If we replace the $\Omega^{2}$ terms with $\omega_{rb}$ and use some averaged surface density within $R$ (labled as $\overline{\Sigma}$) to simplify the integral from 0 to R, we have

[TABLE]

In order for the inner and outer discs to precess together as one ridge body, $|W|$ needs to be a constant or $|\partial W/\partial R\times R|$ needs to be small throughout the disc. In other words, $|\partial W/\partial R\times R|<W$ for all $R$. For orders of magnitude argument, we thus have the rigid body precession condition as

[TABLE]

This equation is most difficult to be satisfied at the radius where $P$ is small. In our planet-disc interaction problem, $P$ is the smallest at the deepest part of the gap (at $\sim R_{p}$) where $P=c_{s}^{2}\Sigma_{gap}$. The rigid rotation condition is then

[TABLE]

Our simulations later will show that the disc cannot break in a highly viscous disc since the gap is shallow. Thus we will limit ourselves to inviscid discs where the breaking condition is then

[TABLE]

or

[TABLE]

Since $\overline{\Sigma}(R_{p})$ represents the averaged disc surface density at the inner disc within $R_{p}$, $\overline{\Sigma}(R_{p})$ is $\sim\Sigma_{0}(R_{p})$ so that $\Sigma_{gap}/\overline{\Sigma}$ can be considered as the gap depth.

If the disc is smooth without a gap so that $\Sigma_{gap}\sim\Sigma_{0}$, Equation 37 becomes the traditional disc breaking condition that the sound crossing time is longer than the precession timescale (Larwood et al., 1996). When a deep gap is formed, Equation 37 can also be written as

[TABLE]

which suggests that the breaking condition is that the sound crossing time is longer than the product of the precession timescale and the square root of the gap depth. The surface density factor reflects that less surface density communicates different parts of the disc less efficiently.

In this section, we derive the gap depth and width for misaligned planets (Equations 12 and 13), the migration and inclination damping rates for planets in gaps (Equations 14 to 19), and the disc breaking condition (Equation 37 and 38). In the next two sections, we will carry out numerical simulations to test these equations.

3 Methods

We solve the compressible Navier-Stokes equations using Athena++ (Stone et al. 2018, in preparation). Athena++ is a newly developed grid based magnetohydrodynamic code using a higher-order Godunov scheme for MHD and the constrained transport (CT) for magnetic fields. Compared with its predecessor Athena (Gardiner & Stone, 2005, 2008; Stone et al., 2008), Athena++ is highly optimized for speed and uses a flexible grid structure that enables mesh refinement, allowing global numerical simulations spanning a large spatial range. Furthermore, the geometric source terms in curvilinear coordinates (e.g. in cylindrical and spherical-polar coordinates) are carefully implemented so that angular momentum is conserved to machine precision, an important feature for studying the disc precession.

We adopt the spherical-polar coordinate system ($r$, $\theta$, $\phi$) in our simulations, with $r$ from 0.2 to 10, $\theta$ from 0.2 to $\pi$-0.2, and $\phi$ from 0 to $2\pi$. In the radial direction, we have 160 grid points uniformly in log($r$). In the $\theta$ and $\phi$ directions, we have 132 and 256 uniform grid points. One level of mesh refinement is used towards the disc midplane, allowing us to accurately simulate the disc evolution but also not to be limited by the small timestep at the disc atmosphere. The grid structure is shown in Figure 1. The outflow boundary condition111This outflow condition is different from the default outflow condition in the code. We follow the standard practice for the outflow boundary condition (Stone & Norman, 1992): if $v_{r}$ in the ghost zones are pointing to the active zones, they are set to be zero to limit the inflow., reflecting boundary condition, and periodic boundary condition are adopted in the $r$, $\theta$, and $\phi$ direction respectively.

Although we adopt spherical-polar coordinates for the simulations, we use cylindrical coordinates to set up the initial condition. In this paper, we use ($R$, $\phi$, $z$) to denote positions in cylindrical coordinates and ($r$, $\theta$, $\phi$) to denote positions in spherical-polar coordinates. In both coordinate systems, $\phi$ always represents the azimuthal direction (the direction of disc rotation).

The initial density profile at the disc midplane is

[TABLE]

where $p$ is set to be -2.25. $R_{0}$ equals the planet’s semimajor axis ($R_{p}$) and $R_{0}=R_{p}=1$ in code units. The temperature is assumed to be constant on cylinders:

[TABLE]

where $q$ is set to be -0.5. With our choice of $p$ and $q$, the disc surface density in our simulations follows $\Sigma\propto R^{-1}$. For our main set of simulations, $T(R_{0})$ is chosen so that the disc scale height $H=c_{s}/\Omega$ at $R_{0}$ is 0.1, where $c_{s}=\sqrt{p/\rho}$ is the isothermal sound speed. This disc scale height at $R_{0}$ is resolved by 10 grid points with out default resolution.

Hydrostatic equilibrium in the $R-z$ plane requires that (e.g. Nelson et al. 2013)

[TABLE]

and

[TABLE]

where $v_{K}=\sqrt{GM_{*}/R}$. For thermodynamics, we adopt the adiabatic equation of state with $\gamma$=1.4, but with a short cooling time as in Zhu et al. (2015). We have chosen the cooling time as 0.01 of the local orbital time, which is motivated by the realistic cooling time for disc regions at $\sim$100 AU (Zhu et al., 2015). Discs with such short cooling time behave very similarly to the disc with the locally isothermal equation of state.

To simulate the interaction between the misaligned planet and the disc, we can either (1) set the disc midplane at the midplane of the grid (the $\theta=\pi/2$ plane) and put the planet on an inclined orbit, as done in Arzamasskiy et al. (2018), or (2) we can set the planet’s orbital plane at the midplane of the grid and tilt the disc with respect to the grid midplane. The setup (1) ensures that the disc’s Keplerian motion is along one of the grid directions ($\phi$ direction) so that it minimizes the grid noise for simulating a Keplerian disc. Thus, we adopted this setup in Arzamasskiy et al. (2018), where the low mass planet does not affect the disc evolution much. However, in current work with a high mass planet, the inner disc can undergo significant nodal precession around the total angular momentum vector of the whole system. Since the disc’s mass is assumed to be negligible compared with the planet’s mass in our simulations, the total angular momentum vector equals the planet’s orbital angular momentum vector. Then, with the setup (1), the disc can move away from the grid midplane during the nodal precession. For example, if the planet is 45 degrees inclined to the disc plane, the disc will move away from the midplane to the polar region of the grid after it precesses for 180o. Although the grid midplane has small numerical diffusion, the polar region has the lowest resolution and largest numerical diffusion. Considering this disadvantage, in this work we adopt the setup (2) where the planet’s orbital plane is at the grid midplane and the disc is tilted with respect to the grid midplane (Figure 1). In this setup, the total angular momentum vector is along the polar direction, and the disc’s angular momentum vector keeps the same angle from the gird polar direction during the disc’s precession so that the tilt between the disc and the grid midplane remains a constant and the numerical error won’t change dramatically with time.

Suppose that the disc is tilted by an angle of $i_{d}$, the coordinates of any point in the coordinate system that has the tilted disc as the midplane ($r^{\prime}$, $\theta^{\prime}$, $\phi^{\prime}$) are related to its coordinates in the numerical grid ($r$, $\theta$, $\phi$):

[TABLE]

where $i_{d}$ is the disc inclination. When $i_{d}>0$, the disc’s angular momentum vector is pointing to the positive $y$ direction in the numerical grid. The Cartesian coordinates ($x$, $y$, $z$) correspond to the spherical coordinates ($r$, $\theta$, $\phi$), while ($x^{\prime}$, $y^{\prime}$, $z^{\prime}$) correspond to ($r^{\prime}$, $\theta^{\prime}$, $\phi^{\prime}$). With $\theta^{\prime}$ calculated for each grid point at ($r$, $\theta$, $\phi$), we can calculate $R^{\prime}$ and $z^{\prime}$ with $R^{\prime}=r\sin\theta^{\prime}$ and $z^{\prime}=r\cos\theta^{\prime}$, and then use Equations 41, 42, and 43 to calculate the temperature, density, and velocities for every grid point ($R$ and $z$ in these equations are now $R^{\prime}$ and $z^{\prime}$).

Our main suite of simulations include 3$\times$3$\times$2=18 simulations with 3 different disc inclinations ($i_{d}$=-0.17, -0.34, -0.68 222So the disc’s angular momentum vector is pointing to the negative $y$ direction in the initial direction, or equivalent to the planet inclination $i_{p}=$0.17, 0.34, 0.68 or 10o, 19o, 39o ), 3 viscosity parameters ($\alpha$=10*-3*, 10*-2*, and 10*-1*), and 2 planet masses ($q=0.003$ and 0.01, or 3 $M_{J}$ and 10 $M_{J}$). The planet is fixed to be on a circular orbit and its potential has a smoothing length of 0.6 $H$. As shown by Rein (2012), dynamical friction has a weak dependence on the smoothing length. We label each simulation in the following manner: P10I10AM1 means a 10 Jupiter mass planet (or $q=0.01$) on a 10o inclined orbit in an $\alpha=10^{-1}$ disc (M1 means minus 1).

We also carry out two additional simulations to complement the main suite of simulations: 1) one inviscid simulation ($\alpha=0$) with $i_{p}$=19o and $M_{p}=10M_{J}$ (P10I19A0) to be compared with the corresponding viscous simulations, 2) a thin disc with $h$=0.05, $\alpha=10^{-3}$, $M_{p}$=10 $M_{J}$, and $i_{p}$=39o(THINP10I39AM3) to study disc breaking in a thin disc.

3.1 Numerical Convergence

To test the numerical convergence of our setup, we also carry out a suite of inviscid simulations using different resolutions and numerical reconstruction schemes. These simulations have the same setup as our fiducial cases with $i_{p}$=19o, except that the planet mass is $q=0.03$ and the radial extend is from $r_{min}=0.1$ to $r_{max}=1$. The larger planet mass makes the disc undergo significant precession during a shorter period of time. The resulting radial profiles of the disc surface density ($\Sigma$), tilt ($i$), and twist ($\gamma$) at 4.2 planetary orbits (4.2 $T_{p}$) are presented in Figure 2. Different resolutions and reconstruction methods (second-order piecewise linear method: PLM, or third-order piecewise parabolic method: PPM) are explored.

The tilt vector at each radius (${\bf l}(R,t)$ in Equation 20) characterizes how a disc is warped. The tilt angle ($i$) and twist angle ($\gamma$) at radius $R$ can be calculated using $\bf{l}$:

[TABLE]

Our definition of twist is to ensure that the initial twist angle is zero. We also define that the twist angle increases in the clockwise direction in the $x-y$ plane.

Figure 2 shows that the simulation with both the PLM reconstruction scheme and 5 grids per scale height at $R_{p}$ has very different $i$ and $\gamma$ profiles compared with simulations having a higher resolution or the PPM reconstruction scheme. Thus, we conclude that, if we use the PLM reconstruction scheme, we need at least 10 grids per scale height. If we use the PPM reconstruction scheme, 5 grids per scale height is sufficient. On the other hand, we notice that the simulation having both 5 grids per scale height and the PPM scheme still has some minor differences from the higher resolution runs. Thus, we decide to use 10 grids per scale height with the PLM scheme as our fiducial setup. On the other hand, with our fiducial setup, the resolution is only 5 grids per scale height for the thin disc simulation THINP10I39AM3. Thus, we use the PPM scheme for this thin disc simulation.

4 Results

We run the simulations for 200-300 planetary orbits. The poloidal cut for the disc density structure of P10I19A0 is shown in Figure 1. The planet, together with the circumplanetary region, can be seen on the right at the grid midplane. Initially, both the inner and outer discs are tilted by 19o. However, the massive planet has carved out a very deep gap so that the inner and outer discs break from each other and precess at different precession rates. After 180 planetary orbits, the inner disc precesses for $\sim$180o around the z-axis and it is now facing the other direction. Thus, the angle between the angular momentum vectors of the inner and outer discs is now 19${}^{o}\times$ 2=38o. The isodensity surface for this snapshot is shown in Figure 3. We can clearly see that the inner and outer discs are misaligned. Spirals are apparent at both the inner and outer discs. Another interesting phenomenon in both Figure 1 and 3 is that the circumplanetary disc (CPD) around the planet is aligned with the outer disc instead of the inner disc. This may be because the CPD forms by accreting material from the outer disc and this incoming material carries the angular momentum of the outer disc.

4.1 Disc Structure

At 200 planetary orbits, the radial profiles of the disc surface density, tilt, and twist are shown in Figures 4 and 5 with a 3 $M_{J}$ and 10 $M_{J}$ planet in the disc respectively. There are several noticeable trends from these figures:

1. The larger $\alpha$ the disc has, the shallower gap the planet induces. This is similar to gap opening by coplanar planets (Fung et al., 2014; Kanagawa et al., 2015). While the planet tries to open a gap by depositing angular momentum into the disc, the disc tries to close the gap due to viscous diffusion. Thus, a higher viscosity leads to a shallower gap. One exception to notice in Figure 4 is that the inviscid case P10I19A0 has a shallower gap than the viscous case P10I19AM3. This is observed in coplanar simulations too (e.g. Zhang et al. 2018). Two effects play roles here: a) a very low viscosity in disks can trigger instabilities (e.g. Rossby Wave Instability, Lovelace et al. 2009) at the gap edge, which leads to turbulence that will close the gap; b) the gap edge can also become eccentric when the viscosity is low (Lubow, 1991a, b; Kley & Dirksen, 2006; Teyssandier & Ogilvie, 2017), in which case azimuthally averaging over an elliptical gap can smear out the density profile of the gap.

2. The more inclined the planet is, the shallower gap the planet induces. This is consistent with our derivation (Equations 10 and 13) that the inclined planet has less gravitational interaction with the disc. We will compare the gap depth with Equation 13 in more detail in §4.3.

3. The tilt ($i$) of both the inner and outer discs does not changed significantly at 200 orbits. When $\alpha$ is large (e.g. $\alpha$=0.1), the tilt seems to decrease slightly towards the inner disc. When $\alpha$ is small (e.g. $\alpha$=$10^{-3}$), the tilt seems to increase slightly towards the inner disc if the gap is not deep (e.g. Figure 5).

4. The deeper the gap is, the larger the difference between the twist of the inner and outer discs is. For the three cases which have the deepest gaps ( P10I10AM3, P10I19AM3, P10I19A0), the twist difference between the inner and outer discs reaches more than 200o. For discs with shallow gaps, the twist difference between the inner and outer discs is limited to a small value.

Figures 6 and 7 show the time evolution of the tilt and twist of the inner (solid curves) and outer (dashed curves) discs. The inner and outer discs are defined as the disc region smaller than $R=0.7$ and larger than $R=1.5$ respectively. Their tilt and twist angles are calculated using the angular momentum vector that has been integrated over the inner and outer discs. For the three cases that have deep gaps ( P10I10AM3, P10I19AM3, P10I19A0), we can see that the inner and outer discs precess at very different rates. The twist ($\gamma$) of the inner disc increases much faster than that of the outer disc, indicating that the two discs break and they are not dynamically connected. The measured precession rates ($\partial\gamma/\partial t$) for the inner discs of P10I10AM3 and P10I19AM3 are $\Omega_{p}/340$ and $\Omega_{p}/300$. Considering that the analytical prediction from Equation 28 is $\Omega_{p}/276$ with i=15o and q=0.01, the simulation results agree well with the theory.

On the other hand, for the shallow gap cases, e.g. P10I19AM2, P3I10AM3, the twist of the inner disc increases fast initially. But it slows down at later times and precesses at the same rate as the outer disc. The inner and outer discs thus have a constant relative twist angle. It seems that the disc tries to break but manages to maintain a twist balance between the inner and outer discs. We will discuss disc breaking more quantitatively in §4.3.

4.2 Migration and Inclination Damping Rates

Although the planet is kept on a fixed circular orbit, we can calculate the planet migration rate and the inclination damping rate using the gravitational force exerted on the planet by the disc. Following Burns (1976), the gravitational force experienced by the planet can be decomposed into

[TABLE]

where $\hat{e_{R}}$, $\hat{e_{T}}$, and $\hat{e_{N}}$ are unit vectors along the radial direction, planet’s velocity direction, and the normal direction to the planet’s orbital plan (the direction of $\hat{e_{R}}\times\hat{e_{T}}$). These components of the force lead to the planet’s orbital evolution,

[TABLE]

where we assume that the eccentricity is zero, $i_{p}$ is the planet’s inclination angle with respect to the outer disc, and $\Omega_{p}t$ represents the planet’s azimuthal angle with respect to the x-axis on the planet’s orbital plane. The quantity $\gamma_{out}(t)$ is the outer disc’s tilt angle at time $t$, and thus it is the angle from the x-axis to the line of nodes between the planet’s plane and the outer disc plane. $(\Omega_{p}t-\gamma_{out}(t))$ is thus the angle from the line of nodes to the planet’s position vector on the planet’s orbital plane. Note that both $r_{p}$ and $i_{p}$ are osculating elements of the orbit. Thus, they vary during one orbital time. To get the long-term orbital evolution of the planet, we integrate $dr_{p}/dt$ and $di_{p}/dt$ from 160 to 200 planetary orbits to calculate the migration and inclination damping rates. Since the outer disc precesses slowly as shown in Figures 6 and 7, we set $\gamma_{out}=0$ for these calculations. We have verified that choosing a slightly different $\gamma_{out}$ based on Figures 6 and 7 has negligible effect on the results.

Figures 8 and 9 show the planets’ migration and inclination damping rates in our simulations. Different colors represent cases with different disc viscosities. The upper panels show the measured rates in code units, while the lower panels show the normalized rates. The normalized rates are defined as

[TABLE]

where $t_{mig}$ and $t_{inc}$ are given in Equations 15 and 17. $\Sigma_{gap}$ is defined as $(\Sigma(R_{p}-H/2)+\Sigma(R_{p}+H/2))/2$. We did not use the surface density at $R_{p}$ as the gap surface density since both the horseshoe material and the circumplanetary region lead to a density spike at $R_{p}$. In Equations 50 and 51, we normalize the rates using the gap surface density. This is because we expect the rates to be roughly proportional to the gap surface density since the dynamical friction process, which determines the orbital evolution of a moderately inclined planet, is the gravitational interaction between the planet and the local disc where the planet is traveling through.

If a gap is not induced by a planet, Equations 50 and 51 are reduced to Equations 14 and 16 except that $t_{mig}$ and $t_{inc}$ are now moved to the left side. These normalized rates without gaps are plotted as the dotted curves in the bottom panels of Figures 8 and 9.

Overall, we can see that the migration and inclination damping rates decrease with a higher planet inclination and a deeper disc gap (lower $\alpha$). The normalized rates agree well with the analytical formulae for most cases. But the measured rates (especially the migration rates) are higher than the analytical formulae for 10 $M_{J}$ planets in $\alpha=10^{-3}$ discs where deep gaps are induced. This indicates that simply scaling the rates with the smallest density in the gap is not adequate. When the planet travels across the gap, its gravity not only interact with the deepest part of the gap but also with the gap edge and even the full disc region which has much higher density. Thus, the relationships between the rates and the gap depth may not be linear.

On the other hand, we can use simulation results to seek empirical relationships between the gap depth and the planet’s orbital evolution. Figure 10 summarizes the planet migration and inclination damping rates with respect to the gap depth. The dashed lines are from the analytical formulae (Equations 18 and 19) using $qs=2/3$ and parameters in the simulations ($\Sigma_{0}$ is set to be 0.07 as inferred from Figures 4 and 5).

If the rates are proportional to $\Sigma_{gap}$, the fudge factor $qs$ in Equations 18 and 19 is 1. Based on the argument above that a misaligned planet can also interact with the disc material beyond the center of the gap, we expect $qs<1$. We found that $qs$=2/3 gives a good fit to most data points (Figure 10). We want to emphasize that $qs$ is the fudge factor which we use to represent complicated physical processes. In reality, $qs$ could vary with both disc and planet parameters, and it can be different between migration rates and inclination damping rates. The two noticeable outliers are the high migration rates for planets with very deep gaps. Thus, something else besides the dynamical friction is pushing the planet inwards. We suspect that the high migration rates are related to the precession of the inner discs in these two cases. An analytical theory which takes into account both the non-uniform gap density profile and the precession of the inner disc is needed in future.

4.3 Gap Opening and Disc Breaking

After studying the planet’s orbital evolution, we study the disc that is under the influence of the planet. We use simulation results to test our analytical estimate on the gap depth and the disc breaking condition. Analytical theory (Equation 12) suggests that the gap depth and width for a misaligned planet are similar to those for a coplanar planet but with a modified planet mass (Equation 11 and 13). In Figure 11, we plot both the measured gap depth from simulations and the prediction from the analytical formula (Equation 13). Good agreements are found for both $q=0.003$ and $q=0.01$ cases, as long as the disc is not undergoing disc breaking (P10I10AM3, P10I19AM3 have disc breaking.).

To test our disc breaking condition (Equation 37), we plot Equation 37 as dashed curves in Figure 11. Whenever the dashed curve is larger than the gap depth curve, the disc should undergo breaking. With a $q$=0.003 planet, the gap is not deep enough for disc breaking even with $\alpha=10^{-3}$. On the other hand, if the planet has $q$=0.01 and a low inclination, the induced gap in the $\alpha=10^{-3}$ disc is deep enough for the disc breaking. The fact that the discs in P10I10AM3 and P10I19AM3 undergo breaking is consistent with our breaking condition.

With the gap depth formula (Equation 12) and the disc breaking condition (Equation 37), we can calculate how massive the planet needs to be in order to break the disc and cause differential precession and large misalignment. For given disc parameters ($H/R$ and $\alpha$) and a given planet inclination ($i_{p}$), we gradually increase the planet mass until the gap depth (Equation 12) reaches the disc breaking condition (Equation 37). We define this planet mass as the disc-breaking planet mass. Planets with masses larger than this mass can cause disc breaking, while planets with less mass cannot break the disc and there will be a constant relative twist between the inner and outer discs. Figure 12 shows the disc-breaking planet mass for different values of $\alpha$ (different colors) and different values of $H/R$ (0.1: solid curves, 0.05: dotted curves, 0.03: dashed curves). We can see that only in extreme conditions, e.g. $H/R\lesssim$0.05 and $\alpha\lesssim 10^{-4}$, a Jupiter mass planet can break the disc. This figure also predicts that, in our thin disc simulation THINP10I39AM3, the thin disc with $H/R$=0.05 and $\alpha=10^{-3}$ (the yellow dotted curve) should break by the $39^{o}$ inclined $10\,M_{J}$ planet. We confirm this by checking the THINP10I39AM3 simulation directly.

5 Observational Signatures

The observational signatures of warped discs (Facchini et al., 2014; Juhász & Facchini, 2017) and broken discs (Facchini et al., 2018) have been studied recently. Although these calculations use circumbinary disc simulations, most of their results can also be applied to warped and broken discs induced by the planet, as shown in this section. We will also highlight some differences between the observational signatures of the broken circumbinary discs and the broken discs induced by the planet.

We calculate the near-IR scattered light images and (sub-)mm velocity channel maps for two simulations (P10I19A0 and THINP10I39AM3) where the gap is deep enough for the disc to break. The inner and outer discs thus have significant misalignment due to their different precession frequencies. We use RADMC-3D333RADMC-3D is an open code of radiative transfer calculations. The code is available online: http://www.ita.uni-heidelberg.de/ dullemond/software/radmc-3d/. for the 3-D radiative transfer calculations, following the similar procedures as in Arzamasskiy et al. (2018). A spherical mesh is adopted. Since one level of mesh refinement is used in our hydrodynamical simulations, we coarsen the fine mesh to generate a mesh with the uniform resolution. To scale the quantities in our simulations, the length unit is set to be 100 au, and the initial gas surface density at 100 au is set to be either 1 g cm*-2* or 100 g cm*-2*. The central star is assumed to be a Herbig Ae/Be star with $T_{eff}=6810$ K, $R_{*}$=1.4 $R_{\odot}$, and $M_{*}$=1.7 $M_{\odot}$. The initial dust to gas mass ratio is 1/100, and we have assumed that 10% of dust mass is in small grains which determine the disc temperature structure and scatter light at near-IR. The dust opacity is calculated using Mie theory for Magnesium-iron silicates (Dorschner et al., 1995) assuming a size of 0.1 $\mu$m.

The radiative transfer calculation with RADMC-3D includes 3 steps. First, the temperature structure of the disc is calculated. Second, using this temperature structure, the near-IR scattered light images are generated including full treatment of light polarization. The Stokes parameters (I, Q, U, V) are computed and the polarized intensity is $\sqrt{Q^{2}+U^{2}}$ (V=0 due to linear polarization from dust scattering). Third, assuming that the number density ratio between 13CO and H2 is 1.75$\times 10^{-6}$, the channel maps for the 13CO (3-2) line are generated.

The images presented below assume a geometry that we are viewing towards the disc along the polar direction in the simulation (so the simulation’s $z$-axis is pointing towards us). Thus, the orbit of the planet is on the plane of sky. The $x$-axis in our simulations points to the right direction in the images and the $y$-axis in the simulations points to the up direction in the images.

5.1 The Near-IR Scattered Light Image

The near-IR scattered light images for P10I19A0 with $\Sigma_{0}=1$ g cm*-2* are presented in Figure 13. The intensity has been multiplied by $R^{2}$ in the image. The colorbar is linear with the maximum chosen to highlight the disc features (the brightest color represents a normalized intensity that is a factor of 5 smaller than the maximum normalized intensity in the image). Since this is an inviscid simulation, we clearly see the gap edge vortex in the images. On the other hand, we also see a dark lane in the horizontal direction in the image. This dark lane is along the nodal lines between the inner and outer discs, and thus is due to the shadow from the inner inclined disc. Since the vortex orbits around the central star at the local Keplerian speed, we can see that the vortex is traveling across the dark lane in these snapshots. When the vortex and the dark lane overlap (T=179.4 $T_{p}$), the vortex looks like it splits into two vortices. After the vortex travels past the dark lane it looks like one vortex again. This is consistent with HD 142527 where the submm vortex seems to split into two parts (Casassus et al., 2015) at the dark lane of the near-IR scattered light image (Marino et al., 2015). Although the split of the vortex in HD 142527 is shown at (sub-)mm wavelengths where the dust thermal emission dominates, it can still be related to the shadow due to the decrease of the dust temperature at the shadow lane. The lower right panel of Figure 13 shows the dust temperature at the disc midplane. Clearly, the casting shadow lowers the disc temperature there, which is also seen in broken circumbinary disc simulations (Juhász & Facchini, 2017; Facchini et al., 2018).

The near-IR scattered light images for THINP10I39AM3 (Figure 14) are quite different from those for P10I19A0. First, due to the finite viscosity with $\alpha=10^{-3}$, vortices do not show up at the gap edge. This is consistent with simulations having coplanar planets in discs (Fu et al., 2014; Zhu & Stone, 2014; Hammer et al., 2017). Second, the shadows are much sharper. This is because, in THINP10I39AM3, the inner disc is both thinner and more misaligned. The larger misalignment is because the planet in THINP10I39AM3 has a higher inclination (39o compared with 19o in P10I19A0) so that the differential precession can lead to a much larger misalignment between the inner and outer discs ($\sim 60^{0}$ in Figure 14 compared with $\sim 35^{o}$ in Figure 13).

We also increase the disc surface density by a factor of 100 (lower panels of Figure 14) to explore how the disc surface density affects the disc shadow. Two effects can be seen in the figure. First, when the disc surface density is higher, the scattering surface at the inner disc is higher so that the inner disc casts a wider shadow on the outer disc. Second, when the disc surface density is higher, the scattering surface at the outer disc is also higher so that the shadow is cast onto a bowl shaped flaring surface instead of a flat surface close to the midplane. As shown in the bottom panels of Figure 14, the two dark lanes concave towards one direction tracing the scattering surface. Thus, we can use the dark lane positions to reconstruct the scattering surface in the disc.

Another interesting phenomenon we observe is that the shadow moves at a slower speed than the inner disc precession speed. By checking the simulation directly, we find that the inner disc’s tilt angle is 60o at $T=100\,T_{p}$ and 120o at $T=190\,T_{p}$. Although the inner disc’s tilt angle changes by 60o, the shadow only rotates $\sim 30^{o}$ by comparing the left and right panels in Figure 14. To understand this difference, Figure 15 shows the geometry of the inner and outer discs. We assign the outer disc midplane as the $x^{\prime}-y^{\prime}$ plane to be different from the planet’s orbital plane that is denoted as the $x-y$ plane. The total angular momentum vector ${\bf L}$ is at the $y^{\prime}-z^{\prime}$ plane and is pointing to the negative $y^{\prime}$ direction. The tilt angle between ${\bf L}$ and the $z^{\prime}$ axis is $i$. The inner disc’s angular momentum vector is ${\bf l_{in}}$ and is also at the $y^{\prime}-z^{\prime}$ plane initially. The angle between ${\bf L}$ and ${\bf l_{in}}$ is $\alpha$. After some time, ${\bf l_{in}}$ precesses around ${\bf L}$ for an angle of $\gamma$ to ${\bf l_{in}^{\prime}}$. We assume that the outer disc does not precess during this short period of time. The line of nodes, which is also the direction of the casting shadow, has an angle of $\theta$ from the $x^{\prime}$-axis in the $x^{\prime}-y^{\prime}$ plane. The relationship between all these angles is

[TABLE]

In our simulation setup, the inner and outer discs are coplanar initially. Thus ${\bf l_{in}}$ and the $z^{\prime}$-axis overlap and $\alpha=i$. So the equation is simplified to

[TABLE]

If we plug in $i=-39^{o}$ and $\gamma=-60^{o}$ at $T=100\,T_{p}$, $\theta$ is $66^{o}$. If we plug in $i=-39^{o}$ and $\gamma=-120^{o}$ at $T=190\,T_{p}$, $\theta$ is $37^{o}$. Considering that $\theta$ is the angle from the positive $x^{\prime}$-axis in the counterclockwise direction, these angles are consistent with the dark lanes in Figure 14. Thus, while the inner disc precesses for $60^{o}$, the shadow only precesses for $\sim 30^{o}$. We can use Equation 52 to explore several extreme scenarios. If the total angular momentum vector ${\bf L}$ is along the outer disc angular momentum vector ($z^{\prime}$-axis), we have $i=0$ and thus ${\rm tg}\theta={\rm tg}\gamma$. The shadow rotates at the same frequency as the inner disc precession. For the other extreme, if $\alpha=i=90^{o}$, $\theta$ is always 90o and the shadow never moves. These indicate that the rotation of shadow not only depends on the inner disc precession but also depends on the geometry of the system. The shadow may or may not move at the same frequency as the disc precession frequency.

Our derivation is similar to Min et al. (2017) where they have pointed out the position angle of the shadow is different from the position angle of the inner disc, except that we choose a different coordinate system that is aligned with the outer disc. This coordinate system is motivated by the disc’s nodal precession. It not only simplifies the calculation, but also makes it much easier to study the time variability of the shadow due to the disc precession.

We can relate the movement of the shadow with the inner disc’s precession by doing time derivative for Equation 52. We then have

[TABLE]

where $d\theta/dt$ is the rotational rate of the shadow while $d\gamma/dt$ is the precession rate of the inner disc. In our simulation setup with $\alpha=i$, we have

[TABLE]

Thus, the shadow can move faster or slower than the inner disc precession, depending on both $i$ and $\gamma$.

5.2 CO Channel Maps

The 13CO channel maps for THINP10I39AM3 are presented in Figure 16. Initially, both the inner and outer discs are oriented in a way that the right side in the image is blue-shifted while the left side is red-shifted. At $T=190$ $T_{p}$, the inner disc precesses around the $z$-axis (the $z$-axis is pointing to us along our line of sight) for 120o in the clockwise direction (or retrogradely compared with the planet’s orbital motion). Thus, the channel maps of the inner disc starts at $\sim$120o away from the positive $x$-axis in the clockwise direction, while the outer disc still starts at the positive $x$-axis. Besides the misalignment between the inner and outer discs, the shadow at $\theta\sim$30o and $\sim$-150o is also visible from the channel maps. This is due to the decrease of the disc temperature within the shadow. One noticeable difference between the $\Sigma_{0}=1$ g cm*-2* and 100 g cm*-2* cases is that, for the latter case, the 13CO emission surface is high in the atmosphere where it is hotter than the midplane so that we can see both the front and back sides of the disc (Rosenfeld et al., 2013).

The moment 0 and 1 maps that are derived from the channel maps in Figure 16 are shown in the middle column of Figure 17. We see less emission at the shadow lane in the moment 0 maps due to the lower temperature there. The different rotational directions between the inner and outer discs are apparent in the moment 1 maps. Even if we view this system at different angles (the other columns in Figure 17), we can still clearly see the different disc orientations between the inner and outer discs. In the -45o maps, we are viewing the outer disc almost edge on, while viewing the inner disc almost face on. In the 45o maps, we are viewing the inner disc edge on while viewing the outer disc face on. If the molecular tracer is optically thin, we can clearly see a wide gap between the inner and outer discs, which is different from disc breaking in circumbinary discs. On the other hand, when the disc surface density is high (the bottom panel block in Figure 17) or the line is more optically thick (e.g. 12CO lines), we may not see the gap between the inner and outer discs. But we can still see that the velocity structure is twisted in the moment 1 maps. Overall, unlike the near-IR scattered light images, the CO (sub-)mm moment maps provide direct information on the orientation of both the inner and outer discs, which can be quite helpful for probing the potential companion.

6 Discussion

6.1 Various Timescales

If the planet is massive enough to break the disc, the inner disc will precess at the precession timescale (Equation 28). Both Figures 12 and 18 show the relationship between the disc precession timescale, the planet inclination, and the planet mass. In Figure 12, the $y$-axis is the planet mass while the color scale is the precession timescale. In Figure 18, the $y$-axis is the precession timescale and the color scale is the planet mass. Combing these two figures, we can easily constrain the disc precession timescale if the disc is known to be broken. For example, for a broken disc with $\alpha=10^{-3}$ and $H/R=0.03$, Figure 12 tells us that $q$ has to be larger than $10^{-3}$ to break the disc. Then, Figure 18 or Equation 28 informs us that with $q>10^{-3}$ the disc precession timescale is shorter than 3000/$\Omega_{p}$.

Although the disc-breaking planet can cause the inner disc to precess (§4.3), the planet’s inclination will be damped at the same time (§4.2). If the inclination damping timescale is shorter than the disc precession timescale, we should not expect misaligned inner and outer discs since the disc will not have enough time to precess before the planet becomes coplanar. Figure 19 shows both the inclination damping timescale (Equation 19) and the disc precession timescale (Equation 28). The precession timescale has a very weak dependence on the planet inclination (cos $i_{p}$), while the inclination damping timescale sensitively depends on the planet inclination. Figure 19 shows that the inclination damping timescale decreases until a gap starts to form, and then the inclination damping timescale increases. If $\Sigma_{0}R_{p}^{2}/M_{*}$=0.001 and $\alpha=10^{-3}$, the inclination damping timescale for a $i=39^{o}$ planet is comparable or longer than the precession timescale, and significant misalignment between the inner and outer discs are possible. For less inclined planets, only massive planets ($q\gtrsim 0.01$) can cause significant misalignment before the planet becomes coplanar. When the disc is more massive or $\alpha$ is larger, the inclination damping timescale gets shorter and an even more massive planet is needed to cause misalignment between the inner and outer discs. On the other hand, if we can constrain that the planet in the misaligned disc is not very massive (e.g. 1 MJ), we can infer that the disc $\alpha$ is not large (e.g. $\alpha\lesssim 10^{-3}$) or the disc mass is small (e.g. $\Sigma_{0}R_{p}^{2}/M_{*}\lesssim 10^{-3}$).

6.2 Inferring Planet Properties From Observations

Based on observations, the geometry of the broken discs can be nicely constrained by studying its near-IR images and (sub-)mm molecular line channel maps (§5). We want to use this broken disc geometry to constrain the mass and inclination of the potential planet in the disk.

The planet’s inclination is relatively easy to derive if the inner and outer discs only experience the nodal precession around the system’s angular momentum vector. Under this condition (which can be violated in some cases as discussed in §6.4), the angle between the angular momentum vectors of the inner and outer discs is between 0 and 2$i$ where $i$ is the angle between the planet’s orbital vector and the angular momentum vectors of either the inner or outer disc (the inner and outer discs have the same $i$ initially and the value of $i$ is a constant during the inner disc’s nodal precession). Thus, if there is an observed system with a misalignment angle of $i_{obs}$ between the inner and outer discs, it indicates that the orbital plane of the potential planet is misaligned with either the inner or outer discs at an angle of $i_{p}\geq i_{obs}/2$.

The minimum planet mass to break the disc is shown in Figure 12. This mass increases with the planet’s inclination since a more inclined planet induces a shallower gap and also a slower precession of the inner disc. Observationally, if we find a broken disc but cannot constrain the misalignment angle between the inner and outer discs, the minimum planet mass is thus the mass at the $i=0^{o}$ limit in Figure 12.

We can derive this disc-breaking planet mass at the $i=0^{o}$ limit analytically since the elliptic integral ($K$) at $i=0^{o}$ is simply one. Combing the gap depth equation at $i=0^{o}$ (Equation 5) with Equations 28 and 38, we can derive a quadratic equation for $q^{2}$ and the solution for $q$ is shown in Figure 20. For a thin disc, a lower mass planet can break the disc and the gap depth is deeper. Since these gaps are normally very deep, Equation 5 can be simplified to $\Sigma_{gap}/\Sigma_{0}=25/K$. Then, the minimum planet mass to break the disc at the small $i$ limit can be further simplified to

[TABLE]

This simple formula agrees with Figure 20 very well.

6.3 Individual Sources

Here we can use our disc breaking condition (Equations 38 and 56, Figures 12 and 20) to study the known protoplanetary discs that have clear signatures of disc breaking.

HD 142527 is a large transitional disc around a 2 $M_{\odot}$ central star (Casassus et al., 2013). It has a large cavity ($\sim$140 au) and a noticeable asymmetric structure at the cavity edge. It also has a clear dark lane at its outer disc, which reveals that the inner disc is $\sim$70o misaligned with the outer disc (Marino et al., 2015). The dark lane coincides with the (sub-)mm intensity dimming inside the asymmetric structure at 140 au (Casassus et al., 2013). This dimming is consistent with Figure 13 where an asymmetric vortex is passing through the inner disc’s shadow and the temperature of the disc under the shadow is decreased. HD 142527 has a known 0.2-0.4 $M_{\odot}$ companion at 13 au (Biller et al., 2012; Christiaens et al., 2018). Such a high mass companion is fully capable of breaking the disc (Figures 12 and 20) and causing differential precession between the inner and outer discs. The companion’s orbit can be very eccentric (Lacour et al., 2016) and it can explain the observed large cavity (Price et al., 2018). If we assume that the companion’s eccentricity is low, the inner disc will precess for a full $2\pi$ angle at a timescale of $\sim$27 planetary orbits (Equation 28), which is $\sim$ 900 yrs with a companion at 13 au. Thus over decades, we may see the change of the disc shadow at the outer disc. Here we have used Equation 28 which has assumed that the inner disc size ($R_{d}$) is the same as the planet’s semi-major axis ($R_{p}$). If $R_{d}$ is much smaller than $R_{p}$, the precession timescale may be much longer.

HD 100453 is a protoplanetary disc with two prominent spiral arms (Wagner et al., 2015). The two spiral arms could be excited by the 0.2 $M_{\odot}$ companion (HD 100453B) which is 108 au away from HD 100453A (Dong et al., 2016; Wagner et al., 2018). On the other hand, this disc also shows dark lanes in the near-IR scattered light images, indicating a misaligned inner disc (Benisty et al., 2017). The misalignment angle between the inner and outer discs is $\sim$ 72o and the gap between the inner and outer discs is from 1 to 20 au. One explanation for the misalignment could be that there is a massive inclined planet inside the gap (at several au) and the planet mass can be constrained by Equation 56. On the other hand, if HD 100453B at 108 au is misaligned with the disc plane of HD 100453A, HD 100453B can also cause precession for both the disc within 1 au and the disc beyond 20 au. Since these two discs have different precession rates, it can lead to misalignment between them. Under this second scenario, massive planets are not needed at several au but HD 100453B should be more than 36o misaligned with either the inner or outer discs and we should not see the movement of the shadow during decades of observations. Using Equation 27 with $R_{d}$=1 au (Benisty et al., 2017), we can derive $\omega_{rb}/\Omega_{b}\sim 4\times 10^{-5}$ where $\Omega_{b}$ is the orbital frequency for HD 100453B. Equation 38 then suggests that $\Sigma_{gap}/\Sigma_{0}$ is $\lesssim 3\times 10^{-6}$ with $h$=0.05. This scenario may also have applications to other systems with both spirals (which may indicate that there is a massive outer companion) and disc shadows.

DoAr 44 is a protoplanetary disc around a 1.4 $M_{\odot}$ star (Bouvier & Appenzeller, 1992). It also has a dark lane at the outer disc and the inferred misalignment between the inner and outer discs is $\sim 30^{o}$ (Casassus et al., 2018). The inner and outer discs are separated by a gap from 5-15 au, based on the detailed modeling for its SED and Near-IR scattered light images. Assuming the disc viscosity $\alpha=10^{-3}$ and $h=0.1$, the minimum planet mass to break the disc is $q\gtrsim$8$\times 10^{-3}$ (Equation 56). The gap depth is thus $\Sigma_{gap}/\Sigma_{0}\lesssim 4\times 10^{-3}$ (Equation 5). This planet should be more than $15^{o}$ inclined with respect to both the inner and outer discs. The inner disc will precess for a full $2\pi$ angle at a timescale of $\sim$300 planetary orbits (Equation 28), which is $\sim$ 6000 yrs with a planet at 8 au. Thus, it will be difficult to detect the position change of the dark lane in the near future, unless the inner and outer discs are oriented in a way which makes $d\theta/dt\gg d\gamma/dt$ (Equation 54).

HD 143006 is a protoplanetary disc with gaps. Near-IR scattered light observations reveal asymmetry at the outer disc (Benisty et al., 2018): half of the disc is fainter than the other half. Such asymmetry is consistent with the shadow cast by a mildly inclined inner disc (Facchini et al., 2018). Assuming that a planet at $\sim$20 au is responsible for the tilt of the inner disc, Figure 12 and Equation 56 suggest that the planet mass needs to be larger than 2 $M_{J}$ in an $\alpha=10^{-3}$ and $H/R=0.04$ disc.

AA Tau shows signatures of a misaligned inner disc in both the dust continuum image and molecular line channel maps (Loomis et al., 2017). If the misalignment between the rings is caused by a misaligned planet, the planet needs to be at least in the Jupiter mass range (Equation 56 and Figure 12). AA tau also exhibits a sudden and long-lasting dimming event at 2011 (Bouvier et al., 2013) after 20 years of a constant brightness. If we assume that 20-200 years is the timescale for the inner disc’s precession due to a potential companion, the companion is massive and close to the central star (Equation 28 suggests that the companion needs to be within several au if $q=0.1$).

Overall, in order for a misaligned companion to break the disc and lead to differential precession between the inner and outer discs, the gap needs to be relatively deep and a massive giant planet ($q\gtrsim 0.001$) is needed. Thus, under this scenario, disc misalignment and shadowing should always be accompanied by deep and wide gaps (e.g. in Pre/Transitional discs), which is consistent with recent observational constraints (Garufi et al., 2018). On the other hand, if we can constrain that the companion mass in the gap of the misaligned disc is small (e.g. $\lesssim M_{J}$), the disc viscosity has to be very low (e.g. $\alpha\lesssim 10^{-4}$) for the disc to break.

6.4 Limitations And Comparison With Previous Works

One caveat in this work is that we do not consider the disc’s gravitational force on the planet and disc’s self-gravity. This excludes the possibility that the planet and disc might undergo secular resonances (Lubow & Martin, 2016; Owen & Lai, 2017; Matsakos & Königl, 2017) to generate large misalignment. On the other hand, secular resonances require a massive companion $\sim 0.01-0.1M_{\odot}$ at $\sim$10-100 au from the star and the circumbinary disc is also massive which dominates the angular momentum budget of the system (Owen & Lai, 2017). Thus, secular resonances may work for HD 142527 which has a known massive companion, but it may not work for ‘Pre/Transitional discs’ where we have not found such massive companions in the disc.

Our work also has implications to secular resonances. With the presence of such a massive planet in the disc to trigger secular resonances, disc breaking can naturally occur (Figure 12) and previous analytical studies which separate the system into the circumprimary disc, companion, and circumbinary disc can be justified.

We also limit the planet’s inclination to be less than 39o, so that the inner disc is not undergoing the Kozai-Lidov oscillation (Martin et al., 2014). If the disc is undergoing the Kozai-Lidov oscillation, the excited eccentric disc can be damped by hydrodynamical processes and the disc will settle to a state with an inclination less than $39^{o}$ and a low eccentricity (Fu et al., 2015). On the other hand, HD 142527 and HD 100453 may be undergoing the Kozai-Lidov oscillation now. These systems have disc misalignment around 70o, implying that the potential planet has a high chance to be misaligned with the inner disc at more than 40o to trigger the Kozai-Lidov oscillation.

We do not couple radiative transfer calculations with hydrodynamical simulations. Instead, we post-process the hydrodynamical simulations with Monte-Carlo calculations. This implicitly assumes that the disc has an infinite amount of time to respond to the stellar radiation. This may not be true in realistic protoplanetary discs, and including both the radiative transfer and the hydrodynamics simultaneously may lead to interesting phenomenon such as generating spiral waves (Montesinos et al., 2016).

Arzamasskiy et al. (2018) found that, if a planet cannot induce a gap in the disc, the planet can most efficiently warp the disc when the planet’s inclination angle is $\sim$2-3 $h$ with respect to the disc. Then, the inclination of the disc warp can reach $\sim 4q$ (their Figure 8). Arzamasskiy et al. (2018) have also carried out Monte-Carlo radiative transfer calculations for an inclined Jupiter mass planet in a $h=0.1$ disc and no disc shadowing is observed. This is expected because the planet mass just reaches the thermal mass and a shallow gap is induced by the planet. In order to break a $h/r=0.1$ disc with a $q=10^{-3}$ planet, Equation 38 suggests that $\Sigma_{gap}/\Sigma_{0}<6\times 10^{-5}$, which requires a very small numerical viscosity and a long simulation time. Since the planet can not break the disc, the planet induces a warp with the $4q$ disc inclination (0.004 in radians or 0.3o with $q$=0.001) which is too small to cause any disc shadow.

Recently, Nealon et al. (2018) have carried out simulations to show that a Jupiter mass planet can cause a $10^{-3}-10^{-2}$ misalignment between the inner and outer discs while a 6.5 Jupiter mass planet can cause a $10^{-2}-10^{-1}$ misalignment between the inner and outer discs. Such warp is roughly consistent with the angle of $4q$ derived in Arzamasskiy et al. (2018). Thus, we should not expect significant disc shadowing from simulations in Nealon et al. (2018) ($0.1$ in radians is only 6o). Disc breaking has not been observed in Nealon et al. (2018) which shows some constant misalignment angle between the inner and outer discs. With a more massive planet (like our simulation THINP10I39AM3), the disc may start to break.

7 Conclusion

Motivated by recent observations that many protoplanetary discs with gaps and holes show signs of disc warping and breaking, we study the interaction between a misaligned massive planet and the protoplanetary disc.

• •

We analytically derive the relationship between the mass of the misaligned planet and the induced gap depth/width (Equations 12 and 13), by replacing the planet mass with an effective planet mass (Equation 11).

• •

We derive the migration and inclination damping rates for misaligned planets that could induce gaps in discs (Equations 14 to 19).

• •

We provide the new disc breaking condition for a disc having a gap (Equations 37 and 38): the sound crossing time through the inner disc is longer than the product of the precession timescale and the square root of the gap depth.

• •

Using 3-D hydrodynamical simulations with mesh-refinement, we carry out a series of simulations to study the interaction between the misaligned planet and the disc. We find that 10 grid points per scale height is needed for studying the disc secular evolution. A higher order reconstruction scheme (e.g. PPM) can lower this requirement to 5 grid points per scale height.

• •

The simulations show that when the gap is deep enough, the inner and outer discs precess independently so that the disc breaks. Even after the disc breaking, the circumplanetary disc around the planet seems to keep its alignment with the outer disc instead of the inner disc.

• •

If a gap is not induced, the planet’s migration and inclination damping rates are similar to the previous study for low mass planets. When a gap is induced, both rates decrease as the gap gets deeper. The simulation results are consistent with analytical expectations (Equations 14 to 19) as long as the disc does not break.

• •

The induced gap depth from simulations is consistent with the analytical expectation very well (Equations 12 and 13). For simulations that show disc breaking, our derived disc breaking condition (Equation 37, Figure 11) also agrees with simulation results very well.

• •

We generate near-IR scattered light images and (sub-)mm CO velocity channel maps for broken discs, by importing the simulated broken disc structure into Monte-Carlo radiative transfer calculations. In the near-IR scattered light images for inviscid simulations, we see that the vortex at the outer gap edge is traveling through the shadow cast by the inner disc. When the vortex and the shadow overlaps, the vortex looks like it splits into two vortices. The temperature of the disc region that is under the shadow decreases, potentially leading to a change of the (sub-)mm intensity there.

• •

The disc shadow in near-IR scattered light images is also affected by the disc surface density. A higher surface density means that the scattering surface is higher and the shadow reveals this bow-shaped flaring surface.

• •

We observe that the shadow rotates at a different orbital frequency from the precession frequency of the inner disc. The relationship between these two frequencies is derived (Equation 54).

• •

CO channel maps clearly reveal the geometry of the broken disc. If the molecular line is optically thick within the gap, the channel maps show a twisted velocity structure.

• •

If observations reveal the broken disc geometry, we can use this geometry to constrain the planet mass (Figures 12 and 20) and the inclination. The minimum planet mass which can break the disc is derived analytically (Equation 56). We use this criterion to explore the potential companions in HD 100453, DoAr 44, HD 143006, and AA Tau.

• •

On the other hand, if we can constrain that the companion mass in the gap of the misaligned disc is small (e.g. $\lesssim M_{J}$), the disc viscosity has to be very low (e.g. $\alpha\lesssim 10^{-4}$) for the disc to break. This provides some indirect constraints on the disc viscosity.

• •

Overall, in order for a misaligned companion to break the disc and lead to differential precession between the inner and outer discs, the gap needs to be relatively deep (e.g. in Pre/Transitional discs), and a planet with at least a Jupiter mass is needed even under some extreme disc conditions. Our study supports the scenario that massive planets are present in Pre/Transitional discs that have disk shadows, where the massive planets (probably larger than $M_{J}$) are responsible for both the wide gaps/cavities and the inner disc misalignment causing disc shadowing. These disks with shadows are thus prime targets for directly searching exoplanets in protoplanetary discs.

Acknowledgements

Z. Z. thank Jaehan Bae for sharing his scripts to generate CO channel maps. Z. Z. thank Steve Lubow for very helpful discussions. Z. Z. acknowledges support from the National Aeronautics and Space Administration through the Astrophysics Theory Program with Grant No. NNX17AK40G and Sloan Research Fellowship. Simulations are carried out with the support from the Texas Advanced Computing Center (TACC) at The University of Texas Austin through XSEDE grant TG-AST130002, and the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.