GTC/CanariCam deep mid-infrared imaging survey of northern stars within 5 pc

In this work we present the results of a direct imaging survey for brown dwarf companions around the nearest stars at the mid-infrared 10 micron range ($\lambda_{c}$=8.7$\mu$m, $\Delta\lambda$=1.1$\mu$m) using the CanariCam instrument at the 10.4 m Gran Telescopio Canarias (GTC). We imaged the 25 nearest stellar systems within 5 pc of the Sun at declinations $\delta>-25^{\circ}$ (at least half have planets from radial velocity), reaching a mean detection limit of 11.3$\pm$0.2 mag (1.5 mJy) in the Si-2 8.7$\mu$m band over a range of angular separations from 1 to 10 arcsec. This would have allowed us to uncover substellar companions at projected orbital separations between $\sim$2 and 50 au, with effective temperatures down to 600 K and masses greater than 30 $M_{Jup}$ assuming an average age of 5 Gyr and down to the deuterium-burning mass limit for objects with ages $<$1 Gyr. From the non-detection of such companions, we determined upper limits on their occurrence rate at depths and orbital separations yet unexplored by deep imaging programs. For the M dwarfs, main components of our sample, we found with a 90% confidence level that less than 20% of these low-mass stars have L and T-type brown dwarf companions with $m \gtrsim 30 M_{Jup}$ and $T_{eff} \gtrsim$ 600 K at $\sim$3.5--35 au projected orbital separations.


INTRODUCTION
Corresponding author: Bartosz Gauza b.gauza@herts.ac.uk High contrast imaging surveys of stars constitute one of the foremost methods to find and study brown dwarfs and extrasolar planets. Their results complement our knowledge on these populations drawn from the objects detected by transit, radial velocity and other techniques. Modern imaging searches extend the accessible range of parameter space of cool companions to wider sepa-rations and longer orbital periods (a > 5 au, P > 10 yr) than those presently explored by radial velocities (RVs) or transits methods (e.g., Lafrenière et al. 2007;Biller et al. 2013;Chauvin et al. 2015a;Meshkat et al. 2015;Galicher et al. 2016;Vigan et al. 2017;Stone et al. 2018; Baron et al. 2019;Vigan et al. 2021;Launhardt et al. 2020). Moreover, directly detected brown dwarfs and planets provide a unique opportunity for spectroscopic characterisation and thereby for detailed study of their fundamental physical properties, in particular their atmospheres (e.g., Gauza et al. 2015a;Skemer et al. 2016;Bonnefoy et al. 2016Bonnefoy et al. , 2018Chinchilla et al. 2021).
A large number of imaging programs have thus far focused on young stars (less than 1 Gyr old) from the solar neighbourhood (d 50 pc), e.g., Biller et al. (2013); Lannier et al. (2016); Nielsen et al. (2019); Vigan et al. (2021). Young nearby stars are ideal targets for direct imaging surveys because their potential substellar companions, at the early stages of evolution are still warm and relatively bright, favouring their detection. As a consequence of both the limitations in sensitivity to the coldest companions at field ages, and enhanced chances of detection at young ages, the most successful searches by means of confirmed detections are so far those carried on known young stars.
The first large imaging program sensitive to substellar companions that targeted the nearest stars was the one led by Oppenheimer et al. (2001). They observed a sample of 163 northern stars in 111 star systems, located within 8 pc from the Sun using the Adaptive Optics Coronograph on the Palomar 1.5 m telescope for the optical imaging and the Cassegrain Infrared Camera on the Palomar 5 m Hale Telescope for the near-IR. For about 80% of the surveyed stars, companions more massive than 40 M Jup at an age of 5 Gyr would have been detected at separations between 40 and 120 au. Among the most sensitive imaging surveys of the nearest stars are those by Dieterich et al. (2012) and Carson et al. (2011). The first group employed the NICMOS on the Hubble Space Telescope to obtain high resolution images of 255 stars within ∼10 pc, while the second used the IRAC on the Spitzer Space Telescope and observed 117 targets at distances from 1.3 to 43.8 pc.
Despite several remarkable discoveries like the methane brown dwarf Gliese 229B (Nakajima et al. 1995) or the planetary system around HR 8799 (Marois et al. 2008), gas giant planets and brown dwarf companions at orbital separations beyond a few astronomical units were found to be rare both by the surveys of young stars and the nearest stars. In contrast to thousands of discoveries from RV and transit methods, only a few tens of companions in or around the planetary mass regime were found by direct imaging surveys sensitive enough to detect them. Dieterich et al. (2012), from analysis of their sub-sample of 138 M dwarfs, calculated a multiplicity fraction of 2.3 +5.0 −0.7 % for L and T-type companions to M dwarfs at orbital separations of 10-70 au. The IRAC/Spitzer search performed by Carson et al. (2011) had the ability to detect 600-1100 K brown dwarf companions at semimajor axes 35 au and 500-600 K companions beyond 60 au. Using Monte Carlo simulations they estimated a 600-1100 K T dwarf companion fraction of < 3.4% at 35-1200 au and < 12.4% for 500-600 K companions at 60-1000 au. Due to limitations in spatial resolution, contrast and sensitivities achieved by available instruments, the orbital separations of less than 10-15 au remained largely unexplored for the presence of massive planets and brown dwarfs by the past imaging surveys.
We used the CanariCam instrument at the Gran Telescopio Canarias (GTC) to carry out deep, high spatial resolution mid-infrared imaging in the 10 micron window, targeting the nearest known stars visible from a northern site (δ > −25 • ), to search for ultra-cool brown dwarfs and massive planets. We aimed to detect companions at 1-10 arcsec separations, which translate to 2.0-50 au, or orbital periods typically longer than 10 years. Therefore, our search extends to orbital separations and periods not explored yet by previous imaging or RV surveys. Since stars in the solar vicinity are typically old, at ages over 1 Gyr, any low-mass substellar companion will have cooled down to T eff well below 1000 K. At such temperatures, the maximum flux emission shifts from the near-to mid-IR. Hence, CanariCam at GTC provided an opportunity to perform a competitive search with respect to direct imaging surveys completed at the optical or near-IR using adaptive optics and coronographic systems, by means of sensitivity to the coolest companions. In this work we present the generic results of the program. The rest of the paper is structured as follows. Section 2 sets out the observed sample of stars, Section 3 describes the observations and data processing steps, and Section 4 presents the analysis and results: relative astrometry of resolved binaries, contrast curves and sensitivities, constraints on physical parameters of detectable companions, and the upper limit estimates on the occurrence rate of companions. Section 5 contains comparison of our results to other surveys and a discussion regarding the stellar binaries and known planets hosts in the sample. Conclusions and final remarks are presented in Section 6.

