VLT/SPHERE survey for exoplanets around young, early-type stars including systems with multi-belt architectures

Dusty debris disks around pre- and main-sequence stars are potential signposts for the existence of planetesimals and exoplanets. Giant planet formation is therefore expected to play a key role in the evolution of the disk. This is indirectly confirmed by extant sub-millimeter near-infrared images of young protoplanetary and cool dusty debris disks around main sequence stars usually showing substantial spatial structures. A majority of recent discoveries of imaged giant planets have been obtained around young, early-type stars hosting a circumstellar disk. In this context, we have carried out a direct imaging program designed to maximize our chances of giant planet discovery and targeting twenty-two young, early-type stars. About half of them show indication of multi-belt architectures. Using the IRDIS dual-band imager and the IFS integral field spectrograph of SPHERE to acquire high-constrast coronagraphic differential near-infrared images, we have conducted a systematic search in the close environment of these young, dusty and early-type stars. We confirmed that companions detected around HIP 34276, HIP 101800 and HIP 117452 are stationary background sources and binary companions. The companion candidates around HIP 8832, HIP 16095 and HIP 95619 are determined as background contamination. For stars for which we infer the presence of debris belts, a theoretical minimum mass for planets required to clear the debris gaps can be calculated . The dynamical mass limit is at least $0.1 M_J$ and can exceed $1 M_J$. Direct imaging data is typically sensitive to planets down to $\sim 3.6 M_J$ at 1 $''$, and $1.7 M_J$ in the best case. These two limits tightly constrain the possible planetary systems present around each target. These systems will be probably detectable with the next generation of planet imagers.


Introduction
How giant planets form and evolve is one of the biggest challenges of modern astronomy and remains the subject of heated debates. This major goal is directly connected to the ultimate search for life over the horizon 2030 to 2040, although several astrophysical (formation, evolution, dynamics, structure and atmosphere), biological (bio-markers) and technical (new technologies developed for next generation of instrumentation) steps must be carried out in that perspective. Understanding how giant planets are formed and structured, how they evolve and interact, is critical as they will completely shape the planetary system architectures and therefore the possibility to form telluric planets capable to host life. More than two decades ago, the only planets we knew were the ones of our Solar System. With the manna of exoplanet discoveries since the 51 Peg discovery (Mayor et al. 1995), the diversities of systems found (Hot Jupiters, irradiated and evaporating planets, misaligned planets with stellar spin, planets in binaries, telluric planets in habitable zone, discovery of Mars-size planet...), the theories of planetary formation have drastically evolved to digest these observing constraints. However, we are still missing the full picture and some key fundamental questions still lack answers like: i/ the physical processes at play to pass the km-size barrier to form planetary cores, ii/ the physics of accretion to form planetary atmospheres, iii/ the formation mechanisms to explain the existence of giant planets at wide orbits, iv/ the physical properties of young Jupiters, v/ the impact of planet-planet and planet-disk interaction in the final planetary system architecture, or vi/ the influence of the stellar mass and stellar environment in the planetary formation processes. Neither, core accretion plus gas capture (CA; Pollack et al. 1996) nor disk fragmentation driven by gravitational instabilities (GI; Cameron 1978) can globally explain all current observables from planet hunting techniques. Alternative mechanisms are then proposed, such as pebbles accretion to enable core accretion to operate at wide orbits (Lambrechts & Johansen 2012), inward/outward migration or planet-planet (Crida et al. 2009;Bromley & Kenyon 2014) or simply the possibility to have several mechanisms forming giant planets (Boley 2009). In this context, each individual discovery of giant planet and young planetary system using direct imaging is rich in terms of scientific exploitation and characterization as these systems offer the possibility: i/ to directly probe the presence of planets in their birth environment, ii/ to enable the orbital, physical and spectral characterization of young massive Jupiters, iii/ to characterize the population of giant planets at all separations in synergy with complementary techniques such as astrometry (GAIA) and radial velocity adapted to filter stellar activity. Dusty debris disks around pre-and main-sequence stars are possible signposts for the existence of planetesimals and exoplanets (Matthews et al. 2014). Numerous T Tauri and Herbig stars indicate that the characteristic timescale for the dispersal of a surrounding dusty, gaseous disk is a few million years (Kennedy & Kenyon 2008). Giant planet formation is therefore expected to play a key role in the evolution of disk. This is indirectly confirmed by extant submillimeter and near-infrared images of cool dusty debris disks around main sequence stars usually showing substantial spatial structure (e.g. Eri, Vega, Fomalhaut, β Pic; see Schneider et al. 2014). It is striking to note that a majority of recent discoveries of imaged giant planets have been obtained around young, dusty, early-type stars. It includes the breakthrough discoveries of Fomalhaut b (3 M Jup at 110 AU, A4V star; Kalas et al. 2008), HR 8799 bcde (5-10 M Jup at 10-64 au, F0V star; Marois et al. 2010), β Pictoris b (8-13 M Jup at 9 au, A5V star; Lagrange et al. 2010), HD 95086 b (3-5 M Jup at 56 au, A8V star; Rameau et al. 2013) and more recently 51 Eri b (2 M Jup at 14 au, F0V star; Macintosh et al. 2015). The presence of dust and the spatial sub-structure (ring, gap, warp and other asymmetries) are possible indirect indicators of the presence of giant planets (Mouillet et al. 1997;Dipierro et al. 2015;Pinte et al. 2020). Direct imaging is here an unique and viable technique to complete our view of planetary system characteristics at wide orbits (≥ 5 au). This technique enables us to directly study the planet-disk connection to constrain the planet and disk physical properties, evolution and formation. In the case of β Pictoris, Lagrange et al. (2012) confirmed that β Pic b was actually responsible for the disk inner warp geometry, perturbing the planetesimals field and shaping the warp up to 40-60 au. HD 95086 and HR 8799 share a com- mon two-components architecture consisting of a warm inner belt (≤ 5 au) and a cold outer disk (100 − 200 au) (see Su et al. 2015). Kennedy & Wyatt (2014) actually showed that the spectral energy distributions of both systems are consistent with twotemperature components compatible with dust emission arising from two distinct radial locations. Such an architecture would be analogous to the outer Solar systems configuration of Asteroid and Kuiper belts separated by giant planets. Therefore, following the strategy of our NaCo DUSTIES (Dusty, yoUng and earlytype STar Imaging for ExoplanetS) survey (Rameau et al. 2013) that led to the discovery of HD 95086 b, we initiated a searching for giant planets with SPHERE at VLT around an newly identified sample of young, early-type stars with indication for some cases of multi-belt architecture to maximize the chances of discoveries. The sample, the observations and the data reduction and analysis are presented in sections 2, 3 and 4 respectively. The results are reported in section 5, and discussed in section 6.

