© ESO 2020

1 Introduction

To improve our knowledge and comprehension of how planets form and evolve in their environments, their close neighborhoods must be investigated. In our own Solar System, planets are often surrounded by satellites and disk-like features, such as dusty rings (Alibert et al. 2005). Moreover, these elements seem to increase in number and extensions with the mass of the central planet. Although there has been a tentative detection of a Saturn-like moon around Kepler-1625 b Teachey & Kipping (2018) and a report of a disk or ring-like system candidate orbiting the pre-main-sequence star 1SWASP J140747.93-394542.6 from a complex photometric dimming event (Mamajek et al. 2012), there has been no exomoon, or exoring, confirmed yet. Another case that is still under debate is that of Fomalhaut b (Lawler et al. 2015), in which the spectral energy distribution was interpreted as the presence of a cloud of dust and planetesimals around a low-mass planet, indicating the presence of an extended source, even if more recent studies have suggested that Fomalhaut b is not a planet but the result of a collision between two large planetesimals (Gaspar & Rieke 2020).

With the current instrumentation and telescopes, the hunt for features comparable to those observed in the Solar System is considerably challenging for any detection technique. However, we would expect the mass of a disk surrounding a substellar object to scale with the mass of the companion itself (Magni & Coradini 2004); satellites and ring systems should also follow such trends (Ward & Canup 2010). Thus, the ideal targets for investigation for bound features would be massive companions spanning the regime from a few Jupiter masses to brown dwarfs. The detectability is further improved for young systems in which planets are still in the thermal emission regime and circumplanetary disks are more likely to be present. Direct imaging instruments such as GPI (Macintosh et al. 2014) and SPHERE (Beuzit et al. 2019) might be used at the extremity of their potential in the search for features around directly imaged planets and brown dwarfs.

Additional advantages arise when considering giant planets and brown dwarfs at wide separations. Indeed, the gravitational influence of a substellar object, as measured by the radius of the Hill Sphere, increases with the mass and the separation from the host star. Then, not only can we expect to find bound features (disks, rings, satellites, etc.) placed at larger distances with respect to such massive and wide companions but also the probability of interaction and pair formation of planets and brown dwarfs is enhanced. For example, the Hill radius of a Jupiter-mass planet at 1 au around a Sun-like star is roughly 0.07 au, whereas, the same planet at 50 au has an Hill radius of more than 3 au. A 50 M_Jup brown dwarf at the same separations has an Hill radius of 0.25 and 12.57 au, respectively. Since typical separations of objects detected with the direct imaging technique are in tens of au, we may find further companions inside their Hill sphere up to few au and this would allow us to disentangle the point spread functions (PSFs) of the two objects with 8-meter class telescopes equipped with xAO. Therefore, even if the terrestrial regime is almost inaccessible for this generation of instruments, we can still constrain upper limits for the mass of the satellites, look for binary planets, and investigate the presence of extended structures such as disks and rings.

Pairs of brown dwarfs in wide orbits around stars have been discovered in several studies, such as ɛ Indi B (King et al. 2010), GJ 569 B (Femenía et al. 2011), GJ 417 B (Kirkpatrick et al. 2001), HD 130948 B (Potter et al. 2002) and AB Dor Ca/Cb (via interferometry, Climent et al. 2019). This kind of configuration might extend to massive giant planets and brown dwarfs close to deuterium burning limits. Isolated brown dwarfs or brown dwarf-giant planet pairs have been already reported as well (e.g., Ophiucus 1622 2405 and 2MASS J1207334393254, Luhman et al. 2007; Chauvin et al. 2004).

On the disk side, there are a number of confirmed cases of substellar companions down to planetary masses showing significant Hα or Paschen β emission. These are considered robust evidence of on-going accretion, as in the case of PDS 70 b (Wagner et al. 2018; Christiaens et al. 2019) and c (Haffert et al. 2019), GQ Lup B (Seifahrt et al. 2007), DH Tau B (Bonnefoy et al. 2014), SR 12 C (Santamaría-Miranda et al. 2018), and USCO1610 B (Petrus et al. 2020), although in certain other cases the observed emission is most likely due to chromospheric activity (Musso Barcucci et al. 2019; Petrus et al. 2020). Some of these candidate circumplanetary disks were also confirmed in the submillimeter with ALMA, (PDS70 c, Isella et al. 2019), while others, such as those surrounding GQ Lup B (MacGregor et al. 2017) and DH Tau B (Long et al. 2019; Wolff et al. 2017; Pérez et al. 2019), escaped detection at this wavelength.

Disk or ring-like features may form at later phases, as a result of collisions that could happen either within a satellite system or between a planet or satellites and an incoming body (Wyatt & Jackson 2016). The latter case may happen more frequently for stars surrounded by debris disks that scatter large bodies inwards and outwards, increasing the chances of collisions with planets or satellites (Hyodo et al. 2017).

In this paper, we present routines based on the fake negative companion technique to suppress the contribution of directly imaged sources and investigate the residuals. These tools were systematically applied to every suitable substellar object, detected during the SHINE/SPHERE survey (Chauvin et al. 2017a) to unveil possible extended structures and to constrain their multiplicity. Since typical contrasts obtained with SPHERE are on the order of 10⁻⁵ − 10⁻⁶, for the youngest systems, the instrument can probe down to a few Jupiter masses. Thus, as mentioned before, with the current instrumentation, we are not sensitive to Solar System-like rocky satellites but only to gas giant companions and to bright and massive accretion disks. In this framework, we present the first indication of a Jupiter-like candidate that orbits a low-mass brown dwarfs bound to DH Tau.

In this paper, we present a description of the method in Sect. 2. In Sect. 3, we investigate the entire SPHERE sample and in Sect. 4, we present the DH Tau Bb candidate companion. In Sect. 5, we show the detection limits and the statistics obtained for the sample of substellarobjects observed with SPHERE. In Sect. 6, we present our conclusions. Finally, in Appendix B we present results obtained by applying our tools to synthetic data and to systems chosen ad hoc. In Appendix C, we describe the method used for the reconstruction of the PSF of the star behind the coronagraphand in Appendix D, we describe the reduction of the LMIRCam/LBTI data for DH Tau B.

2 Methods

2.1 The NEGFC technique

In order to unveil the neighborhood of known companions detected with direct imaging techniques, we need to carefully subtract the contribution of the companion itself. We applied the procedure known as the negative fake companion technique (NEGFC, Lagrange et al. 2010; Marois et al. 2010) which, in short, consists of forward modeling a detected source by injecting a negative model of the instrument PSF into the raw data and minimizing the residuals in the final image after post-processing by adjusting the PSF position and flux.

As model PSF, we used, in general, the off-axis image of the central star, which resemble the shape of the PSF of the companion of interest. The frame containing the model should have dimensions that are sufficient to include several times the full width half maximum (FWHM) of the PSF. This is needed in order to reduce artifacts in the close vicinity of the companion when subtracting the model. In SPHERE observations, for example, the PSF template for IRDIS (see Sect. 3) is an image of 64 × 64 pixels (roughly 16 times the FWHM in the H-band and 12 times in the K-band) of the central star, which is provided by the SPHERE Data Center (Delorme et al. 2017a). Moreover, all our tools rely on the angular differential imaging technique (ADI, Marois et al. 2006a) but they could also be used, in principle, with different post-processing techniques (for example, a median of the frames).

We thus implemented the NEGFC in Python, starting from the already released VIP code (Gomez Gonzalez et al. 2017) and optimizing some of its functions to look for both bound extended features and point-like sources around directly imaged companions. For the faintest objects and companions located in the innermost regions dominated by speckle contaminants, we assumed priors for the (r, θ, F) of the planet, subtracted a model of the latter in each raw frame, and calculated the χ² on a circular aperture centered on the companion in the post processed image obtained with specific techniques, such as classical ADI (cADI, Marois et al. 2006a) and the principal component analysis (PCA, Soummer et al. 2012; Amara & Quanz 2012). The χ² is defined as a sum over the N positive values divided by N − 3, if the function of merit is sum, or as the standard deviation over the values different from zero, if the function of merit is stddev. The (r, θ, F) parameters are then optimized to minimize the χ². Here, we call this procedure the single_fullcube routine.

