Results from The COPAINS Pilot Survey: four new brown dwarfs and a high companion detection rate for accelerating stars

The last decade of direct imaging (DI) searches for sub-stellar companions has uncovered a widely diverse sample that challenges the current formation models, while highlighting the intrinsically low occurrence rate of wide companions, especially at the lower end of the mass distribution. These results clearly show how blind surveys, crucial to constrain the underlying planet and sub-stellar companion population, are not an efficient way to increase the sample of DI companions. It is therefore becoming clear that efficient target selection methods are essential to ensure a larger number of detections. We present the results of the COPAINS Survey conducted with SPHERE/VLT, searching for sub-stellar companions to stars showing significant proper motion differences (Delta mu) between different astrometric catalogues. We observed twenty-five stars and detected ten companions, including four new brown dwarfs: HIP 21152 B, HIP 29724 B, HD 60584 B and HIP 63734 B. Our results clearly demonstrates how astrometric signatures, in the past only giving access to stellar companions, can now thanks to Gaia reveal companions well in the sub-stellar regime. We also introduce FORECAST (Finley Optimised REtrieval of Companions of Accelerating STars), a tool which allows to check the agreement between position and mass of the detected companions with the measured Delta mu. Given the agreement between the values of the masses of the new sub-stellar companions from the photometry with the model-independent ones obtained with FORECAST, the results of COPAINS represent a significant increase of the number of potential benchmarks for brown dwarf and planet formation and evolution theories.


INTRODUCTION
Direct imaging (DI) is the only detection method that provides observations of an exoplanet or brown dwarf (BD) itself, as it captures the thermal emission of self-luminous companions. With the unique opportunity to obtain photometric and spectroscopic observations of substellar objects, this detection method allows for a direct probe of cold companions atmospheres. DI is also necessary to study the outer regions of planetary systems, that cannot be probed by other detection methods. Despite the remarkable efforts that have been invested in the development of new observing technologies and image processing techniques, and a steady increase in the census of wide-orbit companions, only a handful of systems below the deuterium-burning limit have been uncovered around stars in DI programs, and the occurrences of wide companions appear to be intrinsically low (Biller et al. 2007, E-mail: mariangela.bonavita@open.ac.uk 2013; Lafrenière et al. 2007;Nielsen & Close 2010;Vigan et al. 2012Vigan et al. , 2017Vigan et al. , 2021Rameau et al. 2013;Galicher et al. 2016a;Nielsen et al. 2019). In order to empirically constrain the formation, evolution, and atmospheric properties of both isolated and bound sub-stellar companions, we need to uncover a substantial population of these objects, and measure their fundamental properties, such as the effective temperature and mass. However, even when a comprehensive view and an extensive spectro-photometric characterisation is possible, imaging surveys still only provide measurements of an object's luminosity. Mass estimates for imaged planets and brown dwarfs therefore rely entirely on evolutionary models, which currently carry high uncertainties, particularly at young to intermediate ages. An independent determination of masses from dynamical arguments is therefore crucial to overcome the large uncertainties introduced by evolutionary models, and in turn refine the theories. Furthermore, a sample of benchmark objects should ideally span a wide range of properties (e.g., spectral types, masses, ages). As direct imaging is more amenable to very young systems, where low-mass companions are still bright, the majority of the scarce sample of such benchmark objects in the planetary regime orbit relatively young stars (few tens to hundred Myr). On the contrary, most brown dwarf companions with well-defined dynamical masses orbit older hosts with field ages of several Gyrs (see Fig. 1). As known bound sub-stellar companions are even rarer in associations such as Hyades (∼650 Myr Martín et al. 2018), the intermediate age regime remains relatively unexplored, and theoretical models are hence particularly poorly constrained at these ages. Our understanding of the origins and atmospheres of these objects thus remains severely limited by the small number of known systems. In particular, the sparse sample of directly-imaged companions show a large diversity of spectro-photometric characteristics and orbital configurations, and remain challenging to grasp as populations (Bowler 2016). Larger numbers of detections are hence essential to enable a better characterisation and understanding of the wide-orbit companion population, and obtain a clearer picture of their formation patterns. In Fontanive et al. (2019), we presented a new tool COPAINS (Code for Orbital Parametrisation of Astrometrically Inferred New Systems), developed to identify previously undiscovered companions detectable via DI, based on changes in stellar proper motions across multiple astrometric catalogues. A significant proper motion difference (∆µ) between two catalogues for a given star is a good indication of the presence of a perturbing body. For systems showing significant differences between proper motions measured over a long time baseline (e.g., Tycho-2, or Tycho GaiaAstrometric Solution -TGAS; Høg et al. 2000;Michalik et al. 2015), and catalogues that provide short-term proper motions (e.g., Hipparcos, GaiaDR2;ESA 1997;Gaia Collaboration et al. 2018), the tool allows for the computation of secondary mass and separation pairs compatible with the observed trend, marginalised over all possible orbital phases and eccentricities. The resulting solutions are based entirely on dynamical arguments, although a dependence on the adopted (usually model-derived) stellar mass remains in the obtained secondary masses. Compared to the expected sensitivity of an imaging instrument, these predictions can then be used to select the most promising targets for DI searches of low-mass companions. The use of such informed selection processes had already proven to be effective in the stellar regime (Makarov & Kaplan 2005;Tokovinin et al. 2013;Bowler et al. 2021;Steiger et al. 2021), and recently led to the discovery of new brown dwarf companions based on the astrometric signatures induced on their host stars (Currie et al. 2020;Chilcote et al. 2021). Such astrometric systems are particularly valuable, as the combination of relative astrometry from DI information with absolute astrometry from the primary's astrometric signature offers a remarkable opportunity to refine orbital constraints and measure dynamical masses (Calissendorff & Janson 2018;Snellen & Brown 2018;Brandt et al. 2019;Dupuy et al. 2019;Grandjean et al. 2019;Maire et al. 2020;Nielsen et al. 2020;Drimmel et al. 2021). Precise orbital elements for the population of wide-orbit companions can provide key insights into formation mechanisms (e.g., Bowler et al. 2020). Furthermore, model-independent mass measurements for brown dwarfs and giant planets are especially important to bypass the use of mass estimates from theoretical models, which typically carry large uncertainties, both due to the difficulties of determining system ages, and to the systematic uncertainties of the evolutionary and atmosphere models which are particularly pronounced at the lowest masses and youngest ages. Increasing the pool of systems amenable to dynamical mass measurements will therefore be essential to help calibrate theoretical models for substellar objects. In this paper, we present the results of a pilot survey conducted with the SPHERE instrument (Beuzit et al. 2019), an extreme adaptive optics facility at the ESO's Very Large Telescope (VLT), which employed the COPAINS tool for informed target selection. We describe the sample and selection method in Section 2. The observations and data reduction are presented in Section 3. The survey results are reported in Section 4 and discussed in Section 5.