Target Properties
The target selection of the survey has been obtained from the selection of large sample of young, nearby early-type stars with: declination (δ ≤ 25 o ), age (≤ 100 Myr), distance (≤ 100 pc) and R-band brightness (≤ 9.5) to favor good adaptive optics performances. Age selection criteria were applied based on different youth diagnostics (kinematics, isochrones, Lithium, H α emission, X-ray activity, stellar rotation and chromospheric activity). We also used as selection criteria the presence of significant 60 − 70 µm excess from the IRAS and Spitzer missions in the spectral energy distributions (Zuckerman et al. 1995;Zuckerman 2001;Rhee et al. 2007;Zuckerman & Song 2004b,a;Zuckerman et al. 2011;Zuckerman et al. 2013;David & Hillenbrand 2015;Moór et al. 2016) or the existence of multi-belt component analysis from Kennedy & Wyatt (2014). At the end, a total of 30 late-B, A-and early-F-type young stars, observable from the southern hemisphere, were then kept among which 22 were observed between October, 2016 and August, 2019. Their stellar properties are reported in Table 1. The age, distance, spectral type and IR excess properties are shown in Figure 1.

Observations
The SPHERE planet-finder instrument installed at the VLT (Beuzit et al. 2008) is a highly specialized instrument, dedi- cated to high-contrast imaging and spectroscopy of young giant exoplanets. It is based on the SAXO extreme adaptive optics (XAO) system (Fusco et al. 2006;Sauvage et al. 2010;Petit et al. 2014), which controls a deformable mirror with 41 × 41 actuators, and four control loops (fast visible tip-tilt, high-orders, near-infrared differential tip-tilt, and pupil stabilization). The common path optics employ several stress-polished toric mirrors (Hugot et al. 2012) to transport the beam to the coronagraphs and scientific instruments. Several types of coronagraphic devices for stellar diffraction suppression are provided, including apodized pupil Lyot coronagraphs (Soummer et al. 2005) and achromatic four-quadrant phase masks . The instrument has three science subsystems: the infrared dual-band imager and spectrograph (IRDIS, Dohlen et al. 2008), an integral field spectrograph (IFS; Claudi et al. 2008) and the Zimpol rapid-switching imaging polarimeter (ZIMPOL; Thalmann et al. 2008). The sample of young, early-type stars was observed using the IRDIFS-EXT mode, with IRDIS in the dual-band imaging (DBI, Vigan et al. 2010) mode with K 1 K 2 filters (λ K 1 = 2.1025 ± 0.1020, µmλ K 2 = 2.2550 ± 0.1090, µm), and IFS in the Y − H (0.97 − 1.66 µm) mode in pupil-tracking. This combination enables the use of angular and/or spectral differential imaging techniques to improve the contrast performances at the subarcsecond level (Racine et al. 1999;Marois et al. 2006). The choice between IRDIFS mode and IRDIFS-EXT mode is critical to optimize the detection of young, early-T or warm mid-L dwarfs planets considering the primary age and distance. Indeed, it was critical in the case of β Pic b (Lagrange et al. 2009) and HD 95086 b (Rameau et al. 2013) discoveries to properly remove quasi-static speckles that dominate performance detection at close inner angles (0.1 − 2.0 , i.e. 3 − 60 au at 30 pc), but to also maximize the emitted flux by the giant planets. For young ages (10 − 50 Myr), as the potential planets to which we are mostly sensitive are warm, dusty L-type planets with no methane absorption, the choice of the IRDIFS-EXT mode is more appropriate and was chosen for this observing campaign. For the follow-up, as candidates were only detected in the IRDIS field of view, we have opted for the use of the IRDIS the DBI mode with J 2 J 3 filters (λ J 2 = 2.1025±0.1020, µmλ J 3 = 2.2550±0.1090, µm) in pupiltracking. Thus, this second epoch provides, in addition to the possibility to check for common proper motion of the candidates relative to the primary star, the possibility to better discriminate background stars from physical young, early-T or warm mid-L dwarfs planets in the color-magnitude diagram .
The observing sequence used for the survey is the following: PSF flux reference, coronographic centering using the waffle spots, deep coronographic observation of about 70 min in total on target, new coronographic centering using the waffle spots, PSF flux reference, and sky. The PSF flux references are used to estimate the relative photometry of the companion candidates detected in the IRDIS and IFS field of view, as well as the detection limits. The coronographic centering sequence using the waffle spots sequence is critical to obtain the position of the star behind the coronograph and the relative astrometry of the companion candidates. The deep coronographic observation is obtained close to meridian to maximize the field rotation. Finally, the sky background is used to optimize the background subtraction and the flat field correction. The typical observing sequence lasts approximately 90 min including pointing and overheads. The detail of the observations per target is reported in Table 2. As a by-product of the SPHERE observation, one can access the evolution of the different atmospheric parameters seen and registered by the SPHERE XAO system (SAXO). These real-time parameters are good diagnostics of the turbulence conditions (τ 0 , r 0 , integrated wind over the line of sight) and of the XAO correction (Strehl at 1.6 µm) during the observing sequence. The summary of these SAXO parameters over the full survey is reported in Table 2, and shown in Figure 2. Given the brightness of our targets, about 70 % of the survey was obtained under median or good conditions for Paranal with typical Strehl ratio larger than 80 %. A few cases, prior to the UT3 intervention at VLT in 2017, Table 1. Description and properties of the sample. The "Exc." column indicates the presence of an IR excess. The symbol "/" means no IR excess and "Y" means with IR excess. References: (B15) Bell et al. (2015); (D15) David & Hillenbrand (2015)  were affected by the low-wind effect despite good atmospheric conditions.

