# Dust modeling of the combined ALMA and SPHERE datasets of HD 163296. Is HD 163296 really a Meeus group II disk?

G. A. Muro-Arena, C. Dominik, L. B. F. M. Waters, M. Min, L. Klarmann, C. Ginski, A. Isella, M. Benisty, A. Pohl, A. Garufi, et al.

► To cite this version:

G. A. Muro-Arena, C. Dominik, L. B. F. M. Waters, M. Min, L. Klarmann, et al. Dust modeling of the combined ALMA and SPHERE datasets of HD 163296. Is HD 163296 really a Meeus group II disk?. Astronomy and Astrophysics - A&A, 2018, 614, 10.1051/0004-6361/201732299. insu-03635034

HAL Id: insu-03635034

https://insu.hal.science/insu-03635034

Submitted on 8 Apr 2022

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Dust modeling of the combined ALMA and SPHERE datasets of HD 163296

## Is HD 163296 really a Meeus group II disk?

G. A. Muro-Arena<fnref n="1" />, C. Dominik<fnref n="1" />, L. B. F. M. Waters<fnref n="1" />,<fnref n="2" />, M. Min<fnref n="1" />,<fnref n="2" />, L. Klarmann<fnref n="1" />, C. Ginski<fnref n="1" />,<fnref n="3" />,
A. Isella<fnref n="4" />, M. Benisty<fnref n="5" />,<fnref n="6" />, A. Pohl<fnref n="7" />, A. Garufi<fnref n="8" />, J. Hagelberg<fnref n="6" />, M. Langlois<fnref n="9" />,<fnref n="10" />, F. Menard<fnref n="6" />, C. Pinte<fnref n="6" />,
E. Sezestre<fnref n="6" />, G. van der Plas<fnref n="11" />,<fnref n="12" />, M. Villenave<fnref n="6" />, A. Delboulbe<fnref n="6" />, Y. Magnard<fnref n="6" />, O. Möller-Nilsson<fnref n="7" />,
J. Pragt<fnref n="13" />, P. Rabou<fnref n="6" />, and R. Roelfsema<fnref n="13" />

Received 15 November 2017 / Accepted 6 February 2018

## ABSTRACT

**Context.** Multiwavelength observations are indispensable in studying disk geometry and dust evolution processes in protoplanetary disks.

**Aims.** We aim to construct a three-dimensional model of HD 163296 that is capable of reproducing simultaneously new observations of the disk surface in scattered light with the SPHERE instrument and thermal emission continuum observations of the disk midplane with ALMA. We want to determine why the spectral energy distribution of HD 163296 is intermediary between the otherwise well-separated group I and group II Herbig stars.

**Methods.** The disk was modeled using the Monte Carlo radiative transfer code *MCMax3D*. The radial dust surface density profile was modeled after the ALMA observations, while the polarized scattered light observations were used to constrain the inclination of the inner disk component and turbulence and grain growth in the outer disk.

**Results.** While three rings are observed in the disk midplane in millimeter thermal emission at ~80, 124, and 200 AU, only the innermost of these is observed in polarized scattered light, indicating a lack of small dust grains on the surface of the outer disk. We provide two models that are capable of explaining this difference. The first model uses increased settling in the outer disk as a mechanism to bring the small dust grains on the surface of the disk closer to the midplane and into the shadow cast by the first ring. The second model uses depletion of the smallest dust grains in the outer disk as a mechanism for decreasing the optical depth at optical and near-infrared wavelengths. In the region outside the fragmentation-dominated regime, such depletion is expected from state-of-the-art dust evolution models. We studied the effect of creating an artificial inner cavity in our models, and conclude that HD 163296 might be a precursor to typical group I sources.

**Key words.** protoplanetary disks – scattering – techniques: polarimetric – techniques: interferometric – stars: individual: HD 163296

## 1. Introduction

Protoplanetary disks are the sites of planet formation, of which the disks around Herbig Ae/Be are the best-studied examples owing to their brightness and proximity to Earth. The structure and morphology of protoplanetary disks around Herbig Ae/Be stars is as of now not well understood, although various models

have been proposed over the years to link their global geometry to the characteristics of their spectral energy distributions (SEDs). The classification system initially proposed by Meeus et al. (2001) first divided Herbig Ae/Be SEDs into two groups based on the far-infrared (FIR) excess of these sources: those whose SEDs can be fitted by a single power-law component (group II), and those that require an additional broad component

<footnote n="1">¹ Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
e-mail: g.a.muroarena@uva.nl</footnote>

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

<footnote n="3">³ Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="4">⁴ Department of Physics and Astronomy, Rice University, 6100 Main Street, Houston, TX 77005, USA</footnote>

<footnote n="5">⁵ Unidad Mixta Internacional Franco-Chilena de Astronomía, CNRS/INSU UMI 3386 and Departamento de Astronomía,
Universidad de Chile, Casilla 36-D, Santiago, Chile</footnote>

<footnote n="6">⁶ Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="7">⁷ Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="8">⁸ Universidad Autónoma de Madrid, Dpto. Física Teórica, Módulo 15, Facultad de ciencia, Campus de Cantoblanco,
28049 Madrid, Spain</footnote>

<footnote n="9">⁹ CRAL, UMR 5574, CNRS, Université Lyon 1, 9 avenue Charles André, 69561 Saint Genis Laval Cedex, France</footnote>

<footnote n="10">¹⁰ Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, 13388 Marseille, France</footnote>

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

<footnote n="12">¹² Millenium Nucleus Protoplanetary Disks in ALMA Early Science, Universidad de Chile, Casilla 36-D, Santiago, Chile</footnote>

<footnote n="13">¹³ NOVA Optical Infrared Instrumentation Group, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands</footnote>
peaking in the mid- to far-infrared (group I). An additional subdi-
vision was proposed for group I based on the presence (group Ia)
or absence (group Ib) of the 10 µm silicate feature, whereas it
was observed that this silicate feature was present in all group II
sources.

The difference in the SEDs between group I and group II
sources was for a long time attributed to a flaring versus flat
geometry (Dullemond & Dominik 2004), and the larger FIR
excess of group I sources was explained by the absorbed and
reprocessed light in the flaring outer regions of these disks.
Maaskant et al. (2013), on the other hand, pointed at the inner
wall of a flaring outer disk, directly illuminated by the cen-
tral star from the presence of a large gap or inner cavity in
the disk, as the origin of the cold component in group I SEDs.
These authors attributed the near-infrared (NIR) excess to an
optically thin and hot component at the interior of this cav-
ity. We now believe group II sources to be either flat of very
compact objects because they can rarely be resolved in scat-
tered light (Garufi et al. 2017). Although classically classified
as a group II source (Meeus et al. 2001), HD 163296 is some-
times referred to as an intermediate-type source in the literature
(e.g., Dominik et al. 2003), with a FIR excess greater than typi-
cal group II sources yet lower than group I sources. If the view
we hold of group II sources as compact or flat, self-shadowed
disks is correct, then direct observations with the Hubble Space
Telescope (HST; Grady et al. 2000; Wisniewski et al. 2008)
point toward a possible misclassification of the disk around
HD 163296. Dust in the disk was detected in scattered light
by both studies out to radial distances of at least ~450 AU,
which is direct evidence that the disk is not compact and the
outermost regions are flaring and directly illuminated by the
central star. Recently, in a taxonomical study of a sample of
17 group I and group II sources, Garufi et al. (2017) have pro-
posed that HD 163296 might be a precursor to typical group I
disks.

HD 163296 is a ~5.1<sup>0.3</sup><sub>0.8</sub> Myr old (Alecian et al. 2013) Herbig
Ae star located at 122<sup>+17</sup><sub>-13</sub> pc from the Sun (van den Ancker
et al. 1997). The mass of this object is 2.23<sup>0.22</sup><sub>0.07</sub> M<sub>☉</sub> (Alecian
et al. 2013), and it has a luminosity of ~33.1<sup>6.7</sup><sub>2.3</sub> L<sub>☉</sub>
(Alecian et al. 2013). With an effective temperature of
~9200 ± 300 K (Folsom et al. 2012), it is classified as a spectral
type A1Ve star (van den Ancker et al. 1998). The massive disk
surrounding this star was first resolved interferometrically 20
years ago by Mannings & Sargent (1997) using the Owens
Valley Radio Observatory (OVRO) millimeter-wave array; sizes
were derived for both the continuum and gas line emission of
this disk at 1.3 mm of ~100 and 300 AU, respectively. Recent,
more sensitive observations with ALMA have since shown
the gaseous component of the disk extends to distances of at
least 500 AU (Isella et al. 2016; de Gregorio-Monsalvo et al.
2013), while the dust continuum at millimeter wavelengths is
only detected out to a radial distance of ~240 AU from the
central star. The continuum images reveal further structure in
the dust: a bright inner disk component is observed within
the inner ~0.5 arcsec (60 AU) of the image with two bright
rings at ~0.66 arcsec (80 AU) and ~1.0 arcsec (122 AU) and
a third much fainter ring at ~1.6 arcsec (195 AU) (Isella et al.
2016).

