© ESO, 2012

1. Introduction

An extensive population of exoplanets has been discovered down to sub-Jovian masses and at orbital separations below 5 astronomical units (AU) based on large-scale radial velocity (RV) surveys (Mayor & Udry 2008; Marcy et al. 2008) and on transit surveys, which are now identifying hundreds of new candidates (Cabrera et al. 2009; Borucki et al. 2011). These indirectly detected planets provide invaluable information on the distribution of close orbit planets (Cumming et al. 2008) and on their frequency around nearby stars covering a large range of masses (Marcy et al. 2008; Johnson et al. 2010a). Planets at orbital radii larger than 5 AU, comparable to the locations of the giant planets in the solar system, remain outside the range of detection of these methods, and, consequently, little is known about the extrasolar population of wide orbit planetary systems.

Direct detection with high-contrast imaging provides a method to explore the wide orbit planet population. Several adaptive optics (AO) surveys have been conducted to search for low-mass substellar companions around nearby young stars. Because direct imaging needs to overcome the large contrast ratio between the star and a potential substellar companion, the existing surveys were focused on FGKM stars, around which young objects down to a few masses of Jupiter (M_Jup) would be detectable at separation larger than a few tens of AUs. The majority of these surveys reported no (Masciadri et al. 2005; Biller et al. 2007; Lafrenière et al. 2007a; Kasper et al. 2007; Leconte et al. 2010; Ehrenreich et al. 2010; Janson et al. 2011; Delorme et al. 2012) or few (Chauvin et al. 2010) substellar companions to nearby stars and have placed upper limits on the population of massive planets at large orbital separation. Some objects were nonetheless discovered (e.g. Lafrenière et al. 2008; Thalmann et al. 2009; Biller et al. 2010), demonstrating that these objects exist, but are indeed rare (Nielsen & Close 2010).

Recent breakthrough detections around young A-stars – HR 8799 (Marois et al. 2008b, 2010b), β Pictoris (Lagrange et al. 2009, 2010) and Fomalhaut (Kalas et al. 2008; Janson et al. 2012) – have provided new perspectives on the search for companions orbiting more massive stars. The results of RV surveys of the evolved counterparts of A-stars (Johnson et al. 2007) have also provided intriguing results (Johnson et al. 2010b; Bowler et al. 2010a; Johnson et al. 2011), with significant differences in the planet population compared to RV-detected exoplanets around lower mass stars. Statistical analysis of the results clearly shows a higher frequency of planets around more massive stars and larger masses for the detected planets (Johnson et al. 2010a). Moreover, theoretical work on planet formation by core-accretion similarly shows the same kind of correlation between the star and planet mass (Kennedy & Kenyon 2008; Alibert et al. 2011).

Although more technically challenging in terms of observations, these recent results suggest that more massive stars may present more favorable targets in terms of exoplanet detections. And estimations of planet yield show that A-stars will be even more favorable for future large direct imaging surveys (Crepp & Johnson 2011) with dedicated upcoming instruments, such as Gemini/GPI (Macintosh et al. 2008) and VLT/SPHERE (Beuzit et al. 2008). To pursue the possibility of a higher frequency of massive planets around early-type stars, we have constructed a sample of 42 young A−F stars observed at high contrast to be sensitive to massive planets. We report new measurements for 39 stars, and we include 3 A-stars from the literature that have been observed in previous surveys. This survey intends to start defining the population statistics of massive planets and brown dwarfs (BDs) at orbital separations in the tens of AUs from their parent A-stars. This sample is an initial subset of the International Deep Planet Survey (IDPS) spanning M- to B-stars, the results of which will be presented in a forthcoming publication (Galicher et al., in prep.).

The selection of the target sample is explained in Sect. 2. In Sect. 3 the observing strategy, observations and data reduction steps are detailed. The detection limits of the survey are derived in Sect. 4, and we describe the identification of the candidate companions detected in the data. In Sect. 5 we use the detection limits to perform a statistical analysis of the survey, from which we derive a first estimation of the planetary systems frequency at large orbital radii around A-stars. Finally, we discuss our results and conclusions in Sect. 6.

2. Sample selection

The sample includes all A-stars previously observed in high contrast AO imaging surveys sensitive to planets (Nielsen & Close 2010; Chauvin et al. 2010; Janson et al. 2011) combined with the 35 A-stars and 4 F-stars observed for this study, resulting in a total of 42 stars, with 38 A-stars and 4 F-stars. Of the nine A-stars included in previous AO surveys, six were part of our selection (see below); they were re-observed and the new observations provided deeper limits. We used the existing surveys and the new observations to define the target set listed in Table 1, and used the best limit for each target. The sample A-stars span the full A0 − A9 spectral type range. Because Fomalhaut was not observed as part of an AO planet search survey, Fomalhaut b (Kalas et al. 2008) was not included in the sample or the subsequent statistical analysis.

In order to define a representative sample of young A-stars, the A-stars newly observed for this study were selected to have positions on the color − magnitude diagram within or below the band associated with A-stars in the Pleiades. For the A-stars in the sample, the positions of the stars on the color − magnitude diagram are shown in Fig. 1. The target ages directly impact the minimum detectable mass, given the monotonic decline in planet brightness with age (e.g. Fortney et al. 2008; Chabrier et al. 2000). The majority of the A-stars – 22 targets – have age estimates from the literature combining several techniques (Rhee et al. 2007; Moór et al. 2011; Zuckerman et al. 2011; Janson et al. 2011; Moór et al. 2006; Chauvin et al. 2010; Zuckerman et al. 2001; Tetzlaff et al. 2011; Stauffer et al. 1995). These values are reported in Table 1 and used in the analysis of mass detection limits. Of the 22 A-stars with literature ages, 10 are members of nearby moving groups or associations of young stars. For the remaining 17 A-stars, the age is based on a comparison of the position of the targets and the best-fit to the dereddened Pleiades A-stars, a population of  ~125 Myr (Stauffer et al. 1998) A-stars. The individual reddening values were taken from Breger (1986), and Pleiades A-stars on the photometric binary sequence were excluded from the fit. As shown in Fig. 1, the targets without literature ages are located generally very close to the Pleiades fit and within the range of Pleiades members, so an age of 125 Myr is adopted for these 17 stars. To ensure accurate placement on the color − magnitude diagram the A-stars were limited to stars with Hipparcos parallax uncertainty of  <5% (Perryman et al. 1997). A distance cutoff was also imposed to probe angular separations corresponding to orbital separations similar to the outer regions of the solar system and within the radii of protoplanetary disks imaged in scattered light and mm emission (e.g. Fukagawa et al. 2004; Mannings & Sargent 1997) around Herbig Ae stars, the precursors to massive stars. Both debris disk systems and stars without detected excess emission were included. Of the 38 A-stars in the sample, 14 have excess emission from debris disks, as noted in Table 1.

In addition to the A-stars, 4 early F-stars with F0 − F5 spectral types were observed. Since F-stars do not evolve off the Main Sequence as rapidly as A-stars, the color − magnitude diagram is not as a reliable tool for age assessments of young F-stars. Consequently, the F-stars were limited to stars with age estimates from the literature based on space motion linked to nearby moving groups or associations (Moór et al. 2011; Zuckerman et al. 2001; Rhee et al. 2007). Of the 4 F-stars, 2 are debris disk systems. A graphical summary of the sample population incorporating spectral type, age, and distance is given in Fig. 2. In summary, the sample consists of a total of 42 stars, with 38 A-stars and 4 F-stars with a median age of 100 Myr and a median distance of 50 pc; key parameters of the sample are listed in Table 1. All members of the sample have no catalogued visual binary within 5′′ to ensure a detection limit that is uniform with position angle. The targets with wide (>5′′) common proper motion companions from the WDS catalogue (Mason et al. 2001) or close (<0.1′′) speckle companions (Mason et al. 1999) are noted in Table 1.

Fig. 1 MV vs. B − V color − magnitude diagram of the A-stars of our target sample. The stars with published ages are plotted with a green circle. Stars without a published age are noted with blue triangles and have been assigned an age of 125 Myr, based on their positions within the region defined by the Pleiades A-stars, which are indicated with orange stars. An empirical linear fit to the Pleiades is also shown. The older population of 625 Myr Hyades A-stars is plotted with blue squares.

Fig. 2 Age, distance and spectral type of all the stars in our sample. Table 1 gives the individual properties of each target. The size and color of the symbols correspond to the spectral type of each target: bluer colors for earlier spectral types and redder colors for later spectral types. The points grouped at 125 Myr correspond to stars without a published age, as described in Sect. 2.

3. Observations and data reduction