Initial Target List
To select targets for DI campaigns, we searched different catalogues containing proper motions, using as an input a list of known, relatively young, sources, from which the targets protected by the SPHERE Guaranteed Time Observations (GTO) were removed. For the sources with both long-and short-term proper motion information, we first selected stars showing a difference larger than 3σ between two catalogues, in either of the proper motion components. The data presented here were obtained during three ESO periods (P100, P102, and P104) 1 , with several differences in target selection procedure, as a result of the different Gaia data releases available at the time of each selection. The relevant information about the two initial selections are listed here.

Selection 2 (P104)
• Input list: The same compilation as above (Desidera et al. 2015), complemented by bona-fide members of nearby young moving groups . In total, there are about 2200 unique objects in this input list.
• Long-term proper motion catalogue: TGAS, with a proper motion baseline of ∼25 yr.
Following the selection of these so-called ∆µ candidates, stars with known companions from visual and/or radial velocity observations were removed from the list. A star was removed if satisfying any of the following criteria: • Multiplicity flag (MultiFlag) equal to C or O in the Hipparcos catalogue (ESA 1997); • Star appears in the Catalog of Components of Double & Multiple stars (CCDM; Dommanget & Nys 2002); • Star appears in The ninth catalogue of spectroscopic binary orbits (S 9 B ; Pourbaix et al. 2004); • Star has a sub-stellar companion listed in The Extrasolar Planets Encyclopaedia 2 ; • A suffix in the name of the star indicating that it is part of a binary system, or an unresolved binary (e.g. A, B, AB).
• For the Selection 2, we also excluded stars with another source closer than 5 on the sky, and sharing the same parallax (within 3σ), in order to exclude obvious stellar binaries from the survey and correctly focus on sub-stellar companions. For this, we only used the sources with a relatively small parallax uncertainties ( / σ < 0.2).

COPAINS Selection
After removing known binaries as detailed above, we further retained only the targets for which the COPAINS tool (Fontanive et al. 2019) returned that the objects causing the observed astrometric trends could be sub-stellar, if in the parameter space detectable with SPHERE (i.e., within the instrument field-of-view and above the expected sensitivity of observations). COPAINS provides a good indication of the region of the parameter space in which a hidden companion responsible for an observed astrometric offset may be. Based on the formalism from Makarov & Kaplan (2005), the approach assumes that long-term proper motion measurements are representative of a system's centre-of-mass motion, and that short-term measurements correspond to the instantaneous reflex motion of the host star. For a measured ∆µ value, the code allows us to evaluate the possible companion mass and separation pairs compatible with the astrometric data, for a given distance, stellar mass (see Section 2.3), and eccentricity distribution, while assuming face-on orbital inclinations (see Fontanive et al. 2019 for details). We adopted a Gaussian eccentricity distribution centred around e=0 with a width of 0.3 (Bonavita et al. 2013), and the astrometric catalogues listed above were used as long-and short-term proper motion values. Figure 1 shows examples of the resulting trends computed with COPAINS for one target from P100 (left) and one target from P104 (right). In each case, the estimated solutions were compared to the expected sensitivity limit of SPHERE-IFS, shown as a red solid line in Figure 1. The IFS limits were obtained following the approach described by Mesa et al. (2021) and converted to minimum mass limits using the models from Baraffe et al. (2015) and the values of the age and stellar mass available at the time of the selection. Only promising targets, where the intersection of the detection limits and computed regions suggested possible sub-stellar companions detectable in our survey, were kept for the final sample, based on a visual analysis of the obtained plots. Finally, only the targets not previously observed with SPHERE, and observable in the relevant ESO period were kept. The final target list consisting of 25 stars observed with SPHERE is given is Table 1, while Table 2 lists the values of the proper motions from Tycho-II, TGAS, Gaia DR2 and EDR3 and the resulting ∆µ.