In addition to the HST observations, ground-based observa-
tions in polarized scattered light in the *J* band with the Gemini
Planet Imager (Monnier et al. 2017) and in *H* and *Ks* bands with
VLT/NACO (Garufi et al. 2014) have also shown emission in the
NIR out to distances of ~80 AU around HD 163296, that is, much

more compact than the scattered light seen in HST images yet
providing additional evidence that the disk around HD 163296
is not flat and self-shadowed, or at least not at all radii. This
is also smaller than the disk observed in the millimeter con-
tinuum, which poses an interesting question as to the nature of
this and possibly other so-called intermediate sources that might
not fit easily with current interpretations of the group classifi-
cation scheme. The drastically different nature of the observed
disk structure between NIR and millimeter-wavelength obser-
vations shows that neither stands on its own as a complete,
comprehensive picture of this disk. With this in consideration,
we attempted to produce a single, physically consistent, three-
dimensional model capable of simultaneously reproducing the
dust observations in both NIR and millimeter wavelengths.

In this paper we present new polarized scattered light
observations of HD 163296 obtained with SPHERE/IRDIS and
describe the data reduction in Sect. 2. Results of these observa-
tions are presented in Sect. 2.2. In Sect. 3 we detail the procedure
followed to model the surface density of the disk based on the
ALMA observations (Isella et al. 2016), and propose two differ-
ent models capable of explaining some of the features seen in the
SPHERE/IRDIS data. The results are discussed in Sect. 4, par-
ticularly in the context of disk geometry and its role in the Meeus
group classification scheme. Finally, we present our conclusions
and summarize our findings in Sect. 5.

## 2. Observations

### 2.1. Data reduction

HD 163296 was observed as part of the SPHERE (Beuzit et al. 2008) guaranteed time program with the infrared dual band imager and spectrograph (IRDIS; Dohlen et al. 2008) on May 26, 2016. Observations were carried out in dual-polarization imaging mode (DPI; Langlois et al. 2014) in the *J* and *H* bands. The sky was clear during the observations, but seeing conditions were variable with an average seeing of 0.9 arcsec in the optical and seeing spikes up to 1.2 arcsec, especially during the *H*-band sequence. In DPI mode IRDIS uses polarizers inserted into the beam to measure two orthogonal polarization directions simultaneously (ordinary and extraordinary beam). A half-wave plate (HWP) is used to rotate the polarization of the incoming light beam such that the two linear Stokes vector components *Q* and *U* can be measured. The HWP is also used to correct in first order for instrumental polarization by rotating in steps of 22.5°. Such HWP rotation sequence allows the measurement of the four linear polarization components per polarization cycle: *Q*<sup>+</sup>, *Q*<sup>-</sup>, *U*<sup>+</sup> and *U*<sup>-</sup>. The HWP orientation is chosen such that between *Q*<sup>+</sup> and *Q*<sup>-</sup> only the sign of the astrophysical signal of the Stokes *Q* component changes, while the instrumental polarization downstream from the HWP remains unchanged. Thus by subtracting *Q*<sup>-</sup> from *Q*<sup>+</sup> instrumental polarization is canceled out while the astrophysical signal is preserved. For *J* band we recorded a total of 15 polarization cycles with an individual exposure time (DIT) of 16 s and 5 exposures per HWP position (NDIT). This amounted to a total integration time of 80 min. In the *H* band we recorded three polarization cycles: one cycle with a DIT of 8 s and an NDIT of 16 and two cycles with a DIT of 16 s and NDIT of 8. The change in exposure time was performed to adjust to the variable seeing conditions. This amounted to a total integration time of 25.6 min in the *H* band.

The data reduction was performed following largely the
approach layed out in Ginski et al. (2016). We give here a short
summary and highlight differences to the mentioned procedure.
<figure><img src="image_2.png" /><figcaption>Fig. 1. H-band (top) and J-band (bottom) DPI observations of HD 163296 with SPHERE/IRDIS. Left column: $Q_\phi$ images, with the right column showing $U_\phi$. A single asymmetrical ring is observed in scattered light at ~0.6 arcsec from the center of the disk in both $Q_\phi$ images, slightly offset from center in the southwest direction. An additional bright inner component shows up in both bands, the center of which is obscured by the coronograph. Some residual signal shows up in both $U_\phi$ images, mostly constrained to the inner 0.2–0.3 arcsec. The white dashes indicate the position of the major and minor axes labeled a and b, respectively.</figcaption></figure>

After dark, bad pixel, and flat-field correction, we performed a single difference step by subtracting ordinary and extraordinary beams from each other, thereby creating the previously mentioned $Q^+$, $Q^-$, $U^+$, and $U^-$ polarization components. We then performed the so-called double difference step (Kuhn et al. 2001; Apai et al. 2004) by subtracting $Q^-$ from $Q^+$ (and analog for $U$) to generate the Stokes Q and U components for each polarization cycle. Since conditions were variable and thus the tip-tilt correction of the AO system was not always stable, we ensured a good alignment of the frames before the double difference step. For this purpose we used several bright background point sources to which we fitted Moffat functions to measure their positions between frames. After the double difference step, we combined the Q and U frames for all polarimetric cycles. We then performed an instrumental polarization correction as described by Canovas et al. (2011) to subtract the residual instrumental

polarization from the final Stokes Q and U frames. We assumed for this purpose that the stellar light is completely unpolarized. We thus measured the stellar signal at the bright adaptive optics correction cut-off radius in Stokes Q, U, and the intensity image and then subtracted a scaled version of the intensity image from the Q and U images to suppress remaining signal, which should be instrumental in nature. After this final step was performed we computed the polarized intensity (PI) image from the Stokes Q and U components as follows:

$$PI = \sqrt{Q^2 + U^2}. \qquad (1)$$

Since the disk around HD 163296 is not strongly inclined, it is reasonable to assume that the polarized light that we receive from the system is overwhelmingly created in single scattering events (Canovas et al. 2015). We can thus assume that the angle
of linear polarization shows an azimuthal alignment with respect to the central star. It is then convenient to use the polar Stokes components $Q_\phi$ and $U_\phi$ as defined by Schmid et al. (2006), that is,

$$Q_{\phi} = Q \times \cos 2\phi + U \times \sin 2\phi, \qquad (2)$$

$$U_{\phi} = Q \times \sin 2\phi - U \times \cos 2\phi, \qquad (3)$$

wherein $\phi$ is the azimuth with respect to the center of the star. $Q_\phi$ contains all the azimuthal polarized flux as positive signal, whereas $U_\phi$ contains all polarized flux with a 45° offset from azimuthal. Thus $U_\phi$ is in our case expected to contain very little or no signal and can be used as a convenient noise estimator. We show the final reduced $Q_\phi$ and $U_\phi$ images for both bands in Fig. 1.

## 2.2. Observational results

The $Q_\phi$ images of HD 163296 in both J and H bands in Fig. 1 show the same ring structure at approximately 0.6 arcsec from the center of the image with only the north half of the ring clearly visible in the J-band image because of a lower contrast. The ring appears to have an offset from the center of the image in the southwest direction, likely a combined product of the inclination of the disk and the location of the scattering surface at a certain height from the midplane, which is indicative of a geometrically thick disk (Ginski et al. 2016; de Boer et al. 2016). It can be inferred from the direction of the offset that the northeast side of the disk corresponds to its front side.

The polarized flux density in the ring appears asymmetric in both J and H bands. First, there is the usual asymmetry in the direction of the minor axis, likely to be caused by the different scattering angles for the front part of the disk (northeast, brighter because of forward scattering) and back part (southwest, darker), in combination with the scattering phase function. In addition there is also an asymmetry along the major axis in both bands with the brightest side of the ring corresponding to the northwest side and roughly but not exactly coincidental with the major axis. The asymmetry along the major axis cannot be explained by geometric effects and the scattering phase function. Without other factors at play, and assuming the dust is distributed homogeneously in the ring, the measured flux density is expected to be symmetric along the major axis (i.e., mirror symmetric w.r.t. the minor axis). The observed asymmetry is thus probably the result of either an uneven dust distribution in the surface of the disk, an azimuthally variable scale height, or the product of a shadow cast on the surface on one side of the ring by an inner disk component.

Elliptical annuli were fitted to both J- and H-band $Q_\phi$ images to obtain the semimajor axis $a$, position angle PA and offset of the ring ($\Delta x_0$, $\Delta y_0$). The fitting procedure consisted of a Monte Carlo routine that generated $10^6$ elliptical annuli with all four parameters randomized within restricted ranges. The width of each annulus was set to four pixels, and the inclination was set at $i = 46^\circ$ for both H- and J-band images. The total flux was averaged over each aperture, and the best-fit parameters were obtained from the elliptical ring aperture, which maximizes the averaged flux. We chose to fix the inclination at this given value, corresponding to the inclination measured from the ALMA image, as most attempts to measure the inclination from the SPHERE images yielded inclinations larger than $i = 52^\circ$ for the J-band image and similarly values that are too large for the H band. This resulted in apertures that did not properly fit the far side of the ring, most likely due to the variable width and brightness of the ring in both images.

