Shadows and asymmetries in the T Tauri disk HD 143006: Evidence for a misaligned inner disk

While planet formation is thought to occur early in the history of a protoplanetary disk, the presence of planets embedded in disks, or of other processes driving disk evolution, might be traced from their imprints on the disk structure. We observed the T Tauri star HD 143006, located in the 5-11 Myr-old Upper Sco region, in polarized scattered light with VLT/SPHERE at near-infrared wavelengths, reaching an angular resolution of ~0.037"(~6 au). We obtained two datasets, one with a 145 mas diameter coronagraph, and the other without, enabling us to probe the disk structure down to an angular separation of ~0.06"(~10 au). In our observations, the disk of HD 143006 is clearly resolved up to ~0.5"and shows a clear large-scale asymmetry with the eastern side brighter than the western side. We detect a number of additional features, including two gaps and a ring. The ring shows an overbrightness at a position angle (PA) of ~140 deg, extending over a range in position angle of ~60 deg, and two narrow dark regions. The two narrow dark lanes and the overall large-scale asymmetry are indicative of shadowing effects, likely due to a misaligned inner disk. We demonstrate the remarkable resemblance between the scattered light image of HD 143006 and a model prediction of a warped disk due to an inclined binary companion. The warped disk model, based on the hydrodynamic simulations combined with 3D radiative transfer calculations, reproduces all major morphological features. However, it does not account for the observed overbrightness at PA~140 deg. Shadows have been detected in several protoplanetary disks, suggesting that misalignment in disks is not uncommon. However, the origin of the misalignment is not clear. As-yet-undetected stellar or massive planetary companions could be responsible for them, and naturally account for the presence of depleted inner cavities.


Introduction
High angular resolution observations of protoplanetary disks show a wide diversity of structures on different scales. Sub-millimeter (mm) observations show radial structure, for example, bright and dark rings (e.g., ALMA Partnership 2015; Andrews et al. 2016;Fedele et al. 2017;Dipierro et al. 2018), and azimuthal asymmetries (sometimes called "horseshoe") Based on observations performed with SPHERE/VLT under program ID 097.C-0902(A) and 095.C-0693(A).
where the dust continuum emission is much stronger than in the surrounding background disk (e.g., van der Marel et al. 2013;Casassus et al. 2013). Atacama Large Millimeter Array (ALMA) observations at high angular resolution revealed that most, if not all, of the protoplanetary disks are not smooth, and often show ring-like features Cieza et al. 2017;Fedele et al. 2018), suggesting that these features trace universal processes. One of the most interesting scenarios is that they are due to interactions between the planet-forming disk and embedded proto-planets (e.g., Dipierro et al. 2015;Rosotti et al. 2016).
However, other physical mechanisms, such as ice lines, or dead zones, have also been invoked to explain such observations (e.g., Flock et al. 2015;Zhang et al. 2015;Béthune et al. 2016;Okuzumi et al. 2016;Pinilla et al. 2016Pinilla et al. , 2017. Scattered light images show an even more diverse spectrum of structures, including cavities (e.g., Pinilla et al. 2015), rings and gaps, sometimes colocated with sub-mm counterparts (e.g., van Boekel et al. 2017;Pohl et al. 2017a;Avenhaus et al. 2018), and multiple spirals (e.g., Garufi et al. 2013;Benisty et al. 2015). While the continuum mm emission traces the cold disk midplane, scattered light observations trace the small dust particles in the disk upper layers that are directly irradiated by the young star, and hence also depend on the illumination pattern. Therefore, geometrical variations, such as scale-height perturbations due to temperature variations (e.g., Juhász et al. 2015), or due to a warp, can create strong azimuthal asymmetries in scattered light (e.g., Facchini et al. 2018).
Scattered light images of transition disks (disks with an inner dust depleted cavity) show clear evidence for misaligned inner regions, either with narrow shadow lanes in the outer disk Pinilla et al. 2015;Stolker et al. 2016a;Benisty et al. 2017;Casassus et al. 2018) or with low amplitude azimuthal variations observed over time (Debes et al. 2017). Additional evidence for warps in disks are also found in sub-mm observations, through the kinematics of gas lines (Rosenfeld et al. 2012;Casassus et al. 2015;Brinch et al. 2016;Walsh et al. 2017;Loomis et al. 2017;Boehler et al. 2018). Dipper stars provide another example of potentially warped inner disks (e.g., Cody et al. 2014;Bodman et al. 2017). In that case, the warp is thought to be due to a strong dipolar stellar magnetic field that is misaligned with respect to the disk midplane and forces the innermost disk to tilt out of the plane (e.g., AA Tau; Bouvier et al. 2007). The dipper light curves show clear dimming events usually interpreted as the signature of material from a very inclined disk that repeatedly blocks the line of sight. However, recent imaging of moderately inclined outer disks in some dipper stars suggested a strong misalignment between the inner and outer disk regions (Ansdell et al. 2016;Loomis et al. 2017). Such a misalignment can be very large, up to ∼70 • (e.g., Marino et al. 2015;Benisty et al. 2017;Min et al. 2017), and in some cases the resulting shadows can lead to a significant cooling of the outer disk material (Casassus et al. 2018).
In addition to the effect of a strongly inclined stellar magnetic field, warps can result from the gravitational interaction of a massive companion with the disk. If the orbit of the companion is significantly misaligned with respect to the disk midplane, the disk can break, which leads to a significant misalignment between the inner and outer regions of the disk (Nixon et al. 2012;Dogan et al. 2015;Facchini et al. 2018). Other mechanisms, such as secular interactions and precessional resonances Owen & Lai 2017), can also strongly tilt the inner disk. Depending on the location of the companion, these scenarios can lead to a misaligned circumprimary or circumbinary disk. In the first case, such a companion would naturally create a dust depleted inner cavity, as found in transition disks, and explain most of the properties of some of these disks (e.g., HD 142527; Price et al. 2018b). So far, except in HD 142527 and PDS 70, the putative companion inside the cavity of transition disks remains to be clearly detected (Biller et al. 2012;Keppler et al. 2018).
The focus of this study is the protoplanetary disk around HD 143006 (also, 2MASS J15583692-2257153). HD 143006 is a G7 T Tauri star with the following stellar parameters (updated after the Gaia Data Release 2): T eff = 5880 K, L * = 4.58 L , and M * = 1.5 M (Salyk et al. 2013;Garufi et al. 2018). It is located at a distance of 166 ± 4 pc (Gaia Collaboration 2018) and belongs to the Upper Sco star-forming region, which is rather old (5-11 Myr;Preibisch et al. 2002;Pecaut et al. 2012) compared to the typical timescales for disk evolution. Barenfeld et al. (2016) observed ∼100 disk-host candidates in Upper Sco with ALMA, including HD 143006, and showed that, not only does Upper Sco have a lower fraction of disk-host stars than younger regions (such as Taurus or Lupus), but also that these disks have a lower dust mass to stellar mass ratio. The disk around HD 143006 (hereafter only referred to as HD 143006) is resolved at 0.88 mm with an ∼0.35 × 0.3 beam, and its continuum map shows a centrally depleted large cavity (∼84 au; Pinilla et al. 2018b), and a low-contrast brightness asymmetry with an enhanced emission (by a factor 2) in the south-east (Barenfeld et al. 2016). The innermost region, however, remains dust-and gas-rich as indicated by the high near-infrared (IR) excess (∼21%; Garufi et al. 2018), and by near-IR interferometric observations that resolve hot dust at the sub-au scale (Lazareff et al. 2017). A relatively strong H α emission line translates to a mass accretion rate of ∼2 × 10 −8 M yr −1 (Rigliaco et al. 2015).
In this paper we present the first scattered light observations of HD 143006 obtained with the Very Large Telescope (VLT) Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) instrument. Our observations trace the small (suband micrometer-sized) dust grains, well coupled to the gas, in a tenuous surface layer of the disk and show a number of features that share striking similarities with the predictions of a warped disk model. Our paper is organized as follows: in Sect. 2 we present our observations and the data reduction; in Sect. 3 we describe the scattered light images; in Sect. 4 we present the hydrodynamical simulations and radiative transfer predictions; and in Sect. 5 we discuss our findings.