2.2 Evolutions of NEGFC

Since the single_fullcube routine is not able to take into account the frame-by-frame variations of the flux due to changing atmospheric conditions, we propose a first evolution of NEGFC for companions visible in each raw frame (S∕N ≳ 30). The single_framebyframe routine determines (r, θ, F) for each frame of the unprocessed dataset. The single_fullcube routine is then used to suppress the contribution of the faintest or closest companion, whereas the single_framebyframe tool is applied at the brightest or farthest point sources that are detectable in each raw frame and that are clearly distinguishable from speckle contaminants.

For ground-based observations, the shape of PSFs over time is not constant. Indeed, if parameters such as seeing or wind speed change significantly during the observation, the PSF of the companion might be quite different from the first to the last frame of the cube. Also, in order to reach deeper contrasts, the coronagraphic observations are very often the best choice for exoplanets search. With this kind of setting, the star is hidden behind the coronagraph and cannot be used as a model. Thus, we need an image of the off-axis PSF of the star that might nonetheless differ significantly from the PSF of the companion since the two are not taken simultaneously. When subtracting two PSFs with different shapes induced by surrounding conditions, spurious structures may emerge in the residuals. For coronagraphic observations, then, the process of estimating the position and flux of directly imaged companions and the retrievalof features in their residuals works more efficiently for datasets taken in good and stable weather conditions.

When, instead, the observations are taken without the coronagraph, we can use the star itself as a model PSF for each of the frames in the cube. This second evolution of the NEGFC, which uses more than one model PSF, is called the multiple_framebyframe routine here. This tool helps to avoid eventual artifacts that may emerge from changes in the conditions during the observation. Indeed, the model PSF and the PSF of the companion are taken atthe same time and, thus, they have the same shape if their separation is not too large (otherwise, anisoplanatism effects may elongate the PSF and decrease the Strehl ratio). On the other hand, with these kind of settings, the processed images reach poor contrasts compared to a coronagraphic sequence so that in order to distinguish structures around the companion, or even the companion itself, we need to concentrate on brighter planets or brown dwarfs.

One further possibility to use the multiple_framebyframe routine concerns datasets where there are more companions or known background stars in the FoV of the instrument. To be useful as models, such objects should be located in regions of low speckle contaminant noise and with a similar radial separation as the one of the companion under investigation. In any case, the companion used to generate the model PSFs should be very bright with a signal to noise ratio (S/N) on the order of 50 or greater. If this is the case, the peak is expected to be well above the background and detectable in each raw frame. Therefore, suitable point sources in the FoV to be converted into model PSFs are quite rare. We show in Sect. 3.3 an example of the gain in detection limits and improvement of the residuals for the HD 1160 system. This was the only case with a point-like source in the FoV that was bright enough to apply this routine.

It is worth mentioning the four replicas of the PSF of the central star that are often used to measure the accurate position of the star behind the coronagraph. These satellite spots would be, in principle, perfect candidates of model PSFs present in each frame of the cube (Marois et al. 2006b; Langlois et al. 2013), especially for monochromatic observations for which there are no radial elongations. Unfortunately, they are placed very close to the outer working angle of the adaptive optics (AO), where speckles dominate the background. Indeed, based on the preliminary test run on real data, the results demonstrate that the subtraction of the companion using satellite spots with the multiple_framebyframe is worse than the one obtained with the single-PSF routines.

The output of each routine consists of: (1) an image of the companion obtained from a post processed image (median, cADI, or PCA+ADI) of the original cube; (2) an image of the model obtained from the post-processing of an empty cube where the model PSF is added in each frame at the position and with the flux obtained for the companion; (3) an image of the residuals around the companion, obtained by processing (either with a median, cADI, or PCA+ADI) the scientific cube from which the companion was subtracted.

An additional output of each routine is the contrast curve centered on the planet or the brown dwarf. The detection limits are obtained by re-centering the post-processed frame from which the companion was subtracted at the position of the companion itself and calculating the standard deviation in concentric annuli with a width of one FWHM each, starting from a separation of 1/2 FWHM from the center. The limits for the detection are set at the 5σ level and they include the small sample statistics following the discussion presented in Mawet et al. (2014) and the algorithm throughput effects (Jensen-Clem et al. 2018).

2.3 Tests on synthetic and real companions

As a first step, we tested our routines on synthetic companions. A more detailed discussion is presented in Appendix B.1. The most evident result is the gain in contrast close to the PSF of the injected companion once the latter is subtracted. Furthermore, we simulated a disk surrounding the synthetic companion with different inclinations and fixed contrast and separation. The inclination of the system is, indeed, a fundamental parameter for the detection of extended features, together with the luminosity and the separation from the companion (see Appendix B.2). As a result, we found that disk-like sources would also, in principle, be detectable using the method presented in this paper.

We also tested our tools on real companions choosing two ad hoc datasets. The first system we analyzed was the crowded FoV of HIP 87836, with many suitable PSFs that show signs of binarity (see Appendix B.3). The second system was, instead, chosen to prove the efficiency of the method to detect extended features. We analyzed the HST observation of the additional stellar companion that was recently discovered in orbit around GQ Lup at wide separations. GQ Lup C (Alcalá et al. 2020) showed evidence of the presence of a disk based on the infrared excess in its spectral energy distribution.We then applied the single_framebyframe routine to resolve the disk in reflected light. In each of the nine frames (one for each wavelength), we subtracted a model obtained by a nearby star from the PSF of the companion. The disk around GQ Lup C was easily detected with an estimated width of 88 mas and an inclination of 44° (Lazzoni et al. 2020).

3 SPHERE data

3.1 The SPHERE instrument

The components of SPHERE include an extreme adaptive optics system, SAXO (Fusco et al. 2006; Petit et al. 2014), using 41 × 41 actuators, pupil stabilization, and differential tip-tilt control. The instrument has several coronagraphic devices for stellar diffraction suppression, including apodized pupil Lyot coronagraphs (Carbillet et al. 2011) and achromatic four-quadrant phase masks (Boccaletti et al. 2008).

Among the science subsystems available, we used only the Infra-Red Dual-band Imager and Spectrograph (IRDIS, Dohlen et al. 2008) observations, since IRDIS has a wider field of view and most of our targets are placed at few arcseconds. The stars in our sample were observed in IRDIFS mode, with IRDIS in dual-band imaging mode (DBI; Vigan et al. 2010), using the narrow-band H2H3 filters (λ = 1.593, 1.667 μm), or in IRDIFS_EXT mode, using the K1K2 filters (λ = 2.110, 2.251 μm) of the instrument. For one of the observations of HR 2562, the broad band H filter (λ = 1.625 μm) was used.

The IRDIS data were pre-reduced at the SPHERE data center (Delorme et al. 2017a), hosted by OSUG/IPAG in Grenoble, using the SPHERE Data Reduction Handling (DRH) pipeline (Pavlov et al. 2008) and the dedicated Specal data reduction software (Galicher et al. 2018).

3.2 The SPHERE sample

The methods and routines described and tested in the previous Section and in the Appendix have as their primary objective the discovery of (sub)planetary satellites and disks around companions detected with the direct imaging technique. For this reason, we analyzed residuals around substellar companions detected with SPHERE/VLT during the SpHere INfrared survey for Exoplanets (SHINE, Desidera et al., in prep.) guaranteed time of observation (GTO) for 27 young systems hosting one or more companions, with masses spanning from Jupiter to brown dwarfs range. Among the sample we discarded TWA 5 B and HD 284149 B because the central stars are close visual binaries (Macintosh et al. 2001; Bonavita et al. 2017), thus the model PSF obtained from the latter is not suitable for our purposes, along with ROXs42 B and GSC 6214-210 B due to their poor or short datasets.