Data were acquired between 2007 and 2012 with VLT/NaCo and Gemini/NIRI, two instruments located at the focus of 8 m telescopes. NaCo is installed at the ESO Very Large Telescope in Chile. It is the combination of the NAOS adaptive optics system (Rousset et al. 2003) and the CONICA observing camera (Lenzen et al. 2003). It provides diffraction-limited images in the near-infrared on a 1024 × 1024 Aladdin InSb detector array. Our data were acquired using the S13 camera, which provides a spatial sampling of  ~13 mas/pixel (better than Nyquist) and a field of view (FoV) of . The targets were all observed with the broadband K_s filter (2.00−2.35 μm).

NIRI (Hodapp et al. 2003) is installed at the Gemini North telescope in Hawaii. It was used in combination with the ALTAIR AO system (Herriot et al. 2000) to obtain diffraction-limited images on a 1024 × 1024 Aladdin InSb detector array. The f/32 camera provides a sampling of ~22 mas/pixel and a FoV of . To improve the Strehl ratio and provide sharper image quality at large distance from the star, ALTAIR was used with the field lens inside the optical path. The targets were observed either with the CH4-short filter (CH4s, 1.54 − 1.65 μm) or the K′ filter (1.95 − 2.30 μm).

Both instruments were used in pupil-stabilized mode (pupil-tracking mode on NaCo) to reduce the number of independently rotating speckle patterns, improve the stability of the point-spread function (PSF), and provide rotation of the FoV, which is then used for proper implementation of angular differential imaging (ADI; Marois et al. 2006). To optimize the detection of faint companions, the science images were saturated by using detector integration times (DIT) between 15 and 30 s in most cases. These DITs allow the detections to be background-limited at large separation, but since our targets are very bright (H ≤ 7.0), the regions up to 0.4 − 0.8′′ are in general completely saturated and have been masked in the analysis. Each observing sequence started with the acquisition of an unsaturated PSF at five dither positions, obtained with the ND_short neutral density in NaCo (transmission of 1.20 ± 0.04%) or a short DIT combined with narrow-band filters in NIRI (H-continuum in H-band or H₂ in K-band). These PSFs were then followed by 100 − 250 deeply saturated science images. For some observations, sky frames were acquired prior to the science frames following a dithering pattern. Due to restrictions on detector saturation, second epoch data on NaCo were acquired with a much smaller integration time and a larger number of exposures. Some NaCo data was also acquired with a much smaller number of exposures but using co-addition of 5 exposures. All the observations are summarized in Table 2.

The data reduction followed standard procedures. For the unsaturated dithered PSF frames, a sky frame was constructed by taking the median of the images. The individual images were then sky subtracted and divided by the flat field before being registered to a common center. The final unsaturated PSF was obtained by median combining the five images. This PSF was used for the photometric calibration and calculation of the detection limits. For the saturated images, a sky frame was subtracted when available, and each image was divided by the flat field. Bad pixels identified from the flat field were replaced by the median of adjacent pixels using a custom IDL routine. Field distortion was taken into account for NIRI since it reaches several pixels in the outer part of the field: all images were corrected for this effect using a custom IDL routine and an optimized distortion solution that will be presented in Galicher et al. (in prep.). For NaCo, the field distortion was found to be negligible (<1 mas) with the S13 camera in the inner part of the detector (Trippe et al. 2008; Fritz et al. 2010). The center of each saturated image was estimated using the fitting of a Moffat profile (Moffat 1969) on the PSF wings. The saturated region of the images was assigned a zero-weight in the fit to avoid any bias. For NIRI, the center of the saturated images moved by up to 5 pixel on the detector over the full observing sequence. For NaCo, a well-documented problem with the pupil-tracking mode (corrected in October 2011) induced a slow drift of the star during some of the observations. This drift, which could reach several pixels per minute, is proportional to the star elevation. To avoid any important drift, the observing strategy included frequent recentering of the star during the observations and the use short DITs (during the follow-up phase of the survey).

In each observing sequence, a frame selection was performed to reject the worst frames where a sudden increase of seeing or an interruption of the AO loop occurred. The frame selection was based on three parameters: the position of the PSF center in each frame compared to previous and subsequent frames, the overall flux in the frame and a measure of the fraction of flux contained in an annulus covering an area just outside of the saturated region. Frames for which a sudden change of the PSF position or a drop in intensity are observed generally correspond to frames where the AO loop is open, and are thus of poor quality. These frames were inspected visually and in most cases removed from the observing sequence. After the frame selection process, the sequence of images usually contained frames of equivalent quality in terms of AO correction.

For astrometric calibration of the NIRI data, we used archive data for the Θ₁ Ori C field observed on November 30, 2008 with the f/32 camera, K′ filter and without the ALTAIR field lens. We used the coordinates reported by McCaughrean & Stauffer (1994) to determine the mean pixel scale and true North orientation of the NIRI detector. The mean pixel scale and true North orientation measured were respectively 21.38 ± 0.14 mas/pixel and 0.27 ± 0.05 deg. For NaCo observations, we used calibration data collected within our observing programs and within the NaCo Large Program collaboration (NaCo-LP¹; Chauvin 2010). The calibrators used were the Θ₁ Ori C field and IDS 13022+0107, when the former was not observable. The calibration was performed with the S13 camera and the H filter at several epochs between November 2009 and January 2012. The mean pixel scale measured was 13.19 ± 0.07 mas/pixel, and the true North orientation varied from −0.17 ± 0.05 deg to −0.71 ± 0.05 deg between the different epochs. The NaCo PSFs are usually not perfectly centro-symmetric, causing the centers determined on unsaturated and saturated PSFs to be slightly different. The effect has been estimated to  ~6 mas and constitutes the dominant astrometric error. It was estimated using a dedicated calibration where saturated and non-saturated PSFs were taken alternatively, and their center determined using Moffat or Gaussian profile fitting.

Finally, the cleaned registered ADI images were processed with an implementation of the LOCI algorithm (Lafrenière et al. 2007b). All the targets were analyzed using a similar set of parameters: N_δ = 0.75 FWHM, N_A = 300 for NaCo and N_A = 500 for NIRI. For each target, partial subtraction induced by LOCI was measured using fake planets introduced into the data at several position angles and angular separations, in a way similar to that described by Lafrenière et al. (2007b). Temporal smearing induced by the FoV rotation was also measured using fake planets. A typical final reduced image is visible in Fig. 3.

Fig. 3 Example of final image after data reduction and processing with LOCI for target HIP 10670. The central white region corresponds to the saturated area in the science images, which has been masked in the data analysis. The point source 4.7′′ West of the star is a background object 14.6 mag fainter than the central star.

4. Results

4.1. Identification of companion candidates

Detection of the candidate companions (CCs) was performed by visual inspection of the final images after LOCI processing and of the signal-to-noise ratio (SNR) maps. The SNR maps were obtained by normalizing each pixel in the images to the noise measured in an annulus of width λ/D at the same separation as the pixel from the star. The maps were used to identify CC with a detection level above 5σ.

Relative astrometry for each CC was obtained using 2D Gaussian fitting, and the position was refined by maximizing the flux in a 1.5λ/D aperture moved on a 2 × 2 pixel grid by steps of 0.1 pixel around the position found with the Gaussian fitting. The main astrometric error for the NaCo data is the 0.5 pixel uncertainty on the actual center of the PSF in saturated data (see Sect. 3). Other sources of error come from the shift induced by the use of the neutral density (0.1 pixel) or the 2D Gaussian fitting (usually 0.1 pixel). In NIRI data, we assume the dominant error term to be the residual distortion (6 mas) close to the star. At separations larger than 5′′, the distortion correction improves significantly the astrometric measurements for the CCs but not enough to reach a similar accuracy. For these candidates we have assumed a residual error of 1 pixel (21.4 mas). Although there are some uncertainties induced by the distortion in the NIRI field of view, the precision was sufficient to clearly identify background or comoving objects. An updated distortion map for NIRI will be presented in Galicher et al. (in prep.).

Relative photometry was estimated by aperture photometry in 1.5λ/D aperture centered on the best estimation of the position. Contrast difference with respect to the star was obtained by dividing the flux of each CC by the flux of the non-saturated PSF measured in an aperture of the same diameter, scaled by the ratio of exposure time and the transmission of the neutral density (NaCo data) or the correction due to the fact that unsaturated observations were not performed in the same filters as saturated observations (NIRI data). This correction factor for NIRI was calculated using slightly saturated images acquired regularly during the ADI sequence with the same filter as the deeply saturated science images. The non-saturated PSF was normalized to these slightly saturated images using the integrated flux in an annulus outside of the saturated part of the images (usually between 10 and 20 pixels of radius). A calibration performed on data saturated and unsaturated with the same filters indicated an error on the flux normalization of  ~4%. The flux loss induced by LOCI partial subtraction and by the temporal smearing of the CC PSF by field rotation was also included as a correction term. Finally, the photometric error was estimated by including a measure of the local SNR at the separation of each CC and the error on the calibration of the flux normalization between the saturated and unsaturated PSF.

