Combining high contrast imaging and radial velocities to constrain the planetary architecture of nearby stars

TL;DR

This study combines high contrast imaging and radial velocity data to better constrain the planetary architectures of six nearby stars, significantly improving detection limits especially for intermediate semi-major axes.

Contribution

It demonstrates the effectiveness of integrating HCI and RV data, using a Monte Carlo approach to enhance constraints on exoplanet properties around nearby stars.

Findings

HCI data increased completeness by up to 99% for some stars.

Constraints are strongest for semi-major axes of 15-40 AU.

Mass limits reach 5-20 Jupiter masses depending on star age.

Abstract

Nearby stars are prime targets for exoplanet searches and characterization using a variety of detection techniques. Combining constraints from the complementary detection methods of high contrast imaging (HCI) and radial velocity (RV) can further constrain the planetary architectures of these systems because these methods place limits at different regions of the companion mass and semi-major axis parameter space. We aim to constrain the planetary architectures from the combination of HCI and RV data for 6 nearby stars within 6 pc: $\tau$ Ceti, Kapteyn's star, AX Mic, 40 Eri, HD 36395, and HD 42581. We compiled the sample from stars with available archival VLT/NACO HCI data at L$^{\prime}$ band (3.8 $\mu$m). The NACO data were fully reanalyzed using the state-of-the-art direct imaging pipeline PynPoint and combined with RV data from HARPS, Keck/HIRES, and CORALIE. A Monte Carlo approach was used to assess the completeness in the companion mass/semi-major axis parameter space from the combination of the HCI and RV data sets. We find that the HCI data add significant information to the RV constraints, increasing the completeness for certain companions masses/semi-major axes by up to 68 - 99% for 4 of the 6 stars in our sample, and by up to 1 - 13% for the remaining stars. The improvements are strongest for intermediate semi-major axes (15 - 40 AU), corresponding to the semi-major axes of the ice giants in our own solar system. The HCI mass limits reach 5 - 20 $M_{\textrm{Jup}}$ in the background-limited regime, depending on the age of the star. Through the combination of HCI and RV data, we find that stringent constraints can be placed on the possible substellar companions in these systems. Applying these methods systematically to nearby stars will quantify our current knowledge of the planet population in the solar neighborhood and inform future observations.

Figures

Tables

Table 1. Table 1: Stellar properties

Star	Spectral type	Distance	Apparent magnitude	Age [min. - max.]	Mass	References
		(pc)	($L^{\prime }$)	(Gyr)	(M⊙)
$\tau$ Ceti	G8V	3.65	1.7	5.8 [2.9 - 8.7]	0.63	1, 2, 3, 4, 5
Kapteyn’s star	M1VIp	3.93	4.9	1.7 [0.5 - 2.9]	0.27	5, 6, 7, 8, 9
AX Mic	M1V	3.97	2.9	4.8 [1.9 - 7.7]	0.56	5, 6, 8, 10, 11, 12
40 Eri	K0V	4.98	2.4	5.6 [2.8 - 8.4]	0.67	1, 3, 4, 5, 10
HD 36395	M1.5Ve	5.70	3.8	5.8 [2.7 - 8.8]	0.64	5, 6, 7, 9, 13
HD 42581	M1V	5.76	4.0	1.6 [0.3 - 3.0]	0.58	6, 7, 9, 14, 15

Table 2. Table 2: Summary of Archival NACO Imaging Observations

Star	Date	$t_{\text{int}}$/frame	# of frames	# of frames	Total $t_{\text{int}}$	Airmass	Field rotation	Program ID
	(UT)	(sec)	used	removed	(min)		(deg)
$\tau$ Ceti	2008-11-16	0.175	4137	3	12	1.02 - 1.07	19	082.C-0518
Kapteyn’s star	2008-11-16	0.500	2815	20	23	1.10 - 1.14	15	082.C-0518
AX Mic	2009-10-11	0.175	3423	10	10	1.07 - 1.12	14	084.C-0656
40 Eri	2008-11-16	0.500	1268	1693	11	1.05 - 1.07	21	082.C-0518
HD 36395	2008-11-16	0.300	3842	1	19	1.08 - 1.12	18	082.C-0518
HD 42581	2008-11-16	0.300	3723	120	19	1.00 - 1.01	53	082.C-0518

Table 3. Table 3: Summary of RV Observations

Star	# of measurements				$\Delta$t${}_{\text{RV}}$a𝑎aa𝑎a$a$Total time baseline of RV observations.	$a$ for $\Delta$t${}_{\text{RV}}$b𝑏bb𝑏b$b$Semi-major axis corresponding to an orbital period that equals the time baseline of the RV observations given the stellar mass in Table 1.	${\sigma }_{\text{RV}}$
	HARPS	HIRES	CORALIE	total	(days)	(AU)	(m/s)
$\tau$ Ceti	1038	272	314	1624	7353	6.3	2.8
Kapteyn’s star	125	31	0	156	7046	4.7	3.3
AX Mic	0	8	0	8	740	1.3	2.6
40 Eri	0	109	0	109	4766	4.9	12.4
HD 36395	0	41	0	41	6220	5.7	6.0
HD 42581	198	44	0	242	7240	6.1	4.4

