Investigating point sources in MWC 758 with SPHERE

Context. Spiral arms in protoplanetary disks could be shown to be the manifestation of density waves launched by protoplanets and propagating in the gaseous component of the disk. At least two point sources have been identified in the L band in the MWC 758 system as planetary mass object candidates. Aims. We used VLT/SPHERE to search for counterparts of these candidates in the H and K bands, and to characterize the morphology of the spiral arms . Methods. The data were processed with now-standard techniques in high-contrast imaging to determine the limits of detection, and to compare them to the luminosity derived from L band observations. Results. In considering the evolutionary, atmospheric, and opacity models we were not able to confirm the two former detections of point sources performed in the L band. In addition, the analysis of the spiral arms from a dynamical point of view does not support the hypothesis that these candidates comprise the origin of the spirals. Conclusions. Deeper observations and longer timescales will be required to identify the actual source of the spiral arms in MWC 758.


Introduction
While studies of the demographics of exoplanets are well underway, the processes leading to planet formation are still poorly constrained. Understanding how planets form is a vast topic which requires theoretical works, but the most compelling facts could be provided by the direct observation of very young systems, of a few Myrs old in age, which is precisely the moment when planets are still in the formation stage. Regardless of the difficulty in determining the ages of young systems, they are still very much embedded in their envelope. Opacity effects are critical and scale inversely with the wavelength, given that long wavelengths in the thermal regime are more appropriate for reaching the midplane of a disk, whereas short wavelengths in the scattered light regime trace the surface of flared protoplanetary disks.
In such conditions, finding evidence of exoplanets could rely on dynamically induced structures, with observable signatures in the disk morphology. One obvious case is the presence of spiral arms, which -in the case of disks that are not especially massive -could be attributed to density waves launched by planets (Goldreich & Tremaine 1980), but also potentially by other types of perturbations (e.g. vortices). In this context, MWC 758 is of particular interest as a protoplanetay disk with two prominent spiral arms observed in scattered light (Grady et al. 2013;Benisty et al. 2015).
Send offprint requests to: A.
Boccaletti, e-mail: anthony.boccaletti@obspm.fr Based on data collected at the European Southern Observatory, Chile under programs 096.C-0241 and 1100.C-0481 While Benisty et al. (2015) show that reproducing the shape of the spiral arms with planets located within the 50 au cavity would require an unphysical hot disk, Dong et al. (2015) found that outer planets with a few Jupiter masses, orbiting at 100 au, could explain the large pitch angle. More recently, to match the spiral patterns in scattered light, as well as the vortices identified in the submillimeter (Dong et al. 2018), Baruteau et al. (2019) proposed two planets of 1.5 and 5 M J , located at 35 and 140 au, respectively. Finding objects responsible for the propagation of density waves has been attempted in the L band to benefit from a lower star-planet contrast and lower dust extinction. Two independent studies, using Keck telescope and the Large Binocular Telescope (LBT), respectively, have reported two distinct point source candidates but with a low significance: one located at ∼ 20 au from the star (Reggiani et al. 2018) that is residing inside the sub-millimeter cavity, while the second (Wagner et al. 2019) is orbiting at ∼ 100 au, thus showing itself to be potentially in line with the Dong et al. (2015) predictions.
The purpose of this letter is to explore the detection limits provided by SPHERE (Beuzit et al. 2019), the high-contrast imaging instrument at the Very Large Telescope, based on two data sets obtained in 2016 and 2018 in the H and Ks bands.  Table 1. Log of SPHERE observations indicating (left to right columns): the date of observations in UT, the filters combination, the amount of field rotation in degrees, the individual exposure time (DIT) in seconds, the total number of exposures, the total exposure time in seconds, the DIMM seeing measured in arcseconds, the correlation time τ 0 in milliseconds, the variation of the flux during the sequence in %, and the true north (TN) offset in degrees.