For newly detected stellar companions (see Sect. 4.3.1), where a more robust estimation of the relative astrometry and photometry is useful for the determination of their physical parameters, a different approach was taken: instead of estimating the position and signal level of the CC in the final reduced image, we implemented a routine to subtract a “negative CC” in each raw frames before performing the ADI optimization in a way similar to that described in Marois et al. (2010a). The routine was run multiple times for various insertion position of the negative CC in order to minimize the residual noise in the final image at the expected position of the CC and in the surrounding area. When a position minimizing the residuals was found, the routine was used again for different injection flux for the CC, again trying to minimize the residual noise in the final image. Errors were estimated using variations of the best injection point and flux obtained by measuring the residuals in areas of different sizes around the position of the CC.

A total of ~50 CCs were identified above the 5σ level around 22 of the observed targets. Targets with CCs at a projected physical separation less than or equal to 320 AU were followed-up by subsequent observations to verify if they are co-moving with the target star. This cut-off was chosen according to the widest likely known companion at the time of the follow-up, 1RXS J160929.1-210524 b, which is at a projected separation of  ~330 AU (Lafrenière et al. 2008). There is only one target for which no second epoch was obtained although there was one CC within 320 AU: HIP 41307 had a CC at a projected separation of 217 AU but it was confirmed as a background star by Janson et al. (2011). For targets with follow-up observations, the position of the candidates was checked for consistency with either a background or a co-moving object based on the target parallax and proper motion, also taking into account uncertainties on these particular values. Similarly to Chauvin et al. (2010) and Ehrenreich et al. (2010), we performed a χ² probability test on the difference in position Δα and Δδ with respect to the star at two epochs. These probability tests along with the proper motion plots in RA/DEC or in separation/position angle for each CC allowed attributing a status of background, comoving or ambiguous to each of them. The properties of all detected CCs are summarised in Appendix A.

4.2. Detection limits

The ultimate contrast reached by our observations is quantified using 5σ detection limits obtained by measuring at each angular separation the residual noise in an annulus of width λ/D, normalized by the flux of the unsaturated PSF. As mentioned previously, two correction factors were applied to the detection limits. The first one accounts for the flux loss of point sources induced by the LOCI processing. It is calculated by introducing fake companions of known flux at various separation and position angle into the raw data, fully processing the data with LOCI, and measuring the flux of the companions in the final image (see Lafrenière et al. 2007b). This process was repeated 10 times with the companions at various azimuthal positions to obtain the average attenuation as a function of angular separation from the star. The second correction factor takes into account the smearing of the PSF of an object in the field induced by the field rotation during the exposures. The dilution factor was calculated by simulating the smearing of point sources at various separations for every frame. Finally, the detection limits were corrected by dividing them at all separations by the appropriate attenuation for the flux loss and dilution for the smearing.

The best magnitude difference detection limits for angular separations ranging from 0.5′′ to 10.0′′ are listed in Table 3. Example limits representative of good, median and poor contrast performance in K-band are plotted in Fig. 4. The observations are background-limited outside  ~5′′ for all targets. The background limits range between 19 and 23 mag for NIRI in H-band, between 17.5 and 22.5 mag for NIRI in K-band, and between 16.5 and 20.5 mag for NaCo in K-band. These values are within the typical ranges for observations in good conditions with NIRI and NaCo in 0.5 to 1 h of integration time.

A common concern for the estimation of detection limits in high-contrast data is that of the detection threshold. It is known that speckle noise does not follow Gaussian statistics (Goodman 1968; Fitzgerald & Graham 2006). However, Marois et al. (2008a) have performed a detailed analysis of NIRI data processed with ADI (Marois et al. 2006) and demonstrated that the residual noise follows quasi-Gaussian statistics when a sufficient number of images (i.e. of FoV rotation) are combined. For most targets², we obtained a FoV rotation of at least 20 deg (Table 2) over several tens of images, which should be sufficient to reach quasi-Gaussian noise statistics after LOCI processing. For comparison, previous surveys using some form of ADI observations have adopted 5σ or 6σ thresholds to derive their detection limits, either assuming Gaussian statistics for the residual noise (Biller et al. 2007; Lafrenière et al. 2007a; Leconte et al. 2010; Janson et al. 2011) or using a higher threshold to account for the non-Gaussianity of the noise (Kasper et al. 2007).

The 5σ detection limits were finally converted to absolute magnitudes using the Hipparcos parallax and 2MASS H or K_s apparent magnitudes. Since observations were obtained in a filter set different from the 2MASS filters (in particular the CH4s NIRI filter), we checked that for the different spectral types of our targets (from A0V to F5V) there was no significant change in apparent magnitude between the 2MASS filters and the NaCo or NIRI filters. The correction terms were found to be completely negligible (<0.03 mag) compared to the other systematic uncertainties in the photometry.

4.3. Companion detections

Within the 42-star sample, 5 stars have newly identified or previously reported physically associated companions, and we briefly discuss each case of a stellar or substellar companion in this section. Of the companion candidates with follow-up newly reported in this study, three are confirmed as co-moving, and all but one of the remaining candidates have a  <0.01% probability to be comoving. The possible residuals on astrometry are well below the threshold that would change the assessment from background to co-moving. In one case, one of candidate companions to HIP 10670, it cannot be excluded that this object has a proper motion of its own, making it inconsistent with both a stationary background object and a co-moving object. Because it shows a motion much larger than a pixel (13.19 mas on NaCo, 21.38 mas on NIRI) between the two epochs, suggesting a significant motion not compatible with a co-moving object, it is treated as a background object in the statistics.

Fig. 4 Examples of K-band 5σ detection limits representing the best (HIP 69732) and worst (HIP 15648) limits, as well as an intermediate case (HIP 26624). The right axis shows the detection limit in masses of Jupiter for the median target of our survey, a 100 Myr-old A3V star at 50 pc, converted using the COND evolutionary models (Baraffe et al. 2003).

4.3.1. Stellar companions

We have identified three co-moving stellar companions around two stars in our data. Both stars were already known to be in wide multiple systems, but the companions we report were previously unknown due to their large contrast ratio with their parent star. Proper motion plots are given in Fig. 5, and relative positions and fluxes are summarized in Table 4.

HIP 42334 – The star is a 125 Myr old A0V star at 71.07 pc in a double system, with the secondary component at a separation of 20′′ and a magnitude of 13 in the optical (Mason et al. 2001). A point source 8.2 ± 0.1 mag fainter in K_s than the primary was detected at a separation of and position angle (PA) of 11.52° ± 0.12°. The star was re-observed because of the presence of other fainter point sources in the field, but only the brightest candidate proved to be co-moving (Fig. 5, left). Comparison to evolutionary models at 125 Myr (Chabrier et al. 2000) place this object at the bottom of the main sequence, consistent with an M6 − M8 dwarf.

HIP 104365 – The star is 125 Myr old A0V star at 55.12 pc in a quadruple system (Mason et al. 2001). A very tight pair of companions was detected at a separation of  ~9.57′′ and a PA of  ~14.8°. After one year they both clearly showed a common proper motion with HIP 104365 (Fig. 5, right). Given their projected separation of 9.8 AU and their common proper motion, it is extremely likely that they are gravitationally bound together. Their contrast difference with the star in the CH4s filter makes them compatible with M5−M7 dwarfs.

4.3.2. Substellar companions

Three targets in the sample – HR 8799, β Pic, and HR 7329 – have previously reported substellar companions (Marois et al. 2008b, 2010b; Lagrange et al. 2009; Lowrance et al. 2000; Neuhäuser et al. 2011) including four planets and one brown dwarf. A summary of these objects and their properties is given in Table 5. We note that the observations of these three targets were not unusually sensitive and biased, but have detection limits similar to many other targets in the sample.

HR 8799 – This is the first imaged multiple planet system, with four known planetary-mass companions orbiting the star at 68, 38, 24 and 14.5 AU (Marois et al. 2008b, 2010b). It is an A5V star located at 39.4 ± 1.0 pc from the Sun, classified as γ Doradus and λ Bootis. It shows a far-infrared excess emission (Vega-like) consistent with the presence of a warm dust disk within 6 − 15 AU, a massive cold dust disk within 90 − 300 AU and a dust halo extending up to 1000 AU (Rhee et al. 2007; Su et al. 2009). Recent estimations of the system age converge towards 30 Myr (Doyon et al. 2010; Zuckerman et al. 2011; Currie et al. 2011), with HR 8799 being a likely member of the Columba Association. This estimation leads to a mass of 5 − 10 M_Jup for HR 8799 cde, and 3.5 − 7 M_Jup for HR 8799 b. This system has been extensively studied in the literature since its discovery and we do not intend to make an exhaustive review of the numerous results. For detailed studies, see for example Bowler et al. (2010b), Currie et al. (2011), Barman et al. (2011a), Galicher et al. (2011) and Soummer et al. (2011).