In Table A.1, we provide some brief information on the observations of the 23 systems in the sample and in Table 1, we list the characteristics of the 27 objects in the analysis. Several of the datasets have already been published in papers by the SHINE team (references in Table A.1) and others will be presented in detail in forthcoming papers that are currently in preparation. The ages were retrieved either from the literature or from those adopted in SPHERE’s statistical paper (Langlois et al., in prep; Desidera et al., in prep.; for a preliminary report see Meyer 2019) for the objects observed in the first half of the survey.

3.3 Results

Among the objects analyzed in the sample, we distinguish two main categories regarding the residuals we obtained: 15 objects (51 Eri b, β Pic b, GQ Lup B, HD 1160 B, HD 4747 B, HD 95086 b, HIP 65426 b, HIP 107412 B, HR 2562 B, HR 3549 B, PDS 70 b, PZ Tel B, HR 8799 c, d, and e) are located in the innermost part of the system (within 1″), in regions where speckle contaminants are very bright and mostly dominate the residuals once the companion is subtracted. We show in Fig. 1 an example of this kind of system, HIP107412 B, with the brown dwarf placed at 0.27″. The residuals are clearly dominated by bright speckles that cover any fainter feature around the companion.

The other 12 companions are instead orbiting far from the central star and the speckle noise is considerably lower. We make further distinctions between three faint companions (GJ504 B, HD 19467 B, and HR 8799 b), for which the subtraction proves effective (see Fig. 2), and nine bright companions (DH Tau B, CT Cha B, Eta Tel B, HIP 78530 B, HIP 64892 B, AB Pic B, TYC 8047-232-1 B, TYC 7084-794-1 B, and HD 1160 C), for which spurious features may appear in the residuals. Indeed, it may happen that speckles generated around low contrast companion are bright enough to be above the background.Since we are performing post processing techniques with the cube centered on the star, speckles around off-axis PSFs are not cancelled and may contaminate the residuals. This kind of issue does not emerge closer to the star becausethe stellar speckle contaminants dominate there. A second and more frequent bias emerges due to the variation of the PSF during the observation or elongation of the PSF placed close to the edge of the FoV (astigmatism). Thus, we performed further checks, such as testing different model PSF and different apertures.

As an example of the analysis performed to further investigate the nature of the residuals, we show in Fig. 3 the results for HD 1160 C and HD 1160 B. In the residuals of HD 1160 C (last panel of the first row of Fig. 3), a number of structures are present, such as the bright spot in the north of the image (light-blue circle) and the halo close to the center (red ellipse). In order to obtain the second row of Fig. 3, we used the multiple_framebyframe routine on HD 1160 B using as multiple model PSF HD 1160 C. The halo is still present in the residuals, as a negative excess in this case, whereas the bright spot is canceled. The latter is indeed associated to an effect that is, in turn, due to low wind conditions during the observations (Milli et al. 2018), while the presence of the disk-like feature is due to an elongation of the PSF of the C companion, which is as expected given its large separation from the host star. In Fig. 4, we also show the contrast curves (see Sect. 2.2) obtained for HD 1160 B when applying the multiple_framebyframe routine using the other companion in the FoV with respect to the results obtained with the single_framebyframe routine and single_fullcube with the star used as model PSF. The contrasts achieved here are calculated both with respect to the star (left y-axis) and to HD 1160 B (right y-axis).

A first gain is obtained when considering the single_ framebyframe with respect to the single_fullcube due to variable weather conditions during the observation.The further gain in contrast with the multiple_framebyframe routine is evident between 0.05″ and 0.1″, where the bright spot in Fig. 3 is suppressed, whereas the subtraction is slightly worse closer to the central peak due to the difference in shapes between the two PSFs.

We show in Fig. 5 another example of extended structure obtained around, TYC 7084-794-1 B. These results are particularly misleading because they appear in two different observations, and both in H and K, as a disk-like shape as shown in Fig. 5. However, the extended features are produced by a rotational smearing. This happens any time the exposure time for singular frames is too long with respect to the rotation of the objects in the FoV. Usually, this effect is very small and negligible, with the exception of a few systems that have to be treated individually. TYC 7084-794-1 is an example of these since it is very close to the zenith at meridian passage from Paranal and with a detector integration time (DITs) of only 64 s, the brown dwarf, placed at a separation of 2.99″, changes its position by roughly five pixels. For this reason, the PSF of the object is elongated with respect to the PSF of the star. One possible way to overcome this kind of effect might be to elongate the model PSF in the same way as the PSF of the object or to use a reference PSF in the FoV at similar separation, if any. We show in Fig. 6 the sensitively improved subtraction obtained after considering a PSF model elongated in the same direction of the BD.

Based on the analysis performed on a few test companions, it emerged that structures in the residuals must be treated very carefully in order to discriminate between those that are artifacts and those that are real. Starting from the end of 2017, one of the greatest source of artifacts in the residuals, the low wind effect, is almost completely fixed, with a number of affected nights at the level of ~ 3% (Milli et al. 2018).

Useful information about the observational parameters are contained in the SPARTA files that are produced together with the scientific data by the SPHERE Data Center. These files provide information on how the wind velocity and the seeing varied during the observation for each frame. Moreover, non coronagraphics images of the star collected with the differential tip-tilt sensor (DTTS, see Appendix C) in the H band can be used to test the quality of the PSF during the time sequence (see Appendix C for a detailed description). Thanks to this supplementary information, we can determine, with higher confidence, the nature of the residuals that we obtain.

Fig. 1 cADI image of HIP107412, the model, and the residuals.

Fig. 2 cADI image of HR8799 b, the model, and the residuals.

Fig. 3 Upper panel: cADI image of HD 1160 C, the model and the residuals; bottom panel: cADI image of HD 1160 B, using HD 1160 C as model and the processed image of the subtraction of the two.

Fig. 4 Contrast curve around HD 1160 B using one single model PSF with the single_fullcube (blue curve) and single_framebyframe (orange curve) routines, and multiple model PSFs using HD1160 C in place of the star for the multiple_framebyframe routine (green curve). Contrasts are calculated with respect to the central star on the left y-axis and with respect to HD 1160 B on the right one.

Fig. 5 cADI image of TYC 7084-794-1 B, the model and the residuals for the first (upper panel, H band) and the second (bottom panel, K band) epochs.

Fig. 6 cADI image of TYC 7084-794-1 B, with the model taking into account the smearing of the PSF of the companion and the residuals for the first (upper panel, H band) and the second (bottom panel, K band) epochs.

4 DH Tau B

Among the targets we analyzed, the most convincing case is represented by a candidate companion close to DH Tau B.

4.1 The DH Tau system

DH Tau is a young (age ~1.4 Myr, see Table 1) low-mass star in the Taurus star-forming region. DH Tau is still in the accreting phase and forms a wide binary (15″ projected separation) with a similar star, DI Tau. A low-mass brown dwarf or massive planet (DH Tau B) was discovered by Itoh et al. (2005). Its mass is estimated to be close to the deuterium burning limit. DH Tau B was shown to be accreting from prominent Hα and Paschenβ emission (Zhou et al. 2014; Bonnefoy et al. 2014). ALMA observations by Long et al. (2019) resolved a faint and compact disk around the central star, while no disk around DH Tau B was detected (Wu et al. 2017; Pérez et al. 2019).

4.2 Observations