SPHERE imaging
We obtained observations at the VLT located at Cerro Paranal, Chile, on 2016 June 30, using the SPHERE instrument (Beuzit et al. 2008). SPHERE features an extreme adaptive-optics (AO) system (Fusco et al. 2006;Petit et al. 2014;Sauvage et al. 2014) that feeds three science channels and enables high angular resolution and high-contrast imaging at optical (visible and near-IR) wavelengths. HD 143006 was observed with the polarimetric imaging mode of the InfraRed Dual-band Imager and Spectrograph (IRDIS; Dohlen et al. 2008;Langlois et al. 2014), in the J-band (λ 0 = 1.258, ∆λ = 0.197 µm). IRDIS has a plate scale in the J-band of 12.26 mas per pixel (Maire et al. 2016). In addition, to enhance the detection of outer disk features, we used a 145 mas diameter coronagraphic focal mask (N_ALC_YJ_S, with an inner working angle of 0.08 ; Martinez et al. 2009;Carbillet et al. 2011) for one dataset, but removed it for the second dataset, to enable the observation of inner regions otherwise covered by the coronagraph. HD 143006 was observed for ∼35 min on source with the coronagraph, and approximately 5 min on sources without coronagraph, with a seeing of 0.6 -0.8 and 4-5 ms of coherence time. The analysis of the point spread function (PSF) that is estimated from a noncoronagraphic FLUX (short, non-saturated images of the star outside the masked region) measurement shows that the observations reach a 37 mas × 37 mas resolution and a Strehl ratio of 51%. The five inner pixels of the non-coronagraphic polarimetric image (within a radius of ∼60 mas) are saturated. With polarimetric differential imaging (PDI; e.g., Kuhn et al. 2001;Apai et al. 2004), one measures the linear polarization of the light scattered by dust grains in the disk. This technique enables us to efficiently remove the unpolarized stellar contribution and to image with high contrast the outer disk from which we detect polarized scattered light. The instrument splits the beam into two orthogonal polarization states, and a half-wave plate (HWP) is set to four positions shifted by 22.5 • in order to construct a set of linear Stokes images. We reduce the data according to the double difference method (Kuhn et al. 2001), and derive the Stokes parameters Q and U. Assuming only one scattering event for each photon, the scattered light from a protoplanetary disk, at low inclination angle 1 , is expected to be linearly polarized in the azimuthal direction. We therefore describe the Stokes parameters in polar coordinates (Q φ , U φ ; Schmid et al. 2006;Avenhaus et al. 2014), as with φ the position angle of each pixel (x, y) with respect to the star location. The polarized flux along the azimuthal direction appears as a positive signal in the Q φ image, while radial polarization will lead to negative Q φ . If there is only azimuthal and/or radial polarization, the U φ image contains no disk signal and can be used as an estimate of the residual noise in the Q φ image (Schmid et al. 2006). We correct for the instrumental polarization by minimizing U φ and subtracting scaled versions of the total intensity frame from the Stokes Q and U frames (Canovas et al. 2011). The final images were corrected for the true north (by rotating them by 1.775 • in the counterclockwise direction; Maire et al. 2016). The resulting images are shown in Fig. 1.

