Perturbers: SPHERE detection limits to planetary-mass companions in protoplanetary disks

The detection of a wide range of substructures such as rings, cavities and spirals has become a common outcome of high spatial resolution imaging of protoplanetary disks, both in the near-infrared scattered light and in the thermal millimetre continuum emission. The most frequent interpretation of their origin is the presence of planetary-mass companions perturbing the gas and dust distribution in the disk (perturbers), but so far the only bona-fide detection has been the two giant planets around PDS 70. Here, we collect a sample of 15 protoplanetary disks showing substructures in SPHERE scattered light images and present a homogeneous derivation of planet detection limits in these systems. We also estimate the mass of these perturbers through a Hill radius prescription and a comparison to ALMA data. Assuming that one single planet carves each substructure in scattered light, we find that more massive perturbers are needed to create gaps within cavities than rings, and that we might be close to a detection in the cavities of RX J1604, RX J1615, Sz Cha, HD 135344B and HD 34282. We reach typical mass limits in these cavities of 3-10 Mjup. For planets in the gaps between rings, we find that the detection limits of SPHERE are about an order of magnitude away in mass, and that the gaps of PDS 66 and HD 97048 seem to be the most promising structures for planet searches. The proposed presence of massive planets causing spiral features in HD 135344B and HD 36112 are also within SPHERE's reach assuming hot-start models.These results suggest that current detection limits are able to detect hot-start planets in cavities, under the assumption that they are formed by a single perturber located at the centre of the cavity. More realistic planet mass constraints would help to clarify whether this is actually the case, which might point to perturbers not being the only way of creating substructures.


Introduction
Protoplanetary disks (PPDs) are the by-product of the star formation process, and the place where giant planets form before all the gas is accreted onto the star or dispersed over a period of ∼ 3-10 Myr (Fedele et al. 2010). In the last few years, high-resolution observations have opened a new era in our understanding of the gas and dust around young stars. These observations show a plethora of complex substructures in PPDs that are remarkably common when imaged with sufficient angular resolution, including gaps, cavities, rings, vortices, asymmetries, and Article number, page 1 of 26 arXiv:2103.05377v1 [astro-ph.EP] 9 Mar 2021 spiral arms (e.g., van Boekel et al. 2017;Avenhaus et al. 2018;Andrews et al. 2018;Long et al. 2019;Garufi et al. 2020).
The origin of these morphologies remains unclear, and different mechanisms have been proposed to explain them (e.g., Flock et al. 2015;Okuzumi et al. 2016;Takahashi & Inutsuka 2016;Cieza et al. 2016;Gonzalez et al. 2017;Dullemond & Penzlin 2018;Riols et al. 2020). A common interpretation is to describe them as signposts of planetary-mass companions interacting with the disk, which requires a formation of the giant planets and their location within the gaps in less than a few Myrs. In this scenario, the massive planet creates a pressure bump in the gas that stops the radial drift of the dust, which gets trapped in the pressure maxima creating a ring-like structure. In general, mmsized dust is affected by the drag with the gas and it is easily trapped in pressure bumps (Drążkowska et al. 2016;Taki et al. 2016). Small µm-sized dust is however coupled to the gas, following its distribution and possibly populating the gap. As a result, spatial segregation is expected in the distribution of small and large dust particles (Rice et al. 2006;de Juan Ovelar et al. 2013;Pinilla et al. 2015;Hendler et al. 2018). The exact morphology and structure of these regions will eventually depend on planet mass and PPD properties, such as viscosity and temperature (e.g., Whipple 1972;Pinilla et al. 2012).
To understand how these structures are formed and to shed light on the planet formation mechanism, SPHERE highresolution scattered light observations probe the surface layers of the optically thick dust using the polarized differential imaging technique (PDI; de Boer et al. 2020b;van Holstein et al. 2020). This mode detects the polarized light scattered off µmsized grains, and is thus an effective way of removing the stellar halo without altering the underlying disk morphology with postprocessing techniques. To this date, there are ∼ 90 disks imaged in this mode, from individual large disks around intermediate mass stars (e.g., Benisty et al. 2015;Ginski et al. 2016;Stolker et al. 2016;Pohl et al. 2017) to surveys designed to alleviate observational biases (DARTTS; Avenhaus et al. 2018). These results show a ubiquity of substructures in scattered light in IR bright disks, unless the disks are small or too faint.
The indirect detection of protoplanets creating these observed PPD substructures is particularly hindered by the presence of circumstellar material, while stellar activity also limits the performances of radial velocity surveys. As of today, direct imaging is the best technique to detect and characterize at the same time the planet and the disk morphology, although the companion luminosity might also be affected by extinction due to the surrounding disk material (Szulágyi et al. 2018Szulágyi & Garufi 2019;Sanchis et al. 2020). So far, besides the exception of the two planets around PDS 70 Müller et al. 2018;Haffert et al. 2019;Mesa et al. 2019), only not-confirmed planetary-mass candidates have been found with direct imaging in PPDs. A few kinematic detections of embedded protoplanets have also been proposed (Pinte et al. 2018;Teague et al. 2018;Pinte et al. 2019), but there is not yet a final confirmation of these objects either.
Detection limits to low-mass companions via total intensity high-contrast imaging data are available for only a handful of targets that show hints of substructures. These data, however, are analysed with different algorithms that remove the stellar halo, using different mass-luminosity relationships, and sometimes outdated and heterogeneous stellar parameters.
In this paper we present a homogeneous study to determine planet detection limits within the substructures of a sample of PPDs observed in µm-sized dust with SPHERE. Our goal is to perform a systematic analysis to planet sensitivities in these sys-tems, and to understand how far off we are from detecting perturbers with direct imaging. In Section 2 we motivate the PPD sample and present the high contrast imaging reduction of the data. In Section 3 we convert the obtained detection limits to mass sensitivities, and assess the effect of different formation luminosities on the detectability of low-mass companions. In Section 4 we estimate the mass and location of the potential companions carving these substructures, and compare them to the SPHERE detection limits. Finally, in Section 5 we discuss the results and caveats of this paper, and in Section 6 we present the conclusions.

Sample selection and data reduction
We collected PPD systems observed with SPHERE that comply with the following criteria: 1. The PPD shows rings, cavities or spirals in SPHERE/PDI observations. 2. The host star is single or has no stellar companions closer than 3 . 3. The PPD has an inclination below ∼ 60 deg. 4. The PPD has been observed in coronagraphic angular differential imaging mode (ADI, Marois et al. 2006) with SPHERE, suitable for the detection of low-mass companions.
The resulting list of targets includes disks that have substructures observed in the small dust component, and allows us to obtain sensitivities to planets that are not critically affected by projection effects and extinction towards the line of sight (e.g., Szulágyi & Garufi 2019). Moreover, our selection removes close binaries to optimize that the physical and morphological properties of the PPDs are not influenced by the presence of a stellar companion; this allows us to attribute the substructures to intrinsic physical processes related to the disk or host star (e.g., Durisen et al. 2007), or otherwise to the presence of planetarymass objects.
The target selection was mostly based on the compilation by Garufi et al. (2018) and the latest results of the DARTTS-S survey (Avenhaus et al. 2018;Garufi et al. 2020), followed by a cross-match with the ESO archive of ADI data. The final sample is comprised of 15 PPDs, whose substructures are characterized based on their visual appearance in the PDI images. Garufi et al. (2018) classified them in Spirals, Ring, Rim, Giant, Faint, Small and Inclined. This classification is rather subjective and can be affected by projection effects and the observing conditions of the different datasets; this implies that some PPDs will have some overlap, i.e., rims and/or rings (e.g., PDS 66), or giants with spirals (e.g. HD 100546). In this section we further simplify the classification and assign one class to each PPD depending on which feature is most prominent; Spiral when spiral arms in the µm-sized dust are seen, Rings when the PPD shows signs of a series of resolved rings with gaps in between those, and Rim when there is a detection of a large central cavity with a clear bright rim (typically known as transition disks). Some systems reveal different substructures when observed with ALMA in the gas or mm-dust (e.g., Teague et al. 2019), but we do not consider those images for the classification. HD 100546 and HD 34282 are complicated systems with the presence of arcs, spirals and cavities. Here we put them in the Spiral group, but the presence of rims or rings in these systems will also be taken into account (see Appendix A). For instance, HD 34282 hosts an inner cavity with a rim-like structure at ∼ 88 au, and a potential single-armed  (Choi et al. 2016). The resulting stellar parameters are shown in Table 1 and have been derived as explained in Section 2.1 (Right) Mass-distance diagram of the same population. spiral feature farther out better observed after deprojection (de Boer et al. 2020a). These individual substructures within a given PPD (not the global classification made in this section) will be taken into account separately to derive potential planet masses and locations in Section 4.