The DH Tau system was observed across multiple epochs with SPHERE, on the dates 2015-10-26, 2015-12-19, and 2018-12-16 within the guaranteed time of the SPHERE Consortium in IRDIFS mode (H2H3 filters) and on the date 2020-01-05 after the acceptance of an open time proposal (0104.C-0327(A)) in IRDIFS_EXT mode (K1K2 filters). The SPHERE data were reduced as described in Sect. 3. We also analyzed the SPHERE differential tip-tilt sensor (DTTS) images, which show a reconstruction of the PSF of the star behind the coronagraph (see Appendix C). With this supplementary information, we can evaluate whether the PSF was affected by consistent variations in atmospheric conditions and, thus, if potential structures detected in the residuals are more likely to be real features or spurious artifacts.

The system was also observed with LBTI using the LMIRCam instrument in L′ band (λ = 3.6 μm) on the date 2019-09-19. More details about LBTI/LMIRCam and the data reduction performed are given in Appendix D.

Fig. 7 From the top: epoch 1 (2015-10-26), epoch 2 (2015-12-19), epoch 3 (2018-12-16) in H2-band, and epoch 4 (2020-01-05)in K1-band for DH Tau B as observed with SPHERE/VLT, and the observation (2019-09-19) in the L′ -band obtained with LMIRCam/LBTI. The candidate companion is visible in the right panels of the SPHERE observations. There is no detection for the LBT data.

4.3 The candidate companion

The detections were initially retrieved from the first two SPHERE datasets listed in Table A.1. The results are shown in the top two rows of Fig. 7.

To confirm the detection, DH Tau was re-observed with SPHERE with the H2H3 and K1K2 filter pairs. With regard to the H-band observation,the candidate was detected in the residual at a similar separation and position angle with respect to the previous observations (third row of Fig. 7). In the K-band, instead, the candidate is resolved only in K1 (fourth raw of Fig. 7), whereas in K2, a marginal excess is detected in the residuals even if it is not significantly above the background.

To characterize the candidate, which we call DH Tau Bb, we discarded the first epoch due to poor observational conditions and a strong low-wind effect. In Fig. C.1, we show, as an example, the DTTS images of the PSF of the central star at different times during the first and second observations taken with SPHERE (see Appendix C). In the first epoch, thereare constant variations to the shape of the PSF almost throughout the entire observation. During the other observation (and, similarly, during the third and fourth datasets), there are no significant differences between the PSFs, thus strengthening the detection of the candidate companion.

The error bars on the position and contrast of DH Tau Bb were evaluated empirically removing a model of the secondary companion from the raw data, and injecting fake secondary companions at similar contrast and angular separation from the brown dwarf. We then retrieved the contrast and position of the injected fake companions with the tools described in Sect. 2, and evaluated the mean error on the retrieval.

As mentioned above, a disk was detected around the central star with ALMA observations (Long et al. 2019). Since the disk was very faint in the ALMA observations and was not detected in scattered light using our own SPHERE observations, we do not expect it to impact the PSF subtraction. The parameters retrieved for DH Tau Bb are listed in Table 2.

We did notice a drift in the separation of the PSF of the candidate from the H2 to the K1-band. This might indicate that the candidate is not real since it changes its position roughly proportionally to λ∕D. On the other hand, the change in position might be explained by a combination of other factors. First of all, the PSF in the K-band is much broader than in the H-band so that it is more difficult to disentangle the two objects. Moreover, we notice that the position of the candidate is coincident with the position of the first Airy ring of the central companion. Thus, the superimposition of the flux of DH Tau Bb with the flux coming from the Airy ring of the brown dwarf, which is not well-subtracted, might cause a shift in the peak of the candidate PSF. Additionally, the background at increasing wavelengths is more noisy so that the signal to noise ratio is worse in the K-band than in the H-band.

The position angle of the candidate is even more difficult to determine with high precision since, at the separation of the candidate, the uncertainty of one pixel corresponds to ~ 10°. We obtained different values for each epoch and wavelength, with a maximum variation of ~ 15° between the second and third epochs. This might be partially due to the orbital motion of DH Tau Bb around the brown dwarf. In fact, assuming a face-on system and a circular, coplanar orbit for the candidate around the brown dwarf, we would obtain a period of ~ 320 yr for a radial separation of ~0.08″. In the four-year baseline, the candidate companion should then change its position angle by ~ 4°, or more if the orbit is not circular and the DH Tau Bb is close to periastron passage. Thus, taking into account the orbital motion, the unknown inclination of the system and the error bars for each position angle, the values obtained are reasonably consistent. Moreover, the time baseline between the epochs allows us to rule out the (a priori extremely unlikely) possibility of a background object projected by chance very close to DH Tau B, considering the proper motion of the central star (μ_α : 7.07 ± 0.11; μ_δ : − 20.70 ± 0.08 from Gaia DR2). This would imply a shift of about 21 mas toward west and 61 mas toward North of a stationary background source with respect to DH Tau B between epochs 2 and 3.

The last parameter retrieved was the contrast. Adopting the age of 1.4 Myrs, we then calculated the mass of the candidate using the BT-SETTL models. The values obtained for the mass of DH Tau Bb, ~1 M_Jup, are consistent at each epoch and in each filter at which the candidate was detected (H2, H3 and K1). The colors H2–H3 and H2–K1 of the candidate point towards an early T-type, as shown in Fig. 8, where we compare them with those for a collection of other low-mass objects (Bonnefoy et al. 2018).

In order to confirm the detection and better constrain the spectral characteristics of DH Tau Bb, we observed the system in L’-band with the LMIRCam/LBTI instrument (see Appendix D). The candidate was not detected in the residuals, as shown in the last panel of Fig. 7. Therefore, we determined only an upper limit to the contrast of the object in L′ ~ 9.2, which corresponds to an upper limit in mass of 2.6 M_Jup. This is consistent with the mass of DH Tau Bb as obtained in H and K band with SPHERE. The non-detection in L′-band strengthens the hypothesis that the companion, if real, is an early T-type object rather than a late L spectral-type that would have been more common considering the young age of the system. Indeed, from the color-magnitude diagram, an L-type companion is expected to be brighter in the L′ -band, with a gain in Δmag ~ 1−2 with respect to the H-band. The corresponding contrast would have been ~3 − 6 × 10⁻⁴ that is (marginally) detectable in the observations taken with LMIRCAM/LBTI.

The DH Tau system has been observed with many different instruments since its discovery (Itoh et al. 2005). Therefore, we analyzed archive data taken with the WFC3/Hubble Space Telescope (HST) and with the NIRC2/Keck Telescope. For both observations, we did not detect any signs of binarity (the results are omitted from this paper). Indeed, we did not expect to find prominent features either in reflected light or in the poor-quality L′ dataset available.

Fig. 8 Color-magnitude diagrams for the IRDIS/SPHERE photometry for the candidate DH Tau Bb compared to a collection of objects from the literature (Bonnefoy et al. 2018).

4.4 Discussion

The observations described in this paper suggest the presence of a planetary-mass companion (1 M_Jup) at 10 au around DH Tau B. We consider the object as candidate rather than a confirmed companion because of the small wavelength-dependence of the separation and the instability of the PSF affecting some of the observations. On the other hand, the roughly consistent position and flux of the candidate in four different datasets spanning several years argues against an instrumental artifact and it is unique among the targets considered in our study.

This would be the first case of a planetary-mass companion around a brown dwarf bound to a star. If the candidate is confirmed, the DH Tau system will become the first laboratory to test formation mechanisms for binary planet and brown dwarf systems. While a detailed study of the possible formation scenarios is not warranted at this stage, we briefly discuss some possibilities below.

We may confidently exclude the formation of DH Tau Bb by core accretion in a disk surrounding the brown dwarf due to the large masses involved and the wide separation. In order to generate a secondary gas giant, the circumplanetary disk would have to be unrealistically massive with respect to the mass of DH Tau B.