Observations and data reduction
significantly smaller than the HIPPARCOS distance (279 +94 −58 pc) used in Benisty et al. (2015), and slightly larger (but still within the error bars) than the Gaia DR1 value (151 +9 −8 pc) used by Reggiani et al. (2018).
The observing sequences, and data reduction follow the standard procedure implemented in the SPHERE Data Center as described in Delorme et al. (2017). Further details of the data reduction can be found in several papers, such as Boccaletti et al. (2018). The conditions (seeing, coherence time) were far more appropriate for a high contrast in the second epoch in 2018. The variability of the star's intensity (as measured during point spread function exposures at the beginning and at the end of the sequence) is an indicator of good stability (∼4% in 2018 instead of ∼10% in 2016). Beside, the airmass is rather large (Am>1.5) for this target.
We made use of SpeCal (Galicher et al. 2018) to process the data cubes of both IRDIS and IFS with a variety of angular differential imaging (ADI) techniques. In addition, we used custom routines to perform reference differential imaging (RDI) with another star as a calibrator 1 , or an even simpler algorithm to perform spatial filtering.

Global description
The disk of MWC 758 is quite bright, hence, it is already visible in raw data and stacked images, without any ADI processing (referred as no-ADI). Subtracting the stellar residuals with the coronagraphic images of a reference star (RDI) that were obtained in similar conditions reveals the signal of the disk as close as ∼ 0.1 . Several stars were considered as references observed in the same night in the same configuration, however, here we show the one that provides the best result in terms of the detection of the disk. From these data, we measured F disk /F ≈ 0.017 in the H2 and H3 filters integrating the disk signal between 0.12 1 Other targets observed the same night with the same observing mode, which do not necessarily match the color and magnitude of MWC 758 High-pass filtering images of no-ADI and RDI data at two epochs displayed in a large (5 , top) and narrower field of view (1.5 , bottom). North is up, east is left. The intensity scale is arbitrary. and 0.6 from the star. Figure 1 displays the noADI and RDI images, further processed with high-pass filtering to enhance the disk structures. The upper panels show a 5 field of view, where the previously known background object (Grady et al. 2005) is in fact resolved as two components, with their respective separations and position angles of ρ = 2511.4 ± 1.5mas, PA=316.63±0.03 • , and ρ = 2632±32mas, PA=317.9±0.7 • . The error on the second component is significantly larger because of the brightness ratio (∼3.7 mag) between the two stars and their proximity (∼152 mas), which affect the astrometric procedure. The apparent motion between the two epochs is fully consistent with the proper motion of MWC 758.
The images contain several structures that are labeled in Fig.  2, with the spirals S1 and S2 being the most prominent structures. It is relatively difficult to identify precisely where the spirals actually start from the shortest angular separations, but the total intensity data, especially with RDI, definitely provide a higher signal to noise ratio (S/N) of the < 0.25 inner region compared to the DPI image in Benisty et al. (2015) and Ren et al. (2020).
Starting from the shortest separations, S1 visually starts from the northeast (PA=60 • , ρ = 0.18 , although it could even extend  as close as the coronagraphic mask) and can be traced over one and a half rotations (Fig. 2a). S2 originates from the northwest (PA=-30 • , ρ = 0.22 ) and can be traced up to the west, covering a full rotation in total. While we note that it is difficult to identify exactly where the spirals start in the interior as there may be confusion with the residual starlight pattern around the coronagraphic mask, there is some evidence that the spirals extend even closer in. Thus, the RDI process ( Figure 1d) provides a cleaner vision of the central part than the no-ADI case (Figure 1c). At the other extreme, at low intensity levels, we can also identify the trails of the spirals (Fig. 2b). This is by far the deepest observation of these spiral arms. We also recovered some structures in the southwest of S1 (at PA = 235 − 290 • ), labeled S3, which could be reminiscent of a feature identified by Reggiani et al. (2018). The SPHERE data presented here convey the possibility that S3 and S4 could be connected (S4 being in the trail of S3). This possible link could have been already posited from the L band data in Reggiani et al. (2018) and has already been established by Calcino et al. (2020) based on hydrodynamical simulations. Thus, S4 could also be considered as the bottom side of the disk in which case the angular separation between S1 and S4, given the inclination of 21 • is related to the disk opening angle which we estimated to ∼23 • , indicative of a rather large disk thickness.
While S2 is apparently close to an Archimedean spiral, it is not the case for S1, which appears perturbed in the PA range 60-180 • , with a sort of kink near the location: PA≈ 176 • , ρ ≈ 0.25 , possibly splitting into two components (Fig. 2b). Although the RDI processing has lower impact on the disk morphology than ADI, the exact distribution of the structures at such stellocentric distances remain to be confirmed. We also present the IFS RDI processed images (Fig. A.1) which provide a much lower signal to noise ratio than with IRDIS.
When processed with ADI ( Fig. 3), the spirals are no longer visible because of the small field rotation (<30 • ) which induces a quite strong self-subtraction (Milli et al. 2012). Instead, many structures (consistent with the kink and perturbation seen in RDI) are revealed along the spirals, some of which could be confused with the usual signatures of point sources in ADI. These patterns are consistent across wavelengths between IRDIS and IFS (Fig. 3, b and d for instance). Due to the difference in starlight rejection between the two epochs obtained ∼3 years apart, it is still difficult to perform a thorough analysis of the evolution of these structures, especially within the sub-millimeter cavity (< 0.3 ).

