Scattered light shadows in warped protoplanetary discs

TL;DR

This study uses 3D hydrodynamic simulations to explore how a misaligned planet can warp a protoplanetary disc, casting observable shadows and brightness variations, with implications for interpreting observations like those of TW Hya.

Contribution

It demonstrates the impact of a small inclined planet on disc structure and observable shadows, providing a link between planet-induced warps and scattered light images.

Findings

A planet of ~6 Jupiter masses and 2° inclination can cast a significant shadow.

Warped outer discs cause azimuthal brightness variations at large radii.

The model suggests inner disc precession alone may not explain TW Hya's features.

Abstract

Three-dimensional hydrodynamic numerical simulations have demonstrated that the structure of a protoplanetary disc may be strongly affected by a planet orbiting in a plane that is misaligned to the disc. When the planet is able to open a gap, the disc is separated into an inner, precessing disc and an outer disc with a warp. In this work, we compute infrared scattered light images to investigate the observational consequences of such an arrangement. We find that an inner disc misaligned by a less than a degree to the outer disc is indeed able to cast a shadow at larger radii. In our simulations a planet of around 6 Jupiter masses inclined by around 2 degrees is enough to warp the disc and cast a shadow with a depth of more than 10% of the average flux at that radius. We also demonstrate that warp in the outer disc can cause a variation in the azimuthal brightness profile at large radii. Importantly, this latter effect is a function of the distance from the star and is most prominent in the outer disc. We apply our model to the TW Hya system, where a misaligned, precessing inner disc has been invoked to explain an recently observed shadow in the outer disc. Consideration of the observational constraints suggest that an inner disc precessing due to a misaligned planet is an unlikely explanation for the features found in TW Hya.

Equations

$$F(\varphi )=A\cos \left[\varphi -C\right]+1,$$

$${\omega }_{\mathrm{prec}}=-\frac{3Gm}{4a^3}\frac{{\int }_{x_{\mathrm{in}}}^{x_{\mathrm{out}}}\Sigma R^{3}dR}{{\int }_{x_{\mathrm{in}}}^{x_{\mathrm{out}}}\Sigma {\Omega }_k{R}^{3}dR}\cos {\beta }^{\prime }.$$

$$P=\frac{7.21\times 10^7}{{\mu }_c\cos i_{\mathrm{c}}}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{-1/2}{\left(\frac{r_{\mathrm{disc}}}{\text{1 au}}\right)}^{-3/2}{\left(\frac{a_{\mathrm{c}}}{\text{300 au}}\right)}^3\text{yr},$$

$$\Sigma (R)={\Sigma }_0{\left(\frac{R}{R_0}\right)}^{-p}\left(1-\sqrt{\frac{R_{\mathrm{in}}}{R}}\right),$$

$$M_{\mathrm{disc}}={\int }_{R_{\mathrm{in}}}^{R_{\mathrm{p}}}2\pi R\Sigma (R)dR.$$

$$t_{\nu }={\left(\frac{H}{R}\right)}^{-2}\frac{1}{\alpha {\Omega }_{\mathrm{k}}},$$

Full text

††pubyear: 2015††pagerange: Scattered light shadows in warped protoplanetary discs–References

1 Introduction

Evidence is mounting that misalignments and warps in protoplanetary discs are common. Strongly misaligned discs present the most compelling examples, where the inner part of the disc is misaligned to the outer disc. The inner disc casts unique shadows on the outer disc, constraining the relative inclination between the two (Min et al., 2017). Relative inclinations have been measured from  30*∘* for DoAR 44 (Casassus et al., 2018) to 45*∘* for AA Tau (Loomis et al., 2017), 70*∘* for HD 142527 (Marino et al., 2015), 72*∘* for HD 100453 Benisty et al. (2017) and up to 80*∘* for HD 100546 (Walsh et al., 2017). Numerical simulations suggested that such an arrangement may be generated through gravitational interaction with an inclined binary. For example, a circumbinary disc may be broken into two discs and the inner can cast asymmetric shadows onto the outer (Facchini et al., 2013; Facchini et al., 2018; Zhu, 2018). For the specific case of HD 143006, Benisty et al. (2018) demonstrated that such a model reproduces most of the observed features. Using HD 142527 as a case study, Price et al. (2018b) suggested that spiral arms, a cavity and shadows observed in some transition discs could be explained with the presence of an inclined binary. Characteristic disc structures including warps may also be generated through the interaction with an inclined external perturber, as in RW Aurigae system (Dai et al., 2015), MWC 758 (Dong et al., 2018) and as shown by Cuello et al. (2018). Additionally, strongly misaligned discs may also be formed during independent chaotic accretion episodes (Bate, 2018).

Small to moderate misalignments have been invoked at the innermost edge of the disc to explain observations of dipper stars, e.g. Ansdell et al. (2016) and recently by Pinilla et al. (2018). These may be driven by stellar magnetic fields that are misaligned to the disc, as in AA Tau (Donati et al., 2010; Bouvier et al., 1999) and LkCa 15 (Alencar et al., 2018). Further evidence of misaligned discs is motivated by TW Hya, where recent infrared (IR) observations identified a variation of $\sim 20\%$ in the surface brightness of the outer disc that moved with a speed much faster than Keplerian at that radius (Roberge et al., 2005; Debes et al., 2017; Poteet et al., 2018). Such an asymmetry is difficult to explain as TW Hya is considered to be an isolated, single star system (Mamajek, 2005; Gagné et al., 2017)— broadly ruling out interactions with external perturbers or the interaction of the disc with a misaligned binary. Debes et al. (2017) suggested that the feature in TW Hya may instead be caused by a misaligned, inner disc that is casting a shadow on the outer disc. In this interpretation the motion of the inner disc is governed by a misaligned planet that causes the disc to be inclined and rigidly precess, casting a shadow that moves with the precession rate of the inner disc. This model is consistent with existing observations; a small warp or misalignment in the inner region has been suggested (Qi et al., 2004; Pontoppidan et al., 2008; Hughes et al., 2011, $\gtrsim 10$ au) and a planetary companion may explain the inner disc cavity (Calvet et al., 2002; Andrews et al., 2016). Additionally, a small warp is consistent with CO kinematics in the inner part of the disc (Rosenfeld et al., 2012).