THE SAMPLE
The sample of CanariCam targets consists of the nearest known stars from the northern sky, visible from the Roque de los Muchachos Observatory, that is, with declinations δ > −25 • . We used the One Hundred Nearest Star Systems list provided by the Research Consortium On Nearby Stars (RECONS; Henry et al. 1997Henry et al. , 2006Henry et al. , 2018, complemented with the astrometric data from the Gaia DR2 and EDR3 (Gaia Collaboration 2018 where available, starting from the nearest star in the Northern Hemisphere, GJ 699 (Barnard's Star) and moving to more distant ones.
In total we have observed 33 individual stars within 5 pc arranged in 25 systems, five of which are double: GJ 820 A+B, GJ 15 A+B, GJ 65 AB, GJ 725 A+B and GJ 860 AB, and two of which are triple: GJ 866 ABC and GJ 1245 ABC. We count here GJ 866 ABC as two stars, since its individual components AC were not resolved and (AC)B were marginally resolved in our observations. Additionally, two stars, Sirius and Procyon, have known white dwarf companions. The notation "A+B" signifies that the components were observed individually as separate CanariCam targets, and "AB" that both components were observed simultaneously as a single target. The sample includes one A, one F and one G type star, three K stars and 27 M dwarfs. Such distribution of spectral types implies that our statistical results will be significant only for M dwarfs. The sample is a volume limited sample complete up to 4.0 pc. We have imaged all of the 19 known stellar systems at δ > −25 • within this distance, and 6 out of 15 known systems between 4 and 5 pc observable from our site. Due to observational limitations from target brightness, substellar objects were not considered as target primaries. Hence, the Y2-type brown dwarf WISE 0855-0714 located at 2.23 ± 0.04 pc (Luhman 2014;Luhman & Esplin 2016) and the ∼500 K brown dwarf UGPS J072227.51-054031.2 at 4.12 ± 0.04 pc (Leggett et al. 2012) were not included in our sample. The remaining 9 objects between 4 and 5 pc were not observed because of the limited telescope time available for the programme.  Table 1 lists the observed stars, including compiled information on their equatorial coordinates at J2000 epoch (proper motions taken into account), spectral types, trigonometric parallax, distance and proper motions. Table 2 list their near-and mid-IR photometry. Because all known stars in the solar vicinity have large and well determined proper motions, our survey was designed to find common proper motion companions. Any additional source detected within the field of view would have been considered as a potential companion, without any criteria based on photometric colours. Its companionship could be easily verified through second epoch observations.
We have searched through the literature to gather the available information on the planets discovered around our sample stars by RVs, transits and other methods, as well as constraints on the substellar companions from other surveys or signs of RV or astrometric trends indicating a possible presence of a distant companion. Notes with selected essential information regarding each star are compiled in Table A1 in the Appendix.

OBSERVATIONS AND DATA PROCESSING
The program was carried out in queue mode observations, starting in 2012 and completed in 2015. We used the mid-infrared camera CanariCam (Telesco et al. 2008) operating at the Nasmyth-A focal station of the 10.4 m Gran Telescopio Canarias (GTC) at the Roque de los Muchachos Observatory on the island of La Palma (Spain). CanariCam was designed to reach the diffraction limit of the GTC at mid-IR wavelengths (7.5-25 µm). The instrument uses a Raytheon 320×240 Si:As detector with a pixel scale of 79.8 ± 0.2 mas, which covers a field of view of 25.6 × 19.2 arcsec on the sky. We imaged our targets in the 10 micron window, using a medium-band silicate filter centred at λ = 8.7 µm (δλ = 1.1 µm). The choice of this particular bandpass was a compromise between the instrument performance, in particular filter transmissivity, and the sky background contribution, which is significantly higher at the N broad-band and other narrow-band filters than at the Si-2 filter. Si-2 is also favoured by a better spatial resolution, since the diffraction disc is larger at the other available narrow-band filters at longer wavelengths. Observations were executed under the following restricted atmospheric conditions: spectroscopic (clear sky with possible thin cirrus) or better, i.e., photometric/clear sky transparency, precipitable water vapor (PWV) at the level of 5-12 mm and image quality of < 0. 3, corresponding to a seeing of ∼ 0. 8 in the R band.
Observations were performed with the standard chopping and nodding technique used in the mid-IR to remove the sky emission and radiative offset. Chopping consists of switching the telescope secondary mirror at a typical frequency of a few (2-5) Hz between the posi-tion of the source (on-source) and the nearby sky (offsource). This rapid movement of the secondary mirror allows subtraction of the sky background emission that is varying in time at frequencies below the chop frequency. Movement of the secondary mirror changes the optical configuration of the telescope, resulting in two different emission patterns seen by the camera and producing a spurious signal termed the radiative offset seen in the chop-differenced images. To remove the radiative offset, the telescope is moved between two nod positions to swap over on-and off-source positions.  Table 3 continued We used an ABBA nodding sequence and "on-chip" chopping and nodding, with a chop-throw and nod offset of 8 arcsec, a chopping frequency of 1.93 or 2.01 Hz and a nod settle time of about 45 s. On chip method is recommended whenever the scientific target is pointlike, since both on-source and off-source chop positions contain the signal of the target inside the detector field of view and can be aligned and combined. Individual frames of 26 and 19 ms exposures were co-added by Ca-nariCam control software to savesets of 1.6 and 6 s depending on the brightness of the source. We used an on-source integration time of 40 min in total, divided into two observing blocks (OBs) of 20 min. Each block contained three data cube files composed of a set of in- dividual images (savesets) at subsequent chopping and nodding positions. For the two observing blocks we set the instrument at two different position angles to rotate the field of view typically by 30 deg (60 and 90 deg rotations were also used in some cases), and adjusted the configuration of chop and nod position angles so as to maintain the chop/nod parallel to the longer axis of the detector. The use of two different orientations of instrument position angle was a way to initially check the reliability of potential faint sources, distinguish from bad pixels and explore the region along the horizontal axis of detector, where the cross-talk of the star in the 16 channels is more evident and the areas otherwise obscured by the negative off-source chops. A detailed observation log is presented in Table 3. CanariCam images are stored in the standard multiextension FITS files, with a structure of [320, 240, 2, M][N], where 320 and 240 are the image pixel dimensions, 2 is the number of chop positions, M of savesets and N of nod positions. The data were processed using a set of dedicated IRAF/PyRAF 1 scripts developed within our group.
As a first step, the off-source savesets, where the star is not located at the centre of the detector, were subtracted from the corresponding on-source savesets, for respective nod beam position. These chop-/sky-subtracted frames, where the star is located at the centre of the detector, were then aligned to correct for relatively small misalignments (typically of less than five pixels) with respect to the preliminary shifts computed from the chop and nod pointing offsets. Then each pair of frames corresponding to the A and B nod positions was combined, to subtract the radiative offset. The sky-subtracted frames were multiplied by −1 to recover the negative contributions of the star (off-source position of the secondary mirror).
Because the negatives in the A and B nod positions do not overlap, being at opposite sides and at 8 arcsec of the on-source central location, they were also combined to subtract the radiative offset before they were aligned. Residual detector levels constant along single columns or lines but varying across these remained in both the positive and negative chop-and nod-subtracted frames; these were background fitted (masking the target) and subtracted.
The alignment itself was applied at once, to all (positive and negative) images of consecutive repetitions of an OB, relative to a same reference image, and so that the centroids of the target in all images including the reference image were shifted to a same, integer pixel position value. Therefore, the ulterior alignment and combination, even with other epochs, were simplified to integer pixel shifts, obviating the need of reinterpolation. Before aligning, the images were copied into larger ones to avoid the trimming of outer data regions. Then the frames were average-combined per repetition using a shallow sigma upper and lower clipping to discard occasional short transients and sharp outliers.
Each combination involved masking the negative counts of the target. Horizontal patterns of cross-talk features, apparent for the brightest targets were removed. Repetitions from OBs acquired with position angles differing from the North-up East-left orientation were resampled to common orientation using the mscimage task. For the combination of the stacks of the different OBs, the repetitions were flux-scaled according to their zero-point magnitude -as measured on the target -and weighted inversely proportional to the scaled variance of their background noise and the square of the full-width half-maximum (FWHM) of the target.
A final processed image mosaic of one of the target stars -GJ 15A, with a total on-source time of 86 min is displayed in Figure 1. The counts are represented in linear scale and within ±7σ of the background level, where σ is the standard deviation of the background noise. The central region with the highest sensitivity, where the sky area of stacked frames overlap, covers a rectangle of approximately 25 × 19 (∼89×67 au at the distance of GJ 15A). A collection of the processed Ca-nariCam images of all the observed stars is presented in

Relative astrometry of resolved binaries
As part of our CanariCam observations, we resolved the components of the binary stars GJ 65AB and GJ 860AB, and the triple systems GJ 1245ABC and GJ 866(AC)B. Cut-out images of these systems are displayed in Figure 2. We measured the relative angular separations and position angles (ρ, θ) of the resolved components. Observations of GJ 65AB and GJ 866(AC)B were repeated at two separate epochs, about 1.2 and 1.0 yr apart, respectively. In case of GJ 866, the A and C pair is a spectroscopic binary (Delfosse et al. 1999;Woitas et al. 2000) at ρ ∼ 0.01 , and thus only the B component was resolved in the first epoch. In the second epoch the components got closer by their orbital motion, to a separation below the angular resolution of 0.298 achieved on the images. Sirius B was marginally detected at ρ ∼10.3 , θ ∼ 77.7 deg, but its low signal-to-noise (S/N∼4) precludes more accurate measurements.
We obtained centroids of the sources using the IRAF imcentroid task, and transformed the X and Y pixel coordinates with their respective errors into angular separations and North-to-East position angles using the Ca-nariCam pixel scale of 79.8 ± 0.2 mas (http://www.gtc. iac.es/instruments/canaricam/canaricam.php) and the orientation given by the instrument position angle in the image header. We checked that the precision of this orientation is better than 0.3 deg by inspecting the alignment of the chopping-nodding throw with the detector X axis. We did not perform any calibration observations. Therefore, in the determination of the relative astrometry we relied on the default internal calibration of the instrument position angle and pixel scale. Measured values are listed in Table 4.

Contrast curves and achieved sensitivities
Since our observations have a bright star at the centre, by aligning the individual chop-subtracted frames we could improve the resulting FWHM as compared with straight stacking of images done by default by the automatic data reduction pipeline of the instrument provided by the observatory, and to some extent compensate the lack of fast-guiding mode of the telescope. The mean FWHM of the PSFs of all images is of 0.28 (3.5 pix), with the best and the worst values of 0.23 (2.8 pix) and 0.55 (6.9 pix), respectively. The quality of our data is close to the theoretical FWHM of the diffraction-limited PSF, which for GTC is 0.23 at 8.7 µm.
We measured the detection limits on the deepest region of the final images of each target in the survey by using the ratio of the peak counts of the star to 3 times the background noise (σ). To determine the 3σ limiting magnitudes of the images we estimated the magnitudes of the sample stars in the Si-2 filter, using the JHKs photometry from the 2MASS and WISE W1, W2, W3, W4 photometry from the All-Sky and AllWISE Source Catalogs (Wright et al. 2010). In cases when the mid-IR photometry from WISE was highly affected by saturation, we used other measurements available in the literature from Spitzer/IRAC and Akari S09W and L18W bands (Murakami et al. 2007). We converted the 2MASS and WISE magnitudes into fluxes using the corresponding Vega zero points and their errors for each band from Cohen et al. (2003) and Jarrett et al. (2011).    Table 5). Right: the mean contrast curves per given interval of brightness of the stars.    Then, we fitted a power funcion f (λ) = c 1 λ c2 + c 3 to the available measurements via a least squares method and used the obtained parameters to estimate the average flux at λ = 8.7 µm. This approach does not take into account the different spectral types. However, we checked that the effect on the estimated magnitudes at this wavelength is negligible relative to overall uncertainties. The Si-2 magnitude of each star was calculated using the Vega system zero point determined for this CanariCam filter. For Sirius A (GJ 244) we used directly the well-calibrated Si-2 flux within the standard filters from Gemini/T-ReCS observations (Skemer & Close 2011). For tight binary stars in the sample resolved by CanariCam but unresolved by WISE, the magnitudes of individual components were descomposed from the integrated magntiudes and peak-to-peak flux ratios. The values of FWHM, and Si-2 3σ detection limits for individual targets are listed in Table 5.
To measure the sensitivity as a function of angular separation from the central star, we computed the back-ground noise, σ, as a function of radial separation from the star, by measuring the standard deviation in 1 pixel wide concentric annuli around it. The 3σ noise counts were converted to contrast (difference in magnitude between the primary star and the measured quantity, noise in this case) by relating to the peak pixel value of the star's PSF. Results of this method were found to be consistent with a more realistic procedure based on inclusion of artificial sources, as in the analysis of the Barnard's Star data described in Gauza et al. (2015b). Then, the sensitivity limit was calculated using the corresponding Si-2 magnitude of a given star.
The set of graphs of Figs. 3 and 4 comprises the 3σ contrast curves and detection limit curves of the survey. The left plot of Fig. 3 collects the contrast curves (∆Si-2 mag as a function of angular separation ρ) for all observed stars individually. On the right graph are plotted the average achieved contrast versus separation, for a given brightness range of the target stars. For the brightest stars with Si-2 < 2.0 mag the maximum dy- namic range, of about 10 mag, is achieved at ∼3 arcsec. For stars fainter than 4.0 mag (20 stars of the sample) the maximum contrast is reached at 1.0-1.5 arcsec. The detection limit curves for each observed star are plotted on the left panel of Figure 4. On the right panel, we plot the mean detectability curves for stars brighter and fainter than 3.0 mag, and also the best and the worst cases, which were on GJ 699 and GJ 273, respectively. For most of the observed stars (∼80%), i.e., those with Si-2 < 3.0 mag, the detectability limit of Si-2 = 11.3 ± 0.2 mag on average, was reached at ρ 1.0-1.5 arcsec separation. For eight targets (GJ 699, GJ 65AB, GJ 729, GJ 144, GJ 447, GJ 866, GJ 15A and GJ 15B), observations were performed at two epochs separated by ∼1.0-1.7 yr. In these cases, the orbital motion of companions may not be negligible. We estimate, considering a 20-40 M Jup companion in a circular, face-on orbit, that the angular shift induced by orbital motion becomes significant, i.e., exceeds the average spatial resolution of the images, for a 8.0 − 16.5 au orbits, corresponding to ρ 2. 0 − 6. 5 angular separations at d = 2.5-3.5 pc. Stacking of images of well separated epochs would result either in a smearing or point-source doubling for any potential close-in faint companion. In any case, we searched over these closest separations around that stars in the stacks of two epochs observations both combined all together and of each epoch separately. The constrast curves and detection limits reported for these targets were computed on single-epochs images stacks.

Constraints on substellar companions
We did not find any new companion to the stars imaged in this survey. All the final processed images were examined by us for the presence of faint, point-like sources directly by a visual inspection. We searched for  Figure 6. Theoretical absolute magnitudes versus mass of brown dwarfs and giant planets at Si-2 8.7 µm band obtained using the solar abundance Ames-COND models at ages of 0.1, 0.5, 1, 5 and 10 Gyr. Several corresponding isotherms are plotted with dotted lines. The solid and dashed horizontal lines mark the mean 3σ detection limit range of the survey and a 1σ dispersion: MSi2 = 13.6 ± 0.7 mag.
candidate objects both on images combined of all the available CanariCam data, as well as on images stacked separately at different observing epochs or at different orientations of instrument position angle (typically contained in two separate observing blocks of a half of the total on-source time available). We also looked for candidates in the regions masked by the negatives by using the complementing images. We did not employ any PSF subtraction method because adjusting the image display contrast suffice to efficiently inspect the immediate surroundings of the target stars even to one FWHM separation of the stellar PSF.
The only one additional faint source was detected within the field of view on the images of GJ 820A (displayed in Fig. 5). It was detected in both OBs taken with different instrument position angle (North to East sky orientation rotated by 60 deg), with an apparent magnitude of 10.4 ± 0.3 in the Si-2 band, at an angular separation and position angle of ρ = 7.20 ± 0.04 arcsec and θ = 256.2 ± 0.3 deg, respectively. We calculated its equatorial coordinates at the epoch of the Canari-Cam observation based on the precise proper motion of GJ 820A and the measured ρ and θ. Using available archival observations and catalogues of this area of the sky, we found that the source is the background star TYC 3168-590-1 (2MASS J21065820+3845411) with V T =10.737±0.066, J=10.669±0.024, H=10.466±0.016, K s =10.405±0.013 mag. It is also catalogued in Gaia EDR3 with G=11.517±0.003 mag and a parallax of 2.51±0.01 mas. Apart from this source, no other additional sources were identified. In this section we translate the detection limits of each star to constraints on the physical properties of detectable substellar companions, specifically to their masses and effective temperatures. Because the substellar objects continuously cool down as they evolve, it is not possible to determine their mass applying unique relations independent of age, such as the mass-luminosity or mass-effective temperature relations for the main-sequence stars. Therefore, in this case one needs to rely on theoretical models providing a grid of luminosities, effective temperatures and synthetic photometry for different substellar masses as a function of age. In this work, to estimate the minimum masses and temperatures of companions that would have been detected, we used the Ames-COND models (Allard et al. 2001(Allard et al. , 2012Baraffe et al. 2003) for solar metallicity. The COND models are valid up to T eff of 1300 K and extend down to 100 K. They include the formation of dust in the atmospheres, but dust grains are considered to settle below the photosphere and are not included in the photospheric opacity. To compute the synthetic magnitudes for the Si-2 8.7 µm band we used the PHOENIX Star, Brown Dwarf & Planet Simulator available online 2 . For input we used the transmission file of the Si-2 filter and obtained the isochrones for a set of five different ages: 0.1, 0.5, 1, 5, and 10 Gyr.
In Fig. 6 are plotted the synthetic Si-2 absolute magnitudes versus masses obtained using the COND models, for objects at these five ages, jointly with several isotherms in the T eff range between 300 and 1000 K. The solid and two dashed horizontal lines mark the mean detectability level reached in this program, with the 1σ dispersion taken as uncertainty: M Si2 = 13.6 ± 0.7 mag. For an age of 1 Gyr these limits would extend to objects at the deuterium-burning mass limit and, for 10 Gyr, to about 40 M Jup . Assuming an age of 5 Gyr as a typi-2 http://phoenix.ens-lyon.fr/simulator/index.faces cal age expected for stars in the solar vicinity, this sensitivity limit translates to an average minimum mass and temperature of companions that would have been detected of m = 30 ± 3 M Jup and T eff = 600 ± 40 K. The derived constraints on companions masses and temperatures around each observed star, assuming a 5 Gyr age, are listed in Table 5.
We converted each detection limit curve into mass limits using the COND evolutionary models and considering the nominal 5 Gyr age, and counted the number of stars for which a companion at a given mass and projected orbital separation would be detectable. The resulting contour map representing the overall depth of our search for the observed sample over a grid of companion masses and separations is presented in Fig. 7.
In Table A2, we show the estimations of ages for our sample from the literature and those derived using rotation periods obtained from the literature and the gyrochronology relations from Barnes (2007), Mamajek & Hillenbrand (2008), and Angus et al. (2015). Most of the stars have an estimated age in the range of 1-10 Gyr, which justifies the selection of the 5 Gyr isochrone. However, some of them are below the 1 Gyr age and, in these cases, our estimations of detection mass limits were conservative.
We excluded the presence of brown dwarf companions with masses m 40 M Jup and T eff higher than ∼750 K around 24 of the observed stellar systems in the solar vicinity at distances within 4.6 pc, at the range of angular separations from 1.0-3.0 arcsec depending on the brightness of the target star, and up to 10 arcsec. For an average distance of 3.5 pc these angual separations correspond to projected orbital separations in the range Note-M Si2 -absolute Si2 magnitude of an object corresponding to the detection limit, M min,comp , T eff,min,comp -lower limits on mass and effective temperature of detectable companions, s -range of projected physical separations at which detection limits and minimum Mcomp and T eff apply.
from 3.5-10.5 to 35 au. The non-detection of substellar companions throughout this search allowed us to determine an upper limit to the real fraction of companions at the distances and masses mentioned above. We did not assume any shape of underlying distribution of the population of companions in terms of their masses and semi-major axes, which is equivalent to a uniform, linear flat distribution. Considering that the number of objects with companions can be well represented by a Poisson distribution, the probability of having no companions is given by the formula: where λ is the Poisson parameter, k is the number of occurrences, λ = np, with n being the number of events (observed stars) and p the real fraction of companions. For a certain confidence level (γ), where P = 1−γ, the upper limit to the real frequency of substellar companions is p = −ln(P )/n. Considering this general sample of 24 stars, the frequency of such companions is below 9.6% at a confidence level of 90%. Among these 24 stars, 18 stars are M dwarfs. For them, we were sensitive to substellar companions with m 40 M Jup and T eff 750 K at 1-10 arcsec separations, which give 1.8-18 au, ∼5-50 au and 3.5-35 au projected orbital separations for the closest, furthest and average distance of stars in our sample. Excluding the 10% of the nearest and the furthest observed stars, which set the minimum lower limit and the maximum upper limit of the projected physical separations, we covered from 2.5 to 45 au in 90% of this sample, i.e. our sur-vey is complete from 2.5 to 45 au with a confidence level of 90%. Previous imaging programs that targeted the nearest stars (Oppenheimer et al. 2001;Carson et al. 2011;Dieterich et al. 2012) could already detect companions with such masses. However, they explored wider orbital separations, beyond 10-30 au. The CanariCam survey of nearby M dwarfs allowed us to study for the first time the occurrence of substellar companions with m > 40 M Jup at orbits below 10 au. With a 90% level of confidence we set an upper limit of 12.8% on their frequency. This value is consistent with constraints derived by other imaging surveys at larger orbital separations, indicating that brown dwarf companions around lowmass M-type stars of the solar vicinity are rare, both at close orbits between 2.5 and 10 au and at wider ones.
For 11 of the observed M dwarfs we were able to detect companions at 1-10 arcsec separations with m 30 M Jup and T eff 600 K, equivalent to masses and temperatures of the expected T/Y dwarf boundary. We were thus able to establish for the first time a constrain on the upper limit on the frequency of the L and T type companions around M dwarfs, at a range of projected orbital separations between 2-3.5 and 35 au. With a 90% confidence level, we found that frequency of such companions is less than 20% around M stars. Such range of masses and temperatures was explored by imaging surveys so far only in young nearby stars, at wider orbital separations, typically beyond 30-50 au. This imaging search is the first one that probes the presence of L and T companions to M stars at orbital separations around and below 10 au, down to 2.0-3.5 au.

Comparison to other surveys
We attempt to compare the determined frequency limits of stars harboring substellar companions to results of previous works. Of the numerous programs aiming to detect giant planets and brown dwarf companions, we focus mainly on large, high-contrast imaging surveys including general field surveys probing wide separations and the most recent surveys (see references in Table 6) that placed constraints on substellar companions around M dwarfs. The essential outcomes of many of these programs can be found summarized in, e.g., Table 1 of Chauvin et al. (2015b), Bowler (2016) and Vigan et al. (2021). We also consider the results from several RV programs and studies that combine the results from various techniques. Table 6 summarizes the frequency constraints and the explored intervals of masses and orbital separations or periods, determined by each of these surveys. It is not straightforward to compare the results of these surveys to one another and to our results, because the domains of probed parameter ranges are not the same and different teams make different assumptions regarding the distributions laws of the substellar mass companions. In general, all the studies agree that the occurrence rates of substellar objects are below 20%, regardless of the companions masses/separations intervals in question and masses of the primaries. The majority of surveys points to a maximum frequency of 12% or below, typically of a few percent (∼2-5%).
As for the low-mass stars, high-contrast imaging searches targeting specifically the M dwarfs, or those including them as part of the target lists, explored objects at young ages of a few tens to few hundred million years. The most sensitive ones were capable of detecting Jupiter-mass planets at separations down to 10 au, and, most recently, even to 5 au by the SHINE survey (Vigan et al. 2021). On the other side, Doppler measurements restrict the frequencies of planetary companions with minimum masses (m sin i) as small as a few tens Earth masses and are starting to explore orbital periods of up to 10 4 -10 5 days, equivalent to approximately 5-8 au Sabotta et al. 2021).
In this context, our results bridge the explored separation ranges between wide-orbit imaging constraints and those from RVs approaching the snow line. Our limits on the frequency of substellar companions are compatible with previous studies confirming that their presence around M dwarfs is rare. In comparison with recent high contrast imaging programs, our survey is more sensitive to companions of somewhat higher masses, but at more advanced ages of a few Gyr, and hence of significantly lower T eff , extending down to 600 K.

Stellar binaries in the sample
From theory it is expected that a stellar companion alters the formation processes of exoplanets within a protoplanetary disc around the host star (see e.g. a review by Thebault & Haghighipour 2015). The binary component will also introduce a parameter space of dynamical instability in which we would not expect planets or brown dwarf companions to persist on long timescales (Holman & Wiegert 1999;. A non-negligible part of stars in our targets sample has one or more stellar companions known, including five stellar binaries, two triple systems and two stars with a white dwarf companion. We measure a multiplicity frequency (which quantifies the number of multiple systems within our sample) and a companion frequency (which quantifies the total number of companions) of 36±14%

3.5-35 au 90%
Note-a) -re-analysis of archival IRAC data; b) -determined applying DUSTY and COND models, respectively; c) -assuming a hot-start (cold-start) formation scenario; d) -minimum masses, m sin i and 44±16%, respectively, which are consistent within wide uncertainties with the comprehensive determinations by Reylé et al. (2021) in the 10 pc sample and by Ballantyne et al. (2021) for the M dwarfs. For the statistical analysis, we considered the relatively close systems (s 50 au; Sirius, GJ 65, GJ 866, Procyon, GJ 860 and GJ 1245) as individual systems and those having components at wider average orbital separations (s 50 au; GJ 820AB, GJ 725AB and GJ 15AB) as two individual stars.
We add a note of caution that such aproach introduces a potential physical bias in the interpretation of our analysis. Our search probes both circumstellar and circumbinary companions depending on the separation of the binaries, which are very different science cases to one another and to a search around single stars only. We recored a null companion detection in these systems. However, any secondary component introduces a region of instability at certain range of physical separations for both S-and P-type orbits (Holman & Wiegert 1999;Desidera & Barbieri 2007), in which we would not expect an additional body to be found. Thus, an element of bias is inherent in the statistical result derived from a combined single and binary star sample.
Several works have concluded that binarity has a minimal effect on overall planet frequency (Bonavita & Desidera 2007;Bergfors et al. 2013;Piskorz et al. 2015;Southworth et al. 2020). On the other hand, some recent observational studies demonstrate an excess of wide stellar companions to stars which host high-mass hot Jupiters and brown dwarf companions on shortperiod orbits (Ngo et al. 2016;Fontanive et al. 2019;Moe & Kratter 2019;Fontanive & Bardalez Gagliuffi 2021). These conclude that certain types of binaries may support the proccess of formation of close-in highmass planetary and substellar companions. Although early research on this front have been carried, e.g. by the SPOTS (Bonavita et al. 2016;Asensio-Torres et al. 2018) and the VIBES (Hagelberg et al. 2020) surveys, a more detailed analysis of the effect multiplicity has on the occurrence of substellar companions in general requires a larger, dedicated sample that is complete to both single stars and binaries.