Limits of detection to point sources
The limits of detection are estimated with SpeCal, as explained in Galicher et al. (2018). Here, we made use of the KLIP algorithm (Soummer et al. 2012) which produces a higher contrast when compared to other types of ADI implementations in this particular case. The contrast at a given radius is calculated from the standard deviation of the signal contained in an annulus of 0.5×FWHM in width (about 2 pixels) centered on the star. The self-subtraction induced by ADI is estimated with fake planets injected into the datacube (along a spiral pattern to cover a range of separations and azimuths). The presence of the disk's spirals and their residuals after ADI processing inevitably corrupts this measurement. Both the self-subtraction and the radial transmission of the coronagraph (as reported in Boccaletti et al. 2018) are taken into account to produce the final contrast plots in Fig B. 1 (left). The second epoch is significantly better with, for IRDIS, a gain of about one order of magnitude in contrast, in the separation range of 0.1 − 0.2 (as a result of a sub-optimal coronagraph used in K band and poorer conditions), and a more modest gain of 2 to 4 further out. The gain is less pronounced but still noticeable for the IFS sequences.
Two point sources were identified from previous L band observations. The first companion candidate (CC R ) reported by Reggiani et al. (2018) using Keck is located very close to the star inside the cavity at ρ = 0.112 ± 0.006 , PA = 169 ± 4 • , A&A proofs: manuscript no. mwc758_arxiv  (Fig. B.1, right). The K band data are less constraining as the achieved contrast corresponds to a mass range of 1.9 − 4.5 M J . Therefore, the SPHERE observations do not support the presence of a point source at the location of CC W , nor the detection limits are compatible with the expected mass range. Moreover, it is obvious in Fig. C.1b that CC W is not located along the trace of the spiral arm with a radial departure of about 0.08 that is twice the angular resolution of the SPHERE images. Given the level of clumping that is observed from the innermost separations along the spirals and all the way out (Fig. 3), CC W could be a dust feature (although it should be also visible at shorter wavelengths in reflected light), or an artifact given the low signal-to-noise ratio reported in Wagner et al. (2019).
On the contrary, for dynamical reasons, Reggiani et al. (2018) argue that the emission from the closest point source, CC R , cannot originate from its photosphere or it would translate to a mass of 41 to 64 M J ; hence, in the brown dwarf regime and would have induced a noticeable cavity in the small dust grains distribution. Such a massive companion would have been presumably detected with SPHERE, according to Fig. B.1 (right). Instead, the authors postulate that the emission comes for a circumplanetary disk (CPD), the luminosity of which depends on the product of the planet mass (M p ) and the accretion rate (Ṁ). Assuming face-on 1-Myr old disks and constant accretion rate, Zhu (2015) tabulated the absolute magnitudes in near-IR filters for several values of the accretion disk inner radius R in (from 1 to 4 Jupiter radii). Under this theoretical framework, the photometry of CC R in the L band would correspond to an object of ∼ 1 − 8 M J for a moderate accretion rate of 10 −8 M /yr (consistent with the limit of detection in Hα reported by Huélamo et al. 2018). As a comparison, from the L band photometry we can derive the expected absolute magnitude in the H band for such a CPD according to Zhu (2015, Table 1). We found a value ranging from H=7.7 to 10.3 (R in = 1 − 4), while we measured a contrast of ∆H = 9.57, hence, obtaining an absolute magnitude of H limdet = 10.11. Therefore, the non-detection of a point source in SPHERE data at the position of CC R is only marginally consistent with this flux having been produced by a CPD. In fact, this assumption is even less consistent with a non-detection in the K band since K limdet = 8.78 (∆K = 9.00) and the L band photometry of a CPD should correspond to K = 7.08 − 8.35.