Stellar Ages
A full revision of the stellar ages was performed for all the targets following the approaches described in Desidera et al. (2015) and Desidera et al. (2021), considering a variety of indicators and also performing a check for membership to groups using the BANYAN Σ on-line tool 3 . Most of the targets were found to be field objects, while still compatible with a relatively young age, or 3 http://www.exoplanetes.umontreal.ca/banyan/banyansigma. php stars with ambiguous membership. In these cases, we considered indirect age indicators such as the equivalent width of 6708Å Lithium doublet, rotation period, X-ray emission, chromospheric activity, taking as reference the empirical sequences of members of groups and clusters (e.g., Desidera et al. 2015) and isochrone fitting. Several of our targets have ages between Hyades and Pleiades. In this range, we took advantage of the recently derived rotation sequence for Group X at an age of 300 Myr (Messina et al. 2022). When applicable, we considered the age indicators for physical companions outside the SPHERE field of view and, for close binaries, we deblended photometric colors for binarity, to improve the reliability of the derived ages. For the few target found to be belonging to young moving groups, we adopted the values of the ages presented in Bonavita et al. (2016) and Desidera et al. (2021), mostly based on Bell et al. (2015). All targets were found to be younger than 1 Gyr, except for GJ 3346 for which the age is likely to be closer to 5 Gyr as discussed in detail in Bonavita et al. (2020). The age determination process for each target is discussed in Appendix A.