Numerical simulations have confirmed that this arrangement may be possible, as the disc structure can be altered by the influence of an inclined planet (Nealon et al., 2018; Zhu, 2018). In the case that the planet is able to open a gap, the initial disc is separated at the orbit of the planet into an inner and outer disc. A warp forms across the orbit of the planet, where a misalignment develops between the inner and outer disc (Xiang-Gruess & Papaloizou, 2013; Bitsch et al., 2013; Nealon et al., 2018). Xiang-Gruess & Papaloizou (2013) demonstrated that the inner (and, in their case outer) disc tilt toward the orbit of the inclined planet. This behaviour is confirmed by Bitsch et al. (2013) who found that the inner disc can achieve an inclination greater than that of the planet orbit. Both of these works agree that the disc responds most rapidly when the planet is at moderate inclinations (at inclinations $\sim$ the aspect ratio), rather than larger inclinations. Recently, Nealon et al. (2018) conducted similar simulations with a larger outer disc radius (50 au rather than $\sim$20 au). We showed that the larger outer disc moves less under the influence of the misaligned planet, leading to a larger relative misalignment between the inner and outer disc. In the case that protoplanetary discs are larger than $\sim$20 au (Andrews & Williams 2007 and Ansdell et al. 2018, but see Cieza et al. 2019 and van Terwisga et al. 2018), this suggests that the inclinations between the inner and outer disc found in the previous works represent underestimates.

Scattered light observations of protoplanetary discs are the most important diagnostic to investigate the structural properties of protoplanetary discs. The scattering emission from small dust grains at the disc surface layer dominates at short wavelengths, roughly from optical to near-infrared (NIR) part of the spectrum. Small dust grains carry most of the opacity at optical and near-infrared wavelengths, and therefore control the stellar absorption and the disc’s temperature (Hartmann, 2009). Additionally, these dust grains are expected to be well coupled with the gas in the inner and intermediate disc regions and populate the full vertical extent of the gas disc. Therefore, such observations can be considered a fundamental probe of the disc vertical structure (see the review of Watson et al., 2007) and can be used to infer the disc inclination, the shape of the surface and detect features at the surface layer (e.g. induced by gaps carved by planets, Jang-Condell & Turner, 2012). The observational consequences of a low inclination warp in the inner disc have been explored by Arzamasskiy et al. (2018), where they compared the IR scattered light profile in the cases of an aligned and misaligned planet. In their intensity maps, asymmetries due to the warp from the misaligned planet were only observed in the inner regions of the disc. However, the duration of their simulation is shorter than in other works (e.g. Marzari & Nelson, 2009; Xiang-Gruess & Papaloizou, 2013; Bitsch et al., 2013; Nealon et al., 2018), the gap opened by the planet is not fully formed and the planet is kept at on a fixed orbit. This scenario is distinct from that invoked from observations of TW Hya, where the inner disc must be cleanly separated from the outer disc such that it can precess differentially.

In this paper we measure the observational consequences in scattered light for a disc warped by an interaction with a misaligned planet. We use the radiative transfer code mcfost (Pinte et al., 2006; Pinte et al., 2009) to generate scattered light images from the simulations presented in Nealon et al. (2018) (hereafter NDAMN18). In Section 2 we summarise these simulations and our use of mcfost. In Section 3 we present the scattered light results and examine the azimuthal brightness profile in the intermediate and outer disc. In Section 4 we consider the observed features of TW Hya in the context of our model. In Section 5 we discuss these results in the context of the limitations of our simulations and conclude in Section 6.

2 Numerical method

2.1 Three-dimensional hydrodynamic simulations

In NDAMN18 we conducted simulations using the smoothed particle hydrodynamics (SPH) code phantom, a code particularly well suited to examining warps and disc-planet interactions (Price et al., 2018a). The simulations from Section 4 of that work examined the behaviour of an accretion disc with a planet on a misaligned orbit. In the suite of simulations presented the properties of the disc were held constant; the disc extended from 0.1 to 50 au with a disc mass of $5\times 10^{-3}$ M*⊙* and was simulated with $3.4\times 10^{6}$ particles. The surface density profile in each case was initially set with a tapered inner edge and a power-law profile with an index of $p=-1$. The aspect ratio was set to $H/R=0.05(R/R_{0})^{1/4}$ and a vertically isothermal disc assumed with the sound speed profile $c_{\rm s}=c_{\rm s,0}(R/R_{0})^{-1/4}$, where $R$ is the cylindrical radius and $R_{0}$=1 au for the simulations we are considering. As the disc evolves, the temperature of the particles remains a function of radius and is assigned according to $c_{s}(R)$. The viscosity was modelled with the Shakura & Sunyaev (1973) $\alpha$ parametrisation, with $\alpha=10^{-3}$ for all the particles. We refer the interested reader to NDAMN18 for further details.

A planet was placed in each disc at a distance of 5 au from the star with a mass of 0.13, 1.3 or 6.5MJ (corresponding to 0.1, 1.0 and 5.0$\times$ the thermal mass) and an orbital inclination of $2.15^{\circ}$, $4.30^{\circ}$ or $12.89^{\circ}$ (0.5, 1.0 and 3.0$\times$ $H/R$, where $H$ is the scale height). This planet was then allowed to migrate and accrete mass due to the interaction with the surrounding disc for 200 orbits at the initial planet orbital radius. For those planets that were able to open a gap, we found (in agreement with previous studies, see Xiang-Gruess & Papaloizou, 2013) that the initial disc was separated into two discs. The inclination of the planets’ orbit caused a change in the relative orientation of the inner and outer disc, with more massive and strongly inclined planets driving a larger difference.