Fig. 5 Offset position (black circles) with uncertainties measured relative to the primary star for the co-moving stellar companions to HIP 42334 (left) and HIP 104365 (right). The variation expected for a background stationary object is also plotted (grey line) with the expected positions of background objects (grey square).

β Pictoris – This $\mathrm{1}{\mathrm{2}}_{-4}^{\mathrm{+}\mathrm{8}}$ Myr-old (Zuckerman et al. 2001) A5V star at a distance of 19.4 ± 0.2 pc (Crifo et al. 1997) has been well-known for two decades since the imaging of a debris disk (Smith & Terrile 1984) with several asymmetries and a warp at  ~80 AU (Kalas & Jewitt 1995; Mouillet et al. 1997), which could be explained by the presence of a giant planet within 50 AU of the star. Lagrange et al. (2009, 2010) confirmed the presence of a planetary-mass companion with an orbital separation compatible with formation by core-accretion. Using the COND evolutionary models and K-band observations, Bonnefoy et al. (2011) infered a mass range of 7 − 10 M_Jup for the planet. Follow-up of the system with imaging (Lagrange et al. 2012a; Chauvin et al. 2012) and RV measurements (Lagrange et al. 2012b) recently allowed an accurate determination of the orbital parameters and the mass of the planet. Using Markov Chain Monte-Carlo methods, Chauvin et al. (2012) measure a probable range of a = 8 − 15 AU for the semimajor axis, with a low eccentricity e ≤ 0.16 and a high inclination i = 88.5 ± 1.7 deg. Using a dataset where the planet, the main disk and the warp are detected, Lagrange et al. (2012a) also confirmed that the planet is located above the midplane of the main disk, making it aligned with the warp. These observations strongly support that the presence of the planet is responsible for the warped morphology of the disk. Finally, using 8 years of RV measurements on the star, Lagrange et al. (2012b) determine an upper limit of 10 − 12 M_Jup for the mass of the planet, assuming a semimajor axis of 8 − 9 AU.

HR 7329 – The substellar companion to the β Pic moving group member HR 7329 (Zuckerman et al. 2001) was discovered with coronagraphic imaging with the Hubble Space Telescope (Lowrance et al. 2000). It has a magnitude difference of ΔF160W = 6.9 ± 0.1 and separation of (200 AU at 47.7 pc). It was further characterized with spectroscopic observations in the H-band (Guenther et al. 2001). The most comprehensive compilation of astrometric and photometric measurements of this system and a high significance confirmation that the companion is physically associated with the primary is given in Neuhäuser et al. (2011). Like the HR 8799 and β Pic systems, HR 7329 has a debris disk (Smith et al. 2009).

5. Statistical analysis

5.1. Context

Using AO instruments on large ground-based telescopes, a number of deep imaging surveys have been carried out to search for massive planets around young nearby solar-type stars (Chauvin et al. 2003; Neuhäuser et al. 2003; Masciadri et al. 2005; Biller et al. 2007; Kasper et al. 2007; Lafrenière et al. 2007a; Chauvin et al. 2010; Delorme et al. 2012), and more recently around much more massive stars (Ehrenreich et al. 2010; Janson et al. 2011). Despite a large improvement of sensitivity to massive planets in the recent years, due largely to the development of advanced data analysis methods to suppress the speckle noise, none of these surveys have reported the detection of planetary objects potentially formed by core-accretion, i.e. within the circumstellar disk of the primary.

Taking advantage of their null results, the most recent large-scale surveys have placed constraints on the population of long-period massive planets. Using a combination of Monte-Carlo simulations and a Bayesian approach, it is possible to estimate the upper limit of the fraction of stars having substellar companions (f_max) within a certain range of orbital separation and mass. Nielsen & Close (2010) performed a detailed statistical analysis using the data for 118 stars from the surveys of Masciadri et al. (2005), Biller et al. (2007) and Lafrenière et al. (2007a), which allowed them to constrain the presence of planets down to 4 M_Jup around F − M stars. Assuming a flat distribution for the mass and semimajor axis of planets, they rule out at 95% confidence the presence of planets more massive than 4 M_Jup around more than 20% of stars in the range 80−280 AU using the cold-start evolutionary models (see Sect. 5.3 for details on the evolutionary models) and in the range 20 − 500 AU using the hot-start models, which is in agreement with the values previously reported by Lafrenière et al. (2007a). From the observations of 22 very young G − M targets in the Tucana and β Pic moving groups, Kasper et al. (2007) place an upper limit of 5% for planets above 2−3 M_Jup at separations larger than 30 AU. Using a sample of 88 young nearby stars, Chauvin et al. (2010) report with 95% confidence that less than 10% of stars can have planetary-mass companions beyond 40 AU, assuming mass and semimajor axis distributions similar to that reported by Cumming et al. (2008) from RV surveys. Finally, Janson et al. (2011) compared predictions of disk instability models to the observation of 15 B2 − A0 stars, and they showed that less than 30% of massive stars can form and retain companions below 100 M_Jup within 300 AU at 99% confidence.

Previous surveys have also used their non-detections to investigate with what confidence planet distributions obtained from RV surveys can be extrapolated for planets at much larger orbital radii. This is in particular the case for the mass and semimajor axis distributions, which appear to be well fitted by simple power laws (Cumming et al. 2008) up to a certain semimajor axis cutoff. For the power law and planet frequency values reported by Cumming et al. (2008), Nielsen & Close (2010) place the semimajor axis cutoff within 180 AU at 95% confidence with the cold-start models, and within 65 AU with the hot-start models. Similarly, Kasper et al. (2007) rule out the possibility of a positive power law index for a semimajor axis distribution with a cutoff beyond 30 AU. With their data, Chauvin et al. (2010) explored a wide range of parameters for the power law indices and concluded that the non-detection probability is very sensitive to the semimajor power law index, which means that strong constraints on this parameter can be inferred from observations. For the combinations of parameters reported by Cumming et al. (2008), Chauvin et al. (2010) report a value of the semimajor axis cutoff within  ~100 AU with 95% confidence.

Although the reported values are not identical due to the different samples and sensitivities, all previous surveys converge towards a similar conclusion, which is that the frequency of planets above 2−4 M_Jup around F−M stars is below ~10−20% beyond  ~50 AU at high-confidence level. They also show that the distribution of semimajor axis is consistent with a power law of negative index, i.e. with a decrease of the number of planets with orbital separation.

The present work is similar in its aims to previous surveys, but with two major exceptions: it is focused on A-stars, and for the first time it includes the detections of the planetary system around HR 8799 and the planet around β Pic. A-stars are particularly interesting because RV surveys show a strong correlation between planet mass and stellar mass (Johnson et al. 2007, 2010a). Although this correlation may not hold for planets at large radii, it is tempting to test this hypothesis, especially given the planets recently imaged around young A-type stars (Marois et al. 2008b, 2010b; Lagrange et al. 2009, 2010). It is therefore essential that our statistical analysis take into account the confirmed detections for the estimation of the planet frequency around A-stars.

In the following sections, we present our statistical formalism and several applications of our data to place constraints on the frequency of giant exoplanetary systems at large orbital radii around A-stars, and their possible mass and semimajor axis distributions. We perform the analysis for a total of 42 stars: the 39 that we have observed with NaCo and NIRI, and 3 young A-stars included from the survey of young, nearby austral stars by Chauvin et al. (2010) – β Pic (using the data presented in Bonnefoy et al. 2011), HIP 61498 and HIP 98495.

5.2. Statistical formalism

