<figure><img src="image_1.png" /></figure>

# Variable Dynamics in the Inner Disk of HD 135344B Revealed with Multi-epoch Scattered Light Imaging*

Tomas Stolker<fnref n="1" />, Mike Sitko<fnref n="2" />, Bernard Lazareff<fnref n="3" />, Myriam Benisty<fnref n="3" />,<fnref n="4" />, Carsten Dominik<fnref n="1" />, Rens Waters<fnref n="1" />,<fnref n="5" />, Michiel Min<fnref n="1" />,<fnref n="5" />, Sebastian Perez<fnref n="6" />,<fnref n="7" />, Julien Milli<fnref n="8" />, Antonio Garufi<fnref n="9" />, Jozua de Boer<fnref n="10" />, Christian Ginski<fnref n="10" />, Stefan Kraus<fnref n="11" />, Jean-Philippe Berger<fnref n="3" />, and Henning Avenhaus<fnref n="6" />,<fnref n="7" />,<fnref n="12" />

*Received 2017 May 24; revised 2017 August 15; accepted 2017 August 22; published 2017 November 9*

## Abstract

We present multi-epoch Very Large Telescope/Spectro-Polarimetric High-contrast Exoplanet REsearch (VLT/SPHERE) observations of the protoplanetary disk around HD 135344B (SAO 206462). The J-band scattered light imagery reveal, with high spatial resolution (~41 mas, 6.4 au), the disk surface beyond ~20 au. Temporal variations are identified in the azimuthal brightness distributions of all epochs, presumably related to the asymmetrically shading dust distribution in the inner disk. These shadows manifest themselves as narrow lanes, cast by localized density enhancements, and broader features which possibly trace the larger scale dynamics of the inner disk. We acquired visible and near-infrared photometry which shows variations up to 10% in the JHK bands, possibly correlated with the presence of the shadows. Analysis of archival Very Large Telescope Interferometer/Precision Integrated-Optics Near-infrared Imaging ExperIMENT (VLTI/ PIONIER) H-band visibilities constrain the orientation of the inner disk to $i = 18^{\circ}2_{-4.1}^{+3.4}$ and PA = 57$^{\circ}3 \pm 5$ $^{\circ}$7, consistent with an alignment with the outer disk or a minor disk warp of several degrees. The latter scenario could explain the broad, quasi-stationary shadowing in north-northwest direction in case the inclination of the outer disk is slightly larger. The correlation between the shadowing and the near-infrared excess is quantified with a grid of radiative transfer models. The variability of the scattered light contrast requires extended variations in the inner disk atmosphere ($H/r \lesssim 0.2$). Possible mechanisms that may cause asymmetric variations in the optical depth ($\Delta\tau \lesssim 1$) through the atmosphere of the inner disk include turbulent fluctuations, planetesimal collisions, or a dusty disk wind, possibly enhanced by a minor disk warp. A fine temporal sampling is required to follow day-to-day changes of the shadow patterns which may be a face-on variant of the UX Orionis phenomenon.

**Key words:** protoplanetary disks – radiative transfer – scattering – stars: individual (HD 135344B) – techniques: high angular resolution – techniques: interferometric

## 1. Introduction