The errors for each parameter were obtained from the H-band image, as it has a higher signal-to-noise ratio (S/N). The noise was measured from the $U_\phi$ image as the standard deviation $\sigma_{U_\phi}$ of the flux inside the best-fitting aperture. The uncertainty the noise introduces in the flux measured from the $Q_\phi$ image was calculated as $\Delta F = \sigma_{U_\phi} / \sqrt{N}$, where N is the number of pixels in the aperture. We then selected all those apertures from the initial $10^6$ whose average flux is larger than $F_{max} - \Delta F$, and obtained the errors in the four fitted parameters as the standard deviation of the parameters for this subset of apertures.

A semimajor axis of $a = 0.63 \pm 0.045$ arcsec was measured in the H band, which corresponds to a radius of 77 ± 5.5 AU (122 AU per arcsecond, assuming a distance of 122 pc), along with a position angle of PA = 134.8° ± 11°. The ring observed in the SPHERE images thus matches the size and orientation of the first ring observed in the Band 6 ALMA continuum observations, which has a radius of ~0.66 arcsec or 80 AU and a position angle of PA = 134°. The large error measured for the position angle is again likely caused by the uneven width and brightness around the ring, which introduces uncertainty in our fitting procedure by allowing a broad range of position angles to yield similar fluxes inside the apertures. An offset of 105 ± 45 mas was measured in the H band in the southwest direction (PA = 215° ± 22°), in good agreement with the offset measured by Monnier et al. (2017) and corresponding to a height of the scattering surface of 12.7 ± 5.5 AU above the disk midplane. A similar offset of 92 mas was measured in the J-band image, which gives a scattering surface height of 11.2 AU above the midplane. The semimajor axis measured in the J band was $a = 0.61$ arcsec or 74 AU. The position angle measured in this band is PA = 137.5°.

No flux with a S/N above 3 is measured in either the J or H band beyond ~100 AU, in contrast to the second and third rings detected in the ALMA continuum observations out to distances of ~240 AU; but this confirms the nondetections at larger radii by Monnier et al. (2017); Garufi et al. (2014), who also observe a single ring with a similar asymmetry to it. Interestingly, no emission is observed out to larger radii of ~500 AU, whereas HST did detect scattered light in the disk out to these radii (Grady et al. 2000; Wisniewski et al. 2008). This is likely due both to the higher sensitivity of HST to faint, extended features and the fact that this instrument observes total intensity as opposed to the polarized light of the SPHERE images. A bright inner disk component can be seen in the H-band $Q_\phi$ image and less prominently in the J-band $Q_\phi$ image, in both cases partially obscured by the coronograph.

The top panel of Fig. 2 shows the radial profiles of both $Q_\phi$ images along the major axis of the ring. Each data point corresponds to the flux density averaged over all the pixels contained in an aperture with a radius equal to the image resolution, for each band, and then normalized to the maximum value of the profile. The error bars correspond to the standard deviation of the flux density at each aperture over the square root of the number of point spread functions (PSFs) contained in it. The panel shows the brightness asymmetry in the ring along this axis in both photometric bands; positive radii correspond to the northwest direction and negative radii to southeast. The H-band profile shows the northwest side of the ring is brighter than the southeast side by about a factor 3. The southeast side of the ring is not visible in the J-band image, but the position of the northwest side of the ring matches that of the H-band image. The ring width is well resolved in both bands, with the FWHM of the ring approximately four times the image resolution at either end of the major axis. The profile of the inner disk is marred by the
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Top panel: radial polarized intensity profile through the center of the ring along the major axis for H-band (black) and J-band (red) Q<sub>φ</sub> images. Middle panel: radial polarized intensity profile through the star in the direction of the major axis for both bands. Positive radii correspond to the northwest direction, and negative radii to the southeast direction. The vertical dashed lines indicate the approximate location of the ring at 77 AU. The gray region at the center of the figure corresponds to the location of the coronograph. Bottom panel: azimuthal intensity profile of the ring in the H band (black) and J band (red) at a distance of 77 and 74 AU from the center of the disk, respectively, centered at the center of each ring. The position angle PA is measured east from north. All profiles are normalized to the maximum intensity.</figcaption></figure>

small circular (instrumental) artifacts that can be seen close to the center of the image in both bands in Fig. 1: a single one in the *J* band, and three in the *H* band in a triangular formation around the coronograph, two of which are along the ring major axis. The middle panel of Fig. 2 shows the radial profile through the position of the star in the direction of the ring major axis, which provides a better view of the profile of the inner disk. In both images it can be observed that the inner disk appears brighter on the southeast side of the coronograph.

The bottom panel of Fig. 2 shows the azimuthal profile of each Q<sub>φ</sub> image at the location of the ring as a function of the position angle, measured east from north. The data again correspond to the pixel-averaged flux density inside apertures with the same radii as described above, normalized to the maximum of the profile. Error bars correspond to the standard deviation of the flux density at each aperture over the square root of the number of PSFs contained in the apertures. It is observed that the maxima both occur at a position angle of ~350°, about 35° from the semimajor axis toward the front side of the disk. The profiles are in good agreement with each other within the error bars except at position angles between 80° and 180°, where the ring vanishes in the noise in the second quadrant of the *J*-band image whereas the ring still appears clearly visible at this location in the *H*-band image.

### 3. Modeling the dusty disk around HD 163296

The observations with the SPHERE instrument differ greatly from the disk observed in ALMA band 6 images (Isella et al. 2016). In contrast to the single ring at ~77 AU observed in the

J and H bands, the ALMA continuum and line images show a disk that extends at least out to around 500 AU in gas and 240 AU in large dust grains (de Gregorio-Monsalvo et al. 2013; Isella et al. 2016), and presents both a large inner disk extending out to ~50 AU and three outer rings. In this section we present the results of two models that are capable of explaining the differences observed between the midplane thermal emission observations and the scattered light observations of the disk surface. We do this in the context of global disk geometry and the how it ties into the Meeus group classification scheme.

#### 3.1. The initial model

The modeling of the disk was carried out using the three-dimensional Monte Carlo radiative transfer code MCMax3D (Min et al. 2009). The modeling of two independent zones is necessary to fit the disk SED; a first zone composed of an unresolved inner disk (UID) to account for the NIR excess in the SED, and a second zone to model the structure resolved by the polarized scattered light and ALMA continuum observations. The UID has been spatially resolved by NIR interferometric data (Benisty et al. 2010; Lazareff et al. 2017) and modeled independently by several studies (Lazareff et al. 2017; Renard et al. 2010; Tannirkulam et al. 2008). The focus of this work is on the outer disk resolved by the SPHERE and ALMA data, however, and the model of the UID presented here is not tested by the current observations.

The MCMax3D code allows the user to define disk zones with different pressure scale heights, surface densities, grain size distributions, gas-to-dust ratios, and turbulences. The code uses a power-law pressure scale height profile of the form

$$H_g(r) = H_0 \left(\frac{r}{r_0}\right)^p, \qquad (4)$$

where $H_0$ is the scale height of the disk at a radius $r_0$ and $p$ is a positive exponent. A simple power-law surface density profile is defined for the UID, that is,

$$\Sigma(r) \propto r^{-\epsilon}, \qquad (5)$$

where $\epsilon$ is a positive number. The resolved disk is modeled using the following power-law profile with an exponential taper-off:

$$\Sigma(r) \propto r^{-\epsilon} \exp\left(-\left(\frac{r}{r_c}\right)^{2-\gamma}\right), \qquad (6)$$

where $r_c$ is the taper-off radius, and $2-\gamma$ is a positive exponent, while the grain size distribution is modeled as

$$n(a) da \sim a^{-p_a} da, \qquad (7)$$

for a particle size $a$ between a minimum and maximum particle sizes $a_{min}$ and $a_{max}$, respectively, and $p_a$ a positive exponent.

A two-zone model from the DIANA project database (Woitke et al. 2016), obtained from fitting the SED with a genetic fitting algorithm, was used as a starting point for our modeling (P. Woitke, priv. comm.<fnref n="1" />). The parameters characterizing this model can be found in Table 1. A gas-to-dust ratio of a 100 is assumed for the whole disk, while the turbulence is characterized by $\alpha \sim 2.4 \times 10^{-3}$, which is in agreement with the upper limit of $\alpha = 3 \times 10^{-3}$ reported by Flaherty et al. (2017). Both zones are

<footnote n="1">¹ http://www.univie.ac.at/diana/index.php/user/login</footnote>
<table><caption>**Table 1.** Initial *MCMax3D* model parameters.</caption><thead><tr><th></th><th>Zone 1</th><th>Zone 2</th></tr></thead><tbody><tr><td>R<sub>in</sub> (AU)</td><td>4.04 × 10<sup>-1</sup></td><td>1.70</td></tr><tr><td>R<sub>out</sub> (AU)</td><td>1.07</td><td>240</td></tr><tr><td>M<sub>d</sub> (M<sub>☉</sub>)</td><td>6.6 × 10<sup>-9</sup></td><td>4.23 × 10<sup>-4</sup></td></tr><tr><td>ε</td><td>7.35 × 10<sup>-1</sup></td><td>1.035</td></tr><tr><td>γ</td><td>—</td><td>0.585</td></tr><tr><td>r<sub>c</sub> (AU)</td><td>—</td><td>120</td></tr><tr><td>H<sub>0</sub> (AU)</td><td>4.83 × 10<sup>-2</sup></td><td>5.64</td></tr><tr><td>r<sub>0</sub> (AU)</td><td>1</td><td>100</td></tr><tr><td>p</td><td>1.37</td><td>1.15</td></tr><tr><td>a<sub>min</sub> (μm)</td><td>5 × 10<sup>-2</sup></td><td>5 × 10<sup>-2</sup></td></tr><tr><td>a<sub>max</sub> (μm)</td><td>9.75 × 10<sup>-1</sup></td><td>1 × 10<sup>4</sup></td></tr><tr><td>p<sub>a</sub></td><td>3.795</td><td>3.795</td></tr><tr><td>g2d</td><td>100</td><td>100</td></tr><tr><td>α</td><td>2.4 × 10<sup>-3</sup></td><td>2.4 × 10<sup>-3</sup></td></tr></tbody></table>