Figure 1 shows the density cross-section $t=100$ planet orbits of our representative simulation (with a 6.5MJ planet inclined at $12.89^{\circ}$). Here the dashed line represents the initial position of the disc mid-plane. Although the outer disc (with $R>5$ au) has not moved appreciably, the inner has tilted towards the inclined planet orbit. This rotation is accompanied by the inner disc twisting as it undergoes differential precession. The outer disc also warps slightly in response to the misaligned planet (see Figures 8 and 9 of NDAMN18) and traces of the spiral shocks from the planet can be seen, e.g. around $y\sim-8$ au, $y\sim 10$ au in the cross section in Figure 1.

In our simulations, the most massive planet (6.5MJ) was surrounded by a weak tidal stream that was most clearly visible during the initial phases where the gap is not fully carved. As this could potentially complicate any shadows due to the inclination of the inner disc, in this work we only consider snapshots where the planet is on the far side of the disc. We also supplement our parameter study with three new simulations. These have exactly the same properties as the previous set except that we consider a planet mass of 2.6MJ (corresponding to 2.0$\times$ the thermal mass). In line with our previous results the massive planet opens a gap and a warp between the inner and outer disc and driving precession of the inner disc.

For simplicity, in this work we use the $t\sim 100$ orbit snapshots of the simulations presented in NDAMN18. At this time the gap opened by the planet is well established, as the relative tilt between the inner and outer disc is almost at maximum (e.g. see Figure 8, NDAMN18, ). The resolution of the simulations decreases at late times as the particles in the inner disc are accreted onto both the star and planet, but at $t\sim 100$ orbits it has not decreased appreciably — the shell averaged smoothing length to scale-height ratio $\langle h\rangle/H\approx 0.65$ at $R=3$ au and $<0.2$ for $R>10$ au. This choice means we do not take into account any temporal evolution of the discs, noting that the warp amplitude decreases slightly after $t=100$ orbits in all of our simulations.

2.2 Radiative transfer calculations

We use the continuum radiative transfer code mcfost to generate scattered light images from our simulations (Pinte et al., 2006; Pinte et al., 2009). mcfost is particularly well suited for use with SPH simulations because it employs a Voronoi mesh rather than a cylindrical or spherical grid (Camps et al., 2015). This method allows the mesh to be made without interpolation, as it is generated directly from the particle positions. To accommodate regions of the simulation that have lower numerical resolution, mcfost implements a check based on the comparison between the size of the cell and the smoothing length of the particle. In the case that a dimension of any cell is three times larger than the smoothing length, mcfost considers the region outside three smoothing lengths of the particle within that cell to be optically void.

We focus on scattered light observations that are sensitive to the distribution of micron sized dust grains. As these grains are very well coupled to the gas, we can use the gas distribution in our simulations directly. From our simulations we estimate a Stokes number of $\sim 1.2\times 10^{-3}$ in the upper layers for $\sim\mu$m grains at 5 au, confirming this assumption is valid. The total mass of dust in the disc is calculated from the mass of gas in the disc (directly measured from the simulation) and an assumed dust-to-gas ratio. We assume that each of the dust grains are spherical and homogeneous (according to Mie theory), that the dust and gas is in thermal equilibrium and the dust opacities are independent of temperature. We assume a dust-gas ratio of 100, where the dust grains have a power-law size distribution with an exponent of -3.5 (Mathis et al., 1977) between 0.03 and 1000$\mu$m using 100 grain sizes.

To improve efficiency the computation of scattered light is broken into two steps. First, the temperature of the disc is calculated for each of the Voronoi cells. Photon packets are emitted from the central star using a Monte Carlo approach and the path of each photon packet is governed by scattering, absorption and re-emission processes (see Pinte et al., 2006). Second, the scattered light is calculated by emitting photon packets randomly from the disc and central star. In this second step the photon packets are always scattered (i.e. there is no absorption) but the energy of the packet is multiplied by the albedo at each scattering event. Images are then produced with a ray-tracing method, where we integrate the radiative transfer equation using the specific intensities and temperatures computed by the Monte-Carlo runs. We refer to Pinte et al. (2006) and Pinte et al. (2009) for a more detailed description of mcfost.

Unless otherwise stated, we use the following values for our mcfost calculations. To compute the temperature structure and the scattered light we use $10^{8}$ photon packets and our grid is set by the Voronoi mesh. We set our image to have a size of 120 au and define the angular momentum vector of the outer disc as direction of the viewer (that is, the viewer is ‘face-on’ to the outer disc). We do not include any contribution from viscous heating. The luminosity of the star is calculated from 1 Myr Siess isochrone (Siess et al., 2000), assuming a fiducial 1M*⊙* star. This corresponds to an effective temperature of $T_{*}=4278$K and a luminosity of $L_{*}=0.16$L*⊙*. The scattered light images have been computed at a wavelength of $\lambda=1.6\mu$m and convolved with a 2D Gaussian to smooth across 1.8 au (which corresponds to a FWHM of 30mas at a distance of 60 pc, Gaia Collaboration et al. 2018). We chose these values as they represent the angular resolution of TW Hya using a typical instrument, e.g. SPHERE or GPI. The central 5 au of the disc is masked in figures (grey inner circle) to simulate the effect of a coronograph.

2.3 Measuring the surface brightness as a function of angle

Using the NICMOS and STIS instruments on the Hubble Space Telescope (HST), Roberge et al. (2005), Debes et al. (2017) and Poteet et al. (2018) identified subtle fluctuations in the IR surface brightness of TW Hya by considering the flux as a function of position. In order to measure this, they separated the disc into concentric annuli (describing the distance from the star) as well as azimuthal wedges (describing the angle around the disc). This allows the flux in each radial annulus to be described as a function of the azimuthal angle. We follow the method outlined by Debes et al. (2017) in order to describe the flux from our simulations in a comparable manner.