Some disk instability models showed outcomes similar to those observed here (Stamatellos & Whitworth 2009) but with very low probability. However, more recent models of disk instability (Forgan et al. 2018) fail to reproduce the architecture observed for the DH Tau system. Considering the low mass of the central star (0.37M_⊙), this mechanism does not appear as a plausible explanation. Another mechanism that could explain the existence of this system configuration is presented in Ochiai et al. (2014), where the authors assess that the formation of planet-planet pairs should be a quite efficient mechanism, with a 10% success rate for planets that undergo orbital crossing. Therefore, it might be possible that DH Tau B and DH Tau Bb formed independently in the circumstellar disk of DH Tau and then go through orbital crossing and close encounters, forming a binary system. However, the dynamical simulations by Ochiai et al. (2014) predict that pairs of planets formed in this way shrink to a very close separation through tidal dissipation. The wide separation observed might be explained in this framework if we manage to catch the high eccentric phase soon after the encounter.

The latter scenario considers the system to have formed as a low-mass triple system through cloud fragmentation. The mass ratio of 0.027 between DH Tau B and DH Tau and their wide separation is consistent with this mechanism. On the other hand, companions to brown dwarfs are typically found to have mass ratios close to unity, both for isolated brown dwarfs in the field (Fontanive et al. 2018) and for the few cases of bound brown dwarfs (see Sect. 1). The 0.1 mass ratio between DH Tau Bb and DH Tau B would, then, appear quite unusual. A 1 M_Jup object is also expected to be small for cloud fragmentation, being close or below the minimum mass for this process that is quoted in the literature (Boyd & Whitworth 2005; Whitworth & Stamatellos 2006; Zuckerman & Song 2009).

The orbital crossing and cloud fragmentation scenarios, even with the mentioned caveats, seem to be the most likely to explain this peculiar system. Once the candidate is definitively confirmed, dedicated simulations will be needed to investigate these two formation channels in more depth.

An alternative interpretation is that the feature detected in the residuals is a portion of a disk surrounding the brown dwarf. Indeed, accretion through a disc would be consistent with the significant Paschen β emission for DH Tau B. This could explain the shift with the wavelength since our analysis is sensitive to grains of different sizes. Moreover, an extended source would imply large uncertainties in the position angle. However, the disk was not detected in ALMA observations and it was estimated to be very compact around the brown dwarf (Wu et al. 2017; Pérez et al. 2019). This would be inconsistent with our estimate of a separation of ~ 10 au.

5 Statistics of binarity of bound companions

5.1 Detection limits

For each companion presented in Table 1, we took the flux that was derived from the observations listed in Table A.1 and converted it into mass accordingly with regard to the age adopted for the system (for systems with multiple observations available, we used the ones marked with an asterix in Table A.1). We used the BT SETTL models (Allard et al. 2003) in most cases with the only exceptions of 51 Eri b, HD 95086 b,and GJ 504 b for which the ages and contrasts fell out of the boundaries and the AMES COND (Baraffe et al. 2003) model was applied instead. We then calculated the detection limits around 25 out of 27 companions(we excluded HD 1160 C because it falls in the stellar mass range and PDS 70 b because its Hill radius falls inside one FWHM) in order to investigate the contrast reached in their close neighborhoods. The detection limits are obtained from a cADI image of the cube where each companion was already subtracted, as described in Sect. 2.2. The contrasts are then converted into masses using the mean ages of the systems listed in Table 1 and the BT-SETTL models for the most of the sample, with the exception of 51 Eri b, AB Pic B, η Tel B, GJ 504b, HD 95086 b, HIP 65426 b, TYC 7084-794-1 B, and TYC 8047-232-1 B for which AMES-COND models were used instead since contrasts fell out of the range adopted by BT-SETTL conversions. The Hill radii are then calculated adopting the mean value of the masses in Table 1 and the detection limits are interrupted at those separations.

We show in Fig. 9 the detection limits for the close-in (upper panels) and wide (bottom panels) companions. Since HD 1160 B is the only close-in companion with a large Hill radius compared to the other planets and brown dwarfs placed within 1″, we moved it from the first to the second group so that shorter contrast curves are better visualized. From these curves, it emerges how our detection possibilities are limited in the inner part of the system due to the contamination of the star and smaller Hill radii whereas, far from the center, the region of influence around directly imaged companions is background limited and much wider.

5.2 Statistical results

We constrain the binary frequency of bound companions using the statistical code described in Fontanive et al. (2018). This Markov chain Monte Carlo (MCMC) sampling tool was built using the emcee (Foreman-Mackey et al. 2013) python algorithm, enabling the derivation of statistical constraints on stellar or substellar populations (occurrence rates, distributions of companions in the parameter space) based on observed surveys (e.g., Fontanive et al. 2018, 2019).

Binary companions bound to stars have not been studied in depth to date, neither observationally nor theoretically. As a result, the shapes of the companion distributions in mass and separation for such binary systems currently lack any empirical constraints or theoretical predictions. We thus need to rely on assumptions for the underlying population distributions in order to constrain the binary frequency of these bound companions. Given the relatively low masses of our probed planets and brown dwarfs in comparison to isolated, old brown dwarfs in the field, we adopt the mass ratio and separation distributions derived in Fontanive et al. (2018) for late-type T5−Y0 objects in the Solar neighbourhood. This includes a lognormal distribution in separation with a peak at 2.9 au and a logarithmic width of 0.21, and a power law distribution with index 6.1 in mass ratio (see Fontanive et al. 2018 for details).

From the calculated detection limits around each target and assuming that those same distributions hold for our surveyed systems, we constrain the binary rate of bound substellar companions between separations of 0.5−100 au and mass ratios of q = 0.05−1. Figure 10 shows the obtained posterior distributions from the MCMC analyses performed for respective cases of 0 and 1 detection. As expected, the scenario with no detected companion only provides an upper limit on the binary frequency, of f < 8.1% (< 21.2%) at the 1-σ (2-σ) level. In the case where our reported candidate DH Tau Bb is real, we derive a binary fraction of $f={7.3}_{-5.6}^{+8.9}$% on those mass ratio and separation ranges, which is in good agreement with the results for field brown dwarfs from Fontanive et al. (2018) on a comparable region of the parameter space.

We emphasize that the obtained frequencies assume a fixed distribution of companions in the separation-mass ratio space and that these results only hold under those assumptions. While it would be very interesting to explore other possible underlying shapes of companion distributions that may differ from the field brown dwarfs, the existing data do not currently allow for such an analysis. It is worth noting, however, that the one detection present in our observed sample (if confirmed) has a fairly wide separation (~10 au) and a very low mass ratio (q < 0.1) compared to the peaks of the adopted field distributions (~3 au, q ~ 1). Late-type field binaries indeed rarely have separations larger than around 10 au, and typically show strong preferences for near-equal-mass systems. The disparity between the model and observed populations is further highlighted by the fact that low q values like that of the DH Tau B system are only reached around few targets in the SPHERE sample. This suggests that the true binary companion distributions for bound planets and brown dwarfs are likely to be somewhat different than for the field population, although deeper observations of larger sample sizes will be required to confirm this and to further explore formation theories for bound and isolated systems.

Fig. 9 Detection limits for companions placed within 1″ separation from the host star (upper panels) and beyond 1″ (lower panels, HD 1160 B was moved to the second group for a better visualization of the results). Figures on the left: contrasts with respect to the central star as a function of separation (in arcsec); right: detection limits expressed in mass as a function of separation (in au).

Fig. 10 Posterior distributions obtained for the frequency of bound binary companions from the MCMC analysis, in cases of 0 (blue) and 1 (red) detection in our survey. Frequencies were derived over separation ranges of 0.5−100 au and mass ratios of q = 0.05−1. The red dashed line shows the peak of the posterior density in the 1-detection scenario, at 7.3%, while only an upper limit was obtained in the 0-detection case. The shaded areas correspond to the 68% confidence intervals.