Stellar parameters
In Table 1 we show the derived parameters for the stellar hosts present in our sample. We follow the method outlined in Garufi et al. (2018), but in this case we obtain updated values using the new Gaia Early Data Release 3 parallaxes (Gaia ED3, Gaia Collaboration et al. 2016. We build the stellar SEDs using VizieR and fit the wavelength range 0.4 µm -1.3 mm to the PHOENIX photospheric models (Hauschildt et al. 1999), which allows us to calculate the stellar luminosity using literature extinction A V and effective temperature T eff . On the left panel of Figure 1 we show a Hertzsprung-Russell diagram of the sample, where we estimate each individual mass and age using the MIST pre-main-sequence (PMS) tracks (Choi et al. 2016). These isochrones are consistent with other PMS tracks such as Parsec (Bressan et al. 2012), Baraffe (Baraffe et al. 2015) and Dartmouth (Dotter et al. 2008). To calculate uncertainties on the age and mass, the error bars on the luminosity (in turn derived from a 20 % uncertainty on the optical extinction) and on T eff (± 200 K) is propagated through the tracks. The right panel shows the location of the targets in the stellar mass-distance space. Distances to the targets range between 60-300 pc, but most of them lie around 150 pc, as members of star-forming regions such as Sco-Cen and Taurus. Spirals tend to appear in disks around massive stars, while objects classified as Rim, with prominent cavities, are only resolved in the scattered light in < 1.5 M Sun stellar sys-tems (although some of the most massive objects classified as Spiral also count with resolved cavities, see Appendix A).
We note here that the derived isochronal parameters for young PMS stars carry a moderate uncertainty. Stellar ages is the most critical parameter to estimate reliable sensitivities to planet mass; deviations in L/T eff values, the use of different evolutionary tracks, the effect of magnetic fields and the initial position of the star at t = 0 can all contribute to dubious estimations of individual ages (e.g., Asensio-Torres et al. 2019). Our approach allows us to obtain a homogeneous and consistent classification of the different PPDs, while the adopted ages (Table 1) also tend to agree within error bars with other individual studies. For instance, PDS 70 is found to have an age of 5.4 Myr , and TW Hya is most consistent with 8 Myr (Sokal et al. 2018). HD 139614 has an age of 10.75 ± 0.77 Myr, recently derived via astereoseismology (Murphy et al. 2020), in close agreement with our values. For HD 169142, an estimation of 6 +6 −3 Myr was found for the M-type wide binary companion 2M1824 (Grady et al. 2007), for which only the upper uncertainty value is consistent with our estimation, although isochronal ages for low-mass stars seem to be underestimated by a factor of ∼ 2 (Asensio-Torres et al. 2019).

Reduction of the pupil-tracking data
For all the objects in the target list of Table 1, we performed high-contrast imaging ADI reductions to obtain detection limits to the presence of perturbers creating the substructures in the small dust distribution. Table 2 shows the corresponding pupiltracking epochs. When more than one ADI observation of the same system was available, we analysed the epoch that provided a better contrast based on published results, collected sky rota- Notes. Stellar parameters of the PPD hosts considered in this work. References for the stellar temperature and optical extinction can be found in tion and weather conditions (see Table 3). The ADI observing modes are split about equally between IRDIFS and IRDIFS-EXT; IRDIFS was the preferred mode by the SpHere INfrared Exoplanets survey (SHINE; e.g., Chauvin et al. 2017), which mainly targetted 10-200 Myr-old nearby stars. This mode results in dual-band imaging in the H band with the infrared dual-band imager and spectrograph (IRDIS; Dohlen et al. 2008), and Y J spectrophotometry with the near-infrared integral field spectrograph (IFS; Claudi et al. 2008) in the Y-J range with a FoV of 1.7 × 1.7 . With the inclusion of younger stars located beyond 100 pc in the survey (particularly from Sco-Cen) or dedicated follow-up of PPDs, the IRDIFS-EXT mode was used, which extends the dual-band imaging to the K band and the IFS spectrophotometry to H band, as this mode is more adapted for the detection of young red L-type planets.
The reduction and analysis have been conducted homogeneously for the entire dataset. For the IRDIS side we first used the vlt-sphere repository (Vigan 2020) 1 , which is a python-based pipeline to pre-reduce SPHERE data (Beuzit et al. 2019). We sorted out the files and performed static calibrations, such as flat field reduction, dark and sky subtraction, and bad pixel correction. This was followed by star registration, for which we used the waffle pattern to find the location of the host star behind the coronagraphic mask, and flux calibration to get the difference in flux between the unsaturated star and the science frames. After this pre-reduction step we obtained a cube of centred saturated frames, collapsed unsaturated stellar frames for flux calibration, a sky rotation file and the wavelength in which the observations were carried out. To treat the quasi-static speckle noise affecting high-contrast imaging observations, further PSFremoval was conducted with the ANgular DiffeRential OptiMal Exoplanet Detection Algorithm (ANDROMEDA; Cantalloube et al. 2015). This algorithm is based on the maximum likelihood estimation. It first attenuates extended structures by applying a high-pass Fourier filtering that removes 25 % of the lowest frequencies (with an energy loss of ∼ 18 %, see Fig. 1 in Cantalloube et al. 2015). Subsequently, ANDROMEDA performs a model cross-correlation of the signature that a point source leaves after two frames at different rotation angles are subtracted. This processing is thus sensitive to the presence of point sources in the field of view. Before feeding the vlt-sphere output to AN-DROMEDA, we performed a frame selection to remove those that deviate more than n standard deviations of the mean of the 3D cube in a 20 px radius centred on the star. The process is repeated for n=1,2,3... until the rejection fraction is lower than 20 %. Finally, ANDROMEDA provides two 2D maps: (i) the estimated contrast that a point source would have at every location in the image, and (ii) the corresponding uncertainty on this estimated contrast. The limiting contrast is therefore given by the uncertainty map, previously normalised empirically to account for the non-gaussian residual noise at small separation (Cantalloube et al. 2015). This map can be used to derive a 1D 5σ projected detection limit by taking its azimuthal median.
IFS data was pre-reduced via the SPHERE Data Center (Delorme et al. 2017), and the resulting cubes went through the the SpeCal software (Galicher et al. 2018), which further removes the stellar halo by applying an ADI-based method. We selected an Angular and Spectral Differential Imaging (ASDI) reduction following the principal component analysis approach in Mesa et al. (2015) for the spatial and spectral channels. This method is aggressive in terms of speckle subtraction and seems to reach better contrasts than IRDIS-ADI within the IFS FoV (Langlois et al. 2021). Moreover, 1D IFS-ASDI projected contrast curves were computed as the azimuthal standard deviation of the noise residuals, corrected for self-subtraction via fake planet injections. To obtain 2D IFS-ASDI limiting magnitude maps that ac-  count azimuthally for the presence of disk residuals, we evalu-ated, around each pixel, the standard deviation within a ring of inner and outer radius of 0.7 and 2.5 FWHM, respectively.

Detection limits
In Figure 2 we show the projected 1D contrast limits achieved by ANDROMEDA/IRDIS and SpeCal-ASDI/IFS for the entire PPD sample. In general, ANDROMEDA/IRDIS achieves magnitude differences of ∼ 11-14 mag at 0.5 in the contrast-limited regime, and down to 16 mag at larger separations limited by background sensitivity. With the use of ASDI, IFS contrasts seem to improve the IRDIS limits at close separations, within 0.4 . Besides the presence of disk residuals, variations in Strehl ratio, weather conditions, collected sky rotations and magnitudes of the host stars all contribute to the final detection limits.
To convert these contrast to detectable companion masses, a post-formation luminosity of the perturber needs to be assumed. This is known as the 'initial luminosity' with respect to the cooling phase, and is set by the radiative transfer and thermodynamics during the accretion phase (e.g., Mordasini et al. 2017;Szulágyi & Mordasini 2017;Marleau et al. 2019a). These differences can lead to divergences between the hot (more luminous) and cold (fainter) starts of about 3 orders of magnitude in luminosity during the first few Myrs (e.g., Marley et al. 2007;Spiegel & Burrows 2012). Although it is likely that a smooth range in planet formation luminosities between the extreme hot and cold models exist, the coldest starts do not seem to be a valid representation of the physical conditions of the current population of imaged planets; for instance, dynamical masses of the two planets around β Pic seem to fall close to the predictions of a hot start, either formed through disk instability or a hot core accretion process (Brandt et al. 2020;Nowak et al. 2020 four giant planets in the HR 8799 system (Wang et al. 2018). Although these results could be affected by an observational bias, where we simply have not yet detected cold-start planets because they are fainter (see Section 5.1), theoretical studies of the accretion shock also favour hotter start models (Marleau et al. 2019b).
Here we rely on the AMES-DUSTY atmospheric models (Chabrier et al. 2000) to reproduce the hottest, most luminous outcome of planet formation, usually associated with disk instability. This model includes dust absorption and scattering, since young planets seem to have condensates in their atmospheres (e.g., Müller et al. 2018). AMES-DUSTY has been widely used by the community, and starts with a completely formed planet at arbitrary initial entropy for a given mass.
We also use the Bern EXoplanet cooling curves (BEX), which are based on the population synthesis models by Mordasini et al. (2017). They provide luminosities of young giant planets formed through the core accretion paradigm. Depending on whether the accretion shock luminosity is radiated away or deposited into the planet, the initial luminosities of the newborn planets are classified in diminishing order of brightness from 'hottest' to 'coldest'. These are then expanded into evolutionary tracks at constant mass, reproducing the cooling under different atmospheric conditions (Marleau et al. 2019a). We make use of the BEX-COND tracks, where the initial planet luminosities are coupled to the boundary conditions given by the condensatefree AMES-COND atmospheric models (Allard et al. 2001), and use these resulting BEX-COND gravities and temperatures as a function of time to calculate the magnitudes via DUSTY. We adopt the 'hot' and 'warm'-start initial conditions from BEX-COND, which we call BEX-HOT and BEX-WARM. These two relations correspond to a fit extension of the cold nominal (BEX-HOT) and cold classical (BEX-WARM) population in Marleau et al. (2019a); see Equation 1 in Vigan et al. (2020).
On the right panel of Figure 2, we see the effect of these different post-formation luminosities on the median projected sensitivity to planet masses around a stellar host located at the typical distance of our sample (145 pc). Down to ∼ 15 au (0.1 ), assuming that the AMES-DUSTY models are valid, the SPHERE observations would be sensitive to planets below 10 M Jup , and down to 4 M Jup at > 0.5 . This is different by a factor ∼ × 1.5 if the BEX-HOT start models are considered, and planets below 11 M Jup would not be detectable if the cold BEX-WARM models are a good representation of planet formation.
There are however a few caveats in the use of evolutionary models. If the perturber is formed via core accretion, there will be a delay time with respect to the formation of its host star (e.g., Fortney et al. 2005;Bonnefoy et al. 2014). For gas giants the delay cannot be longer than the typical PPD lifetime of ∼ 3 Myr (e.g. Ribas et al. 2015); most giant planets around an A-type star will have acquired their final mass after 2 Myr (Williams & Cieza 2011;Marleau et al. 2019a). Moreover, if the planet is still forming, its intrinsic luminosity could be augmented by accretion onto the planet surface and/or shock and thermal emission from a potential circumplanetary disk (CPD) (Szulágyi et al. 2014 Aoyama et al. 2021), while at the same time also reddened by the presence of material in the circumplanetary environment (Szulágyi et al. 2018Szulágyi & Ercolano 2020). In these cases where formation is ongoing, these evolutionary models of post-formation luminosities are probably not a very solid representation of the physical parameters of the perturbers. However, without comprehensive spectroscopic analysis it is extremely difficult to derive the mass of a forming protoplanet from the expected magnitudes, as the contribution of accretion can be very important, creating a spread in luminosity of up to ∼ 2 dex for a given planet mass (see Figure  . 5σ mass sensitivity to companions in the PDS 66 system and the IRDIS H2 filter. The removal of the stellar halo has been performed with ANDROMEDA, and AMES-DUSTY initial conditions have been assumed for the perturber's luminosity. The black dashed circumferences show the outer radius of the bright compact region at 25 au and the location of the outer ring at 85 au, respectively, as seen in scattered light by Avenhaus et al. (2018). The white thick circumference shows the proposed location of a perturber, at 55 au, creating the ring-like substructure, which coincidentally corresponds to a further depleted region in the scattered light data. The black circle shows the location of the star and has a radius of 0.1 . The 2D map has been deprojected based on the disk parameters of Table 1. dasini et al. (2017)). Here we suppose that the perturber starts its post-formation cooling at the same time as the star reaches its pre-main sequence phase, and assume that any potential time delay would be included within the estimated error bars on the age of the host star. Extinction effects are discussed in Section 5.2.
Finally, we convert the projected 2D limiting contrasts to mass limits for each pixel in the image. These are then used to (i) obtain a detection probability map for each target (included in Appendix C) via the Multi-purpose Exoplanet Simulation System (MESS, Bonavita et al. 2012) and (ii) obtain 1D deprojected mass sensitivities. For the latter, we deprojected the 2D maps by inclination and position angle to be in the disk plane (see an ANDROMEDA/IRDIS example in Figure 3, and in 4 for the ASDI/IFS mode), and at each radial separation we take the median of the limiting mass values to construct 1D sensitivities, which are shown for each individual system in Appendix B. We also overplot the location of the substructures seen in scattered light and in the mm continuum as explained in Appendix A. In the following section we will use this method to constrain the detectability of the perturbers that may be creating the PPD substructures.

Population of perturbers in SPHERE PPDs
To approximate the mass and location of the perturbers causing the substructures in the PPDs of our sample, we first rely on the simulated results in the literature where gap morphology in the µm-sized dust is taken into account   . 5σ mass sensitivity to companions in the HD 135344B system achieved by the ASDI/IFS reduction. As in Figure 3, AMES-DUSTY initial conditions have been assumed for the perturber's luminosity, and the 2D map has been deprojected based on the disk parameters of Table  1. Disk residuals caused by the spiral pattern are clearly seen. The black dashed circumference shows the location of the cavity in scattered light as observed by SPHERE/PDI . The inner red thick circumference shows the proposed location of a perturber in the middle of the cavity, very close to the coronagraphic mask, and the outer one represents the proposed location of a 5-10 M Jup companion creating the spirals (Dong & Fung 2017a Baruteau et al. 2019). However, these derivations make up only a subset of the substructures seen in our target list. Given the intricate disk-planet interaction and the large numbers of gap structures we are dealing with, in addition to the simulated perturbers obtained from the literature we also consider two different ways of estimating their masses: 1) Gap width proportional to the planet Hill radius and 2) Difference in cavity size in near infrared scattered light images and ALMA mm continuum.

Gap Opening Mass from Hill Radii
The gap width that a planet carves in the disk can be approximated by a proportionality factor applied to its Hill sphere of influence, defined as R Hill = R p (M p /3M ) 1/3 . A simple estimation that can relate the observed width in scattered light to the perturber's mass is given by: where the width is defined as the distance between the gap minimum R PDI,gap and the peak of the ring R PDI,peak seen in PDI, and k is a scaling factor. This formulation is similar to the latest estimations of gap-carving planet masses in the millimetre continuum (e.g., Lodato et al. 2019), which are in agreement with the individually-modelled DSHARP planet population and the broader collection from Bae et al. (2018).
To estimate the scaling factor in the µm-sized dust, we consider the gap located at ∼ 127 au in the PPD around HD 97048. The edges of this gap are well resolved in SPHERE/PDI and total intensity, and the surface brightness of the bottom is detected at 5σ background sensitivity (Ginski et al. 2016). First, we apply the analytical gap-opening criterion from Crida et al. (2006), which does not take into account gap morphology nor dust dynamics and evolution. This solution gives the minimum planet mass that clears 90 % of the gas disk when the following criterion is met: where H is the local scale height and Re the Reynolds number defined as r p 2Ω p /ν, with r p being the radius of the planet orbit, Ω p the angular velocity and ν the viscosity ν = α c 2 s /Ω p . The shaded area of Figure 5 shows the allowed mass values according to this approach for α = 10 −3 and assuming H/r p = 0.056 (same value as found by Dong & Fung (2017b) at the gap location through Equation 4, see later). A perturber with a mass ≥ 1.3 M Jup is necessary to create a well-depleted gap in the gas.
To take into consideration the morphological features that a planet causes in the disk, hydrodynamical simulations have investigated the gap parameters that a planet embedded in a PPD with a given physical properties would impose on the gas (e.g., Fung et al. 2014;Duffell 2015;Kanagawa et al. 2015Kanagawa et al. , 2016. The associated analytical criteria approximate these simulations and link planet masses to gap depth and width. Following Dong & Fung (2017b), the gap depth in gas is given by: and the normalised and dynamical gap width is found to scale only with the local disk aspect ratio as: where q, α and H/r are the mass ratio between the companion and the host star, the accretion disk viscosity parameter (Shakura & Sunyaev 1973), and the disk aspect ratio, respectively. The locations (r in, and r out, ) are the edges of the gap, defined as the positions where the surface density reaches the mean of the minimum (r min, ) and undepleted 0 (r min, ) surface density. These morphological parameters also need to be connected to the observed ones in scattered-light images. 2D and 3D hydrodynamical simulations coupled to radiative transfer calculations show that this translation depends mostly on the inclination of the disk and the angular resolution and sensitivity of the observations (Szulágyi et al. 2018). Using this formulation, Dong & Fung (2017b) obtained a value for the mass of the perturber of 1.3 M Jup for α = 10 −3 and H/r = 0.056 (see Figure 5). This is in agreement with the derivation using Crida's criterion, although it does not always have to be the case, as equations 3 and 4 are accommodated to each specific gap morphology, unlike equation 2. When comparing the two criteria at the depth assumed by Crida, this later prescription tends to derive a higher planet mass (see, e.g., Figure 12 in Bergez-Casalou et al. 2020).
We finally apply the gap-opening Hill radius criterion of equation 1 to this gap. A scaling factor of k = 8 is necessary to simulate the mass derived by Dong & Fung (2017b). We find that for LkCa 15, the ∼90 au gap in TW Hya and HD 169142, we are also able to approximate the simulations of Dong & Fung (2017b) with k ∼ 8-10. A k >10 scaling is however necessary to account for their derived masses in the ∼ 78 au gap around RX J1615.3-3255 and the ∼ 22 au gap in TW Hya, which is related to the individual gap morphology that the Hill radius approximation does not account for. To give an insight on the uncertainties involved when applying this basic Hill radius approximation, variations in the k factor from 8 to 12 can be translated to a perturber's mass estimate between ∼1.5 to 0.5 M Jup for the gap we are considering in HD 97048. Likewise, another source of uncertainty is the viscosity α assumed by the models; Dong & Fung (2017a) find that the mass of the perturber changes by a factor of 10 between the two extreme values they considered, from α = 10 −4 to α = 10 −2 .
Here we take from the literature the location of all the gaps between rings and within cavities found in the SPHERE scattered light images of our PPD sample, and assume that one planet at that location opens a gap that scales with k = 8. This also includes PDS 70, for which we ignore the two giant planets detected in its cavity for the rest of the analysis in order to be used as a point of reference and comparison with other systems. An individual brief discussion of each system and its substructures is provided in Appendix A, together with a very simple representation of their scattered light appearance and the locations at which we locate the perturbers (Figure A.1). If the gap minimum location is not explicitly calculated in the literature, we place the putative planet in the centre between a pair of rings. When the region between the star and the inner ring is depleted down to the coronagraphic IWA in the µm-sized dust, we treat it as a cavity. In these cases, we locate the planet at half the distance between the star and the rim, which does not necessarily need to be the case in a real scenario. We also treat the outer rings of PDS 70 and LkCa 15 as cavities, although non-resolved inner disk signals in the small dust are detected close to the mask. With this prescription we seek to account for a global estimate of planet masses in real observations of non fully-depleted gaps.  6. Underlying population of potential perturbers derived from PPD substructures in the SPHERE scattered light images. For planets in the gaps between rings (circles with inner dots) and planets in cavities (circles), these values are estimated using a proportionality factor of 8 R Hill to account for the mass that would create the observed width in the scattered light gap (Equation 1). When the gap morphology in scattered light was taken into account, planet masses are taken from detailed simulations in the literature for α = 10 −3 , and are marked with an asterisk. Simulated planets creating spiral features are also included as spiral symbols (see references in the Appendix A and see also Figure A.1). The distribution of confirmed planets in the exoplanets.eu catalogue are shown as green (radial velocity discoveries), orange (transits) and blue (imaging) small dots. Figure 6 shows the planet masses derived through Hill radius and k = 8. The perturbers are located at the gap minimum within the cavities and in between rings of the PPDs in our sample. When derivations from simulations exist in the literature, we adopt those ones. They are marked with an asterisk in the figure. We also include the latest simulated planet masses and locations that explain the spirals seen in HD 36112  and HD 135344B (Dong & Fung 2017a). We do not consider in this section three systems: CQ Tau, HD 100546 and HD 139614, because they have complicated structures with no clear gaps, or no simulated planets creating the spirals exist. From the figure, two different populations of perturbers seem to arise; more massive 3-10 M Jup objects that create big cavities, and a more sparse population of planets creating gaps within rings in the ∼ 0.01-1 M Jup regime. The latter are all massive enough to potentially create detectable gaps in PPDs (e.g., Bitsch et al. 2018). We note that these derived perturbers are objects that could migrate and evolve (e.g., Lodato et al. 2019), and so the mass-position di-  agram of Figure 6 is not directly comparable to the population observed by the indirect methods. In Figure 7 these masses are compared to the derived 5σ mass detection limits at the proposed location of the perturber, computed for the different starting luminosities after deprojecting the 2D mass sensitivity maps. These two figures suggest that:

Underlying population of perturbers
-Our sensitivities do not reach masses for planets creating the gaps in the disk population with rings.
-Planet masses needed to create the observed cavities are in general larger than those in the gaps between rings, unless the assumption that only one planet induces the cavity is incorrect. -If hot-start models are a good representation of the initial planet luminosities, we could be close to a detection in several systems, omitting the extinction due to the potential surrounding circumplanetary and circumstellar material. So far, two giant planets in PDS 70 are the only bona-fide detections, which might imply the presence of currently undetected lower-mass companions in these systems, or that other physical mechanisms may be responsible for the formation of the cavity, such as dead-zones, MHD winds, instabilities or grain growth at snow-lines (Flock et al. 2015;Pinilla et al. 2016;Takahashi & Inutsuka 2016;Cieza et al. 2016).
-Planet masses creating gaps within scattered-light rings are consistent with the population of undetected companions that is emerging from the high-resolution ALMA data (see Figure  21 in Zhang et al. 2018). Figure 8 shows a histogram version of the results in Figure 7. We can see that there are 8 cases (i.e., the first bin) where these putative planets with hottest starts could be detectable; 6 in cavities (HD 34282, Sz Cha, HD 135344B, RX J1604, PDS 70 and RX J1615) and 2 creating spiral patterns (HD 135344B and HD 36112). Additionally, 3 perturbers might be close to a detection (defined here as less than half an order of magnitude away in mass, forming the second bin) in the cavity of HD 169142 and the gaps between rings in HD 97048 and PDS 66. Table 4 shows these most promising systems where to look for perturbers based on this analysis. For the perturbers in cavities very close to the coronagraphic mask, spatial resolution and projection effects can also affect the detection. Overall, these results might point towards the fact that planets may not be the only way of creating substructures, given the current number of detections, unless cold initial conditions represent a better description of their luminosity, which is unlikely according to Marleau et al. (2019b).

Gap opening mass from scattered light-mm continuum comparison
In this section we derive planet masses for the systems in which resolved rings and/or cavities exist both in SPHERE scattered light and in ALMA mm continuum observations (see Table 5). SPHERE PDI observations trace the small dust coupled to the gas in a PPD, while ALMA data probe mm-sized dust grains settled in the disk midplane. For this reason, the pressure bump created by the presence of a planetary-mass companion can leave its imprint in the different appearance of the disk at different wavelengths. In general, large grains will accumulate more effectively at the peak of the gas pressure, as no drag force exists in the pressure bump, making the gas to rotate near or exactly at Keplerian velocity, while µm-sized dust is well-coupled to the gas and moves along with it. This effect creates a spatial difference in the cavity sizes observed in the mm continuum with ALMA and in polarised scattered light, which is dependent on the planet mass and disk location. de Juan Ovelar et al. (2013) simulated the radial distribution of the dust due to the presence of 1, 9, and 15 M Jup giant planets located at 20, 40 and 60 au around a solar-mass star, after 3 Myr of evolution and α viscosity of 10 −3 . We use their functional form of the ratio in the cavity size seen in µm and mm dust for a given planet mass: with f (M p ) defined as the ratio between the location of the scattered light inner edge of the cavity seen with SPHERE/ZIMPOL, and the peak in 850 µm flux with ALMA, c ∼ 0.85 and Γ ∼ [−0.22, −0.18, −0.16] for planet orbital radii R p = [20, 40, 60], respectively.
The location of the scattered light inner edge is usually not provided from the observations. Following Villenave et al.
(2019), we derive for each gap a minimum planet mass (M p (min)), assuming that the scattered light edge is at the half distance between the peak of the disk in PDI and the gap minimum, and a maximum planet mass (M p (max)), taking the scattered light edge as the position of the peak in PDI. We use as R p the location of the gap in scattered light. In the cases where f (M p ) is close to 1, the predicted mass is below 1 M Jup . This regime is outside the range of the simulations, and observational uncertainties in the derivation of the peak locations would be of critical importance. We thus limit ourselves to the 1-15 M Jup mass range, and treat these values as the upper and lower thresholds of the resulting planet masses outside these limits. We also extrapolate the value of Γ to the location of the gaps in our sample, assuming that it behaves linearly at distances >60 au within the 1-15 M Jup mass range.
From this approximation, we see in Figure 9 that perturbers within large cavities are close to detection under AMES-DUSTY conditions (RX J1604, PDS 70 and HD 135344B). Again, planets in the gaps between a pair of rings such as HD 97048 and HD 169142 seem to be at least about half an order of magnitude lower in mass than currently detectable. The derived masses for these systems are also in good agreement to those derived in the last section from Hill radii or estimated from simulations that relate the mass to the morphological features of the gaps (see Tables 4 and 5). In Figure 10 we show this difference for those substructures carved by 1-15 M Jup planets. For most of these perturbers, the mass derived using the ALMA observations is typically within a factor of ∼ 0.4 and 1.5 of that obtained via the Hill radius criterion. We do not show here the potential planet creating the gap at ∼ 45 au around HD 97048, as only an upper limit can be constrained with Equation 2, but both approaches yield a planet mass well below 1 M Jup . Notes. Perturbers that, according to the estimations of Figure 7, are the most promising to be found. PDS 70 treated as planet-less here. Planet masses in gaps and creating the spiral patterns are taken as described in Figure 6. The columns AMES-DUSTY and BEX-HOT are our derived SPHERE detection limits (best of ANDROMEDA/IRDIS and ASDI/IFS) using the AMES-DUSTY and BEX-HOT models (see Section 3).
References for the locations of the perturbers can be found in Appendix A. See also Figure A.1 for a sketch of these morphologies. Notes. Minimum and maximum mass of the perturbers, M p (min) and M p (max), respectively, calculated using Eq. 5. The derived detection limits at the perturber's location using AMES-DUSTY and BEX-HOT tracks (best of ANDROMEDA/IRDIS and ASDI/IFS) are shown in the last two columns. References for the location of the substructures can be found in Appendix A.

Why is the current detection rate so low?
These homogeneously-derived detection limits in mass have clear implications on the detectability of potential companions in very young PPDs. Indeed, our results suggest that, if planets creating substructures in PPDs have a very hot start, the majority of them are probably less massive than ∼ 4 M Jup (see Figure 2). The mass detection limits are however affected by a factor ∼ 3 if planets were formed via cold core-accretion models, which is a consequence of the very young age of our sample. If we compare our sensitivities to the SPHERE/SHINE survey (see  (Gratton et al. 2019), but remain uncertain (e.g., Rameau et al. 2017;Currie et al. 2019).
An explanation for the low detection rate of giant planets in PPDs could be that the perturbers' masses are overestimated. Indeed, the problem of estimating masses from the observed morphologies in the dust is still complex; current analytical equations describing the gap can only be used after making several assumptions, such as the disk temperature, viscosity, gas accretion effects, dust evolution or migration. Simulations of planetdisk interactions also need to account for more complex and realistic thermodynamics that account for the cooling efficiency of the disk, which would impact the morphology and the number of gaps that a planet can produce (e.g, Szulágyi 2017; Miranda & Rafikov 2020). Including more physics in the theoretical simulations in fact suggests that the derived masses by locally isothermal relations might be underestimated. Under the isothermal assumption, for a given planet mass the gap is wider than in the case where the disk does not cool rapidly. This implies that the derived masses by current isothermal relations might be underestimating the real mass of the perturbers, as even small planets could open up wide gaps. The same conclusion is found if gas accretion is taken into account (Bergez-Casalou et al. 2020).
If the planet is not on a circular or coplanar orbit, the dynamical mass estimates (as Equations 3 and 4) would also be mod- ified, e.g., a planet on an eccentric orbit may open a wider but shallower gap (Muley et al. 2019). This assumption on the planet orbit is mostly justified on the basis of parameter restriction on the theoretical simulations, and it is not clear how this could affect when applied to a large sample. For instance, PDS 70b is found to be on an eccentric orbit of e ∼ 0.17 (Wang et al. 2021), and the presence of misaligned inner disks, such the one around J1604 in this sample, might point to the presence of a companion in a highly-inclined orbit (Mayama et al. 2012;Pinilla et al. 2018). Nevertheless, our masses derived via gap morphology and those calculated with ALMA data seem to agree, and are probably representative of the order of magnitude of the perturbers' mass. For the special case of PDS 70, which we treated as a planet-less PPD, the derived mass of a single perturber carving its cavity at 27 au has a mass of 4.9 M Jup using the Hill radii prescription, or between 2.1-8.5 M Jup following equation 5. This is comparable to the mass of PDS 70b or c (both of ∼ 5 M Jup ) at 23 and 30 au, respectively (Mesa et al. 2019). If accretion luminosity is taken into account as a contributor to the observed luminosity, the mass of PDS 70b decreases to about 1 M Jup (Stolker et al. 2020). Assuming that PDS 70c follows a similar pattern, the combined mass of both planets will be in the range 2-3 M Jup . This fact may indeed be affecting the current number of detections in cavities; as claimed by Dodson-Robinson & Salyk (2011), several undetected lower-mass planets of mass ≥ 1 M Jup might be creating the large holes in optically thin transition disks, instead of a single very massive compan-ion. Assuming the presence of more than one perturber inside of a cavity would greatly affect the mass estimates of Equation 1, which lead to the assumption that Rim systems have more massive planets than Ring systems, with the only exception of LkCa 15. There is not enough information in scattered light yet to derive the mass of multi-planet systems in cavities from their observed morphology, but simulated scattered light images of cavities opened by multiple Jovian planets seem to resemble the observations (see Figures 5 and 7 in Dong et al. 2015).
Another possibility is that some (or most) of these planets reradiate the energy influx generated during the gas accretion shock, and consequently are better described by cold starts, which is not theoretically supported (Szulágyi & Mordasini 2017;Marleau et al. 2019b). In this scenario, in direct imaging surveys we would have detected only those that are born with high entropies, but the bulk of the planets that create substructures, which would be formed by a cold process, remains undetected after the PPD phase (Vigan et al. 2020). However, planets of <10 Myr may still be actively accreting the material that resides within the gap, inside the circumplanetary disk from which the planet feeds. In this case of a cold formation pathway, the resulting accretion luminosity radiated away by these planets would be easily detectable, and could be even brighter than the planet irradiation itself, especially at near-IR wavelengths (Zhu 2015;Sanchis et al. 2020). Brittain et al. (2020) circumvent this apparent inconsistency by proposing a episodical accretion of circumplanetary disks, in which accretion outbursts where planets could be detected only occur during 2 % of the runaway accretion phase.
Nonetheless, the number of systems with resolved depleted cavities or spiral patterns in scattered light which also have high-contrast total intensity observations is still low, as shown in this work. Detection probabilities are also hampered by projection effects in PPDs that are not seen face-on, especially at close semimajor axes (see the detection probability maps in Appendix C), and the emission of disk residuals in the ADI data, which are sometimes difficult to distinguish from real companions (Sissa et al. 2018;Gratton et al. 2019). Further observations of young cavity-hosting disks in both polarised and ADI modes would help constrain all these different scenarios, and to understand whether other theoretical interpretations not related to the presence of perturbers might be necessary to explain the various morphologies, such as effective dust growth at snowlines ), variations in the ionisation level of the gas (Flock et al. 2015), or photoevaporation winds (Owen et al. 2012).

Extinction effects
Throughout this study we have assumed that the observed emission of these young planets will not be affected by the presence of disk material along the line of sight. However, all forming planets have to be surrounded by circumplanetary material regardless of their mass, since they still accrete. As long as they do it, there has to be material around to accrete from, which means that planetary atmospheres can well be extincted, shifting the emission to longer wavelengths. If an extended cloud of material around the planet exists, observations might not even be detecting the planet itself, but the reprocessed emission of the dust (e.g., Stolker et al. 2020). As seen for PDS 70b, even with 1-5 µm photometry and spectroscopy, the protoplanet nature of these objects makes it difficult to disentangle the intrinsic planet emission from that of the circumplanetary environment, and whether the dust reddening the spectrum is located in the planet's atmosphere, in the circumplanetary envelope or in the circumstellar disk Haffert et al. 2019;Mesa et al. 2019;Christiaens et al. 2019;Stolker et al. 2020;Wang et al. 2020).
If the bulk of the perturbers creating the scattered-light gaps between rings have masses of a fraction of that of Jupiter as the results suggest, they might not open deep enough gaps to clear all the disk material, and the extinction could indeed completely impede their detection in the near-infrared. According to 3D high-resolution hydrodynamical simulations of disk-planet interaction by Sanchis et al. (2020), their emission (intrinsic and accretion shock luminosity) can be highly attenuated by the surrounding dust (of about ∼13 mag in H band and ∼4 in L at locations of ∼100 au for a 1 M Jup companion), depending on disk aspect ratio, surface densities, viscosities, dust processing and planet formation parameters. However, when the undepleted disk surface density is low, even low-mass objects seem to suffer little extinction at any distance in the PPD. For instance, in PDS 70 the infrared extinction has an incidence in H band only for planets < 2 M Jup at distances of < 40 au, and due to the very low surface density of the PPD around CQ Tau (Ubeira Gabellini et al. 2019), extinction is only somewhat relevant for 1 Jupitermass planets within 20 au, with an L-band extinction of 0.3 mag (Sanchis et al. 2020). For the three potential sub-Jupiter planets carving the gaps in TW Hya, van Boekel et al. (2017) estimated an extinction in H band of about 2 magnitudes during the late detached phase, and for the same system Sanchis et al. (2020) increased the mass upper limits obtained by Ruane et al. (2017) only from 1 to 2 M Jup , when extinction is taken into account. Maire et al. (2017) also found that the detection mass limits beyond the spirals in HD 135344B might be underestimated by not more than 2 M Jup . These results suggest that, depending on the individual system, the derived mass sensitivities to low-mass perturbers depleting the µm-sized dust may be too optimistic, with a critical effect in dense PPDs, and a low or moderate influence in PPDs with a lower surface density.
Perturbers with masses higher than 4 M Jup will not suffer substantially from extinction in any band or distance, according to Sanchis et al. (2020), even for very dense disks, because the material is cleared within the gap. However, those simulations do not count with a well-resolved planet vicinity, which means that the density of material is lower than it would be with a fullyresolved CPD (see ). In addition, they assume a face-on configuration of the disk, whereas the circumstellar disk could distort the vision of the planet for inclined systems (e.g., Szulágyi & Garufi 2019). For these reasons, these values can be better considered as lower limits for extinction.
In summary, we expect that the emission originating in the potential > 4 M Jup planets carving the cavities will not be heavily obscured, as it has been the case for PDS 70b (Mesa et al. 2019). Hence, in the majority of situations perturbers susceptible of being importantly extincted (>2 mag) in the H and K bands are in any case too low mass (≤ 2 M Jup ) to be detected by the SPHERE detection limits, while more massive planets that are within SPHERE's reach would probably suffer less obscuration. A more detailed luminosity estimate of the forming planets would consist in coupling evolutionary models of 'naked' planets (as in Figure 7) with individual high-resolution hydrodynamical simulations of the disk.

Potential for future detection
Based on the results of the two previous sections, giant planets in cavities seem to be within the reach of current high-contrast imagers such as SPHERE, and constitute the most promising places SPHERE/IRDIS 1.5h Fig. 11. Contrast magnitudes for 8 Myr-old planets in the NG76 population synthesis model (Emsenhuber et al. 2020a,b). Two populations of 0.05-0.5 M Jup and 3-10 M Jup planets are artificially located at 80 and 20 au, respectively. Estimated extinction coefficients have been taken from Sanchis et al. (2020), see text. The METIS contrast curve has been obtained from Carlomagno et al. (2020), and the SPHERE limits correspond to the derived contrast in this work for RX J1604, which we adopt as a representative host star. Small random shifts in separation are applied for a better visibility of the populations, and dot sizes are connected to planet radius.
where to image a giant planet. Planets carving the gaps between rings however have masses of about ∼0.1 M Jup , similar to those of the undetected planetary-mass companions responsible for the mm-dust gaps seen in high-resolution ALMA observations in the DSHARP survey (Zhang et al. 2018) and in the Taurus starforming region (Long et al. 2018). Even without considering the presence of extinction, this level of contrast is unattainable for current instruments (see Figures 7 and 8).
If the surface density depletions are indeed caused by the presence of embedded sub-Jupiter mass companions, their detection might be more favorable at longer wavelengths. A future instrument that would enable their detection is the Mid-Infrared E-ELT Imager and Spectrograph (METIS; Brandl et al. 2014). The top panel of Figure 11 shows the detection limit in L band of METIS after 1 h of integration, compared to two different populations of synthetic planets; 0.05-0.5 M Jup companions that might potentially be creating gaps in the outer disks, which we locate at 80 au, and more massive 3-10 M Jup objects opening up cavities at 20 au. Planet parameters have been obtained from the synthetic population model NG76 2 (Emsenhuber et al. 2020a,b) for a median age of our sample of 8 Myr. We considered objects with an atmosphere and in the detached evolutionary phase or later, and calculated their approximate brightnesses in L band assuming that the planets radiate as blackbodies, following van Boekel et al. (2017). The simulated contrast curve was taken from Carlomagno et al. (2020), and uses a ring-apodized vortex coronagraph (RAVC). To calculate the contrast between the host star and the synthetic population, we have assumed a host star such as RX J1604, an early K dwarf located at the median distance of our sample. We also include extinction coefficients as derived by Sanchis et al. (2020) for an unperturbed surface density of Σ = 127 g/cm 2 , consistent with the densest disks in star forming regions; we take their values for a 5 M Jup planet at 20 au (see their Table A.2) and for a 1 M Jup planet at 100 au (see their  Table A.4), and apply them to the in-cavity and in-rings planets, respectively. The sizes of the dots are correlated with the planet radius.
We see that METIS will be able to detect giant planets forming the resolved cavities in scattered light with masses > 3 M Jup , and possibly also a small fraction of those in the gaps between rings in the systems where extinction is low. On the lower panel of Figure 11 the same calculations are made for SPHERE/IRDIS in K band. As expected, it shows that the current instrumentation is in the limit of providing giant planet images within cavities, but the bulk of the sub-Jupiter-mass planets is out of reach.
This implies that the use of complementary pipelines and the weather conditions during the observations may now have a critical effect in detecting giant planets for which instruments such as SPHERE are close to the detection limis (see, e.g., Fig A.1 in Keppler et al. 2018). Obtaining the best possible observing conditions for the most promising systems in Table 4, together with detailed and varied reduction schemes, could be the most favorable case for a new protoplanet detection in a PPD in the near future.
Direct coronagraphic imaging at long wavelengths (> 3 µm) with the James Webb Telescope (JWST) will also be useful for detecting and characterizing young Saturn-like planets at projected distances > 1 (e.g., Beichman et al. 2010). For nearby and large PPDs, JWST may be able to provide the most constringent mass limits of planets creating substructures at hundreds of au, down to ∼ 0.2 M Jup around dwarf K and M host stars (Carter et al. 2020), and also to characterize their atmospheric compositions and planet parameters, which would in turn shed light on their formation pathway (Linder et al. 2019).
Finally, another way to look for these perturbers is to detect hydrogen emission lines, such as Hα, originating in the shock-heated gas during the accretion of disk material (Aoyama & Ikoma 2019). If detected, this emission constitutes the main observational signature of the ongoing formation of a substellar object (Eriksson et al. 2020). Observations of PPDs with SPHERE/ZIMPOL and MUSE have so far been unable to detect planetary-mass companions around PPDs hosts other than the pair around PDS 70, despite reaching potential sensitivities down to 1 M Jup (Haffert et al. 2019;Zurlo et al. 2020).

Conclusions
In this work we have presented a homogeneous ADI data reduction of 15 PPDs observed with SPHERE which are characterised by rings, cavities or spirals in scattered light. We provided detection limits to low-mass companions creating these substructures, and studied the difference between the current SPHERE mass limits and the estimated sensitivity that would be needed to achieve a detection in these systems. The main results of this paper can be summarised as: 1. We find that the detection limits in our PPD sample (K and earlier-type host stars of 10 Myr) are greatly affected by the assumed initial luminosity of the perturber, with median sensitivities that vary for the hottest AMES-DUSTY starts from ∼ 9 to 4 M Jup , respectively for perturbers in close (∼10 au) and wide (∼100s of au) orbits. If a cold start based on core accretion formation is considered, these detection limits are deteriorated by a factor of ∼ 3, and sensitivities to perturbers barely reach the planetary-mass regime. An assessment of the detection limits in mass for each individual system can be found in Appendix B.
2. We performed an estimation of the mass of the currently undetected perturbers via three approaches: 1) Literature simulations based on the scattered light morphology of the disk; 2) Gap width in scattered light proportional to the Hill radius of the planet; 3) Location of the scattered light cavity compared to its radius in ALMA mm continuum observations. Assuming that one perturber is responsible for creating one substructure, the underlying population of planets is located between ∼ 10-400 au and within a mass range of 0.01-1 M Jup for gaps in between rings, and more massive companions of up to ∼ 3-10 M Jup in cavities and spirals. They represent the potential objects that would be needed to create the current substructures seen in young disks, which have so far no direct correspondence to any detected population of exoplanets.
3. We compared the estimated masses to the obtained detection limits at the perturbers' positions. We find that SPHERE is about one order of magnitude away in mass from detecting the majority of planets creating the gaps in between scattered-light rings. If the presence of material attenuating the planet emission is important, this difference could be more pronounced. Perturbers creating cavities tend to be more massive, and current SPHERE detection limits seem to be good enough for potential detections, consistent with the discovery of two giant planets in PDS 70.
4. Besides PDS 70, we find that the cavities in RX J1604, RXJ1615, Sz Cha, HD 135344B and HD 34282 are the most promising systems for a future direct detection of giant planets. Perturbers of the order of a Jupiter mass might also be found in between the rings of PDS 66 and HD 97048.
5. Future imaging instruments such as ELT/METIS and the JWST may be able to explore the hottest population of planets depleting the µm-sized dust between rings. These surveys will help to constrain the initial luminosities of the perturbers and link planet masses to gap morphologies.
6. Current dynamical mass constraints assume the presence of a single planet on a coplanar and circular orbit in an isothermal disk. This premise might not be universally realistic, which could lead to overestimated mass derivations. Further work in this area exploring a larger parameter space in the properties of planets would help reduce these uncertainties on the planet masses creating substructures in the scattered light.
PA are taken from Hashimoto et al. (2012). Follow-up high resolution observations with ALMA emission found the outer cavity peaking at 75 au in the dust continuum, and an inner ring at 10 au Francis & van der Marel 2020).
Sz Cha: This disk has been observed in IRDIS H-band DPI mode. It is characterised by three different rings at semi-major axis between 0.2 and 0.6 . We adopt here the ellipse fitting parameters derived by Hagelberg et al. (priv. comm.), and classify the disk as Ring, treating the inner one as a cavity. We are not aware of current observations of this object in the mm emission with enough resolution to observe substructures.
HD 34282: Scattered light image from de Boer et al. (2020a). Two rings seems to be well-fitted, the inner one (rim R2) at semi-major axis of ∼ 89 au. The second one after scale height deprojection seems to be more coincident with a single-armed spiral. There seems to be no simulations of planets causing the spiral in scattered light yet. Inclination and PA are from this paper. ALMA data from van der Plas et al. (2017a) see a ring in the continuum extending from 75 au to 359 au, but peaking at the deprojected radius of 138 au. According to their work, a 50 M Jup object at 0.1 (∼ 30 au) could explain the ALMA ring.
We thus classify this object as Spiral, and use the inner rim to derive a planetary mass inside the scattered-light cavity via Hill radius, which will be similar to the ALMA planet candidate.
RX_J1604.3-2130A: This transition disk shows variable dips and a ring peaking at 66 au in scattered light . The dips might be originated by the presence of an inner misaligned ring. The cavity is wider in ALMA observations, which reveal an elliptical fit to the outer disk peaking at 90 au in the dust continuum (Mayama et al. 2018).
HD 36112: This transitional disk with an inner eccentric cavity is resolved in the mm wavelengths out to ∼ 50 au in deprojected distance (Dong et al. 2018), but not in PDI with SPHERE down to ∼ 15 au (Benisty et al. 2015). Scattered light observations however show a pair of spirals with a large opening angle from ∼0.25-0.45 . The object is thus classified as Spiral, and the inclination and position angle are taken from Isella et al. (2010). Baruteau et al. (2019) have simulated planet masses and locations that create these substructures. They take into account previous upper limits on planet sensitivities (e.g., Reggiani et al. 2018), and reproduce the spirals in scattered light and the sub-mm crescent substructure with two planets of mass 1.5 M Jup and 5 M Jup located at 35 au and 140 au (for a distance of 160 pc), inside and outside the spirals, respectively. We assume these simulated planets as the ones creating the substructures and compare them to our sensitivities. As no resolved inner cavity is seen with SPHERE, we do not use this object for the ALMA-Scattered light comparison.
TW Hya: SPHERE shows a face on disk with three wide albeit shallow gaps in the µm-sized dust at approximately 7, 22 and 90 au (van Boekel et al. 2017). This work derives the planet masses carving these non fully-depleted gaps from their depths following Duffell (2015). We scale these values for α = 10 −3 and obtain 0.04 M Jup for the innermost gap. We rely on Dong & Fung (2017b) for the masses of the perturbers causing the two outer gaps, 0.15 and 0.08 M Jup , respectively. We treat the three gaps as rings, as ALMA observations in the sub-mm show a ring at ∼ 2 au. The ∼22 au depression is also observed by ALMA in the surface brightness of the continuum (e.g., Huang et al. 2018), but no clear correspondence with brightness peaks can be done between the millimetre and scattered light. Even though van Boekel et al. (2017) assumed a face-on disk, here we take the inclination and position angle as derived by Qi et al. (2004) from CO imaging of the outer disk.
CQ Tau: SPHERE PDI data shows a spiral pattern extending up to ∼0.4 , but a cavity is not resolved outside the coronagraphic mask (Benisty et al. in prep.). NIRC2 observations in L band excluded the presence of giant planetary-mass companions down to 5 M Jup outside the spiral region (Uyama et al. 2020). Ubeira Gabellini et al. (2019) reported ALMA data showing a depleted cavity with a rim at ∼ 57 au. Inclination and position angle are taken from this work. We include only the gaps and substructures that were considered as to be created by perturbers in Section 4. The adopted locations of these perturbers are shown as green dashed curves. Scattered light rings are illustrated with thick rings, while arms and spiral patterns are shown as flat curled curves. The disks are projected according to their inclination and position angle in Table 2. We also include HD 100546, HD 139614 and CQ Tau for completeness, although no clear gaps or simulated planets creating the spiral pattern exist.  14. 5σ sensitivity to planet masses in RXJ1615 (see Figure B.1).
As the minimum age for this system is below the AMES-DUSTY lower limit, the minimum mass estimate for this model is not included here.