We discretise the outer disc into wedges where the annuli are delimited by 5, 15, 30 and 50 au and the azimuthal slices have an angular size of $5^{\circ}$. In each single wedge, the flux from all pixels is averaged and the standard deviation quantifies the noise in the flux at that location in radius and azimuthal angle. The flux across each radial annulus is averaged and we use it to normalise the flux in each wedge in that radial annulus. This means that the amplitude of the surface brightness profiles plotted as a function of azimuthal angle represent the variation from the average brightness in that radial annulus. Following the method of Debes et al. (2017), we assume this azimuthal profile can be fitted with a sinusoidal function of the form

[TABLE]

where $A$ and $C$ are unconstrained and $\phi$ is the angle measured from the positive $x$-axis in an anti-clockwise direction. This form is scaled to the average brightness around the annulus, where $A$ determines the magnitude of the variation compared to the average and the position angle $C$ is the angular location of the strongest variation.

The fit is achieved using a least-squares curve fitting tool scipy.optimize (Jones et al., 2001). When conducting the fit, each measurement of the flux at a specific azimuthal wedge is weighted by the standard deviation of the flux in that wedge. Each parameter ($A$ and $C$) is returned with an associated standard deviation. We use the standard deviation of the parameter $A$ as an estimate for our accuracy of the fit on the amplitude of the surface brightness profile.

3 Synthetic scattered light observations

We computed the scattered light images for our 12 simulations described above using mcfost. Figure 2 shows an example of the resulting total intensity image for the simulation with a 6.5MJ planet initially inclined at $12.89^{\circ}$. It is apparent that the upper, right side of this image is more strongly illuminated than the lower, left side. This brightness feature is consistent with the orientation of the misaligned inner disc casting a shadow towards the left, lower side of the disc in Figure 2. We examine the properties of this shadow in the following sections.

3.1 Azimuthal brightness profile

The left panels of Figure 3 show the density rendering of our simulation (upper), the corresponding scattered light image (middle, as in Figure 2) and the azimuthal brightness profile (lower). The surface brightness has been measured across 5-15 au, 15-30 au and 30-50 au as indicated on the scattered light image. The surface brightness profile shows asymmetries are present at all radii outside of the planet orbit, with the deepest shadow occurring near the inner disc. The profile of the shadow is quite noisy for $\phi<180^{\circ}$, likely due to low resolution in the particle distribution close to the star being strongly illuminated.

In order to test the effect of the inner disc, we remove all the gas inside of 15 au as a post-processing step and re-run the radiative transfer calculation. The density rendering, scattered light image and surface brightness profile are shown in the right panels of Figure 3. We note that the same profile shown in Figure 3 is achieved by only removing the gas inside of 5 au, but with 15 au we can be confident there are no effects due to warping at the inner edge of the outer disc. By eliminating all gas in the inner region, we remove any structures that could cast a shadow such that the outer disc is uniformly illuminated by the star. The surface brightness still shows significant but much smaller variations in the azimuthal profile (lower, right panel). This variation is due to the influence of the misaligned planet on the structure of the disc; the planet induces a warp in the outer disc and this causes one side of the disc to be illuminated more strongly. Because the illumination of the disc depends on the differential twist between the source of the illumination (e.g. the star) and the radius of the disc that is being illuminated, this affect appears stronger at larger radii. The lower, right panel of Figure 3 demonstrates this, with the variation between 30-50 au stronger than in the intermediate disc. Due to warp propagation, we note that at later times in our simulation the warp amplitude is lower in the outer disc.

To confirm this variation in brightness in the outer disc is exclusively caused by differential twist in the disc driven by the planet, we first repeat the radiative transfer but force the scattering to be isotropic. Here we recover the same variation in the surface brightness, confirming that this effect is due to the preferential illumination of the disc. Second, we repeat this calculation at $t=40$ orbits when the inner disc is inclined but the outer disc has not yet been able to twist up (see Figures 8 and 9, NDAMN18, ). At this earlier time when the outer disc is not yet warped, we find no variation in the azimuthal brightness profile.

Our simulations with a misaligned planet therefore demonstrate that a non-axisymmetric azimuthal brightness profile can be driven by two mechanisms. First, the inner misaligned disc can cast strong shadows that emanate from the inner region. Second, the outer disc warp causes some part of the disc to be preferentially illuminated with the strongest variation occurring at large radii. Understanding the shadow in TW Hya requires us to consider both of these effects together.

3.2 Amplitude of the brightness variation

We quantify the variation in the surface brightness profile observed in Figure 3 using the curve fitting procedure outlined in Section 2.3. Figure 3 suggests that a variation could be caused by either shadowing from the inner disc or warping of the outer disc, so we present fits for both of these. Figure 4 shows the azimuthal brightness for both cases with the variation due to the misaligned inner disc in the left panel and the variation measured in the outer disc in the right panel. Both are shown with the standard deviation in the flux measurement as a shaded region. We measure the variation from the misaligned inner disc between 5-15 au, where there is little confusion from the outer disc warp and the shadow is strongest. We also measure the variation due to the outer disc warp between 30-50 au, where it is strongest (and after eliminating the inner disc, as described above). Both measurements thus represent an upper limit on the amplitude of the shadow for $t=100$ orbits and should be considered as demonstrations of the effect of the warp.

For a 6.5MJ planet inclined at $12.89^{\circ}$, the intensity of the shadow cast by the misaligned disc between 5-15 au is best fit with an amplitude of 78%$\pm$3% (compared to the mean in this annulus). For the same simulation, in the outer disc between 30-50 au we measure a variation in the surface brightness caused by the outer disc warp of 25%$\pm$9%. It is clear that the fit to the shadow in the left panel of Figure 4 is not a good representation of the actual profile cast by the shadow. This is particularly apparent for azimuthal angles $\lesssim 120^{\circ}$, where the profile is dominated by noise likely due to the strong illumination of low resolution particles. However we chose a simple sinusoid function because of its simplicity and also to allow direct comparison to the observations in Debes et al. (2017) and Poteet et al. (2018), as they used this form.

