© ESO 2018

1 Introduction

High-resolution and high-contrast imaging capabilities provided by the new generation of adaptive optics based instruments such as the Spectro-Polarimeter High-contrast Exoplanet REsearch (SPHERE) instrument (Beuzit et al. 2008) and the Gemini Planer Imager (GPI; Macintosh et al. 2014) open the path to directly image new protoplanetary disks in scattered light. Protoplanetary disks are optically thick in the optical and near-infrared (NIR), so that scattered light imaging probes micron-sized dust grains in the disk’s surface layer, while mm observations trace larger, mm-sized grains located close to the disk’s mid-plane. Transition disks (TD) with gas and dust gaps (e.g. Brown et al. 2007; Merín et al. 2010; van der Marel et al. 2016; Ansdell et al. 2016; van Boekel et al. 2017) are particularly interesting since they might correspond to the stage where planet forming processes are active, thus we expect to see signposts of planet–disk interactions. The Atacama Large Millimeter Array (ALMA) now provides sufficient sensitivity and resolution at sub-mm wavelengths to detect these gaps (Fedele et al. 2017) and to estimate the mass of protoplanetary disks around YSOs (Ansdell et al. 2016; Trapman et al. 2017), although the bulk mass (H₂ and He) is not directly observable. In recent observational studies of transition disks, giant gaps and cavities have also been directly imaged in scattered light (e.g. Thalmann et al. 2010, 2015; Hashimoto et al. 2012; Avenhaus et al. 2014; Follette et al. 2015; Ohta et al. 2016; Stolker et al. 2016; Benisty et al. 2017; Pohl et al. 2017; van Boekel et al. 2017. In addition, the detection of non-axisymmetric disk features that could be related to the presence of planets is fundamental in order to improve our current understanding of the disk evolution and the planet formation process. However, the small angular separations involved pose great observational challenges.

RY Lup has been classified as a T Tauri G-type star showing type III variability located in the Upper Centaurus Lupus (UCL) association with an estimated age of 10–20 Myr (e.g. Manset et al. 2009; Pecaut et al. 2012). However, a colder effective temperatureof 5000 K and a best fit spectral type of K2 have been estimated by Biazzo et al. (2017) in the most recent literature based on X-shooter spectroscopy. Its distance has been estimated by recent measurements using Gaia (first GAIA data release, Gaia Collaboration 2016), which provided a distance estimate of 151 pc ± 1 pc. RY Lup is a TD, as confirmed by Ansdell et al. (2016) using the ALMA data, but was not previously identified assuch from its SED (Manset et al. 2009). The spectral energy distribution of RY Lup shows an infrared excess, as well as a modest UV excess (Evans et al. 1982; Gahm et al. 1989). Its strong 10 μm silicate emission feature, seen in its IRS spectrum (Kessler-Silacci et al. 2006), washes out the mid-IR dip in its broad-band fluxes.

The photometric variability of the star has been revealed by Woods using Harvard plates in 1921. The star was classified as a typical RW Aur star (a variable eruptive T Tauri-type star) on the basis of its light curve. From a long-term study of the optical photometry of the star, Gahm et al. (1989) reported significant photometric variations (up to 1.5 mag at visible wavelengths) on timescales of hours, a constant variability with a 3.75 day period over 50 yr, and a decline in the mean and maximum photometric magnitudes over decades. This photometric behaviour of RY Lup is typical of type III variability defined by variable stellar obscuration by circumstellar dust. These stars are characterized by an increase in the degree of linear polarization as the star becomes faint, supporting the idea that variable obscuration is responsible for the fading of the star. Historically the star has been observed to be redder when it is at a fainter magnitude (Evans et al. 1982; Liseau et al. 1987; Hutchinson et al. 1989; Covino et al. 1992), and the amplitude of the variations also decreases with wavelength. The intrinsic position angle of T Tauri star polarization is generally a function of both wavelength and time (Bastien 1981, 1985). These authors also noticed remarkably large and rapid variations in both polarization and position angle in RY Lup. More recently, Manset et al. (2009) has investigated these photometric and polarimetric variabilities using simultaneous BV polarimetric and UBV photometric observations. They showed that the polarization is high (3.0%) when the star is faint and red (V = 12.0 mag, B − V = 1.3 mag), and low (0.5%) when it is bright and bluer (V = 11.0 mag, B − V = 1.1 mag). The photometric and polarimetric variations share a common period of 3.75 d. They concluded that linear polarization is produced by dust scattering in an asymmetric (flat) circumstellar envelope, and they explained the photometric and polarimetric variations by an almost edge-on circumstellar disk that is warped close to the star, where it interacts with the star magnetosphere. They claim that the inhomogeneous disk matter contained in the warp could be corotating with the star and could partially occult it during part of the rotation period, which explains the dips in luminosity and the accompanying increase in polarization. However, the most recent estimate of v_{sin i} = 16.3 ± 5.3 km s⁻¹ from (Alcalá et al. 2017) is smaller than earlier measured. This implies that the stellar rotation axis has an inclination likely <70 degrees. Therefore, the disk is probably misaligned with respect to the star’s equatorial plane.

Recent ALMA high-resolution sub-mm observations of RY Lup in the 890 micron dust continuum (Ansdell et al. 2016) provide a detailed map of the spatial distribution of large, mm-sized dust, at a linear spatial resolution of 35 au. They show a clear signature of an inner mm dust cavity with a diameter of 0.8″ (60 au) and a clearly resolved dust ring. These continuum emission measurements constrain M_dust to less than 0.2 M_⊙ and additional CO isotopologueemission (₁₃CO and C₁₈O 3–2 lines) constrain M_gas to less than 2.6 M_⊙.The gas-to-dust ratio has been estimated to be 5–50 when using their parameterized model framework. Although the gas mass estimation may be underestimated due to their assumption of an ISM-like abundance, the weak CO isotopologue emission indicates rapid disk evolution, either directly in the gas-to-dust ratio or chemically via permanent loss of volatiles to solids. All previous works on RY Lup are consistent with a system comprising a star surrounded by a nearly edge-on disk. Manset et al. (2009) estimate the mass of RY Lup to M_⋆ ∕M_⊙ = 1.4, with an age of 12 Myr. This region’s proximity and age makes it ideal for a baseline study of disk properties.

Dodson-Robinson & Salyk (2011) demonstrate that transitional disks with wide cavities (>15 AU) and high accretion rate greater than 10⁻⁸ M_⊙ yr⁻¹ presents unambiguous evidence for tidal clearing by multiple planets. The cavity formation can only occur once evaporation rates exceed accretion rates, as otherwise the hole will be replenished by accretion. Photoevaporation models do not appear to be able to explain the existence of disks with such characteristics. The RY Lup disk clearly falls in that category with an estimated mass accretion rate 10^−8.2 M_⊙ yr⁻¹ from Gahm et al. (1993).

As part of the SPHERE guaranteed time observations dedicated to a large survey to search for planets and to study circumstellar disks around members of young and nearby associations, we recently recorded high-contrast images of RY Lup. The data resolve the disk for the first time in intensity in the near-infrared and in polarized light. This article presents the observational results, and aims to develop qualitative arguments on the detected disk. In parallel to the observation, detailed numerical hydrodynamical simulations in the context of planet–disk interaction in combination with radiative transfer calculations were performed to study the disk morphology. We first describe the observations and the data analysis (Sect. 2), then the results obtained on the disk and the search for planets in this system (Sect. 3). This is followed by detailed numerical modelling of the disk (Sect. 4).

2 Observations and data reduction

2.1 Observations

All observations were part of the SPHERE consortium guaranteed time programme (SPHERE High-Contrast Imaging Survey for Exoplanets and SPHERE Disk Survey). The Infra-Red Dual-beam Imager and Spectrograph (IRDIS, Dohlen et al. 2008) and Integral Field Spectrograph (IFS, Claudi et al. 2008) were used. The first RY Lup datasets were recorded on 16 April 2016 with the IRDIFS instrumental configuration, where observations are simultaneously carried out with the IRDIS subsystem dual-band imaging (DBI) mode using H2, H3 filters respectively centred at 1.59 and 1.67 μm (Vigan et al. 2010), and the Integral Field Spectrometer operating in Y J bands (0.95–1.35 μm) with a spectral resolution of 42 in Y band to 40 in J band. Following the discovery of the disk in intensity, we obtained Dual-band Polarimetric Imaging (DPI, Langlois et al. 2010) observations in field stabilized mode on 27 May 2016. IRDIS provides a 10–12′′ square field of view (1500–1800 square au, given the star’s distance 151 pc, using 12.25 mas/pixel platescale i.e. 1.85 au per pixel). The IFS dataset consists of 21 000 spectra each spread over 5.1 ×41 pixels on the detector. After extraction, the FoV is 1.7′′ square and the spaxel size is 7.46 × 7.46 mas, i.e. 1.13 au per pixel. We used an apodized Lyot coronagraph (Soummer 2005; Boccaletti et al. 2008; Carbillet et al. 2011; Guerri et al. 2011), including a 185 mas focal mask, as well as a pupil mask (N-ALC-YJH-S). The IRDIFS observations were performed in pupil stabilized mode in order to perform Angular Differential Imaging (ADI) post-processing (Marois et al. 2006). The field rotation for this set of data is 73 degrees. The observing sequences for both IRDIS and IFS can be summarized as follows: (1) point spread function (PSF) imaging, including a small star offset to displace the star away from the coronagraphic mask in order to record the unsaturated PSF which provides relative photometric calibration (we note that a neutral density filter was also insertedto avoid saturation); (2) imaging of the star behind the coronagraphic mask, with four crosswise faint replicas of the star artificially generated using the deformable mirror (Langlois et al. 2013) used for fine monitoring of the star centre; and (3) science coronographic sequence. At the end of the observing sequence, we repeated the first and second steps followed by a calibration of the sky background, with detector integration times (DITs) corresponding to the DITs of the coronagraphic observations. Finally, the true north (TN) and pixel plate scales were measured using the astrometric calibrator NGC3603 (Khorrami et al. 2016) observed as part of the standard SPHERE calibration routines performed for the Guaranteed Time Observations (GTO) survey (Maire et al. 2016). More details on the datasets are presented in Table 1. The observing conditions were very good for all observing sequences. The Strehl ratio estimation provided in Table 1 is based on an extrapolation of the phase variance deduced from the reconstruction of the SPHERE Adaptive Optics system open-loop data using deformable mirror and tip-tilt voltages, and wavefront sensor closed-loop data (Fusco et al. 2004). The observing conditions and the different data reduction methods for each dataset are described in more detail in Sects. 2.2 and 2.3.

2.2 Intensity images using IRDIS-DBI (H2, H3 bands) and IFS (YJ band)

The IRDIS data were first corrected for cosmetics and sky background using the SPHERE Data Reduction and Handling (DRH) pipeline (Pavlov et al. 2008) implemented at the SPHERE Data Centre. The outputs include cubes of left and right images recentred onto a common origin (the centre of the field rotation) using the satellite spots. This preprocessing also includes background subtraction, bad-pixel interpolation, flat-field correction, distortion correction, and wavelength calibration (for IFS). After these first steps, the best frames were selected according to their quality, leading to the use of 72 frames out of 80 for IRDIS and of all the available frames for IFS. The IRDIS frame selection criteria is based on the central coronagraphic spot flux being within a ± 2σ deviation from the median value of the full sequence.

Following these steps the data were processed to remove the stellar halo (i.e. in order to achieve high-contrast) with the SpeCal pipeline (Galicher et al., in prep.) developed to reduce and analyse our GTO data from the SPHERE High-contrast Imaging Survey for Exoplanets (SHINE) survey. This pipeline implements a variety of ADI-based algorithms from which we used classical Angular Differential Imaging (cADI, Marois et al. 2006), Template Locally Optimized Combination of Images (TLOCI, Marois et al. 2014), and reference differential imaging (RDI) using a compatible reference star from our SPHERE GTO target library. In parallel, the IFS data were also processed using angular and spectral principal component analysis (PCA-ASDI) as describedin Mesa et al. (2015). In the following we discuss the results based on the noADI (classical averaging with no subtraction), cADI (see Fig. 1), and RDI for the disk morphology. On the other hand, the photometric analysis of the point sources is based on TLOCI since it delivers more accurate photometry (Galicher et al., in prep.)

2.3 IRDIS-DPI (BBH band)

The IRDIS-DPI observations of RY Lup were carried out on 27 May 2016 with the BBH filter (1.625 μm, width = 290 nm) using the same apodized pupil Lyot coronagraph (N-ALC-YJH-S) as for the IRDIFS observations. Sixty-four polarimetriccycles were taken, each consisting of one data cube for each of the four half wave plate (HWP) positions (0, 45, 22.5, 67.5 degrees). Dedicated coronagraphic images were taken at the beginning and at the end of the science sequence to determine accurately the star centre behind the coronagraph using the four crosswise faint replicas created by the deformable mirror. The IRDIS data were first corrected for cosmetics using the DRH pipeline implemented at the SPHERE Data Centre and following the same steps as the DBI data. The outputs include cubes of left and right images (parallel and perpendicular polarized beams, respectively) recentred onto a common origin. The data were then reduced following the prescriptions of Avenhaus et al. (2014), using the azimuthal Stokes formalism (Q_φ, U_φ) and relative polarimetric measurements. The parallel and perpendicular polarized images are combined to produce Q and U Stokes images. To obtain cleanStokes Q and U images, that is to correct for instrumental polarization downstream of the HWPs position in the optical path, Q+ and Q− (0 and 45 degrees), and U+ and U− (22 and 67.5 degrees) are included in the image combination and used to compute the azimuthal Stokes components. The azimuth is defined with respect to the star position as shown by Canovas et al. (2015) and the U_φ signal is very small for centrally symmetric disks, but not for very inclined disks as shown by Pohl et al. (2017).

However, there might still be an instrumental polarization left upstream of the HWP in the final Q and U images, which is assumed to be proportional to the total intensity image as shown in Canovas et al. (2011). Since the RY Lup disk is inclined, the first standard reduction did not include U_φ minimization of the signal in an annulus around the central star. We also reduced the IRDIS-DPI data with a new method based on an inverse approach method (Denneulin et al., in prep.). This method also enables the use of any set of frames with incomplete polarimetric cycles. The method is based on an optimization using the electric field (Jones Matrix) rather than the intensity (Mueller matrix) and relies on the inverse approach (i.e. fitting a model of the data to the observed dataset) which is significantly less biased and more efficient at minimizing instrumental artefacts. No assumptions about the angle of linear polarization of the source are made to correct for the instrumental polarization. The method is based on the use of an instrumental model describing the effect of the instrument on the incident electromagnetic field. The electric field is divided into non-polarized and polarized components, $$I_k=\Vert v_k{\Vert }^2\rho +|\langle v_k,m\rangle {|}^2+\text{noise},$$(1)

where I_k is the flux measured on the detector,ρ is the non-polarized intensity (unknown), m is the electric field linearly polarized (unknown), and v describes the effect of the instrument/telescope for a given half wave plate polarizer position. The model describes every polarimetric measurement and the polarized intensity is obtained by solving for both polarized and unpolarized components of the electrical field for every pixel individually using the criteria $$\underset{m,\rho }{min}\sum _k\frac{{\omega }_k}{2}{\Vert \Vert v_k{\Vert }^2\rho +|\langle v_k,m\rangle {|}^2-I_k\Vert }^2,$$(2)

where ω_k is a weighting factor calculated from the inverse of the noise variance. The derived Q_φ and U_φ estimations using this latest method is similar to the standard reduction we have performed, but the estimation of the polarization angles is more accurate. In the RY Lup case the disk is very bright which makes the level of instrumental polarization left in the Q_φ and U_φ images very small compared to the disk signal, and we conclude that the disk morphology is not affected by these polarimetric residuals. The Q_φ and U_φ images are displayed in Figs. 2 and 3. Figure 4 shows the polarized intensity overplotted with polarization vectors representing the angle of linear polarization. This strengthens the hypothesis that there is a clear departure from azimuthal polarization in inclined disks such as RY Lup.

Fig. 1 IRDIS cADI (top left) and noADI (top right) H2+H3 image normalized to the maximum of the unsaturated non-coronagraphic PSF showing the close-in environment of RY Lup. cADI (bottom left) and noADI (bottom right) IFS YJ image normalized to the maximum of the unsaturated non-coronagraphic PSF showing the close-in environment of RY Lup. The dark central region corresponds to the area masked by the coronagraph. The intensity scale is logarithmic. The orientation is standard: north is up and the east is towards the left

Fig. 2 IRDIS DPI Qφ (left) and Uφ (right) images showing the close environment of RY Lup. The dark central region corresponds to the area masked by the coronagraph. The intensity scale follows the inverse hyperbolic sine.

Fig. 3 Top: polarized intensity image showing the close-in environment of RY Lup. Bottom: unsharp masking of the polarized intensity image showing the close-in environment of RY Lup. The intensity scale follows the inverse hyperbolic sine.

3 Results

We obtained a very clear detection of the RY Lup disk in all datasets presented in this paper. The analysis of the morphology of the detecteddisk primarily focuses on the IRDIS and IFS intensity images and IRDIS-DPI (see Sects. 2.2 and 2.3). Furthermore, the cADI IFS and IRDIS images are used to search for point-sources focusing on non-polarized companions because the intensity images reach higher contrast (Sect. 3.2). Our physical modelling of the polarimetricand total intensity images allows a qualitative analysis of the disk, and is presented in Sect. 4.

3.1 Disk properties

The IRDIS and IFS datasets reveal clearly the inclined disk around RY Lup as seen in Fig. 1 in intensity and in Fig. 2 in polarimetry. Figure 1 for IFS is a median image over wavelength range from Y to J band. In our images, the disk appears as a dominating double-arch structure in the SE–NW direction extending to 0.53′′ (80 au) in radius. Our observations support a high disk inclination with respect to the line of sight, with a position angle of 107 ±1 degrees. The position angle was estimated by considering the brightest parts of the disk. In addition to these bright components the images reveal the presence of two fainter spiral-like features on the SE direction as seen on the display of the Q_φ image (Fig. 2). A spiral-like feature is also visible on the opposite side of the disk along its major axis. This either means that the m = 1 and m =2 modes are basically seen on opposite sides of the disk or, together with the inner SE spiral arm, these two components could be approximately traced by an ellipse. In Fig. 3 we also applied an unsharp mask to the polarized intensity in order to enhance the disk structure. Other double-arch structures in TDs have been recently reported with high-contrast imaging instruments such as SPHERE (Janson et al. 2016; Garufi et al. 2016; Pohl et al. 2017), but spiral-like features such as the one observed here have not yet been detected in very inclined disks. The spiral-like features detected in the RY Lup disk images could be created by an interaction of the disk with a planet located in its surrounding (Dong et al. 2016). We demonstrate in Sect. 5 that it is a reliable hypothesis.

We note that the ADI processing of these images may have been biased and are not a faithful representation of the true intensity and geometry. In this particular case the disk is bright and these artefacts are negligible, as shown by the similarities between the cADI, noADI, and the DPI images.

The reduced U_φ image signal can be interpreted as radial polarization. The high peak-to-peak ratio between U_φ and Q_φ (30%) can be attributed to the high inclination of the disk (~70 degrees as determined from the total intensity image). In a highly inclined disk, multiple scattering can occur (i.e. where light that has already been polarized is scattered). In this case, this is the prime contributor to the U_φ signal. This is consistent with a theoretical study by Canovas et al. (2015) and Pohl et al. (2017), who found that multiple scattering can produce significant non-azimuthal polarization. They showed that the peak-to-peak ratio between U_φ and Q_φ can even be as high as 50% for a disk inclination of 70 degrees depending on the mass and grain size distribution of the disk model.

We used the RDI IFS, RDI IRDIS, and DPI IRDIS images to derive the relative surface brightness (SB − M_⋆) profiles of the disk along its major axis (Fig. 5), assuming that the major axis has a position angle of 107 degrees. These profiles are calculated using an average width of 5 pixels along the midplane with the following formula: SB − M_⋆ = 2.5 × log(normalized Intensity/pixel−area²) − M_⋆. The conversion of the intensity into mag arcsec⁻² was performed using the 2MASS stellar magnitudes (J = 8.54 mag, H = 7.69 mag, Cutri et al. 2003) and the ratio of the maximum to the total flux of the measured unsaturated non-coronagraphic PSF. With this normalization choice, the profiles provide information on the scattering efficiency of the dust grains. We note some brightness asymmetry between the west and east disk wings which is slightly more pronounced in the polarimetric intensity (represented with arbitrary surface brightness scale) than on the intensity cADI profiles (IFS and IRDIS). We also note that the radial profiles are very similar from the Y to H band, which suggests a shallow dependence efficiency of the scattering with wavelength. This is confirmed by the RY Lup disk spectrum ranging from 5.05 to 5.25 mag arcsec⁻² − M_⋆ across the 0.9–2.3 μm wavelength range. This spectrum was estimated by computing the integrated flux from the IFS image at each wavelength in a constant area, defined by the area where the flux is above 1.5σ of the collapsed image over all wavelengths. The observed brightness and colour of scattered light images put strong constraints on the scattering properties of protoplanetary dust as shown by Mulders & Dominik (2012). The very flat spectrum of the RY Lup disk suggests a minimum grain size of ~ 3 μm assuming astrosilicates and a grain size distribution slope of − 3.5.

Fig. 4 Polarized intensity image showing the close-in environment of RY Lup. The white stripes represent the angle of linear polarization (fixed length, not scaled with the degree of polarization). The intensity scale follows the inverse hyperbolic sine.

3.2 Search for planets

Twenty candidate companions (CC) were detected in the IRDIS field of view (Fig. 7) using DBI with H2H3 filters, whereas no point-like source was found in the IFS image. The photometry and astrometry of the companions were measured using the TLOCI algorithm applied to each spectral band separately as presented in Table 2. We divided each science frame into annuli of 1.5 full width at half maximum (FWHM, where FWHM= 41 mas). For each science frame and annulus, we computed the stellar residuals using the best linear combination of the most correlated frames for which the self-subtraction of mock point sources was at maximum 20%. The astrometry and the photometry were estimated by a fit of the point spread function template onto the data. The photometric errors were estimated conservatively by including the variations in the stellar flux during the sequence (estimated from the fluctuations of the stellar residuals), the accuracy of the fitting procedure, and the PSF variability between the beginning and the end of the sequence. The astrometry of the detected companions was calibrated using pixel scales of 12.255 mas and a TN angle offset of −1.742 degrees (Maire et al. 2016). The astrometric errors include the accuracy of the fitting procedure and the star centring error.

Because of the photometric bias from the TLOCI algorithm, we estimated the throughput of the technique at each position in the field using fake planet injection and produced throughput-corrected final images. The contrast curve for each spectral channel is estimated by the azimuthal standard deviation of the throughput corrected image in an annulus of 0.5 FWHM width. No small statistical correction at short separation was applied (Mawet et al. 2014). The azimuthally averaged contrast curves shown in Fig. 6 were estimated for the TLOCI reduction for IRDIS and PCA ASDI reduction for IFS because these reductions provide the best compromise for contrast, stellar rejection, and self-subtraction correction for the point source detection. The disk is bright, but strongly attenuated by the TLOCI algorithm (which is optimized for point-sourcedetection). As a consequence, the disk residuals increase the measured noise level by a small amount leading to slightly pessimistic contrast curves at very short distances. For the conversion of the contrast limits to mass limits in Fig. 6, we used the COND atmospheric and evolutionary models of Baraffe et al. (2015) and Baraffe et al. (2003). All the point sources were classified as probable background stellar objects using the SPHERE consortium tools developed for the classification and the ranking of the companion candidates discovered in SPHERE/SHINE. With the same suite of tools, we derived the colour magnitude diagram (CMD) of all point sources, as shown in Fig. 8, to compare their H2–H3 colour to those of field and young dwarfs covering theM, L, T, and Y types. Further details about the derivation of this colour magnitude diagram are provided in Zurlo et al. (2016). We also computed a background probability for the point sources by assuming the Besancon model predictions (Robin et al. 2003). As seen in Table 2, there is only one candidate with low probability of contamination, but its position on the CMD diagram corresponds to a background star. All the detected point sources also have very similar offsets from the M–L sequence and their separations are also greater than a hundred au, so we can conclude that they have a very low likelihood of being bound companions.

Fig. 5 RY Lup disk surface brightness profile measured along the midplane on the IFS intensity image (asterisk line) expressed in mag arcsec−2 − M⋆, IRDIS intensity image (red dotted line) expressed in mag arcsec−2 − M⋆, and IRDIS polarimetry image (blue dashed line) expressed in arbitrary units. The grey area represents the separation range masked by the coronagraph. The intensity is normalized to the unsaturated non-coronagraphic PSF intensity peak. The error bars are within 0.1 mag.arcsec−2.

4 Physical disk modelling

Our physical disk modelling aims to reproduce the polarized intensity, the orientation of the polarization vectors, and the two faint spiral arms observed on the SE and SW sides in the SPHERE intensity and polarimetric images. This requires the loss of disk symmetry, which could be explained by planet–disk perturbation dynamics. Thus, we want to test whether the spiral pattern triggered by a Jupiter-like planet could eventually result in the scattered light features that we detect around RY Lup, assuming that its disk is seen close to edge-on. For this purpose, we carry out 3D global numerical hydrodynamical simulations to study the morphology of spiral arms excited by an outer planetary perturber. Using detailed 3D follow-up radiative transfer models, we produce synthetic scattered light images for a highly inclined disk configuration at NIR wavelengths. This modelling approach is motivated by recent work on spiral arms in the context of scattered light by Dong et al. (2015, 2016), Juhász et al. (2015), Pohl et al. (2015).

4.1 Hydrodynamical planet–disk simulations

We use the Three-dimensional RAdiation-hydrodynamical Modelling Project (TRAMP, Klahr et al. 1999) code in its extended version (Klahr & Kley 2006), which applies a flux-limited diffusion approximation for the radiation. We consider a giant planet embedded in a fully viscous disk, where the consistent release of accretion energy fromthe planet via radiation is included. The disk model considers a spherical polar coordinate system (r, θ, φ), where the disk midplane coincides with the θ = 0 plane. The number of grid cells in the radial (r), polar (θ), and azimuthal (φ) directions are 64, 31, and 89, respectively. The computational disk grid has a radial extent from ~18 au to 235 au. The complete set of vertical disk structure equations is solved self-consistently, including both viscous dissipation and heating by irradiation from the central star (cf. D’Alessio et al. 1998). This results in detailed profiles of the density and temperature structure with vertical height and disk radius. Thus, it allows us to determine the 3D temperature structure of the disk at each step of the planet–disk interaction processes. More precisely, the initial 3D dataset, before adding the planet, is determined from a set of 1D vertical structure models for a given accretion rate (which will result in the desired total disk mass). Thus, the dynamical viscous evolution in the 3D radiation hydro simulation of the disk will not significantly alter the initial global gas distribution, but any development of structure in the 3D simulations is the result of the torques exerted by the planet. We make sure that the 1D and 3D models use the same parameters including the dust opacities and alpha viscosity. We assume a central stellar mass of M_⋆ = 1.4 M_⊙, a radius of R_⋆ = 1.72 R_⊙, and an effective temperature of T_eff = 5669 K. These parameters are consistent with the properties of RY Lup derived by a piecewise linear interpolation of the values from Manset et al. (2009) assuming an object’s distance of 151 pc (Gaia Collaboration 2016). Furthermore, the input parameters for the viscosity parameter and the mass accretion rate are chosen to be α = 10⁻³ and M_accr,0 = 3.3 × 10⁻¹⁰ M_⊙ yr⁻¹, respectively. The latter is rather low compared to typical values for T Tauri disks, but a natural outcome of the disk mass considered. The accretion rate is calculated such that the disk gas mass corresponds to M_disk = 2.5 × 10⁻³ M_⊙, distributed between an inner disk radius of r_in = 18 au and an outer disk radius of r_out = 123 au, as measured by Ansdell et al. (2016).

The position of a possible planetary perturber is uncertain, but spiral arms exterior to a planet’s orbit are unlikely to explain the observations as they are too tightly wound given typical disk scale height values (Juhász et al. 2015). However, a location further out in the disk is compatible with the outer disk radius in the ¹³ CO ALMA emission map. Furthermore, the mm dust continuum suggests a truncation at ~120 au and a depletion of large dust beyond (Ansdell et al. 2016). Following the approach by recent spiral arm models (Dong et al. 2015), the planet can be roughly located at a distance three times the location of the inner arms. Thus, we consider a planet located at 190 au, that is outside of the disk radius in both scattered light and in thermal emission. Given that the planet carves out a wide gap in the disk, this location is a justified assumption. We assume a fixed circular orbit, although spiral density waves may be also excited by a companion on an eccentric orbit. The mass of the embedded giant planet is considered to be M_pl ∕M_⋆ = 1.4 × 10⁻³ (~2 M_jup). This mass is restricted by the mass detection limits for planetary companions obtained from our SPHERE contrast curves. Figure 6 shows that the upper mass limit at an angular distance of ~ 1.1′′ is ~ 2 M_jup. It is worth noting that the mass limit is model dependent and slightly higher if a warm start model is assumed (~5 M_jup) (Marley et al. 2007; Spiegel & Burrows 2012). It is also worth mentioning that a planet might be at a location behind the disk or might be geometrically close to the central star, where it could remain undetected. Hence, it is geometrically likely to have a projected separation less than ~ 1.1′′, at which a higher planet mass would be possible according to Fig. 6. However, our choice of a lower planet mass is thus conservative, which is reasonable given the uncertainties on the age of the system and on the atmospheric/evolutionary models used to estimate themass detection limits. The planet–disk features, especially the contrast between spirals and the background disk, would be even stronger for higher planet masses (cf. Dong et al. 2015; Pohl et al. 2015). The simulation is run for 90 orbits at 190 au (~ 1.8 × 10⁵ yr), which is sufficiently long for the inner spiral arms to reach a quasi steady state, although this evolutionary timescale would be not enough to carve out a full gap around its orbit. This is, however, reasonable for our modelling purposes since the planet is located outside the disk and we are mainly interested in the observational appearance of the spiral arms.

A map of the surface density Σ_gas after an evolutionary time of 90 orbits at 190 au is shown in Fig. 9. While the planet always excites inner and outer waves, we concentrate here on the m = 2 spiral arm configuration inward of the planet’s position at 190 au. The primary arm originates from the planet location, the secondary is shifted by ~180° in the azimuthal direction. The density contrast of the secondary spiral wave with respect to the background disk is weaker than the primary one.

Fig. 6 Contrast curves and companion mass limits derived for IFS (black curve) using the ASDI PCA reduction method, and for IRDIS H2 and H3 bands (blue and red curves, respectively) using the TLOCI reduction method. The coronagraph inner working angle (radius where the coronagraphic transmission is 50%) is 0.095 arcsec in radius.

Fig. 7 IRDIS TLOCI H2+H3 image normalized to the maximum of the unsaturated non-coronagraphic PSF showing the detected point sources in the environment of RY Lup. The intensity scale is linear.

Fig. 8 Colour magnitude diagram displaying the discovered candidate companions, which are marked in red including photometric error bars, compared to known substellar field (coloured symbols) and young objects. This plot assumes that the candidate companions are at the same distance as the star. The candidates were classified as probable background objects based on their position in the CMD. Several candidates lying near the L–T transition could still be considered as potential companions.

Fig. 9 Simulated 2D vertically integrated density map for a planet-induced spiral wave considering a planet-to-star mass ratio of 1.4 × 10−3. The green dot highlights the position of the planet.

4.2 Radiative transfer models

The resulting 3D disk density and temperature structure presented in Sect. 4.1 is subsequently fed into the Monte Carlo (MC) radiative transfer code RADMC-3D developed by Dullemond et al. (2012)¹. To convert the gas density to the dust density used in the radiative transfer calculations, we adopt a gas-to-dust ratio of 50, which is lower than the canonical ISM value as suggested by the ALMA Lupus survey of Ansdell et al. (2016). For the dust opacity calculation we consider a species mixture of silicate (60%), carbonaceous (10%), and icy (30%) material, where the optical constants are taken from Draine (2003), Zubko et al. (1996), Warren & Brandt (2008). We assume two different grain size populations, small (0.01–1 μm) and large (1 μm–1 mm) grains, where each size distribution is a smooth power law with exponent − 3.5. The mass ratio between the two populations is determined such that the number density follows n(a) ∝ a^−3.5. Thus, the large fraction of mass is in the large grains (~95% of total dustmass), but the small grains dominate the opacity at NIR wavelengths. We assume that the gas and the small dust is well mixed, thus their dust scale height is equal to the pressure scale height. We are aware that this choice does not reproduce the scattering colours of the disk, but we only focus on the morphology of the image in our model approach.

The stellar parameters for the radiation source as well as the radial disk extension are taken to be the same as in the hydrodynamical calculations. For the polar and azimuth we use a finer grid sampling with N_φ = 192 and N_θ = 64, and interpolate the density and temperature array values accordingly. All scattering MC simulations are run with 5 × 10⁸ photon packages. We consider multiple scattering effects and produce images for all four Stokes components, from which the total intensity and polarized intensity images are calculated at H band. The disk orientation is set by the inclination angle (0° is face-on) and the position angle (PA, from north to east). The theoretical radiative transfer images are convolved by the instrument PSF with a FWHM of 0. ′′ 04 so that it mimicks the angular resolution of the data.

5 RT Modelling results

5.1 Synthetic scattered light images

The convolved H-band polarized intensity image for our planet–disk interaction model in the RY Lup system is shown in Fig. 10. The disk is inclined by 75 degrees with respect to the line of sight and rotated by a PA of 107 deg relative to the disk major axis. We note that we ran a grid of models for inclinations between 70 and 80 degrees in steps of 1 deg, where the data is best represented by an inclination of 75 degrees. This is consistent with the measurement from total intensity given the uncertainties. For such a highly inclined system, the symmetric two-spiral arms produced by the external planetary perturber and known from face-on disk configurations can have a completely different morphology in scattered light (cf. Dong et al. 2016). The spiral arms (m = 1 and m = 2 modes) instead appear as two arms located along the major axis at the top (illuminated) side of the disk (see sketch in Fig. 11). There is also a substructure branching off the SE spiral, possibly a second winding of the spiral density wave, but in the convolved synthetic image it is not as pronounced as in the observed PDI image. Both spiral arms have an almost equal brightness contrast in the scattered light image. This result is consistent with a parametric study by Dong et al. (2016) on how spirals driven by companions can appear in scattered light at arbitrary viewing angles.

5.2 Polarimetric comparison

Our disk model leads to a qualitatively good match with the IRDIS polarimetric image from Fig. 2. The bright double-arch structure in the SE–NW direction extending to 0.53′′ in radius to the NW is well reproduced and corresponds to scattered light from the bottom side of the disk. The model also confirms that the disk must have a high inclination with respect to the line of sight. In addition to these bright components the two fainter spiral-like features in the SE direction, which break off, are also reproduced. Due to the multiple disk and planet parameters to adjust in these models and due to the high inclination of the disk, there is no perfect match of the brightness ratios between the east and west sides, and the bright and fainter features of the disk. For such a highly inclined system, the symmetric two spiral arms generated by the external planetary perturber in fact have very different morphology in scattered light with the viewing angle.

Figure 10 shows our polarized intensity model for RY Lup overlaid with linear polarization vectors. We note that all vectors have the same length as we are especially interested in the polarization orientation. As also detected in the IRDIS polarimetric image (Fig. 3), it is noticeable that the polarization direction along the minor axis is radial, which means that there is a significant non-azimuthal polarization. Since the polarization angle is indeed more complicated than expected for single scattering, the U_φ signal must be physical and connected to multiple scattering, which is included in our radiative transfer calculations as well.

Fig. 10 Convolved polarized intensity H-band image considering a planet-induced spiral model with Mpl∕M⋆ = 1.4 × 10−3. The image is produced at an inclination of 75 deg, and a PAs of 107 deg. The inner 0.18″ are masked (black circular area) to mimic the coronagraph region. The white stripes represent the angle of linear polarization (fixed length, not scaled with the degree of polarization).

Fig. 11 Sketch of the disk model suggested by our observations and modelling. North is up and east is towards the left. The major axis of the disk is indicated by the dashed line and divides the disk into a near and far side. The top (illuminated) and bottom (obscured) halves are separated by the disk midplane.

5.3 Spectral energy distribution

The SED model based on our hydrodynamical simulation compared to the observed photometry is shown in Fig. 12. The stellarspectrum is taken from the Castelli–Kurucz atlas (Castelli & Kurucz 2004), choosing the model for a stellar type of G8V (T_eff = 5500 K), and a surface gravity of log g = 4.5. The flux measurements from the VLT Interferometer (VLTI)/PIONIER, SPHERE/IRDIS/IFS, 2MASS, IRAS, WISE, and AKARI hint at a strong NIR and mid-IR (MIR) excess expected for a disk around a young star. Although our model lies within the measurement uncertainties of the SPHERE data points, there is a clear trend that it underestimates the fluxes between ~5 and ~50 μm. The explanation is that our model only considers the dust disk from our hydrodynamical simulations with a radial grid starting at R_in = 18 au. The inner disk region is hidden behind the coronagraph for our SPHERE data, so there is no strong constraint on the geometry and orientation of an inner dust belt or halo. This makes the modelling of the inner disk regions very degenerate, and thus, in our model there is no dusty material close to the star. In general, the inner radius regulates the maximum temperature of the dust grains at the inner rim. For larger R_in the luminosity excess shifts from the NIR to the MIR (see e.g. Woitke et al. 2016). Exploring a larger grid of dust and disk structure parameters is beyond the scope of the paper as we did not attempt to do a detailed fit of the SED. The photometric point in the mm regime is nicely reproduced with our simulated flux. Due to the absence of additional photometric points in the (sub-)mm, the mm slope fit cannot be evaluated. For example, the dust size distribution power-law index (typical value of 3.5 in our case for both populations) changes, in particular, the mm and cm slopes. Although our current model indicates a good disk mass estimate, this situation might change when additional dust is included in the inner disk regions in order to fit the NIR/MIR fluxes. The underestimate of this emission could result in a slight error in dust temperature and requires a higher disk mass. This is certainly compatible with the disk mass estimate from Ansdell et al. (2016) as this value only gives a lower limit.

Fig. 12 Synthetic SED from our radiative transfer model compared with direct measurements obtained with SPHERE IRDIS and IFS in violet, VLTI/PIONIER (Anthonioz et al. 2015), 2MASS, IRAS, WISE, and AKARI2 and ALMA (Ansdell et al. 2016) in grey.

6 Conclusions

We observed the large circumstellar disk around RY Lup with SPHERE/IRDIS in scattered light with high angular resolution using polarimetric and angular differential imaging and uncovered directly for the first time an inclined disk with spirals. The disk position angle is 107 degrees and its inclination is 70 degrees. Our observations show the complementarity between polarimetric and angular differential imaging. Polarimetric observations allow the study of the disk without the self-subtraction effects inherent to ADI-processed data, while the ADI observations provide deeper detection limits further away from the star, enabling us to study the disk and search for exoplanets. We also retrieved a high signal-to-noise ratio spectrum for the disk using the IFS. The disk colour is grey, which is an indication that dust grains larger than the wavelength dominate the scattering opacity in the disk surface.

We have studied the morphology and surface brightness of the disk and formulated a hypothesis on the origin of the spiral arms. An explanation for the spiral arms could not be uniquely determined due to the high inclination of the disk. Our numerical model shows that the spiral arms could be explained by one low-mass planet (2 M_jup) located at 190 au orbiting exterior to the spiral arms. Although this model describes the formation of spiral arms qualitatively similar to the features in the NIR scattered light observations of RY Lup, it is only a suggested configuration for the system, and not a best-fitting model. However, given the high inclination of the disk it is unlikely that the planet’s thermal radiation is directly detectable. While spirals can be excited by the tidal interaction with the companion, they can also be triggered by the close proximity of shadows in the disk as discussed in Montesinos et al. (2016), Benisty et al. (2017). Also, we cannot exclude that the 0.8′′ diameter gap detected in sub-mm continuum emission could host additional massive planets that could play a role in the disk morphology.The innermost gap discovered in the ALMA observation cannot be detected in scattered light due to the high inclination of the disk. The scattered light flux shows a small profile asymmetry which does not coincide with the symmetry of the sub-mm continuum emission. This is likely the result of surface density perturbation related to the presence of the spiral arms and could be strongly related to the viewing orientation.

Our observations are compatible with the hypothesis made by Manset et al. (2009) of photometric and polarimetric variations created by an almost edge-on circumstellar disk that is warped (or inclined) close to the star where it interacts with the star magnetosphere. This configuration could lead to a partial shadowing of the outer disk and a brightness asymmetry that can remain undetected due to the high inclination of the disk. We also performed a detailed radiative transfer study, which reproduce well our scattered light observations. The planet–disk interaction scenario is in agreement with the ALMA high-resolution dust continuum image being truncated at ~120 au. This modelling effort cannot fully constrain all the parameters in this case, but helps to understand the disk structure and distribution of small grains in particular.

Acknowledgements

We acknowledge our anonymous reviewers for the very careful reading and constructive suggestions that contributed to improving this paper. 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 astronomique de l’universite de Geneve (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France), and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland), and NOVA (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). This work was supported by the Programme National de Planétologie (PNP) and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU co-funded by CNES. This work has also been supported by a grant from the French Labex OSUG2020 (Investissements d’avenir - ANR10 LABX56) and by the ANR grant ANR-14-CE33-0018 (GIPSE). This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS /LAM/CeSAM (Marseille), OCA/Lagrange (Nice), and Observatoire de Paris/LESIA (Paris). This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France. Q. Kral also acknowledges support from STFC (consolidated grant) via the institute of Asronomy, Cambridge.