Known planets hosts
As many as 12 of the 33 stars in the sample have at least one small planet of a few to a few tens Earth masses on a close orbit at around 1 au or less. All of these cases, except the G8.5V-type τ Cet, are M dwarf stars. In contrast, only two of these stars were found to host more massive planets -K2. There are no theoretical premises indicating that such close-in planets will preclude the formation of more massive companions in wider orbits. Instead, observational studies showed that hot Jupiters host stars tend to have far-away companions (e.g., Knutson et al. 2014;Lodieu et al. 2014;Wang et al. 2015). The dynamical interactions between multiple planets or a distant brown dwarf companion may affect the final orbital configuration of the system. A variety of mechanisms, such as planet-planet scattering, the Kozai-Lidov effect or secular gravitational interactions (Hansen & Murray 2015;Naoz 2016), have been invoked for smaller planets to explain inwards migration to short orbital periods (Masset & Papaloizou 2003).

CONCLUSIONS AND FINAL REMARKS
We completed a deep, high spatial resolution imaging search of substellar companions around the nearest northern stars in the mid-IR at 8.7 µm using the Ca-nariCam instrument at the 10.4 m GTC telescope. Our target sample included 25 stellar systems composed of 33 individual stars, with declinations δ > −25 • within 5 pc of the Sun. No previously undetected companions were identified in our survey.
We explored the angular separations between 1-3 and 10 arcsec with sensitivities sufficient to detect companions with masses and temperatures higher than 40 M Jup and 750 K for 24 of the observed stars, and as low as 30 M Jup and 600 K for 11 M-type stars. Considering an average distance of 3.5 pc of our sample, 3.5-35 au projected orbital separations were probed for faint companions. The non-detections enabled us to determine upper limits for the fraction of substellar companions. At a 90% confidence level, we found that less than 9.6% of the nearby stars have companions with m 40 M Jup and T eff 750 K, and that less than 20% of the closest M dwarfs have L and T companions with m > 30 M Jup and T eff 600 K within the range of explored physical separations. This is one the first imaging programs capable to detect mature substellar companions in this range of masses and temperatures below 10 au separations and provides evidence that substellar companions to lowmass M dwarfs are rare also at such closer orbits. Concurrently, extending the constraints beyond 5 au, our results are complementary to the evidences brought by RV programmes (e.g., Bonfils et al. 2013;Winn & Fabrycky 2015;Morales et al. 2019;Sabotta et al. 2021), which find occurrence of giant planets around M stars to be less than 3% at orbital periods up to 25-30 yr (∼5 au).
This work demonstrates that the modern groundbased mid-IR imaging instruments operating on 10 m class telescopes can reach angular resolutions and sensitivity limits as good as and, in certain cases (e.g. nearby, relatively old stars), better than adaptive optics systems in the optical or near-IR or space telescopes. This technique presents a high potential for direct imaging detection and studies of brown dwarfs and exoplanets. Our survey was not extensive enough to determine more precisely the true fraction of L and T brown dwarf companions at close orbits. Nonetheless results on the observed sample provide valuable constraints for next generation facilities, such as the James Webb Space Telescope or the Extremely Large Telescope, which will allow for detection and accurate characterization of the coldest companions of stars Danielski et al. 2018).
We thank the anonymous referee for a careful review of our manuscript and his/her constructive comments which substantially helped improving the quality of the paper. We are grateful to the GTC staff for performing the CanariCam observations. Based on observations made with the Gran Telescopio Canarias (GTC), installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. BG and DJP acknowledge support from the UK Science and Technology Facilities Council (STFC) via the Consolidated Grant ST/R000905/1. VJSB, MRZO and JAC acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia, Innovación y Universidades and the European Regional Development Fund through projects PID2019-109522GB-C5[1,3] We acknowledge the use of Carmencita, the CARMENES input catalogue (Caballero et al. 2016). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.  (Bonfils et al. 2013), HIRES (Butler et al. 2017) and CARMENES (Alonso-Floriano et al. 2015). Relatively active flare star (Reiners et al. 2007 (Bonfils et al. 2013); on the CARMENES GTO target list . No planets reported yet.
GJ 725 A HD 173739 Monitored using HARPS-N (Berdiñas et al. 2016). On the CARMENES GTO target list