We repeat this measurement of the variation in the surface brightness due to either the misaligned inner disc or warp in the outer region for each of our 12 simulations. Figure 5 shows the amplitudes from each of these fits with their associated uncertainties, plotted as a function of the relative inclination between the inner and outer disc (i.e. the maximum projected inclination, calculated as the angle between their average angular momentum vectors). The depth of the shadow driven by the inner disc (left panel, Figure 5) is largest for more massive planets. This agrees with our previous findings where the massive planets are able to open a gap more cleanly, allowing a larger difference in the orientation between the inner and outer disc (NDAMN18). For those planets not able to cleanly separate the inner and outer disc (in our simulations, planets with $0.13$MJ), we find that the amplitude of the shadow depths measured are consistent with zero.

Figure 5 thus suggests that an appreciable shadow (i.e. more than $\sim$10%) is only cast by planets much more massive than the gap-opening mass. For an aligned planet, the gap-opening mass is approximately the same as the thermal mass of 1.3MJ (Crida et al., 2006; Baruteau et al., 2014). With the parameters used in our simulations, the limit for generating a brightness variation corresponds to $\gtrsim 6$MJ planets (at inclinations $\gtrsim 2^{\circ}$) that drive misalignments of only a few degrees between the inner and outer disc. We find a similar limit exists on the amplitude of the variation driven by the warp in the outer disc (right panel, Figure 5); planets that are more massive than the gap clearing mass are able to drive appreciable variations in the brightness profile in the outer disc. Additionally, in NDAMN18 we showed that the warp created across the planet orbit for a massive planet at low/moderate inclinations was mainly caused by differential precession of the inner disc. These shadows therefore represent the misalignment due to the precession of the inner disc, rather than the inner disc being strictly tilted away from the outer disc.

3.3 Effect of viewing inclination

In the above results we have assumed the viewing inclination is parallel to the angular momentum vector in the outer disc (i.e. the outer disc is face on). However, as the variation in surface brightness is a geometric effect it is likely to change when viewed at other inclinations. A cursory examination of the disc at different inclinations shows that this is indeed the case with shadows only clear for viewing inclinations less than about $30^{\circ}$. At inclinations larger than this, brightness profiles are complicated by the dust scattering function and flaring structure.

We thus consider the extreme case of viewing the disc edge-on (i.e. the viewing inclination is $i=90^{\circ}$ and the angular momentum vector of the outer disc is perpendicular to the line of sight). As our simulation does not go for a full precession of the inner disc, we emulate precession by rotating the disc through $90^{\circ}$, $180^{\circ}$ and $270^{\circ}$, measured clock-wise from the orientation in Figure 3. When the disc is subsequently turned edge on, the inner disc is then at four different orientations, as it would be across a full precession.

Figure 6 shows the inclined view of the entire disc. There is a flux asymmetry between the upper and lower disc surfaces depending on the rotation of the disc, with the lower disc surface brighter with an viewing rotation of $0^{\circ}$ and the upper disc surface at a rotation of $180^{\circ}$. To confirm whether this asymmetry is caused by the misaligned inner disc or the warp in the outer disc, we repeat this calculation but remove the gas inside of 15 au. As before, this removes any shadowing that could be caused by the inner misaligned disc such that any changes in flux should only be caused by the warp in the outer disc. We find that the disc is brighter in absolute flux (as the remaining disc is now more strongly illuminated) but that the same flux asymmetry pattern appears. This comparison suggests that a flux asymmetry between the upper and lower surfaces of an edge-on disc is a potential signature of a disc warp. If the warp is steady, in observations these asymmetries would evolve according to the precession time-scale of the warp in the outer disc. We note that similar asymmetries have already been observed on much shorter time-scales, e.g. HH 30 (Stapelfeldt et al., 1999).

3.4 Temperature structure

While our simulations in NDAMN18 assume the disc is vertically isothermal, mcfost calculates the temperature structure due to the illumination by the star. To quantify the effect of our assumption, Figure 7 compares the temperature calculated by mcfost for each Voronoi cell in the disc mid-plane with the temperature predicted in our simulation. Here we have discretised the Voronoi mesh into annuli (defined by cylindrical radius) and by comparing with the azimuthally averaged profile of the angular momentum, extracted and averaged the temperature of only those cells that are within a distance of $|z|<0.05$ from the mid-plane. Thus the profile in Figure 7 represents the temperature of the (sometimes misaligned) disc mid-plane. Outside of the planet orbit (at $R\gtrsim 10$ au) the radial temperature gradient is consistent with that assumed in the hydrodynamical simulations (i.e. $T(R)\propto R^{-1/2}$). Around the planet orbit (at 5 au) the spiral density waves created by the planet are identified in the temperature structure. The disc inside the planet orbit and around the planet is warped and not uniformly illuminated by the star, hence we do not expect good agreement in this inner region.

The offset in Figure 7 demonstrates that the temperature assumed in the simulation is hotter than that calculated using mcfost, by a factor of $\sim 4$. This offset corresponds to a difference in the sound speed in the disc of $\sim 2$ such that communication (via wave propagation) is occurring about twice as fast than it would if radiative transfer were included. A change in the sound speed (or equivalently, the thickness of the disc) may prevent the planet carving a gap or change the precession behaviour. However, for the parameters we have simulated in NDAMN18 we estimate that such a small change does not affect the relationship between the sound crossing and precession time-scales. As this relationship determines the evolution of the disc (in the event the planet is massive enough to open a gap, whether the inner disc will precess with respect to the outer disc), we thus expect that a cooler disc will exhibit similar dynamical evolution. Additionally, the disc temperature calculated by mcfost is primarily determined by the luminosity of the star, which in turn depends on the input parameters of its age and mass. A higher luminosity star would reduce the discrepancy (our adopted value is $L_{*}=0.16$ L*⊙*).

Figure 8 shows the azimuthal temperature profile between $9.9<R<10.1$ au, with the cells coloured by their azimuthal angle (such a narrow radial annulus is only chosen for plotting convenience, we find the same results with larger annuli). Significant departures from the locally isothermal approximation occur at only more than two scale-heights from the mid-plane.