Data reduction and analysis
To calibrate the IRDIS and IFS dataset on sky, the platescale and true north solution at each epoch were corrected based on the long-term analysis of the SPHERE Guaranteed Time Observation astrometric calibration described by Maire et al. (2016). The rotation correction considered to align images to the detector vertical in pupil-tracking observations is −135.99 ± 0.11 • . Anamorphism correction is obtained by stretching the image Y-direction with a factor of 1.0060 ± 0.0002. All IRDIS and IFS datasets were reduced using the SPHERE Data Reduction and Handling (DRH) automated pipeline (Pavlov et al. 2008) and additional IDL routines for the IFS data reduction (Mesa et al. 2015) at the SPHERE Data Center (Delorme et al. 2017) to correct for each data cube for bad pixels, dark current, flat field, and sky background. After combining all data cubes with an adequate calculation of the parallactic angle for each individual frame of the deep coronagraphic sequence, all frames are shifted at the position of the stellar centroid calculated from the initial star center position. For independent check, two pipelines were then used to process the data in angular differential imaging (ADI), and in combined spectral and angular differential imaging (ASDI): the IPAG-ADI pipeline , and the SpeCal (Galicher et al. 2018). These routines allow to reduce the data cubes with almost the same set of algorithms (classical ADI, Marois et al. 2006;LOCI, Lafrenière et al. 2007; PCA, Soummer et al. 2012;Andromeda, Cantalloube et al. 2015), and to exploit the spectral diversity of the IRDIS and IFS observations using ASDI techniques in addition to ADI only. Following the principles described in Galicher et al. (2018), SpeCal (and IPAG-ADI) delivers for various algorithms and observing techniques (ADI, ASDI) contrast curves, SNR maps and the possibility to locally characterize the astrometric, photometric and spectroscopic signal of any companion candidate using either a template or negative fake planet injection approach. As consistent results were found with both pipelines, the full set of observations was reduced with SpeCal (routinely used with the SPHERE GTO) using the TLOCI algorithm (in ADI and ASDI) for IRDIS, and the PCA algorithm (in ASDI) for IFS. A spatial filtering for each data cube was automatically applied to the deep coronographic observations and the reference PSFs before the use of SpeCal.
The TLOCI algorithm is implemented in SpeCal, to attenuate the background signal. The TLOCI algorithm locally subtracts the stellar speckle pattern for each frame in annuli of 1.5 × FWHM further divided in sectors. The subtraction is based on a linear combination of the best 20 (N parameter) correlated reference images calculated in the optimization region and selected to minimize the self-subtraction at maximum 20% (τ parameter), see Galicher et al. 2011 andMarois et al. 2014 for further description of the reference frame selection, and the subtraction and optimization regions. For IFS, in the PCA version we used each frame is subtracted from its average over the field of view to estimate the principal components. The spectral diversity is exploited after proper rescaling and renormalization of the IFS data cubes as detailed by (Mesa et al. 2015). Considering the significant field rotation of our observations, the first 100 principal components were subtracted.