References

All Tables

All Figures

Fig. 1 IRDIS cADI (top left) and noADI (top right) H2+H3 image normalized to the maximum of the unsaturated non-coronagraphic PSF showing the close-in environment of RY Lup. cADI (bottom left) and noADI (bottom right) IFS YJ image normalized to the maximum of the unsaturated non-coronagraphic PSF showing the close-in environment of RY Lup. The dark central region corresponds to the area masked by the coronagraph. The intensity scale is logarithmic. The orientation is standard: north is up and the east is towards the left

In the text

	Fig. 2 IRDIS DPI Qφ (left) and Uφ (right) images showing the close environment of RY Lup. The dark central region corresponds to the area masked by the coronagraph. The intensity scale follows the inverse hyperbolic sine.
In the text

	Fig. 3 Top: polarized intensity image showing the close-in environment of RY Lup. Bottom: unsharp masking of the polarized intensity image showing the close-in environment of RY Lup. The intensity scale follows the inverse hyperbolic sine.
In the text

	Fig. 4 Polarized intensity image showing the close-in environment of RY Lup. The white stripes represent the angle of linear polarization (fixed length, not scaled with the degree of polarization). The intensity scale follows the inverse hyperbolic sine.
In the text

Fig. 5 RY Lup disk surface brightness profile measured along the midplane on the IFS intensity image (asterisk line) expressed in mag arcsec−2 − M⋆, IRDIS intensity image (red dotted line) expressed in mag arcsec−2 − M⋆, and IRDIS polarimetry image (blue dashed line) expressed in arbitrary units. The grey area represents the separation range masked by the coronagraph. The intensity is normalized to the unsaturated non-coronagraphic PSF intensity peak. The error bars are within 0.1 mag.arcsec−2.