Additionally, the striking pattern on either side of the mid-plane demonstrates the effect that shadowing has on the temperature structure in the disc. The upper and lower surface of the disc have opposing temperatures at the same azimuthal position, such that one side of the disc is warm while the other is cool. This temperature differential is driven by the shadowing due to the misaligned inner disc and will affect the dust thermal emission and disc scale-height of the disc if taken into account during our simulation.

4 Casting shadows: an application to TW Hya

Recent observations of TW Hya have identified an azimuthal surface brightness variation that is proposed to be created by an inclined inner disc driven by a misaligned planet. From observations, this shadow has two defining features; the speed at which it rotates ($22.7^{\circ}$/yr) and the radial location (outside of $\sim$50 au) (Roberge et al., 2005; Debes et al., 2017; Poteet et al., 2018). In Section 3 we have shown that an inner, misaligned disc generated by an inclined planet is capable of casting shadows with similar properties. We thus seek to test our model by applying it to the shadow observed in TW Hya. Indeed, if this model is consistent with the observations it will place constraints on the location and mass of the misaligned planet. Based on the features of the shadow observed in TW Hya, we consider the location, speed, depth of the shadow and secular evolution of the disc in turn.

4.1 Location of the shadow

The shadow identified in TW Hya appears only in the outer disc, where $R\gtrsim 53$ au. Debes et al. (2017) suggest that the absence of the shadow at intermediate radii may be caused by a few different mechanisms. First, they suggest the relative geometry between the inner and outer disc may result in a shadow only being cast in the outer disc. With the arrangement proposed by Debes et al. (2017), the shadow is likely to appear as in the left panels of Figure 3 where we find the shadow is weaker in the outer disc. Second, they propose that the inner rim of the outer disc is puffed in such a way as to cast a shadow that masks any azimuthal variation (as in Dullemond et al., 2001; Dong, 2015). However, the azimuthal temperature profile in Figure 8 demonstrates that if there is a misaligned inner disc it will drive an asymmetry in the temperature structure at the inner edge of the outer disc. This suggests that the puffing would also be non-axisymmetric. Third, Debes et al. (2017) suggest that one of the currently observed gaps in the inner region is capable of interfering with the shadow. While shadows from the puffed edges of a planet gap have been considered (Jang-Condell, 2008, 2009), our simulations do not include a planet at larger radii and so cannot address this.

Our scattered light image in the right panels of Figure 3 suggests an additional option: warping of the outer disc by a misaligned planet produces comparable fluctuations in the azimuthal surface brightness profile. Crucially, the strength of this effect depends on the radius due to the radial dependence of the warp. Figure 5 also demonstrates that a planet massive enough to separate an inner disc (in our simulations, $\gtrsim$2MJ) is also massive enough to warp the outer disc. Thus any planet that can drive an inner, precessing disc must cause warping of the outer disc which would affect the profile of the shadow as a function of the radius. In this scenario, a misaligned planet could self consistently explain two of the key features of the shadow in TW Hya.

4.2 Depth of the shadow

The depth of the shadow (measured as the amplitude of the variation in the surface brightness profile, $A$ in Equation 1) is likely to have contributions from both the shadow cast by the misaligned inner disc and the warp in the outer disc. How these effects sum depends on their relative phase; as the warp in the outer disc precesses much more slowly than the precessing inner disc. Thus there must exist particular arrangements where they add constructively and the shadow appears deeper than it otherwise would (and vice versa). For the parameters chosen in our simulations (in particular, the location of the planet orbit and the aspect ratio), Figure 5 indicates that either mechanism in isolation could produce variations that are consistent in amplitude to TW Hya for planet masses $\sim 6$MJ ($\sim 20\%$, Debes et al., 2017).

4.3 Speed of the shadow

The rate of precession of the inner disc depends strongly on the mass of the planet and its orbital inclination to the precessing disc. In the case that the inner disc is precessing differentially relative to the outer disc, we have additional constraints on the relationship between the planet mass, radial extend of the inner disc and the radial extent of the outer disc.

The frequency of precession of a rigidly precessing circumprimary disc in an inclined binary is given by Larwood et al. (1996),

[TABLE]

Here $x_{\rm in}$ and $x_{\rm out}$ are the inner and outer edges of the disc respectively, $m$ is the planet mass, $a$ is the semi-major axis of the orbit, $G$ is the gravitational constant and $\beta^{\prime}$ is the angle between the disc and the inclined planet. This expression can be used to calculate the precession time-scale, as $t_{\rm prec}=2\pi/\omega_{\mathrm{prec}}$. If we assume $\Sigma\propto R^{-1}$ and a very small inner edge, this reduces to the form used by Debes et al. (2017)111In deriving Equation 3 from Lai (2014), we find a different leading coefficient to that in Equation 2 of Debes et al. (2017) by a factor of $50^{3/2}$. (originally from Lai, 2014),

[TABLE]

where $\mu_{\rm c}$ is the mass ratio of the secondary to the primary, $i_{\rm c}\equiv\beta^{\prime}$ in Equation 2, $r_{\rm disc}\equiv x_{\rm out}$ and $a_{\rm c}$ is the semi-major axis of the secondary. Following Debes et al. (2017), we can solve either of these equations for the position and mass of the planet that will drive precession at the rate observed in TW Hya.

The solution to both Equations 2 and 3 for $t_{\rm prec}=P=15.9$ years is shown in Figure  9. In this representation, a planet must lie in the vicinity of the solid and dashed lines in order to cause precession of a disc that agrees with the rate observed in TW Hya. Here we have assumed the disc is truncated at 0.05 au (although it may be smaller than this, Johnstone et al., 2014), the disc extends to 90% of the planet’s orbit (as in Debes et al., 2017), the star has a mass of $0.7$M*⊙* (Huélamo et al., 2008) and the inclination between the disc and the planet is small. Additionally, we assume a $\Sigma$ profile in Equation 2 that has a taper at the inner edge such that $\Sigma(R_{\rm in})=0$. Both of these relations agree that in order to drive the observed precession rate, such a planet must be massive ($\gg$ the gap clearing mass) and relatively close in. For example, a planet of $12.5$MJ located at $\sim 0.2$ au would drive a disc between 0.05 and 0.2 au to precess with a period of 15.9 years. These planet masses are broadly consistent with previous predictions by Facchini et al. (2014).