<table><caption>**Table 2.** Inner and outer radii of zones 2a, 2b, and 2c in our model.</caption><thead><tr><th></th><th>Zone 2a</th><th>Zone 2b</th><th>Zone 2c</th></tr></thead><tbody><tr><td>R<sub>in</sub> (AU)</td><td>1.70</td><td>63</td><td>105</td></tr><tr><td>R<sub>out</sub> (AU)</td><td>63</td><td>105</td><td>240</td></tr></tbody></table>

assumed to have the same inclination and position angle. For this model we used the values measured from the first ring in the ALMA image, that is, PA = 132° and i = 46°.

A graphical representation of the model zones is shown in Fig. 3. The UID is shown in red (zone 1), with the resolved disk shown in shades of blue (zone 2). Zone 2 was additionally subdivided into zones 2a, 2b, and 2c. These zones account for the RID seen in the ALMA image (zone 2a), the first ring seen in both ALMA and SPHERE images (zone 2b), and the two outer rings seen only in the ALMA image (modeled jointly in zone 2c). The inner and outer radii for these zones are listed in Table 2. Below we discuss the different steps used to arrive at a full model of the HD 163296 disk. Starting from the initial model, we model the surface density in order to match the ALMA observations (Sect. 3.2). While this is successful, it turns out that in order to match the SPHERE observations as well, modifications are required in zone 2. These modifications are described in Sects. 3.3–3.5. All initial parameters for zones 2a, 2b, and 2c in our models correspond to the parameters of zone 2 listed in Table 1. For consistency, we refer to the different zones in the disk instead of disk components from now onward when referring to the models.

## 3.2. Surface density

The SPHERE/IRDIS images are an example of how polarized scattered light images do not constitute the best tracer for the dust mass in the disk. Compared to the millimeter emission seen in the ALMA image, the disk appears to be much reduced in size with all of the structure outside the ring at 0.63 arcsec (77 AU) completely invisible under the noise level. This is likely due in part to the low percentage of the total intensity that is polarized. The emission in the RID is affected by leftover instrumental artifacts, which limit the fidelity with which we can reconstruct the surface density of the disk in this region.

Zone 1

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

Zone 2

**Fig. 3.** Graphical representation of the 2-zone *MCMax3D* model. Zone 1, in red, corresponds to a UID component with an outer radius of 1 AU. Zone 2 corresponds to the outer disk resolved by the ALMA observations. Zone 2 has been subdivided into zones 2a, 2b, and 2c, corresponding to the resolved inner disk (RID, zone 2a), inner ring (zone 2b), and second and third rings (zone 2c). This subdivision was created for the purpose of locally altering some model parameters. For the purpose of modeling the surface density of the disk, however, zone 2 is treated as a single region. The inner and outer radii for each zone are given in Tables 1 and 2.

The thermal emission we observe at 1.3 mm in the ALMA continuum image delivers a more robust view of the bulk of the dust mass. It is not subject to shadowing effects and it is not affected by the instrumental artifacts caused by the coronograph and the bright point source at the center of the image, although it does depend on the temperature structure of the disk and its optical depth, and the radial extent of the disk is limited by radial drift. The disk is extremely bright at millimeter wavelengths and due to the high sensitivity of ALMA, the S/N is also very large. With these considerations, we can use the spatial brightness distribution of the disk as observed in the ALMA image to model the radial surface density profile of the dust to high fidelity, limited only by the spatial resolution of the image and dependent on the parameters used for the initial model.

This modeling of the surface density is carried out iteratively since the RID is optically thick even at millimeter wavelengths, which translates to nonlinear variations of the surface brightness of the disk as a function of the surface density. Also, changing the surface density influences the temperature in the disk, requiring iteration. The process followed can be summarized in the following steps, assuming a surface density distribution that is perfectly azimuthally symmetric:

- The radiative transfer model for the initial set of parameters (described in Sect. 3.1) is obtained with *MCMax3D* for a given surface density profile, starting with the power-law surface density characterized by the parameters in Table 1.

- The obtained model is ray traced to obtain an image of the disk at 1.3 mm.

- The ray-traced image is convolved with a Gaussian kernel to simulate the ALMA reconstructed image. The kernel shape and orientation are chosen to match the characteristics of the
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** *Left panel:* radial profiles of the ALMA image (black), ray-traced model image (blue), and ray-traced model image after convolution with a 2-dimensional Gaussian kernel with dimensions of 0.22 arcsec by 0.15 arcsec and a PA of -88° (green). Each point was obtained by integrating the flux density in elliptical annuli of varying radii and normalizing first by the aperture area, and then by the maximum of the ALMA image profile. The error bars for the ALMA profile were obtained as the rms of the image divided by the number of synthesized beams contained in each annuli. *Right panel:* side-by-side comparison of the ALMA image (left) and the ray-traced image of one of our models (right) after Gaussian blurring.</figcaption></figure>

synthesized beam of the ALMA observations: beam dimensions are 0.22 × 0.15 arcsec with a PA of -88° measured north to east.

- The radial profile of the ALMA image is obtained by integrating the disk flux in elliptical annuli of increasing size centered on the center of the image, and normalized by the annulus area. The elliptical annuli have a PA of 132°, eccentricity of 0.719 (corresponding to an inclination of 46°) and a width of 5 pixels.

- The radial profile of the convolved model image is obtained in the same manner.

- The radial surface density profile used as input for the model is then scaled by the point-by-point ratio between the ALMA brightness profile and the model brightness profile. In this way mass is added to the disk at radii at which the model is too faint, and mass is subtracted at radii at which the model is too bright (i.e., gaps).

- This scaled surface density is then used as the starting point for the next iteration of the model.

By iterating the steps above to model the surface density of zone 2, a radial intensity profile can be obtained that converges to that of the ALMA image profile, as seen in the left panel of Fig. 4. This plot shows an excellent agreement between the radial intensity profile of the ALMA observation and our model (using the parameters in Table 1) after the iteration of the surface density. We have called this model M0. A side-by-side comparison of the ALMA observations and the convolved image of the model can be seen in the right panel of Fig. 4, which shows excellent overall agreement between the ALMA data and model M0. Both the data and model appear to show an RID with a slightly lower PA than the rest of the disk, but this is an effect of the orientation of the synthesized beam and Gaussian kernel that was used to convolve the ray-traced model image. This method is only applicable to azimuthally symmetric sources, and the results are also of course model dependent. A different grain size distribution, turbulence, chemical composition, or scale height of one of the disk components, for example, produces a different resulting surface density profile. In addition, it can only be used as far as the larger grains extend out in the disk. The ALMA observations provide us with the ability of measuring the surface density out to radii of ~240 AU, as no thermal

emission is detected in the continuum image above the noise level beyond that point. The model we present here is thus truncated at 240 AU, despite the observational evidence of both gas and small dust grains on the disk surface extending out to over 500 AU.

The width and depth of the gaps observed in the disk structure cannot be uniquely determined from the observations, as these are limited by resolution. Isella et al. (2016) have studied the gap width-depth degeneracy in detail and offer a more comprehensive view on the extent of this degeneracy based on the assumption that the surface density of the disk can be approximated as a power law with square gaps in it. The method described here offers more degrees of freedom in that the overall dust surface density of the disk is not constrained to a radial power law, and in that it allows for a smooth and continuous radial surface density profile, which is probably a more natural approximation to the real surface density of the disk than that of sharp gap edges. This is supported by the radial profile observed in the *J*- and *H*-band images in Fig. 2, which have an angular resolution ~4 times higher than that of the ALMA image and in which no sharp edges are detected for the ring at 77 AU.