6 Conclusions

In this paper, we present tools based on the NEGFC technique to suppress the flux of known faint planets and brown dwarfs in order to detect potential companions or disks around them. The different routines are optimized in order to determine, with the best possible precision, the position and the flux of the companion, to subtract it, and to look into the residuals for additional faint features bound to it. We distinguished among two kinds of routine: the single_fullcube and the single_ framebyframe routines, which use a single model PSF, and the multiple_framebyframe routine, which uses more PSFs as models to subtract the companion. The single-model PSF tools are usually used for coronagraphic observations, whereas if a suitable PSF is visible in each raw frame the multiple-PSFs routine can be applied (both in coronagraphic and non-coronagraphic datasets). The latter tool is very helpful to suppress spurious features that may appear in the residuals due to the variation of the PSF during the observation. However, PSFs that are bright enough and at a good-enough separation to be used as models are rarely found in the FoV of substellar companions. In non-coronagraphic observations, instead, the central star can be used as a model for each frame but such settings are not common for planets hunting due to the poor contrast that can be achieved.

These tools were first tested on synthetic point sources and disks and on selected on-sky test cubes. We then applied it to the substellar companions imaged by SPHERE within the SHINE survey. This is the first methodical study that aims to look for circumplanetary disks and satellites or, more specifically, for the kind of instrument that we have, pairs of bound planets. Among the objects in the sample, we identified a candidate companion around DH Tau B, with a projected separation of ~ 10 au from the brown dwarf and with mass a of ~ 1 M_Jup. The DH Tau system was observed with SPHERE (H2–H3 and K1–K2 bands) in four different epochs and with LMIRCam/LBTI (L′ band). The candidate was detected in each of the H-band datasets and in the K1-band. The K2 and L′-band observations did not lead to a detection but they were used to obtain only the upper limits on the mass of the candidate, which is consistent with the previous results. The color-magnitude diagrams point towards an early T-type object.

We derived the contrast curves for each of the 25 suitable substellar companions observed with SPHERE and we compared the statistics of the multiplicity of our bound companions to the one published in Fontanive et al. (2018) for isolated objects. If the companion around DH Tau B is real, the frequency of binary substellar companion is ~ 7.3%, which is consistent with the result obtained for field brown dwarfs.

The techniques presented in this paper will serve as a useful tool for the hunt of companions and extended sources around directly imaged companions in the future. Indeed, in the long-term, we plan to adapt the tools described here to prepare for the next-generation telescopes, such as JWST and ELT.

MIRI@JWST will provide high-contrast imaging and IFS imaging from 5.6 to 28.8 μm. The greatest advantage of this telescope is the stability of the PSF over time. For MIRI, the expected contrast goes from 10⁻⁴ to 10⁻⁶ at separations > 2 arcsec (Danielski et al. 2018). Even if the contrast were lower than what is achievable with SPHERE, it will be more efficient for colder and smaller objects due to the longer wavelength coverage.

A dedicated high-contrast imager on an ELT (e.g., Kasper et al. 2011) is expected to reach contrasts of 10⁻⁸ at 30 mas to 10⁻⁹ at larger separations. These contrasts will allow for the detection of Earth-like planets. With this instrument, which will also provide a very stable PSF over time, exomoons and exorings will be much more efficiently detected. These detections will increase the probability of finding habitable zones even at far distances from the host star, much the same way as the moons of our Solar System are best candidates for the research of life forms.

Acknowledgements

SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF-Osservatorio di Padova (Italy), Observatoire de Geneve (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France) and ASTRON (Netherlands) in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland) and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004-2008), grant number 226604 for FP7 (2009-2012) and grant number 312430 for FP7 (2013-2016). This work has made use of the SPHERE Data Centre, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CESAM (Marseille), OCA/Lagrange (Nice) and Observatoire de Paris/LESIA (Paris). We thank P. Delorme and E. Lagadec (SPHERE Data Centre) for their efficient help during the data reduction process. This work is based on observations collected at LBT Observatory. The LBT is an international collaboration among institutions in the United States, Italy and Germany. LBT Corporation partners are: The University of Arizona on behalf of the Arizona university system; Istituto Nazionale di Astrofisica, Italy; LBT Beteiligungsgesellschaft, Germany, representing the Max-Planck Society, the Astrophysical Institute Potsdam, and Heidelberg University; The Ohio State University, and The Research Corporation, on behalf of The University of Notre Dame, University of Minnesota and University of Virginia. We acknowledge financial support from the “Progetti Premiali” funding scheme of the Italian Ministry of Education, University, and Research. A.Z. acknowledges support from the FONDECYT Iniciación en investigación project number 11190837. We acknowledge financial support from the Programme National de Planetologie (PNP) and the Programme National de Physique Stellaire (PNPS) of CNRS-INSU. This work has also been supported by a grant from the French Labex OSUG020 (Investissements d’avenir – ANR10 LABX56). The project is supported by CNRS, by the Agence Nationale de la Recherche. This work has been supported by the project PRIN-INAF 2016 The Cradle of Life – GENESIS-SKA (General Conditions in Early Planetary Systems for the rise of life with SKA). We also acknowledge the use of the Large Binocular Telescope Interferometer (LBTI) and the support from the LBTI team. T.H. acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. We also want to thank the referee for the useful comments.

Appendix A Additional table

Appendix B Data analysis

B.1 Synthetic companions

In order to understand how much we can gain from the subtraction of the companion in term of detection limits and which kind of point sources we may expect to detect, we simulated a companion (C1) in the FoV of TYC 8047-232-1, observed with SPHERE in date 2015-09-25 (see Table A.1). This C1 was injected at a separation of r = 3.9″ and with contrast of 10⁻³ with respectto the central star.

To obtain more realistic residuals, we used two different PSFs, one for the injection of the fake companion and one as a model. At this stage, the differences between the single_ framebyframe and single_ fullcube routine are not appreciable since the injected PSF does not vary. The detection limits around C1 are calculated on a cADI image obtained before and after the subtraction of the companion using single_framebyframe routines and following the procedure described in Sect. 2.

Residuals and contrast curves centered on C1 are shown in Fig. B.1. The gain in contrast when we get very close to the PSF of the companion appears once the latter is subtracted. This region is also the most interesting to investigate when looking for bound features since, in most cases, it is inside the Hill radii of substellar companions.

In order to simulate two bound substellar objects, we injected very close to C1 a second fake companion (C2) with increasing separations and with contrasts corresponding to the values at detection limits as shown on the contrast curves in Fig. B.1. The residuals for three different positions are shown in Fig. B.2. In the upper row, the second companion is placed at ~ 0.05″ from the main peak with a contrast of 5 × 10⁻⁵, in the mid row at ~0.075″ and with contrast of 4 × 10⁻⁵ and in the bottom panel at ~0.1″ with contrast of 2 × 10⁻⁵. From these preliminary tests, it emerges that our procedure is, in principle, quite efficient in unveiling additional companions placed very close to the PSF of directly imaged planets and brown dwarfs. However, artifacts that are present in the residuals, especially at few λ∕D, may hide or simulate the presence of faint objects. In the upper raw, for example, the fake companion is detected in the residuals even if it could be interpreted as an additional artifact attributed to the PSF subtraction. For such a scenario, more datasets would be needed. Once, instead we move farther from the central peak, the PSF of the secondary is more evidently disentangled by the surrounding residuals.

B.2 Synthetic disks

Since the method described in this paper is optimized to subtract the best point-like model to a given PSF, any conspicuous deviation from such a profile should be visible in the residuals whenever it is above the background noise. For this reason, we should be able not only to detect point-like sources but also extended ones. In order to understand which kind of residuals we should expect from extended sources, we inject a circular annulus around the mock planet C1 (see previous section) with different inclinations.