Stellar Masses
The masses for all the stars in the sample were also revised, using the Manifold Age Determination for Young Stars ( , Squicciarini & Bonavita in preparation; see Squicciarini et al. 2021, for a description of the tool) and the updated values of the stellar ages. M retrieved and cross-matched photometry from Gaia EDR3 (Gaia Collaboration et al. 2021) and 2MASS (Skrutskie et al. 2006) for all our targets and then applied a correction for interstellar extinction by integrating along the line of sight the 3D extinction map by (Leike et al. 2020); the derived A(G) were turned into the photometric band of interest using a total-to-selective absorption ratio R=3.16 and extinction coefficients A λ from Wang & Chen (2019). The derived absolute magnitudes were then compared with a grid of isochrones with an age range based on the minimum and maximum age values included in Table 1 to yield a mass estimate. M can use several available grids, but in this instance the PARSEC isochrones (Marigo et al. 2017) were used, due to their large dynamical range spanning the entire stellar regime. A constant solar metallicity, appropriate for most nearby star-forming regions, was assumed (D'Orazi et al. 2011) for all targets except for HIP 21152 and HIP 21317, for which we assumed [Fe/H]=+0.13 based on their membership to the Hyades (see Appendix A for details). For each star, a sample of mass estimates was constructed by computing the best-fit mass at different ages within the given age range; its median was taken as the final mass estimate, while the reported errors represent the 16th and the 84 percentile; photometric uncertainties were naturally propagated on the final result via a Monte Carlo approach, i.e. by randomly varying, in a Gaussian fashion, photometric data according to their uncertainties while building the sample of mass estimates.

OBSERVATIONS AND DATA REDUCTION
All observations were performed with VLT/SPHERE (Beuzit et al. 2019) with the two Near Infra-Red (NIR) subsystems, IFS (Claudi et al. 2008) and IRDIS (Dohlen et al. 2008) observing in parallel (IRDIFS Mode), with IRDIS in dual-band imaging mode (DBI; Vigan et al. 2010). For all targets we used the IRDIFS-EXT mode, which enables covering the Y-, J-, H-, and K-band in a single observation, which is meant to provide a high-level of spectral content for subsequent analyses. A summary of the observing parameters and conditions is given in Figure 1. Examples of the solutions computed with COPAINS for the observed astrometric trends of one target selected through the Selection 1 process (using Tycho-2 and TGAS; left) and one target from the Selection 2 process (using TGAS and GaiaDR2; right). The solid black lines corresponds to the median curves of solutions, and the dark and light shaded areas represent the 1 and 2-σ confidence intervals, respectively. The dashed grey lines show the expected detection limit used for the survey selection, computed as detailed in the text. Figure 2. IFS detection limits expressed in contrast (left) and minimum companion mass (right) vs projected separation, for all the targets in our sample. The dashed vertical line marks the coronographic radius. Note that objects with multiple epochs will appear more than once. The COND models (Baraffe et al. 2003) were used for the magnitude to mass conversion, using the adopted age from Table 1. Table 3. The observing sequence adopted was similar to those designed for the SHINE Guaranteed time survey (see e.g. Chauvin et al. 2017) and consisted of: • One PSF sub-sequence composed of a series of off-axis unsaturated images obtained with an offset of ∼0.4 relative to the coronagraph center (produced by the Tip-Tilt mirror). A neutral density filter was used to avoid saturation 4 and the AO visible tip-tilt and high-order loops were closed to obtain a diffraction-limited PSF.
• A star center coronagraphic observation with four symmetric satellite spots, created by introducing a periodic modulation on the deformable mirror (see Langlois et al. 2013, for details), in order 4 www.eso.org/sci/facilities/paranal/instruments/sphere/ inst/filters.html to enable an accurate determination of the star position behind the coronagraphic mask for the following deep coronagraphic sequence.
• The deep coronagraphic sub-sequence, for which we used here the smallest apodized Lyot coronagraph (ALC-YH-S) with a focalplane mask of 185 mas in diameter.
• A new star center sequence, a new PSF registration, as well as a short sky observing sequence for fine correction of the hot pixel variation during the night.
IRDIS and IFS data sets were reduced using the SPHERE Data Reduction and Handling (DRH) automated pipeline (Pavlov et al. 2008) at the SPHERE Data Center (SPHERE-DC, see 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. In order to calibrate the IRDIS and IFS data sets on sky, we used images of the astrometric reference field 47 Tuc observed with SPHERE at a date close to our observations. The plate scale and true north values used are based on the long-term analysis of the GTO astrometric calibration described by Maire et al. (2016).

Detection Limits
In order to evaluate our sensitivity to stellar companions, we determined detection limits for point sources. We used the standard procedure to derive detection limits outside the coronagraphic field masks that makes use of the SPECAL software as described in Galicher et al. (2018) Table 4. The corresponding values of the minimum companion mass limits, obtained using the evolutionary models from Baraffe et al. (2015) for the magnitude to mass conversion, are shown in the right panel Fig. 2. The achieved average limits appear to be significantly worse than the average of the expected limits (dashed-dotted line in Fig. 2) used for the target selection. This is most likely due to the fact that, since our program was executed in service mode and as filler, most of our targets were observed in sub-optimal conditions and with very small field rotation, with a strong negative effect on the quality of the high contrast imaging performances, especially at short separations. This is also confirmed by the fact that the average limit achieved for the first 150 targets of the SHINE survey (gray solid line, from Vigan et al. (2021)), where all target were observed in the best possible conditions, is instead much better than both the measured and estimated COPAINS limits.

RESULTS
We detected a total of 14 candidate companions, 4 of which were found to be background sources thanks to additional epochs available in the literature or obtained in our program (HIP 63862, HD 57852, HIP 33690). The 10 co-moving companions are shown in Fig. 3. Eight of the co-moving companions have separations below 0.9 arcseconds, and are therefore in the IFS field of view, while the remaining 2 were only observed with IRDIS. Five are new discoveries, including a white dwarf companion at ∼ 3.6 from GJ 3346 (already presented in Bonavita et al. (2020)) and four new sub-stellar companions: HIP 21152 B 5 , HIP 29724 B, HD 60584 B and HIP 63734 B. Their properties are discussed in details in Sec. 4.5, 4.6 and 4.7.

SPHERE astrometry and photometry
The astrometry and photometry measurements from all our SPHERE observations, are listed in Table 5. For each epoch we report the 5 HIP 21152 B was independently discovered as part of two other surveys targeting accelerating stars, as detailed in Kuzuhara et al. (2022) and Franson et al. 2022 (in preparation). Both works include an in-depth characterisation of the system, the latter also including a full spectral and orbital analysis combining all available data sets. projected separation and position angle, and the contrast (expressed as apparent magnitude difference) in the IFS Y and J filters, as well as the IRDIS K 1 and K 2 for the IRDIFS-EXT observations 6 . The probability that the source is a background star, evaluated as described in Sec. 4.2, is also listed.

Common proper motion confirmation
Multiple epochs, used to clarify the bound or background nature of our candidates, were available for 9 of the program stars. Except for HIP 15247 and HD 60584, which were re-observed with SPHERE as part of the program, all additional epochs were retrieved from other surveys, catalogues (including Gaia) or papers dedicated to specific objects. The complete list of astrometric measurements for all our systems is presented in Table 7, together with the references used for each entry. Figure 4 shows the resulting common proper motion analyses for both the co-moving and background interlopers. Note that for HD 60584 we detected two sources, one at 0.5 arcsecs (CC1, discussed in Sec. 4.7) and one at ∼3 arcsecs (CC2), shown in Fig. 4 and confirmed to be background. The remaining candidates are bright companions at very small separation and are then very likely physically related as the probability of having such bright background stars at these separations is very low. To confirm this, we used the code described in Section 5.2 of Chauvin et al. (2015) and adapted it to our results to estimate the probability of finding a background contaminant at the given separation and contrast as a function of galactic coordinates by comparison with the prediction of the Besançon galactic model (Robin et al. 2012). As expected all the resulting probabilities (reported in the last column of Table 5) are below 10 −4 . A further confirmation of common proper motion is also provided by the agreement between the properties of the imaged companions with observed ∆µ, as discussed in Sec. 4.4.

Companion spectra
Spectra for the bright companions (providing relatively high SNR data) were obtained using the IFS data. We used the SpeX Prism Library Analysis Toolkit (SPLAT; Burgasser & Splat Development Team 2017) to estimate the spectral classification of each companion. We used the built-in spectral fitting function in SPLAT to compare, by minimising the χ 2 value, the observed spectra to both the SPLAT near-infrared spectral standards and to the full library of templates available in SPLAT. For bright companions, we simply used the spectrum obtained by rotating and summing the images. For HIP 21152 B, which was not detectable without removing the speckle pattern, we injected negative point spread functions (derived from the flux calibration) on the individual monochromatic images at the average companion position, and changed its intensity minimising the root mean square of the residuals in area of 9 × 9 pixels centred on this mean position. Figure 5 shows the spectral standard (red) and template (blue) providing the lowest χ 2 values to the observed spectrum of each detected companion. Note that the images of HIP 15247 B was saturated at most wavelengths, so it was not possible to obtain an usable spectrum. Among the five remaining companions, four have an estimated mid-to late-M spectral type, while HIP 21152 B is clearly a substellar object and compatible with a late-L or early-T spectrum. 6 see www.eso.org/sci/facilities/paranal/instruments/ sphere/inst/filters.html for a full description of the SPHERE filters. The strong absorption feature visible in the spectrum of HIP 78549 B is most likely spurious and the result of the extremely low quality of the wavelength calibration available for this object, particularly between 1.32 and 1.45 µm. We therefore masked this region of the spectrum while performing the fit. The resulting mid M spectral type is in agreement with the properties estimated from the photometry.

FORECAST (Finely Optimised REtrieval of Companions of Accelerating STars)
Since the ∆µ can be considered as an approximation of the instantaneous acceleration due to an unseen companion, it can be represented as a vector in the plane of the sky directed toward the position of the companion -at the epoch of the latest astrometry observation. The position angle of an imaged companion compatible with the ∆µ should then be along the same direction of the acceleration, plus or minus the change in angle due to the orbital motion of the companion between the latest astrometry epoch and the imaging one. Based on these considerations, it is possible to highlight a region on the plane of the sky where the companion compatible with the ∆µ should lie, based on the ∆µ orientation, and also associate a value of the companion mass corresponding to each point in the resulting 2D map, based on the ∆µ absolute value. An example of the ∆µ maps obtained with this method using the ∆µ at the epoch of Gaia EDR3 (2016.0) is shown in Fig. 6, while a complete description of the method, and of the FORECAST (Finely Optimised REtrieval of Companions of Accelerating STars) code used to produce the maps will be the subject of a dedicated companion publication (Bonavita et al. 2022, in preparation). In each image, the position of the detected companion is marked as a star. As all companions lie within the allowed region (as shown in the left panel of Fig. 7, which compares the PAs of the companions and those of the ∆µ vectors), we could also obtain a first estimate the expected companion mass at that position and compare it with the values of the value of the mass derived from the photometry. The dynamical masses (M ∆µ ) were estimated using the method described in Kervella et al. (2019), which allows to also take into account the effect of the unknown orbital eccentricity and inclination, as well as that of the observing window smearing, which are instead not taken into account by the COPAINS tool. The comparison with the photometric masses (M Phot ) is shown in the right panel of Fig. 7, with the error bars on the values of M ∆µ taking into account the errors on the ∆µ and companion position, and those on M Phot mainly arising from the uncertainties on the stellar ages. The good agreement between the values shows the potential of our method in providing a high number of potential new benchmark objects.

HIP 21152 B: a new bound substellar companion in the Hyades
The companion to HIP 21152 is clearly the most interesting of our high SNR detections. The IRDIS photometry and the IFS spectra, shown in Fig. 5, point towards an early T spectral type. The  Table 7). The blue line shows the motion of a background object relative to the target, based on the EDR3 parallax and proper motion of the primary over the same time frame. In all panels, the filled circles show the measured separation and position angle of the companions at each epoch, colour coded as explained in the legends. The crosses indicates the expected position of a background object at the same epoch, following the same colour code. The companions for HIP 28474, HIP 78549, HIP 15274, HIP 71899 and GJ 3346 (updated from the dedicated work published in Bonavita et al. (2020), which only included GaiaDR2 epoch), are clearly found to be co-moving with our targets.
mass of the object, as judged from its spectral type and age, and adopting the models from Baraffe et al. (2015) 7 , is estimated to be 0.032 ± 0.005 M . This is also in good agreement with the mass derived from the FORECAST analysis, which is 0.021 ± 0.007 M . The Lithium-depletion boundary, an observational limit that separates low-mass stars from brown dwarfs based on their ability to burn Li in their cores, has been estimated to occur around spectral type L3.5 -L4 in Hyades (Martín et al. 2018). Currently, about a dozen objects with spectral type L3.5 or later claimed as members of the cluster. Lodieu et al. (2014) published spectra of 12 L-dwarf candidates from Hogan et al. (2008), and confirmed one of them as a potential brown dwarf member of Hyades, with spectral type 7 New atmospheric models have been recently developed by Phillips et al. (2020). However, such models (known as ATMO2020) are mostly valid for relatively old ultra-cool objects (late-T and Y companions). Therefore, given the age and brightness of HIP 21152 B as well as the estimated spectral type, we decided that the models from Baraffe et al. (2015) would be more suited for the mass estimate in this case. The same applies to the other BD companions described in Sec. 4.6 and 4.7 L3.5. The same objects has later been listed as a L4 member of Hyades, with Li absorption detected in the spectra by Martín et al. (2018). Four more L-type members (two L5 and two L6) have been confirmed in Schneider et al. (2017); Pérez-Garrido et al. (2017,2018). Bouvier et al. (2008) were the first to report the existence of T-type candidate members in Hyades. The membership of one of them (CFHT-Hy-21) has been confirmed by Lodieu et al. (2014), whereas the other one (CFHT-Hy-20) is still considered a candidate (Lodieu et al. 2014;Zhang et al. 2021), given the lack of radial velocity measurement. Four T-type probable members have been reported in Zhang et al. (2021), spanning the spectral type range T2 to T6.5. The multiplicity of Hyades stars has been assessed through imaging (Reid & Gizis 1997;Patience et al. 1998;Siegler et al. 2003;Duchêne et al. 2013), high-resolution spectroscopy (Reid & Mahoney 2000), a combination of the two (Guenther et al. 2005), and Gaia astrometry (Deacon & Kraus 2020). While several binary systems consisting of two substellar objects, or objects close to the (sub)stellar border, were reported in Siegler et al. (2003); Reid & Mahoney (2000); Duchêne et al. (2013), no substellar companions to FGK stars have so far been reported in Hyades, making HIP 21152 B Figure 5. Comparison between the spectra of five of the companions detected at high SNR (black bullets with error bars) with the best fit from SPLAT (Burgasser & Splat Development Team 2017), showing in red the best fit among the spectral standards provided in SPLAT, and in blue the overall best fit among all available templates.
the first of its kind. HIP 22152 B thus occupies a unique place in the luminosity-age parameter space compared to the known population of directly-imaged systems, making it a highly valuable benchmark for empirical constraints to theoretical models. BD companions such HIP 22152 B have been shown in the past to be themselves close pairs of substellar objects (see e.g. Martín et al. 2000;Potter et al. 2002) Table 2. The position of the companions is marked with a red star. Colours are according to the dynamical mass responsible for the ∆µ at a give distance; the same logarithmic scale was used for all stars, according to the colour scale shown on the bottom of the figure. The empty area at center is the area covered by the coronagraphic mask.
corresponds to a projected separation of roughly 82 mas. While this target has good potential for the detection of additional (possibly less massive) companions, an analysis like the one described by Lazzoni et al. (2020) will require higher quality data than those presented in this work.

HIP 29724 B: a new high mass brown dwarf
We also detected a very close (99.9±1.5 mas) companion to HIP 29724. The photometry suggests a mass of 0.063±0.008 M , thus placing HIP 29724 B also in the sub-stellar regime. The FORE-CAST analysis confirms that the companion is compatible with the observed TGAS-EDR3 ∆µ, both in terms of position and corresponding mass (0.067 ± 0.021 M ). With a projected separation of just 6.3 au, RV monitoring could significantly contribute to the refinement of the BD orbit and dynamical mass. However, there are no high-precision RV measurements available up to now. Sparse RV data over several decades (Nordström et al. 2004;Torres et al. 2006 Table 8 for all companions. In both panels blue symbols are high SNR detections; red symbols are low SNR candidates that needs further confirmation.

HD 60584 B and HIP 63734 B: two more potential low mass brown dwarf candidate companions
The ∆µ maps obtained with FORECAST are not only useful to confirm the nature of the candidates found in objects with highly significant ∆µ, but could in principle also be used as finding charts to highlight and retrieve possible additional companions appearing in the imaging data at a SNR lower than the threshold required for a confirmed detection, which would otherwise be overlooked. We had in fact noticed that all our high SNR detections were around stars for which the ratio between the ∆µ absolute value and its error (hereafter S NR ∆µ ) is higher than 10. So we decided to produce the FORECAST maps also for the targets with low or intermediate S NR Max (HD 60584, HIP 63734, and HIP 22506) as a test for their possible use as finding charts. While the data for HIP 63734 were taken in fairly good atmospheric conditions, the target was observed quite far from meridian passage so that field rotation is limited (only 7.13 degrees). The analysis of the IFS data with ASDI-PCA revealed a candidate companion with S NR ∼ 8.5, at separation 555 ± 2 mas and PA=329.6 ± 0.2 degree, with a contrast of dJ=12.08 and dH=11.26. The comparison with Baraffe et al. (2003) models yields an evolutionary mass of 0.011 ± 0.005M . While the PA of the object is fully compatible with that of the ∆µ, the corresponding companion mass is a bit higher (0.032 ± 0.020M ) is a bit higher than the evolutionary one. However, this detection is uncertain because the small rotation angle makes the noise distribution quite different from a Gaussian. The case of HD 60584 is slightly more complicated. Although there were four available epochs, the resulting IFS data were all of relatively poor quality. The second and third epochs were taken at about one month interval; since we did not expect a large orbital motion between them, we combined them as a single observation. This allowed us to identify a point source compatible with the FORECAST predictions at a separation 543±5 mas and PA=52.5±0.5 degrees, with a SNR of 4.9. At that position, the mass of a companion compatible with the ∆µ would have to be below 0.01 M which, at the relatively high age of the targets would mean that the companion would have to be fairly faint. While the IFS photometry points towards a higher mass (0.028 ± 0.009M ), given the high uncertainties due to the image and calibration quality and the very field small rotation angle we still deemed it acceptable. As for the other BDs, no evident elongation pointing towards possible additional companions was observed for the PSF of HIP 63734 B and HD 60584 B. However, given the low SNR of the detections, we estimate that the only potentially detectable companions within the Hill Radius of these BDs (103 and 80 mas, respectively) would have been equal-mass ones. Although further investigation is required to confirm the nature of both HD 60584 B and HIP 63734 B, mostly due to the poor quality of the imaging data, these additional detections clearly show the power of the approach, which pushes the ∆µ method towards companions with smaller masses, whose detection is more uncertain and likely to be below the usual 10 σ threshold used for automatic retrieval of point sources in imaging data. Although repeated observations will be needed to confirm low SNR imaging detections, with the use of FORECAST a single observation is potentially enough to confirm that a companion observed in high contrast imaging is the one responsible for the ∆µ, and therefore co-moving. Other low SNR point sources were retrieved around HIP 22506, but were discarded mostly based on the strong disagreement between their brightness and the value of the mass compatible with the position within the FORECAST map. This also shows the power of FORECAST in terms of vetting of possible background or spurious sources.

DISCUSSION
The COPAINS selection method has proven very successful in ensuring a high detection rate in both the stellar and sub-stellar regime.
We detected a total of 14 candidate companions. Two were known binaries (HIP 15247 B and HIP 78549 B) and four were identified as background sources thanks to additional available epochs found in the literature, which also allowed us to confirm the common proper motion nature of four additional new stellar companions, including the white dwarf companion to GJ 3346, described in Bonavita et al. (2020). The masses of the remaining four candidates, derived using the available photometry and the evolutionary models from Baraffe et al. (2015), place them in the sub-stellar regime. Such high sub-stellar companion detection rate confirms the efficiency of the COPAINS selection method, as well as the new FORECAST code used to confirm their nature, based on the agreement between their properties with the predictions based on the measured ∆µ.

Comparison with blind surveys
To quantify the improvement in terms of detection rate compared to blind surveys, we compared our results with those from the first 150 targets from the SHINE survey (Vigan et al. 2021). Only 93 of the SHINE-150 targets are included in TGAS and DR2, and only 13 have ∆µ more significant than 3 σ, and would have therefore been selected for COPAINS. Three of these SHINE ∆µ stars have sub-stellar companions, including HIP 65426, one of the two new SHINE detections. The resulting sub-stellar companions frequency is then ∼ 25% which is significantly higher than what obtained with the blind approach (∼ 9% without any correction due to prior knowledge about the companions, see Vigan et al. (2021) for details), once again showing the efficiency of the COPAINS selection method, especially when combined with the FORECAST maps.

Limitations
We have demonstrated with our campaign the power of using informed target selection processes such as the COPAINS tool (Fontanive et al. 2019) to identify new directly-imaged companions. The high detection rate obtained here strongly validates the use of such approaches in survey designs, despite the numerous assumptions and limitations of the work conducted in this pilot survey. Indeed, our original sample selection considered catalogue proper motions taken at face value. Instead, a more accurate approach would require placing all measurements in the same reference frame and at the same epoch, such as the Hipparcos-Gaia Catalog of Accelerations (HGCA) defined by Brandt (2018Brandt ( , 2021. These catalogues provides Hipparcos and Gaia DR2/EDR3 proper motions, as well as a Gaia-Hipparcos scaled positional differences (close to the TGAS proper motions), placing all proper motions at the epochs of GaiaDR2 and EDR3, respectively, with recalibrated uncertainties. Nonetheless, given the number of approximations made in COPAINS, these differences were found to be negligible, especially given the use made of the resulting computed trends (i.e., visual selection based on comparisons to expected detection limits). As COPAINS considers long-term proper motions as representative of the system's center-of-mass motion, and short-term measurements as the instantaneous reflex motion of the host, the considered ∆µ are only good approximations for systems with orbital periods that roughly match the long-and short-term timescales of the considered astrometric catalogues (see Fontanive et al. 2019). These limitations are further added to the adopted eccentricity distribution, primarily impacting the width of the computed solutions, and the fact that the approach assumes face-on orbits, which implies that estimated trends actually provide lower mass limits for a given separation. While information from the now available HGCA catalogues (Brandt 2018(Brandt , 2021 or proper motion anomalies measured by Kervella et al. (2019Kervella et al. ( , 2022 would therefore provide more robust and reliable astrometric trends to use for selection purposes (and should consistently be used for orbital and dynamical mass constraints), our results were not impacted by the use of uncorrected catalogue values, and we have nonetheless shown that the idea behind our method offers a highly promising pathway for future observing programs.

CONCLUSIONS
We presented the results of the COPAINS survey, a search for companions to 25 stars selected using the COPAINS tool by Fontanive et al. (2019), and the discovery of four new brown dwarf companions: HIP 21152 B, HIP 29724 B, HD 60584 B and HIP 63734 B. The value of blind surveys of course lies in their ability to constrain the underlying planet and sub-stellar companion population, but comes with a high cost in terms of telescope time. On the other hand, surveys like COPAINS offer an undeniably efficient selection method, providing a much higher success rate with a considerably smaller time commitment. Moreover, the possibility to derive model independent mass estimates for the companions to accelerating stars also means that each new detection arising from surveys like CO-PAINS can be added to the currently scarce number of much needed benchmark objects, providing new crucial constraints the evolutionary models. Given the good agreement between the values of the masses of these companions obtained from the photometry with the model-independent ones based on the ∆µ, this work represents a considerable addition to the current benchmarks sample. The combination with the FORECAST maps ensures a high detection rate even when the quality of the imaging data is not ideal, while at the same time further enhancing the sensitivity to lower mass companions. Fig. 8 shows the mass ratio (left panel) and companion mass (right panel) vs separation of the COPAINS sub-stellar candidates compared to those of the known DI companions (data from exoplanet.eu updated on April 4 th , 2022). An estimate of the ∆µ was possible for about 40% of the known companions and about half of these (shown as coloured dots in Fig. 8) would have been selected using COPAINS. This plot once again shows how a selection like the one provided by COPAINS can already lead to the detection of companions well in the planetary mass regime. Catalogues like those by Brandt et al. (2019) and Kervella et al. (2022) provide a more robust estimate of the accelerations, and therefore their use with COPAINS is likely to lead to an even more effective selection. Moreover, these acceleration catalogues can also be used to retrieve a time series of absolute astrometry which, combined with the imaging data, allows for a more detailed characterisation of the orbit, and thus of the dynamical mass (see e.g. Drimmel et al. 2021). The availability of Gaia-only accelerations in the upcoming full third GaiaRelease (DR3) will allow for another great step further. It will in fact not only free the method from the boundaries so far imposed by the use of external catalogue as source of long term proper motion, but more importantly will allow for a significant improvement in terms of uncertainties, thus allowing future survey to truly focus on targets with accelerations caused by planetary mass companion.   ID  pmRA  pmDE  ∆µ  TYCHO-2  TGAS  DR2  EDR3  TYCHO-2  TGAS  DR2  EDR3  TYC-TGAS TGAS-DR2  TGAS-DR3  HIP 15247 77.      Table 6. Gaia astrometry and photometry of companions retrieved in EGDR3, including additional companions outside SPHERE FoV. Separations and position angles were derived using the positions from Gaia EDR3, when available. Further details about the known systems can be found in Appendix A. The ID of the target is reported in the first column to be consistent with the rest of the tables in the paper. The Gaia ID, and the values of the parallax and proper motion reported are those retrieved in EDR3 for the companions. Although the companions of HIP 15247 and HD 129501 were detected by Gaia, no astrometric solution was available in EDR3 (hence the blank fields).    Bonavita et al. (2020), the masses of the primaries (M A ) were derived as described in Sec. 2.3. The same is true for M B (phot), which is the value of the secondary mass obtained from the photometry, using the COND models (Baraffe et al. 2003

DATA AVAILABILITY
The data used for this work are available through the ESO Science Archive Facility (http://archive.eso.org/cms.html).
(2020). It was included in the sample as it appears as a young star with moderately high activity level, due to accretion of material lost by the WD progenitor.
HIP 27441 = HD 39126 = TYC 7601-371-1 Field object. Age indicators consistently provide an age of 250±100 Myr. The star has a wide common proper motion companion, 2MASS J05483951-3956087, at 34" (1400 au at the distance of the star). From colors and absolute magnitudes, the companion is expected to be a M2 star.

HIP 28474 = HD 41071 = RT Pic
The star is classified as member of Columba association in several works and membership is confirmed by our analysis with updated Gaiainputs. Age indicators support a young age although Li EW and rotation period are in mild disagreement with the typical values of Tuc-Hor and Columba members and more similar to Pleiades and AB Dor ones. For this reason, we adopt Columba age but with upper limit at the age of the Pleiades. The classification as eclipsing binary with 2.6 mag deep eclipses (Malkov et al. 2006) appears spurious, as there are no indications of them in the TESS light curve (4 sectors) and the RV results constant over decades at km/s level (Nordström et al. 2004;Gaia Collaboration et al. 2018) and at few tens of m/s level over three months from HARPS data. A close stellar companion was identified by Chauvin et al. (2010) and further observed by Tokovinin & Briceño (2020). It is also detected in our data (mass from photometry 0.16 Msun) and is the responsible of the ∆µ signature.
HIP 29724 = HD 43976 Young field object, not associated to any known moving group. Li EW and X ray emission are close to the mean locus of the Pleiades while the rotation period from TESS suggests a slightly older age. A new substellar companion is detected in this study (see Sec. 4.6 for details).
HIP 33690 = HD 53143 Star with resolved debris disk (Kalas et al. 2008). The star is flagged as a possible member of IC2391/Argus in some literature works (e.g., Nielsen et al. 2019). However, BANYAN analysis with updated Gaiaparameters return 0% membership probability for this group and no significant probability for other known MGs. Furthermore, the non detection of Lithium (Torres et al. 2000) clearly rules out the young age of IC2391/Argus. The Li non-detection and the other age indicators are fully compatible with an age similar to the Hyades. The companion candidate seen in our images and previously identified by Wahhaj et al. (2013) is a background object.
HD 57852 = HIP 35564 = HR 2813 = TYC 8132-2112-1 Bona-fide member of Carina-Near MG, confirmed with BANYAN Σ analysis using the updated kinematic parameters. The star has highly significant RV variability (Andersen & Nordstrom 1983;Pribulla et al. 2014;Desidera et al. 2015;Gaia Collaboration et al. 2018), although without orbital solution. Peak-to-valley RV range for the available data is 13 km/s. The longest-period plausible orbital solution consistent with the data yields a period of 994 d with moderately high eccentricity. Shorter periods are also possible considering the poor time sampling. The spectroscopic component is likely the responsible of the ∆µ signature and of the large RUWE in GaiaEDR3 (8.69), and it is expected to be at a separation too close to the star to be detected in SPHERE images. HD 57852 has a wide companion (HD 57853) at 9", which is itself a spectroscopic binary with three components with total mass 2.4 Msun (Saar et al. 1990;Desidera et al. 2006) The system is then a hierarchical quintuple. The source detected at 4.9" and previously identified by Chauvin et al. (2015) is a background object.
n Pup A = HD 60584 = TYC 6539-3802-1 Mid-F type star with a nearly twin wide companion (n Pup B = HD 60585 at 9.9"). Katoh et al. (2018) found the primary to an SB1 with period 366 days, eccentricity 0.49, RV semi-amplitude 12 km/s, corresponding to a minimum mass of 0.5 Msun. This orbital solution and variability above 1.6 km/s (peak-to-valley) are not supported by HARPS RV (Lagrange et al. 2009;Trifonov et al. 2020). The HARPS RV variability appears to be dominated by short-term variations (intra-night and from contiguous nights) with limited variability on longer timescales. The star is an X-ray source, but considering the very limited age dependence of X-ray emission for mid-F stars we rely only on isochrone fitting for our age estimate. The identification of a faint BD candidate, compatible with the ∆µ signature, is discussed in Sect. 4.7. The source at 4.4" is a background object.
HIP 48341 = HD 85364 = HR 3899 = 6 Sex = TYC 4899-1904-1 Early type star, identified as a possible member of UMa group by Soderblom & Mayor (1993); King et al. (2003). Chupina et al. (2006) classified it as a member of corona surrounding the UMa nucleus. Kinematic analysis using BANYAN returns a field object classification. Independently on the kinematic assignment, isochrone fitting yields an age of 660±260, formally compatible with Ursa Major. The star has also IR excess indicating the presence of a debris disk (Chen et al. 2014).
HIP 52462 = HD 92945 = V419 Hya Star with spatially resolved debris disk (Golimowski et al. 2011;Marino et al. 2019). The stellar properties and the presence of companions combining imaging, RV, and astrometric signatures were recently investigated by Mesa et al. (2021).
HIP 59726 = HD 106489 Field object with age similar to the Hyades. The RV results constant within 0.5 km/s from few literature measurements spanning several decades.

HIP 61804
Young star, not associated to any known moving group. The rotation period from TESS and chromospheric activity are similar to Pleiades members.  Chauvin et al. (2015), is a background object.
HIP 71899 = HD 129501 = TYC 1483-1030-1 F8 star. A bright source is detected within IRDIS field of view and is also seen (without full astrometric solution) in Gaia Collaboration et al. (2018Collaboration et al. ( , 2021. The combination of SPHERE and Gaiadata confirms its physical association, with a mass of 0.39 Msun. Indications of the young age of the system comes from the prominent X-ray emission (intermediate