While this model is successful at reproducing the data at milimeter wavelengths, it fails in the NIR. This is seen clearly in Fig. 5, which shows a comparison in both *H* and *J* bands between the SPHERE data (top row), and the ray-traced images of M0 (second row). The figure also shows the ray-traced images of models M1 and M2, to be discussed below in Sects. 3.3–3.5 (third and bottom rows, respectively). Reproducing the structure seen in the SPHERE/IRDIS images requires the localized modification of some parameters to account for three major differences between these images and the ray-traced images of model M0 at NIR wavelengths. In first place, the *H*- and *J*-band images of M0 show clearly the presence of the second ALMA ring, which is not visible in our NIR data. Secondly, the produced images of model M0 show a more extended RID than observed in the data. And lastly, since model M0 is azimuthally symmetric, it fails to reproduce the observed asymmetry along the major axis of the first ring, described in Sect. 2.2. The modification of parameters required to reproduce the SPHERE/IRDIS images, and the two models resulting from these modifications, are discussed in Sect. 3.3 through 3.5.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Side-by-side comparison of the *Q*<sub>Φ</sub> scattered light images for both *H* band (**left**) and *J* band (**right**) for our observations (**top**), model M0 (second row), model M1 (third row), and model M2 (**bottom**). Model M0 consists of the model characterized by the parameters in Table 1. Model M1 consists of a UID inclined by 3° with respect to the outer disk, characterized by an inclination *i* = 45° and position angle PA = 136°, along with a small grain depletion characterized by a minimum grain size *a*<sub>min</sub> = 3 µm in zone 2c. Model M2 consists of inclined UID and RID components inclined by 1.3° with respect to the rest of the disk. This inclination is characterized by *i* = 45° and PA = 133.2°. Increased settling in the outer disk of this model is characterized by *α* = 1 × 10<sup>-5</sup> in zone 2c.</figcaption></figure>

### 3.3. Absence of outer rings in polarized scattered light: Increased settling or grain growth in the outer disk?

Compared to the ALMA image, the polarized scattered light images in both the *H* and *J* band show a curious absence of emission beyond the ring at ~0.63 arcsec. While the rings seen in the ALMA image are bright and emission is measured out to almost three times the radial extent of the scattered light disk, no structure is seen in polarized scattered light in the outer disk

with a S/N above ~2.5. Considering scattered light has also been detected out to 500 AU with roll subtraction (Schneider & Silverstone 2003; Fraquelli et al. 2004) by Grady et al. (2000) using HST STIS and by Wisniewski et al. (2008) with HST ACS, we are left to wonder the reason for the absence of polarized intensity observed in our SPHERE images beyond the first ring. While the sensitivity of our ground-based observations is lower than that of HST, model M0 shows that after modeling the surface density of the outer disk the second ring should still be
clearly visible at 1.25 and 1.6 µm even considering only ~30% of scattered light is polarized by dust. Considering the properties detailed for our model in Sect. 3.1, and assuming no abrupt radial variations in these parameters, a sharp decrease in surface density would be required at the location of the second gap, around ~0.85 arcsec or 105 AU, to obscure the outer rings. However this is not observed in the results obtained from the modeling of the surface density based on the ALMA continuum image.

We can infer from this that some of the input parameters for our model are in fact not radially constant. Since we know there is a large dust mass at radii larger than 100 AU, it follows that there are no, or few, small grains present in the surface of the disk directly illuminated by the star in this region. This can be explained by the small grains settling closer to the midplane of the disk, in the shadow cast by the first ring, or by a local depletion of the smallest grains starting at around 100 ± 15 AU from the center of the disk.

Our next model (M1) attempts to explain the absence of outer rings in the *J* and *H*-band images using a modified grain size distribution in zone 2c. The optical depth needs to be decreased to hide the outer rings in this zone while maintaining turbulence constant throughout the disk. We do this by depleting this region of the smallest dust grains, which are responsible for the optical depth in the NIR. We implement this in our first model by increasing the minimum particle size while keeping *a*<sub>max</sub> and *a*<sub>pow</sub> unchanged. With this modified setup we then execute the iterative procedure described in Sect. 3.2. We find that depleting this zone of grains up to a size of 3 µm is necessary to decrease the optical depth in this region enough to obscure the second and third rings in the *H*-band and *J*-band *Q*<sub>Φ</sub> images. We find a depletion of all small grains up to 2.5 µm still allows for the second ring to show up clearly in the ray-traced images of our models. In Sect. 4.2 we present a dust evolution model that shows it is possible to obtain a rapid decrease in the surface density of small dust grains, and we show this partial depletion is strong enough to almost completely obscure the presence of a second ring at 125 AU.

Our last model (M2) uses increased settling of the dust grains at radii larger than 105 AU. This is achieved by decreasing the turbulence in zone 2c to bring the small grains on the surface of the disk closer to the midplane and into the shadow cast by the first ring (zone 2b), while keeping the grain size distribution in this zone equal to the rest of zone 2. With this modified setup we then execute the iterative procedure described in Sect. 3.2, just as we did for model M1. We find that the *α* parameter needs to be lowered to ~1×10⁻⁵ to achieve this effect. While we model this as a sharp decrease in the turbulence at 105 AU for simplicity, a more gradual decrease across the gap separating zones 2b and 2c would yield similar results. The global distribution of the small dust grains in these two models is summarized in Fig 6. The top panel shows a flaring disk with a localized depletion of small grains in a ring around the disk, starting in our case somewhere in between zones 2b and 2c, and extending out indefinitely, possibly all the way to the outer edge of the disk, corresponding to the geometry of model M1. The bottom panel shows the geometry of the disk for model M2, consisting of a flared structure that extends out to the first ring at 77 AU; a decrease in the height of the scattering surface is caused by localized increased settling in zone 2c. While both models are on their own sufficient to explain the absence of scattered light observed in the HD 163296 disk beyond ~100 AU, they are not completely independent. Even a partial depletion of small grains in the outer disk would cause the height of the τ = 1 surface to decrease in some degree, resulting in a geometry that is effectively very similar to that

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Simplified graphical representation (not to scale) of the two possible scenarios that might explain the lack of observed polarized scattered light in the SPHERE/IRDIS images in the outer disk. Top panel: small grain depletion scenario is shown, in which small grains are present throughout the inner regions of a flared disk but are largely depleted on the surface of the disk at the location of the second ring. The larger grains are shown on the disk midplane, while the locations of the gaps are indicated by the white vertical bands. Bottom panel: increased settling scenario shows a double flaring structure for the dust and therefore for the scattering surface. A flaring outer edge with a small number of small grains remaining is necessary in both scenarios to account for the scattered light observed in HST images.</figcaption></figure>

of the first model, albeit with a different mechanism at work behind it. Additionally, increased settling and depletion of small grains in this zone are not mutually exclusive, and the absence of the outer dust rings in the scattered light images might well be caused by a combination of these two effects, or in combination with a decreased pressure scale height or dust chemistry. The obtained surface density profiles for these models are shown in Fig. 7, along with the total mass of zone 2c obtained for each model in Table 3. We can see the effect of modifying the grain size distribution and turbulence in zone 2c: by increasing the minimum particle size in the outer disk (M1), we are effectively increasing the mass of this component of the disk by ~50% during the modeling of the surface density either as a product of a decreased temperature or optical depth. Since the smallest dust grains provide only a small fraction of the optical depth at millimeter wavelengths (a few percent), we can attribute this to a decreased temperature caused by the diminished absorption of stellar light in this zone. Similarly, the mass in this zone is approximately doubled if instead the turbulence is decreased to $\alpha = 1 \times 10^{-5}$, also as a product of a lowered midplane temperature (M2).

### 3.4. The resolved inner disk

One of the most evident differences between the observations at millimeter and NIR wavelengths is the apparent size of the RID (zone 2a). The RID appears about 50% larger in the ALMA image than in the *H*-band polarized scattered light image and even larger compared to the *J*-band image (see Fig. 5).
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Surface density profiles for models M0 (dashed), M1 (dotted), and M2 (solid). Zones 2a, 2b, and 2c are indicated in the figure as progressively lighter shades for the background color.</figcaption></figure>

<table><caption>**Table 3.** Total mass for each zone for all three models, in solar masses, after the iterative modeling of the surface density.</caption><thead><tr><th></th><th>M0</th><th>M1</th><th>M2</th></tr></thead><tbody><tr><td>Zone 1</td><td>6.44 × 10<sup>-9</sup></td><td>6.44 × 10<sup>-9</sup></td><td>6.44 × 10<sup>-9</sup></td></tr><tr><td>Zone 2a</td><td>1.3 × 10<sup>-4</sup></td><td>2.3 × 10<sup>-3</sup></td><td>2.25 × 10<sup>-3</sup></td></tr><tr><td>Zone 2b</td><td>1 × 10<sup>-4</sup></td><td>1.1 × 10<sup>-4</sup></td><td>9.8 × 10<sup>-5</sup></td></tr><tr><td>Zone 2c</td><td>1.4 × 10<sup>-4</sup></td><td>2.1 × 10<sup>-4</sup></td><td>2.8 × 10<sup>-4</sup></td></tr><tr><td>Total</td><td>3.7 × 10<sup>-4</sup></td><td>2.6 × 10<sup>-3</sup></td><td>2.6 × 10<sup>-3</sup></td></tr></tbody></table>

The scale height of zone 2a was modified to reduce the apparent size of the RID in models M1 and M2. The scale height was fixed at the inner radius of zone 2a, but the flaring exponent *p* was reduced from *p* = 1.15 to *p* = 1.07, making this zone close to wedge-shaped. Decreasing the scale height of this zone in this manner has the desired effect of reducing the apparent size of the RID, but it also modifies its temperature structure. By decreasing the scale height we are also reducing the amount of stellar light being absorbed by the dust in this component of the disk, which makes the temperature in this region drop. The drop in temperature translates directly to a drop in the surface brightness, and therefore to a larger dust mass obtained in this zone as a result of the iterative modeling of the surface density described in Sect. 3.2. The surface density in this zone is found to increase by over an order of magnitude in models M1 and M2 as a product of the reduced scale height, yielding a total dust mass for the whole disk of ~2.6 × 10⁻³ M☉, which seems unusually high; see Fig. 7, which shows the surface density profiles in zone 2 for models M0, M1, and M2. Accurately estimating the mass of zone 2a is a complicated matter, considering this zone is optically thick in our models. However, because of the large dust masses required in this zone for models M1 and M2, it seems unlikely that a reduced scale height of this zone is solely responsible for its smaller apparent size in the scattered light images.