Therefore, each raw frame was initially de-rotated into its real position, a one pixel-width annulus with uniform contrast was inclined andinjected around the companion. For the forward-modelling process, the planet+annulus were then convolved with the PSF of the instrument. Finally, the frame was rotated again in the original configuration. We must stress that such modelsare only qualitative and not representative of the physical processes involved (see discussion below). However, more accurate simulations of disks around substellar companions are beyond the purposes of this paper.

Fig. B.1 Upper panel: cADI images of, from the left, the companion, the model PSF, and residuals. Bottom panel: contrast curves centered on the position of the companion as obtained in the post-processed frame with (blue curve) and without (orange curve) the companion. Contrasts are calculated with respect to the star (left y-axis) and with respect to the companion (right y-axis).

Fig. B.2 Injection of a second fake companion at the detection limits placed at 0.05″ (upper panel), 0.075″ (mid panel), and 0.1″ (bottom panel) from the peak of the primary.

We then applied the single_framebyframe routine to the cube with C1 and the synthetic disk. We show in Fig. B.3 the results obtained. The circular annulus was injected at a separation of ~ 35 mas, with an integrated flux ratio with respect to the companion of F_disk∕F_pl = 1∕10. Moreover, we simulated an inclination of the disk moving from a face-on system (upper row), through 45° of inclination (mid row), to an nearly edge-on configuration (85°, bottom row). From Fig. B.3 it emerges how much the inclination influences the detectability of extended features even in an ideal case like this.

Fig. B.3 Example of a companion with a synthetic disk with inclination, from top to bottom, of 0°, 45°, 85°.

Fig. B.4 cADI image of HIP 87836. Encircled in blue, the two stars used for the analysis are described in Appendix B.3.

Other factors should be taken into account with regard to the detectability of extended features. Apart from the inclinations that, as shown before, form a critical aspect, also the brightness of the disk and its extension are crucial parameters. In particular, disks around sub-stellar companions are less massive than circumstellar ones. However, their brightness may not be proportional to the quantity of dust grains since they are illuminated both from the inside, by the companion, and from the outside, by the star. Even ifdisks around planets and brown dwarfs may have higher contrast than expected, we still may not be able to detect them due tothe small Hill radii involved. Such compact disks are not resolved and may only be detected by their emission lines.

Fig. B.5 Images of the object P1 around HIP 87836, model and residuals of, from the top, the primary P1, the secondary S1, and the tertiary T1.

Fig. B.6 Contrast curve (upper figure) with and without the P2 point source (S2 was removed before calculating the contrasts and is represented by the red dot) and images of the object P2 (lower figures) around HIP 87836, model and residuals of, from the top, the primary P2, and the secondary S2.

B.3 HIP 87836 data

We also tested the routines described in Sect. 2 on an ad hoc dataset, HIP 87836. This system is particularly suitable for our purposes since it is placed almost on the galactic plane (Galactic coordinates 001.9927977373370 -01.6464052166128); thus the FoV of the instrument is very crowded with PSFs, as shown in Fig. B.4. Due to the projection on the sky plane, a few of these stars are detected very close to each other so that they may resemble multiple bound systems, similar to the ones that we are investigating.

The dataset we analyzed comes from observations taken with SPHERE/VLT in May 2018 with the coronagraph on (see Sect. 3 for the description of the observations and their pre-processing). To test the efficiency of our routines in detecting further companions to directly imaged sources, we chose two apparent multiple systems with different flux ratios and distances between the primary and the secondary that are highlighted by blue circles in Fig. B.4.

The first point source (P1) is placed at r_P1 = 4.11″ and θ_P1 = 209.1° and it is very bright, exhibiting a contrast with respect to the central star F₁ = 2.5 × 10⁻⁴. A second fainter object is placed nearby P1. Our aim is to visualize it in the residuals and, subsequently, to determine its precise position and flux. In order to do that, we ran the single_framebyframe routine to optimize the position and contrast of P1 and, then, to subtract it to obtain the final image shown in the upper panel of Fig. B.5. As expected, the secondary is clearly visible in the upper right corner of the residuals. Moreover, another companion, very close to the primary and not visible by eye in the original image, was unveiled after the subtraction. In order to determine contrasts and positions for the secondary (S1) and the tertiary (T1) we applied twice our routines to subtract each object, in turn. Since S1 and T1 are too faint to optimize (r, θ, F) in each of the unprocessed frames of the scientific cube, we used the single_fullcube routine. Results of the two subtractions are shown in the second and third rows of Fig. B.5, respectively. The contrast with respect to the primary are very similar, F₂ ∕F₁ = 1.7 × 10⁻² and F₃ ∕F₁ = 1.8 × 10⁻², and the spots are placed at r_S1 = 0.13″, θ_S1 = 305.8°, and r_T1 = 0.09″, θ_T1 = 72.8° with respect to P1.

The second point source (P2) is placed at r_P2 = 4.07″ and θ_P2 = 267.1° and is much fainter than P1, with a contrast of F₁ = 1.4 × 10⁻⁵. The secondary (S2) is not clearly distinguishable from the primary. Thus, this system is a perfect test bench for our routines with regard to fainter objects with close companions. We ran the single_framebyframe routine to optimize (r, θ, F) for P2 and to subsequently subtract it, as shown in the upper row of Fig. B.6. In the image of the residuals, S2 is clearly visible and, similarly to what we have done for S1 and T1, we ran the single_fullcube routine to determine position and flux for this object. Results of the second subtraction are shown in the bottom panel of Fig. B.6. The contrast obtained with respect to the primary is F₂ ∕F₁ = 2.1 × 10⁻¹ and the spot is placed at r_S2 = 0.07″, θ_S2 = 99.3° with respect to P2. Together with the images of the two point sources, we show in Fig. B.6 the contrast curves centered on the position of P2 before (blue curve) and after (orange curve) subtracting P2. Such values were obtained subtracting S2 from each raw frame and, then, using the same procedure described in Sect. 2. The red dot in the plot represents the position andcontrast obtained for S2.

Unfortunately, we were not able to find any background star that could be used as a model for each frame. Indeed, the only two bright point sources are placed at the very edge of the FoV so that anisoplanatism affects such PSFs. We tested both objects as model PSFs for the multiple_framebyframe routine but the residuals are higher than the ones obtained with the PSF of the star.

Appendix C DTTS data

The differential tip-tilt sensor of SPHERE is a second-stage tip-tilt sensor located just before the near-infrared coronagraphic wheel (Beuzit et al. 2019). At this location, a grey beam splitter separates a few percent of the light, which is sent to a technical camera sensing precisely the position of the focal spot at a rate of 1 Hz. This enables the compensation of any slow movement between the coronagraph focus and the visible wave-front sensor due to thermal movements, residual differential dispersion, etc. The differential tip-tilt loop of the AO system ensures the final and accurate centering of the star on the coronagraph. Its precision has been assessed in the integration phase and on-sky to be < 0.5 mas, (Baudoz et al. 2010) which means 1/80th of the diffraction width in H-band, for bright stars. For stars as faint in the H-band as DH Tau, the accuracy is probably reduced.

During the observations, the DTTS images are regularly recorded as part of the AO system telemetry. The system records a 32 × 32-pixel, 1-s average image every 30 s. This series of images enables monitoring different information: flux variation, residual jitter of the PSF, influence of the low-wind effect on the PSF (Sauvage et al. 2015; Milli et al. 2018). It is, however, important to note that for faint stars, only the tip of the PSF will be visible in the DTTS images due to the faintness of the star combined with the small fraction of the flux that is taken for the tip-tilt sensing. Therefore, mild low-wind effect may not necessarily be detected using the DTTS images.