Companion candidate detection and characterization
Using the IRDIS and IFS SNR maps provided by SpeCal, we identified by eye a total of 8 companion candidates at relatively large separation (≥ 3.0 ) in the IRDIS field of view of 6 targets (HIP 16095, HIP 95619, HIP 101800, HIP 34276, HIP 117452 and HIP 8832) of the complete survey. One companion candidate was identified at relatively close separation in the IFS field of view of HIP 41307, but latter flagged as a bright quasi-static speckle through various processing tests and therefore discarded.   Fig. 4. Left: absolute magnitude in K 1 -band versus K 1 -K 2 color for brown dwarfs with discovered companions. Right: Same but for absolute magnitude in J 3 -band versus J 3 − J 2 color. The targets from our survey are noted in red. Figure 3 shows the IRDIS imaged reduced in TLOCI ADI of HIP 16095 (bright companion located at 3.3 in the K 1 and K 2 combined filters), and IFS image reduced in PCA ASDI of HIP 41307 (quasi-static speckle discarded located at 0.5 in the combined YJH-bands) as illustration of the detection process. All companion candidates were then characterized using SpeCal with the TLOCI algorithm in ADI only and a template approach. The relative astrometry and photometry are reported in Table 3. As first diagnostics, we reported in Figure 4 Bonnefoy et al. (2018). We used here the most recent parallaxes of the young objects from Greco & Brandt (2016), and added additional companions (Gauza et al. 2015;Stone et al. 2016;De Rosa et al. 2014) at the L/T transition. At first look, we see that all detected companion candidates fall on the expected sequence of possible bound companions from the early-M spectral type for the candidates around HIP 117452, late-L spectral types for HIP 95619, HIP 16095 and HIP 8832, to early-T for HIP 34276. The companion around HIP 101800 was detected only in K 1 -band during the first epoch. After a verification of the public archive, the companion candidates around HIP 34276 (cc1 and cc2) and HIP 101800 (cc1 and cc2) were previously known and characterized as stationary background sources by Wahhaj et al. (2013) as part of the NICI Campaign around de-  Matthews et al. (2018) as physically bound, confirming that this system was actually a quadruple system with an A0 primary (HIP 117452 A), orbited by a close binary pair Ba and Bb also resolved in this survey, and additionally by a K-type star at about 75".
Follow-up observations of the candidates were automatically scheduled and obtained using the DBI mode with J 2 J 3 filters of IRDIS, well adapted to discriminate background stars from physical young, early-T or warm mid-L dwarfs planets and offering an additional epoch for a proper motion test. Follow-up observations were then processed using SpeCal with the TLOCI algorithm in ADI only and a template approach as before. All companion candidates were re-detected except the one around HIP 8832 falling outside the IRDIS field given its large separation and an observing sequence not perfectly centered with the meridian passage. The results are reported in Table 3. The use of a different pair of filters enabled us to explore the companion candidate photometric properties in the J 2 -band and J 3band based CMD this time for which we have also reported the distribution of background stars observed in previous crowded fields (see Figure 4, Right). One can directly see that most of our late-L to early-T potential companion candidates, including the previous ones identified as stationary background stars around HIP 34276 (cc1 and cc2) and HIP 101800 (cc1 and cc2), fall onto the background contaminant sequence indicating that they are most likely background stars. For a further check, we used the relative astrometry obtained at two epochs to estimate the proper motion of the companion candidates relative to their primary stars. Figure 5 shows the proper motion plots of each candidate and confirms that the companions candidates around HIP 34276, HIP 101800, and HIP 95619 are not co-moving with their primary star. The distance and proper motion of the stars, with their uncertainties, are taken from the Gaia Data Release 2 catalogue (Gaia Collaboration et al. 2018). For HIP 16095, given the relatively low proper motion of the star, the status of the companion candidate HIP 16095-cc1 remains ambiguous. However, the J 2 -band and J 3 -band based CMD still supports a background contamination. If bound, this candidate would have an estimated mass between 7 and 12 M Jup at the system age (≤ 100 Myr) and distance (88 pc) illustrative of the SPHERE detection performances around young, nearby stars beyond 10 au.
For HIP 117452 Ba and Bb, the colors and magnitudes in K 1 and K 2 compared to the predictions of the evolutionary models of Siess et al. (2000) suggest that Ba and Bb are a pair of M1 and M2 low-mass stars considering an age of 40 Myr at a distance of 42 pc. Combining our relative astrometry with the one reported by Matthews et al. (2018) and shown in Table 3, we performed a first orbit fitting of the pair. Following the method developed by Chauvin et al. (2012), we used a Markov-Chain Monte-Carlo (MCMC) Bayesian analysis technique (Ford & Gregory 2007), which is well suited for observations covering a small part of the whole orbit (for large orbital periods). We did not consider any prior information using the proximity of the primary star. The results are reported in Figure A.1 and favor a relatively inclined orbit i ∼ 98 +8 −5 deg, longitude of ascending node fairly wellconstrained at Ω = 20±2 deg, tight semi-major axis a ∼ 14 +7 −4 au, but surprisingly large eccentricities e ≥ 0.4. These large values of eccentricity are not dynamically expected given the proximity of the primary star located at a physical projected separation of ∼ 150 au, although the orbit of the binary companion around HIP 117452 is not known. Fitting solutions using a least squares Levenberg-Marquardt (LSLM) algorithm (Press et al. 1992) to search for the model with the minimal reduced chi 2 are also reported for comparison. Further dynamical study of the global system considering the debris disk architecture around HIP 117452 and the binary companion HIP 117452 BaBb configuration is be needed.