The total mass of the disk and its individual zones are detailed for each of the three models in Table 3. Interestingly, the surface density in zone 2b remains approximately equal in all three models, showing that the temperature structure on this zone remains independent of the modified scale height in zone 2a.

### 3.5. *Brightness asymmetry: Evidence of self-shadowing?*

Another predominant feature of the scattered light images is the asymmetry observed in the flux density of the ring (zone 2b) along the major axis, which cannot be explained by either the scattering phase function or the angle dependence of the polarization efficiency. Such an asymmetry could be caused by an accumulation of dust on the north side of the disk, or by the shadow of an inclined inner disk on the south side of the ring, similar to what has been proposed for HD 142527 (Marino et al. 2015), HD 135344B (Stolker et al. 2016), and HD 100453 (Benisty et al. 2017). Considering this asymmetry is not observed in the ALMA image, we favor the shadowing scenario.

Taking into account the structure of HD 163296 two possibilities are available, in which either zone 1 (UID) is inclined by a few degrees with respect to the rest of the disk, or in which both zones 1 (UID) and 2a (RID) are misaligned by similar degree. We find that an inclination of 3° with respect to the rest of the disk for zone 1, characterized by an inclination of 45° and a position angle of 136°, produces a brightness asymmetry qualitatively very similar to the asymmetry observed in the scattered light images. This inclination however also has the effect of casting a shadow on the south side of zone 2a (RID). If both zones 1 and 2a of the disk are slightly inclined, however, this shadow is avoided. We find that if both components are inclined, an inclination relative to the outer disk of just ~1.3°, characterized by an inclination of 45° and a PA of 133.2°, is required to reproduce to a high degree the asymmetry seen in the scattered light images.

Figure 5 shows the *H*- and *J*-band *Q*<sub>φ</sub> images of the SPHERE data (top) and the obtained ray-traced *Q*<sub>φ</sub> images of models M0, M1, and M2 (second, third, and bottom row, respectively). These images were convolved with Gaussian kernels with a FWHM of 41 and 50 mas to simulate the resolution of the *J*- and *H*-band images, respectively, measured from the PSF of the flux images. The models were scaled down such that the integrated flux density in an elliptical annulus containing the ring, with a width of four times the image resolution, equals that of the observations. The noise was measured from the *H*-band and *J*-band *U*<sub>φ</sub> images in apertures with a radius of 4 and 3.4 pixels, respectively, assuming they contain no remaining signal, and added artificially to the model images pixel by pixel using a Gaussian random number generator to simulate the noise. In model M1 zone 1 was inclined with respect to the rest of the disk to show the produced asymmetry. In model M2 both zones 1 and 2a were inclined as described to show the resulting asymmetry. The choice of inclining different zones in these two models was motivated by purely illustrative purposes, and has no larger significance. The purpose is to show that similar results can be obtained by inclining only zone 1, or both zones 1 and 2a. Either choice of inclined components would work for either model.

The resulting azimuthal profiles of models M1 and M2 are shown for both the *H* and *J* bands in Fig. 8 along with the SPHERE profiles for comparison. Both models show a similarly good agreement with the *J*-band profile in the third and fourth quadrant of the disk with some deviations in the first and second quadrant. However, the S/N in this band is much lower than in the *H* band and therefore should be interpreted with care. For the *H* band, we can see an excellent agreement with the data for both models, which manage to reproduce the observed asymmetry to a high degree except at position angles between ~50° and 130°. The reason for this difference is currently unclear. It might be due to a local disturbance in the scale height of the disk, either in the ring or in the inner regions causing a localized shadow on the ring. It is also possible that the dust properties or
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Azimuthal profiles of the ring for the SPHERE data (black), model M1 (blue), and model M2 (green) for both *H* band (top) and *J* band (bottom). The profiles were measured at a distance of 77 AU and 74 AU for *H*- and *J*-band images, respectively.</figcaption></figure>

<table><caption>**Table 4.** Parameters modified in models M1 and M2 compared to the parameters used in model M0.</caption><thead><tr><th></th><th>M0</th><th>M1</th><th>M2</th></tr></thead><tbody><tr><td>α in zone 2c</td><td>2.4 × 10<sup>-3</sup></td><td>2.4 × 10<sup>-3</sup></td><td>1 × 10<sup>-5</sup></td></tr><tr><td>a<sub>min</sub> (μm) in zone 2c</td><td>5 × 10<sup>-2</sup></td><td>3</td><td>5 × 10<sup>-2</sup></td></tr><tr><td>H<sub>0</sub> (AU) in zone 2a</td><td>5.64</td><td>0.073</td><td>0.073</td></tr><tr><td>r<sub>0</sub> (AU) in zone 2a</td><td>100</td><td>1.7</td><td>1.7</td></tr><tr><td>p in zone 2a</td><td>1.15</td><td>1.07</td><td>1.07</td></tr><tr><td>i zone 1</td><td>46°</td><td>45°</td><td>45°</td></tr><tr><td>PA zone 1</td><td>132°</td><td>136°</td><td>133.2°</td></tr><tr><td>i zone 2a</td><td>46°</td><td>46°</td><td>46°</td></tr><tr><td>PA zone 2a</td><td>132°</td><td>132°</td><td>133.2°</td></tr></tbody></table>

the abundance of small grains are different at this location, but a structure like this would be expected to wash out within a few orbital timescales.

A summary of the parameters used for each model is found in Table 4.

# 4. Discussion

## 4.1. Disk geometry

While typically classified as a group II source, HD 163296 is often referred to in the literature (e.g., Dominik et al. 2003) as an intermediate-type source. It is tempting therefore to assume it might be an example of a type of protoplanetary disk that constitutes an evolutionary link between typical group I and group II sources. While this is indeed a possibility, the true nature of the evolution of Herbig disks might be different.

Two separate evolutionary processes should be distinguished in noncompact disks: the formation of an extended gap or cavity in the outer disk leading to group I-like SEDs, and the deflaring or vertical collapse of the disk, leading to flat group II sources, assuming all disks start out as both continuous and flaring. It follows from observational evidence that the disk around HD 163296 is neither compact nor flat or – at least not

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Spectral energy distribution of HD 163296. The points and their corresponding error bars correspond to the photometry of the source. The total SED of model M2 is shown in black, along with its individual components: stellar spectrum (red), zone 1 (blue), zone 2a (green, solid), zone 2b (green, dashed), and zone 2c (green, dotted).</figcaption></figure>

completely – self-shadowed, as it can be observed in scattered light images as far out as ~500 AU. Nevertheless, observations of the source so far have failed to resolve a cavity in the outer disk, either in millimeter dust emission or scattered light, although cavities have been detected in $^{13}$CO and C$^{18}$O with ALMA out to a few tens of AU (Isella et al. 2016), which are comparable to the cavity sizes of group I sources. We may tentatively conclude then, until further observational evidence, that HD 163296 constitutes part of a subcategory of protoplanetary disks that has undergone neither of the processes that characterize either group I or group II sources. It follows that the disk either constitutes part of a separate group of protoplanetary disks or that the disk is in an evolutionary stage anterior to either group I or group II sources.

With this in mind we computed the SED for each of the components of model M2, which can be seen in Fig. 9, by ray tracing our model at 200 different wavelengths between 0.1 µm and 7.5 mm. The SED obtained for model M1 is qualitatively very similar. The figure shows the photometry of HD 163296, along with the total SED of the model and the SED contribution of each of its individual components. The SEDs of models M0 and M1 are very similar, since the initial model used fits the SED and the modifications implemented in the models do not have a big effect on it. We can see the bulk of the NIR flux comes from the thermal emission of zone 1, and only a small fraction is attributed to zones 2a and 2b in scattered light. In the FIR the thermal emission from zone 2a dominates and has a weaker, but still comparable, contribution from the thermal emission of zone 2b, that is smaller only by about a factor 2 up to ~100 µm and about equal from there up to millimeter wavelengths. This seems to suggest zones 2a and 2b combined might be playing a role in the SED of HD 163296 that is similar to that of the bright inner wall of the outer disk in group I sources.

To analyze this further, we generated four different models identical to M2 in all but the total dust mass in zones 1 and 2a (UID and RID). The dust in these zones of the first, second, and third of these models was depleted by a factor $1 \times 10^2$, $1 \times 10^3$, and $1 \times 10^4$, respectively. In the fourth model the dust in zones 1 and 2a was depleted completely. The resulting SED for all these models is shown in Fig. 10. The figure illustrates the effect in the SED of progressively removing the inner disk
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Spectral energy distribution of HD 163296. The points and their corresponding error bars correspond to the photometry of the source. The total SED of model M2 is shown as a thick solid red line. The solid, dashed, dot-dashed, and dotted black lines show the SEDs obtained by depleting zones 1 and 2a of the model by a factor 10², 10³, 10⁴, and complete depletion, respectively.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Reconstructed grain size distribution for a 2.5 Myr disk with an initial mass of 0.2 $M_\odot$, initial dust-to-gas ratio of 0.02, turbulence $\alpha = 2.4 \times 10^{-3}$, and fragmentation velocity $v_{frag} = 3$ m s⁻¹.</figcaption></figure>