Appendix C: Detection probability maps
Here we present the detection probability maps to low-mass companions in the PPDs present in our sample (see Table 1). These maps represent the companion mass-semimajor axis parameter space of the SPHERE pupil-tracking data. To obtain them, we first converted the projected 2D limiting magnitude maps from the reduced IRDIFS data to planet masses, using the AMES-DUSTY (Chabrier et al. 2000), BEX-HOT and BEX-WARM (Marleau et al. 2019a) models. We then ran the MESS code, a Monte Carlo tool for the predictions of exoplanet search results (Bonavita et al. 2012) for each PPD system, which calculates the probability of detecting a planet of a given mass at a given separation. This is achieved by the injection of a set of test companions on different orbits, which parameters are drawn from probability distributions. The resulting projected separation of each companion is then calculated and compared to the mass sensitivity of our data at the same location, determining whether or not the planet would have been found by SPHERE. For each target we generated a uniform grid of mass and semi-major axis in the interval [0.5, 200] M Jup and [1, 1000] au with a sampling of 0.5 M Jup and 0.5 au, respectively. We then generated 10 4 orbits for each point in the grid, randomly oriented in space from uniform distributions in ω, e and M, corresponding to the argument of periastron with respect to the line of nodes, eccentricity, and mean anomaly, respectively, while the inclination and position angles were fixed to be in the disk plane (Table 2).
Article number, page 23 of 26 A&A proofs: manuscript no. aanda  Table 1 and AMES-DUSTY conditions have been assumed to convert ANDROMEDA contrast limits to companion masses. The blue curves circumscribe the probability of detecting a planet of a given mass-semimajor axis space.   Table 1 and AMES-DUSTY conditions have been assumed to convert the ASDI contrast limits to companion masses. The blue curves circumscribe the probability of detecting a planet of a given mass-semimajor axis space.