Extinction by the disk
Implicitly, the assessment of detection limits in the previous section assumes the same optical depth at all wavelengths. In a more realistic case, the surrounding disk is expected to produce some extinction and it is necessary to estimate the variation of optical depth between the L band, and the H and Ks bands to derive meaningful conclusions. The optical depth can be assessed using an estimate of the surface density at the position of the companion candidates and a model of the dust opacity. Boehler et al. (2018) considered the disk surface brightness of MWC 758 measured in the continuum with ALMA, and coupled to a radiative transfer model, to derive azimuthally averaged surface density values of the order of 4.4 10 −3 (3.3 10 −3 g/cm 2 resp.) at the position of CC W (CC R resp.). Using a dust model in which the composition and the size distribution are consistent with parameters used by Boehler et al. (2018) 2 , we computed dust mass-opacities in the filters of interest, in the H, K, and L' bands. Combined with the half surface densities estimates (to consider only the amount of dust mass from the disk midplane to the surface) at CC R and CC W locations, we calculated the difference of magnitudes, which are displayed in Tab. 2. However, since water ice is absent from many dust models of protoplanetary disks (Woitke et al. 2019, and this is particularly relevant for the surface layers of disks where most of the opacity in the H/K/L' bands is produced), we also computed these magnitude differences when the water ice component is removed. The differences of magnitudes (H-L' and K-L') can even become negative due to some extinction resonances at specific wavelengths. In the case of the "standard" particle size distribution starting at a min = 0.05 µm, the magnitude differences between the H/K and L' bands are small. In the case of a larger minimum size (a min = 1 µm), which would correspond to a more evolved/processed dust population, we observe more dispersion in the magnitudes differences. Yet, all  Table 2. Magnitude differences with respect to the L' band estimated at the positions of CC R and CC W for two assumptions on the minimum grain size.
the magnitudes extinctions remain below one magnitude, which shows that a large amount of differential extinction is unlikely -an argument with the potential to reconcile present SPHERE observations and previous ones. Furthermore, CC R and CC W inferred companions are massive enough to carve gaps in the disk, which would reduce the dust opacity effects even more. The detection limits derived in the previous section can be revisited given the extinctions estimated in Table 2 and in considering the worst-case scenario (one of the two dust compositions). Assuming the H band detection limit would suffer from differential extinction of 0.83 magnitude (resp., 0.34 magnitude in K band) at the location of CC W for the smallest minimum grain size, the mass limit increases to 0.9 − 2.0 M J (resp., 2.2 − 5.3 M J in K band), which is still lower than (resp., comparable to) the mass inferred by Wagner et al. (2019, 2 − 5 M J ). Larger minimum grain sizes imply either a lower extinctions in H (-0.13 magnitude) or slightly larger in K (0.51 magnitude), which does not, therefore, change drastically the outcome.
In what concerns the object CC R , the extinction moves the limit of sensitivity from H limdet = 10.11, to 9.49 for a min = 0.05 µm, respectively 10.53 for a min = 1 µm. In the K band, the sensitivity reaches 8.52 or 8.39 instead of 8.78 depending on the minimum grain size. Comparing with the extrapolated magnitudes of a CPD in the H band (H = 7.7 − 10.3) and the K band (K = 7.08 − 8.35), these limits of sensitivity are sufficient to rule out -either partially (H band) or totally (K band) -the presence of a CPD.
We note that at the position of CC W and CC R , the local surface density in sub-micronic/micronic particles could be different from that extrapolated from Boehler et al. (2018), with ALMA measurements tracing the larger ones. Also, additional accretion along planet poles could also make the extinction higher than evaluated, especially in a face-on geometry (Fung & Chiang 2016).