In the text

	Fig. 6 Contrast curves and companion mass limits derived for IFS (black curve) using the ASDI PCA reduction method, and for IRDIS H2 and H3 bands (blue and red curves, respectively) using the TLOCI reduction method. The coronagraph inner working angle (radius where the coronagraphic transmission is 50%) is 0.095 arcsec in radius.
In the text

	Fig. 7 IRDIS TLOCI H2+H3 image normalized to the maximum of the unsaturated non-coronagraphic PSF showing the detected point sources in the environment of RY Lup. The intensity scale is linear.
In the text

Fig. 8 Colour magnitude diagram displaying the discovered candidate companions, which are marked in red including photometric error bars, compared to known substellar field (coloured symbols) and young objects. This plot assumes that the candidate companions are at the same distance as the star. The candidates were classified as probable background objects based on their position in the CMD. Several candidates lying near the L–T transition could still be considered as potential companions.

In the text

	Fig. 9 Simulated 2D vertically integrated density map for a planet-induced spiral wave considering a planet-to-star mass ratio of 1.4 × 10−3. The green dot highlights the position of the planet.
In the text

	Fig. 10 Convolved polarized intensity H-band image considering a planet-induced spiral model with Mpl∕M⋆ = 1.4 × 10−3. The image is produced at an inclination of 75 deg, and a PAs of 107 deg. The inner 0.18″ are masked (black circular area) to mimic the coronagraph region. The white stripes represent the angle of linear polarization (fixed length, not scaled with the degree of polarization).
In the text

	Fig. 11 Sketch of the disk model suggested by our observations and modelling. North is up and east is towards the left. The major axis of the disk is indicated by the dashed line and divides the disk into a near and far side. The top (illuminated) and bottom (obscured) halves are separated by the disk midplane.
In the text

	Fig. 12 Synthetic SED from our radiative transfer model compared with direct measurements obtained with SPHERE IRDIS and IFS in violet, VLTI/PIONIER (Anthonioz et al. 2015), 2MASS, IRAS, WISE, and AKARI2 and ALMA (Ansdell et al. 2016) in grey.
In the text