Radial velocity (RV) observations of pre-main-sequence stars are challenging, due to the high levels of stellar activity, but if we treat the observations of Setiawan et al. (2008) and Huélamo et al. (2008) as upper limits on the mass of any companion in TW Hya then they yield a combined upper limit to the RV amplitude of 238$\pm$9 m s*-1*; any putative planet in the system must induce a stellar reflex velocity below this limit. Assuming a viewing inclination of 7*∘, we solve for the planet orbit and mass that will cause a radial velocity measurement of 238$\pm$9 m s-1*. The relationship between the planet mass and radius is shown in Figure 9 by the dotted line and represents an upper limit for the mass of the planet — a planet that is more massive than this should have been revealed in existing observations. This comparison suggests that it is possible for a planet that is massive enough to cause a precessing, inner disc to remain undetected by current radial velocity measurements.

4.4 Long term evolution

Our simulations in NDAMN18 are limited in that they do not complete a full precession of the inner disc. We instead consider the viscous and alignment time-scales to infer the long term evolution.

The mode of differential precession invoked to cast the moving shadow requires that the planet and inner disc have some mutual misalignment (Equation 2). This was evident in NDAMN18, where we see differential precession cease after $t\sim 50$ orbits, when the planet and inner disc show the same tilt (Figures 8 and 9 of that paper). When the inner disc and planet are aligned, a shadow will still be cast (the inner disc will still be misaligned with respect to the outer disc) but it will move with the much slower precession rate of the entire disc. The rate at which this alignment occurs between the planet and the inner disc is on the order of the precession time-scale (Scheuer & Feiler, 1996) — this is $\sim$15.9 years in the case of TW Hya. However, this azimuthal feature has been identified in observations spanning about 17 years (Roberge et al., 2005; Debes et al., 2017; Poteet et al., 2018). Given the length of time this shadow should be visible compared to the disc lifetime, it is thus unlikely that the inner disc is driven by a single misaligned planet.

Independent of the mechanism that causes the disc to precess, its lifetime may be limited by how quickly it is accreted onto the star. The accretion rate in TW Hya is estimated to be $\sim 1.8\times 10^{-9}$ M*⊙*/yr (Ingleby et al., 2013). We relate the mass of the disc inside the planet orbit to this accretion rate to find the expected disc life-time, optimistically estimating that the disc could extend to 0.3 au (from Figure 9). Assuming a surface density profile of (Pringle, 1981)

[TABLE]

where $\Sigma_{0}$ is calibrated to the total disc mass, $R_{0}=1$ au and $p=1$, the disc mass from the inner edge of the disc $R_{\rm in}=0.05$ au to the planet orbit at $R_{p}=0.3$ au is

[TABLE]

Assuming a total disc mass (i.e. between 0.1 au and 50 au) of $5\times 10^{-3}$ M*⊙* as in our simulations, this suggests the inner disc life-time of $\sim 6.3\times 10^{3}$ years or less than $\sim 0.1$% of the total disc life-time of $\sim 10^{7}$ years (e.g. Weinberger et al., 2013). This measurement is similar to the estimation using the viscous life-time, expressed as

[TABLE]

where $H/R$ is the aspect ratio, $\alpha$ the Shakura-Sunyaev viscosity parameter and $\Omega_{\rm k}$ the Keplerian angular velocity. Here we set $H/R=0.05(R/R_{0})^{1/2-q}$, with $q=1/4$, $R_{0}=1$ au and $\alpha=10^{-3}$. At $R=0.3$ au the viscous time-scale is then $\sim 2.3\times 10^{4}$ years. It is thus highly unlikely that this disc can be sustained long enough to be observed unless it is being fed by the outer disc, but such replenishment is not expected across the gap that must exist in order for the precession of the inner disc to occur.

We note that the above estimate is directly proportional to the total mass of the disc but this is difficult to measure (Thi et al., 2010; Bergin et al., 2013; Trapman et al., 2017; Teague et al., 2018). For the largest realistic total disc mass of $0.1$ M*⊙* the inner disc life-time increases to $1.2\times 10^{5}$ years or $\sim 1.2$% of the disc life-time. The life-time of this inner, precessing disc is thus much shorter than the estimated disc life-time, even taking into account the uncertainty in the total disc mass.

5 Discussion

Observations of a moving azimuthal brightness pattern in TW Hya have been proposed to be caused by a misaligned, precessing inner disc casting a shadow on the outer disc (Debes et al., 2017; Poteet et al., 2018). Although our results in Section 3 show that the observational signatures of a disc misaligned by a planet may provide a consistent explanation, our discussion in Section 4 shows that this is unlikely to be the complete picture. We thus consider other mechanisms that may result in shadowing from the inner disc and the relevant limitations of our simulations from NDAMN18.

5.1 Alternative mechanisms to misalign the inner disc

While precession of the inner disc due to a single misaligned planet is unlikely, the interaction of multiple bodies may be able to sustain precession. For example, the presence of another massive planet may lead to planet-planet scattering interactions that could continue to perturb the inner planet (Thommes & Lissauer, 2003). Indeed, there is circumstantial evidence of planets outside of 1 au in TW Hya (Calvet et al., 2002; Debes et al., 2013; Andrews et al., 2016; Tsukagoshi et al., 2016; van Boekel et al., 2017; Mentiplay et al., 2018). Additionally, if the outer disc is particularly massive it may affect the precession rate of the inner disc. While our simulations only consider a disc mass of $5\times 10^{-3}$ M*⊙*, the observed disc mass is poorly constrained (Thi et al., 2010; Bergin et al., 2013; Trapman et al., 2017; Teague et al., 2018) and may be massive enough to influence the inner disc precession.