Detection limits and survey completeness
To exploit the information from the actual non-detection in IFS and IRDIS observations of the survey, the detection limits of each individual observations were then estimated. Based on SpeCal results, we derived a standard pixel-to-pixel noise map for each observing sequence corrected from the flux loss related to the ADI or ASDI processing by injecting fake planets. The detection limit maps at 5σ are then obtained using the pixel-topixel noise map divided by the flux loss and normalized by the relative calibration with the primary star (considering the different exposure times, the neutral density and the coronograph transmission). These detection limits are finally corrected from small number statistics following the prescription of Mawet et al.  (2014) to adapt our 5σ confidence level at small angles with IRDIS and IFS. The 5σ contrast curves, resulting from the azimuthal average of the detection maps, are reported for IFS and IRDIS in Figure 6.
To convert the detection limits in terms of survey completeness of the mass and semi-major axis parameter space explored with SPHERE, we used the MESS (multi-purpose exoplanet simulation system) code, a Monte Carlo tool for the statistical analysis and prediction of exoplanet search results (Bonavita et al. 2012). This code has been extensively used in previous direct imaging surveys for that same purpose (Chauvin et al. , 2015(Chauvin et al. , 2018Vigan et al. 2012Vigan et al. , 2017Rameau et al. 2013;Lannier et al. 2016). With MESS, we then generated a uniform grid of mass and semi-major axis in the interval counting the number of detected planets over the number of generated ones by simply comparing the on-sky projected position (separation and position angle) of each synthetic planet with the SPHERE 2D detection limit maps at 5σ converted in masses based on the COND (hot-start) model predictions (Baraffe et al. 2003). The primary age, distance and magnitude, reported in Table 1, are considered for the luminosity-mass conversion. The resulting detection probability map of the complete survey is then reported in Figure 7. This result shows that, despite the relatively wide spread in age (20 to 120 Myr) and distance (10 to 102 pc) of our sample, we achieved a relatively good detection probability larger than 50 % for giant planets with masses larger than 5 M Jup and semi-major axes between 10 and 500 au, sufficient for the detection of system analogues to HR 8799 or HD 95086. In principle, the degeneracy between mass and initial entropy could change these limits considerably (e.g., Marleau & Cumming 2014;Brandt et al. 2014). In practice, however, taking more realistic post-formation entropies into account strongly mitigates this problem, as shown for instance in the case of HIP 65426 b by Marleau et al. (2019).

Discussion
Our survey is composed by relatively old systems which are gas-free. Therefore, some of these systems contain debris disks. We assumed that planets are a valid explanation for the formation of debris structure, as shown in the case of the solar system where planets are known to reside between two belts of debris and in the case of HR 8799 and HD 95086 where planets are known to reside in two-temperature debris disks. The analysis of our survey follows the work by Matthews et al. (2018). The temperature values of debris belts are found in Chen et al. (2014), in which these temperatures are estimated by using a two-temperature blackbody model and a Bayesian parameter estimation to select the best model to fit the SED. The disk radii are calculated following Pawellek & Krivov (2015) assuming that dust are composed by 50% astrosilicate and 50% ice. In addition we constrained our SPHERE/IRDIS observations with dynamical arguments on the possible planetary systems hiding within the debris gaps (Shannon et al. 2016).
Mass limits are calculated with the MESS code as described in Section 6 and shown in Figure 8. The theoretical mass for a single planet to clear the observed gap is large ≥ 25M J (Nesvold & Kuchner 2014). Therefore, in our cases, we infer that the systems must be in multi-planet configuration as in HR8799, in which several planets with lower mass clear the gap. We plotted in this Figure 8, the minimum masses of planets required to clear the debris gaps, as well as their location, and their number 'Np' based on the N-body simulations of Shannon et al. (2016). This model considers only planets with low eccentricity. The mass and the number 'Np' will change if the eccentricity is larger. The mass that can be read off is the minimum mass per planet, with uncertainties based on the age of the system and on the belts radius. The minimum mass calculation assumes that planets are spaced by a typical value of ∼ 20 mutual Hill radius (R H ) which is consistent with the value of 21.7 ± 9.5R H predicted by Fang & Margot (2013).
By combining the observational upper and theoretical lower mass constrains, only a small region of parameter space is unconstrained. For 12 targets in our survey for which temperatures values are found in Chen et al. (2014), we infer a multi-planet system based on the large theoretical clearing masses. In such a multi-planet system, the widest separation planet will have a physical separation close to that of the outer debris belt, where our direct imaging limits are relatively tight. In main cases planets must be at least ∼ 0.1M J to clear the observed gap based on dynamical arguments, and in some cases the dynamical mass limit exceeds 1M J . In Figure 8, for the target HIP7345, the mass limit, ∼ 1.3M J at 90% in the gap, is close to the dynamical mass limit ∼ 0.9.
Among our 12 targets for which we have the presence of two debris belts, no exoplanetary mass companions were detected. Our sample is too small for a detailed statistical analysis. However, a non-detection in a sample of 12 stars is not inconsistent with the debris disks occurrence rate of 6.27% in a debris disk sample of planets between 5 − 20M J and 10-1000 au (Meshkat et al. 2017), since we would expect some companions might be below our detection limits. Our non-detections are also consistent with the lower occurrence rate of ∼ 1% found in Bowler (2016) and Galicher et al. (2016). The results of this 12 targets sample are not incompatible with the theory that planets are carving wide debris gaps, since in each case our direct imaging mass limits are higher than the theoretical mass limits that we calculate.
Future observations combining radial velocity, astrometry with GAIA for the inner parts (≤ 5 au), but also with the next generation of planet imagers from the ground (SCExAO, KPIC, SPHERE+, GPI2.0 on 10m-class Telescopes, then with the ELTs) and space (JWST, WFIRST) to explore will help shedding light on the existence of the planetary perturbers beyong 5 au, and potentially sculpting these architectures.

Conclusions
We have reported the observations and analysis of a survey of 22 stars with VLT/SPHERE with IRDIS in the DBI mode with K 1 K 2 filters and J 2 J 3 for the follow-up observations, and IFS in the Y − H filters, with the goal to detect and characterize giant planets on wide orbits. The selected sample favors young, i.e. ≤ 100 Myr, nearby, ≤ 100 pc, dusty, and early-type stars to maximize the range of mass and separation over which the observations are sensitive. The optimized observation strategy with the angular differential imaging in thermal bands and a dedicated data reduction using various algorithms allow us to reach a typical contrast 12.5 mag at 0.25 and 14 mag at 1.0 in IRDIS. These contrasts are converted to mass limits for each tar- Fig. 8. Constraints on the planetary systems for 12 targets in our survey. The positions of the inner and the outer debris belts are shown in orange by shading the regions inside the inner and beyond the outer. Our mass contrast limits are based on SPHERE/IRDIS and COND model predictions. Dynamical mass constraints for a slightly closer planet spacing of 20 mutual Hill radii from Shannon et al. (2016) are shown in green with masses below this value shaded. Np is the number of planets with the mass, indicated in green, required to open the gap. The uncertainties for debris belt position and dynamical mass limit are calculated based on the uncertainty in debris belt temperature, and indicated with hatching.
get. Despite the good sensitivity of our survey, we do not detect any new giant planets. We confirmed that the sources detected around HIP 34276, HIP 101800, HIP 16095 and HIP 95619 are stationary background sources by analyzing K 1 -band, K 2 -band, J 2 -band and J 3 -band images and their relative motion. The status of the candidate around HIP 8832 still requires further follow-up. HIP 117452 BaBb is resolved and confirmed as a binary companion (De Rosa et al. 2011;Matthews et al. 2018). For 12 targets of our survey, where we determined the debris belts radii, we derived upper and lower mass limits. We used Monte-Carlo simulations to estimate the sensitivity survey performance in terms of planetary mass and semi-major axis to perform the upper limit. We additionally calculate the minimum required mass for planets in the system to have cleared the observed debris gap to perform the lower mass limit. Combining our upper and lower mass limits, we are able to tightly constrain the unexplored parameter space around these systems: typically, planets must be at least ∼ 0.1M J in main cases to clear the observed gap based on dynamical arguments, and in some cases the dynamical limit exceeds 1M J . Direct imaging data from VLT/SPHERE is sensitive to planets of ∼ 3M J for a typical target in our survey. Several of the planetary systems will likely be detectable with the next generation of high contrast imagers.