Table 4. Table 4: PCA PSF subtraction parameters

Star	Stackinga𝑎aa𝑎a$a$Number of frames that were stacked prior to PSF modeling and subtraction, or none if no stacking was performed.	PCs usedb𝑏bb𝑏b$b$Quantity of principal components used to model the stellar PSF, expressed as a percentage of the total number frames or as the number of components used.		Mask sizec𝑐cc𝑐c$c$Radius of the central mask.
		(%)	(#)	(arcsec)
$\tau$ Ceti	none	10	414	0.30
Kapteyn’s star	none	1	28	0.37
AX mic	10	2.5	9	0.40
40 Eri	none	0.5	6	0.27
HD 36395	20	7.5	14	0.32
HD 42581	none	5	186	0.10

Full text

*Aims. *We aim to constrain the planetary architectures from the combination of HCI and RV data for 6 nearby stars within 6 pc: $\tau$ Ceti, Kapteyn’s star, AX Mic, 40 Eri, HD 36395, and HD 42581. We explored where HCI adds information to constraints from the long-term RV monitoring data for these stars.

Methods. We compiled the sample from stars with available archival VLT/NACO HCI data at L′ band (3.8 $\mu$m), where we expect substellar companions to be brighter for the typically older ages of nearby field stars (¿1 Gyr). The NACO data were fully reanalyzed using the state-of-the-art direct imaging pipeline PynPoint and combined with RV data from HARPS, Keck/HIRES, and CORALIE. A Monte Carlo approach was used to assess the completeness in the companion mass/semi-major axis parameter space from the combination of the HCI and RV data sets.

*Results. *We find that the HCI data add significant information to the RV constraints, increasing the completeness for certain companions masses/semi-major axes by up to 68 – 99% for 4 of the 6 stars in our sample, and by up to 1 – 13% for the remaining stars. The improvements are strongest for intermediate semi-major axes (15 –- 40 AU), corresponding to the semi-major axes of the ice giants in our own solar system. The HCI mass limits reach 5 –- 20 $M_{\textrm{Jup}}$ in the background-limited regime, depending on the age of the star.

*Conclusions. *Through the combination of HCI and RV data, we find that stringent constraints can be placed on the possible substellar companions in these systems. Applying these methods systematically to nearby stars will quantify our current knowledge of the planet population in the solar neighborhood and inform future observations.

1 Introduction

Nearby stars (¡20 pc) are particularly interesting targets for exoplanet searches and characterization using a variety of techniques. These stars are bright enough for precise radial velocity (RV) measurements, and many have been the targets of long-term RV monitoring with high-resolution spectrographs such as HARPS (Pepe et al. 2002; Mayor et al. 2003) and Keck/HIRES (Vogt et al. 1994). The proximity of these stars also improves the companion masses that can be detected through direct or high contrast imaging (HCI) in the background-limited regime by increasing the apparent magnitude of their thermal emission (e.g., Thalmann et al. 2011; Quanz et al. 2012; Mawet et al. 2019). The current space missions TESS (Ricker et al. 2015) and Gaia (see Perryman et al. 2014) are specifically focused on searches for planets around nearby stars using the transit and astrometry techniques, respectively. Future large-aperture missions now under study, such as HabEx and LUVOIR, or a space-based mid-infrared interferometer will also target nearby stars to directly image and characterize many Earth-like planets (Mennesson et al. 2016; Pueyo et al. 2017; Kammerer & Quanz 2018). Constraining the planet population in the solar neighborhood using current data and facilities will help prioritize targets for these missions as well as future ground-based instruments and telescopes such as VLT/ERIS (Davies et al. 2018), ELT/METIS (Quanz et al. 2015; Brandl et al. 2018), and ELT/HiRES (Marconi et al. 2018).

To quantify our current constraints on the planetary architecture of nearby stars, it is advantageous to combine the various exoplanet detection methods since these probe complimentary parameter space in companion mass and semi-major axis and can give information on different exoplanet parameters. The RV method is most sensitive to close-in planets and yields estimates of the planet’s minimum mass and some of its orbital parameters, while HCI is most sensitive to massive planets at wide separations and measures photons from the planet itself that can be used to characterize the planet’s atmosphere and estimate its mass. When a planet is detected using both RV and HCI, and more recently by adding information from Hipparcos and Gaia Data Release 2, it is possible to break the planet mass/inclination degeneracy from RV data alone and the resulting dynamical mass can then be compared to the masses derived from evolutionary and atmospheric models to further refine these predictions (e.g., Crepp et al. 2012; Bowler et al. 2018; Cheetham et al. 2018; Snellen & Brown 2018; Brandt et al. 2018; Dupuy et al. 2019). Even for cases in which no planet is detected, the complimentary nature of these different detection methods can put substantial constraints on the possible masses and semi-major axes of substellar companions in individual systems.

In this work, we focus on combining constraints from HCI and RV data to build a census of planets in the solar neighborhood. The sample presented here consists of 6 stars that are all within 6 pc, which have been targets of long-term RV monitoring campaigns and have available archival HCI data. We specifically consider stars with available HCI data at the L*′* band ($\lambda_{0}$ = 3.8 $\mu$m) from the NACO instrument at the VLT (Rousset et al. 2003; Lenzen et al. 2003). Cool substellar companions are expected to be brighter at these thermal infrared (IR) wavelengths than at shorter wavelengths and thus these data can reach better mass constraints even with weaker contrast limits compared to extreme adaptive optics (AO) near-infrared instruments such as VLT/SPHERE and Gemini/GPI (e.g., Heinze et al. 2010; Vigan et al. 2015; Beuzit et al. 2019; Macintosh et al. 2014). Here we consider the constraints on a system’s planetary architecture that can be derived from the combination of HCI and RV data when no massive companion is detected, in contrast to the works mentioned above that combined HCI and RV data to measure dynamical masses of known companions.

This paper is organized as follows. In Sec. 2, the archival HCI data and the RV data are described. In Sec. 3, we derive mass limits from the HCI data and combine these with companion mass/semi-major axis constraints from the RV data. The resulting limits on the planetary architecture of these systems are presented in Sec. 4 and discussed in Sec. 5.

2 Sample and data reduction

2.1 High contrast imaging data

The sample of stars for this work was constructed by searching the ESO archive for programs targeting very nearby stars ($<6$ pc) with the NACO instrument at the VLT. NACO consists of NAOS, an AO system, and CONICA, a 1 – 5 $\mu$m imager and spectrograph. As discussed above, we selected stars with data taken in the L*′* filter ($\lambda_{0}$ = 3.80 $\mu$m, $\Delta\lambda$ = 0.62 $\mu$m), since planets are expected to be brighter at these thermal infrared wavelengths for the old ages ($>1$ Gyr) of solar neighborhood stars. It was also required that the data were taken in pupil tracking mode to enable angular differential imaging (ADI; Marois et al. 2006) and in “cube” mode, so that each individual short (0.175 – 0.5 sec) exposure is written out. The properties of the resulting sample of 6 stars is detailed in Table 1.

The NACO HCI data for all 6 stars in the sample are summarized in Table 2. The data primarily come from the ESO 082.C-05182 program (PI: Apai), which targeted all sources within 6 pc and RV planet hosts within 13 pc. The data on the star AX Mic are from the ESO 084.C-0656 program (PI: Apai), which targeted the HR 8799 planets. All the archival data analyzed in this work have not previously been published. The data were taken with the AO system using the target star as the natural guide star and without a coronagraph. These observations also employed a dither pattern that moves the star from the center of one detector quadrant to another every few data cubes to facilitate the subtraction of the strong sky background at L*′* . Due to the lack of coronagraph and the brightness of these nearby stars, all data are saturated in the stellar core.

The raw NACO HCI data were downloaded from the ESO archive and reduced using the state-of-the-art direct imaging data pipeline PynPoint (Amara & Quanz 2012; Stolker et al. 2019, v0.5.3). PynPoint includes modules for calibration steps as well as modules to subtract the stellar point spread function (PSF) using Principal Component Analysis (PCA). All data were dark, flat, and bad pixel corrected, and then the sky was subtracted using an average of the preceding and subsequent data cubes in which the star was at a different dither position. Automatic frame selection was performed to remove any frames in which the AO loops were open. This frame selection was done by performing aperture photometery of the star in each frame using a circular aperture with a radius of 0$\aas@@fstack{\prime\prime}$2, and then removing frames in which the measured stellar counts differed from the median by more than 4 times the standard deviation. Manual frame selection was also performed for 40 Eri because these data suffered from a tracking issue that caused the star to be close to the edge of the detector in over half of the data set. Next, the star was located in each frame and an area of 4″ by 4″ around the star was cropped out. The resulting images were aligned relative to each other using cross correlation of the full 4″ by 4″ images. The aligned images were then centered by fitting a Gaussian to a central circular region of 1″ in radius in the average of all images. The cross correlation and the Gaussian fit were performed over large regions of the images so that the alignment and centering are accurate despite the core of the stellar PSF being saturated. Finally, the stellar PSF was modeled using PCA after masking out a central region corresponding to minimum separation at which a planet could be detected as determined by the field rotation of the data (see Sec. 3.1 for details). The final reduced and PSF-subtracted image for each star is shown in Fig. 1.

To search for substellar companions at projected separations outside of 2″, classical ADI was performed (specifically, the mean of all images was subtracted from each image before de-rotating and averaging; Marois et al. 2006). At these larger separations from the star, the subtraction of the stellar PSF is less critical as the data are in the background-limited regime and no longer the contrast-limited regime. Prior to the ADI PSF subtraction, the individual frames were cropped to 13″ by 13″ and the same central region as above was masked out. This 13″ field of view (FOV) is set by the size of one quadrant of the NACO detector at the plate scale of $\sim$ 27 mas. The classical ADI subtraction was not performed for 40 Eri because the star was too close to the edge of the detector to probe separations larger than 2″ for all position angles. The final classical ADI images are shown in Fig. 2.

2.2 RV data

RV data were taken from the Data Analysis Center for Exoplanets222https://dace.unige.ch/ (site accessed on 11 January 2019). (DACE) hosted by the University of Geneva, which compiles available RV measurements from a variety of sources. The RV data for the targets in this work consist of new and previously published HARPS data (for $\tau$ Ceti, 40 Eri, HD 42581, and Kapteyn’s Star), published Keck/HIRES data (Butler et al. 2017, for all stars), and CORALIE data (for $\tau$ Ceti). CORAVEL data for $\tau$ Ceti are also available in the DACE database, but these data were removed due to their significantly higher scatter compared to the other data sets (270 m/s versus 2.8 m/s).

Using the analysis tools available on DACE, the RV data were binned on a nightly basis. A nightly binning was chosen because this work focuses on companions with orbital periods much longer than a 1 day timescale. The binned RV time series were then downloaded directly from DACE and then the secular acceleration of each star was corrected given the star’s parallax and proper motion from the parallax reference listed in Table 1. The RV offsets of the different instruments were then corrected. For this step, each data release from each instrument (e.g., COR98, COR07, and COR14 from the CORALIE instrument) was treated independently, since upgrades or other changes to the instrument can introduce different RV offsets. Each RV offset was determined by adding the data sets one by one from the earliest to the latest starting time and setting the mean of the temporarily overlapping data points to 0. The HARPS15 data sets for Kapteyn’s star and HD 42581 do not overlap in time with any other RV data set, so for these stars, the HARPS15 data were instead anchored to the HARPS03 data using an average of the offsets measured for the two RV standard stars with an early M spectral type (M2.5V) similar to Kapteyn’s star and HD 42581 in Lo Curto et al. (2015). For these HARPS15 data, the dispersion of the measured RV offsets for the standard stars (2 m/s) is added in quadrature to the HARPS15 error. We note that the uncertainty in the HARPS15 RV offset is likely even larger for a given star, since the value depends on both the spectral type and the $v\leavevmode\nobreak\ \textrm{sin}\leavevmode\nobreak\ i$ of the star. Finally, the overall standard deviation of the calibrated RV data was taken for use in the subsequent analysis.

The resulting calibrated RV data are shown in Fig. 3 and Table 3 summarizes the properties of these RV data for each target star. The time baseline of the RV observations ranges from just over 2 years for AX Mic to $\sim$20 years for $\tau$ Ceti. The scatter in the RV measurements also spans a wide range from 2.8 to 12.4 m/s. The larger RV scatter values are likely due to stellar activity. Instrumental effects could also play a roll in the case of 40 Eri, which exhibits a jump in the measured RV in the most recent HIRES epochs that is uncorrelated with activity index changes and that drives this star’s high RV scatter.

3 Analysis

3.1 High contrast imaging contrast and mass limits

We first searched the PSF-subtracted images for any companion candidates. A point source is detected in the HD 42581 data to the SSE whose position is consistent with the known brown dwarf companion HD 42581 B (Nakajima et al. 1995; Oppenheimer et al. 1995). We also identified a point source to the west of HD 36395 with a separation of 3.5″, a PA of 276 degrees, and a S/N of 16.5 (where S/N is defined as in Mawet et al. 2014). We searched for possible background sources coinciding with this candidate in a range of published catalogues and found that the position of this detection is consistent with that of a faint optical source (B mag = 19.9) in the Guide Star Catalog V2.3 within the reported errors (GSC2.3 identifier: S1LI000310; Lasker et al. 2008). The position is also consistent with a source in the GAIA DR2 catalogue (BP = 19.2 mag, GAIA DR2 identifier: 3209938267881341056; Gaia Collaboration et al. 2018). This GAIA source has a measured parallax of $0.7\pm 0.4$ mas, which corresponds to a distance of 0.9 – 3 kpc. Based on these catalogue identifications, we conclude that this candidate companion is a background star. Finally, a second possible point source to the SE of HD 36395 is visible at a separation of 0.84″and a PA of 145 degrees, but this detection is not significant (S/N = 3.2). Having found no new significant point sources in the data, we focus the rest of the analysis and discussion on the HCI detection limits.

Contrast limits for the final HCI images were derived by planting artificial companions in the individual frames, performing the full PCA or ADI PSF subtraction, and then computing the resulting detection significance. This method accounts for the self-subtraction of the planet signal that occurs when the stellar PSF is modeled using PCA or subtracted using ADI. We tested a range of values for the number of principal components (PCs) and the number of images that were temporally stacked prior to PSF subtraction to determine the optimal parameters for each data set. For these tests, the limiting contrast was computed at 4 evenly-spaced separations from the edge of the central mask (0$\aas@@fstack{\prime\prime}$1 – 0$\aas@@fstack{\prime\prime}$4; see Table 4) out to the edge of the 4 by 4″ field of view and along a position angle of 0 deg for a number of PCs ranging from 0.5 to 40% of the total frames and a number of stacked frames from 1 to 40. For each stacking value, the number of PCs that yielded the two best contrast limits were then selected and tested at 3 additional position angles. Finally, the optimal combination of stacking and number of PCs was chosen by prioritizing the wider 2 separations, since these separations are more likely to contribute additional constraints to those from the RV data alone. We note that by optimizing for the contrast at wide separations, the closest separations can have substantial self-subtraction and in some cases a detection with the specified significance could not be achieved. For the classical ADI PSF subtraction, we also tested the same range of stacking values and found that changing the number of stacked frames did not affect the final contrast limits. In the case of HD 42581, we masked out the region around HD 42581 B when determining the noise at the widest separation. The final parameters for each star used in the PSF subtraction are presented in Table 4.

The final contrast curves with the optimal PSF subtraction parameters were computed at steps of 0$\aas@@fstack{\prime\prime}$1 (a total of 40 – 60 separations for stars other than 40 Eri) by averaging over the contrasts at 3 positions angles (0, 120, and 240 degrees). We used an average of 3 position angles due to constraints on computational time: determining the contrast for each image position required running the full PSF subtraction multiple times to find the contrast with the specified significance. We computed the contrast for a constant false positive fraction (FPF) of $2.9\times 10^{-7}$, while accounting for the small sample statistics at close separations (Mawet et al. 2014). This FPF is equivalent to a 5-sigma detection threshold for an underlying Gaussian distribution and a large number of samples. We note that contrast curves with a constant FPF differ from the traditional contrast curves that fix the detection threshold as a function of separation and thereby allow the FPF to vary. Any image positions for which the specified FPF could not be reached due to self-subtraction were subsequently ignored.

In the end, the contrast curves from the PCA and ADI PSF subtractions were combined and then converted to mass limit curves assuming the star’s age and L*′* magnitude and using the COND evolutionary models (Baraffe et al. 2003). The COND model was chosen primarily for convenience: many of the available models do not go down to low enough effective temperatures for the L*′* contrast limits reached in this work. We note that although the COND model assumes a hot start for the planet formation, the different predictions of the hot and cold start models have converged for the ages of the stars in the sample and for the mass limits we achieve (see Marley et al. 2007).

3.2 Radial velocity mass limits

To determine whether a planet with a given mass and semi-major axis would be detected by the RV data, we consider the overall scatter of the RV data and the epochs of the RV observations. For a given planet, we predict the stellar RV at each observation epoch and find the maximum difference in RV measured across the observations. If this maximum difference in the predicted RV at the observations epochs is larger than 5 times the standard deviation of the measured RVs (reported in Table 3), then the planet signal is considered detectable by the RV data. We use this criterion for an RV detection instead of one derived from the typical Lomb-Scargle periodogram (Lomb 1976; Scargle 1982) because the semi-major axes of interest correspond to orbital periods that are generally longer than the RV time baseline (for $a$ ¿ 1 – 6 AU), and thus these RV signals would not be periodic in the RV data. For the largest semi-major axes explored here of 50 – 70 AU only 1 – 5% of the orbit is covered by the RV time baseline, making a periodogram a poor method for detecting these RV signals.

Prior to determining the standard deviation of the measured RVs, we consider the impact of known planets on the measured RV signal. Three of the six stars in the sample have had planets detected in their RV data with minimum masses in the super Earth/mini Neptune range: $\tau$ Ceti with 4 planets (1.8 – 3.9 $M_{\oplus}$, Tuomi et al. 2013; Feng et al. 2017), Kapteyn’s star with 2 planets (4.8 and 7 $M_{\oplus}$, Anglada-Escude et al. 2014; but Robertson et al. 2015 found that the lower mass planet signal is due to stellar activity) and 40 Eri with 1 planet (8.5 $M_{\oplus}$, Ma et al. 2018). All the measured RV semi-amplitudes for these planets are a factor of 3 – 6 lower than the overall scatter measured in the data, so we do not expect these signals to change the detectabilty of the long orbital period companions as determined by the criterion above. Tuomi et al. (2014) also reported a planet with a mass of 1.9 times the mass of Neptune around HD 42581 based on combined UVES and HARPS RV data; however, we cannot reproduce the planet signal with the HARPS and Keck/HIRES RV data used in this work. We therefore analyze the HD 42581 data without removing a planetary signal, but we note that if this planet is real then the RV constraints for this target would improve.

The brown dwarf companion to HD 42581 (mass = 29 – 39 $M_{\textrm{Jup}}$; Nakajima et al. 2015) could also have an effect on the RV signal of this star. This effect would most likley be in the form of a linear RV trend given the companion’s wide separation and the time baseline of the RV data, with the slope of this trend highly depending on the companion’s orbital parameters. Considering only circular orbits and given the approximate separation of HD 42581 B in the NACO data ($\sim$40 AU), the linear trend could be as high as $\sim$2 m/s/year; however, we do not find evidence for such a trend in the RV data. The non-detection of an RV signal is consistent with companion orbits that have semi-major axes larger than about 80 – 100 AU when using the same detection criterion used above (i.e., the maximum difference in RV being larger than 5 times the overall RV standard deviation).

3.3 Combining constraints from high contrast imaging and RV

To combine the constraints on planetary architecture from HCI and RV, we employ a Monte Carlo approach to find the percentage of companions that would be detected for a given mass and semi-major axis ($a$). We explored a grid of companion masses ranging from 1 to 50 $M_{\textrm{Jup}}$ and semi-major axes from 1 to 45 AU (37 AU for 40 Eri, adjusted for the smaller FOV of this data). For each combination of mass and semi-major axis values, we randomly sample the other orbital parameters and then determine whether a planetary/substellar companion with this mass and orbit would have been detected with each method. Only circular orbits are considered ($e$ = 0), and the longitude of the ascending node ($\Omega$) is fixed because this orbital parameter does not affect the RV signal or the projected separation between the star and the companion in the HCI data. The time of closest approach ($T_{0}$) is also fixed to an arbitrary value, since this parameter is degenerate with the argument of periapse ($\omega$) for circular orbits. The mass of the star is set to the literature value given in Table 1. For each combination of planet mass and semi-major axis, we sample the remaining orbital parameters of inclination ($i$; from a sin($i$) distribution for randomly oriented orbital planes) and $\omega$ (from a uniform distribution spanning 0 to 360 deg) 10,000 times.

For each of the 10,000 random orbits, the predicted projected separation at the epoch of the HCI observation and the predicted RV signal at the epochs of all RV measurements are found. To determine whether a companion on a given orbit would have been detected in the HCI data, we compare the companion mass to the mass limit curve at the projected separation. If the mass is higher than the mass limit, this planet is considered detectable by the HCI data. We apply the criteria outlined in Sec. 3.2 to determine if the same planet is detectable by the RV data alone. Finally, we find the percentage of planets detectable by each method to measure the completeness for each combination of mass and semi-major axis.

4 Results

4.1 Contrast and mass limits from HCI data

The limiting L*′* contrast curves for all stars in the sample are shown in Fig. 4. We constrain the presence of companions out to projected separations of 2.0″ for 40 Eri and to 6.5″ for the remaining targets. These angular separations correspond to physical separations of 10 AU and 25 – 35 AU, respectively, for the very nearby distances of these stars. The limiting contrasts for a fixed FPF of $2.9\times 10^{-7}$ go down to $\sim$14 magnitudes in the background-limited regions, which are at angular separations of 1 – 3″ depending on the brightness of the target star. To explore the completeness of the NACO HCI data alone, we also present performance maps in Fig. 5, following the method described in Jensen-Clem et al. (2018). These performance maps present the true positive fraction (TPF) as a function of separation and contrast for the fixed FPF of $2.9\times 10^{-7}$, allowing a more thorough exploration of the trade-off between FPF and TPF. The 50% completeness curve is highlighted, which is the same completeness shown in Fig. 4 by definition, as well as the 95% completeness. In both the final contrast curves and performance maps, separations for which the specified FPF could not be reached due to self-subtraction are not plotted.

Fig. 6 shows the 50% completeness contrast limits converted into mass limits given the star’s age, distance, and apparent magnitude. The ages of the stars are quite uncertain with typical errors of $\sim$50 – 70%, and this results in a spread in the limiting mass of about 10 $M_{\textrm{Jup}}$ . The mass limits reach down to 2 – 10 $M_{\textrm{Jup}}$ for the stars with youngest ages (¡3 Gyr; Kapteyn’s star and HD 42581). The remaining stars are likely to be older than 2 – 3 Gyr, resulting in mass limits of 5 – 10 $M_{\textrm{Jup}}$ in the background limit, assuming the youngest possible ages, or 15 – 20 $M_{\textrm{Jup}}$ , assuming the oldest possible ages up to 9 Gyr.

4.2 Planetary architecture constraints from combining HCI and RV

The percentage of the 10,000 orbits for each mass/semi-major axis combination that can be detected by HCI, RV, and by either HCI or RV are plotted in Figs. 7 – 12. Semi-major axes out to 2 times (4 times for 40 Eri) the maximum projected separation probed by the NACO FOV are explored, since orbits with $a$ greater than the FOV may also fall into the NACO frame depending on the orbital phase at the time of the observation. We focus on large semi-major axes that correspond to longer orbital periods than are fully covered by the RV data, since this is the parameter space where HCI could improve constraints when combined with RV. Over these semi-major axes, we find that the limiting mass of the RV data increases quickly as a function of semi-major axis, whereas the HCI constraints are relatively flat for semi-major axes corresponding to the projected physical separations of the background limit. For 4 out of 6 of the stars, HCI adds significant information to the RV constraints (AX mic, Kapteyn’s star, HD 36395, HD 42581). The additional constraints from HCI are concentrated at low companion masses near the HCI detection limit (5 – 20 $M_{\textrm{Jup}}$ ) and large semi-major axes out to the maximum separation probed by the NACO FOV ($\sim$25 – 40 AU). When including the information from HCI, the percentage of detectable substellar companions increases by up to 50 – 90% for Kapteyn’s star, HD 36395, and HD 42581 depending on the age assumed. For AX mic, which is the star with the fewest RV measurements (8 measurements after nightly binning) and the shortest RV time baseline (741 days), the percentage of companions that can be detected increases by up to 99% irrespective of the stellar age.

For the stars $\tau$ Ceti and 40 Eri, the addition of HCI data increases the detectable orbits by at most 10%. The HCI mass limits for 40 Eri only go out to projected physical separations of $\sim$10 AU, and the RV data provide strong constraints for semi-major axes smaller than this value despite the relatively large scatter in the data (12.4 m/s). $\tau$ Ceti is well known as an extremely stable RV star, and of the stars in this sample it has the longest RV time baseline (7,354 days or $\sim$19 years) and lowest RV scatter (2.8 m/s). Thus the RV data of this star place very powerful constraints on the planetary architecture even at wide semi-major axes up to 40 AU (orbital periods of over 300 years). Still, for 40 Eri and $\tau$ Ceti, HCI can detect a few tens of the 10,000 sampled orbits that RV misses for small companion masses and large semi-major axes. The orbits at which HCI can detect planets tend to have small inclinations (less than 20 – 30 deg) and projected separations close to the maximum possible value at the time of the HCI observations. It is possible to detect these companions in HCI but not in RV because at these large semi-major axes, the RV time baseline only covers a small fraction of the orbit (down to 13% for $a$ ¡ the NACO FOV) and thus the smallest changes of RV over the baseline occur around the times of maximum and minimum RV value, where the derivative of the RV is 0. This is also the section of the orbit where the projected separation is the highest, allowing some companions to be detectable by HCI.

HCI can even add information for semi-major axes beyond the extent of the image for lower inclinations, due to the fact that there is a fraction of orbits for which the projected separation at the time of the HCI observation still falls in the FOV. These constraints are as high as 40 – 80% for the semi-major axes just larger than the FOV and steadily drop with increasing semi-major axis. The triangular shape of the additional HCI constraints is set by the mass limit in the HCI background-limited regime at the lower end and by the RV mass limit at the upper end. The power of HCI to add information for semi-major axes beyond the observed projected separations is amplified by the fact that higher inclinations are more likely given randomly-oriented orbital planes. For the star 40 Eri, HCI can even detect companions with semi-major axes up to 4 times the extent of the NACO FOV, depending on the planet’s orbital phase at the time of the observations, despite the smaller projected separations probed in the imaging data of this target.

5 Discussion

We have found that by combining constraints from NACO L*′* HCI with long-term RV monitoring data, the constraints on the planetary architecture of nearby stars are greatly improved, primarily for substellar masses of 5 – 20 $M_{\textrm{Jup}}$ and for semi-major axes similar to those of the ice giants in our own solar system (15 – 40 AU). The amount of improvement from including the HCI data depends primarily on the age of the star, with ages ¡5 Gyr leading to mass limits in the planetary mass regime. Differences in the RV data time baseline and scatter can also significantly change the constraints from these data and thus how much information HCI can add. For example, the difference between the contributions from HCI for HD 36395 and HD 42581 for the age of $\sim$3 Gyr are driven by the longer time baseline (7,240 days versus 6,220 days) and the lower RV scatter (4.4 versus 6.0 m/s) of HD 42581, resulting in HCI detecting up to 96% more companions for HD 36395 compared to up to 56% for HD 42581.

The fact that we did not detect any additional companions in the archival HCI is perhaps not surprising given the current statistics on giant planets and substellar companions at wide separations. HCI surveys have found that giant planets on wide orbits are relatively rare: planets of 5 – 13 $M_{\textrm{Jup}}$ and $a$ = 30 – 300 AU are found around less than $\sim$4% of FGK and M stars (Bowler 2016). Combined statistics on M stars from RV, HCI, and microlensing surveys yield similarly low giant planet occurrence rates of 1.5 – 3% for planets of 1 – 10 $M_{\textrm{Jup}}$ and $a$ = 10 – 1000 AU (Meyer et al. 2018). However, we note that the semi-major axes probed for the nearby stars in this sample are smaller than those probed by the typical HCI survey, since those surveys focus on young stars that are generally found in moving groups at around 30 – 50 pc compared to the ¡6 pc distances explored in this work. The intermediate semi-major axis range of 15 – 40 AU probed by HCI here is also beyond those semi-major axes well constrained by RV surveys, which give statistical constraints out to a few AU (e.g., Johnson et al. 2010; Bonfils et al. 2013). While this sample of 6 stars does not allow us to draw any significant conclusions about the statistics of substellar companions at these intermediate semi-major axes, it does highlight the new information gained by direct imaging searches around nearby stars both now and with future facilities, despite the generally older ages of field stars.

The HCI and RV limits reported here assume that all possible orbital planes are equally likely and thus use no prior information on the inclination of the possible planetary system. However $\tau$ Ceti, for example, hosts a debris disk that was found to have an inclination of 35 $\pm$ 10 deg by Herschel Space Observatory imaging (Lawler et al. 2014). Assuming this inclination for any possible planetary companions would reduce the strength of their RV signal significantly (57% of the RV signal from the same planet with an edge-on orbit) and would also increase the likelihood that a companion would fall at a wide separation in the single epoch of HCI data. This shows that for future HCI observations of other nearby stars, prior information on the inclination of the system such as debris disk inclination or stellar rotation axis could be taken into account to improve the information gained by adding HCI to existing long-term RV monitoring (Janson 2010).

In addition to being favorable targets for HCI and RV observations, nearby stars such as those in this sample will also be prime targets for astrometric searches for planets with Gaia (e.g., Perryman et al. 2014). With an astrometric precision on the order of tens of $\mu$as for V band magnitudes down to 12555https://www.cosmos.esa.int/web/gaia/science-performance gives the nominal astrometric performance of Gaia (site accessed on 17 April 2019)., Gaia will be especially sensitive to planets on wide orbits around the nearest stars for orbital periods well covered by the nominal Gaia mission time baseline of 5 years. Gaia will be less sensitive to planets on orbital periods longer than the mission time baseline, since it is then more challenging to disentangle the planet signal from the parallax and proper motion signals (see Ranalli et al. 2018). Therefore astrometry with Gaia will probe the parameter space of a ¡ 5 – 6 AU for masses down to a few Neptune masses for the very nearby stars in this sample. However, further analysis will be required to determine what planet masses will detectable for much longer orbital periods where HCI currently adds the most information. The astrometry for very bright stars (G ¡ 6 mag, $\tau$ Ceti and 40 Eri in this sample) in the Gaia DR2 is also currently limited by calibration issues from CCD saturation effects; however, the calibration is expected to improve in future data releases (Lindegren et al. 2018).

Recent works have also combined proper motion and/or position information from the Hipparcos and Gaia catalogues, leveraging the $\sim$24 year time baseline between these two missions to constrain the masses of known companions (Calissendorff & Janson 2018; Snellen & Brown 2018; Brandt et al. 2018; Dupuy et al. 2019) or search for new companions (Kervella et al. 2018). The latter work specifically targeted nearby stars (¡ 50 pc) to search for ”proper motion anomalies”: the difference between the long-term proper motion obtained from the positions in the Hipparcos and Gaia catalogues and the proper motion from either catalogue due to the presence of a companion. That work presents proper motion anomaly measurements for all stars in this sample except HD 36395, and found a significant (¿5$\sigma$) change in proper motion only for HD 42581, which is due to its brown dwarf companion HD 42581 B. We note that the error bars on the measured proper motion anomalies may be underestimated in that work given the results of Brandt (2018), which found that combining the Hippcarcos and Gaia catalogues in the same reference frame requires inflating the proper motion and position errors.

The upcoming JWST will be a powerful tool for HCI exoplanetary science thanks to the coronagraphs available in the NIRCam and MIRI instruments (Boccaletti et al. 2005; Beichman et al. 2010). We explored the mass limits that can be reached with JWST/NIRCam for the stars in this sample, specifically for the separations at which the sensitive background limit should be reached (approximately 3 – 5″ or 10 – 30 AU). Based on the evolutionary and atmospheric models presented in Linder et al. (2019) and assuming an effective integration time of $\sim$30 min, we find that planets with masses down to 1 – 2 $M_{\textrm{Jup}}$ would be detected by JWST for these nearby field stars, roughly an order of magnitude better than the mass limits presented here. With these mass limits, the JWST HCI data could detect Jupiter-mass planets currently hidden in the RV data, even for stars with very long time baselines and small overall RV scatter like $\tau$ Ceti.

In the near future, new HCI observations of these targets from the ground could also improve the imaging constraints presented in this work. The NACO data analyzed here have relatively short integration times of $\sim$10 – 25 min, and increasing the total integration time would improve the mass limits at wide separations by a few $M_{\textrm{Jup}}$ . The upcoming ERIS instrument for the VLT will provide improved 1 – 5 $\mu$m HCI capabilities compared to NACO when it sees first light in 2020 (Davies et al. 2018; Kenworthy et al. 2018). ERIS will use the AO Facility of the UT4 telescope with its deformable secondary mirror, resulting in fewer warm reflections before the light enters the cryogenic camera compared to NACO and improved sensitivity in the background limit. Finally, having a $\sim$12 years time baseline between these NACO observations and any new HCI observations would improve completeness due to the orbital motion of planets between the two epochs. Planets on certain orbits will have moved significantly on the sky, such that a planet initially too close to the host star would be further away and more easily detected in the later epoch. This effect is stronger for nearby stars like those in this sample because the HCI data can probe smaller semi-major axes and shorter orbital periods compared to more distant systems.

6 Summary

We quantified the limits on possible planetary architectures around 6 very nearby stars (d ¡ 6 pc) by combining limits on the mass and semi-major axes of possible companions from HCI and from RV measurements. We find that HCI adds information to the long-term RV monitoring constraints for intermediate semi-major axes (15 – 40 AU) and masses down to the background limit of the archival L*′* NACO imaging data reanalyzed here (5 – 20 $M_{\textrm{Jup}}$ ). The completeness was improved significantly for 4 of the 6 stars in the sample, and even by up to a few percent for $\tau$ Ceti, the star with the longest RV time baseline and lowest scatter, specifically for companions on nearly face-on orbits ($i$ ¡ 30 deg) and for cases where the RV time window included the time when the derivative of the RV signal is 0. These results show the power of combining constraints from different and complementary exoplanet detection methods. In the future, combining information from these the HCI and RV methods as well as from Gaia astrometry and JWST HCI observations will allow us to fully characterize our knowledge of the substellar companion population in the solar neighborhood at intermediate semi-major axes. The resulting information can be used to inform the target lists for future space missions like HabEx and LUVOIR, which will look for Earth analogues around the nearest stars.

Acknowledgements.

This work has been carried out within the framework of the National Center for Competence in Research PlanetS supported by the SNSF.