components, effectively creating a large cavity in the disk inside the first ring. Not only does the NIR component of the total disk SEDs decrease as dust is removed from both the inner disks, but the thermal emission of the outer disk increases as the temperature in this zone rises, which is a consequence of progressively increasing the amount of stellar radiation received by this zone. Lost NIR excess notwithstanding, our models show that artificially creating a dust cavity in the disk, even a partially depleted cavity, can produce a larger FIR excess reminiscent of typical group I sources. This leads us to believe that the lack of such a cavity might be the largest discriminator between HD 163296 and other group I sources.

## 4.2. Absence of multiple rings in scattered light

Three bright rings at ~75, 125, and 200 AU from the center of the disk can be seen in the thermal emission from the mid-plane, while our polarized scattered light images with IRDIS in the J and H bands show a clear lack of the two outermost

of these rings. Two models were presented in Sect. 3.3 capable of reproducing the observations, explaining the lack of small dust grains on the disk surface beyond ~100 AU by either an increased settling (model M2) or a depletion of small grains (model M1) in the outer disk (zone 2c). While we modeled this depletion by completely removing all dust grains smaller than 3 µm in zone 2c, we aim to show that a strong but partial depletion of the smallest dust grains in the outer disk is not only capable of reproducing the absence of the outer rings in scattered light, but is also consistent with dust evolution models.

As a proof of concept we present the results from a dust evolution simulation obtained using the dust evolution code twopoppy (Birnstiel et al. 2012). Our goal was not to produce an accurate model of the surface density of the disk, but rather to show there may be regions in a disk where a depletion of small grains can be achieved naturally as a product of dust evolution. The code was used to evolve a disk with an initial mass of 0.2 $M_\odot$. A turbulence parameter of $2.4 \times 10^{-3}$ was used for consistency with our MCMax3D model, along with a low fragmentation velocity of 3 m/s and an initial dust to gas ratio of 0.02. An initial gas surface density characterized by the parameter $\gamma = 0.8$ and a characteristic radius $r_c = 165$ AU was used, that is,

$$ \Sigma(r) \propto r^{-\gamma} \exp\left(-\left(\frac{r}{r_c}\right)^{2-\gamma}\right). \qquad (8) $$

A high initial disk mass and dust-to-gas ratio are required to obtain a total dust mass of $5 \times 10^{-4} M_\odot$, in the lower end of the $5-17 \times 10^{-4} M_\odot$ dust mass range found in the literature (Natta et al. 2004; Tannirkulam et al. 2008; Mannings & Sargent 1997; Isella et al. 2007), and a total disk mass of $\sim 7 \times 10^{-2} M_\odot$ after 2.5 Myr of evolution. Owing to the rate at which the disk loses dust mass it is not possible to achieve a dust mass of $5 \times 10^{-4} M_\odot$ after 5 Myr of evolution without a much higher initial disk mass or, alternatively, a much higher initial dust-to-gas ratio. A way of achieving the larger dust mass observed in the disk around HD 163296 without requiring a very high initial dust-to-gas ratio or initial mass is with a disk that has continued accreting gas and dust from the circumstellar environment for a couple million years after its formation, replenishing the dust content of the disk and allowing this way for the large dust mass we observe today. If this were the case for HD 163296, this would mean the disk we are observing is effectively younger than the 5 Myr estimated for the system. Evolving such a disk for 2.5 Myr produces the grain size distribution seen in Fig. 11, obtained using the grain size distribution reconstruction module of twopoppy (Birnstiel et al. 2015).

The figure shows a grain size distribution that is limited by fragmentation in the inner ~100 AU of the disk, and mostly by grain growth in the outer disk, except for a small transition zone in between that is radial drift-limited. The small grains in the inner region grow quickly to the fragmentation limit and are then replenished by the collision and fragmentation of the larger grains. In the outer disk dust grains have grown quickly to the drift limit and drifted inward, leaving behind a population of small grains with a density that is too low to grow efficiently to larger sizes. This means the smaller grains in this region can still grow slowly, but are not efficiently replenished as they do not reach the fragmentation limit. In the intermediate transition region, dust is still able to grow to millimeter sizes, but the drift limit is reached before they can fragment, making the replenishing of small grains very inefficient.
<figure><img src="image_13.png" /><figcaption>Fig. 12. Comparison between the H-band Q<sub>φ</sub> SPHERE image (left) and model M2 (right). A second ring was artificially added to the image of the model with a semimajor axis of 125 AU, a thickness of 24 AU, and a surface brightness 46 times lower than that of the first ring.</figcaption></figure>

Since the outer disk is growth-limited, we observe no large amount of large grains beyond 100–200 AU. This is qualitatively consistent with the ALMA observations, which show a sharp decrease in optical depth around ~125 AU between the second and third rings. The remaining intermediate size grains beyond ~125 AU could also provide the low optical depth responsible for the faint thermal emission seen from the third ring in the ALMA image. The low abundance of small grains remaining in the outer disk in this model might be enough to explain the emission seen in the HST scattered light images, while at the same time this abundance is low enough to explain the absence of the second and third ALMA rings in the polarized scattered light images.

To determine if this depletion of small grains in the outer regions of the disk is enough to explain the nondetection of the second ring in the polarized scattered light images, we estimate the radial optical depth contrast between the location of the first ring, at 77 AU, and the second ring at 125 AU. The optical depth per unit length in the radial direction at a distance *r* from the star in the *J* and *H* bands is given by the sum of the optical depths contributed by each grain size population,

$$ \frac{\Delta\tau(r)_{J,H}}{\Delta l} = \sum_{a_i} \kappa(a_i)_{J,H} \rho(r, a_i), \qquad (9) $$

where $κ(a_i)_{H,J}$ is the opacity of dust grains of size $a_i$ at the *J* or *H* band, obtained from *MCMax3D*, $ρ(r, a_i)$ the density of grains of size $a_i$ on the surface of the disk at a distance *r* from the star, and $\Delta l$ a radial distance. Since both rings are unresolved (or barely resolved) in the ALMA image, which has a resolution of 0.2 arcsec or ~24 AU, we cannot estimate their thickness. However, the FWHM for the first ring was measured from the *H*-band $Q_φ$ image and found to be approximately 24 AU. We shall assume for simplicity this is the thickness of both rings, although it is more likely an upper limit for the width of the second ring. This assumption allows us to compare the optical depth contrast, and therefore the surface brightness contrast, between the two rings. We approximate $ρ(r, a_i)$ by

$$ \rho(r, a) \sim \frac{\Sigma(r, a)}{H(r)}, \qquad (10) $$

where $\Sigma(r, a)$ is the surface density of particles of size $a_i$ at a distance *r* from the center of the star obtained from the grain size reconstruction of our dust evolution model, and $H(r)$ the *MCMax3D* pressure scale height at a distance *r* from the

center of the disk. The optical depth contrast between the first and second ring can then be expressed as

$$ \frac{\Delta\tau(R_1)_{J,H}}{\Delta\tau(R_2)_{J,H}} \sim \frac{H(R2) \sum_{a_i} \kappa(a_i)_{J,H} \Sigma(R1, a)}{H(R1) \sum_{a_i} \kappa(a_i)_{J,H} \Sigma(R2, a)}, \qquad (11) $$