We show in Fig. C.1 the DTTS images retrieved for DH Tau for two epochs obtained with SPHERE. While in the first epoch large PSF variations and significant low-wind effect can be seen, in the second epoch, the PSF as resulting from the DTTS images is remarkably stable along the full sequence of observations.

Appendix D LBT data reduction

The observations at the LBT were taken on the 2019-09-19 with the L- and M-band InfraRed Camera (LMIRCam) behind the cryogrenically cooled beam combiner of the Large Binocular Telescope Interferometer (LBTI, Hinz et al. 2016). The LBTI’s primary purpose is to perform nulling interferometry (Ertel et al. 2018, 2020), but the instrument also provides the possibility to obtain among other modes sensitive, high-contrast adaptive optics imaging observationsusing one of the LBT apertures or both independently (e.g., Stone et al. 2018). Here, we used the latter mode.

The observations were taken without a coronagraph so that we could apply the multiple_framebyframe routines using the PSFs of the star as model. For this purpose, short (0.3 s, unsaturated frames) and long (0.9 s, saturated frames) exposure times were alternated during the 2h40min time of the observation, with a stable seeing of ~ 0.65″ and a FoV rotation angle of ~37°. A total of 6620 frames were acquired, with 40 unsaturated frames every 400 saturated ones. The images of the star from the two apertures of the LBT were placed next to each other on the detector at a distance sufficient not to influence the discovery space around the DH Tau system as seen by the left and right mirrors. Nodding in the down-up direction onthe detector (NOD-A and NOD-B position, respectively) was done to subtract the variable sky and telescope background and detector artifacts. During each sequence of 40 unsaturated or 400 saturated frames, the star was placed in the NOD-A position for the first half of the frames and in the NOD-B position for the second half. Due tothe presence of clouds and internal errors, a total of 460 frames were discarded before starting the reduction.

These frames are 2048 × 2048 pixels, joining the images of the star, either in NOD-A or NOD-B position, gathered from the right and left mirrors. Each frame was actually composed of two images, the first of which with a minimum integration time of 27 ms, so that it could be used as a quasi-instantaneous bias. Thus, as a first step of the data reduction, we subtracted the two frames to remove bias drifts at frequencies higher than the nodding frequency. Then, we subtracted to each frame in one NOD position the mean of the closest sequence taken in the opposite NOD. For example, considering the sequence of frames from 1 to 400, we subtracted to each frame in NOD-A (from 1 to 200) the mean of the frames in NOD-B (from 201 to 400). This procedure removes all the remaining bias, the majority of detector artifacts, and the sky- and telescope background.

Fig. C.1 Reconstruction of the PSF of the central star during the observations for the first (left) and second (right) epochs.

The LMIRCam observations are strongly affected by distortion effects, depending on the position of the star in the FoV. We then applied a distortion correction provided by the LBTI team (Spalding & Stone 2019) and available on the LBTO webpage¹.

We divided the frames into four datasets, one for each mirror (left or right) and NOD position, recentering the images on the position of the star and, subsequently, cropping them to a FoV of 7.5″ × 7.5″ (701 × 701). In order to detect the brown dwarf in each frame where the central star was saturated, we stacked together each 20 consecutive images and gathered them into two separate cubes: the scientific one composed of long exposure time frames and the one with unsaturated frames to be used as model PSFs. We thus obtained four scientific cubes, each with 140 frames, and four model PSFs cubes, each with 14 frames.

To further improve the contrast, we also removed the bad pixels with a sigma-clipping and the “striping” that occurs in LMIRCam images due to the vertical detector channels. Finally, we applied the PSF subtraction using the multiple_framebyframe routine, considering one model PSF every ten frames of the scientific cube.

Unfortunately, the NOD-B images for both apertures were obtained with the target close to the edge of LMIRCam’s FoV and we have discovered that uncorrected image distortion appears to affect the image quality in this case. Thus, these frames could not be used in our analysis. The NOD-A cubes were, instead, good-quality observations and we used them both individually and jointly to determine the parameters presented in the text.

References

All Tables

All Figures

	Fig. 1 cADI image of HIP107412, the model, and the residuals.
In the text

	Fig. 2 cADI image of HR8799 b, the model, and the residuals.
In the text

	Fig. 3 Upper panel: cADI image of HD 1160 C, the model and the residuals; bottom panel: cADI image of HD 1160 B, using HD 1160 C as model and the processed image of the subtraction of the two.
In the text

	Fig. 4 Contrast curve around HD 1160 B using one single model PSF with the single_fullcube (blue curve) and single_framebyframe (orange curve) routines, and multiple model PSFs using HD1160 C in place of the star for the multiple_framebyframe routine (green curve). Contrasts are calculated with respect to the central star on the left y-axis and with respect to HD 1160 B on the right one.
In the text

	Fig. 5 cADI image of TYC 7084-794-1 B, the model and the residuals for the first (upper panel, H band) and the second (bottom panel, K band) epochs.
In the text

	Fig. 6 cADI image of TYC 7084-794-1 B, with the model taking into account the smearing of the PSF of the companion and the residuals for the first (upper panel, H band) and the second (bottom panel, K band) epochs.
In the text

	Fig. 7 From the top: epoch 1 (2015-10-26), epoch 2 (2015-12-19), epoch 3 (2018-12-16) in H2-band, and epoch 4 (2020-01-05)in K1-band for DH Tau B as observed with SPHERE/VLT, and the observation (2019-09-19) in the L′ -band obtained with LMIRCam/LBTI. The candidate companion is visible in the right panels of the SPHERE observations. There is no detection for the LBT data.
In the text

	Fig. 8 Color-magnitude diagrams for the IRDIS/SPHERE photometry for the candidate DH Tau Bb compared to a collection of objects from the literature (Bonnefoy et al. 2018).
In the text

	Fig. 9 Detection limits for companions placed within 1″ separation from the host star (upper panels) and beyond 1″ (lower panels, HD 1160 B was moved to the second group for a better visualization of the results). Figures on the left: contrasts with respect to the central star as a function of separation (in arcsec); right: detection limits expressed in mass as a function of separation (in au).
In the text

Fig. 10 Posterior distributions obtained for the frequency of bound binary companions from the MCMC analysis, in cases of 0 (blue) and 1 (red) detection in our survey. Frequencies were derived over separation ranges of 0.5−100 au and mass ratios of q = 0.05−1. The red dashed line shows the peak of the posterior density in the 1-detection scenario, at 7.3%, while only an upper limit was obtained in the 0-detection case. The shaded areas correspond to the 68% confidence intervals.

In the text

	Fig. B.1 Upper panel: cADI images of, from the left, the companion, the model PSF, and residuals. Bottom panel: contrast curves centered on the position of the companion as obtained in the post-processed frame with (blue curve) and without (orange curve) the companion. Contrasts are calculated with respect to the star (left y-axis) and with respect to the companion (right y-axis).
In the text

	Fig. B.2 Injection of a second fake companion at the detection limits placed at 0.05″ (upper panel), 0.075″ (mid panel), and 0.1″ (bottom panel) from the peak of the primary.
In the text

	Fig. B.3 Example of a companion with a synthetic disk with inclination, from top to bottom, of 0°, 45°, 85°.
In the text

	Fig. B.4 cADI image of HIP 87836. Encircled in blue, the two stars used for the analysis are described in Appendix B.3.
In the text

	Fig. B.5 Images of the object P1 around HIP 87836, model and residuals of, from the top, the primary P1, the secondary S1, and the tertiary T1.
In the text

	Fig. B.6 Contrast curve (upper figure) with and without the P2 point source (S2 was removed before calculating the contrasts and is represented by the red dot) and images of the object P2 (lower figures) around HIP 87836, model and residuals of, from the top, the primary P2, and the secondary S2.
In the text

	Fig. C.1 Reconstruction of the PSF of the central star during the observations for the first (left) and second (right) epochs.
In the text