Spirals as tracers of point sources
We measured the spine of the two main spirals (S1 and S2) following the method developed in Boccaletti et al. (2013), along with an adequate numerical mask to isolate each spiral (data points are reported in Fig. C.1, left). The spirals are relatively well fitted with an archimedean function (ρ = a × θ + b) in the PA range [−70 • , 330 • ] for S1 (green circles, χ 2 ν = 1.27, a = 0.29 mas/ • , b = 159 mas) and [−30 • , 210 • ] for S2 (blue circles, χ 2 ν = 1.20, a = 1.11 mas/ • , b = 128 mas). Planetary mass objects in the MWC 758 system can be responsible for the launching of density waves materializing into spiral arms in and out of the planet orbital radius. An obvious test would be to check whether we can associate the spirals S1 and S2 with CC R and CC W . However, a thorough investigation will require dedicated hydrodynamical simulations in 3D which is beyond the scope of this paper. Previous works extensively used the linear approximation to the density waves theory (Rafikov 2002) to interpret the geometry of spiral arms observed in scattered light (Muto et al. 2012;Boccaletti et al. 2013;Benisty et al. 2015) although this approach has been found unreliable for massive planets and for cases when the spiral wake is localized away from the disk plane, at the disk surface, as in the case for scattered light images . Nevertheless, the pitch angle of the outer spiral arm is reduced by the disk 3D structure which compensates for the wake broadening. Hence, the linear theory still provides a decent match for the outer spiral even for a massive planet (several Jupiter masses). With these limitations in mind we performed a qualitative analysis following the prescription of Muto et al. (2012): where h c is the disk aspect ratio at the planet location (r c , θ 0 ), α and β are the power-law exponents of the angular frequency and temperature profile dependence with r (in the standard case α = 1.5, β = 0.25). In addition, we assumed the disk inclination and position angle: i = 21 • and PA = 62 • , together with h c = 0.18 (Boehler et al. 2018;Andrews et al. 2011). We also accounted for a standard flaring index γ = 1.2. Imposing the locations of the two candidate companions CC R and CC W , we calculated the associated spirals launched from these sites (respectively the green and blue lines in Fig. C.1, right).
Focusing on the outer arms of the spirals, beyond the planets orbital radii, we found no obvious match between these spiral models and those observed in the data. This may suggest that the perturber(s) could be much less massive or they could be located at smaller physical distances (or both). Again, we caution that the density wave linear theory does not capture the complexity of this system and cannot be taken as a definitive result in this case. A more elaborate modeling is proposed in Calcino et al. (2020), which qualitatively reproduces the main disk features assuming a 10 M J planet at 33.5 au (about 0.2 ) with some eccentricity (e = 0.4).
As long as planet-induced spirals are co-rotating with the planets, we can also consider dynamical arguments to evaluate the amount of expected rotation of the spirals. For an object at the position of CC R , and CC W , we expect an angular motion between the two epochs (2.96 years apart) of ∼ 17.9 • , ∼ 1.3 • respectively, for a central star with a mass of 1.7M (average value in the literature). Therefore, we can safely reject CC R as being the source of the spiral arms, as also derived by Ren et al. (2020). To obtain a more precise estimate of the spiral rotation, we isolated the spiral S2 with a numerical mask in the high-pass filtered RDI images from Fig. 1 (we ignore S1 as being significantly contaminated by diffraction residuals in Jan. 2016 data), deprojected the images, multiplied by the square of the stellocentric distance to give more weight to the outer parts of the spiral, and then we solved for the rotation angle which minimizes the difference between the two epochs with a least square metric. While the exact value depends on the image processing (masking and  Spines of the spirals S1 (green circles) and S2 (blue circles) superimposed on RDI processed image from Dec. 2018 together with the fit of archimedean spirals (left panel). The right panel shows two spiral models with CC R (green line) and CC W (blue line) as the perturbers, assuming linear density wave theory.