One limitation of our simulations in NDAMN18 is that they do not take into account any magnetic field effects. At the predicted location of the misaligned planet (Figure 9) the stellar magnetic field is likely to interact with the disc. An example of how this may affect the inner edge of the disc is found in AA Tau, which has a stellar magnetic field that is misaligned by about 20*∘* (Bouvier et al., 1999; Donati et al., 2010). Non-ideal MHD simulations by Romanova et al. (2004) demonstrated that the inner edge of a disc warps under the influence of the misaligned magnetic field. This type of warping represents an additional method to cast a shadow from the inner disc to the outer disc, but here the time-scale for the shadow to move would be related to the stellar magnetic field structure. In applying such a model to TW Hya, we note that measurements of the magnetic field geometry suggest it is misaligned by $\lesssim 10^{\circ}$ (Donati et al., 2011) and the inner radius of the disc is 0.05 au (Eisner et al., 2006; Johnstone et al., 2014). Thus if the shadow cast from the misaligned inner disc was due to magnetic fields it would likely have a similar time-scale to the stellar rotation period, which is observed to be days in TW Hya. It is not clear how time-scales of years required by the observed shadow in TW Hya can be related to the much shorter time-scales due to stellar magnetic field effects at the inner edge of the disc.

Our simulations are also limited by their radial extent and duration. In NDAMN18 we demonstrated that the misalignment of the inner disc and the warp structure in the outer disc both evolve dramatically across the 200 planet orbits simulated. Both the variation due to the misaligned inner disc and the warped outer disc decrease across the last 100 orbits of our simulations. We thus caution that the amplitudes presented in Figures 3 and 5 are optimistic, but note that the other features of the shadow (i.e. radial location) do not change appreciably for $t>100$ orbits. The radial profile of the warp under the influence of multiple planets also remains unexplored.

5.2 Alternative mechanisms to cast a shadow

The shadow in TW Hya may alternatively be the cast by the dense tidal stream when the planet is on the near side of the disc. We find hints of this in the earlier stages of our simulations, when the planet is still rapidly accreting. In this case the speed of the shadow may be explained by the Keplerian orbit of a misaligned planet at 5.6 au. Additionally, it implies that the shadow would only be observed for the half of the orbit when the planet is on the near side of the disc. This is consistent with current observations (e.g. Figure 7 of Debes et al. 2017) and suggests that the shadow will not be visible after $\sim$2022.

The radial dependence of the shadow observed in TW Hya may be self-consistently explained if the disc is warped by the presence of a misaligned planet. However, warping of the outer disc could also be caused by a smaller, misaligned planet at larger radii or by a previous interaction with an external perturber. Recent simulations by Cuello et al. (2018) showed that the disc may be warped by such an episode and retain the warp long after the interaction (see also Nixon & Pringle, 2010). For the relatively isolated TW Hya, such an interaction appears unlikely. In any case, we note that the time-scale of the shadow cannot be matched exclusively by a warp in the outer disc.

6 Conclusions

We have examined the signatures in scattered light that are caused in a protoplanetary disc under the influence of a misaligned planet. Using our simulations from NDAMN18, we have conducted radiative transfer calculations using mcfost. We have then examined the azimuthal profile of the surface brightness from both a shadow cast from an inner, misaligned disc and the variation due to a warp in the outer disc. Our results are summarised as:

An inner disc misaligned by more than a degree can cast non-axisymmetric shadows on the outer disc which are observable in IR scattered light (with a shadow depth of $\gtrsim 20\%$ compared to the average flux at that radius). This shadow is strongest closest to the inner edge of the outer disc and decreases in intensity at larger radii. The azimuthal location of this shadow is consistent with a shadow cast from the inner disc. In the case that the disc is created by a misaligned planet, our simulations suggest the planet must be $\gtrsim 6.5$MJ and inclined by more than $\sim 2^{\circ}$ to cast an appreciable shadow. These values are specific to our model, but correspond to a planet with 5.0$\times$ the thermal masses at an inclination of 0.5$\times H/R$. In contrast to previous works, a only small relative misalignment between the inner and outer disc is required to produce shadowing (and this inclination is less than the initial inclination of the planet). 2. 2.

The presence of a warp can additionally cause measurable variations in the azimuthal surface brightness. The amplitude of this variation depends on the relative warp between the star and the disc, and thus can change the appearance of this variation as a function of radius. In our simulations with a misaligned planet, the variation is strongest at large radii. Our simulations suggest that planets that are massive enough to cause a misaligned inner disc also drive a warp in the outer disc. 3. 3.

Our simulations suggest that a misaligned inner disc may be consistent with observations of the moving shadow in TW Hya. Despite this, we argue that it is difficult to sustain this configuration over observationally relevant time-scales, making this an unlikely arrangement for TW Hya.

Although the scattered light images from a misaligned inner disc are consistent with observations, we suggest that a single misaligned planet is unlikely to be the mechanism driving such a disc. Alternative scenarios include a misaligned planet in the outer disc providing the outer disc warp, the inner disc warped by a misaligned stellar magnetic field, Keplerian motion of a misaligned inner planet at $5$ au or precession driven by more than a single planet.

Even the case that the misaligned disc is not driven by a planet, this work demonstrates that small to moderate misalignments in the disc are capable of casting shadows onto the outer disc. In combination with previous studies looking at strongly misaligned discs (Facchini et al., 2018; Zhu, 2018), this confirms that warps and misalignments can be used to inform on the underlying dynamics in protoplanetary discs.

Acknowledgements

The authors warmly thank Daniel Price, Cassandra Hall, Chris Nixon and Jim Pringle for useful discussions. We also thank John Debes and the referee for constructive comments. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 681601). DM is funded by a Research Training Program Stipend from the Australian government. This work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. Figure 1 and the upper panels of Figure 3 were plotted with splash (Price, 2007). All other figures were produced using the community open-source Python package Matplotlib (Hunter, 2007).