ALMA observations
HD 143006 was observed by ALMA in Cycle 2 (2013.1.00395.S) with a synthesized beam of 0.35 × 0.30 , PA = −73 • . For more details on the data reduction and calibration, we refer the reader to Barenfeld et al. (2016), who first presented the data. In this paper, we will only discuss the optically thick 12 CO J = 3-2 observations that trace the surface layers as our scattered light observations. 1 In inclined disks, multiple scattering effects lead to a strong signal in U φ (e.g., T Cha; Pohl et al. 2017b).

Scattered light images
Our scattered light images, with and without coronagraph, are shown in Fig. 1, and present distinct features that are annotated in left panel of Fig. 2. These features are more evident when the image is scaled by r 2 to account for the drop of the stellar illumination with radius. Figure 3 provides radial and azimuthal profiles.
As indicated in the left panel of Fig. 2, we detect the following features, from outside in: (a) An outer disk extending from ∼0.3 to 0.5 (∼50 to 83 au) in scattered light. We refer to this region as Ring #1. It shows a clear azimuthal asymmetry, and is only detected between 0 and 185 • .
(b) A region with less polarized signal than the surrounding disk, between ∼0.24 and 0.3 (∼40 to 50 au), which we call Gap #1. We note that this region is not devoid of scattered light signal at our angular resolution (see Fig. 3, left panel). We measure a ratio of polarized surface brightness, after azimuthal averaging, of ∼70% between radii of 0.3 and 0.4 .
(c) A ring-like feature (Ring #2) between ∼0.11 and 0.18 (∼18-30 au) that presents a strong azimuthal asymmetry. At PAs of ∼190 • and 340 • , two dark regions that we will refer to as shadows can be seen, and at PAs between ∼110 • and 170 • , we observe an overbrightness. This range of PAs is the same as the one over which the outer disk is the brightest. We also note that the peak signal of that ring is at different separations on the west and east side.
Along Ring #2, we measure a ratio of brightness of ∼0.3/1 between PA ∼ 50 • and 140 • , respectively, and a similar brightness along east/west. Finally, we measure a radially averaged polarized surface brightness of ∼7 and ∼12% of the max brightness at the shadows' locations (see Fig. 3).
(d) A dark region, called Gap #2, inside 0.11 (∼18 au). The non-coronagraphic image ( Fig. 1, right) provides a better view of that region and supports the detection of a gap, which we marginally detect right outside the limits set by the saturated pixels (∼0.06 , i.e., ∼10 au). The presence of these saturated inner pixels prevents us from determining the inner edge, if any, of this gap, and whether there are other rings inside 10 au. It indicates the presence of an inner cavity that extends beyond the coronagraph (defined with the vertical dotted line). The error bars are defined as the standard deviation of Q φ in the corresponding bin, divided by the square root of the number of PSF contained in the bin. Right panel: azimuthal profile of the same image after deprojection and averaging across the Ring #2 and Ring #1 widths (0.12 -0.17 , blue curve; 0.4 -0.48 , red curve, respectively). Each curve is normalized to its maximum value, and the red profile is offset by 1.5 for clarity. The error bars are indicated with the shaded region, defined as the standard deviation of Q φ in the corresponding bin, which reflects the large variations of the surface brightness.
also prevents a direct detection of an inner disk inside this area. However, the near-IR excess observed in the spectral energy distribution (SED) of HD 143006 indicates the presence of hot dust at the sublimation radius, which was also spatially resolved with near-IR interferometry (at ∼0.1 au; Lazareff et al. 2017). All the features are very apparent in the polar mapping of the r 2 -scaled image. Figure 2, right panel, presents such an image, obtained after deprojection using i = 17 • and PA = 170 • (see Sect. 4.1). It clearly shows an east/west asymmetry, as well as the bright area between PA ∼ of 110 • -170 • . This overbrightness on Ring #2 does not appear to be co-radial and has some contribution from larger radii. The overbrightness along this range of PA is also evident in the outer disk. The two shadow lanes, the inner gap inside ∼0.11 (Gap #2), and the outer disk (Ring #1 and Gap #1) are also clearly visible.
The presence of an inner gap (Gap #2), depleted in small dust grains, as well as of a second, much shallower gap (Gap #1), clearly appears in the radial profile of the r 2 -scaled Q φ image (Fig. 3, left). This plot is produced after deprojection and azimuthal averaging. The right panel presents the azimuthal profiles of Rings #1 and #2 after deprojection and an average over a width of [0.13 -0.19 ] and [0.32 -0.47 ], respectively. The large azimuthal brightness asymmetries are clearly visible.
The images shown in Figs. 1 and 2, left panel, are not deprojected. The middle panel of Fig. 1, and both panels of Fig. 2 are shown with r 2 scaling to account for the drop-off in stellar illumination and enhance the visibility of faint outer features. This procedure does not take into account the effects of inclination and PA and of the non-planarity of the surface layers that scatter the stellar light (see such a method in Stolker et al. 2016b). Considering the low inclination of the object (as determined in Sect. 4.1), these effects are expected to be small and should not dramatically affect the shape of the features described in this section.

Inner and outer disk misalignments
The characteristics of shadows observed in scattered light images of disks have been well reproduced with a significant misalignment (up to ∼70 • ) between inner and outer disk regions (e.g., Min et al. 2017). To investigate if this could apply to HD 143006, we discuss in this section the inclination and PA values inferred for HD 143006 using various tracers of the inner and outer disk.
Outer disk. Among other transition disks, Pinilla et al. (2018c) modeled the continuum ALMA observations of HD 143006, which trace the outer disk regions, using an asymmetric Gaussian ring model, and found i ∼ 30 • and PA ∼ 148 • . Their continuum dataset (first presented in Barenfeld et al. 2016) resolves a large millimeter dust cavity (∼84 au or ∼0.5 ). It shows a low-contrast asymmetry in the south-east that might affect their values of i and PA.
However, the dust continuum, assuming it is optically thin, traces the disk midplane, while the scattered light signal comes from the surface layers of the disk. We therefore also consider the kinematics of the 12 CO line, an optically thick tracer of the surface layers. To derive the inclination and PA of the outer disk, we consider the Moment 1 map of the 12 CO line published by Barenfeld et al. (2016) and use a simple analytical model of a planar disk in Keplerian rotation around a 1.5 M star. We compute a projected velocity map, convolve it with a two-dimensional (2D) elliptical Gaussian beam (inferred from the ALMA data, 0.35 × 0.30 ) and fit it to the Moment 1 map. We only consider the disk regions where the intensity is above 2σ in the integrated intensity map. We perform our fit with the Markov chain Monte Carlo (MCMC) method using emcee . We note that if the flaring of the disk is large, the thin disk model that we use is not accurate but considering the rather low inclination of the system, it is a reasonable approximation.
The discrepancy between the two estimates based on the ALMA data might be due to the two kinds of observations tracing different regions of the disk, or to the large size of the ALMA beam (∼0.35 , i.e., 60 au) and the complex structure of the disk in the continuum, which makes the continuum-based estimate likely to be less reliable than the one based on CO kinematics.
Inner disk. Inclination and PA measurements of the inner disk are very challenging due to the very high angular resolution that is required to spatially resolve the inner au. HD 143006 was observed with the VLTI H-band instrument PIONIER (Precision Integrated-Optics Near-infrared Imaging ExpeRiment) in the context of a large program focused on Herbig AeBe stars by Lazareff et al. (2017). The H-band visibilities and closure phases trace the thermal emission of the hot dust located in a narrow region at the sublimation radius. In the case of HD 143006, the H-band emission appears to be very compact (∼0.1 au), and is only marginally resolved. Hence, all the analytical models considered for HD 143006 (ellipsoid, ring with and without azimuthal modulation) fit the data equally well and are not well constrained. The inferred inclination and PA values are the following: i = 27 • ± 3 • , PA = 1 • ± 13 • for the ellipsoid model, respectively; i = 23 • ± 5 • , PA = 168 • ± 15 • for the ring model with m = 1 modulation; and, i = 31 • ± 4 • , PA = 148 • ± 21 • for the ring model with m = 2 modulation. Considering the very small extent of the region probed by PIONIER, and the limited angular resolution of the observations, the data can only provide a rough estimate of the inner disk geometry, as indicated by the large error bars and the strong model dependence of the results.
The discrepancy between the various estimates of inclination and PA, as well as the features detected in the scattered light image, suggest that the inner and outer disks might be misaligned. Min et al. (2017), the location and shape of the shadows seen in scattered light depend on the morphology of the inner disk, and on the shape and height of the scattering surface of the outer disk at the cavity edge, hereafter z scat . For a given orientation of the outer disk, and a given z scat , at the location of the shadows, the inclination and PA of the inner disk can be obtained by solving the equations that define the PA of the line connecting the shadows (α) and the offset in declination of this line with respect to the star (η): tan(α) = sin(i 1 )cos(i 2 )sin(PA 1 ) − cos(i 1 )sin(i 2 )sin(PA 2 ) sin(i 1 )cos(i 2 )cos(PA 1 ) − cos(i 1 )sin(i 2 )cos(PA 2 ) η = z scat · cos(i 1 ) cos(i 2 )sin(i 1 )sin(PA 1 ) − cos(i 1 )sin(i 2 )sin(PA 2 ) ,

Misalignments. As shown by
where the indices 1 and 2 refer to the inner and outer disks, respectively. For an inclined disk, at PA = 0 • , the near side of the disks is in the west. These equations lead to two families of solutions, depending on which side of the outer disk is the closest to us. The two possible configurations for the geometry of the inner and outer disks can lead to similar misalignment angles. Because of the complex morphology of the object, and the fact that part of its surface is shadowed, we cannot directly infer from the scattered light observations which side of the outer disk (east or west) is the nearest to us. We therefore present the two families of solutions for the orientation of the inner disk that would reproduce the location of the shadows (in red and blue in Fig. 4). These solutions are shown as a function of z scat /R, with z scat the height of the scattering surface of the outer disk (which differs from the pressure scale height H p by a factor of ∼2-4). The right panel shows the corresponding misalignment β , for the inner disk inclinations and position angles provided in the left and middle panels, respectively. As an example, assuming an outer disk inclination and PA of i = 17 • and PA = 170 • , a misalignment of 30 • (which is the value that we will use in Sect. 4.2) is obtained when the inner disk inclination and PA are ∼13 • and ∼11 • , respectively, and z scat /R is 0.12 at 0.11 (blue curve). If, instead, the near side of the outer disk is opposite (i.e., i = 17 • and PA = 350 • ), such a misalignment is obtained when the inner disk inclination and PA are ∼47 • and ∼356 • , respectively, and z scat /R = 0.17 at the outer disk rim location (red curve). Assuming that the scattering surface corresponds to ∼2-4 H p , we find that the disk aspect ratio H p /R is ∼0.03-0.08 at the outer disk rim.

Hydrodynamic and radiative transfer model
In earlier studies, a parametric approach was used to determine the disk geometry and density structure in the inner and outer disks that would lead to the observed shadowing pattern seen in scattered light observations of protoplanetary disks (e.g., Marino et al. 2015;Benisty et al. 2017). To model HD 143006, we use 3D hydrodynamical simulations, first presented in Facchini et al. (2018). Our observations (Fig. 1) are strikingly similar to their predictions (see their right panel of Fig. 10), in particular regarding the east/west brightness asymmetry.
We provide here a summary of the simulations, and for more details, we refer the reader to Facchini et al. (2018). The 3D simulations have been performed with the Smoothed Particle Hydrodynamics (SPH) code PHANTOM (Price et al. 2018a), using 10 6 particles. We consider a protoplanetary disk and an equal mass binary with a semi-major axis a 0 , and inclined by 60 • with respect to the disk. The disk has an initial surface density scaling with r −1 , with r being the radial coordinate. The temperature profile in the simulation is taken to be vertically isothermal, with the temperature scaling as T ∝ r −1/2 , and an aspect ratio of H p /r = 0.041 at r = 1.7a 0 . After a few binary orbits the circumbinary disk breaks into two separate annuli, driven by the tidal torques generated by the binary on an inclined orbit. Once the inner disk disconnects from the outer disk, it precesses freely around the binary angular momentum vector ). The inclined ring extends from 1.7a 0 to 5a 0 . For this specific setup, the inner and outer disks can show a mutual misalignment between 10 • and 110 • , with the angle varying as the inner disk precesses.
To generate synthetic observables from the SPH hydrodynamic simulations, we assume that the small dust grains, which scatter light efficiently, and the gas, are dynamically coupled, and use the 3D radiative transfer code RADMC-3D 2 . We note that the temperature profiles of the hydrodynamical and radiative transfer simulations are not computed self-consistently. We first scale the hydrodynamic simulations such that the binary orbital separation a 0 is 5.2 au. The disk aspect ratio being H p /r = 0.041 at 8.8 au, with a flaring index of 0.25, implies that at the outer disk location (0.11 , i.e. ∼18.3 au), H p /r = 0.05, consistent with the estimate obtained from the relative misalignment of the disks and the location of the shadows (see Sect. 4.1). After the scaling of the simulations, we interpolate the particle-based density distribution in the SPH simulations to the 3D spherical mesh used in the radiative transfer calculations, using the standard cubic spline kernel. The spherical mesh consists of N r = 220, To compare to our SPHERE images, we post-process a snapshot corresponding to 245 binary orbits, for which the misalignment between the inner and outer disks is ∼30 • (see Fig. C.1). In this snapshot, the outer disk inclination is ∼16 • , and its PA is ∼170 • , while for the inner disk the inclination and PA are ∼16 • and 14 • , respectively. As explained in Sect. 4.1, since we do not know which side of the outer disk is closer to us, we consider a second solution with ∼14 • and ∼350 • as the inclination and PA of the outer disk, respectively, and ∼44 • and ∼355 • for the inclination and PA of the inner disk. In the first solution, while the inclination with respect to the line of sight is the same, the orientation of the two disks is almost opposite: the east side of the outer disk is closer to us, while the near side of the inner disk is in the west. In the second solution, the near side of both the inner and outer disks is the west side, but the inclinations differ by 30 • . These values are close to the estimates derived in Sect. 4.1, based on the location of the shadows. We note however, that these values are model-dependent as they depend on the density considered in the inner disk and within the gap separating it from the outer disk. We therefore expect that other types of models (e.g. with a small circumprimary disk instead of a circumbinary disk) would provide slightly different values as long as the misalignment is moderate (∼20-30 • ).
We use the stellar parameters mentioned in Sect. 1. Barenfeld et al. (2016) derived a dust mass of ∼24.3 Earth masses for the disk, by converting the sub-millimeter continuum flux using d = 145 pc. We scale this value using the new Gaia distance (d = 166 pc), and consider 9.5 × 10 −3 M . To be consistent with the hydrodynamical simulations, we used two identical stars that reproduce the total stellar luminosity of HD 143006. The dust opacity and scattering matrix elements were calculated from A171, page 6 of 14 the optical constants of astronomical silicates (Weingartner & Draine 2001) for a grain size distribution of n(a) ∝ a −3.5 for grain sizes between a = 0.1 µm and a = 1 mm. We use 10 8 photon packages to calculate the dust temperature in a thermal Monte Carlo simulation as well as to calculate scattered light images in the J-band (1.2 µm). Once the synthetic images are computed, we convolve the scattered light predictions with a FLUX image from the dataset, to reduce their resolution to that of the observations. The model predictions are shown in Fig. 5. The inner disk is inclined and casts a shadow onto the outer disk at two points (similar to the cases of HD 142527 and HD 100453), but as the inclination is only moderate, the shadow also darkens half of the outer disk. Because we do not know which side of the disk is closer to us, we provide two solutions in Fig. 5 (models A  and B). The upper panels show a sketch of the 3D structure of the disk, with a color coding that indicates which part of the disk is above or below the plane perpendicular to our line of sight. Our model reproduces most of the features observed in the scattered light observations: a clear east/west brightness asymmetry, two narrow shadows, and two bright arcs tracing Ring #2. We note that the circumbinary disk clearly appears in the synthetic image, while it is not detected in the observations up 60 mas (i.e., ∼10 au). This supports the presence of a small, misaligned circumprimary disk rather than a circumbinary disk. This is discussed further in Sect. 5. We note that the model, and the consequent shadowing due to a misaligned inner disk, cannot reproduce the bright region observed along PA ∼ 110-170 • , in particular along Ring #2, nor the outermost gap (Gap #1).

Discussion
Warps have been inferred in many protoplanetary disks, with various observational tracers. Shadows in scattered light appear as steady low brightness regions (Stolker et al. 2016a;Benisty et al. 2017;Casassus et al. 2018), for which a moderate to large misalignment between the inner and outer disks up to ∼70 • was suggested from radiative transfer modeling. In one of them, HD 142527, a stellar companion on an eccentric orbit is thought to be responsible for the misalignment (Price et al. 2018b). In other objects, such as SAO 206462 and RXJ1604.3-2130A, the shadows appear very variable in amplitude, width, and location (Stolker et al. 2016b;Pinilla et al. 2018a). For example, in RXJ1604.3-2130A, an object known to be an aperiodic dipper (Ansdell et al. 2016), we find that the timescale for the variations is shorter than a day, indicating a very complex and dynamic inner disk (Pinilla et al. 2018a). In AA Tau, a strong, inclined magnetic field induces a warp at the disk inner edge that periodically rotates with the stellar period (Bouvier et al. 2007). However, recent observations by Loomis et al. (2017) indicate that the inner disk is also perturbed, with an additional warp and evidence for a radial inflow, possibly due to gap-crossing streamers. Another case of a perturbed inner disk with a warp is V354 Mon, which presents a low gas-to-dust ratio in the inner disk, and dimming events that could be due to small dust particles that result from the fragmentation of larger particles that drift from the outer disk (Schneider et al. 2018). In all these objects, the presence of a companion in the stellar or sub-stellar mass regime, at a separation of a couple of tens of au, could explain some characteristics of the observations.

Origin of the warp
In this paper, we consider a model of a broken and misaligned circumbinary disk due to the gravitational influence of a stellar companion on an inclined orbit. It successfully reproduces the general characteristics of the scattered light image, supporting the idea that the disk of HD 143006 hosts a warp inside the observed cavity. We stress that our observational predictions hold for any warped disk with a moderate misalignment, independent of what is causing the torque that leads to the misalignment. In particular, as HD 143006 shows a near-IR excess, indicating the presence of hot dust grains very close to the star, it is likely that it hosts a circumprimary disk that is possibly tilted, rather than a circumbinary disk, as in the HD 142527 system. However, since we do not have any direct image of the innermost regions, we cannot directly determine the outer extent of the inner disk, the location of the warp, or the location of the putative companion inducing it.
Upper limits on the binarity. Kraus et al. (2008) led a survey on the stellar binarity of 82 young stars from Upper Sco, using non-redundant aperture masking interferometry, a technique that allows us to search for companions at the diffraction limits. Combining their results with the ones from the literature, they report a frequency of binary companions of ∼33 +5 −4 % at separations of 6-435 au, using d = 145 pc as the distance of Upper Sco. For HD 143006, they estimate lower limits on the K-band contrast of ∼3.5 for separations within 20-40 mas, and of ∼5.1 within 40-80 mas, that is, detection limits of companions with K-band apparent magnitude of 10.6 and 12.2, respectively. Using the BT-SETTL models (Baraffe et al. 2015), considering d = 166 pc, and assuming an age of 10 Myrs, these detection limits translate into companion masses of 0.45 and 0.16 M , corresponding to a mass ratio of q = 0.3 and 0.1, respectively. Using archival near-IR (H-band) interferometric (VLTI/ PIONIER) observations of HD 143006, we performed an analysis of the closure phase data, assuming that the inner disk is point-symmetric and that any departure from point-symmetry detected in the closure phase would be due to the presence of a binary companion. We used the Companion Analysis and Non-Detection in Interferometric Data algorithm (CANDID; Gallenne et al. 2015), which allows an estimate of the 3-σ detection limit for any companion at varying separation from the primary star. Similar to Kraus et al. (2008), we find that a companion with a H-band contrast lower than 3.7 magnitude would have been detected at more than 3-σ within 150 mas (i.e., 25 au, bandwidth smearing limitation), which translates into a mass ratio of q = 0.2 using the parameters mentioned above.
A warp induced by a companion. Considering these detection limits, the presence of an equal mass binary can be excluded, but a low mass stellar object, or a massive planet, could still be present in the cavity of HD 143006. For example, assuming that HD 143006 has a stellar mass of 1.5 M , a 10 M Jup planet (q = 0.006) would not be detected by the interferometric observations. Such a massive object could possibly be responsible for the misalignment of a circumprimary disk. If the inner disk angular momentum is lower than that of an inclined planet, the latter can tilt the disk inside its orbital radius (e.g., Matsakos & Königl 2017). If the planet is massive enough to carve a gap, the inner and outer disks separate and can be both tilted with respect to the midplane (e.g., Bitsch et al. 2013;Nealon et al. 2018). For a 6-M Jup planet with initial inclinations ranging from 10 to 80 • , the misalignment can be up to 15 • (Xiang-Gruess & Papaloizou 2013). Owen & Lai (2017) recently investigated a mechanism that can lead to large misalignments with a companion-star system that is originally coplanar. They find that for a stellar/companion mass ratio of 0.01-0.1 and separation of 10-100 au, secular precession resonances can generate large misalignments between the inner (circumprimary) and outer (circumbinary) disks. The resonance between the inner disk and the companion can lead to a wide range of misalignment angles between the inner and outer disk, some as large as the ones observed for HD 142527 and HD 100453, while the companion remains in the plane of the outer disk. Finally, a massive planet could become misaligned via the Kozai effect , when an additional external companion is present. However, this scenario does not seem plausible, since there is no hint of a binary companion orbiting outside the circumstellar disk.
In general, it is unclear if the disk would be broken and misaligned at a specific stage of its evolution. In the conditions of our hydrodynamical simulations, disk breaking is favored for a low disk aspect ratio, therefore a misaligned companion could more easily break the disk when the star gets older and colder (which in turn implies a colder disk and lower aspect ratio). This would support the findings of Garufi et al. (2018), that shadows observed in scattered light, smoking guns of misalignments, are primarily found around rather old objects. Interestingly, these are often also the objects with higher photospheric metallicity (Banzatti et al. 2018), showing spiral arms, and with the highest near-IR excess, like in the case of HD 143006 (∼21%), which indicates a strongly inflated inner disk with micrometer-sized dust grains reprocessing stellar light at high altitude above the disk midplane. This might be a result of the interaction of the inner disk with an inclined companion. If such a companion creates a gap or cavity, which is depleted with time, the inner disk will lose mass without efficient replenishment from the outer disk, and at some point the angular momentum of the companion will become larger than the one from the inner disk, favoring a late misalignment between the inner and outer disks.
The misalignment of the companion(s) could also occur very early in a disk's history and be an imprint from the early stages of star formation. Misaligned planets could form by fragmentation in the accretion phase of the protostellar envelope (Terquem & Papaloizou 2002) and while most are ejected from the system, some could remain in an inclined orbit around the star. In this framework, Teyssandier et al. (2013) find that massive planets (approximately the mass of Jupiter) would be circularized and align in the disk midplane but that Neptune-mass planets would remain on inclined orbits over the disk lifetime. They note, however, that if the disk mass steeply decreases with time, more massive objects could remain misaligned. Interestingly, a Neptune-mass planet would allow a continuous replenishment of small dust in the inner disk from the outer disk (Pinilla et al. 2015;Rosotti et al. 2016), which is needed in the case of HD 143006.
Other scenarios. Apart from being induced by a misaligned planet, inner disk warps can be due to a strong misaligned dipolar magnetic field as in the well-documented case of AA Tau. Lavail et al. (2017) find a magnetic field of 1.4 kG in HD 143006 (called V1149 Sco in their paper), but the data are not sufficient to obtain a topology of the magnetic field. If the magnetic field is strongly inclined, and warps the inner disk edge, such a warp would rotate with the stellar period (approximately a few days), and lead to fast changing shadows. This scenario can easily be tested with multi-epoch observations. Unfortunately, our two J-band epochs (with and without coronagraph) were obtained on the same observing night but a low signalto-noise ratio optical image, obtained a year before with the Zurich IMaging POLarimeter (ZIMPOL) instrument is shown in Fig. D.1. Although the low quality of the image prevents a detailed analysis of the disk image, we find that the brightest features (east/west sides of Ring #2) have not moved significantly, and hence that it is likely that the shadows do not move at the stellar period.
Another scenario that would lead to misalignments of inner and outer disk regions, is if they are primordial ones, due to a late accretion of material that has an angular momentum misaligned with that of the star (Bate et al. 2010;Dullemond et al. 2018). Recent hydrodynamical simulations carried out by Bate (2018) lead to such a case, with a circumbinary disk whose inner and outer disk planes differ, but this is not a common outcome. Instead, many of the simulations with multiple stars lead to circumstellar and/or circumbinary disks that can be misaligned with each other (Bate 2018).

Rings, gaps, and brightness asymmetries
Our scattered light image shows two rings and two gaps. The outer ring (Ring #1) might be tracing the outer disk up to what our sensitivity allows and might therefore not be tracing any ring-like perturbation (in density or scale height) of the disk. Gap #1 shows significant scattered light signal, and does not appear empty of small grains, in particular in the region located between PAs 110 • and 170 • , as seen in Figs. 2 and 3, right panels, in which the disk appears almost continuous. Gap #1 could be tracing a marginal depletion in small dust, in a gap opened by a planet (e.g., Dong & Fung 2017), while Ring #2 would be its inner edge. It is also possible that Gap #1 is due to selfshadowing by the inner edge of the outer disk (Ring #2). In that case, the denomination of Ring #1 is artificial, in the sense that it would trace the outer disk illumination beyond the shadow of Ring #2 rather than a depleted region. By contrast, Gap #2 appears to be quite depleted in dust, so Ring #2 could be directly irradiated by the star and puffed up. This would naturally lead to Ring #2 casting a shadow on the outer disk.
Modeling of the visibility data of ALMA observations reveals that there is a large cavity in the millimeter emission, surrounded by an asymmetric ring peaking at 84 au (0.5 ; Pinilla et al. 2018b). This location coincides with the outer edge of Ring #1 that extends from ∼0.3 to 0.5 and peaks approximately at ∼0.45 . If this outer ring is a real density enhancement (instead of being only an effect of self-shadowing), Ring #1 and the asymmetric ALMA ring may have a common origin due to dust trapping in pressure bumps possibly induced by a planetary companion. On the other hand, if Ring #1 is not a density enhancement and if it only results from self-shadowing, the very large radial segregation between the inner ring (Ring #2) and the ring observed with ALMA cannot be explained by a single giant planet, and instead, a stellar-mass companion would be required.
The near-IR excess indicates the existence of an optically thick dusty belt located close to the dust sublimation radius within the first few astronomical units. At these distances from the star, the gas density is expected to be high enough that the millimeter dust particles too should be coupled to the gas. Therefore, with very high angular resolution (∼0.01-0.02 ) observations with ALMA, it may be possible to detect a misaligned inner disk in the millimeter emission as well, and confirm our findings based on the shadows observed in the scattered light image.
Even though a misaligned disk model is successful in reproducing the east/west asymmetry, it does not account for the narrow overbrightness between PA 110 • and 170 • (see Fig. 5). Azimuthal asymmetries in scattered light can be due to the scattering angle and polarization efficiency of dust grains located on the disk surface. Large grains (≥5 µm) are efficient forward scatterers, often leading to one side of the disk (the near side) being brighter than the other. Assuming that the PA of the outer disk is ∼170 • , we would expect to see the overall east or west side brighter than the other, and not over such an azimuthally narrow region. The polarization efficiency being maximum for 90 • scattering angles can in turn lead to bright lobes along the major axis of the disk for inclined disks (see, e.g., the radiative transfer model in Benisty et al. 2017). The disk of HD 143006 being only very slightly inclined, the latter possibility cannot account for the strong brightness asymmetry. It is possible that the disk presents a local overdensity, for example, due to the formation of a vortex, an eccentric ring, or a spiral arm at the outer edge of the cavity. If the overbrightness traces the tip of a spiral arm, we would expect it to not be co-radial, as observed for Ring #2. Interestingly, the continuum ALMA data (Barenfeld et al. 2016;Pinilla et al. 2018b) show a low contrast asymmetry (contrast less than a factor of two) along the same angle as the one observed in the scattered light data, supporting a density enhancement. However, it is located close to the outer ring (Ring #1), at large radii (in the ALMA image, peaking at 0.5 or 84 au), and does not coincide with the inner ring (Ring #2). From the current ALMA observations, the asymmetry is not resolved and its exact morphology is still an open question. As shown in the polar map, not only Ring #2 shows the overbrightness, but the outer disk does too. This suggests that along this range in PAs, the disk is more strongly illuminated than the rest of the disk. An equally possible scenario is that most of the disk (between PAs of 0 and 110 • ) lies in a partial shadow caused by the inner disk, while the bright region (PAs 110-170 • ) is unshadowed. This additional shadow could be due to asymmetric features from the inner disk that have a moderate radial optical depth, likely from tenuous surface layers. In either scenario, the physical cause of this phenomenon is not clear.

Conclusions
In this paper, we present the first scattered light observations of the circumstellar disk around the T Tauri star HD 143006. Our observations reveal two rings and two gaps, a strong east/west brightness asymmetry, an overbrightness along a narrow range of position angles (110-170 • ), and two dark narrow lanes. Such azimuthal brightness variations are indicative of shadowing effects, in particular due to a misaligned inner disk.
We analyze the kinematics of the 12 CO line, as observed with ALMA, using an analytical model of a razor-thin disk in Keplerian rotation and derive an inclination of ∼17 • and a position angle of 170 • for the outer disk. Combined with the constraints derived from inner disk observations with near-IR interferometry, this suggests that the inner and outer disk regions are moderately misaligned. We provide two possible solutions for the inner disk orientation to reproduce the location of the shadows, depending on the side of the outer disk that is nearest to us.
Our scattered light image shares a striking resemblance with synthetic predictions based on hydrodynamical simulations of a protoplanetary disk warped by an inclined equal mass binary . In these simulations, the circumbinary disk breaks into two distinct annuli (inner and outer disk, both circumbinary), and the inner disk precesses freely around the binary angular momentum vector. To compare with our observations, we post-process a snapshot of these simulations for which the relative misalignment between the two annuli is ∼30 • , considering the stellar parameters of HD 143006. This model reproduces the east/west asymmetry, but does not account for the additional overbrightness along a narrow range of position angles which might be due to an overdensity not included in our model.
Although our model uses an equal mass binary, which can be ruled out by current detection limits, we stress that our predictions hold for any warped disk, independently of the cause of the misalignment. In particular, a massive planet (e.g., with a mass ratio of 0.01-0.1) might break the disk, or alternatively, an inclined magnetic field could misalign the innermost disk edge as in AA Tau, although our marginal evidence that the shadows have not rotated within a year does not support the latter scenario.
Further observations of this system, in particular with ALMA at high resolution, will allow us to constrain the orientation of the inner disk in large (mm) grains, and could confirm or contradict that the features observed in our scattered light images are due to a misaligned inner disk. Comparison between ALMA and SPHERE observations at similar angular resolution will also allow us to constrain whether the disk shows evidence for cooler regions in the mm, due to the shadows, as in DoAr 44 (Casassus et al. 2018).
Scattered light shadows have now been found in a handful of objects, often in transition disks with large cavities that could host a high planetary mass or a low stellar mass companion (as in HD 142527; Biller et al. 2012), which would still be below the current detection limits. It is therefore possible that all transition disks host stellar or planetary-mass companions with a mass ratio of ∼0.01, some with inclined orbits. It is unclear however, if the misalignments that the shadows observed in scattered light trace could be the origin of the relative inclinations between the stellar rotation axis and orbit orientation found in many exoplanetary systems.