Our formalism for the statistical analysis is based on previous works by Carson et al. (2006) and Lafrenière et al. (2007a). We consider the observation of N stars enumerated by j = 1...N. We note f the fraction of stars that have at least one companion of mass and semimajor axis in the interval [m_min,m_max] ∩ [a_min,a_max], and p_j the probability that such a companion would be detected from our observations. The analysis does not differentiate between a single planet system and a multiple planet system, i.e. f represents the frequency of giant planetary systems at large orbital radii. With these notations, the probability of detecting such a companion around star j is (fp_j) and the probability of not detecting it is (1 − fp_j). Denoting  { d_j }  the detections made by the observations, such that d_j equals 1 if at least one companion is detected around star j and 0 otherwise, the likelihood of the data given f is $\mathit{L}\left(\mathrm{\left\{}{\mathit{d}}_{\mathit{j}}\mathrm{\right\}}|\mathit{f}\right)\mathrm{=}\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{f}{\mathit{p}}_{\mathit{j}}\mathrm{)}}^{\mathrm{1}\mathrm{-}{\mathit{d}}_{\mathit{j}}}\mathrm{·}{\mathrm{\left(}\mathit{f}{\mathit{p}}_{\mathit{j}}\mathrm{\right)}}^{{\mathit{d}}_{\mathit{j}}}\mathit{.}$(1)The determination of the probability that the fraction of stars having at least one companion is f is obtained from Bayes’ theorem: $\mathit{p}\left(\mathit{f}|\mathrm{\left\{}{\mathit{d}}_{\mathit{j}}\mathrm{\right\}}\right)\mathrm{=}\frac{\mathit{L}\left(\mathrm{\left\{}{\mathit{d}}_{\mathit{j}}\mathrm{\right\}}|\mathit{f}\right)\mathrm{·}\mathit{p}\mathrm{\left(}\mathit{f}\mathrm{\right)}}{{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\mathit{L}\left(\mathrm{\left\{}{\mathit{d}}_{\mathit{j}}\mathrm{\right\}}|\mathit{f}\right)\mathrm{·}\mathit{p}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathrm{d}\mathit{f}}\mathit{,}$(2)where p(f) is the a priori probability density of f, or prior distribution, and p(f |  { d_j } ) is the probability density of f given the observations  { d_j } , or posterior distribution. Since we have no a priori knowledge of the wide-orbit massive planets frequency, we adopt a “maximum ignorance” prior, p(f) = 1. This approach is similar to previous studies for deep imaging surveys, allowing a direct comparison of the results.

Given a confidence level (CL) α, we can use the posterior distribution p(f |  { d_j } ) to determine a confidence interval (CI) for f using $\mathit{\alpha }\mathrm{=}{\mathrm{\int }}_{{\mathit{f}}_{\min }}^{{\mathit{f}}_{\max }}\mathit{p}\left(\mathit{f}|\mathrm{\left\{}{\mathit{d}}_{\mathit{j}}\mathrm{\right\}}\right)\mathrm{d}\mathit{f,}$(3)where f_min and f_max represent the bounding values of f, i.e. the minimal and maximal fraction of stars with at least one planetary companion. Previous studies have always assumed a value of f_min = 0 since no detection was included. In that case, Eq. (3) becomes an implicit equation on f_max, which can be solved numerically to estimate f_max. In our case, we want to use our detections to infer a value for both f_min and f_max. Following Lafrenière et al. (2007a), an equal-tail CI  [f_min,f_max]  is found by solving: $\begin{array}{ccc}\frac{\mathrm{1}\mathrm{-}\mathit{\alpha }}{\mathrm{2}}& \mathrm{=}& \int \begin{array}{}\mathrm{1}\\ {\mathit{f}}_{\max }\end{array}\mathit{p}\left(\mathit{f}|\mathrm{\{}{\mathit{d}}_{\mathit{j}}\mathrm{\}}\right)\mathrm{d}\mathit{f}\\ \frac{\mathrm{1}\mathrm{-}\mathit{\alpha }}{\mathrm{2}}& \mathrm{=}& {\mathrm{\int }}_{\mathrm{0}}^{{\mathit{f}}_{\min }}\mathit{p}\left(\mathit{f}|\mathrm{\{}{\mathit{d}}_{\mathit{j}}\mathrm{\}}\right)\mathrm{d}\mathit{f}\mathit{.}\end{array}$Again, when a value of α is fixed, Eqs. (4) and (5) become implicit equations on f_min and f_max that can be solved numerically.

5.3. Monte-Carlo simulations

The critical step of the statistical analysis is the determination of p_j, the completeness, for all observed targets. This completeness represents the fraction of companions that would have been detected in the interval [m_min,m_max] ∩ [a_min,a_max] with the observations. It is directly related to the sensitivity of the observations to the physical parameters of the planets, which in turn is a function of the detection limits, age and distance of each star.

Fig. 6 The detection limit (dashed line) for wide orbit planets from the observations of the target HR 8799 in the NIRI CH4s filter (left) and K′ filter (right) compared to the DUSTY (red line) and COND (black line) evolutionary models (Chabrier et al. 2000; Baraffe et al. 2003). The DUSTY model mass-magnitude grid does not extend to magnitudes as faint as the detection limits. The current estimates of the HR 8799 planet masses (Marois et al. 2010b) are also plotted, and the empirical masses are consistent with the COND magnitudes within the error bars.

A classical approach to estimate p_j is to use Monte-Carlo (MC) simulations to generate large populations of planets with random physical and orbital parameters, and check their detectability given the sensitivity limits of the survey. For the MC simulations, we used the Multi-purpose Exoplanet Simulation System (MESS, Bonavita et al. 2012), a tool specifically designed for the statistical analysis of exoplanet surveys, which has already been used for the statistical analysis of several surveys (Chauvin et al. 2010; Janson et al. 2011; Delorme et al. 2012). The mass and semimajor axis of the planets are generated either from a linear grid or from specific distributions like power laws extrapolated from RV surveys of close-in planets (e.g. Fischer & Valenti 2005; Cumming et al. 2008; Johnson et al. 2010a). Then for each combination of mass/semimajor axis, a large number of planets are generated by sampling randomly the other orbital parameters: the inclination has a uniform distribution in sini, the longitude of the ascending node and the argument of periastron have a uniform distribution between 0 and 2π, and the time of passage at periastron has a uniform distribution between 0 and P, with P the orbital period. The eccentricity distribution is uniform between 0 ≤ e ≤ 0.8 as suggested by surveys of long period exoplanets (Cumming et al. 2008). However, given the large number of planets generated at each mass/semimajor axis point, the effect of the eccentricity distribution does not impact significantly the results. Once all the parameters are chosen for a planet, its projected separation on the sky plane is calculated taking into account the distance of the star, all the orbital parameters and the date of the observations.

The other important observable needed for comparison to the detection limits is the luminosity of the planet. This quantity is directly related to the planet mass and age, and can be estimated using evolutionary models. The predictions from evolutionary models can vary significantly depending on which model is considered. In particular, the cold-start (Marley et al. 2007; Fortney et al. 2008) models predict much fainter planets at young ages for a given mass than the hot-start models (Chabrier et al. 2000; Baraffe et al. 2003; Burrows et al. 2003). Hot-start models generally use arbitrary initial conditions, independent of the outcome of the formation phase, assuming that they will be “forgotten” after a few millions of years (Baraffe et al. 2002). In contrast, cold-start models have been developed to use initial conditions calculated from a core-accretion mechanism for the formation of the planets. However, the simplified treatment of the core-accretion shock results in significant uncertainties on the predictions of these models (Fortney et al. 2008). Existing observations are more consistent with predictions from the hot-start models rather than the cold-start models (see Janson et al. 2011, Fig. 4), although the much fainter magnitudes of cold-start planets would not have been detectable by most observations. Only hot-start models are considered for this study, since the brightness levels of cold-start planets are too faint for the typical sensitivity of current observations.

Fig. 7 Mean probability of detection of our survey as a function of semimajor axis and planet mass (left) and for specific masses, given above or below each curve in MJup, as a function of semimajor axis (right). The mean probability of detection is obtained by averaging the completeness maps of all targets of the survey. β Pic b (orange square) and the 4 planets of the HR 8799 system (red circles) have been overlaid for reference using the current best estimates of their mass and separation. Note that the HR 8799 planets are represented at their projected physical separation because the true physical separation is not yet known precisely.

More precisely, the predicted magnitude from evolutionary models (Chabrier et al. 2000; Baraffe et al. 2003) – calculated in the CH4s and K′ filters of NIRI or in the Ks filter of NaCo (Baraffe & Allard, priv. comm.) – are used to estimate the planet magnitudes as a function of mass for the age of each target star. Only two extensive evolutionary model grids are available – DUSTY and COND (Chabrier et al. 2000; Baraffe et al. 2003). Of the two evolutionary models, only the COND models extend to sufficiently faint magnitudes and low planet masses to reach the detection limits of our observations, so the analysis in this paper uses the COND models. As an example, the DUSTY and COND mass-magnitude relations at the age of 30 Myr are shown in Fig. 6 which illustrates the limitation of the DUSTY models compared to the sensitivity limit of the observations. Figure 6 also plots the estimates of the masses of the outer three HR 8799 planets based on dynamical stability limits and evolutionary models (Marois et al. 2010b) and compares the masses to the predictions of the DUSTY and COND models. Within the current capacity to perform empirical comparisons, the H-band and K-band magnitudes are consistent with the COND models, as the range of estimated masses for the HR 8799 planets encompasses both the predicted CH4s and K-band magnitudes for that mass range. Of the 42 stars in the sample, 9 were observed with the CH4s filter and 33 were observed with K-band filter. Spectroscopic and photometric observations of the HR 8799 planets have measured differences from the theoretical DUSTY and COND atmospheres and highlighted the need for updated models incorporating non-equilibrium chemistry and an improved treatment of clouds (Janson et al. 2010; Bowler et al. 2010b; Barman et al. 2011a; Galicher et al. 2011), however grids of models including these effects are not currently available.

In the MC simulations, the masses of the generated planets are converted to absolute magnitudes in the filter matching the observations using the evolutionary models, and the signal of each planet is then compared to the 5σ detection limit of the observations at the appropriate angular separation. Since we are using the 1D detection limits, the position angle of the simulated planets is simply ignored. A planet is considered as detected if its observed signal lies above the detection limit at its projected position. The completeness is then calculated as the fraction of detected planets over the number of generated planets.

The two new binaries discovered in the course of our survey, HIP 42334 and HIP 104365, were included into the analysis taking into account the possible disruption induced by the stellar binarity. The limiting values for the semimajor axis of the planets were computed using the analytical expressions of Holman & Wiegert (1999). They define the limiting values that the semimajor axis of a planet can reach and still maintain its orbital stability, as a function of the mass-ratio and orbital elements of the binary. For HIP 42334 and HIP 104365, the ranges of semimajor axes where long-term stability is unlikely are respectively 85−720 AU and 198−1750 AU, so no planets were generated within these intervals.

Finally, for the 13 targets where candidates were detected outside of 320 AU but not followed-up by subsequent observations, planets generated in the MC simulation with a projected separation larger than the closest candidate detected were considered as non-detectable.

5.4. Mean probability of detection

The overall sensitivity of our observations can be estimated in the mass interval 1−14 M_Jup and the semimajor axis interval 5−320 AU. For this analysis, we perform MC simulations for each target assuming no particular distribution for the mass and semimajor axis of the simulated planets. At each point on a linear grid of semimajor axis and mass, 10⁴ planets are simulated following the procedure described in Sect. 5.3, and their detectability is tested with the detection limits to calculate a completeness map for each target of the sample. The mean survey detection thresholds are calculated to illustrate the sensitivity levels, but the completeness of each target is used in the statistical analysis rather than an aggregate sample sensitivity.

The mean detection probability of the full survey is estimated by averaging all the completeness maps for each of the 42 targets. Only the deepest observation of each target was used in the MC simulations, since the second epoch observations were generally only deep enough to redetect the candidates. For the particular case of HR 8799, which is the only target with five very deep observations, the completeness maps for the five epochs were averaged together before being included with the maps of the other targets to avoid giving disproportionate weight to that specific target.

The mean probability of detection for the full survey is represented in Fig. 7 with contour plots as a function of mass (1−30 M_Jup) and semimajor axis (10 − 2500 AU), and cuts at specific masses. For a planet of given mass and semimajor axis, Fig. 7 describes the probability that we would have detected it from our survey data. Observing bright A-stars necessarily results in a greater minimum-mass detection limit, due to the increased contrast compared to solar-type stars. This is clearly visible in Fig. 7, which shows that our data is most sensitive above 10 M_Jup and at separations of 50 AU or more, reaching the highest probabilities of detection between 75 and 300 AU. This is consistent with the median distance of our sample (50 pc), a saturation radius of 0.3 − 0.5′′ for the data, and the fact that the sensitivity limits usually reach a plateau at separations larger than 3 − 4′′. The peak sensitivity reaches 45%, 80% and 94% respectively for 5, 10 and 14 M_Jup. The overall sensitivity to planets similar to HR 8799 bcd and β Pic is low, with a detection probability ranging from 1% for β Pic b to 30% for HR 8799 b, due to the small number of targets covering the appropriate range of orbital separations. The capacity to detect planets like β Pic b is particularly low, and we consider this complication in the frequency calculation in Sect. 5.5. As expected, the sensitivity of the survey to BDs is much higher, with a detection probability above 50% between 75 and 300 AU.

Fig. 8 Representation of the probability density of f given the observations  { dj }  (posterior distribution) as a function of the planetary systems frequency f, in the semimajor axis interval 5−320 AU and mass interval 3−14 MJup, for two detections around HR 8799 and β Pic (see Sect. 5.5 for details). The posterior distribution is obtained from Eq. (2) assuming a linear-flat prior for the frequency and a flat mass distribution for the planets. An equal-tail CI has been calculated using Eqs. (4) and (5) for confidence levels α of 68% (light grey area) and 95% (dark grey area). The dashed line gives the binomial distribution for two detections out of 42 targets, which represent the probability density that would be obtained in an ideal case where the observations would be sensitive to all planets in the range 5 − 320 AU and 3 − 14 MJup.

5.5. Estimation of the planet and brown dwarf systems frequency

In this section we publish the first estimation of the frequency of exoplanetary systems in the mass interval 3−14 M_Jup and the semimajor axis interval 5−320 AU for A-stars. For each target in j = 1...42, the individual target completeness map is used to estimate the mean value of p_j over the 3−14 M_Jup and 5−320 AU intervals. The values of p_j for each target is then used in Eq. (1) to calculate the likelihood for values of f between 0 and 1. Using numerical integration of the likelihood with respect to f, the posterior distribution p(f |  {d_j}) is derived using Eq. (2). Finally, 95% and 68% confidence intervals are calculated using numerical integrations of Eqs. (4) and (5).

Figure 8 presents the probability density function of f calculated from our observations, using a linear-flat prior, compared to a binomial distribution that represents an ideal case in which the observations would be sensitive to all planets in the intervals 5−320 AU and 3−14 M_Jup. We have seen from Sect. 5.4 that our sensitivity to low-mass, close-in planets is limited, so the calculated probability density departs significantly from the ideal binomial distribution. When considering both detections around HR 8799 and β Pic, we estimate the planetary systems frequency around A-stars to $\mathrm{8.}{\mathrm{7}}_{-2.8}^{\mathrm{+}\mathrm{10.1}}\mathrm{\%}$ at a confidence level of 68% (1σ). This value relies on the assumption made in the Monte-Carlo simulations (see Sects. 5.3 and 5.4) that the mass distributions of the simulated planets is flat. We investigate constraints on other possible distributions in Sect. 5.6. Although this probability is relatively high, it is a direct consequence of our decreased sensitivity to short-period planets with low mass. Upcoming large surveys with increased sensitivity to planets in ~5−50 AU orbits will provide a more accurate estimation of the frequency over the full 5 − 320 AU range.

We note that the inclusion of β Pic in the sample must be considered with care. Although it was included as part of the previous deep imaging survey by Chauvin et al. (2010), it is important to note that their observations were not sensitive enough to detect the planet. For that particular target, we reanalyzed the data presented by Bonnefoy et al. (2011), which is more sensitive than the data of Chauvin et al. (2010). This could potentially cause a bias, because it is unknown whether there are similarly sensitive unpublished datasets for the other A-stars included from the Chauvin et al. (2010) survey. An additional bias could occur because β Pic may have been selected among all potential A-stars targets due to indicators correlated with the presence of a planet, such as the known warp in the disk. Because of these complications, we also present an estimation of the frequency considering β Pic as a non-detection³. In that case, we estimate the frequency to $\mathrm{4.}{\mathrm{3}}_{-1.3}^{\mathrm{+}\mathrm{9.1}}\mathrm{\%}$ at 68% confidence level ($\mathrm{4.}{\mathrm{3}}_{-3.3}^{\mathrm{+}\mathrm{17.7}}\mathrm{\%}$ at 95% confidence level).

Finally, since our data is very sensitive to the brown dwarf regime, we can make an estimation of the brown dwarf systems frequency around A-stars including the detection of HR 7329 B in the survey. In the intervals of 14 − 75 M_Jup and 5 − 320 AU, we estimate the frequency to $\mathrm{2.}{\mathrm{8}}_{-0.9}^{\mathrm{+}\mathrm{6.0}}\mathrm{\%}$ at a confidence level of 68%.

5.6. Constraints on giant planets population distributions

Another interesting aspect of our sample is the possibility to place some initial constraints on the distributions of mass and semimajor axis for the giant planet population. The approach taken by previous studies was to try extrapolating the distributions coming from RV surveys, which are usually parameterized with power laws. One of the most complete RV study on this topic for FGKM stars is the one performed by Cumming et al. (2008) using the sensitivity limits from the Keck Planet Search. They determine that the planet frequency for solar-types stars is 10.5% for 0.3 − 10 M_Jup planets with periods less than 1826 days, and that the mass and semimajor axis distributions are best modeled with power laws respectively of index⁴α =  −1.31 and β =  −0.61, i.e. dN ∝ M^-1.31dM and dN ∝ a^-0.61da. For early-type stars, Johnson et al. (2010a) find a positive correlation between planet frequency and primary mass, and determine a frequency of 11 ± 2% for stars in the range 1.3−1.9 M_⊙ (see Table 6 for their measured frequency in different bins of stellar mass). For the distributions of mass and semimajor axis around A-stars, Bowler et al. (2010a) use the detections of their survey to place constraints on the possible values of α and β. They find that the values of Cumming et al. (2008) are inconsistent with their observations, with a probability below 1%, and consequently they favor positive values for both parameters.

For our study we assume the frequency f is known over a certain range of mass and semi-major axis from RV surveys, and we model the mass and semimajor axis distributions as power laws. A cutoff on the semimajor axis distribution (a_cutoff) is added to define an upper limit to the separation at which planets are present. To obtain meaningful results it is important to define the range of mass and semimajor axis over which the frequency f is known and valid. In our case, we use the value reported by Bowler et al. (2010a) and Johnson et al. (2010a) in their survey of giant planets around old A stars: f = 11 ± 2% measured for planets in the ranges 0.5−14 M_Jup and 0.1 − 3.0 AU (median sensitivity of their RV data).

Fig. 9 Contour plots showing the confidence with which we can exclude a model with specific parameters for the mass distribution power law (α, from  −1.5 to 1.5, top to bottom), the semimajor axis distribution power law index (β, from  −2.5 to 0.2) and the semimajor axis distribution power law cutoff (acutoff, from 10 to 320 AU). The hatched grey area corresponds to the values of acutoff within the separation of the known planet HR 8799 b (dashed line). Given the existence of the planet at a separation of  ~68 AU, values of acutoff below that separation are de facto excluded. The best fit values of Cumming et al. (2008) for solar-type stars in terms of mass and semimajor axis are α =  −1.31 and β =  −0.61. Bowler et al. (2010a) report a 50% confidence for values of β higher than zero.

However, our observations are not sensitive to the same range of mass and semimajor axis as the RV data, so the value of f needs to be scaled to be valid inside the range of parameters where we are sensitive, [m_min,m_max] for the planet mass and [a_min,a_max] for the semimajor axis. Once the indices α and β are chosen, the power laws can be integrated over [0.5, 14] M_Jup and [0.1, 3.0] AU for the RV, and over [m_min,m_max] and [a_min,a_max] for the direct imaging. Finally, the value of f is scaled by the ratio of the direct imaging integrated values over the RV integrated values. This normalization ensures that the new frequency, f′, is valid over [m_min,m_max] and [a_min,a_max], and that its value over [0.5, 14] M_Jup and [0.1, 3.0] AU is always equal to f = 11% to match the RV data. In our case, we used the current knowledge of the properties of the detected companions around HR 8799 and β Pic to define the appropriate range of mass and semimajor axis: intervals of [3.5, 12] M_Jup and [8, 68] AU were used to encompass these currently known detections.

For each of the targets in our sample, we simulate different populations of planets with distributions generated on a grid of values for α (−1.5 to 1.5 by steps of 0.1), β (−2.5 to 0.5 by steps of 0.1) and the semimajor axis distribution cutoff, a_cutoff (10 to 320 AU by steps of 10 AU). For each point of the grid, 10⁵ planets are simulated and the fraction of detectable planets in the intervals [3.5, 12] M_Jup and [8, 68] AU is measured. If we assume that all stars of the sample have one planet in that range, the sum of these fractions over our 42 stars would give the number of detection that would be expected from the observations (taking into account the sensitivity around each target). However, we have assumed that the frequency of planetary systems is equal to f′, so the sum needs to be multiplied by f′ to obtain the expected number of detections over [3.5, 12] M_Jup and [8, 68] AU, taking into account the known frequency from RV surveys. Finally, we can use Poisson statistics to obtain the probability that this expected number of detections matches our real detections. In Fig. 9 we represent contour plots showing the confidence with which we can exclude different combinations of parameters given the detections in the range [3.5, 12] M_Jup and [8, 68] AU. Before describing the results, we remind that increasing the value of α is equivalent to increasing the proportion of high-mass planets, and that increasing the value of β is equivalent to increasing the proportion of planets at large orbital separations. We first note that given the detection of the planet HR 8799 b at an orbital separation of 68 AU, we assume that the cutoff for the semimajor axis distribution, a_cutoff, lies at higher orbital separations. In the three contour plots, the combination of parameters that can be excluded with 95% confidence defines a vertical area for a_cutoff ≥ 68 AU. The general trend when increasing α, i.e. when increasing the proportion of high-mass planets, is to move the allowed range of parameters towards smaller values for β.

The top plot of Fig. 9 is close to the value of α =  −1.31 reported by Cumming et al. (2008). With this value of α, their value of β =  −0.61 is not in the range of parameters allowed by our observations. When exploring higher values of α, as suggested by Bowler et al. (2010a), a value of  − 0.61 for β becomes even less likely given a cutoff of a_cutoff ≥ 68 AU. We also note that Bowler et al. (2010a) report higher confidence for positive values of both α and β. However, when increasing α in our simulations, the positive values of β become less likely, suggesting a possible change in the population of exoplanets between the short separation planets detected by RV and the planets at very wide orbital separations around A-stars.

6. Discussion and conclusions

The present survey results on A-star wide orbit planetary systems frequency and population distribution can be compared and combined with previous results to investigate the effects of host star mass, the predictions of theoretical models, the impact of debris disks, and the current state of the planet population across a range of orbital separations. Table 6 summarizes the frequencies and survey sensitivities reported in both AO (current work and Nielsen & Close 2010; which combines the results of Masciadri et al. 2005; Biller et al. 2007; and Lafrenière et al. 2007a) and RV surveys (Johnson et al. 2010b). These studies were selected for comparison, since each covers a large sample and considers the same 68% confidence interval. For the sample with the most similar target mass range, the evolved A-stars, the AO measurement of the wide orbit planetary systems frequency range of 5.9−18.8% encompasses the close orbit planet frequency of 11 ± 2%. Direct imaging searches covering similar orbital separations to this study, but targeting lower mass stars (Nielsen & Close 2010) find an upper limit of  < 20%, while RV searches for close orbit planets measure a declining frequency for the successively lower mass samples (Johnson et al. 2010a). From the current statistics on the wide orbit imaged planet population, it is not possible to determine how the planetary systems frequency scales with host star mass. Theoretical models have predicted a rising (Kennedy & Kenyon 2008), peaked (Ida & Lin 2005), or declining (Kornet et al. 2006) planet frequency on mass, with models incorporating different treatments of factors such as the location and evolution of the snow line, the initial disk size and its dependence on host star mass, and the orbital migration of planets.

Since the survey was designed to detect planets, the full sample is extremely sensitive to brighter brown dwarf companions in wide orbits (probability of detection  >90% in ~75−300 AU, see Fig. 7), enabling an initial assessment of the brown dwarf systems frequency around stars more massive than the Sun. Among the 42 targets, one brown dwarf was previously identified around HR 7329 (Lowrance et al. 2000; Neuhäuser et al. 2011), yielding a brown dwarf systems frequency of $\mathrm{2.}{\mathrm{8}}_{-0.9}^{\mathrm{+}\mathrm{6.0}}$% (see Sect. 5.5). This low level of brown dwarf companions is consistent with the $\mathrm{3.}{\mathrm{2}}_{-2.7}^{\mathrm{+}\mathrm{3.1}}$% frequency measured for a larger scale survey of young (3−3000 Myr), solar-type (F5−K5) stars covering a wider range of separations (28−1590 AU) with an AO survey (Metchev & Hillenbrand 2009).

Theoretical models of the formation of giant planets by gravitational instability have proposed that the directly imaged planets may be the lowest mass products of the same process that forms binary stars and brown dwarfs (Kratter et al. 2010), and these simulations predict a higher proportion of brown dwarfs than planets. While the statistics of the current study are limited, the similarity of the brown dwarf and planetary systems frequencies suggests that the imaged planets may have formed through a different process such as core accretion (e.g. Pollack et al. 1996). An investigation into the conditions for gravitational instability for the specific case of the HR 8799 planets (Nero & Bjorkman 2009) found that the likelihood of this formation mechanism increased with radius towards the outermost planets. Other possible explanations for the wide orbits of some of the directly imaged planets include scattering to larger orbital radii from interactions in multiple planet systems (e.g. Veras & Armitage 2004), outward migration in a resonance with multiple planets (e.g. Crida et al. 2009), and outward migration through a planet-disk interaction (e.g. Veras & Armitage 2004).

The full sample of 42 young, A- and F-stars includes 17 systems with dusty debris disks sustained by the ongoing collisional grinding of planetesimals into smaller particles (Backman & Paresce 1991) or an event such as a catastrophic collision of planets (Cameron 1997; Melis et al. 2010). The 3 targets in this sample with imaged planets or brown dwarfs all reside in systems encircled by dust disks. Additionally, Fomalhaut (not included in this study) would provide another example of a planet-disk system (Kalas et al. 2008) if the presence of a planet was confirmed (see e.g. Janson et al. 2012). While not all A-stars with debris disks harbor massive giant planets in wide orbits, the frequency of planetary systems appears higher among the targets with excess emission from dust. The combination of resolved disk maps and planet imaging is an especially powerful tool to understand the planet-disk interactions that may sculpt young planetary systems (e.g. Liou & Zook 1999; Kuchner & Holman 2003; Wyatt 2006; Quillen & Faber 2006). The structure of the Fomalhaut disk inner edge has been compared with dynamical models of planets with different masses (Chiang et al. 2009), and the size and shape of the HR 8799 disk have been compared with models of the dynamically cleared zones and orbital migration history (Patience et al. 2011). Advances in disk imaging with ALMA and planet imaging with extreme AO will enable more detailed studies to investigate planet-disk interactions which may generate structures such as asymmetries and clumps.

By combining the results and sensitivities of existing direct imaging A-star planet searches with those of RV programs targeting retired A-stars, a comprehensive summary of the currently known A-star planet population is given in Fig. 10. Analysis of the current data already shows distinct differences. The close orbit A-star planets are best fit by a distribution that is flat or rising with increasing orbital separation (Bowler et al. 2010a) and rising with increasing planet mass (Bowler et al. 2010a; Johnson et al. 2010a). In contrast, the wide orbit planet data, incorporating the outer cutoff implied by the outermost HR 8799 planet, are consistent with a distribution that is declining with increasing orbital separation. This result is independent of whether the distribution of planets is rising, flat, or declining with planet mass, as shown in the three panels of Fig. 9 (the allowed region always has a negative power law index). Thus, the existing surveys have identified the boundaries of planets around A-stars – from  ~0.6 AU (Bowler et al. 2010a) to  ~70 AU – and have indicated that the data cannot be fit by a single distribution of properties.

Fig. 10 Comparison of the mean probability of detection for the present survey (plain orange contours) with a simulated survey with GPI of the same 42 targets (dashed blue contours). The mean probability of detection is obtained using Monte-Carlo simulations as described in Sect. 5.4. In addition to β Pic b (orange square) and the HR 8799 planets (red circles), the giant planets detected around old A stars by Johnson et al. (2010b, 2011) have been overplotted (purple triangles). The dotted purple line shows the median detection threshold of the Bowler et al. (2010a) RV survey around old A-type stars. The HR 8799 planets are represented at their projected physical separation because the true physical separation is not yet known precisely. Note that Fomalhaut b has not been included due to its uncertain nature (see e.g. Janson et al. 2012).

Figure 10 also presents the sensitivity of a simulated survey with the upcoming high-contrast imager Gemini/GPI (Macintosh et al. 2008). Using simulated contrast curves for that instrument, we performed Monte-Carlo simulations similar to the ones described in Sect. 5.4 on the 42 targets of our sample. They show that extreme AO imaging will provide the crucial coverage to connect the close orbit and wide orbit populations and to reveal the full distribution of planets as a function of separation and planet mass. Extreme AO data will determine if the transition from a flat or rising population of close orbit planets to a declining population of wide orbit planets is smooth or discontinuous and may discover a peak in the separation distribution. The total number of planets detected in extreme AO surveys, combined with trends measured for close orbit A-star planets, can be used as a test of formation models (e.g. Crepp & Johnson 2011). Since the upcoming extreme AO instrument include integral field units, it will be possible to investigate the nature of the exoplanet atmospheres in addition to the population statistics. Spectra of the currently imaged planetary mass companions have revealed differences with brown dwarfs of similar temperatures (e.g. Janson et al. 2010; Patience et al. 2010; Barman et al. 2011a,b; Skemer et al. 2011), and upcoming AO observations will further explore the architectures and atmospheres of exoplanets around A-stars.

Online material

Appendix A: Point sources properties

Table A.1 summarizes the properties of all point sources detected with our observations.

Acknowledgments

A.V. and J.P. acknowledge support from a Science and Technology Facilities Council (STFC) grant (ST/H002707/1). J.P. acknowledges support from the Leverhulme Trust through a research project grant (F/00144/BJ). The authors would like to thank G. Chauvin for providing his published detection limits, as well as E. Nielsen, M. Viallet and F. Pont for fruitful discussions on the statistical analysis, and J. Johnson for information on comparisons with RV studies. We thank the ESO and Gemini staff for performing the observations. This research made use of the SIMBAD database, operated at CDS, Strasbourg, France.

References

All Tables

All Figures

Fig. 1 MV vs. B − V color − magnitude diagram of the A-stars of our target sample. The stars with published ages are plotted with a green circle. Stars without a published age are noted with blue triangles and have been assigned an age of 125 Myr, based on their positions within the region defined by the Pleiades A-stars, which are indicated with orange stars. An empirical linear fit to the Pleiades is also shown. The older population of 625 Myr Hyades A-stars is plotted with blue squares.

In the text

	Fig. 2 Age, distance and spectral type of all the stars in our sample. Table 1 gives the individual properties of each target. The size and color of the symbols correspond to the spectral type of each target: bluer colors for earlier spectral types and redder colors for later spectral types. The points grouped at 125 Myr correspond to stars without a published age, as described in Sect. 2.
In the text

	Fig. 3 Example of final image after data reduction and processing with LOCI for target HIP 10670. The central white region corresponds to the saturated area in the science images, which has been masked in the data analysis. The point source 4.7′′ West of the star is a background object 14.6 mag fainter than the central star.
In the text

	Fig. 4 Examples of K-band 5σ detection limits representing the best (HIP 69732) and worst (HIP 15648) limits, as well as an intermediate case (HIP 26624). The right axis shows the detection limit in masses of Jupiter for the median target of our survey, a 100 Myr-old A3V star at 50 pc, converted using the COND evolutionary models (Baraffe et al. 2003).
In the text

	Fig. 5 Offset position (black circles) with uncertainties measured relative to the primary star for the co-moving stellar companions to HIP 42334 (left) and HIP 104365 (right). The variation expected for a background stationary object is also plotted (grey line) with the expected positions of background objects (grey square).
In the text

Fig. 6 The detection limit (dashed line) for wide orbit planets from the observations of the target HR 8799 in the NIRI CH4s filter (left) and K′ filter (right) compared to the DUSTY (red line) and COND (black line) evolutionary models (Chabrier et al. 2000; Baraffe et al. 2003). The DUSTY model mass-magnitude grid does not extend to magnitudes as faint as the detection limits. The current estimates of the HR 8799 planet masses (Marois et al. 2010b) are also plotted, and the empirical masses are consistent with the COND magnitudes within the error bars.

In the text

Fig. 7 Mean probability of detection of our survey as a function of semimajor axis and planet mass (left) and for specific masses, given above or below each curve in MJup, as a function of semimajor axis (right). The mean probability of detection is obtained by averaging the completeness maps of all targets of the survey. β Pic b (orange square) and the 4 planets of the HR 8799 system (red circles) have been overlaid for reference using the current best estimates of their mass and separation. Note that the HR 8799 planets are represented at their projected physical separation because the true physical separation is not yet known precisely.

In the text

Fig. 8 Representation of the probability density of f given the observations  { dj }  (posterior distribution) as a function of the planetary systems frequency f, in the semimajor axis interval 5−320 AU and mass interval 3−14 MJup, for two detections around HR 8799 and β Pic (see Sect. 5.5 for details). The posterior distribution is obtained from Eq. (2) assuming a linear-flat prior for the frequency and a flat mass distribution for the planets. An equal-tail CI has been calculated using Eqs. (4) and (5) for confidence levels α of 68% (light grey area) and 95% (dark grey area). The dashed line gives the binomial distribution for two detections out of 42 targets, which represent the probability density that would be obtained in an ideal case where the observations would be sensitive to all planets in the range 5 − 320 AU and 3 − 14 MJup.

In the text

Fig. 9 Contour plots showing the confidence with which we can exclude a model with specific parameters for the mass distribution power law (α, from  −1.5 to 1.5, top to bottom), the semimajor axis distribution power law index (β, from  −2.5 to 0.2) and the semimajor axis distribution power law cutoff (acutoff, from 10 to 320 AU). The hatched grey area corresponds to the values of acutoff within the separation of the known planet HR 8799 b (dashed line). Given the existence of the planet at a separation of  ~68 AU, values of acutoff below that separation are de facto excluded. The best fit values of Cumming et al. (2008) for solar-type stars in terms of mass and semimajor axis are α =  −1.31 and β =  −0.61. Bowler et al. (2010a) report a 50% confidence for values of β higher than zero.

In the text

Fig. 10 Comparison of the mean probability of detection for the present survey (plain orange contours) with a simulated survey with GPI of the same 42 targets (dashed blue contours). The mean probability of detection is obtained using Monte-Carlo simulations as described in Sect. 5.4. In addition to β Pic b (orange square) and the HR 8799 planets (red circles), the giant planets detected around old A stars by Johnson et al. (2010b, 2011) have been overplotted (purple triangles). The dotted purple line shows the median detection threshold of the Bowler et al. (2010a) RV survey around old A-type stars. The HR 8799 planets are represented at their projected physical separation because the true physical separation is not yet known precisely. Note that Fomalhaut b has not been included due to its uncertain nature (see e.g. Janson et al. 2012).

In the text