Spatially resolved observations provide detailed insight into the physical and chemical processes occurring in protoplanetary disks. The disk around the intermediate-mass star HD 135344B (SAO 206462) is a suitable target to be observed with high-resolution due to its proximity ($156 \pm 11$ pc; Gaia Collaboration et al. 2016), spatial extent (~1"15 in the scattered light; Grady et al. 2009), low inclination (16°; van der Marel et al. 2016), and brightness from visible to millimeter wavelengths (e.g., Carmona et al. 2014). The disk is classified as a transition disk (e.g., Espaillat et al. 2014) with a large dust cavity resolved at continuum (sub)millimeter wavelengths ($R_{cav} = 51$ au at 156 pc; Andrews et al. 2011). In scattered light, two symmetric spiral arms

have been detected (Muto et al. 2012), which might indicate the presence of a massive gas giant in the outer disk (Dong et al. 2015; Fung & Dong 2015; Dong & Fung 2017). So far, only upper limits on planet masses have been derived from direct imaging observations ($3 M_{Jup}$ at 0"7, assuming hot-start evolutionary models; Maire et al. 2017). The cavity radius in scattered light ($R_{cav} = 27$ au at 156 pc; Stolker et al. 2016), tracing micron-sized grains in the disk surface, is located inward with respect to the large grains in the midplane, which can be explained by planet-induced dust filtration (Garufi et al. 2013). The scattered light cavity coincides with a region in which the surface density of CO gas is significantly reduced (van der Marel et al. 2016).

Pre-main-sequence stars are commonly variable at optical and near-infrared wavelengths on various timescales, for example due to rotational modulation by stellar spots, variable accretion, dust obscuration, and structural changes in the inner disk (e.g., Eiroa et al. 2002). Variability also occurs at mid-infrared wavelengths, for example, *Spitzer/IRS* spectra show a

<footnote n="1">¹ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands; T.Stolker@uva.nl</footnote>

<footnote n="2">² Department of Physics, University of Cincinnati, Cincinnati OH 45221, USA</footnote>

<footnote n="3">³ Université Grenoble Alpes, IPAG, F-38000 Grenoble, France CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="4">⁴ Unidad Mixta Internacional Franco-Chilena de Astronoma (CNRS UMI 3386), Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="5">⁵ SRON Netherlands Institute for Space Research, Sorbonnelaan 2, 3584 CA Utrecht, The Netherlands</footnote>

<footnote n="6">⁶ Departamento de Astronomía, Universidad de Chile, Casilla 36-D, Santiago, Chile</footnote>

<footnote n="7">⁷ Millennium Nucleus “Protoplanetary Disks,” Chile</footnote>

<footnote n="8">⁸ ESO, Alonso de Córdova 3107, Vitacura, Casilla 19001, Santiago de Chile, Chile</footnote>

<footnote n="9">⁹ Universidad Autónoma de Madrid, Dpto. Física Teórica, Módulo 15, Facultad de Ciencias, Campus de Cantoblanco, E-28049 Madrid, Spain</footnote>

<footnote n="10">¹⁰ Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="11">¹¹ University of Exeter, Astrophysics Group, Stocker Road, Exeter EX4 4QL, UK</footnote>

<footnote n="12">¹² Institute for Astronomy, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland</footnote>

<footnote>* Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 087.C-0702(A.B), 087.C-0458(B.C), 087.C-0703(B), 088.C-0670(B), 088.D-0185(A), 088.C-0763(D), 089.C-0211(A), 091.C-0570(A), 095.C-0273(A), 097.C-0885(A), 097.C-0702(A), and 297.C-5023(A).</footnote>
typical anti-correlation between the amplitude of the near- and mid-infrared emission, indicating changes in the height of the inner disk at sub-au distance and consequent shadowing of the disk further outward (Espaillat et al. 2011). The spectral energy distribution (SED) of HD 135344B contains a large near-infrared excess ($F_{NIR}/F_{*} = 0.27$; Garufi et al. 2017) due to the presence of hot dust in the innermost disk region (Brown et al. 2007). It was noted by Garufi et al. (2017) that HD 135344B belongs to a sub-category of group I protoplanetary disks with both a large near-infrared excess and spiral arms, as well as shadows in most cases. The near-infrared continuum emission is variable up to 20%–30%, while the 10 µm flux exhibits fluctuations of 60% (Grady et al. 2009; Sitko et al. 2012). Furthermore, Grady et al. (2009) observed an anti-correlation between the strength of the *J*- and *L'/M'-band* fluxes that was linked to geometrical changes of the inner disk.

Multi-wavelength polarimetric differential imaging observations by Stolker et al. (2016) revealed three shadow lanes in the *J*-band and a broader shaded region bound by two of the shadow lanes. A comparison with optical images from a month earlier showed that the southern *J*-band shadow lane was not present in the *RI* bands, pointing toward a transient or variable origin. Those shadows are presumably cast by dust in the inner disk, which is asymmetrically perturbed and/or misaligned with respect to the outer disk (Stolker et al. 2016). Similarly, Wisniewski et al. (2008) found that the scattered light flux from the protoplanetary disk around HD 163296 showed variations between different imagery epochs obtained with the *Hubble Space Telescope* (HST). This might indicate a time-variable shadowing of the outer disk by scale height variations of the inner disk wall (Wisniewski et al. 2008), in line with the photometric variability in the near-infrared due to structural disk changes near the dust sublimation zone (Sitko et al. 2008). Furthermore, Ellerbroek et al. (2014) reported enhanced extinction in the optical for HD 163296, lasting from a few days up to a year, which was interpreted as caused by a dusty disk wind.

In this paper, we present multi-epoch, polarized scattered light imagery of the protoplanetary disk around HD 135344B that was obtained with the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE; Beuzit et al. 2008) instrument. We aim to detect and characterize brightness variations caused by shading dust in the inner disk. The shadow patterns and their variability allow us to probe the physical processes occurring in the innermost disk region, which are not directly accessible by high-contrast imaging instruments. The scattered light images are complemented with multi-epoch visible and near-infrared photometry that we aim to link to the scattered light variations, and near-infrared interferometry, allowing us to place constraints on the orientation of the inner disk. Furthermore, we will use radiative transfer simulations to quantify the correlation between the scattered light contrast and near-infrared excess in order to estimate the extent of the inner disk variations.

## 2. Observations and Data Reduction

### 2.1. SPHERE/IRDIS Dual-polarimetric Imaging (DPI)

Imaging polarimetry data sets were obtained on 2016 May 04, May 12, June 22, and June 30 with the near-infrared imager (IRDIS; Dohlen et al. 2008) of SPHERE at the European Southern Observatory's Very Large Telescope (VLT). Observations were carried out with the broadband *J*-filter (BB_J, 1.245 µm) in DPI (Langlois et al. 2014) mode. The pixel scale

of the detector is 12.26 mas pix⁻¹ (Maire et al. 2016). An apodized Lyot coronagraph was employed (N_ALC_YJH_S, 185 mas mask diameter), allowing for an integration time of 32 s. The four standard half-wave plate orientations were cycled with two or four subsequent integrations per half-wave plate orientation. The extreme adaptive optics system (SAXO; Fusco et al. 2006) provided a typical Strehl ratio of ~75% in the *H* band (see Table 1).

Seeing conditions were mostly good (<0.′′7), except on 2016 May 12 when the observations were executed with an average seeing of 2.′′1 and the presence of strong winds. Nonetheless, the AO loop remained closed for a total integration time of 17 minutes. The integration time was 34 minutes or longer for the other observations. The point-spread function (PSF) of the sequence on 2016 May 04 was partly affected by low-order aberrations caused by low wind speeds (1.4 m s⁻¹ on average). The PSF quality was evaluated from the non-coronagraphic images recorded by the differential tip-tilt sensor after which 13 polarimetric cycles were removed, leaving a total of 11 cycles (47 minutes).

Standard calibration procedures were applied with the SPHERE data reduction and handling (DRH) pipeline (v0.18.0; Pavlov et al. 2008) which included sky subtraction, flat field correction, and bad pixel interpolation. The frames with horizontally and vertically polarized flux were separated and subsequently processed by a custom pipeline for DPI data. We obtained coronagraphic images with four symmetric satellite spots, induced by a periodic modulation applied to the deformable mirror, before and after the science sequence. These frames were used to determine the position of the star behind the coronagraph mask. To center the coronagraphic DPI data, we interpolated linearly between the start and end position of the star. Stokes *Q* and *U* images were obtained with the double-difference method (Hinkley et al. 2009) and subsequently collapsed with a mean stacking.

To correct for instrumental polarization, we used the method described by Canovas et al. (2011) which assumes that the central star is unpolarized. The $U_φ$ signal, which provides an estimate of the noise level in the single-scattering limit, was minimized by stepwise changing the inner and outer radius of the annulus used to measure the (assumed to be unpolarized) signal close to the star. The azimuthal counterparts of the Stokes *Q* and *U* images, $Q_φ$ and $U_φ$, were calculated with an additional minimization applied on the $U_φ$ image by correcting for a minor rotational offset of the half-wave plate (Avenhaus et al. 2014) for which the optimized values ranged from −2.°0 to −1.°3. We note that the procedure of minimizing the $U_φ$ signal is not strictly valid because part of the scattered light flux from the disk will be present in the $U_φ$ image due to multiple scattering (Canovas et al. 2015). The effect will be small when the disk inclination is low; however, the high signal-to-noise ratio (S/N) of the disk detection around HD 135344B might reveal a real signal in the $U_φ$ image (Stolker et al. 2016). Finally, the images were rotated by −1.°8 toward the true north orientation (Maire et al. 2016).

Flux frames were obtained at the start of each sequence by shifting the star away from the coronagraph, with a shorter integration of 0.87 s and an additional neutral density filter (ND1.0) in the optical path to avoid saturation. The flux frames are used to determine the angular resolution of the images (see Table 1), as well as the scattered light contrast of the disk (see Section 3.1). The data reduction procedure
<table><caption>**Table 1**
SPHERE/IRDIS Observations</caption><thead><tr><th>UT Date</th><th>Integration<sup>a</sup><br>(minutes)</th><th>Airmass</th><th>Seeing<sup>b</sup><br>(arcsec)</th><th>Wind Speed<sup>c</sup><br>(m s<sup>-1</sup>)</th><th>Coherence Time<sup>d</sup><br>(ms)</th><th>Strehl Ratio<sup>e</sup><br>(%)</th><th>PSF FWHM<sup>f</sup><br>(mas)</th></tr></thead><tbody><tr><td>2015 May 03<sup>g</sup></td><td>76.8</td><td>1.11–1.37</td><td>0.69(0.07)</td><td>6.1(0.3)</td><td>2.4(0.4)</td><td>73(4)</td><td>41.3 × 38.4</td></tr><tr><td>2016 May 04</td><td>102.4</td><td>1.06–1.38</td><td>0.52(0.16)</td><td>1.4(0.7)</td><td>13.1(3.4)</td><td>78(2)</td><td>42.0 × 41.9</td></tr><tr><td>2016 May 12</td><td>17.1</td><td>1.14–1.18</td><td>2.14(0.30)</td><td>14.1(0.8)</td><td>3.0(0.6)</td><td>60(7)</td><td>43.1 × 44.8</td></tr><tr><td>2016 June 22</td><td>34.1</td><td>1.09–1.15</td><td>0.63(0.07)</td><td>8.1(0.3)</td><td>8.3(1.4)</td><td>72(2)</td><td>40.1 × 43.6</td></tr><tr><td>2016 June 30</td><td>34.1</td><td>1.02–1.03</td><td>0.37(0.05)</td><td>8.0(0.5)</td><td>9.5(1.5)</td><td>79(1)</td><td>40.2 × 38.1</td></tr></tbody></table>

**Notes.** Values in parentheses provide the standard deviation of the average measured value.

<sup>a</sup> Multiplication of the integration time (DIT), the number of integrations (NDIT), the number of polarimetric cycles (NPOL), and the number of half-wave plate orientations (NHWP = 4).

<sup>b</sup> Seeing measured by the differential image motion monitor (DIMM) at 0.5 µm.

<sup>c</sup> Wind speed at ground level.

<sup>d</sup> Coherence time measured by the multi-aperture scintillation sensor (MASS), except the first epoch which is estimated from the DIMM.

<sup>e</sup> H-band Strehl ratio estimated by SAXO.

<sup>f</sup> PSF FWHM fitted to the non-coronagraphic flux frames.

<sup>g</sup> Archival data from Stolker et al. (2016).

included a dark-frame subtraction, flat-field correction, and bad-pixel interpolation. To remove any residual background and bias from the images, we calculated the mean pixel value (with the inner 1′′5 masked) for each detector column and subtracted that value from all pixels in that column. The PSF of HD 135344B was fitted with a 2D Gaussian profile, which yielded a typical FWHM of 41 mas. More details on the observations and conditions are provided in Table 1.

## 2.2. Rapid Eye Mount (REM) Visible and Near-infrared Photometry

HD 135344B was observed with the REM at La Silla, Chile, in June 2016. The visible camera, ROSS2, is a simultaneous multi-channel imaging camera that delivers the $g'r'i'z'$ bands onto four quadrants of the same CCD detector. A dichroic enables simultaneous observations with the infrared camera, REMIR, with a similar field of view of $10' \times 10'$. Observations were executed with 3 s exposures in the $g'r'i'z'$ bands and 1 s exposures in the JHK-bands. For the JHK photometry, a standard five-position dither pattern was used and additional sky frames were obtained. Data was acquired during a total of twelve nights, but the data of two nights were rejected due to thick clouds.

The photometry of HD 135344B was measured differentially with respect to HD 135344A (SAO 206463), an A0V star with a separation of $21''$ from HD 135344B, which appears to be photometrically stable (Sitko et al. 2012). Differential photometry allowed us to measure with high precision the absolute fluxes of HD 135344B, also in variable conditions or with the presence of thin clouds. The JHK magnitudes of HD 135344A were retrieved from the 2MASS catalog (Skrutskie et al. 2006) while the $g'r'i'z'$ magnitudes were calculated through a transformation of the BVRI magnitudes with the relations from Jordi et al. (2006). The BVRI magnitudes of HD 135344A were obtained from Sitko et al. (2012): $B = 7.879 \pm 0.003$ mag, $V = 7.756 \pm 0.003$ mag, $R = 7.708 \pm 0.006$ mag, and $I = 7.662 \pm 0.004$ mag.

## 2.3. PIONIER Interferometry

We retrieved all of the available archival near-infrared interferometric data of HD 135344B from the Precision Integrated-Optics Near-infrared Imaging ExperIMENT (PIONIER)

instrument (Le Bouquin et al. 2011) at the Very Large Telescope Interferometer (VLTI). The data were taken during multiple epochs from 2011 to 2013.<sup>13</sup> The instrument recombines the four auxiliary telescopes that were positioned in the short (A1-B2-C1-D0), intermediate (D0-G1-H0-I1), and long (A1-G1-I1-K0) baseline configurations. The projected baseline, *B*, ranged from 7 to 135 m, enabling a maximum angular resolution of $\lambda/2B = 1.2$ mas across the seven spectrally dispersed channels in the H-band ($R \simeq 40$). Each observation of HD 134453B was preceded and followed by an observation of a calibration star to characterize the instrumental and atmospheric contribution to the visibilities and closure phases (i.e., the transfer function). Calibration stars were identified with the *SearchCal* tool (Bonneau et al. 2006, 2011). The data were reduced with the *pndrs* package, described in detail by Le Bouquin et al. (2011). The 27 calibrated OIFITS files (Pauls et al. 2005) are available in the *Optical interferometry DataBase* (OiDB) (Haubois et al. 2016).

The observing conditions were best during the two epochs in 2011 April, with a coherence time of $\sim 5$ ms and a seeing of $0''5-0''7$, while the conditions were significantly poorer when the other data sets were taken. We rejected low S/N measurements by selecting only the data points where the error estimates were in the first quartile of the total distribution of errors. The dominating factor in the error estimate is the stability of the transfer function during the night, which is determined by the temporal scatter of the calibrator visibilities. For each data set, the quality assessment was done independently for the six visibilities, four closure phases, and seven spectral channels. The selected measurements were retrieved from 11 out of 13 epochs, with 2011 May 26 and 2012 April 27 excluded.

# 3. Results

## 3.1. Multi-epoch Polarized Light Imagery

Scattered light images are displayed in chronological order (top to bottom) in Figure 1. Besides the newly obtained data sets from 2016, we also show the 2015 J-band imagery from Stolker et al. (2016). The first column shows the unscaled $Q_\phi$ images, defined such that positive values correspond to

<footnote n="13"><sup>13</sup> UT dates: 2011 April 27 and 29, 2011 May 26 and 27, 2011 June 3, 2011 August 7 and 8, 2012 March 6, 28, 29 and 30, 2012 April 27 and 2013 May 15.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Multi-epoch polarized scattered light images of HD 135344B in the J-band. The columns show from left to right the unscaled $Q_\phi$ images, unscaled $U_\phi$ images, $r^2$-scaled $Q_\phi$ images, and an unsharp-masked version of the $r^2$-scaled $Q_\phi$ images. The field of view of each image is $1''4 \times 1''4$ with north and east in the upward and leftward direction, respectively. The surface brightness of the images has been normalized to the integrated $Q_\phi$ flux (see the main text for details). The dynamical range of the color stretch is fixed in each column, except for the unsharp-masked images. The dynamical range of the $U_\phi$ images is a factor 10 smaller than the partner $Q_\phi$ images. Orange corresponds to positive values, blue to negative values, and black is the zero point. The extent of the coronagraph has been masked out. The major axis position angle of the outer disk, PA = 63° (purple line; van der Marel et al. 2015), and the inner disk, PA = $57.3 \pm 5.7$° (yellow line; see Section 3.4), are shown in the top row.</figcaption></figure>
azimuthally polarized flux. The images are normalized to the disk-integrated $Q_{\phi}$ flux, measured with an annulus aperture between 0.′′1 and 1.′′0, and shown with identical dynamical range. The second column shows the corresponding $U_{\phi}$ images, containing flux with a ±45° rotational offset of the direction of polarization with respect to the $Q_{\phi}$ flux. The third column contains the $Q_{\phi}$ images with a stellar irradiation correction applied (i.e., $r^2$-scaling). The fourth column shows unsharp-masked images that were obtained by smoothing the $r^2$-scaled $Q_{\phi}$ images with a Gaussian kernel ($\sigma = 200$ mas) and subtracting the smoothed images from the original $r^2$-scaled images. This procedure enhances the contrast of small-scale features by removing low spatial frequencies. The dynamical range of the unsharp-masked images is limited to positive values, and for each image separately normalized to the peak intensity.

The SW direction (PA = 180°–270°) of the disk, located around the major axis (see Figure 1), appears relatively bright in all $r^2$-scaled images. In contrast to the shadowing variations, the origin of that brightness enhancement is presumably intrinsic as the bright wedge in scattered light coincides with the (sub)millimeter emission peak of the crescent-shaped dust continuum (Pérez et al. 2014; Stolker et al. 2016). An enhancement of the surface density and/or midplane temperature will elevate the height of the scattering surface; therefore, a larger geometrical cross section of the disk surface is irradiated, which increases the scattered light flux. In this work, we focus on local brightness variations between the five epochs. We refer the reader to Muto et al. (2012), Garufi et al. (2013), and Stolker et al. (2016) for a detailed analysis and discussion of the scattered light detection of the spiral arms and cavity edge.

The $U_{\phi}$ images in Figure 1 contain contributions from noise, residuals of instrumental polarization, and possibly multiple scattered light from the disk. The residual signal could not be removed with the minimization steps explained in Section 2.1, but a more detailed analysis is required to determine if the remaining signal is an instrumental artifact. We distinguish between two different type of signals in the $U_{\phi}$ images that appear to be related to the total integration time and therefore the S/N. First, the relative contribution of noise is particularly well visible at separations ≳300 mas from the star in the last three epochs for which the total integration time was shortest, while the relative contribution is lower in the first two epochs. Second, the $U_{\phi}$ images show an enhanced signal within 300 mas of which the relative strength is larger in the first two epochs. The inner signal reaches only mildly above the background noise level in the remaining epochs. The inner $U_{\phi}$ signal appears to be variable between epochs, revealing a variety of brightness patterns. For example, the first epoch shows a complex pattern of multiple positive and negative lobes, whereas the second epoch displays an anti-symmetric signal which is bisected in a positive and negative side. The relative strength of the $U_{\phi}$ flux with respect to the $Q_{\phi}$ flux is quantified in Section 3.3.

Although we deem it likely that the remaining $U_{\phi}$ signal is an uncorrected instrumental artifact, we speculate that the increasing strength of the innermost $U_{\phi}$ signal with increasing S/N could be a result of multiple scattered light from the cavity edge. A minor fraction of the scattered light will be non-azimuthally polarized because of the small inclination of the outer disk (16°–20°; van der Marel et al. 2015, 2016).

Furthermore, a small fraction of the stellar light will scatter in the extended inner disk atmosphere (see Section 4.3) before it scatters from the outer disk. The $U_{\phi}$ signal from those photons will get modulated by the temporal variations in the inner disk, such that each epoch may show a different $U_{\phi}$ image. However, more detailed analysis and modeling is required in order to determine if the $U_{\phi}$ signal is real and if the temporal variations could be caused by the inner disk.

## 3.2. Asymmetric $Q_{\phi}$ Brightness Variations

A comparison of the $r^2$-scaled $Q_{\phi}$ images in Figure 1 shows epoch-to-epoch brightness variations. Azimuthal brightness minima are visible in all images with variations in their location, shape, and strength. The locality of the brightness minima (e.g., the shadow lane at PA = 169° on 2015 May 03) points toward a shadowing effect, likely caused by dust in the (unresolved) inner disk (Stolker et al. 2016). Furthermore, the brightness variations occur on a timescale similar to the dynamical timescale of the inner disk (the finest temporal resolution is eight days) and the variability timescale of the near-infrared photometry. Minor brightness variations in the unscaled $Q_{\phi}$ images are also visible in the region between the cavity edge and the coronagraph, possibly related to the flow of gas and small dust grains from the outer disk.

Azimuthal brightness variations, and their changes between epochs, are more evidently revealed with polar projections of the scattered light images, which are displayed in Figure 2. For clarity, we choose the unsharp-masked images for the identification of the shadow features in the polar projections. We caution that applying an unsharp mask may introduce a bias in the identification of brightness variations. However, the shadow features, which we will discuss below and are marked with arrows in Figure 2, are also visible in the regular $r^2$-scaled $Q_{\phi}$ images, but the contrast between shadowed and non-shadowed regions is smaller.

Asymmetric illumination/shadowing variations are visible in the scattered light imagery of all five epochs. Here we list the main characteristics of the shadow features, in consonance with the locations that are pointed out in Figure 2.

1. Epoch 1, 2015 May 03—Three narrow shadow lanes are present at position angles of 34°, 169°, and 304°, and an azimuthally broader dimming is visible in a north-northwest direction that is bound by two of the shadow lanes (Stolker et al. 2016).

2. Epoch 2, 2016 May 04—The eastern half of the disk is mildly shadowed, approximately in the position angle range of 10°–170°. Deeper shadows are superimposed near the edges of the global shadow. The deepening is particularly well visible in the north, extending at the cavity edge from PA = −50° to PA = 50°. Radially outward, the location of the shadow shows an azimuthal gradient. There is a hint of a localized shadow lane at PA = 170°, approximately colocated with the southern shadow lane in the first epoch.

3. Epoch 3, 2016 May 12—The bisection of the brightness distribution from the second epoch seems to have disappeared, although the poor observing conditions and the short total integration time (see Table 1) make the identification of the shadows challenging. The broad, northern shadow from the previous epoch is still present.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Polar projections of the r²-scaled, unsharp-masked $Q_{\phi}$ images shown in chronological order (top to bottom). North corresponds to PA = 0° and positive position angles are measured east from north. Localized and broad shadow features are indicated with solid and dashed arrows, respectively.</figcaption></figure>

There is also a hint of the broad, southern shadow, whereas the narrow southern shadow has disappeared.

4. Epoch 4, 2016 June 22—A broad shadow is present between PA = −90° and PA = 30° upon which finer shadow variations are superimposed, including narrow shadow lanes at the boundary with the non-shadowed region, similar to the northern shadow features in the first epoch. The shadow lane at PA = 30° possibly coincides with the location of the northeast shadow lane detected in first epoch.

5. Epoch 5, 2016 June 30—The cavity edge is shadowed between PA = −45° and PA = 40°, while shadowing of the exterior spiral arm only occurs from PA = 0° onwards. The radial extent of the broad shadow increases with increasing position angle, similar to the broad shadows in the first and second epoch. The shadow covers the full radial extent of the disk between PA = 0° and PA = 40°, similar to the shadow feature at the same location in the second and third epoch. The narrow shadow lanes from the fourth epoch seem to have disappeared.

The brightness depth of the shadow lanes appears typically larger than the broader shadows (e.g., epoch 1 in Figure 1), indicating larger density enhancements in the inner disk. Indeed, there is no evident radial dependence in the depth of the shadow lanes, even though the height of the outer disk increases with radius. The brightness depth of the broad shadows is typically deepest at the cavity edge but weakens toward larger radii (e.g., epoch 4 in Figure 2), presumably an effect of the increasing outer disk height. The broad shadow in north-northwest direction (PA ≃ −80°–30°) of all epochs shows an azimuthal gradient which is possibly a result of the asymmetric spiral arm perturbation of the disk surface.

The absence of strong azimuthal gradients, caused by a light travel time effect, in the location of the narrow shadows provides a lower limit on the radius from where the shadows are cast. For example, the position angle of a shadow will change by 10° between the inner disk and 80 au (500 mas) if the responsible dust clump is located at 0.15 au (see also Kama et al. 2016). Arguably, some of the shadow lanes show very minor tilts in Figure 2, but the precision and angular resolution of the observations challenge the identification of light travel time effects. The best candidate is the shadow lane at PA = 169° in the 2015 epoch. Stolker et al. (2016) speculated that its azimuthal tilt could be caused by orbital motion in the inner disk from which an orbital radius of 0.06 au (at 140 pc) was estimated. For the other shadow lanes, we may conclude that the dust clumps responsible are presumably located at distances ≳0.15 au.

A quantification of the azimuthal brightness variations is shown in Figure 3. The mean $Q_{\phi}$ flux is measured in position angle bins of 10° wide across a radial separation of 0′′1–0′′7 and divided by the angular area of a pixel. The polarized surface brightness (in counts s⁻¹ arcsec⁻²) is normalized to the total Stokes I flux (in counts s⁻¹), which is measured with a circular aperture on the unsaturated, non-coronagraphic flux frames after a correction for the integration time and response of the neutral density filter. The optimal aperture size (1′′5) was determined by measuring the photometric flux with a large range of aperture sizes (up to 3′′0), from which it was established that the total encompassed flux flattened for apertures larger than ~1′′5. The mean error bar is calculated from the standard error on the individual contrast points.

The polarized surface brightness contrast in Figure 3 shows typical values in the range of (2–6) × 10⁻³, except in the southwest direction where the contrast goes up to 8 × 10⁻³ at PA ≃ 240°. The integrated disk brightness consists mainly of signal from the
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Polarized brightness of the disk normalized to the total Stokes *I* flux. The plot shows the contrast from four of the epochs (colored dashed lines), the mean contrast (black solid line), and the relative optical depth variation between the minimum and maximum contrast (black dashed line, right y-axis). The gray shaded area covers the total variation of the contrast between the epochs. The mean error on the contrast, across all epochs and position angles, is shown on the bottom of the figure (see the main text for details). The image from 2016 May 12 has been excluded as it was affected by the poor observing conditions.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Visible and near-infrared photometry obtained with the REM during 10 nights in 2016 June. The mean and standard deviation of the fluxes are provided by the horizontally dashed lines and shaded regions, respectively. Photometric monitoring overlapped with the two most recent SPHERE epochs of which the UT dates are indicated with vertically dotted lines. The arrows point in the direction of an increasing and decreasing near-infrared variation.</figcaption></figure>

cavity edge and the spiral arms (see unscaled $Q_φ$ images in Figure 1), which are both asymmetrically shadowed resulting in dimming variations of 20%–45%. The contrast variation is minimal (10%) in the south-southwest direction, where the large-scale density and/or scale height enhancement affects the disk surface. The opening angle of the scattering surface is larger in that direction, such that shadowing by the inner disk requires dust to be located at higher altitude above the midplane. The contrast variations provide an upper

limit on the local changes in optical depth through the inner disk atmosphere. The stellar radiation that is transmitted through the atmosphere will be attenuated by a factor $e^{-\tau}$; therefore, the relative change in optical depth, $\Delta\tau$, can be calculated from the minimum and maximum contrast (see Figure 3). However, this only provides an upper limit on the optical depth variations because (i) the total flux does also affect the contrast and is correlated with the shadowing of the outer disk (see Section 4.3), (ii) part of the thermal emission from
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Top: integrated polarized scattered light contrast (black crosses) of the five J-band SPHERE epochs. Bottom: the integrated Stokes *I* flux (purple squares) and the disk-integrated *Q*<sub>φ</sub> flux (red circles), shown in arbitrary normalized units. The disk-integrated *U*<sub>φ</sub> flux (green circles) is computed from the absolute pixel values and shown relative to the *Q*<sub>φ</sub> flux with a factor 10 enhancement. The uncertainties are given at a 5σ level.</figcaption></figure>

the inner disk may illuminate the outer disk without being affected by
any local obscurations, and (iii) PSF smearing will lower the
brightness contrast between shadowed and non-shadowed regions.

3.3. Photometry and Scattered Light Contrast

Multi-epoch photometry in the g'r'i'z' and JHK-bands are displayed in Figure 4, covering 10 nights between 2016 June 07 and 2016 July 04 with nearly a daily sampling from June 16 till June 22. The error bars reflect the uncertainty on both the science and calibration star. The early June photometry in the g'r'i'z'-bands starts with a minor increase of 3%, remains approximately constant in mid-June, and increases with 2% by the end of June.

The temporal course of the JHK photometry appears more irregular with variations up to 10% with respect to the mean (horizontally dashed lines in Figure 4) in the JHK-bands. During the first half of 2016 June, the fluxes increased with approximately 6%, 9%, and 10% in the J-, H-, and K-bands, respectively, following the trend of the g'r'i'z' fluxes but with a larger fractional increase. In mid-June, the JHK fluxes also remained approximately constant but they increased further from June 22 onwards, in contrast to the g'r'i'z' fluxes. The final epoch shows a decline in the J- and H-bands, while the K-band photometry remained constant. Although the JHK variability in the first part of June seems correlated with the g'r'i'z' photometry, in the second half no correlation is apparent; that is, the JHK fluxes increased up to 10% with respect to the mean while the g'r'i'z' photometry remained constant. The total temporal coverage of the photometry is too short to reveal any trends and the sampling is too sparse to resolve possible variations on timescales less than one day.

In addition to the absolute REM photometry, we measured
the relative disk photometry of the SPHERE data, which is
presented in Figure 5. The disk-integrated polarized flux
was determined from the Q<sub>φ</sub> and U<sub>φ</sub> images with an annulus
aperture (0"1–1"0) centered on the star. The Stokes I flux was
measured with a circular aperture (1"5 radius) from the
dedicated flux frames. Absolute pixel values were used for the

$U_{\phi}$ photometry. The photometric contrast (top panel in Figure 5) is calculated as the ratio of the $Q_{\phi}$ and Stokes *I* flux after correcting for the difference in integration time and the response of the neutral density filter. Relative photometry allows for an epoch-to-epoch analysis without requiring an absolute flux calibration, assuming that both the coronagraphic sequence and the total flux data were obtained during similar observing conditions. The $Q_{\phi}$ and $U_{\phi}$ surface brightness errors are computed as the standard deviation within an aperture centered on each pixel with a radius of 62 mas (i.e., 1.5 resolution elements) and propagated accordingly to an integrated error (3%–5%). The uncertainty on the total flux (1%–2%) is computed in a background-limited region as the standard error of the sum, $\sigma_{\sqrt{N_{\text{pix}}}}$, with an annulus aperture equal in size to the aperture used for the Stokes *I* photometry.

The integrated contrast in Figure 5 varies between
$5.4 \times 10^{-3}$ and $7.2 \times 10^{-3}$. The photometry is not calibrated,
so only relative variations of $Q_{\phi}/\text{Stokes I}$ and $U_{\phi}/Q_{\phi}$ are
meaningful. The second epoch shows a consistent decrease of
both the $Q_{\phi}$ and Stokes *I* flux with respect to the first epoch. In
the third epoch, the contrast decreased by 20% possibly due to
the poor observing conditions (see Section 2.1), particularly
affecting the $Q_{\phi}$ photometry. During the last two epochs, the
relative increase of the total flux is large compared to the $Q_{\phi}$
flux, resulting in relatively low contrast. The increase of the
total flux during the last epochs seems consistent with the REM
photometry (see June 22 and 30 in Figure 4). The relative $U_{\phi}$
photometry shows an increase in the second epoch due to the
residual within 200 mas (see Section 3.1), possibly caused by
multiple scattered light from the inner and/or outer disk, while
the relative $U_{\phi}$ photometry in the third epoch is larger due to
the enhanced noise residual.

3.4. Parametric Model Fitting of the Visibilities

The normalized, squared visibilities, $V^2$, across the spectrally
dispersed H-band channels of the multi-epoch PIONIER
observations are displayed in the left panel of Figure 6. The
visibilities decrease continuously with increasing spatial
frequency (i.e., $B/\lambda$), indicating that the region from which
the H-band flux originates is resolved. At the longest baselines,
there appears no turnover point to an asymptotic value of the
stellar flux, so the circumstellar emission is not over-resolved.
The closure phases, $\Delta\phi$, are consistent with zero within the
error bars ($|\Delta\phi| \lesssim 3^\circ$), therefore, the brightness distribution of
the inner disk is point symmetric within the uncertainties and at
the spatial resolution probed by the observations. There is no
significant dispersion visible for visibility points obtained with
baselines of similar lengths but different position angles, which
implies that the inclination of the inner disk is small. The
diagonal scatter of the visibilities is an effect of chromatic
dispersion (Lazareff et al. 2017). Coverage of the (u,v)-plane is
shown in the right panel of Figure 6.

The orientation and characteristic radius of the inner disk H-band emission is determined by fitting, in Fourier space, a parametric model to the visibilities, following the procedure described in detail by Lazareff et al. (2017). This allowed us to apply a $\chi^2$ minimization and to assign formal error bars to the inferred parameter values. We parameterized the H-band emission with an elliptical brightness distribution that is radially parameterized by a weighted combination of a Gaussian and pseudo-Lorentzian profile. The inner rim is not
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Left: squared visibilities ($V^2$) of in the VLTI/PIONIER H-band channels with $R \simeq 40$ (top) and the fitting residuals of the ellipsoidal brightness model (bottom). Right: coverage of the (u,v)-plane, shown with the same color coding as the visibilities.</figcaption></figure>

fully resolved by the longest baselines, which justifies an ellipsoidal distribution instead of a broadened ring.

The best-fit model ($\chi^2 = 1.10$) corresponds to an inclination and major axis position angle of $18^\circ 2_{-4.1}^{+3.4}$ and $57^\circ 3 \pm 6^\circ 3$, respectively. The values are, within the uncertainties, very similar to those of the outer disk (see the discussion in Section 4.1.1). The best-fit position angle is shown in Figure 1 in comparison with the outer disk value from van der Marel et al. (2015). The half-flux semimajor axis is $0.71 \pm 0.03$ mas (=0.11 au), which implies that a significant fraction of the H-band emission originates from within the silicate sublimation radius ($R_{\text{sub}} = 0.2$ au; Carmona et al. 2014). A detailed overview of the fitting results is provided in the Appendix, where the best-fit values of all parameters are listed with their dependence on the cutoff level of the selection criterion for the (u,v) points.

We caution that the visibilities were combined from multiple epochs, while the inner disk is variable on a timescale of days or less (see Section 3.3). The visibilities in Figure 6 depend on the absolute H-band flux and the relative contributions of the star and disk. Therefore, an additional uncertainty has been introduced by combining the visibilities from multiple epochs. While the absolute H-band flux is variable and not measured at the nights of the observations. Also, we made the assumption that the orientation of the inner disk did not change due to precession between those epochs.

# 4. Discussion

In Section 4.1 we will discuss the available constraints on the orientation of the inner and outer disk, as well as the presence of the broad, quasi-stationary shadow. In Section 4.2, we will provide an observational perspective on the variability and discuss several processes that may affect the inner disk dynamics. In Section 4.3, we will estimate the extent of the inner disk variations by quantifying the correlation between the thermal emission and the scattered light flux with a grid of radiative transfer models.

## 4.1. Constraints on the Inner Disk (mis)Alignment

### 4.1.1. Near-infrared and Submillimeter Interferometry

The relative orientation of the inner and outer disk is determined by their inclination with respect to the plane of the sky and their position angle, with the misalignment defined by the angle between the normal vectors of the two midplane orientations. The outer disk’s orientation of HD 135344B has been determined in several studies (see Carmona et al. 2014, for an overview). Here we list those values that were derived from spatially resolved (sub)millimeter CO observations: $i = 11^\circ 5 \pm 0^\circ 5$ and PA = $64^\circ \pm 2^\circ$ (Lyo et al. 2011), $i = 20^\circ$ and PA = $63^\circ$ (van der Marel et al. 2015), $i = 16^\circ$ and PA = $63^\circ$ (van der Marel et al. 2016). In Section 3.4, we found that the inner disk inclination and position angle are $i = 18^\circ 2_{-4.1}^{+3.4}$ and PA = $57^\circ 3 \pm 5^\circ 7$, respectively. The values are consistent with the orientation of the outer disk within the uncertainties of the model fitting and the available values for the outer disk. This result contrasts several earlier studies that suggested a significant misalignment between the inner and outer disk (Fedele et al. 2008; Grady et al. 2009; Stolker et al. 2016). However, we can not exclude a minor disk warp given the uncertainties on both the inner and outer disk orientation.

A possible misalignment of the two disk components relies also on the identification of the near and far side of both the inner and outer disk. Spatially resolved observations of CO gas show a redshifted and blueshifted velocity in the southwest and northeast direction, respectively (e.g., Pérez et al. 2014). This implies that the near side of the outer disk is along the southeast direction of the minor axis if we assume that the spiral arms in scattered light follow a trailing motion. For the inner disk, there is no direct constraint on the near and far side at the angular resolution of the PIONIER observations. However, the misalignment will be ~38° if the near side of the inner disk is in a northwest direction, a scenario that can be excluded from the absence of two stationary shadow lanes similar to HD 142527 ($\Delta\theta = 70^\circ$; Marino et al. 2015) and HD 100453 ($\Delta\theta = 72^\circ$; Benisty et al. 2017). This means that the near side
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Raytraced $Q_{\phi}$ images of the DIANA radiative transfer model of HD 135344B (Woitke et al. 2016). Left: the inner disk is aligned with the outer disk. Right: the misalignment of the inner disk is $2^{\circ}6$, by adopting the best-fit inclination and position angle from Section 3.4. Images have been convolved with a Gaussian kernel (FWHM = 41 mas) to match the angular resolution of the SPHERE imagery. Surface brightness values are provided along the major and minor axis direction of the outer disk.</figcaption></figure>

of the inner disk is located in a sourtheast direction and a possible misalignment will be minor.

### 4.1.2. Quasi-stationary Shadowing by a Minor Disk Warp?

The location and width of the shadows provide further constraints on the orientation of the inner disk with respect to the outer disk. Three scenarios remain possible: (i) the inner disk and outer disk are aligned, such that shadowing occurs only through uplifting of dust in the inner disk atmosphere, (ii) a minor misalignment ($\Delta\theta \sim 1^{\circ}-2^{\circ}$) is present that might not cast a stationary shadow, but additional variations in the inner disk will cause preferential shadowing in the direction where elevation of the inner disk above the outer disk midplane is largest, (iii) an intermediate misalignment ($\Delta\theta \sim 2^{\circ}-10^{\circ}$) is present which casts a broad, stationary shadow when the misalignment of the inner disk is similar to the opening angle of the scattering surface of the outer disk. An example of the third scenario is the broad shadow detected with HST on the TW Hya disk ($\Delta\theta = 8^{\circ}$; Rosenfeld et al. 2012).

The inclination of the inner disk could be either smaller or larger than the outer disk given the uncertainties on its orientation (see Section 2.3) that could result in a broad shadow in an approximately northwest or southeast direction, respectively. Interestingly, a broad shadow seems present in all of the images in a north-northwest direction (see Figure 2), although with small variations in its precise location, shape, and depth. This may imply that the inclination of the inner disk is slightly smaller (i.e., more face-on) than the outer disk. Variable fine structure is present in the north-northwest shadow, which requires additional optical depth variations through the atmosphere of the inner disk as will be discussed in more detail in Section 4.2.

A broad dimming was also present in the northwest direction of the Subaru/HiCIAO H-band imagery by Muto et al. (2012) and the VLT/NACO H and $K_s$-bands imagery by Garufi et al. (2013). This was interpreted in both studies as a depolarization effect that can occur when the disk is inclined because the

polarization efficiency peaks around the major axis where scattering angles are close to $90^{\circ}$ (e.g., Murakawa 2010; Min et al. 2012). The inclination of the outer disk is relatively small, so a strong depolarization effect might not be expected. We speculate that the broad dimming in the HiCIAO and NACO imagery might be the same quasi-stationary shadowing effect that is seen in the SPHERE imagery, possibly related to a minor disk warp.

To illustrate this scenario, we have adopted the DIsc ANAlysis (DIANA; Woitke et al. 2016) radiative transfer model of HD 135344B. The left image in Figure 7 displays the raytraced $Q_{\phi}$ image for a setup in which the inner disk is aligned with the outer disk ($i = 20^{\circ}$, $PA = 63^{\circ}$; van der Marel et al. 2015). A mild depolarization effect is seen in both directions of the minor axis, but the effect is slightly stronger on the near side (southeast) due to the flaring geometry of the disk surface (Min et al. 2012). In the right image of Figure 7, we changed the orientation of the inner disk to the best-fit values from Section 3.4. The misalignment here is $2^{\circ}6$; that is, comparable to the disk warp seen in the debris disk around $\beta$ Pic ($\Delta\theta = 4^{\circ}6$; Heap et al. 2000). The atmosphere of the inner disk casts a mild shadow on the outer disk in a northwest direction, reminiscent of the broad shadow on the cavity edge of the SPHERE images in Figures 1 and 2. In the opposite direction, the outer disk becomes more strongly irradiated, which introduces an azimuthal brightness modulation along the cavity edge (see also Rosenfeld et al. 2012), although intertwined with the modulation by the polarization efficiency of the outer disk.

## 4.2. Origin of the Shadows and their Variability

### 4.2.1. Observational Perspective and Variability Timescale

The most prominent azimuthal brightness variations were identified in Section 3.2 and we noticed that the shadows can be classified into two categories, (i) localized shadow lanes and (ii) broader shadows that are tens of degrees wide. The shadow variations are presumably caused by small density
enhancements in the atmosphere of the inner disk with the optical depth variations being approximately smaller than unity (see Figure 3). Although a broad shadow is present in the north-northwest direction of all epochs, none of the shadows are fully stationary. This indicates that the atmosphere of the inner disk is a dynamical environment in which the gas and the micron-sized dust grains (which are dynamically coupled) are subjected to processes that change their vertical distribution on a fast timescale.

Variability of the shadows might be caused by orbital and/or vertical motion of dust enhancements in the inner disk. The vertical response of the disk occurs approximately on the Keplerian timescale, meaning that a local perturbation of the disk will settle to an equilibrium state within approximately one orbit. However, the orbital timescale provides only a lower limit on the response timescale as it depends on the heating and cooling timescales of the gas. The H-band flux is emitted from a characteristic radius of 0.11 au (see Section 3.4), which corresponds to a Keplerian timescale of 10 days, similar to the finest temporal sampling of the SPHERE imagery (eight days). However, the shadows are presumably cast further outward as inferred from the absence of significant light travel time effects (see Section 3.2). The temporal sampling of the SPHERE observations is too sparse to determine the timescale by which the shadow features appear and disappear, therefore, disentangling variability due to orbital motion and the (dis)appearance of shadows is not possible.

The JHK fluxes in Figure 4 show variations up to 10% on a timescale of days to weeks (see also Grady et al. 2009; Sitko et al. 2012), while the $g'r'i'z'$ fluxes varied only by 1%–2%, indicating that mainly the thermal emission from the inner disk is affecting the near-infrared variability. For reference, the J-band flux of HD 135344B consists of 69% stellar radiation and 31% inner disk emission (see Section 4.3). The fast variability timescale of both the shadows and the near-infrared photometry may point toward a common origin in the inner disk.

Small variations of the visible photometry might be related to episodic accretion events, photospheric/chromospheric activity, stellar pulsations, minor attenuation variations by an (optically thin) dust envelope, or increased scattering from the inner disk. The scattered light flux from the inner disk is in the optical ~1% of the total flux, so structural changes in the inner disk may cause a change both in thermal emission and scattered light. In that case, the photometric variations in the visible should be correlated with variations in the near-infrared, which seems the case in the first half of the REM photometry but not the second half.

### 4.2.2. Inner Disk Processes Affecting the Dust Dynamics

Several processes may have an effect on the dynamics and distribution of the gas and dust in the inner disk, possibly causing shadow variations on the outer disk. For example, hydrodynamical fluctuations, such as turbulent eddies and filaments, may produce short-lived obscuration events (Dullemond et al. 2003; Flock et al. 2017). Or, catastrophic collisions between planetesimals, possibly stirred-up by a planet (Kenyon & Bromley 2004), will locally enhance the dust density, although a gaseous environment will affect the dust dynamics differently than in a debris disk. Dedicated simulations are required to determine if such processes could produce disk perturbations that are localized and strong enough to

explain the narrow shadows. Precession of a disk warp, for example driven by a companion (Lai 2014), will result in a variable location of the casted shadow (Debes et al. 2017). The broad shadow in the north-northwest direction of HD 135344B appears approximately stationary, so a fast precession can be excluded if the shadow is cast by a minor disk warp.

Variations in the inner disk might also be related to star-disk interactions. HD 135344B is an F4V-type star with a weak magnetic field ($\langle B_z \rangle = 32 \pm 15$ G; Hubrig et al. 2009), so a magnetic coupling to and warping of the inner disk seems unlikely. Indeed, with the absence of a breaking mechanism, the star has been able to spin-up to a near break-up rotational velocity (Müller et al. 2011). The fast rotation might drive a viscous decretion disk, similar to classical Be stars (Rivinius et al. 2013), by which gas from the stellar atmosphere spreads outward, possibly creating disturbances in the inner disk by an interaction with the inward accretion flow. Also, accretion may play a role in the distribution of material in the inner disk. Fairlamb et al. (2015) measured a rate of $10^{-7.4} M_{\odot} \text{yr}^{-1}$, and Sitko et al. (2012) determined a factor of two variation during the course of a few months. Therefore, accretion of gas and dust from the outer disk might be an irregular process, possibly mediated by one or multiple companions inside the large dust cavity, such that the inner disk gets asymmetrically replenished.

The near-infrared emission, hydrogen line fluxes, and the He I line profile are all variable on various timescales, related to the processes occurring in the inner disk and near the stellar surface (Grady et al. 2009; Sitko et al. 2012). The variable P Cygni profile of the He I line is likely related to a wind with an orientation that changes on timescales of a day or less (Sitko et al. 2012), which is launched in the star-disk interaction region (e.g., Edwards 2009). Micron-sized dust grains in the inner disk atmosphere could be entrained by a photoevaporative wind that is driven by the UV radiation of the star (Owen et al. 2011; Hutchison et al. 2016), or dust could be uplifted by a centrifugally driven disk wind (Bans & Königl 2012). An extended low-density atmosphere could be supported by the magnetic field of the inner disk, which may explain a large near-infrared excess and possible shadowing (Turner et al. 2014). Alternatively, the central star could drive a wind from the circumplanetary disk of a planetary companion (Tambovseva et al. 2006), or disk perturbations by a companion on an inclined orbit may also cause an asymmetric illumination of the disk (Demidova et al. 2013). Three-dimensional radiation nonideal magnetohydrodynamical simulations show turbulent velocities in the inner disk up to 10% of the sound speed and a nonaxisymmetric shadow on the outer disk cast by a dead zone-induced vortex (Flock et al. 2017). An inner disk vortex will orbit the star with a Keplerian velocity, so it is likely not responsible for the broad, quasi-stationary shadow.

### 4.2.3. Face-on Variant of the UX Orionis Phenomenon?

Although the origin of the shadows remains uncertain, we notice a resemblance between the photometric variations of UX Orionis stars (UXORs) and the spatially resolved shadow variations on the disk surface around HD 135344B. UXORs are a subclass of Herbig Ae/Be stars that are characterized by sudden declines in brightness up to several magnitudes in the optical, associated with increased extinction and polarization, which suggests changes of the column density in the line of sight toward the star (Waters & Waelkens 1998). It has been proposed that such photometric variations could be caused by
orbiting dust clouds when the disk is observed almost edge-on (Grinin et al. 1994), although alternative explanations involving disk winds (Grinin & Tambovtseva 2003) and turbulent filaments exist (Dullemond et al. 2003). Similar processes could be invoked to explain the variability of UXORs and HD 135344B, which may suggest that HD 135344B is a face-on variant of the UXOR phenomenon.

### 4.3. Radiative Transfer Models of a Shadowed Outer Disk

Shadow variations on the outer disk depend on changes in the vertical distribution of dust in the inner disk. More specifically, the strength of the shadowing is set by the inner disk radius, where the transition region from optically thin to thick occurs highest above the midplane. The height of this transition region is mainly determined by the pressure scale height and the surface density if the dust opacities are constant throughout the inner disk. For a flaring geometry, the thickest part will be close to the outer edge of the inner disk. Alternatively, a puffed-up inner rim may shadow the outer disk in case the exterior of the inner disk is fully shadowed by the rim (Dullemond & Dominik 2004; Dong 2015).

#### 4.3.1. Parametric Model Setup

To understand quantitatively the effect of the inner disk on the near-infrared excess and the scattered light contrast, without making an assumption about the origin of the shadows, we have constructed a grid of 20 × 20 radiative transfer models that are meant as a proof of concept rather than an accurate fit to the data. We adopted the DIANA model setup of HD 135344B, which provides a multiwavelength fit to the SED and several other gas and dust observables (Woitke et al. 2016). The inner disk is aligned with the outer disk and ranges from 0.16 au to 0.21 au (at 140 pc), slightly beyond the characteristic radius inferred from the H-band visibilities (see Section 3.4), with a surface density profile parameterized as $\Sigma \propto r^{-1.8}$ and a sharp inner rim. The pressure scale height profile is parameterized as $H \propto r^{-0.08}$ with the normalization provided by a reference aspect ratio, $H_0/r_0$, at $r_0 = 0.2$ au. The negative flaring index implies a decrease of the scale height with increasing radius, which is expected due to direct heating of the inner rim by the star. The radial extent of the inner disk is small, so the flaring index has only a minor impact on the disk geometry. The grid covers values of the inner disk dust mass in the range of $10^{-11}-10^{-7} M_\odot$, and the inner disk aspect ratio, $H_0/r_0$, in the range of 0.005–0.25. The flaring index of the outer disk is 1.14, with a reference aspect ratio of 0.16 at 50 au. The radiative transfer and image raytracing were done with `MCMax3D`, a Monte Carlo continuum radiative transfer code (Min et al. 2009).

#### 4.3.2. Polarized Scattered Light versus Thermal Emission

For each model, we computed the disk-integrated $Q_\phi$ intensity of the outer disk and the total J-band flux, which are displayed in the top panel of Figure 8. Several effects of the inner disk mass and the aspect ratio on the polarized intensity are evident in the colored map. Increasing the aspect ratio from very small values up to ~0.05 results in a larger fraction of the stellar light being reprocessed by the inner disk and reemitted in the J-band, thereby increasing the scattered light flux from the outer disk. Similarly, the increase of the polarized intensity with increasing dust mass is also the result of a larger fraction

<caption>**Figure 8.** Top: radiative transfer simulations of the $Q_\phi$ intensity of the outer disk (colored map) and the J-band photometry (white contours). Center: scattered light contrast (colored map) with the mean and $1\sigma$ of the REM J-band flux (black dashed lines) and the contrast of the SPHERE imagery (white dashed lines). The white cross denotes to the DIANA model of HD 135344B (Woitke et al. 2016). Bottom: spectral energy distributions (SEDs) of all models (black solid lines), superimposed by SED of the DIANA model (yellow dashed line) and the photometry (red points; Carmona et al. 2014). Error bars of the photometry have been excluded when they are smaller than the symbol.</caption>

of thermal radiation from the inner disk scattering from the outer disk. In this regime of the aspect ratio, the opening angle of the inner disk is too small to shadow the outer disk.

Shadowing of the outer disk by the inner disk starts to have an effect for $H_0/r_0 \gtrsim 0.05$, such that the opening angle of the inner disk atmosphere is comparable to or larger than the opening angle of the scattering surface of the outer disk. The polarized intensity decreases with increasing aspect ratio when the inner disk dust mass is smaller than $\sim 10^{-9} M_\odot$ because the optical depth through the inner disk atmosphere toward the outer disk increases. Consequently, a larger fraction of the stellar light is attenuated by extinction in the inner disk. For inner disk masses larger than $\sim 10^{-9} M_\odot$, the effect of

<figure><img src="image_9.png" /></figure>
enhanced reprocessing by the inner disk starts to dominate over the shadowing effect. More specifically, even though increas- ing the aspect ratio beyond 0.05 will increase the amount of shadowing of the outer disk, enhanced reprocessing by a larger amount of dust in the inner disk results in an increased irradiation of the outer disk in the J-band. Since the exponent is negative for both the surface density and the scale height profile, the inner disk is thickest at the location of the inner rim (i.e., the radius where the vertical τ = 1 surface is located highest above the midplane). As a result, reprocessed stellar light by the inner rim can directly irradiate the outer disk, in contrast to an inner rim which is partly shielded from the outer disk by a flaring scale height profile.

The J-band photometry of the radiative transfer models is superimposed on the polarized surface brightness in the top panel of Figure 1. The inner disk dust mass and aspect ratio have a similar effect on the J-band flux because both parameters affect the amount of reprocessing by the inner disk. While the stellar J-band flux remains constant, the thermal emission from the inner disk is positively correlated with both the inner disk dust mass and the reference aspect ratio. The SEDs of all models are shown in the bottom panel of Figure 1. As expected, the near- and mid-infrared excess varies greatly in the explored parameter space.

### 4.3.3. Observed and Simulated Scattered Light Contrast

An absolute flux calibration of the SPHERE data is not possible with high-precision, so we use a relative brightness measurement to compare the results from the grid of models with the observations. The central panel of Figure 8 displays the polarized surface brightness contrast that has been computed with the same procedure described in Section 3.2, here simply given by the ratio of the mean polarized intensity of the outer disk and the total J-band flux. The J-band flux and shadowing of the outer disk are interlinked in the lower part of the disk mass regime; therefore, the contrast distribution across the grid appears similar to the individual components in the top panel. In the regime of the large inner disk mass and aspect ratio, the contrast balances between the increasing scattered light flux and the increasing J-band photometry. The latter dominates over the former, such that the contrast decreases toward the top right of the distribution, yet the polarized intensity has a significant effect on the contrast. This means that using a relative flux measurement instead of an absolute one may lead to a degeneracy in the interpretation because a low contrast could correspond to either a high or low polarized intensity.

The azimuthally and temporally averaged scattered light contrast of the SPHERE observations, as well as the REM J-band photometry, are superimposed on the contrast from the radiative transfer models in the central panel of Figure 8. We rejected the third quadrant because the brightness of the southwest part of the disk is mainly determined by the enhanced midplane density, which would otherwise bias the characteristic contrast variations by the shadowing. Significant changes in the total dust mass of the inner disk are not likely to occur, as the local viscous timescale will be multiple orders of magnitude longer than the timescale probed with the SPHERE observations. For a fixed dust mass of $9.86 \times 10^{-10} M_{\odot}$ from the DIANA model, the observed scattered light contrast corresponds to aspect ratio values in the range of 0.11–0.22, likely beyond what is expected from a hydrostatically

supported disk. On the other hand, the observed J-band fluxes cover a lower aspect ratio regime (0.07–0.1). A less extended atmosphere will be required when a small misalignment is present between the inner and outer disk (see Section 4.1), in which case the estimated value of $H_0/r_0$ will be smaller.

### 4.3.4. Evidence for Extended Variations in the Inner Disk

We infer from Figure 8 that extended variations in the inner disk atmosphere ($H/r \lesssim 0.2$) can explain the variations of the scattered light contrast and the related shadowing. However, the predicted J-band flux from an extended atmosphere is too large compared to the REM photometry. The discrepancy between the scattered light contrast and J-band flux possibly points to the uncertain origin of the near-infrared excess. We want to stress that the goal of the radiative transfer models is to provide a quantitative estimate of the required extent of the inner disk atmosphere without making an assumption about the origin of the near-infrared excess (e.g., super-refractory grains, disk wind, magnetically supported atmosphere). Nonetheless, alternative disk structures or dust properties could be considered to dissolve the discrepancy, for example, an inner disk consisting of an optically thick but geometrically thin component and an optically thin but geometrically extended component. The relative strength of the near-infrared excess and the outer disk shadowing is also sensitive to the abundances of amorphous silicates, amorphous carbon, and metallic iron. Amorphous laboratory silicates have a single scattering albedo close to unity (Dorschner et al. 1995), while the absorption cross-section of carbon and iron grains is larger (Zubko et al. 1996). Therefore, increasing the relative amount of silicate grains will lower the near-infrared excess, while the extinction through the atmosphere may remain large enough to cast a shadow on the outer disk.

We also tested the dependence of the turbulence mixing strength, set by the dimensionless viscosity parameter $\alpha$, on the scattered light contrast and near-infrared excess. The parameter controls the dust settling in the radiative transfer models by assuming an equilibrium between upward turbulent mixing and downward gravitational settling (Woitke et al. 2016). The turbulence/settling parameter mainly affects the brightness contrast when its value is $\lesssim 10^{-4}-10^{-5}$, such that even the micron-sized grains slightly settle. The maximum $\alpha$ value affecting the brightness contrast increases with increasing aspect ratio. The dependence on the flaring index of the inner disk was also tested. As expected, for a given reference aspect ratio there appeared no significant dependence on the flaring index because of the narrow radial extent of the inner disk.

## 5. Summary and Conclusions

We have presented a multi-epoch scattered light study of the protoplanetary disk surrounding HD 135344B, an isolated pre-main-sequence F4V-type star in the Scorpius OB2-3 association. Polarimetric differential imaging obser- vations with VLT/SPHERE in the J-band revealed, with a spatial resolution of ~6.4 au, azimuthal shadowing variations on the outer disk related to the vertical dust distribution of dust in the inner disk. The imagery shows irregularly changing shadow patterns between all epochs, although similarities have been identified. Shadows appear both as localized lanes and broader structures, typically colocated, which likely trace small-scale perturbations and large-scale
<table><caption>**Table 2**
Visibility Fitting Results</caption><thead><tr><th>Cutoff (%)</th><th>n<sub>V</sub></th><th>χ<sup>2</sup><sub>r</sub></th><th>i (degree)</th><th>PA (degree)</th><th>a (mas)</th><th>k<sub>c</sub></th><th>f<sub>c</sub></th><th>f<sub>L</sub></th></tr></thead><tbody><tr><td>25</td><td>199</td><td>1.10</td><td>18.2<sup>+3.4</sup><sub>-4.1</sub></td><td>57.3 ± 5.7</td><td>0.71 ± 0.03</td><td>-3.34 ± 0.57</td><td>0.61 ± 0.04</td><td>0.20 ± 0.02</td></tr><tr><td>35</td><td>278</td><td>1.44</td><td>19.9<sup>+3.1</sup><sub>-3.7</sub></td><td>57.3 ± 6.3</td><td>0.76 ± 0.03</td><td>-3.04 ± 0.57</td><td>0.56 ± 0.03</td><td>0.17 ± 0.03</td></tr><tr><td>45</td><td>358</td><td>1.67</td><td>19.9<sup>+3.1</sup><sub>-3.7</sub></td><td>56.1 ± 5.7</td><td>0.74 ± 0.03</td><td>-2.81 ± 0.60</td><td>0.57 ± 0.04</td><td>0.18 ± 0.03</td></tr><tr><td>55</td><td>437</td><td>1.85</td><td>21.6<sup>+2.9</sup><sub>-3.4</sub></td><td>56.7 ± 6.3</td><td>0.74 ± 0.05</td><td>-2.79 ± 0.68</td><td>0.58 ± 0.04</td><td>0.18 ± 0.03</td></tr><tr><td>65</td><td>517</td><td>2.07</td><td>18.2<sup>+3.4</sup><sub>-4.1</sub></td><td>56.7 ± 7.4</td><td>0.72 ± 0.05</td><td>-2.70 ± 0.59</td><td>0.58 ± 0.04</td><td>0.18 ± 0.03</td></tr><tr><td>75</td><td>597</td><td>2.28</td><td>19.9<sup>+3.1</sup><sub>-3.7</sub></td><td>51.6 ± 9.7</td><td>0.78 ± 0.05</td><td>-2.64 ± 0.62</td><td>0.55 ± 0.04</td><td>0.15 ± 0.03</td></tr><tr><td>85</td><td>676</td><td>2.76</td><td>18.2<sup>+3.4</sup><sub>-4.1</sub></td><td>48.1 ± 10.3</td><td>0.74 ± 0.05</td><td>-2.39 ± 0.70</td><td>0.57 ± 0.05</td><td>0.18 ± 0.03</td></tr><tr><td>95</td><td>756</td><td>2.99</td><td>18.2<sup>+3.4</sup><sub>-4.1</sub></td><td>45.3 ± 9.1</td><td>0.76 ± 0.05</td><td>-2.66 ± 0.77</td><td>0.57 ± 0.06</td><td>0.16 ± 0.04</td></tr><tr><td>100</td><td>796</td><td>3.07</td><td>19.9<sup>+3.1</sup><sub>-3.7</sub></td><td>49.8 ± 9.7</td><td>0.74 ± 0.05</td><td>-1.95 ± 0.76</td><td>0.58 ± 0.05</td><td>0.20 ± 0.04</td></tr></tbody></table>

Note. Table columns (from left to right): cutoff level in the cumulative distribution of the error estimates, number of (u,v) points, reduced χ², inclination, position angle of the major axis, half-flux semimajor axis, spectral index of the circumstellar component, fractional flux of the circumstellar component, and weighting factor of the pseudo-Lorentzian profile. The uncertainties are provided at a 1σ level. Most parameter values appear stable over a broad range of quality cutoff levels.

dynamics, respectively. A broad, quasi-stationary shadow is present in the north-northwest direction of all scattered light images, in particular well visible at the cavity edge where the outer disk scattering surface is lowest. This might suggest that the inner disk is misaligned by several degrees, requiring additional optical depth enhancements through the inner disk atmosphere to shadow the disk further outward. However, fitting of a parametric brightness model to VLTI/PIONIER H-band visibilities provided a best-fit inclination and position angle which is, within the uncertainties, consistent with the outer disk. The photometry showed only minor variations (1%–2%) in the g'r'i'z'-bands, while the JHK fluxes varied up to 10%, indicating significant changes in the amount of reprocessing of stellar light by the inner disk. Variability of the photometry and shadowing appear to be correlated, and both are related to structural changes in the inner disk which we have quantified with a grid of radiative transfer models. The observed variations in scattered light contrast require extended variations in the inner disk atmosphere (H/r ≲ 0.2), beyond what is inferred from the near-infrared excess alone, highlighting the uncertainty about the origin of the near-infrared excess. The variability of the shadows is likely related to the structure and dynamics of the inner disk, the H-band emission of which originates from a characteristic radius of 0.11 au, that is, inside the silicate sublimation zone. Asymmetric shadowing variations by the inner disk might be caused by mechanisms such as turbulent fluctuations, planetesimal collisions, a dusty disk wind, or asymmetric and episodic accretion. Simultaneous observations of scattered light, photometry, spectroscopy, and/or near-infrared interferometry will provide more stringent constraints on the driving processes in the inner disk for which an approximate daily sampling is required to trace the fast disk dynamics.

We would like to thank Doug Lin, Mario Flock, and the other attendees of the Ringberg workshop on “The Atmospheres of Disks and Planets” for insightful discussions. We also thank the referee for providing constructive comments and suggestions. This research has made use of the OiDB<fnref n="14" /> and SearchCal<fnref n="15" /> services at the Jean-Marie Mariotti Center.

M.B. acknowledges funding from ANR of France under contract number ANR-16-CE31-0013 (Planet Forming Disks). S.P. acknowledges support from CONICYT-Gemini grant 32130007. S.K. acknowledges support from an STFC Rutherford fellowship (ST/J004030/1) and an ERC Starting Grant (Grant Agreement No. 639889).

*Facilities:* VLT:Melipal (SPHERE), VLTI (PIONIER).

*Software:* MCMax3D (Min et al. 2009), SPHERE DRH (Pavlov et al. 2008), SearchCal (Bonneau et al. 2006, 2011), Astropy (Astropy Collaboration et al. 2013).

## Appendix
### Details on the Visibility Fitting

The PIONIER H-band visibilities can be expressed as (Lazareff et al. 2017),

$$V(u, v, \lambda) = \frac{f_s(\lambda_0/\lambda)^{k_s} + V_c(u, v)f_c(\lambda_0/\lambda)^{k_c}}{f_s(\lambda_0/\lambda)^{k_s} + f_c(\lambda_0/\lambda)^{k_c}}, \quad (1)$$

where $V_c$ is the visibility of the circumstellar component, $k_s$ and $k_c$ are the spectral index of the stellar and circumstellar component, respectively, $f_s$ and $f_c$ are the fractional flux of the stellar and circumstellar component at a reference wavelength $\lambda_0$, and $\lambda$ is the wavelength of the spectral channel. The flux fractions of the star and circumstellar disk are normalized such that $f_s + f_c = 1$.

As described in Section 3.4, the visibilities are fitted with an ellipsoidal brightness distribution that is parameterized by a weighted combination of a Gaussian, $\mathcal{F}_G(r)$, and a pseudo-Lorentzian, $\mathcal{F}_L(r)$, radial distribution that allows for some freedom in the steepness of the asymptotic decay. The brightness distributions are derived from

$$\mathcal{F}_G(r) = \frac{\ln 2}{\pi a^2} \exp \left[ - \left( \frac{r}{a} \right)^2 \ln 2 \right], \quad (2)$$

$$\mathcal{F}_L(r) = \frac{a}{2\pi\sqrt{3}} \left( \frac{a^2}{3} + r^2 \right)^{-3/2}, \quad (3)$$

followed by an anamorphosis along the minor axis, and where r is the emission radius and a is the semimajor axis of the half-light isophote. The Hankel transforms of the brightness distributions take a simple analytical form. The visibility of

<footnote n="14">¹⁴ Available at http://oidb.jmmc.fr.</footnote>

<footnote n="15">¹⁵ Available at http://www.jmmc.fr/searchcal.</footnote>
the circumstellar component is

$$V_c = (1 - f_L)V_G + f_L V_L, \quad (4)$$

where $f_L$ is a weighting factor for the contribution of the pseudo-Lorentzian component, and $V_G$ and $V_L$ are the Hankel transforms of the Gaussian and pseudo-Lorentzian brightness distribution, respectively.

The best fit results in Section 3.4 were obtained by selecting the first quartile of the visibility data, on the basis of individual error estimates. Table 2 provides an overview of all the fitting results of the PIONIER visibilities to show the effect of the cutoff level on the best-fit values. The cutoff level in the cumulative distribution of the error estimates is stepwise loosened from 25% up to 100%, that is, the first quartile and all available data points, respectively. The error estimates for individual data points are derived by the PIONIER reduction pipeline from the internal dispersion on a short timescale. However, experience shows that other errors, not captured by the estimate, become increasingly prevalent as the observing conditions degrade. Our choice of the first quartile is a compromise between a severe selection that discards valid information and a loose selection that includes corrupted data, as is reflected in the increase of the reduced $\chi^2$ in Table 2.

## ORCID iDs

Tomas Stolker https://orcid.org/0000-0002-5823-3072

Mike Sitko https://orcid.org/0000-0003-1799-1755

Myriam Benisty https://orcid.org/0000-0002-7695-7605

Carsten Dominik https://orcid.org/0000-0002-3393-2459

Sebastian Perez https://orcid.org/0000-0003-2953-755X

Julien Milli https://orcid.org/0000-0001-9325-2511

Stefan Kraus https://orcid.org/0000-0001-6017-8773

## References

Andrews, S. M., Wilner, D. J., Espaillat, C., et al. 2011, ApJ, 732, 42

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Avenhaus, H., Quanz, S. P., Schmid, H. M., et al. 2014, ApJ, 781, 87

Bans, A., & Königl, A. 2012, ApJ, 758, 100

Benisty, M., Stolker, T., Pohl, A., et al. 2017, A&A, 597, A42

Beuzit, J.-L., Feldt, M., Dohlen, K., et al. 2008, Proc. SPIE, 7014, 701418

Bonneau, D., Clausse, J.-M., Delfosse, X., et al. 2006, A&A, 456, 789

Bonneau, D., Delfosse, X., Mourard, D., et al. 2011, A&A, 535, A53

Brown, J. M., Blake, G. A., Dullemond, C. P., et al. 2007, ApJL, 664, L107

Canovas, H., Ménard, F., de Boer, J., et al. 2015, A&A, 582, L7

Canovas, H., Rodenhuis, M., Jeffers, S. V., Min, M., & Keller, C. U. 2011, A&A, 531, A102

Carmona, A., Pinte, C., Thi, W. F., et al. 2014, A&A, 567, A51

Debes, J. H., Poteet, C. A., Jang-Condell, H., et al. 2017, ApJ, 835, 205

Demidova, T. V., Grinin, V. P., & Sotnikova, N. Y. 2013, AstL, 39, 26

Dohlen, K., Langlois, M., Saisse, M., et al. 2008, Proc. SPIE, 7014, 70143L

Dong, R. 2015, ApJ, 810, 6

Dong, R., & Fung, J. 2017, ApJ, 835, 38

Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJL, 809, L5

Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503

Dullemond, C. P., & Dominik, C. 2004, A&A, 417, 159

Dullemond, C. P., van den Ancker, M. E., Acke, B., & van Boekel, R. 2003, ApJL, 594, L47

Edwards, S. 2009, in AIP Conf. Proc. 1094, Stellar Systems, and the Sun, ed. E. Stempels (Melville, NY: AIP), 29

Eiroa, C., Oudmaijer, R. D., Davies, J. K., et al. 2002, A&A, 384, 1038

Ellerbroek, L. E., Podio, L., Dougados, C., et al. 2014, A&A, 563, A87

Espaillat, C., Furlan, E., D'Alessio, P., et al. 2011, ApJ, 728, 49

Espaillat, C., Muzerolle, J., Najita, J., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 497

Fairlamb, J. R., Oudmaijer, R. D., Mendigutía, I., Ilee, J. D., & van den Ancker, M. E. 2015, MNRAS, 453, 976

Fedele, D., van den Ancker, M. E., Acke, B., et al. 2008, A&A, 491, 809

Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2017, ApJ, 835, 230

Fung, J., & Dong, R. 2015, ApJL, 815, L21

Fusco, T., Rousset, G., Sauvage, J.-F., et al. 2006, OExpress, 14, 7515

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2016, A&A, 595, A2

Garufi, A., Meeus, G., Benisty, M., et al. 2017, A&A, 603, A21

Garufi, A., Quanz, S. P., Avenhaus, H., et al. 2013, A&A, 560, A105

Grady, C. A., Schneider, G., Sitko, M. L., et al. 2009, ApJ, 699, 1822

Grinin, V. P., & Tambovseva, L. V. 2003, in Conf. Proc. Towards Other Earths: DARWIN/TPF and the Search for Extrasolar Terrestrial Planets, ed. M. Fridlund, T. Henning, & H. Lacoste (Noordwijk: ESA), 429

Grinin, V. P., The, P. S., de Winter, D., et al. 1994, A&A, 292, 165

Haubois, X., Mella, G., Duvert, G., et al. 2016, Proc. SPIE, 9907, 99073V

Heap, S. R., Lindler, D. J., Lanz, T. M., et al. 2000, ApJ, 539, 435

Hinkley, S., Oppenheimer, B. R., Soummer, R., et al. 2009, ApJ, 701, 804

Hubrig, S., Stelzer, B., Schöller, M., et al. 2009, A&A, 502, 283

Hutchison, M. A., Laibe, G., & Maddison, S. T. 2016, MNRAS, 463, 2725

Jordi, K., Grebel, E. K., & Ammon, K. 2006, A&A, 460, 339

Kama, M., Pinilla, P., & Heays, A. N. 2016, A&A, 593, L20

Kenyon, S. J., & Bromley, B. C. 2004, AJ, 127, 513

Lai, D. 2014, MNRAS, 440, 3532

Langlois, M., Dohlen, K., Vigan, A., et al. 2014, Proc. SPIE, 9147, 91471R

Lazareff, B., Berger, J.-P., Kluska, J., et al. 2017, A&A, 599, A85

Le Bouquin, J.-B., Berger, J.-P., Lazareff, B., et al. 2011, A&A, 535, A67

Lyo, A.-R., Ohashi, N., Qi, C., Wilner, D. J., & Su, Y.-N. 2011, AJ, 142, 151

Maire, A.-L., Bonnefoy, M., Ginski, C., et al. 2016, A&A, 587, A56

Maire, A.-L., Stolker, T., Messina, S., et al. 2017, A&A, 601, A134

Marino, S., Perez, S., & Casassus, S. 2015, ApJL, 798, L44

Min, M., Canovas, H., Mulders, G. D., & Keller, C. U. 2012, A&A, 537, A75

Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, A&A, 497, 155

Müller, A., van den Ancker, M. E., Launhardt, R., et al. 2011, A&A, 530, A85

Murakawa, K. 2010, A&A, 518, A63

Muto, T., Grady, C. A., Hashimoto, J., et al. 2012, ApJL, 748, L22

Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 411, 1104

Pauls, T. A., Young, J. S., Cotton, W. D., & Monnier, J. D. 2005, PASP, 117, 1255

Pavlov, A., Möller-Nilsson, O., Feldt, M., et al. 2008, Proc. SPIE, 7019, 701939

Pérez, L. M., Isella, A., Carpenter, J. M., & Chandler, C. J. 2014, ApJL, 783, L13

Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&ARv, 21, 69

Rosenfeld, K. A., Qi, C., Andrews, S. M., et al. 2012, ApJ, 757, 129

Sitko, M. L., Carpenter, W. J., Kimes, R. L., et al. 2008, ApJ, 678, 1070

Sitko, M. L., Day, A. N., Kimes, R. L., et al. 2012, ApJ, 745, 29

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Stolker, T., Dominik, C., Avenhaus, H., et al. 2016, A&A, 595, A113

Tambovseva, L. V., Grinin, V. P., & Weigelt, G. 2006, A&A, 448, 633

Turner, N. J., Benisty, M., Dullemond, C. P., & Hirose, S. 2014, ApJ, 780, 42

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2016, A&A, 585, A58

van der Marel, N., van Dishoeck, E. F., Bruderer, S., Pérez, L., & Isella, A. 2015, A&A, 579, A106

Waters, L. B. F. M., & Waelkens, C. 1998, ARA&A, 36, 233

Wisniewski, J. P., Clampin, M., Grady, C. A., et al. 2008, ApJ, 682, 548

Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103

Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321