where $R_1$ and $R_2$ are the radii of the first and second rings, respectively. We integrated the product of $κ(a_i)\Sigma(r, a_i)$ over all grain sizes at both radii, at a wavelength of 1.6 µm, and found this value to drop by a factor 14 from the first to the second ring. Accounting also for the spread of dust grains over a larger vertical extent (first term on the right-hand-side of Eq. (11) using the scale height of our *MCMax3D* model (see Table 1), we find the optical depth at 125 AU to drop by an additional 1.8 factor relative to the optical depth of the first ring, yielding a total optical depth contrast of about 25.

Assuming the radial optical depth at the location of the first ring is $τ \sim 1$, this amounts to a relative intensity about 16 times lower for the second ring as a consequence of the drop in surface density and the increasing scale height. Taking into account the $r^{-2}$ dependence of the stellar radiation field as well, the total intensity of scattered light we could expect for the second ring in this model can be as low as 46 times fainter than the intensity of the first ring.

To test whether this contrast in intensity is high enough to account for the absence of the second ring in the scattered light images, we manually added a second ring to the *H*-band image of our model M2. This was carried out by measuring the average flux of the first ring at each position angle, scaling it down by a factor 46, and adding it to the image at a distance of 125 AU from the center of the first ring to recreate a similar offset. The flux was measured at each position angle in circular apertures of radius *r* = 8 pixels, and then added in circles of the same radius at the expected location of the second ring, and at the corresponding position angle. This allowed us to simulate a second ring with a similar width to that of the first ring and a similar asymmetry. The resulting image is shown in the right side of Fig. 12. We see that, at ~2% of the intensity of the first ring, the artificially added second ring is almost completely invisible at the noise level of our polarized scattered light *H*-band image. Only the brightest side of the ring shows up, at a very low S/N, at the northwest side of the image.

Even though this scenario can successfully explain the nondetection of the outer rings in the scattered light images, a depletion of the smallest grains in the outer disk to this degree is not strictly necessary to produce this effect. As discussed in
Sect. 3.3 an increased settling in zone 2c, characterized by a small turbulence $\alpha$ parameter, can also explain why the two outer rings observed in the ALMA image are not visible in polarized scattered light. With the small dust grains responsible for the opacity in the NIR settling closer to the midplane, the scattering surface can be pulled down and into the shadow cast by the inner ring. The radiative transfer models obtained using an $\alpha = 1 \times 10^{-5}$ in zone 2c (see Fig. 5, bottom row) show that this increased settling can be successfully used to explain the absence of the outer rings in the scattered light images.

While these two simplified cases are enough, on their own, to explain why we observe only the innermost ring in the scattered light images, a more complex model lying somewhere in between these two extremes is possible. A lower turbulence, coupled with a partial depletion of the smallest grains on the surface of the disk, can be expected to produce comparable results. An incomplete depletion of smaller grains in the outer disk is for us the favored scenario, since a small remaining number of small grains in the outer disk out to several hundred AU is required to explain the observations of scattered light out to ~500 AU observed with HST. Shadowing of the outer disk could also be caused by increasing the scale height in the region of the first ring. However, this would influence the SED and change the observed offset of the scattered light ring from the star.

## 5. Conclusions

We have presented in this paper two independent models that are capable of reproducing simultaneously the millimeter continuum observations of the midplane of HD 163296 with ALMA and polarized scattered light observations of the disk surface with SPHERE/IRDIS.

The observations of the midplane allowed us to obtain a surface density profile of the disk, for the parameters defined in Sect. 3.1, out to a truncation radius of 240 AU. This truncation radius marks the outer edge of the dust disk at 1.3 mm wavelength. The gas component of the disk, on the other hand, is seen to extend out to ~500 AU in both gas line observations and scattered light. An inner disk component and a single ring at 77 AU characterize the polarized scattered light images in both *J* and *H* bands. The ring presents an offset of 105 mas, corresponding to a height of 12 AU above the disk midplane in the *H* band and 11.5 AU in the *J* band. This observation of scattered light is further evidence of a flaring disk. However, a simple flaring disk model is unable to account for the nondetection of any signal beyond the second gap at 100 AU in the SPHERE images.

The nondetection in scattered light of the two outer rings observed in the ALMA dust continuum image is evidence of the absence of small dust grains on the disk surface in this zone. This points to either the smaller grains settling closer to the midplane, under the shadow cast by the first ring, or to a partial-to-complete depletion of the small grains in this region. We find with model M2 that a settling characterized by $\alpha \sim 1 \times 10^{-5}$ is required to settle the small grains responsible for the opacity at *J* and *H* bands into the shadow cast by the first ring. In contrast, if settling is not increased with respect to the rest of the disk, a total depletion of grains smaller than 3 µm can have the same effect (model M1). The two mechanisms, while different, are approximately equivalent in that the produced effect is the lowering of the scattering surface under the surface of the disk directly illuminated by the central star. A dust evolution model has been provided as proof of concept that shows the feasibility of the depletion

of small grains in the middle-disk, supporting the small grain depletion scenario.

The evidence of HD 163296 possessing a global flaring geometry points to its misclassification as a group II source, insofar as our current understanding of these sources and their geometry goes. The lack of a largely depleted inner cavity sets it apart from group I sources as well. With an intermediate-type SED, we find from both of our models that the bulk of the FIR excess comes in similar measure from the RID and the first ring. With the RID and UID in the way, there is no large inner wall directly illuminated by the star, which likely explains why the FIR region of its SED is lower than a typical group I source. At the same time, the gap between the RID and the first ring is reminiscent of the extended dust-depleted gaps or cavities commonly found in group I sources. Our models show depleting the inner regions of this disk brings the properties of this object much closer to a group I source. Whether the disk is currently in the process of forming a dust-depleted inner cavity and transitioning toward a group I stage or not, both our models and observational evidence point to HD 163296 having more in common with group I sources than with group II sources.

The brightness asymmetry observed in the IRDIS images is likely the product of three separate effects: the scattering phase function, angle-dependent polarization efficiency, and self shadowing. We provided one model that shows a small misalignment of ~1.3° of the RID and UID suffices to reproduce this asymmetry to a high degree in all but the southeast quadrant, and a second model with the UID misaligned by 3° that produces similar results.

We may conclude from this work that multiwavelength observations are fundamental for a thorough understanding and proper modeling of protoplanetary disks as a consequence of the pronounced differences that arise at different wavelengths caused by the complex structures of disks. Some of the conclusions presented on this paper would not have been possible without the analysis of the disk observations at multiple wavelengths.

*Acknowledgements.* SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF – Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (The Netherlands), ONERA (France), and ASTRON (The Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland), and NOVA (The Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012), and grant number 312430 for FP7 (2013–2016). G. M.-A. acknowledges funding from the Netherlands Organisation for Scientific Research (NWO) TOP-1 grant as part of the research program “Herbig Ae/Be stars, Rosetta stones for understanding the formation of planetary systems”, project number 614.001.552. M.B. acknowledges funding from ANR of France under contract number ANR-16-CE31-0013 (Planet Forming Disks).

## References

Alecian, E., Wade, G. A., Catala, C., et al. 2013, MNRAS, 429, 1001

Apai, D., Pascucci, I., Wang, H., et al. 2004, in *Star Formation at High Angular Resolution*, eds. M. G. Burton, R. Jayawardhana, & T. L. Bourke, IAU Symp., 221, 307

Benisty, M., Natta, A., Isella, A., et al. 2010, A&A, 511, A74

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

Beuzit, J.-L., Feldt, M., Dohlen, K., et al. 2008, in *Ground-based and Airborne Instrumentation for Astronomy II*, Proc. SPIE, 7014, 12

Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148

Birnstiel, T., Andrews, S. M., Pinilla, P., & Kama, M. 2015, ApJ, 813, L14
Canovas, H., Rodenhuis, M., Jeffers, S. V., Min, M., & Keller, C. U. 2011, A&A, 531, A102

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

de Boer, J., Salter, G., Benisty, M., et al. 2016, A&A, 595, A114

de Gregorio-Monsalvo, I., Ménard, F., Dent, W., et al. 2013, A&A, 557, A133

Dohlen, K., Langlois, M., Saisse, M., et al. 2008, in Ground-based and Airborne Instrumentation for Astronomy II, Proc. SPIE, 7014, 70143L

Dominik, C., Dullemond, C. P., Waters, L. B. F. M., & Walch, S. 2003, A&A, 398, 607

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

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150

Folsom, C. P., Bagnulo, S., Wade, G. A., et al. 2012, MNRAS, 422, 2072

Fraquelli, D. A., Schultz, A. B., Bushouse, H., Hart, H. M., & Vener, P. 2004, PASP, 116, 55

Garufi, A., Quanz, S. P., Schmid, H. M., et al. 2014, A&A, 568, A40

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

Ginski, C., Stolker, T., Pinilla, P., et al. 2016, A&A, 595, A112

Grady, C. A., Devine, D., Woodgate, B., et al. 2000, ApJ, 544, 895

Isella, A., Testi, L., Natta, A., et al. 2007, A&A, 469, 213

Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101

Kuhn, J. R., Potter, D., & Parise, B. 2001, ApJ, 553, L189

Langlois, M., Dohlen, K., Vigan, A., et al. 2014, in Ground-based and Airborne Instrumentation for Astronomy V, Proc. SPIE, 9147, 91471R

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

Maaskant, K. M., Honda, M., Waters, L. B. F. M., et al. 2013, A&A, 555, A64

Mannings, V., & Sargent, A. I. 1997, ApJ, 490, 792

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

Meeus, G., Waters, L. B. F. M., Bouwman, J., et al. 2001, A&A, 365, 476

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

Monnier, J. D., Harries, T. J., Aarnio, A., et al. 2017, ApJ, 838, 20

Natta, A., Testi, L., Neri, R., Shepherd, D. S., & Wilner, D. J. 2004, A&A, 416, 179

Renard, S., Malbet, F., Benisty, M., Thiébaut, E., & Berger, J.-P. 2010, A&A, 519, A26

Schmid, H. M., Joos, F., & Tschan, D. 2006, A&A, 452, 657

Schneider, G., & Silverstone, M. D. 2003, in High-Contrast Imaging for Exo-Planet Detection, ed. A. B. Schultz, Proc. SPIE, 4860, 1

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

Tannirkulam, A., Monnier, J. D., Harries, T. J., et al. 2008, ApJ, 689, 513

van den Ancker, M. E., de Winter, D., & Tjin A Djie, H. R. E. 1998, A&A, 330, 145

van den Ancker, M. E., The, P. S., Tjin A Djie, H. R. E., et al. 1997, A&A, 324, L33

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