The Gemini Planet Imager Exoplanet Survey: Giant Planet and Brown Dwarf Demographics From 10-100 AU

TL;DR

This study analyzes the demographics of giant planets and brown dwarfs around stars at 10-100 au using GPIES data, revealing correlations with host star mass and distinct distribution patterns for planets and brown dwarfs.

Contribution

It provides the first statistical constraints on the occurrence and distribution of substellar companions at wide separations, highlighting differences between planets and brown dwarfs.

Findings

Higher planet occurrence around stars >1.5 M_sun

Power-law distribution for giant planets with α = -2.4, β = -2.0

Brown dwarfs are less common and have different mass and semi-major axis distributions

Abstract

We present a statistical analysis of the first 300 stars observed by the Gemini Planet Imager Exoplanet Survey (GPIES). This subsample includes six detected planets and three brown dwarfs; from these detections and our contrast curves we infer the underlying distributions of substellar companions with respect to their mass, semi-major axis, and host stellar mass. We uncover a strong correlation between planet occurrence rate and host star mass, with stars M $>$ 1.5 $M_\odot$ more likely to host planets with masses between 2-13 M$_{\rm Jup}$ and semi-major axes of 3-100 au at 99.92% confidence. We fit a double power-law model in planet mass (m) and semi-major axis (a) for planet populations around high-mass stars (M $>$ 1.5M$_\odot$) of the form $\frac{d^2 N}{dm da} \propto m^\alpha a^\beta$, finding $\alpha$ = -2.4 $\pm$ 0.8 and $\beta$ = -2.0 $\pm$ 0.5, and an integrated occurrence rate of $9^{+5}_{-4}$% between 5-13 M$_{\rm Jup}$ and 10-100 au. A significantly lower occurrence rate is obtained for brown dwarfs around all stars, with 0.8$^{+0.8}_{-0.5}$% of stars hosting a brown dwarf companion between 13-80 M$_{\rm Jup}$ and 10-100 au. Brown dwarfs also appear to be distributed differently in mass and semi-major axis compared to giant planets; whereas giant planets follow a bottom-heavy mass distribution and favor smaller semi-major axes, brown dwarfs exhibit just the opposite behaviors. Comparing to studies of short-period giant planets from the RV method, our results are consistent with a peak in occurrence of giant planets between ~1-10 au. We discuss how these trends, including the preference of giant planets for high-mass host stars, point to formation of giant planets by core/pebble accretion, and formation of brown dwarfs by gravitational instability.

Tables

Table 1. Table 1: Known Spectroscopic Binaries

Name	Type	$P$	$q$	$M_1$	$M_2$	Reference
		(days)		(M⊙)	(M⊙)
CC Eri	SB2	$1.56$	$0.54\pm 0.02$	${0.54}_{-0.03}^{+0.02}$	$0.29\pm 0.02$	Amado et al. (2000)
HR 784	SB1	$37.09$	$\cdots$	${1.21}_{-0.04}^{+0.02}$	${0.60}_{-0.13}^{+0.11}$	Gorynya & Tokovinin (2018)
$\zeta$ Hor	SB2	$12.93$	$0.88\pm 0.03$	${1.82}_{-0.13}^{+0.12}$	${1.58}_{-0.30}^{+0.04}$	Sahade & Hernández (1964)
HR 3395	SB1	$14.30$	$\cdots$	${1.24}_{-0.05}^{+0.05}$	${0.84}_{-0.10}^{+0.07}$	Abt & Willmarth (2006)
4 Sex	SB2	$3.05$	$0.945\pm 0.002$	${1.42}_{-0.02}^{+0.01}$	${1.35}_{-0.01}^{+0.01}$	Torres et al. (2003)
HD 142315	SB1	$1.264$	$\cdots$	${2.15}_{-0.06}^{+0.08}$	${2.01}_{-0.08}^{+0.06}$	Levato et al. (1987)
$\zeta$ TrA	SB1	$12.98$	$\cdots$	$1.13\pm 0.01$	${0.40}_{-0.03}^{+0.05}$	Skuljan et al. (2004)
V824 Ara	SB2	$1.68$	$0.909\pm 0.014$	$1.18\pm 0.01$	$1.08\pm 0.01$	Strassmeier & Rice (2000)
V4200 Sgr	SB2	$46.82$	$0.983\pm 0.001$	$0.88\pm 0.01$	$0.87\pm 0.01$	Fekel et al. (2017)
Peacock	SB1	$11.75$	$\cdots$	${5.31}_{-0.27}^{+0.29}$	${4.45}_{-0.63}^{+0.36}$	Luyten (1936)
$\iota$ Del	SB1	$11.04$	$\cdots$	${1.93}_{-0.24}^{+0.10}$	${0.87}_{-0.21}^{+0.76}$	Harper (1935)
42 Cap	SB2	$13.17$	$0.727\pm 0.005$	$1.94\pm 0.03$	$1.41\pm 0.02$	Fekel (1997)

Table 2. Table 2: Detected companion properties

Name	Epoch	$\rho$11Projected separation at the first detection in the GPIES campaign.	$a_{\mathrm{proj}}$11Projected separation at the first detection in the GPIES campaign.	$\Delta H$	Mass22Mass derived from absolute $H$ magnitude.	Astrometry reference
		(arc sec)	(au)	(mag)	(MJup)
51 Eri b	20141218	$0.451\pm 0.002$	$13.43\pm 0.09$	$14.4\pm 0.2$	$2.6\pm 0.3$	De Rosa et al. (2015)
$\beta$ Pic b	20151106	$0.249\pm 0.001$	$4.84\pm 0.02$	$9.3\pm 0.1$	$12.9\pm 0.5$	Wang et al. (2016)
HD 95086 b	20160229	$0.621\pm 0.005$	$53.65\pm 0.46$	$13.7\pm 0.2$	$2.6\pm 0.4$	Rameau et al. (2016)
HR 8799 c	20160919	$0.944\pm 0.001$	$38.99\pm 0.15$	$12.0\pm 0.1$	$8.3\pm 0.6$	Wang et al. (2018b)
HR 8799 d	20160919	$0.674\pm 0.001$	$27.85\pm 0.11$	$12.0\pm 0.1$	$8.3\pm 0.6$	Wang et al. (2018b)
HR 8799 e	20160919	$0.385\pm 0.002$	$15.89\pm 0.09$	$11.6\pm 0.1$	$9.2\pm 0.6$	Wang et al. (2018b)
HD 984 B	20150830	$0.216\pm 0.001$	$9.90\pm 0.04$	$6.5\pm 0.1$	$48\pm 3$	Johnson-Groh et al. (2017)
HR 2562 B	20160125	$0.619\pm 0.003$	$21.07\pm 0.11$	$11.7\pm 0.1$	${26}_{-13}^{+9}$	Konopacky et al. (2016)
PZ Tel B	20150730	$0.502\pm 0.001$	$23.65\pm 0.08$	$5.3\pm 0.1$	$64\pm 5$	this work

Table 3. Table 3: Constraints on substellar object populations from the GPIES 300-star sample, between 10-100 au. The power-law model has a number of substellar companions per star of f 𝑓 f within the given mass and semi-major axis limits, with power law indexes of α 𝛼 \alpha , β 𝛽 \beta , and γ 𝛾 \gamma for companion mass, semi-major axis, and stellar host mass, respectively.

	$f$	$\alpha$	$\beta$	$\gamma$
Planets around high-mass stars: 5 MJup $m$ 13 MJup, $1.5$MMM⊙
Median	8.9%	-2.37	-1.99	–
68% CI	5.3 – 13.9%	-3.16 – -1.51	-2.47 – -1.44	–
95% CI	2.9 – 20.6%	-3.90 – -0.56	-2.92 – -0.81	–
Planets around all stars: 5 MJup $m$ 13 MJup, $0.2$MMM⊙
Median	3.5%	-2.277	-1.68	2.03
68% CI	2.1 – 5.4%	-2.98 – -1.48	-2.18 – -1.10	1.05 – 3.07
95% CI	1.1 – 8.0%	-3.65 – -0.54	-2.63 – -0.38	0.09 – 4.16
Brown dwarfs around all stars: 13 MJup $m$ 80 MJup, $0.2$MMM⊙
Median	0.8%	-0.47	-0.65	-0.85
68% CI	0.3 – 1.6%	-1.53 – 0.67	-1.69 – 0.89	-2.33 – 0.59
95% CI	0.1 – 2.8%	-2.60 – 1.96	-2.52 – 3.86	-3.76 – 2.01
Substellar objects around all stars: 5 MJup $m$ 80 MJup, $0.2$MMM⊙
Median	4.7%	-2.09	-1.44	1.15
68% CI	3.2 – 6.7%	-2.49 – -1.71	-1.92 – -0.91	0.35 – 1.96
95% CI	2.0 – 9.1%	-2.89 – -1.33	-2.36 – -0.28	-0.48 – 2.75

Table 4. Table 4: 300-star sample properties and observing log

Name	R.A.	Dec.	SpT	ref	$d$	ref	$I$	$H$	Moving	Age	ref	Mass	UT Date	$t_{\exp }\times n_{\mathrm{coadd}}$	$t_{\mathrm{int}}$	$\Delta$PA
									Group
	hh:mm:ss.ss	dd:mm:ss.s			pc		mag	mag		Myr		M⊙		sec	min	deg
BU Psc	00:02:02.59	-02:45:58.2	A8/9V	1	$64.70\pm 0.38$	2	6.5	6.3	$\cdots$	554	3	1.49	20160920	$59.6\times 1$	34.8	21.3
HD 105	00:05:52.54	-41:45:11.0	G0V	4	$38.845\pm 0.084$	2	6.9	6.2	TucHor	45	5	1.14	20141110	$59.6\times 1$	38.8	31.0
HR 9	00:06:50.08	-23:06:27.1	F3V	4	$39.960\pm 0.098$	2	5.8	5.3	$\beta$ Pic	26	6	1.49	20141109	$59.6\times 1$	34.8	44.4
HD 987	00:13:52.83	-74:41:17.5	G8V	4	$45.847\pm 0.051$	2	8.0	7.1	TucHor	45	5	0.94	20151105	$59.6\times 1$	28.8	16.6
HD 984	00:14:10.25	-07:11:56.8	F7V	1	$45.91\pm 0.12$	2	6.8	6.2	Columba	42	5	1.28	20150830	$59.6\times 1$	21.9	15.1
													20160919	$59.6\times 1$	13.9	9.5
HD 1237 A	00:16:12.70	-79:51:04.0	G8V	4	$17.5649\pm 0.0073$	2	5.9	5.0	Carina-Near	200	7	1.01	20150831	$59.6\times 1$	34.8	13.6
HD 1466	00:18:26.10	-63:28:39.0	F8V	4	$42.970\pm 0.055$	2	6.9	6.2	TucHor	45	5	1.19	20150703	$59.6\times 1$	42.7	21.0
9 Cet	00:22:51.79	-12:12:34.0	G2.5V	8	$21.310\pm 0.023$	2	5.8	5.0	$\cdots$	300	9	1.14	20150901	$59.6\times 1$	36.8	31.2
$\kappa$ Phe	00:26:12.12	-43:40:47.7	A5IVn	10	$23.809\pm 0.085$	11	3.7	3.6	$\cdots$	578	3	1.83	20150702	$52.4\times 1$	41.0	47.5
${\lambda }^1$ Phe	00:31:24.98	-48:48:12.7	A1Va	12	$52.97\pm 0.67$	11	4.5	4.8	$\cdots$	473	3	2.18	20150704	$59.6\times 1$	41.8	35.3
${\beta }^1$ Tuc	00:31:32.67	-62:57:29.6	B8/A0	13	$41.41\pm 0.34$	11	4.4	4.7	TucHor	45	5	2.44	20151106	$59.6\times 1$	37.8	18.6
HD 3221	00:34:51.20	-61:54:58.0	K4Ve	4	$44.499\pm 0.050$	2	8.4	6.7	TucHor	45	5	1.06	20160918	$59.6\times 1$	35.8	19.2
HD 3375	00:36:00.80	-59:43:01.9	F2V	13	$57.099\pm 0.095$	2	6.4	6.1	$\cdots$	30	9	1.51	20151105	$59.6\times 1$	33.8	22.5
HD 3821 A	00:40:47.48	-07:13:56.6	G1V	4	$26.496\pm 0.033$	2	6.4	5.6	$\cdots$	300	9	1.05	20160921	$59.6\times 1$	34.8	24.5
HD 3888	00:40:51.58	-53:12:35.7	F7V	13	$44.177\pm 0.060$	2	6.8	6.2	$\cdots$	70	9	1.24	20150801	$59.6\times 1$	12.9	8.6
													20150831	$59.6\times 1$	36.8	25.6
ksi Phe	00:41:46.30	-56:30:05.2	A3VpSrCr Ksn	14	$68.43\pm 0.42$	2	5.5	5.6	$\cdots$	767	3	1.91	20150729	$59.6\times 1$	17.9	12.0
HD 5294	00:54:59.27	+24:06:01.0	G5	15	$29.917\pm 0.045$	2	6.7	5.9	$\cdots$	300	9	1.03	20160921	$59.6\times 1$	35.8	12.3
HIP 6276	01:20:32.27	-11:28:03.7	G9V	4	$35.332\pm 0.060$	2	7.7	6.7	AB Dor	149	5	0.94	20160920	$59.6\times 1$	35.8	30.4
HD 8558	01:23:21.14	-57:28:50.5	G7V	4	$45.314\pm 0.053$	2	7.8	6.9	TucHor	45	5	1.01	20160919	$59.6\times 1$	37.8	23.8
CC Phe	01:28:08.66	-52:38:19.1	K1V	4	$39.844\pm 0.044$	2	8.3	6.9	TucHor	45	5	0.86	20160918	$59.6\times 1$	35.8	24.9
48 Cet	01:29:36.10	-21:37:45.5	A0V	16	$71.38\pm 1.48$	11	4.8	5.0	$\cdots$	553	3	2.27	20141112	$59.6\times 1$	39.8	63.3
49 Cet	01:34:37.72	-15:40:34.9	A1V	16	$57.07\pm 0.33$	2	5.4	5.5	Argus	50	17	2.02	20141111	$59.6\times 1$	16.9	17.4
EX Cet	01:37:35.37	-06:45:36.7	K0/1V	1	$24.044\pm 0.025$	2	6.9	5.9	TucHor	45	5	0.88	20141215	$59.6\times 1$	6.0	3.6
													20150901	$59.6\times 1$	30.8	24.5
HD 10472	01:40:24.06	-60:59:56.6	F2IV/V	13	$71.15\pm 0.15$	2	7.2	6.7	TucHor	45	5	1.44	20151106	$59.6\times 1$	37.8	19.5
DK Cet	01:57:48.91	-21:54:05.0	G4V	4	$41.411\pm 0.092$	2	7.4	6.6	TucHor	45	5	1.08	20141109	$59.6\times 1$	34.8	10.1
HD 13183	02:07:17.98	-53:11:56.4	G7V	4	$49.535\pm 0.081$	2	7.9	7.0	TucHor	45	5	1.05	20141215	$59.6\times 1$	38.8	33.8
HD 13246	02:07:26.02	-59:40:45.8	F7V	4	$45.610\pm 0.056$	2	7.0	6.3	TucHor	45	5	1.23	20141218	$59.6\times 1$	38.8	27.4
$\varphi$ Eri	02:16:30.60	-51:30:44.0	B8IV	18	$47.12\pm 0.27$	11	3.5	4.0	TucHor	45	5	3.58	20150831	$59.6\times 1$	37.8	27.0
HIP 11152	02:23:26.64	+22:44:06.7	M1V	19	$27.193\pm 0.045$	2	9.3	7.6	$\beta$ Pic	26	6	0.50	20160918	$59.6\times 1$	37.8	12.5
HD 15115	02:26:16.20	+06:17:33.6	F4IV	20	$49.00\pm 0.10$	2	6.3	5.9	Columba	42	5	1.45	20141216	$59.6\times 1$	62.6	35.2
CC Eri	02:34:22.60	-43:47:47.0	K7V	4	$11.5373\pm 0.0049$	2	7.1	5.1	$\cdots$	30	9	0.84†	20141109	$59.6\times 1$	34.8	32.0
HD 16699 A	02:38:44.29	-52:57:03.1	F8V	4	$57.876\pm 0.093$	2	7.2	6.7	$\cdots$	27	9	1.28	20160919	$59.6\times 1$	33.8	22.6
$ϵ$ Hyi	02:39:35.23	-68:16:00.8	B9Va	18	$46.55\pm 0.20$	11	3.9	4.4	TucHor	45	5	2.70	20150729	$59.6\times 1$	47.7	28.5
HR 789	02:39:47.93	-42:53:29.8	A1Vb	10	$40.47\pm 0.64$	2	4.5	4.6	$\beta$ Pic	26	6	2.17	20141110	$59.6\times 1$	39.8	43.6
HR 784	02:40:12.42	-09:27:10.4	F8VFe-0.4	10	$21.93\pm 0.12$	2	5.3	4.7	$\cdots$	125	9	1.81†	20141214	$59.6\times 1$	40.8	32.7
$\zeta$ Hor	02:40:39.60	-54:33:00.0	F6V	4	$48.94\pm 0.28$	2	4.7	4.3	$\cdots$	30	9	3.34†	20150831	$59.6\times 1$	21.9	15.7
84 Cet B	02:41:13.78	-00:41:42.0	–	–	$22.59\pm 0.43$	11	6.6	4.5	$\cdots$	124	9	0.73	20160920	$59.6\times 1$	37.8	25.5
$\iota$ Hor	02:42:33.47	-50:48:01.1	F8V	4	$17.327\pm 0.016$	2	4.8	4.3	$\cdots$	300	9	1.24	20141111	$59.6\times 1$	39.8	29.6
HD 17207	02:45:01.13	-22:09:58.7	F7V	10	$38.860\pm 0.069$	2	6.7	6.0	$\cdots$	125	9	1.20	20141112	$59.6\times 1$	38.8	58.9
HR 826	02:46:45.10	-21:38:22.2	F3IV/V	16	$48.198\pm 0.084$	2	6.1	5.6	$\cdots$	300	9	1.53	20141109	$59.6\times 1$	30.8	69.5
HD 17332 B	02:47:27.21	+19:22:21.2	G5V	21	$33.065\pm 0.056$	2	7.0	5.7	AB Dor	149	5	0.97	20150830	$59.6\times 1$	35.8	14.4
VZ Hor	02:51:53.00	-61:37:05.0	K1V	4	$40.126\pm 0.048$	2	7.9	6.7	$\cdots$	300	9	0.99	20141216	$59.6\times 1$	38.8	18.6
EP Eri	02:52:32.13	-12:46:11.0	K1V	16	$10.3584\pm 0.0047$	2	5.1	4.2	$\cdots$	126	9	0.92	20160920	$59.6\times 1$	37.8	26.8
HD 18809	03:00:19.70	-37:27:16.0	G4V	4	$48.89\pm 0.20$	2	7.8	6.9	$\cdots$	125	9	1.05	20141111	$59.6\times 1$	39.8	58.4
HD 20232	03:13:50.19	-38:48:33.4	A2/3III/IV	22	$79.47\pm 0.19$	2	6.5	6.4	$\cdots$	429	3	1.66	20141215	$59.6\times 1$	39.8	53.9
${\kappa }^1$ Cet	03:19:21.70	+03:22:12.7	G5V	8	$9.140\pm 0.023$	11	4.2	3.0	$\cdots$	300	9	1.05	20141218	$59.6\times 1$	38.8	20.5
$ϵ$ Eri	03:32:55.84	-09:27:29.7	K2V	8	$3.2161\pm 0.0017$	11	2.8	1.9	$\cdots$	500	9	0.86	20141110	$14.5\times 4$	35.0	29.2
HR 1139	03:43:33.84	-10:29:08.2	A5III-IV	14	$41.81\pm 0.22$	2	5.3	5.2	$\cdots$	534	3	1.62	20160918	$59.6\times 1$	37.8	22.2
HR 1189	03:48:35.48	-37:37:19.2	A1Va	12	$52.64\pm 0.33$	2	5.3	5.4	TucHor	45	5	2.01	20151103	$59.6\times 1$	37.8	60.4
HR 1190	03:48:35.88	-37:37:12.6	B9.5Van	12	$52.70\pm 0.48$	2	4.6	4.6	TucHor	45	5	2.53	20151103	$59.6\times 1$	36.8	38.7
HD 24916	03:57:28.70	-01:09:34.1	K4V	1	$15.301\pm 0.011$	2	7.0	5.5	$\cdots$	125	9	0.80	20141214	$59.6\times 1$	37.8	21.8
HR 1249	04:02:36.74	-00:16:08.1	F7/8V	1	$18.771\pm 0.041$	2	4.9	4.3	AB Dor	149	5	1.28	20141217	$59.6\times 1$	32.8	18.9
39 Tau	04:05:20.26	+22:00:32.0	G5V	23	$16.895\pm 0.020$	2	5.3	4.5	$\cdots$	300	9	1.11	20150829	$59.6\times 1$	37.8	12.6
HD 25953	04:06:41.52	+01:41:02.9	F6V	1	$56.97\pm 0.15$	2	7.3	6.7	AB Dor	149	5	1.25	20151221	$59.6\times 1$	36.8	17.6
$\gamma$ Dor	04:16:01.59	-51:29:11.9	F1V	10	$20.46\pm 0.15$	11	3.8	3.5	$\cdots$	894	3	1.56	20141110	$59.6\times 1$	37.8	28.6
HD 28287	04:26:38.57	-28:57:06.5	G9.5Vk:	10	$35.645\pm 0.039$	2	8.0	6.9	AB Dor	149	5	0.89	20141112	$59.6\times 1$	38.8	145.3
													20141214	$59.6\times 1$	9.9	149.3
51 Eri	04:37:36.13	-02:28:24.8	F0IV	14	$29.78\pm 0.12$	2	4.8	4.8	$\beta$ Pic	26	6	1.58	20141218	$59.6\times 1$	37.8	23.8
													20150131	$59.6\times 1$	62.6	36.5
													20150831	$59.6\times 1$	44.7	29.9
V834 Tau	04:41:18.86	+20:54:05.4	K4V	24	$13.2381\pm 0.0100$	2	6.7	5.3	$\cdots$	82	9	0.71	20160919	$59.6\times 1$	37.8	12.7
HD 30447	04:46:49.50	-26:18:09.0	F3V	22	$80.54\pm 0.25$	2	7.4	7.0	Columba	42	5	1.45	20160920	$59.6\times 1$	31.8	93.0
HD 31392	04:54:04.20	-35:24:16.2	G9V	10	$25.765\pm 0.020$	2	6.9	5.9	$\cdots$	200	9	0.99	20141216	$59.6\times 1$	40.8	62.1
HR 1597	04:54:53.04	-58:32:51.5	F5V	10	$30.468\pm 0.028$	2	5.7	5.1	$\cdots$	125	9	1.33	20160125	$59.6\times 1$	36.8	21.4
${\pi }^1$ Ori	04:54:53.71	+10:09:04.0	A0Va lB	12	$35.66\pm 0.32$	11	4.4	4.5	$\cdots$	424	3	1.85	20151221	$59.6\times 1$	23.9	9.9
HIP 23309	05:00:47.13	-57:15:25.5	M0Ve	4	$26.901\pm 0.020$	2	8.1	6.4	$\beta$ Pic	26	6	0.79	20151106	$59.6\times 1$	37.8	23.1
$\zeta$ Dor	05:05:30.69	-57:28:22.7	F9VFe-0.5	10	$11.646\pm 0.024$	11	4.1	3.4	$\cdots$	300	9	1.18	20151105	$18.9\times 3$	35.3	20.8
AS Col	05:20:38.02	-39:45:18.0	F6V	4	$47.663\pm 0.073$	2	6.9	6.2	Columba	42	5	1.28	20141111	$59.6\times 1$	38.8	43.6
HD 35650	05:24:30.14	-38:58:10.2	K6V	4	$17.4793\pm 0.0086$	2	7.9	6.1	AB Dor	149	5	0.71	20141110	$59.6\times 1$	39.8	24.7
HIP 25434	05:26:22.98	-43:22:36.4	G0	25	$92.41\pm 0.26$	2	8.3	7.9	Columba	42	5	1.23	20151221	$59.6\times 1$	36.8	45.0
HD 35841	05:26:36.59	-22:29:23.7	F3V	16	$103.68\pm 0.30$	2	8.4	7.8	Columba	42	5	1.30	20160228	$59.6\times 1$	46.7	46.1
AF Lep	05:27:04.76	-11:54:03.5	F8V(n)k:	10	$26.877\pm 0.019$	2	5.7	5.1	$\beta$ Pic	26	6	1.15	20141214	$59.6\times 1$	28.8	25.8
HD 36435	05:27:39.35	-60:24:57.6	G5V	4	$19.478\pm 0.012$	2	6.3	5.3	$\cdots$	125	9	0.97	20141217	$59.6\times 1$	39.8	22.0
HD 36329	05:29:24.10	-34:30:56.0	G3V	4	$72.46\pm 0.15$	2	7.8	7.1	Columba	42	5	1.35	20141111	$59.6\times 1$	35.8	40.4
$\zeta$ Lep	05:46:57.35	-14:49:19.0	A2IV-V(n)	10	$21.608\pm 0.075$	11	3.3	3.3	$\beta$ Pic	26	6	2.06	20160229	$29.1\times 2$	37.8	37.5
$\beta$ Pic	05:47:17.08	-51:04:00.2	A6V	10	$19.440\pm 0.045$	11	3.6	3.5	$\beta$ Pic	26	6	1.73	20151106	$59.6\times 1$	36.8	28.0
HD 40781	05:59:13.38	-38:42:38.1	G0V	4	$64.19\pm 0.71$	2	7.7	7.0	$\cdots$	127	9	1.21	20160227	$59.6\times 1$	37.8	60.4
HD 41278	06:00:55.34	-54:57:05.0	F5V	13	$57.201\pm 0.083$	2	6.9	6.4	AB Dor	149	5	1.37	20151217	$59.6\times 1$	36.8	27.0
HR 2234	06:15:29.66	-04:54:52.7	A3Vas	26	$73.28\pm 0.45$	2	5.8	5.9	$\cdots$	440	3	1.91	20150130	$59.6\times 1$	40.8	26.3
AO Men	06:18:28.21	-72:02:41.5	K4Ve	4	$39.263\pm 0.039$	2	8.4	7.0	$\beta$ Pic	26	6	0.89	20151221	$59.6\times 1$	35.8	16.6
V1358 Ori	06:19:08.04	-03:26:20.0	G0V	1	$52.01\pm 0.15$	2	7.3	6.6	Columba	42	5	1.24	20150129	$59.6\times 1$	29.8	18.1
													20151219	$59.6\times 1$	30.8	18.4
AB Pic	06:19:12.89	-58:03:15.8	K1V(e)	4	$50.120\pm 0.073$	2	8.3	7.1	Carina	45	5	1.00	20151106	$59.6\times 1$	35.8	19.5
V435 Car	06:21:50.05	-51:14:15.8	A5V	27	$71.43\pm 0.17$	2	6.4	6.3	$\cdots$	439	3	1.64	20141217	$59.6\times 1$	39.8	29.9
HD 45270	06:22:30.96	-60:13:07.7	G1V	4	$23.890\pm 0.014$	2	5.9	5.2	AB Dor	149	5	1.15	20141216	$59.6\times 1$	39.8	22.7
HIP 31878	06:39:50.06	-61:28:42.2	K7V(e)	4	$22.241\pm 0.014$	2	8.4	6.6	AB Dor	149	5	0.69	20160128	$59.6\times 1$	33.8	17.4
HD 48524	06:42:05.52	-38:00:13.7	G0V	22	$67.84\pm 0.49$	2	7.8	7.1	$\cdots$	125	9	1.17	20141214	$59.6\times 1$	33.8	74.7
26 Gem	06:42:24.33	+17:38:43.1	A2V	15	$43.63\pm 1.28$	11	5.1	5.1	Columba	42	5	1.87	20160226	$59.6\times 1$	37.8	14.0
HD 48370	06:43:01.02	-02:53:19.3	K0V(+G)	1	$36.072\pm 0.068$	2	7.2	6.4	$\cdots$	125	9	1.06	20150201	$59.6\times 1$	39.8	24.2
HD 49855	06:43:46.32	-71:58:35.8	G6V	4	$55.59\pm 0.10$	2	8.2	7.4	Carina	45	5	0.93	20160320	$59.6\times 1$	36.8	15.3
HR 2562	06:50:01.01	-60:14:57.9	F5VFe+0.4	10	$34.041\pm 0.048$	2	5.7	5.1	$\cdots$	300	9	1.43	20160125	$59.6\times 1$	33.8	19.4
HR 2558	06:51:42.45	-36:13:48.5	A3V	22	$55.41\pm 0.17$	2	5.6	5.5	$\cdots$	687	3	1.69	20151217	$59.6\times 1$	35.8	74.6
													20160316	$59.6\times 1$	56.7	99.4
HD 51797	06:56:23.55	-46:46:55.1	K0V(e)	4	$95.65\pm 0.24$	2	8.9	7.9	Columba	42	5	1.20	20160321	$59.6\times 1$	30.8	30.9
HD 53143	06:59:59.66	-61:20:10.2	G9V	4	$18.3602\pm 0.0079$	2	6.1	5.1	IC2391	50	17	1.00	20151219	$59.6\times 1$	36.8	20.3
HR 2720	07:12:04.08	-30:49:17.0	A8III/IV	22	$47.423\pm 0.087$	2	5.7	5.5	$\cdots$	724	3	1.54	20160327	$59.6\times 1$	36.8	170.6
$\lambda$ Gem A	07:18:05.61	+16:32:25.7	A4IV	26	$30.93\pm 0.19$	11	3.4	3.5	$\cdots$	638	3	2.12	20160319	$29.1\times 2$	34.9	14.6
HR 2813	07:20:21.40	-52:18:41.0	F5V	10	$33.15\pm 0.24$	2	5.5	5.1	Carina-Near	200	7	1.43	20160320	$59.6\times 1$	35.8	27.0
HR 2814	07:20:21.86	-52:18:33.3	F9V	10	$36.65\pm 0.59$	2	5.9	5.3	Carina-Near	200	7	1.35	20141217	$59.6\times 1$	38.8	28.8
HD 57969	07:20:22.98	-56:17:41.3	A1V	13	$72.20\pm 0.17$	2	6.3	6.4	$\cdots$	217	3	1.71	20151221	$59.6\times 1$	31.8	32.3
HD 58192	07:23:53.76	-17:24:48.2	F7V	16	$36.594\pm 0.053$	2	6.3	5.7	$\cdots$	300	9	1.24	20160324	$59.6\times 1$	29.8	48.2
HD 58556	07:25:57.22	-02:14:54.5	G1V	1	$33.925\pm 0.052$	2	6.4	5.7	$\cdots$	300	9	1.19	20160328	$59.6\times 1$	36.8	21.1
HD 59704	07:29:31.40	-38:07:21.0	F7V	22	$53.370\pm 0.097$	2	7.2	6.5	Carina-Near	200	7	1.26	20160326	$59.6\times 1$	22.9	62.5
HR 2882	07:30:42.57	-37:20:22.1	G3V	4	$21.770\pm 0.013$	2	6.0	5.3	$\cdots$	200	9	1.06	20141111	$59.6\times 1$	38.8	76.3
HD 61033	07:34:28.03	-52:58:05.4	G7V(e)	4	$28.93\pm 0.38$	11	6.8	5.8	$\cdots$	300	9	1.05	20160324	$59.6\times 1$	36.8	22.9
HD 61005	07:35:47.50	-32:12:14.7	G8Vk	10	$36.486\pm 0.043$	2	7.5	6.6	Argus	50	17	0.98	20150130	$59.6\times 1$	40.8	109.5
HIP 37288	07:39:23.04	+02:11:01.2	M0.0Ve	28	$15.058\pm 0.012$	2	8.1	6.1	$\cdots$	300	9	0.67	20141218	$59.6\times 1$	38.8	22.6
HD 62850	07:42:36.10	-59:17:51.0	G3V	4	$32.912\pm 0.028$	2	6.6	5.9	Carina-Near	200	7	1.15	20160318	$59.6\times 1$	37.8	19.5
HD 63008	07:44:12.51	-50:27:24.2	F9VFe-0.5	10	$30.084\pm 0.022$	2	6.1	5.4	$\cdots$	300	9	1.19	20160322	$59.6\times 1$	35.8	28.0
HD 63008 B	07:44:16.48	-50:27:59.8	G5V	10	$30.080\pm 0.027$	2	6.9	6.0	$\cdots$	300	9	1.01	20160322	$59.6\times 1$	35.8	27.2
YZ CMi	07:44:40.17	+03:33:08.8	M4.0Ve	29	$5.9874\pm 0.0021$	2	8.3	6.0	$\cdots$	70	9	0.29	20160225	$59.6\times 1$	32.8	22.5
HD 62883	07:45:50.88	-07:31:46.4	F5V	1	$58.23\pm 0.23$	2	6.8	6.2	$\cdots$	30	9	1.26	20150201	$59.6\times 1$	40.8	43.3
HD 63608	07:46:17.00	-59:48:34.0	K0IV-V(k)	10	$33.806\pm 0.030$	2	7.5	6.5	Carina-Near	200	7	0.97	20160125	$59.6\times 1$	32.8	21.0
HR 3070	07:49:12.90	-60:17:01.0	F4V	10	$33.767\pm 0.074$	2	5.2	4.9	Carina-Near	200	7	1.56	20150408	$59.6\times 1$	41.8	22.4
HD 64122	07:49:29.35	-54:54:04.4	F6/7V	13	$80.66\pm 0.19$	2	7.9	7.3	$\cdots$	30	9	1.32	20141214	$59.6\times 1$	32.8	25.1
TYC 206-419-1	08:07:09.09	+07:23:00.1	K8	15	$41.469\pm 0.068$	2	8.9	7.3	AB Dor	149	5	0.82	20160319	$59.6\times 1$	26.8	12.2
HD 68143	08:11:15.15	+01:16:36.4	F3/5+F6/7	1	$80.15\pm 1.37$	2	8.1	7.4	$\cdots$	125	9	1.25	20160226	$59.6\times 1$	34.8	19.2
q Pup	08:18:33.31	-36:39:33.4	A8V	10	$28.63\pm 0.15$	11	4.2	4.1	$\cdots$	894	3	1.71	20160321	$59.6\times 1$	37.8	76.6
HD 70703	08:21:00.48	-52:13:40.9	A0/1V	27	$72.71\pm 0.14$	2	6.3	6.3	Carina	45	5	1.74	20151217	$59.6\times 1$	34.8	30.9
V478 Hya	08:22:49.98	+01:51:34.0	G1/2V	1	$59.28\pm 0.16$	2	8.1	7.3	$\cdots$	200	9	1.08	20141216	$59.6\times 1$	13.9	6.8
													20151220	$59.6\times 1$	9.9	8.7
													20160328	$59.6\times 1$	35.8	21.7
HR 3300	08:22:55.18	-52:07:25.6	A0V	27	$73.25\pm 0.38$	2	5.7	5.9	Carina	45	5	2.21	20150407	$59.6\times 1$	41.8	29.7
V592 Pup	08:25:17.72	-34:22:01.2	G1IV	4	$66.98\pm 0.15$	2	7.1	6.4	$\cdots$	80	9	1.34	20160227	$59.6\times 1$	44.7	86.6
HR 3341	08:26:25.24	-52:48:27.0	A0V	27	$69.33\pm 0.27$	2	5.8	5.9	AB Dor	149	5	1.96	20151221	$59.6\times 1$	37.8	23.1
HD 72687	08:33:15.36	-29:57:23.8	G5V	4	$45.436\pm 0.085$	2	7.6	6.8	$\cdots$	125	9	1.05	20160229	$59.6\times 1$	35.8	172.2
HR 3395	08:35:50.98	+06:37:12.8	F8V	30	$26.280\pm 0.090$	2	5.4	4.8	$\cdots$	125	9	2.08†	20150404	$59.6\times 1$	32.8	19.5
V401 Hya	08:37:50.46	-06:48:25.2	G8/K0(IV)	1	$24.342\pm 0.026$	2	6.1	5.3	$\cdots$	300	9	1.11	20150406	$59.6\times 1$	36.8	29.0
HD 74341 A	08:40:43.54	-57:32:43.6	A4V	31	$72.63\pm 0.19$	2	6.4	5.9	$\cdots$	320	3	1.72	20151219	$59.6\times 1$	37.8	24.6
$\eta$ Cha	08:41:19.51	-78:57:48.1	B8V	13	$94.97\pm 1.44$	11	5.5	5.7	eta Cha	11	5	3.20	20151218	$59.6\times 1$	37.8	14.7
Asellus Borealis	08:43:17.21	+21:28:06.9	A1IV	15	$55.56\pm 0.65$	11	4.5	4.8	$\cdots$	527	3	2.22	20160301	$59.6\times 1$	36.8	14.7
HD 74576	08:43:18.00	-38:52:57.0	K3V	4	$11.1860\pm 0.0032$	2	5.7	4.4	$\cdots$	125	9	0.86	20141214	$59.6\times 1$	38.8	58.1
HD 75519	08:49:05.70	-39:57:16.0	G8V	4	$36.509\pm 0.049$	2	7.2	6.3	$\cdots$	125	9	1.06	20141215	$59.6\times 1$	19.9	38.8
HD 76378	08:54:57.24	-24:23:39.4	K2VCN+1	10	$33.684\pm 0.042$	2	7.8	6.6	$\cdots$	125	9	0.96	20150129	$59.6\times 1$	40.8	65.4
HD 77056	08:59:18.26	-24:43:43.9	G8V	16	$37.881\pm 0.063$	2	7.6	6.7	$\cdots$	300	9	0.97	20150501	$59.6\times 1$	18.9	52.7
V405 Hya	09:04:20.76	-15:54:51.0	K2V	16	$27.382\pm 0.039$	2	7.8	6.5	$\cdots$	125	9	0.85	20150130	$59.6\times 1$	40.8	42.2
$\beta$ Car	09:13:11.98	-69:43:01.9	A1III-	10	$34.70\pm 0.13$	11	1.2	1.5	$\cdots$	324	3	3.31	20160319	$2.9\times 10$	18.2	22.0
HD 80846	09:22:17.72	+06:17:51.3	G0	15	$71.90\pm 0.32$	2	7.9	7.3	$\cdots$	125	9	1.23	20160225	$59.6\times 1$	30.8	37.7
HD 82282	09:27:57.54	-66:06:07.7	K0V	4	$38.209\pm 0.036$	2	8.1	6.9	$\cdots$	300	9	0.89	20150504	$59.6\times 1$	41.8	19.1
HIP 46634	09:30:35.04	+10:36:00.0	–	–	$35.364\pm 0.094$	2	7.9	6.8	$\cdots$	125	9	0.92	20150201	$59.6\times 1$	40.8	17.1
LQ Hya	09:32:25.57	-11:11:04.7	K1Vp	1	$18.290\pm 0.022$	2	6.9	5.6	$\cdots$	70	9	0.82	20141218	$59.6\times 1$	37.8	36.8
HD 82943	09:34:50.73	-12:07:46.3	F9VFe+0.5	10	$27.612\pm 0.029$	2	5.9	5.2	$\cdots$	500	9	1.23	20160226	$59.6\times 1$	33.8	28.6
HD 84075	09:36:18.03	-78:20:42.0	G2V	13	$64.098\pm 0.091$	2	7.9	7.2	Argus	50	17	1.16	20160228	$59.6\times 1$	36.8	14.2
HD 83946	09:38:54.10	-64:59:26.7	F5V	13	$90.55\pm 0.22$	2	8.2	7.7	Volans-Carina	80	2	1.25	20160228	$59.6\times 1$	36.8	18.3
HD 84040	09:41:41.92	-25:57:51.8	F7V	16	$71.61\pm 0.24$	2	8.0	7.3	$\cdots$	70	9	1.20	20160229	$59.6\times 1$	35.8	107.0
HD 84330 B	09:42:40.61	-55:49:53.5	G3V	4	$43.700\pm 0.071$	2	7.3	6.0	$\cdots$	122	9	1.10	20151218	$59.6\times 1$	37.8	24.2
													20160318	$59.6\times 1$	18.9	13.1
HD 84330 A	09:42:40.80	-55:49:55.0	G1V	4	$43.695\pm 0.072$	2	6.6	5.9	$\cdots$	300	9	1.15	20150404	$59.6\times 1$	33.8	25.6
HD 84273	09:43:20.34	-29:48:14.7	G5V	22	$45.235\pm 0.094$	2	7.6	6.9	$\cdots$	300	9	1.00	20160324	$59.6\times 1$	36.8	168.9
4 Sex	09:50:30.08	+04:20:37.1	F5V	30	$47.21\pm 0.35$	2	5.8	5.2	$\cdots$	300	9	2.69†	20150501	$59.6\times 1$	31.8	17.3
6 Sex	09:51:14.02	-04:14:35.8	A4V	1	$65.24\pm 0.44$	2	5.7	5.7	$\cdots$	710	3	1.79	20160316	$59.6\times 1$	35.8	22.4
HD 87099	10:02:00.46	-34:10:25.0	A8III	22	$76.37\pm 0.27$	2	6.7	6.4	$\cdots$	925	3	1.54	20150408	$59.6\times 1$	39.8	117.0
AN Sex	10:12:17.67	-03:44:44.4	M2V	32	$7.7048\pm 0.0050$	2	7.0	5.3	AB Dor	149	5	0.55	20160316	$59.6\times 1$	32.8	21.9
HD 298936	10:13:14.78	-52:30:54.0	K3Ve	4	$53.350\pm 0.091$	2	8.6	7.4	Carina	45	5	0.97	20160430	$59.6\times 1$	34.8	24.5
HR 4013	10:13:24.73	-33:01:54.2	G0V	10	$22.442\pm 0.016$	2	5.8	5.1	$\cdots$	300	9	1.12	20160320	$59.6\times 1$	34.8	126.0
q Vel	10:14:44.16	-42:07:18.9	A2Va	10	$31.08\pm 0.15$	11	3.6	3.7	Argus	50	17	2.23	20151219	$59.6\times 1$	7.0	6.2
													20160319	$52.4\times 1$	37.5	49.7
HR 4068	10:21:07.92	-17:59:05.6	F3IV	16	$62.58\pm 0.22$	2	6.0	5.6	$\cdots$	70	9	1.57	20160322	$59.6\times 1$	35.8	30.7
HD 90712	10:27:47.78	-34:23:58.1	G0VCH-0.3	10	$37.892\pm 0.071$	2	6.9	6.1	$\cdots$	300	9	1.13	20160321	$59.6\times 1$	35.8	48.7
HD 90885 B	10:28:17.90	-52:33:40.0	K1V	4	$33.373\pm 0.049$	2	7.8	6.6	$\cdots$	300	9	0.89	20160323	$59.6\times 1$	36.8	26.1
HD 90885 A	10:28:18.50	-52:33:42.0	K0.5V(k)	10	$33.441\pm 0.051$	2	7.6	6.5	$\cdots$	300	9	0.93	20160323	$59.6\times 1$	36.8	26.2
HD 90905	10:29:42.23	+01:29:28.0	G1V	1	$31.038\pm 0.091$	2	6.3	5.6	$\cdots$	70	9	1.19	20150130	$59.6\times 1$	37.8	20.6
V419 Hya	10:43:28.32	-29:03:51.1	K1V	4	$21.540\pm 0.020$	2	6.9	5.8	AB Dor	149	5	0.90	20160226	$59.6\times 1$	26.8	147.5
HD 95086	10:57:03.09	-68:40:02.6	A8III	13	$86.44\pm 0.24$	2	7.0	6.9	LCC	15	33	1.61	20160229	$59.6\times 1$	37.8	16.9
HR 4296	11:00:08.34	-51:49:04.1	A3III/IV	27	$58.84\pm 0.13$	2	5.8	5.8	$\cdots$	456	3	1.68	20160318	$59.6\times 1$	36.8	27.3
TW Hya	11:01:51.91	-34:42:17.0	K6Ve	4	$60.09\pm 0.15$	2	9.1	7.6	TWA	10	5	0.77	20160128	$59.6\times 1$	34.8	79.4
b Leo	11:02:19.78	+20:10:47.1	kA1VmA3V	34	$38.87\pm 0.27$	11	4.2	4.3	$\cdots$	425	3	2.04	20160427	$59.6\times 1$	35.8	12.8
HR 4334	11:08:44.05	-28:04:50.2	A1V	22	$53.08\pm 0.53$	2	5.2	5.4	TWA	10	5	2.08	20160324	$59.6\times 1$	40.8	133.3
HR 4372	11:17:12.06	-38:00:51.8	A1V	35	$74.21\pm 0.29$	2	6.0	6.0	$\cdots$	431	3	1.85	20160322	$59.6\times 1$	37.8	41.3
HD 99827 A	11:25:17.70	-84:57:16.0	F5V	13	$90.28\pm 1.12$	2	7.2	6.6	$\cdots$	70	9	1.46	20160229	$59.6\times 1$	78.5	26.4
													20160229	$59.6\times 1$	78.5	26.4
HD 100546	11:33:25.51	-70:11:41.2	B9Vne	36	$110.02\pm 0.62$	2	6.2	6.0	LCC	15	33	2.21	20141217	$59.6\times 1$	32.8	12.9
89 Leo	11:34:21.95	+03:03:36.6	F5.5V	10	$27.062\pm 0.060$	2	5.3	4.8	$\cdots$	300	9	1.38	20150201	$59.6\times 1$	37.8	21.3
HR 4502	11:41:19.86	-43:05:44.4	A0V	27	$64.30\pm 0.42$	2	5.3	5.5	Argus	50	17	2.25	20160322	$59.6\times 1$	36.8	43.4
HD 102458	11:47:24.55	-49:53:03.0	G4V	4	$113.40\pm 0.43$	2	8.2	7.6	LCC	15	33	1.37	20160228	$59.6\times 1$	35.8	36.0
HD 103599	11:55:43.58	-54:10:50.4	F3IV	13	$108.03\pm 0.52$	2	8.0	7.4	LCC	15	33	1.37	20160318	$59.6\times 1$	32.8	22.6
HD 103673	11:55:57.12	-77:00:30.2	G5V	13	$46.788\pm 0.065$	2	7.9	7.0	$\cdots$	300	9	1.01	20160321	$59.6\times 1$	36.8	14.2
HD 104125	11:59:23.78	-57:10:04.7	A2V	13	$99.15\pm 1.35$	2	6.4	6.3	LCC	15	33	2.00	20160126	$59.6\times 1$	35.8	23.7
HR 4597	12:02:37.76	-69:11:32.2	B9V	13	$104.86\pm 1.01$	2	5.7	6.1	LCC	15	33	2.87	20160227	$59.6\times 1$	54.7	23.1
HD 105613	12:09:38.83	-58:20:58.7	A3V	13	$104.61\pm 0.44$	2	7.1	7.0	LCC	15	33	1.79	20160319	$59.6\times 1$	37.8	21.0
3 Crv	12:11:03.88	-23:36:08.5	A1V	16	$58.81\pm 0.66$	2	5.2	5.4	$\cdots$	388	3	1.99	20150201	$59.6\times 1$	18.9	34.9
HD 105994	12:11:58.86	-50:46:12.4	F3V	27	$112.80\pm 0.68$	2	8.2	7.6	LCC	15	33	1.48	20150129	$59.6\times 1$	35.8	35.2
HR 4669	12:17:06.36	-65:41:34.6	A0V	13	$105.94\pm 0.60$	2	5.8	6.0	LCC	15	33	2.49	20150504	$59.6\times 1$	41.8	17.8
													20160125	$59.6\times 1$	30.8	14.9
													20160430	$59.6\times 1$	48.7	19.7
HD 106906	12:17:53.23	-55:58:31.8	F5V	13	$103.33\pm 0.46$	2	7.4	6.8	LCC	15	33	1.51	20150504	$59.6\times 1$	41.8	25.3
HD 107146	12:19:06.50	+16:32:53.9	G2V	37	$27.471\pm 0.032$	2	6.4	5.6	$\cdots$	125	9	1.08	20160225	$59.6\times 1$	37.8	21.9
HR 4692	12:20:28.28	-65:50:33.5	B9V	13	$97.36\pm 0.59$	2	6.0	6.3	LCC	15	33	2.46	20160125	$59.6\times 1$	32.8	17.5
													20160430	$59.6\times 1$	14.9	6.7
HD 108857	12:30:46.32	-58:11:16.7	F7V	13	$104.52\pm 0.93$	2	7.9	7.2	LCC	15	33	1.39	20160427	$59.6\times 1$	33.8	22.6
HR 4796 A	12:36:01.08	-39:52:10.0	A0V	22	$72.78\pm 1.75$	11	5.6	5.8	TWA	10	5	2.23	20160318	$59.6\times 1$	36.8	52.7
TYC 9412-59-1	12:39:21.26	-75:02:39.2	K3Ve	4	$103.67\pm 0.29$	2	9.0	8.0	eps Cha	4	38	1.19	20160229	$59.6\times 1$	30.8	15.2
HD 110058	12:39:46.22	-49:11:55.4	A0V	27	$129.98\pm 1.33$	2	7.6	7.6	LCC	15	33	1.70	20160319	$59.6\times 1$	35.8	34.2
$\rho$ Vir	12:41:53.06	+10:14:08.3	A0Va lB	35	$36.27\pm 0.28$	11	4.7	4.8	$\cdots$	198	3	1.83	20160226	$59.6\times 1$	36.8	16.3
LO Mus	12:47:19.00	-66:14:14.8	K1Vk	10	$30.44\pm 0.25$	2	7.7	6.7	$\cdots$	70	9	0.83	20160128	$59.6\times 1$	34.8	17.1
HD 111520	12:50:19.75	-49:51:48.8	F5/6V	27	$108.94\pm 0.65$	2	8.4	7.8	LCC	15	33	1.26	20150702	$59.6\times 1$	39.8	34.7
HD 112383	12:57:26.25	-67:57:38.4	A2IV/V	13	$101.65\pm 0.45$	2	6.5	6.7	LCC	15	33	2.07	20160321	$59.6\times 1$	36.8	16.6
													20160427	$59.6\times 1$	11.9	5.4
TYC 9245-617-1	12:58:25.58	-70:28:49.2	K0Ve	4	$94.80\pm 0.26$	2	8.8	7.7	LCC	15	33	1.29	20160228	$59.6\times 1$	36.8	26.7
HD 113414	13:03:39.01	-16:20:11.7	F7/8V	16	$54.06\pm 0.15$	2	7.2	6.5	$\cdots$	70	9	1.29	20150130	$59.6\times 1$	22.9	24.5
HD 113524	13:04:59.47	-47:23:48.4	F6/8	27	$119.33\pm 0.97$	2	8.5	7.9	LCC	15	33	1.29	20150504	$59.6\times 1$	40.8	34.7
HR 4979	13:12:03.20	-37:48:11.0	G3V	22	$20.67\pm 0.12$	11	4.0	3.3	$\cdots$	125	9	1.78	20160322	$29.1\times 2$	12.6	14.1
													20160528	$18.9\times 3$	37.8	71.1
HR 5008	13:17:13.92	-43:58:46.0	A2mA4-A9/F0	27	$71.93\pm 0.96$	2	5.5	5.5	$\cdots$	907	3	1.89	20160323	$59.6\times 1$	18.9	24.3
HD 115820	13:20:26.84	-49:13:25.0	A7/8V	27	$115.90\pm 0.69$	2	7.6	7.4	LCC	15	33	1.65	20150704	$59.6\times 1$	40.8	33.2
$\iota$ Cen	13:20:35.82	-36:42:44.3	kA1.5hA3mA3Va	10	$18.021\pm 0.055$	11	2.5	2.7	Carina-Near	200	7	2.14	20150705	$23.3\times 2$	35.7	72.6
HR 5121	13:37:23.50	-46:25:40.2	B8V	27	$141.39\pm 4.53$	2	5.7	6.2	LCC	15	33	3.42	20150501	$59.6\times 1$	40.8	40.4
													20160427	$59.6\times 1$	64.6	56.6
HD 118878	13:40:37.67	-44:19:48.7	A0V	27	$122.91\pm 3.63$	2	6.3	6.3	UCL	16	33	2.21	20150405	$59.6\times 1$	42.7	41.6
HD 118991 A	13:41:44.81	-54:33:33.7	B8.5Vn	12	$86.82\pm 1.29$	2	5.1	5.4	LCC	15	33	3.21	20150404	$59.6\times 1$	36.8	25.8
HD 119727	13:46:14.54	-54:40:59.9	A1V	13	$103.76\pm 0.94$	2	6.2	6.2	LCC	15	33	2.18	20150404	$59.6\times 1$	38.8	29.8
HD 120326	13:49:54.53	-50:14:23.7	F0V	27	$113.93\pm 1.54$	2	7.9	7.6	LCC	15	33	1.49	20160229	$59.6\times 1$	36.8	30.2
HD 122973	14:05:10.69	-09:02:54.6	G2V	1	$51.85\pm 0.19$	2	7.4	6.8	$\cdots$	300	9	1.19	20160226	$59.6\times 1$	34.8	32.2
HD 123058	14:07:29.28	-61:33:43.9	F3/5V	13	$71.69\pm 0.17$	2	7.1	6.8	Argus	50	17	1.41	20150705	$59.6\times 1$	41.8	21.3
													20160321	$59.6\times 1$	21.9	11.1
HR 5357	14:19:23.91	-37:00:10.1	A2Vn	14	$68.94\pm 0.38$	2	5.7	5.7	UCL	16	33	2.12	20150406	$59.6\times 1$	8.0	47.2
HD 126488	14:27:30.48	-52:31:30.2	F2V	27	$160.44\pm 8.23$	2	8.0	7.6	UCL	16	33	1.56	20150703	$59.6\times 1$	40.8	28.7
HD 128400	14:41:52.50	-75:08:22.0	G6V	4	$20.255\pm 0.010$	2	6.0	5.2	$\cdots$	300	9	1.01	20160427	$59.6\times 1$	37.8	14.1
HD 129060	14:44:14.14	-69:40:26.8	F7V	4	$36.331\pm 0.042$	2	6.4	5.8	$\cdots$	125	9	1.24	20160301	$59.6\times 1$	58.7	32.4
													20160430	$59.6\times 1$	24.9	10.5
HD 129679	14:44:16.82	-11:37:02.5	G3V	1	$59.75\pm 0.20$	2	7.7	7.1	$\cdots$	57	9	1.17	20160225	$59.6\times 1$	36.8	34.6
HD 129926 A	14:46:00.08	-25:26:35.4	F0VSr	31	$30.28\pm 0.85$	11	4.8	4.3	$\cdots$	498	3	1.64	20150702	$59.6\times 1$	41.8	58.7
HD 130697	14:50:58.73	-42:49:20.8	A2V	27	$145.60\pm 1.70$	2	6.5	6.5	UCL	16	33	2.17	20150405	$59.6\times 1$	34.8	54.7
HD 131435	14:54:54.80	-41:21:53.2	F5V	27	$43.46\pm 0.16$	2	6.5	6.0	$\cdots$	131	9	1.29	20160526	$59.6\times 1$	36.8	49.7
HD 131835	14:56:54.48	-35:41:43.4	A2IV	22	$133.66\pm 3.55$	2	7.6	7.6	UCL	16	33	1.77	20150504	$59.6\times 1$	40.8	95.2
HD 132094	14:58:24.28	-37:21:44.7	B9V	22	$176.67\pm 4.15$	2	7.0	7.3	UCL	16	33	2.48	20150704	$59.6\times 1$	35.8	84.2
HD 133803	15:07:14.96	-29:30:15.9	A9V	22	$110.71\pm 0.79$	2	7.7	7.4	UCL	16	33	1.55	20150701	$59.6\times 1$	27.8	164.7
NY Aps	15:12:23.44	-75:15:15.6	G9V	4	$52.536\pm 0.062$	2	8.4	7.5	Argus	50	17	0.84	20160227	$59.6\times 1$	19.9	16.3
													20160430	$59.6\times 1$	36.8	13.8
HR 5663	15:15:53.66	-48:04:25.1	A2IIIs	27	$52.37\pm 0.17$	2	5.6	5.5	$\cdots$	547	3	1.66	20150729	$59.6\times 1$	41.8	29.4
$\beta$ Cir	15:17:30.85	-58:48:04.3	A3Va	10	$30.55\pm 0.18$	11	3.8	3.8	$\cdots$	536	3	1.93	20150408	$59.6\times 1$	39.8	21.2
$\gamma$ TrA	15:18:54.69	-68:40:46.1	A1V	15	$56.37\pm 0.38$	11	2.6	2.5	$\cdots$	291	3	3.41	20160317	$18.9\times 3$	35.9	17.0
HD 136482	15:22:11.27	-37:38:08.0	B8/9V	22	$148.27\pm 2.42$	2	6.5	6.8	UCL	16	33	2.71	20160228	$59.6\times 1$	37.8	71.8
HD 137919	15:30:21.31	-41:55:08.3	B9.5V	27	$154.16\pm 1.41$	2	6.2	6.5	UCL	16	33	2.87	20150406	$59.6\times 1$	41.8	35.1
HR 5751	15:32:04.19	-38:37:21.2	A3III	22	$55.32\pm 0.15$	2	5.9	5.8	UCL	16	33	1.69	20150704	$59.6\times 1$	20.9	54.8
V343 Nor	15:38:57.54	-57:42:27.3	K0V	4	$40.11\pm 0.10$	2	7.2	6.0	$\beta$ Pic	26	6	1.18	20150703	$59.6\times 1$	34.8	23.9
g Lup	15:41:11.38	-44:39:40.3	F3/5V	27	$17.437\pm 0.049$	11	4.1	3.7	$\cdots$	200	9	1.42	20160322	$52.4\times 1$	39.3	43.6
HD 139883	15:41:53.23	-34:53:19.6	F2V	22	$111.03\pm 0.64$	2	7.9	7.5	UCL	16	33	1.51	20150504	$59.6\times 1$	38.8	97.9
HD 141254	15:49:39.65	-38:46:38.7	F3V	22	$152.83\pm 1.01$	2	8.7	8.1	UCL	16	33	1.35	20160427	$59.6\times 1$	32.8	59.0
HD 141569	15:49:57.76	-03:55:16.2	A2V(e)kB9mB9 lB	39	$110.63\pm 0.55$	2	6.8	6.9	$\cdots$	5	40	2.04	20150404	$59.6\times 1$	40.8	24.4
$ϵ$ Ser	15:50:48.97	+04:28:39.8	kA2hA5mA7V	41	$21.598\pm 0.089$	11	3.5	3.4	$\cdots$	553	3	1.76	20150801	$59.6\times 1$	30.8	18.6
NZ Lup	15:53:27.29	-42:16:00.7	G2	4	$60.34\pm 0.18$	2	7.2	6.4	UCL	16	33	1.29	20150408	$59.6\times 1$	41.8	47.9
HD 142315	15:54:41.60	-22:45:58.3	B9V	16	$145.34\pm 1.11$	2	6.6	6.7	US	10	33	4.17†	20150501	$59.6\times 1$	25.8	47.6
HD 142705	15:56:47.86	-23:11:02.5	A0V	16	$144.32\pm 1.59$	2	7.4	7.3	US	10	33	2.05	20150730	$59.6\times 1$	29.8	55.5
HD 142851	15:57:59.36	-31:43:43.9	A0V	22	$144.05\pm 1.24$	2	6.8	7.1	UCL	16	33	2.43	20160526	$59.6\times 1$	37.8	156.6
48 Lib	15:58:11.38	-14:16:45.5	B8Ia/Iab	16	$143.46\pm 4.96$	11	4.7	4.8	US	10	33	4.68	20160326	$59.6\times 1$	36.8	57.4
HD 142678	16:02:16.06	-73:30:54.8	A6III	13	$69.25\pm 0.22$	2	6.3	6.1	$\cdots$	666	3	1.61	20160918	$59.6\times 1$	37.8	15.6
HD 143811	16:03:33.44	-30:08:13.2	F5V	22	$137.14\pm 0.82$	2	8.4	7.8	US	10	33	1.35	20160430	$59.6\times 1$	35.8	176.0
HR 6051	16:14:28.89	-21:06:27.2	B9V	16	$112.18\pm 2.47$	2	6.2	6.4	US	10	33	2.56	20150703	$59.6\times 1$	39.8	58.3
HD 145689 A	16:17:05.52	-67:56:28.0	A3V	13	$55.55\pm 0.18$	2	5.7	5.7	Argus	50	17	1.80	20150705	$59.6\times 1$	40.8	18.5
d Sco	16:18:17.90	-28:36:50.5	A1Va	12	$41.29\pm 0.38$	11	4.6	4.9	$\beta$ Pic	26	6	2.04	20150702	$59.6\times 1$	41.8	148.9
V371 Nor	16:19:15.90	-55:30:17.0	K3Ve	4	$28.113\pm 0.029$	2	8.2	6.7	$\cdots$	125	9	0.78	20150504	$59.6\times 1$	41.8	23.9
													20150730	$59.6\times 1$	35.8	23.0
													20160430	$59.6\times 1$	10.9	6.4
HR 6094	16:24:01.30	-39:11:35.0	G5V	4	$12.908\pm 0.023$	2	4.8	4.0	$\cdots$	300	9	1.10	20150501	$59.6\times 1$	39.8	63.3
$\zeta$ TrA	16:28:28.14	-70:05:03.8	F9V	10	$12.177\pm 0.025$	2	4.4	3.6	$\cdots$	125	9	1.54†	20150406	$59.6\times 1$	41.8	17.6
HD 149813	16:39:47.00	-51:59:39.0	G3V	4	$31.428\pm 0.044$	2	6.9	6.1	$\cdots$	300	9	1.01	20160328	$59.6\times 1$	35.8	30.2
HD 152555	16:54:08.14	-04:20:24.7	F8/G0V	1	$45.316\pm 0.088$	2	7.1	6.5	AB Dor	149	5	1.16	20150703	$59.6\times 1$	40.8	22.5
HR 6297	17:00:06.24	-54:35:49.2	A5IV/V	13	$53.25\pm 0.36$	2	5.3	5.3	$\cdots$	765	3	1.75	20150704	$59.6\times 1$	41.8	25.2
HD 155114	17:11:51.24	-45:52:27.6	G3V	27	$32.480\pm 0.049$	2	6.9	6.1	$\cdots$	300	9	1.05	20150408	$59.6\times 1$	40.8	40.5
													20150730	$59.6\times 1$	9.9	13.0
V824 Ara	17:17:25.50	-66:57:03.7	K1Vp	13	$30.506\pm 0.025$	2	6.0	4.9	$\beta$ Pic	26	6	2.26†	20150404	$59.6\times 1$	38.8	17.9
HD 156751	17:22:47.89	-58:28:23.7	A7II/III	13	$68.25\pm 0.21$	2	6.5	6.3	$\cdots$	396	3	1.56	20160321	$59.6\times 1$	36.8	20.8
TYC 8728-2262-1	17:29:55.07	-54:15:48.6	K1V	4	$67.76\pm 0.16$	2	8.5	7.5	$\beta$ Pic	26	6	1.04	20160430	$59.6\times 1$	37.8	25.7
$\pi$ Ara	17:38:05.52	-54:30:01.6	A5IV/V	13	$44.56\pm 0.52$	11	4.9	4.9	$\cdots$	872	3	1.72	20160427	$59.6\times 1$	37.8	25.1
$\gamma$ Oph	17:47:53.58	+02:42:26.9	A1VnkA0mA0	41	$31.52\pm 0.21$	11	3.5	3.7	$\cdots$	520	3	2.14	20150702	$59.6\times 1$	41.8	16.8
HD 164249 A	18:03:03.36	-51:38:55.7	F6V	4	$49.61\pm 0.12$	2	6.5	6.0	$\beta$ Pic	26	6	1.38	20150501	$59.6\times 1$	41.8	29.6
													20150829	$59.6\times 1$	23.9	17.5
HR 6748	18:06:23.70	-36:01:11.0	G1V	4	$17.204\pm 0.034$	2	5.3	4.6	$\cdots$	300	9	1.09	20150405	$59.6\times 1$	32.8	98.7
HD 167425 A	18:19:40.13	-63:53:11.6	F9.5V	10	$23.486\pm 0.022$	2	5.6	4.9	$\cdots$	300	9	1.21	20150405	$59.6\times 1$	38.8	19.7
$\varphi$ Oct	18:23:36.44	-75:02:39.6	A0V	13	$59.71\pm 0.46$	2	5.2	5.4	$\cdots$	254	3	2.12	20150831	$59.6\times 1$	37.8	14.5
TYC 8757-314-1	18:30:11.98	-58:16:27.7	Mwl	42	$14.379\pm 0.014$	2	7.9	6.2	AB Dor	149	5	0.63	20160430	$59.6\times 1$	37.8	19.3
HR 6948	18:33:00.92	-39:53:31.3	F5V	10	$37.055\pm 0.058$	2	5.9	5.3	$\cdots$	125	9	1.46	20150701	$59.6\times 1$	41.8	50.4
V889 Her	18:34:20.16	+18:41:24.7	G2V	15	$35.375\pm 0.059$	2	6.7	5.9	$\cdots$	10	9	1.02	20150730	$59.6\times 1$	39.8	16.1
HR 7012	18:45:26.90	-64:52:16.5	A7V	10	$28.55\pm 0.15$	11	4.4	4.3	$\beta$ Pic	26	6	1.70	20150408	$59.6\times 1$	41.8	20.8
111 Her	18:47:01.27	+18:10:53.5	A5III	15	$28.89\pm 0.38$	11	4.2	4.4	$\cdots$	420	3	1.78	20150829	$59.6\times 1$	37.8	13.6
PZ Tel	18:53:05.90	-50:10:50.0	G9IV	4	$47.13\pm 0.13$	2	7.6	6.5	$\beta$ Pic	26	6	1.10	20150730	$59.6\times 1$	38.8	31.8
													20160920	$59.6\times 1$	13.9	4.7
V1285 Aql	18:55:27.36	+08:24:09.6	M3Ve	43	$11.1359\pm 0.0061$	2	7.8	5.7	$\cdots$	70	9	0.32	20160919	$59.6\times 1$	36.8	16.9
HD 175726	18:56:37.20	+04:15:55.2	G0V	20	$26.621\pm 0.038$	2	6.1	5.4	$\cdots$	300	9	1.11	20160921	$59.6\times 1$	37.8	21.0
$\rho$ Tel	19:06:19.92	-52:20:26.3	F6V	27	$58.86\pm 0.50$	2	4.7	4.0	TucHor	45	5	2.64	20150504	$59.6\times 1$	41.8	29.4
$\alpha$ CrA	19:09:28.34	-37:54:16.1	A2Va	10	$38.43\pm 0.37$	11	3.9	3.9	$\cdots$	620	3	2.10	20150801	$59.6\times 1$	41.8	82.7
$\eta$ Tel A	19:22:51.21	-54:25:26.2	A0V	44	$48.22\pm 0.49$	11	4.8	5.1	$\beta$ Pic	26	6	2.15	20150501	$59.6\times 1$	40.8	26.7
HD 181327	19:22:58.90	-54:32:17.0	F6V	4	$48.21\pm 0.13$	2	6.6	6.0	$\beta$ Pic	26	6	1.39	20150703	$59.6\times 1$	41.8	26.8
Rukbat	19:23:53.18	-40:36:57.4	B8V	27	$55.74\pm 0.68$	11	3.8	4.2	AB Dor	149	5	2.92	20150704	$59.6\times 1$	41.8	55.8
c Aql	19:29:00.99	+01:57:01.9	A0IVp	14	$60.68\pm 0.27$	2	5.5	5.6	$\cdots$	448	3	1.84	20150831	$59.6\times 1$	36.8	17.7
HD 183414	19:35:09.70	-69:58:32.0	G3V	4	$38.184\pm 0.055$	2	7.2	6.4	$\cdots$	125	9	1.06	20150504	$59.6\times 1$	41.8	16.6
V4200 Sgr	19:54:17.75	-23:56:27.9	K3VaCN1	8	$14.1069\pm 0.0079$	2	5.4	4.2	$\cdots$	125	9	1.75†	20150702	$59.6\times 1$	39.8	73.7
$ϵ$ Pav	20:00:35.55	-72:54:37.8	A0Va	10	$32.22\pm 0.18$	11	3.7	3.8	Argus	50	17	2.30	20150729	$59.6\times 1$	33.8	18.5
HD 191089	20:09:05.21	-26:13:26.5	F5V	22	$50.14\pm 0.11$	2	6.8	6.1	$\beta$ Pic	26	6	1.35	20150831	$59.6\times 1$	36.8	6.9
Peacock	20:25:38.85	-56:44:05.6	B2IV	15	$54.82\pm 1.57$	11	1.9	2.5	TucHor	45	5	9.69†	20150704	$11.6\times 4$	34.1	24.4
HD 194188	20:26:04.75	-46:39:35.8	A9III/IV	27	$66.93\pm 0.26$	2	6.4	6.2	$\cdots$	699	3	1.57	20150501	$59.6\times 1$	28.8	27.4
HD 195005	20:28:49.92	-01:44:04.1	F7V	1	$32.129\pm 0.053$	2	6.3	5.6	$\cdots$	300	9	1.20	20150830	$59.6\times 1$	29.8	23.1
o Cap B	20:29:52.56	-18:35:10.3	A7/8V	16	$71.11\pm 0.27$	2	6.4	6.3	$\cdots$	387	3	1.65	20160921	$59.6\times 1$	35.8	45.5
$\iota$ Del	20:37:49.10	+11:22:39.7	A1IV	14	$59.47\pm 0.49$	2	5.2	5.4	$\cdots$	66	3	3.34†	20160919	$59.6\times 1$	37.8	17.5
AU Mic	20:45:09.53	-31:20:27.2	M1VeBa1	8	$9.7248\pm 0.0046$	2	6.9	4.8	$\beta$ Pic	26	6	0.64	20141110	$59.6\times 1$	43.7	4.7
													20150703	$59.6\times 1$	36.8	166.6
HR 8013	20:56:47.27	-26:17:46.4	F6V	10	$21.287\pm 0.068$	2	5.2	4.7	$\cdots$	200	9	1.24	20150901	$59.6\times 1$	36.8	97.8
HD 200968 A	21:07:10.40	-13:55:23.0	K1IV	4	$17.411\pm 0.016$	2	6.3	5.1	$\cdots$	300	9	0.91	20160919	$59.6\times 1$	37.8	33.1
HD 200798	21:07:51.19	-54:12:58.8	A5/6IV/V	13	$69.68\pm 0.30$	2	6.3	6.1	TucHor	45	5	1.76	20150702	$59.6\times 1$	39.8	26.3
HD 202917	21:20:49.92	-53:02:02.4	G7V	4	$46.849\pm 0.085$	2	7.9	7.0	TucHor	45	5	1.02	20150704	$59.6\times 1$	39.8	28.1
V390 Pav	21:21:24.49	-66:54:57.4	K2V	4	$31.980\pm 0.044$	2	7.9	6.5	$\beta$ Pic	26	6	0.81	20150504	$59.6\times 1$	41.8	17.8
HD 203244	21:24:40.64	-68:13:40.2	G5V	4	$20.81\pm 0.15$	2	6.2	5.3	$\cdots$	300	9	0.99	20150831	$59.6\times 1$	36.8	16.6
IK Peg	21:26:26.61	+19:22:32.2	kA6hA9mF0+DA	45	$47.33\pm 0.32$	2	5.7	5.5	$\cdots$	439	3	1.53	20150829	$59.6\times 1$	37.8	12.3
HD 205905	21:39:10.20	-27:18:24.0	G1.5IV-V	8	$26.302\pm 0.041$	2	6.2	5.4	$\cdots$	300	9	1.14	20150901	$59.6\times 1$	34.8	61.0
42 Cap	21:41:32.86	-14:02:51.4	G2IV	8	$33.23\pm 0.35$	11	4.5	3.6	$\cdots$	300	9	3.34†	20150830	$59.6\times 1$	34.8	38.8
													20160918	$59.6\times 1$	24.9	29.0
HN Peg	21:44:31.33	+14:46:19.0	G0V	15	$18.128\pm 0.020$	2	5.3	4.6	$\cdots$	200	9	1.10	20150829	$59.6\times 1$	37.8	14.5
HD 207575	21:52:09.72	-62:03:08.5	F6V	13	$47.039\pm 0.079$	2	6.8	6.1	TucHor	45	5	1.32	20150801	$59.6\times 1$	35.8	21.0
HIP 108405	21:57:41.20	-51:00:22.0	M2Ve	4	$14.8604\pm 0.0087$	2	8.3	6.1	$\cdots$	30	9	0.40	20160918	$59.6\times 1$	37.8	27.4
$ϵ$ Ind	22:03:21.70	-56:47:10.0	K5V	4	$3.6224\pm 0.0037$	11	3.5	2.3	$\cdots$	500	9	0.81	20150703	$11.6\times 5$	38.8	26.5
Alnair	22:08:13.99	-46:57:39.5	B6V	10	$30.97\pm 0.20$	11	1.7	2.0	AB Dor	149	5	3.51	20150702	$11.6\times 4$	32.6	35.7
HR 8547	22:28:37.67	-67:29:20.6	A3/5V	13	$43.93\pm 0.18$	2	5.2	5.1	$\cdots$	508	3	1.73	20150801	$59.6\times 1$	41.8	18.9
TYC 9340-437-1	22:42:48.92	-71:42:21.3	K7Ve	4	$36.661\pm 0.032$	2	8.7	7.1	$\beta$ Pic	26	6	0.80	20160526	$59.6\times 1$	35.8	30.8
TYC 1701-642-1	22:44:41.54	+17:54:18.3	K0	15	$49.74\pm 0.28$	2	8.5	7.2	AB Dor	149	5	0.97	20160920	$59.6\times 1$	35.8	13.5
TW PsA	22:56:24.10	-31:33:56.0	K4Ve	8	$7.6082\pm 0.0050$	2	5.5	3.8	$\cdots$	200	9	0.81	20150704	$59.6\times 1$	41.8	165.9
Fomalhaut	22:57:39.05	-29:37:20.0	A4V	10	$7.704\pm 0.028$	11	0.8	0.9	$\cdots$	749	3	1.80	20150829	$2.9\times 10$	13.6	1.7
													20150830	$2.9\times 10$	19.9	166.1
HD 217343	23:00:19.29	-26:09:13.5	G5V	4	$31.780\pm 0.041$	2	6.8	6.0	AB Dor	149	5	1.06	20141111	$59.6\times 1$	40.8	38.7
HR 8799	23:07:28.71	+21:08:03.3	F0+VkA5mA5	41	$41.29\pm 0.15$	2	5.6	5.3	Columba	42	5	1.55	20160919	$59.6\times 1$	59.6	20.9
HR 8843	23:16:57.69	-62:00:04.3	F6V	10	$20.460\pm 0.026$	2	5.1	4.6	$\cdots$	200	9	1.25	20150702	$59.6\times 1$	41.8	21.4
$\kappa$ Psc	23:26:55.96	+01:15:20.2	A2VpSrCrSi	14	$47.06\pm 0.64$	11	4.7	5.0	AB Dor	149	5	2.07	20150901	$59.6\times 1$	35.8	20.6
HIP 116003	23:30:13.44	-20:23:27.5	M2Ve	4	$15.956\pm 0.023$	2	8.6	6.6	Columba	42	5	0.46	20160920	$59.6\times 1$	37.8	51.4
HD 221231	23:31:02.85	-69:04:36.3	G1V	4	$30.313\pm 0.025$	2	6.5	5.7	$\cdots$	300	9	1.10	20150801	$59.6\times 1$	32.8	15.5
HD 221853	23:35:36.15	+08:22:57.4	F0	46	$65.41\pm 0.21$	2	7.0	6.4	$\cdots$	125	9	1.46	20151105	$59.6\times 1$	8.0	3.9
													20151106	$59.6\times 1$	21.9	8.7
HD 224228	23:56:10.67	-39:03:08.4	K2V	4	$21.927\pm 0.046$	2	7.3	6.0	AB Dor	149	5	0.84	20141110	$59.6\times 1$	39.8	64.7
$\eta$ Tuc	23:57:34.97	-64:17:53.1	A1V	13	$47.44\pm 1.10$	11	4.8	4.9	TucHor	45	5	2.14	20150703	$59.6\times 1$	40.8	19.1

Equations

$$P(\delta |\text{data})\propto {\int }_0^{\infty }{P}_{HM}(f|\text{data})P_{LM}(f-\delta |\text{data})df$$

$$\frac{d^{2}N}{dmda}=f{C}_1{m}^{\alpha }{a}^{\beta }$$

$$C_1={\left[{\int }_{m1}^{m2}{\int }_{a1}^{a2}{m}^{\alpha }{a}^{\beta }dmda\right]}^{-1}$$

$$C_1=\frac{\alpha +1}{m_2^{(\alpha +1)}-m_1^{(\alpha +1)}}\frac{\beta +1}{a_2^{(\beta +1)}-a_1^{(\beta +1)}}$$

$$P(f,\alpha ,\beta |\mathrm{data})\propto P(\mathrm{data}|f,\alpha ,\beta )P(f)P(\alpha )P(\beta )$$

$$P(\mathrm{data}|f,\alpha ,\beta )=\prod _i\prod _j\frac{e^{-E_{i,j}}{E}_{i,j}^{O_{i,j}}}{O_{i,j}!}$$

$$\frac{d^{2}N}{dmda}=f{C}_1{m}^{\alpha }{a}^{\beta }{\left(\frac{M_{\ast }}{1.75M_{\odot }}\right)}^{\gamma }$$

$$P(f,\alpha ,\beta ,\gamma |\mathrm{data})\propto P(\mathrm{data}|f,\alpha ,\beta ,\gamma )P(f)P(\alpha )P(\beta )P(\gamma )$$

Full text

††software: Astropy (Astropy Collaboration et al., 2013), emcee (Foreman-Mackey et al., 2013), GPI DRP (Perrin et al., 2016), IDL Astronomy Library (Landsman, 1993), pyKLIP (Wang et al., 2015)

1 Introduction

Since giant extrasolar planets ($\gtrsim$1 MJup) are typically the easiest to detect for most search techniques, there is a long record of studies that investigate the characteristics of the underlying population. A principal scientific goal is to distinguish between different planet formation and evolution theories by measuring planet occurrence rates with respect to the properties of the planets (e.g., mass, semi-major axis, eccentricity) and host stars (e.g., mass and metallicity). A correlation between giant planet frequency and stellar metallicity for FGK stars was among the first empirical findings using radial velocity (RV) data (Gonzalez, 1997; Fischer & Valenti, 2005). Tabachnik & Tremaine (2002) and Cumming et al. (2008) analyzed datasets of RV-detected planets and completeness to derive occurrence rates for giant planets, fitting power law distributions to the population in mass and semi-major axis. Johnson et al. (2007) and Johnson et al. (2010a) observed a correlation between RV giant planet occurrence rate and stellar host mass, with more massive stars more likely to host giant planets with $a<2.5$ au, though there is uncertainty in deriving a main sequence mass for subgiants (Lloyd, 2011; Johnson & Wright, 2013; Lloyd, 2013). Results from asteroseismology suggest spectroscopically-derived stellar masses are too high (Johnson et al., 2014; Campante et al., 2017), but not significantly so (North et al., 2017; Stello et al., 2017). Thus the underlying correlation appears robust: for RV-detected giant planets the occurrence rate is proportional to $M_{*}^{\gamma}$, with $\gamma=1.05^{+0.28}_{-0.24}$ (Ghezzi et al., 2018).

Direct imaging in the infrared targets young, nearby stars to detect self-luminous giant planets at high contrast at separations of $\lesssim 1^{\prime\prime}$ from their parent stars. This technique is sensitive to giant planets at wider separations ($\gtrsim 5$ au), and has thus far been detecting planets primarily around higher-mass stars ($\sim$1.5-1.8 M⊙, e.g., Marois et al. 2010, Lagrange et al. 2010), despite the lower intrinsic luminosity of lower-mass stars making it easier to image planetary-mass companions. While higher-mass stars hosting most directly imaged planets is suggestive that the RV trend of rising giant planet fraction with increasing stellar mass holds at wider separations, it has yet to be robustly demonstrated (Wahhaj et al., 2013; Galicher et al., 2016; Bowler, 2016). Lannier et al. (2016) found that while imaged intermediate mass-ratio companions (e.g., brown dwarfs) appeared equally common around higher-mass stars (spectral types A and F) and lower-mass stars (M dwarfs), for low mass ratio companions (such as giant planets) there was a 74.5% probability that the distributions around the higher-mass and lower-mass stars were different. With new observations and a larger statistical sample we reexamine this question.

Brown dwarfs, intermediate in mass between stars and giant planets, have typically been detected in greater numbers by imaging surveys compared to giant planets, despite relatively low brown dwarf occurrence rates (Metchev & Hillenbrand, 2009). Recent advances in high contrast imaging instrumentation (Macintosh et al., 2014; Beuzit et al., 2008), observing techniques (Marois et al., 2006), and data reduction (Lafrenière et al., 2007; Soummer et al., 2012; Pueyo, 2016; Wang et al., 2015) have led to detections of planets orbiting HR 8799 (Marois et al., 2008, 2010), $\beta$ Pictoris (Lagrange et al., 2009), HD 95086 (Rameau et al., 2013b, a), 51 Eridani (Macintosh et al., 2015), and HIP 65426 (Chauvin et al., 2018), among others. Brandt et al. (2014) proposed that such wide separation (10–100 au) giant planets and brown dwarfs are part of the same underlying population. This work tests this hypothesis with new observations.

2 Survey Description

The Gemini Planet Imager (GPI) is an instrument optimized to directly detect planets at high contrast within $\sim 1^{\prime\prime}$ radius from their parent star. Coupling a woofer/tweeter extreme adaptive optics system (Poyneer et al., 2014; Bailey et al., 2016) with an apodized Lyot coronagraph (Soummer et al., 2007; Savransky et al., 2014) and an integral field spectrograph (Chilcote et al., 2012; Larkin et al., 2014), GPI routinely reaches contrasts between $\sim 10^{-5}$ and $10^{-6}$ at $0.4^{\prime\prime}$ from bright stars (m${}_{I}\lesssim$ 8.0 mag) after post-processing (Macintosh et al., 2014; Ruffio et al., 2017b). The coronagraphic inner working angle (IWA) is 2.8$\lambda/D$, or 0.12*′′* in $H$-band, though in practice the contrast near the IWA after postprocessing depends on additional factors such as the total sky rotation achieved during the observations of each star.

The Gemini Planet Imager Exoplanet Survey (GPIES) utilizes GPI for a 600-star survey, begun in late 2014 and set to conclude in 2019 (Macintosh et al., 2018). Based at Gemini South at Cerro Pachon, Chile, GPIES is targeting the youngest, closest stars in the sky in a systematic survey for giant planets between 10–100 au. In particular, the goal of the planet-finding campaign is to directly measure the planet occurrence rate, determine how planet occurrence correlates with stellar properties, probe the atmospheres of giant planets and brown dwarfs, and determine the orbital properties of substellar companions. The debris disk portion of GPIES is described in depth in Esposito et al. (2019, in prep.).

Here we present the results from the first 300 stars observed by GPIES. As the first half of the planned survey, 300 stars represents the natural midcourse point to consider the occurrence rate of giant planets partway through the survey. The 300th star was observed on UT 2016-09-21, and in this work we consider only GPIES campaign observations on or before this date. Candidate companions from direct imaging observations can be reduction artifacts, background stars, or bound companions. Second epoch observations confirm that a source is astrophysical and not spurious, and common proper motion indicates whether the star is a bound companion or a background star (e.g., De Rosa et al. 2015, Nielsen et al. 2017). For these 300 stars we have completed second epoch observations of all but seven of the 55 candidate companions, and most of these remaining seven have spectra consistent with background stars; follow-up observations are ongoing for the stars observed after this date. Thus the 300-star sample represents an essentially complete sample where candidate companions have been classified into PSF reduction artifacts, background stars, or bound companions.

2.1 Target Selection and Properties

A sample of young, nearby stars was constructed by combining members of kinematic associations (Zuckerman et al., 2001, 2011; de Zeeuw et al., 1999) with X-ray selected FGKM stars within 100 pc that also satisfy the observing limits (m${}_{I}\lesssim 9.0$ mag, $\delta\leq 25\arcdeg$). Candidate B and A stars were selected for youth using evolutionary tracks (Siess et al., 2000). We obtained echelle spectra for approximately 2000 candidate stars to further restrict the sample to the $\sim$800 most promising candidates by measuring lithium abundance and chromospheric activity indicators. Apparent binaries with angular separation $0.02$–$3.0\arcsec$ and $\Delta$mag $\leq 5$ (including those discovered during the course of the campaign) were then removed because the presence of the companion degrades closed-loop performance. The target list was updated in late-2016 to replace newly-resolved binaries with stars from the original sample, and newly-identified moving group members. The final GPIES sample consisted of 602 nearby stars. Almost half of the final sample are members of kinematic associations; 152 are from nearby moving groups and 104 from the more distant Scorpius-Centaurus OB2 association. Of the remaining 346 stars not in known moving groups or associations, 92 were B and A-type stars selected based on their position on the color-magnitude diagram, and 254 were solar-type stars selected based on the strength of chromospheric activity indicators (e.g., Ca H&K). A few stars in the sample were selected for other indicators of planetary systems such as debris disk structure. For example, even though the age of Fomalhaut is $>$400 Myr (Mamajek et al., 2013), the stellocentric offset of the outer belt (Kalas et al., 2005) and eccentric orbit of Fomalhaut b could be the dynamical outcome of a yet-to-be-detected Fomalhaut c in the system (Kalas et al., 2013).

The distributions of various target properties for the first 300 stars observed during the GPIES campaign are shown in Figure 1. Parallax measurements were obtained from the Gaia DR2 catalog for 251 stars (Gaia Collaboration et al., 2018). We used Hipparcos parallaxes from van Leeuwen (2007) for the remaining 49 stars because these were either too bright for Gaia or had worse astrometric precision in the DR2. Spectral types were obtained from SIMBAD, compiled from a number of literature sources. The properties for each target are given in the observing log in Table A.

The ages of the stars in the 300-star sample were estimated using a number of methods: 145 based on membership of a nearby kinematic moving group or association (Zuckerman et al., 2001, 2011), 109 based on activity indicators in their optical spectra or measured X-ray luminosity, and 48 based on their position on the color-magnitude diagram. For these 109 stars, the process to derive their ages from X-ray and spectroscopic indicators, as well as the new high-resolution spectra obtained for a subset, will be discussed in more detail in an upcoming paper. The ages of the kinematic groups are derived from isochrone fitting (e.g., Naylor & Jeffries, 2006), lithium depletion boundary analysis (e.g., Binks & Jeffries, 2014), trace-back analysis (e.g., Riedel et al., 2017) and eclipsing or resolved spectroscopic binaries (e.g., Nielsen et al., 2016). Applying these various techniques to many coeval stars simultaneously provides tight constraints on the ages of these groups, and as such these ages are preferred.

For the stars not in a known kinematic group, we instead estimated their ages by the measured decrease in lithium abundance and chromospheric activity indicators for lower-mass stars, and the evolution across the color-magnitude diagram for higher-mass stars. We used a combination of lithium abundance, Ca II H/K and H$\alpha$ emission line strength, and X-ray luminosity to assign an age to the lower-mass stars (F-type and later) based on a comparison to similar measures of stars within open clusters and kinematic associations (Mamajek & Hillenbrand, 2008; Soderblom, 2010). These methods cannot be applied to the more massive stars within the sample (B through early F-type). For example, lithium depletion only occurs efficiently in stars that have convective envelopes deep enough to transport lithium from the photosphere to regions where the temperature is sufficient for lithium burning to occur. As the depth of the convective envelope is inversely proportional to the mass of the star, this process only occurs efficiently in lower-mass stars (Castro et al., 2016). Early-type stars also do not possess the magnetic field necessary to generate the chromospheric activity indicators observed in later-type stars. Instead, we rely on the rapid evolution of these more massive stars across the color-magnitude diagram to determine their ages. The position of early-type (B and A-type) stars on the color-magnitude diagram is primarily dictated by their age and mass. These stars rapidly evolve onto the zero-age main sequence (ZAMS) where they spend roughly one-third of their main sequence lifetime at their ZAMS location on the HR diagram, before the rising helium abundance in their cores causes them to become cooler and more luminous. Secondary effects include metallicity, rotation, binarity, and extinction, which may need to be considered when attempting to determine the age of a star from its position on the color-magnitude diagram. Table A presents our best estimate for the ages and masses of each target star, the median of the distribution if posteriors are available. We are currently exploring the production of posteriors for the 109 stars that rely on X-ray, calcium, and lithium age indicators, and will present the results of this analysis and age and mass posteriors for the full sample in an upcoming paper.

2.1.1 Higher-mass stars: ages and masses

We combined the stellar evolutionary model MESA (Paxton et al., 2010), used to predict how the fundamental properties of a star evolves with time, with the ATLAS9 model atmospheres (Castelli & Kurucz, 2004), used to predict the emergent flux from the stellar photosphere, to estimate both the ages and masses of higher-mass stars and the masses of lower-mass stars based on their position on the color-magnitude diagram. We first generated a grid of evolutionary models spanning a range of masses ($0.6\leq M/{\rm M}_{\odot}\leq 9.8$), metallicities ($-0.5\leq[M/H]\leq 2.0$), and initial rotation rates ($0.0\leq\Omega/\Omega_{c}\leq 0.7$ where $\Omega/\Omega_{c}$ denotes the rotation rate as a fraction of the critical rotation rate, see Paxton et al., 2013), assuming a solar abundance of $X_{\odot}=0.7154$, $Y_{\odot}=0.2703$, and $Z_{\odot}=0.0142$ (Asplund et al., 2009). These models were generated using the same parameters used to create the MESA Isochrones and Stellar Tracks (MIST) model grid (Dotter, 2016; Choi et al., 2016), and our models are consistent with their grid at the rotation rates they computed. Rotation is treated by initializing solid-body rotation at the ZAMS, and we exclude the pre-main sequence (PMS) portion of the evolutionary track given the discontinuity seen at the ZAMS for the rapidly-rotating stars. For higher-mass stars, the PMS is very brief ($\sim 10$ Myr), and it is far more likely that a star consistent with the location of a PMS A-type star on the color-magnitude diagram is actually an older star evolving away from the ZAMS. As with the MIST model grid, stars with $M<1.2$ M*⊙* were only computed without rotation.

We collected optical photometry ($G_{B}$, $G$, $G_{R}$) from the Gaia catalogue (Gaia Collaboration et al., 2018), combining the catalog uncertainties with the systematic uncertainties reported in Evans et al. (2018). For the handful of stars too bright for Gaia, we used Tycho2 photometry ($B_{T}$, $V_{T}$; Høg et al., 2000) instead. We fit these measurements to synthetic photometry derived from a combination of our evolutionary model and the ATLAS9 model atmospheres. We used the emcee affine-invariant Markov Chain Monte Carlo (MCMC) sampler (Foreman-Mackey et al., 2013) to sample the posterior distributions of star mass $M$, age $t$, metallicity $[M/H]$, initial rotation $\Omega/\Omega_{c}$, inclination of the rotation axis $i$, and parallax $\pi$. We used the following priors:

• •

uniform for $t$ and $\cos i$,

• •

a Kroupa (2001) initial mass function for $M$,

• •

a normal distribution for $[M/H]$ ($-0.05\pm 0.11$; Nielsen et al., 2013),

• •

a mass-dependent Maxwellian distribution for the rotation rate $v$ (Zorec & Royer, 2012),

• •

a combination of a normal distribution (Gaia or Hipparcos measurement and uncertainty) and an assumed uniform space density with an exponential drop-off in distance (Bailer-Jones, 2015) for $\pi$.

If the star is part of a known kinematic group we used a normal distribution described by the age and corresponding uncertainty of that group as the age prior. We also use $v\sin i$ measurements as an additional prior, when available. At each step in the MCMC chain, the rotation rate and radius of the star are combined to give equatorial velocity. This is combined with the inclination angle, also an MCMC parameter, to give $v\sin i$ for the star at that step in the chain. The prior on this derived quantity has a form of a gaussian centered on the measured value with a sigma of the reported error. For known spectroscopic binaries, each stellar component of the system was included in the fit. Three parameters ($M$, $\Omega/\Omega_{c}$, $i$) were added for each additional star in the system, with a single age, metallicity, and parallax describing the system. The mass ratio measured from the orbit fit to the radial velocities of both components in the double-lined systems was used as an additional prior. For systems that either eclipsed or contained three stars, the mass measurements for each star were used as priors instead. The masses reported in Table A for these spectroscopic systems are the total system mass, rather than the mass of the brightest component.

At each step in the MCMC chain the evolutionary model grid was interpolated to determine the luminosity, equatorial radius, and current rotation rate of the star. A two-dimensional model of the star was constructed using the formalism described in Lovekin et al. (2006) and Espinosa Lara & Rieutord (2011) to account for the latitudinal dependence of the effective temperature and surface gravity, and for the angle between the rotation axis and the observer. We sampled the surface of the star equally in latitude and longitude ($90\times 180$) to model the star. The emergent flux at each point within each bandpass was computed from the ATLAS9 model atmosphere given the the local effective temperature and surface gravity, the metallicity of the star, and the viewing angle between the surface normal and the observer. The total flux of the star in each bandpass was calculated by weighting the flux from each point by the projected surface area of each segment as viewed from Earth; segments not visible from Earth have zero projected surface area. The total flux was scaled by the distance to the star and compared to the measured photometry. The MCMC chains were advanced until the median and $1\sigma$ credible region derived from the last half of the chain were no longer evolving.

The ages we derive for these higher-mass stars are consistent with previous works using this method. For example, we derive an age of $750^{+170}_{-190}$ for Fomalhaut, consistent with $520\pm 100$ Myr from Nielsen et al. (2013) and $520^{+100}_{-60}$ from David & Hillenbrand (2015). Mamajek (2012) finds an age for the system of $440\pm 40$ Myr, based largely on an analysis of the age of the companion K4V star TW PsA (Fomalhaut B). This is similar to the two previous analyses of Fomalhaut itself, and about $\sim$1.5 $\sigma$ lower than our estimate. This age discrepancy is largely a result of different metallicity values used in these analyses, with the Siess et al. (2000) models used by Nielsen et al. (2013) assuming a solar metallicity of $Z_{\odot}=0.02$, David & Hillenbrand (2015) using the PARSEC model grid which assumes $Z_{\odot}=0.0152$, while Mamajek (2012) used a $Z_{\odot}=0.017$ and the Bertelli et al. (2008) model grid. We use a lower value of $Z_{\odot}=0.0142$ along with the model grid discussed above. By increasing the metallicity from 0.0142 to 0.017 to match Mamajek (2012), we recover a more consistent age of $550\pm 70$ Myr. Similarly, increasing the metallicity further to $Z_{\odot}=0.02$ gives a still lower age, of $440\pm 70$ Myr. As we note above, the Asplund et al. (2009) value of $Z_{\odot}=0.0142$ more closely matches modern estimates of the protosolar abundance.

For stars in young moving groups, this method tends to overestimate the ages of these higher-mass stars, due to the fact that there is essentially no movement across the HR diagram during the first third of a star’s main sequence life. For example, the absolute $V$ magnitude of a 1.8 $M_{\odot}$ star will change less than 0.05 mags from the ZAMS ($\sim$15 Myr) until $\sim$240 Myr, given the models of Siess et al. (2000). Thus $\beta$ Pictoris, despite being the namesake of the $26\pm 3$ Myr $\beta$ Pic moving group, has a poorly constrained age based on its color-magnitude diagram position alone, with a 1-$\sigma$ upper limit of $<345$ Myr, and $<510$ Myr at 2-$\sigma$. This is similar to the age derived by Nielsen et al. (2013) of $110^{+20}_{-30}$ Myr, and $520^{+400}_{-170}$ Myr from David & Hillenbrand (2015).

Thus, while there is variation between the different literature age determinations using this method, the precise set of models used and especially assumed protosolar metallicity affects the final shape of the age posterior. Nevertheless, the ages used for the higher-mass stars in our sample not in moving groups are derived in a consistent manner, rather than relying on a heterogeneous compilation of literature estimates.

2.1.2 Lower-mass stars: masses

The above technique was used for stars bluer than $B-V=0.35$ (earlier than F1). We used a similar framework to estimate the masses of the lower-mass stars within the sample. Rotation of these stars can safely be ignored; observed $v\sin i$ values are significantly lower than for early-type stars (Gallet & Bouvier, 2013), and as such there is a negligible dependence of the effective temperature and surface gravity on stellar latitude. For these stars, the emcee sampler was used to sample the posterior distributions of star mass $M$, age $t$, metallicity $[M/H]$, and parallax $\pi$. As for the higher-mass stars, an additional mass term was included in the fit for each stellar companion within known spectroscopic binary systems. A prior on age was applied for these stars, either for the group age for members of moving groups or associations, or the age distribution derived by activity indicators. At each step the effective temperature, surface gravity, and radius of the star was calculated by interpolating the evolutionary model. The flux from the star in each bandpass was then computed from an interpolation of the ATLAS9 model atmosphere grid, scaled by the square of the ratio of the radius of the star and the distance. These lower-mass stars not in moving groups or associations are assigned a gaussian age uncertainty with $\sigma$ of 20%. We are in the process of deriving more realistic age uncertainties for individual FGKM stars based on calcium emission, X-ray emission, and lithium absorption, and will present updated age and mass posteriors in a future paper. For this work, then, we assume a single value of age and mass for each star when deriving completeness. For a future paper detailing the final description of the GPIES sample, however, we will incorporate full posteriors in age and mass.

Extinction is not treated for either the higher-mass or lower-mass stars due to the relative proximity of the stars within the Local Bubble. The overall sample is also relatively free of binaries by construction, although there are undoubtedly spectroscopic binaries that remain undiscovered in the sample. The effect of a binary companion would be to make a star appear more luminous and redder, and therefore typically appear older within our analysis. An overestimate of the age of the star would lead to an underestimate of the sensitivity of the observations of that star to planetary-mass companions, and an overestimate in the overall planet frequency in the final analysis. This bias would be most significant for unresolved binaries with a near equal flux ratio, the type of systems that are most amenable for detection via radial velocity. The effect of this bias on the analysis presented below is expected to be small given that the B and A-type stars with ages determined using this technique only constitute 16% of the sample, and only a small fraction of them are expected to be undiscovered equal-mass spectroscopic binaries.

There are a total of 12 spectroscopic binaries remaining in the sample and their properties and estimates of the component masses from the SED fit described previously are given in Table 1. There are several physical binaries and multiples wider than $3\arcsec$ (the limit used in constructing the sample), for example GJ 3305 orbits 51 Eridani at $66.7\arcsec$ (Feigelson et al., 2006). We expect these very wide systems to have less of a dynamical influence on the disk and planet formation. A more detailed analysis of the effect binarity has on planet formation requires a larger, dedicated sample that is complete to both single stars and binaries, with early work on this front being carried out by the SPOTS survey (Asensio-Torres et al., 2018). In this work, we consider very close SBs and wide-orbit binaries the same as single stars, though for SBs we adopt the combined mass of the system as the stellar host mass.

2.2 Data Acquisition

GPIES has been operating in priority visitor mode, where members of the GPIES team operate GPI to take campaign observations, since the start of the campaign in late 2014. An automated target-picker was used to select the best star to observe given the ages and distances of the remaining unobserved stars, and the current hour angle and environmental conditions (McBride et al., 2011). Targets with known substellar companions were not prioritized. The instrumental setup was similar for each target. A typical sequence consists of 38 one-minute coronagraphic exposures taken with GPI’s integral field spectrograph at a low spectral resolution ($\lambda/\Delta\lambda=50$) in the H-band (1.5–1.8µm). The Cassegrain rotator was fixed for all GPI observations causing astrophysical signals to rotate in the IFS images as the target being observed transits overhead. Observations were conducted under program codes GS-2014B-Q-500, GS-2014B-Q-501, GS-2015A-Q-500, GS-2015A-Q-501, GS-2015B-Q-500, and GS-2015B-Q-501. Observations in 2016 were conducted under the 2015B semester program codes. An observing log is given for the first 300 stars observed in the campaign—and any subsequent observations taken to confirm or reject candidate companions—in Table A.

3 Data reduction and contrast curves

GPIES data are processed automatically, with the reduction steps described in detail in Wang et al. (2018a); we briefly describe them here. Each frame of integral field spectroscopy data from GPI is processed with the GPI Data Reduction Pipeline (DRP; Perrin et al., 2014, 2016). The GPI DRP produces dark frames and wavelength solutions for each lenslet using calibration data taken during the daytime. Immediately before each observing sequence, a snapshot argon arc lamp is taken and is processed by the GPI DRP to measure the instrument flexure between the observing sequence and the daytime calibration data. For each science frame, the GPI DRP subtracts the dark frame, corrects for bad pixels, compensates the wavelength solution for flexure, constructs 3D spectral data cubes, corrects for optical distortion, and measures the locations and fluxes of the satellite spots that are used for calibration.

A typical observing sequence consists of 38 one-minute exposures, each of which is reduced into a spectral datacube with 37 wavelength channels, resulting in a total of 1406 individual spectral slices, with some correlation between adjacent wavelength channels. The stellar PSF (i.e., speckles) is subtracted from each of these spectral slices. With GPI, we are able to utilize both angular differential imaging (ADI; Marois et al., 2006) and spectral differential imaging (SDI; Racine et al., 1999; Marois et al., 2000; Sparks & Ford, 2002) to disentangle the potential planets from the speckles. The speckle subtraction and the planet detection are performed using an open-source Python package pyKLIP (Wang et al., 2015), which includes an implementation of the Karhunen-Loève (K-L) Image Projection algorithm (KLIP; Soummer et al., 2012; Pueyo et al., 2015). We follow the forward model matched filter approach described in Ruffio et al. (2017b) to which a few improvements, which are described in the following, were made.

The first step of the reduction consists in removing the time-averaged uncorrected atmospheric speckles, which result in a smooth stellar halo steeply rising near the edge of the coronagraphic mask. When the wind is strong, the halo becomes lopsided and aligned with the direction of the wind, sometimes denoted wind-butterfly. Consequently, the halo cannot be efficiently subtracted with a classical spatial high-pass filter and a Gaussian kernel, which was the approached taken in Ruffio et al. (2017b). In this work, we first subtract the radial intensity profile in the image, which removes the symmetric component of the stellar halo. Then, the lopsided component is subtracted from its projection onto a set of principal components computed from a library of GPI images containing a strong wind-butterfly.

The second processing step consists in subtracting the speckles individually in small sectors of each image. The results presented in Ruffio et al. (2017b) suffered from artifacts related to the edge of the sectors. In this work, we used tight overlapping sectors centered at each separation in the image therefore mitigating the edge-effects.

Finally, flux, standard deviation and Signal-to-Noise ratio (S/N) maps are computed according to Ruffio et al. (2017b) using a matched filter accounting for the heteroskedasticity (i.e., non-uniformity) of the noise and the bias in the planet signal induced by the speckle subtraction (Pueyo, 2016). The matched filter assumes a specific spectrum for the planet; however, it was demonstrated that two templates (a cool T-type and a warm L-type spectral template) were sufficient to recover most spectral types (See Ruffio et al., 2017a). The entire survey was therefore reduced twice, once for each template.

The speckle subtraction has several free parameters that were optimized in Ruffio et al. (2017a) for GPI. We used 30 K-L modes to model the stellar PSF. For each subtraction zone, the K-L modes were generated from the 150 most correlated reference images selected from the subset of the observing sequence in which a planet in that zone would move by at least 0.7 pixel in H band due to a combination of ADI and SDI

Planet detection thresholds are defined from the standard deviation map after calibration (in unit of planet-to-star flux ratio). The standard deviation map calibration is described in Ruffio et al. (2017b). It guarantees that the S/N map has a standard deviation of one and corrects for any algorithm related flux losses (i.e., algorithm throughput) using extensive planet injection and recovery. Then, the 1D standard deviation as a function of separation is computed as the root-mean-square of the standard deviation map in 4-pixel concentric annuli. Detection thresholds are generally derived from a false positive rate chosen a priori. However, the estimation of planet occurrence rates, the goal of this work, requires that all the candidates above the detection threshold be vetted. The vetting process generally requires follow-up observations, which is a limited resource due to telescope time constraints. As a result, we define the detection threshold of $\sigma$ as the boundary used over the GPIES campaign for candidate follow-up. This 1D detection threshold is the contrast curve. Unfortunately, while the stellar halo subtraction generally improved GPIES’ sensitivity at intermediate separations, it introduced a large number of false positives at the very edge of the coronagraphic mask when using the T-type spectral template. We therefore increased the detection threshold to $8\sigma$ below $<0.3^{\prime\prime}$ in these reductions. Full contrast curves for the entire GPIES campaign, including the subset discussed here, will be presented along with the final statistical analysis of the campaign in a future paper.

4 Companion Detections

Nine substellar companions around seven host stars have been detected within the first half of the GPIES survey, three brown dwarfs and six planetary-mass companions. The measured position, contrast, and derived mass for each companion is given in Table 2. The first GPIES campaign image and corresponding sensitivity curves are given in Figure 2 for the imaged planets, and Figure 3 for the brown dwarfs. 51 Eridani b and HR 2562 were new discoveries from these images, while $\beta$ Pictoris b, HD 95086 b, HR 8799 cde, HD 984 B, and PZ Tel B were detections in our campaign images of known substellar companions.

4.1 51 Eridani

51 Eridani b was discovered early in the GPIES campaign (Macintosh et al., 2015) around the F0IV $\beta$ Pictoris moving group member (Bell et al., 2015) 51 Eridani (hereafter 51 Eri). The H-band spectrum measured in the discovery epoch exhibited strong methane and water absorption, consistent with both predictions for low-mass exoplanets around young stars, and also the observed spectra of low-temperature T-type brown dwarfs. Follow-up astrometric observations confirmed that 51 Eri b was in a bound orbit around 51 Eri (De Rosa et al., 2015), with a negligible probability of it being an interloping field T-dwarf in chance alignment with 51 Eri.

Subsequent spectroscopic observations with GPI (Rajan et al., 2017) and VLT/SPHERE (Samland et al., 2017), combined with thermal-IR photometry obtained with Keck/NIRC2 (Macintosh+15, Rajan+17), have been used to characterize the atmospheric properties of the planet. The observed spectral energy distribution is consistent with a photospheric cloud fraction between 75-90% (Rajan et al., 2017) and 100% (Samland et al., 2017), depending on the specifics of the cloud model used in the analysis. Both fits find a similar effective temperature ($T_{\rm eff}=$ 600–760 K), but a range of values for the surface gravity ($\log g=$ 3.5–4.5 [dex]). The estimated mass of the planet is strongly dependent on the assumptions made regarding its formation. Under an assumption of a high-entropy (so-called ‘hot-start’) formation scenario (Marley et al., 2007), the luminosity of the planet corresponds to a mass of 1–2 MJup (Rajan et al., 2017). Uniquely for young directly-imaged companions, the luminosity of 51 Eri b is also consistent with predictions of low-entropy, ‘cold-start’ formation scenarios (Marley et al., 2007; Fortney et al., 2008), with a significantly higher model-dependent mass of 2–12 MJup. Measurements of the reflex motion induced by the orbiting planet, either via spectroscopy or astrometry, are required to differentiate between these two scenarios.

4.2 $\beta$ Pictoris

The exoplanet $\beta$ Pictoris b ($\beta$ Pic b) was discovered around the A6V namesake of the $\beta$ Pictoris moving group using VLT/NaCo in 2003, later confirmed to be a bound companion in follow-up observations taken with the same instrument in 2008 (Lagrange et al., 2009, 2010). The planet is orbiting interior to a large circumstellar disk (Smith & Terrile, 1984), and has a measurable effect on the morphology of the inner disk (Lagrange et al., 2012; Millar-Blanchaer et al., 2015). $\beta$ Pic b was observed several times with GPI; as a first light target during commissioning, as a part of an astrometric monitoring program to constrain the orbit of the system, and finally as part of the GPIES Campaign. Due to its proximity and youth, the primary was a highly ranked GPIES target, and was selected by our automated target-picker early in the campaign, resulting in the detection of $\beta$ Pic b. Milli-arcsecond astrometry obtained with GPI was critical for constraining the orbital parameters of the planet (Wang et al., 2016). Observed variability in the light curve of the host star $\beta$ Pic was interpreted as being caused by a transit of either the planet or by material within its Hill sphere (Lecavelier Des Etangs & Vidal-Madjar, 2009). Wang et al. (2016) used GPI astrometry to tightly constrain the inclination of the orbit of the planet and excluded a transit at the 10-$\sigma$ level, although a transit of the Hill sphere of the planet was still a possibility. Precise timings on the transit of the planet’s Hill sphere were made, allowing for a unique opportunity to probe the circumplanetary environment of a young exoplanet.

The low contrast between $\beta$ Pic b and its host star compared to other directly-imaged exoplanets makes it an ideal target for atmospheric characterization. Chilcote et al. (2017) used a combination of GPI spectroscopy and literature photometry (compiled and calibrated by Morzinski et al., 2015) to measure the spectral energy distribution of the planet between 0.9 and 5 µm, wavelengths spanning the bulk of the emergent flux of the planet. The data were consistent with models that incorporated the condensation of dust into thick vertically-extended clouds within the photosphere, and matched the observed spectral energy distributions of isolated low-surface gravity brown dwarfs (Chilcote et al., 2017). Unlike 51 Eri b, the bolometric luminosity of $\beta$ Pic b is only consistent with the predictions of the ‘hot-start’ luminosity models, and inconsistent with the Fortney et al. (2008) ‘cold-start’ models. However, this is not a criticism of core accretion; see, e.g., Berardo et al. (2017); Owen & Menou (2016) for how core accretion can lead to hot/warm starts, and also our Section 6.3. Chilcote et al. (2017) derived a model-dependent mass of $12.7\pm 0.3$ MJup, consistent with the model-independent mass derived from Hipparcos and Gaia astrometry of the host star of $11\pm 2$ MJup (Snellen & Brown, 2018) and 13$\pm$3 MJup (Dupuy et al., 2019).

4.3 HD 95086

HD 95086 b was discovered around the A8III member of the 15 Myr (Pecaut & Mamajek, 2016a) Lower Centaurus Crux subgroup of the Scorpius-Centaurus OB2 association (Rizzuto et al., 2011) using VLT/NaCo in 2011 (Rameau et al., 2013b, a), one of a growing number of exoplanets that have been resolved around Sco-Cen members (Bailey et al., 2014; Chauvin et al., 2017). The planet was noted for its unusually red infrared colors that were ascribed to a photosphere dominated by thick clouds (Meshkat et al., 2013). HD 95086 was observed as a part of the GPIES campaign to characterize the atmosphere of the planet, monitor its orbital motion, and to search for additional interior companions. Like $\beta$ Pic, the target star was selected for observations for our planet-search program by our automated target-picker, resulting in a detection of the planet. The low luminosity and very red infrared colors of the planet combined with its distance of 86 pc make it a challenging target for spectroscopic observations, especially at shorter wavelengths. De Rosa et al. (2016) reported the first spectroscopic measurement of the planet obtained using GPI, a spectrum covering the blue half of the $K$ band. This was combined with literature photometry at $J$ and $L^{\prime}$ to reveal a spectral energy distribution consistent with model atmospheres incorporating significant amounts of photospheric dust, a result confirmed in a later analysis incorporating SPHERE observations of the planet (Chauvin et al., 2018).

Monitoring the orbital motion of the planet can provide insight into the architecture of the HD 95086 system. The presence of two large circumstellar dust rings was inferred from the spectral energy distribution of the star (Moór et al., 2013; Su et al., 2017); the outer ring was subsequently resolved with millimeter interferometric observations with ALMA (Su et al., 2017). HD 95086 b lies between these two rings and, despite only sampling a small fraction of the full orbit of the planet with GPI and literature astrometry, Rameau et al. (2016) were able to constrain the eccentricity of the orbit to be near-circular, under the assumption of co-planarity with the outer ring. While the dynamical impact of the planet on the outer ring is still uncertain (Su et al., 2017), Rameau et al. (2016) demonstrated that it is unlikely that the planet is responsible for clearing the gap in to the predicted outer radius of the inner disk, and posited this as indirect evidence of the existence of additional lower-mass companions that are below present detection limits, a finding also confirmed by Chauvin et al. (2018).

4.4 HR 8799

The HR 8799 system is currently the only example of a multiple planet system detected via direct imaging (Marois et al., 2008, 2010; Konopacky & Barman, 2018). The host star is a $\lambda$ Bootis star that exhibits $\gamma$ Doradus variability (Zerbi et al., 1999; Gray & Kaye, 1999). It is a probable member of the $42^{+6}_{-4}$ Myr (Bell et al., 2015) Columba moving group (Zuckerman et al., 2011; Malo et al., 2013) and has a spectral type of either A5 and F0, depending on which set of spectral lines are used for typing (Gray et al., 2003). The star is also host to two circumstellar dust rings, at $\sim$10 and $\sim 100-300$ au, with the planets occupying the gap between them, and a large diffuse dust halo at $\sim$300-1000 au (Su et al., 2009). While the inner ring is only inferred from the spectral energy distribution of the star, the outer ring has been resolved with ALMA and the SMA (Booth et al., 2016; Wilner et al., 2018).

The HR 8799 system has been observed multiple times with GPI; as an early science target during the commissioning of the instrument (Ingraham et al., 2014), and during the campaign to both characterize the atmospheres of the known planets (Greenbaum et al., 2018) at multiple bands, and as part of the planet search in $H$ when chosen by the automated target-picker. The campaign observations resulted in the detection of three of the four planets, HR 8799 cde, with HR 8799 b outside the field of view of GPI in our nominal search mode. Astrometry obtained with GPI was also used in conjunction with measurements obtained with Keck/NIRC2 to investigate the dynamical properties and long-term stability of the system (Wang et al., 2018b). By utilizing N-body simulations, the combination of astrometry and dynamical limits placed tighter constraints on the orbital parameters and the masses of the individual planets, which when combined with evolutionary models give a mass for HR 8799 b of $5.8\pm 0.5$ MJup, and $7.2^{+0.6}_{-0.7}$ MJup for planets c, d, and e.

4.5 Brown Dwarf Companions

In addition to the four planetary systems discovered during the first half of the GPIES campaign, three brown dwarf companions were resolved; PZ Tel B, HD 984 B, and HR 2562 B. The companions to PZ Tel (Mugrauer et al., 2010; Biller et al., 2010) and HD 984 (Meshkat et al., 2015) were known before the GPIES observations of these stars, and were detected when their host stars were selected by our automated target-picker, while HR 2562 B was first discovered with GPI (Konopacky et al., 2016) as part of our planet search campaign. PZ Tel and HD 984 are probable members of the $\beta$ Pictoris and Columba moving groups, respectively (Malo et al., 2013), while HR 2562 is thought to be $<$1 Gyr based on photometric (Casagrande et al., 2011) and spectroscopic (Pace, 2013) analyses. The three host stars have a later spectral type than the planet host stars discussed previously; G9IV (1.1 M*⊙) for PZ Tel, F7V (1.3 M⊙) for HD 984, and F5V (1.3 M⊙*) for HR 2562.

A brown dwarf companion was later detected around the campaign target HD 206893 (initially discovered with the VLT/SPHERE by Milli et al., 2016), but these observations were obtained after the 300th star was observed and the companion is therefore not included in any of the analyses presented below for the first half of the GPIES campaign.

5 Survey Modeling and Inference

Having observed 300 stars and detected six planets and three brown dwarfs, GPIES is one of the largest, deepest direct imaging surveys, with a high yield for substellar companions. Using this rich dataset, we proceed to measure the overall occurrence rate of giant planets and trends in the underlying population.

5.1 Survey Completeness to Date

We determine the completeness to planets for each of our observations using the Monte Carlo procedure described in Nielsen et al. (2008), Nielsen & Close (2010), and Nielsen et al. (2013). For each target star, completeness to substellar companions is determined over a grid that is uniform in log space in mass and semi-major axis. At each grid point, $10^{4}$ simulated substellar companions, each with the same value of mass and semi-major axis, are randomly assigned orbital parameters: inclination angle from a $\sin(i)$ distribution, eccentricity from a linear distribution fit to wider-separation RV planets (Nielsen et al., 2008), $P(e)=2.1-2.2e$ with $0\leq e\leq 0.95$, and argument of periastron and epoch of periastron passage from a uniform distribution. Since contrast curves are one dimensional, position angle of nodes is not simulated. The projected separation (from the orbital parameters and mass and distance of the host star) and the $\Delta H$ (from the mass of the companion, the age and absolute magnitude of the host star, and luminosity models) for each companion are then computed. A simulated companion that lies above the contrast curve is considered detectable, while one that lies below (or outside the field of view) is undetectable. If more than one contrast curve is available for a single star, the simulated companions generated at the first epoch are advanced forward in their orbits to the next epoch, and compared to the next contrast curve; a companion is detectable if it is above at least one contrast curve.

As noted above, for the forward model matched filter (FMMF) contrast curves described in Section 2.2, we use an 8$\sigma$ threshold inside of 0.3*′′, and a 6$\sigma$ threshold outside. FMMF contrasts were produced for both a T-type and L-type template for each star (Ruffio et al., 2017b). The T-type reduction assumes significant methane absorption in the $H$ band, so that images in redder channels may be subtracted from images in bluer channels without subtracting the planet from itself. To stay consistent, the same evolutionary models that produce $H$ band flux are used to give temperature as a function of mass and age, and planets hotter than 1100 K are assumed undectable using the T-type curve. The L-type curve, which does not suffer from the possibility of self-subtraction, is used for simulated companions of all temperatures. These FMMF contrast curves are defined out to 1.7′′, though beyond 1.1′′* some parts of the field go beyond the edge of the detector, and observations are not complete over all position angles. We account for this with the same method as Nielsen et al. (2013), by recording the fractional completeness as a function of separation, and using a uniform random variable to reject simulated planets that would not fall on the detector.

We utilize the CIFIST2011 BT-Settl atmosphere models (Baraffe et al., 2015; Allard, 2014; Caffau et al., 2011)111https://phoenix.ens-lyon.fr/Grids/BT-Settl/CIFIST2011/COLORS/colmag.BT-Settl.server.2MASS.Vega, tied to the COND ‘hot start’ evolutionary model grid (Baraffe et al., 2003), to give absolute $H$ magnitude and temperature as a function of mass and age. From these, we produce completeness maps for each star. These maps give the fractional completeness to companions as a function of mass and semi-major axis (completeness at a given value of mass and semi-major axis is not binary since orbital parameters will set whether a given companion is detectable for a certain contrast curve). These completeness maps are shown for stars with detected companions in Figures 2 and 3. Summing the maps across all target stars produces a map representing our depth of search (e.g. Lunine et al. 2008) for the survey, as given in Figure 4. The contours give the number of stars to which the survey is complete to companions of a given mass and semi-major axis. Since target stars cover a large range of distances (3.2–176.7 pc), the limited field of view of GPIES ($\sim$0.15–1.1*′′, beyond which not all position angles lie on the detector) means that completeness plots for individual stars do not precisely line up. Individual completeness plots never reach 100% for a similar reason, as some simulated planets at each semi-major axis can fall within the IWA or outside the outer working angle, and so go undetectable. As such the maximum value reached on the plot is 246 stars, though all 300 stars have some fractional sensitivity to substellar companions. The detected companions are indicated by red dots, plotted at their inferred mass and projected separation at the first GPIES epoch, as given in Table 2. In this figure, and in the subsequent analysis, we assume masses for each companion inferred from the BT-Settl models, the age of the system, and the $H$-band magnitude, even if a more robust mass is available from full SED fitting, to maintain consistency in the completeness calculation. Most of these assumed masses are close to the published mass, with the notable exception of HD 95086, where the $H$-band magnitude and age suggest a mass of 2.6 $\pm$ 0.4 MJup*, compared to 4.4 $\pm$ 0.8 MJup inferred from the $Ks$-band magnitude and a bolometric correction (De Rosa et al., 2016). HD 95086 is a particularly unusual case, as its dusty atmosphere leads to a redder spectrum than most models predict. We note that the analysis below does not change significantly whether a mass of 2.6 or 4.4 MJup is used.

We also investigated the effect of utilizing full age posteriors for each star, rather than using a single age for each target. We considered the AB Dor moving group target star AN Sex, and computed the completeness plot for a single age of 149 Myr (as used in this work), for a disjoint gaussian, and for a symmetric gaussian. The disjoint gaussian has a $\sigma$ of 19 Myr below 149 Myr, and a $\sigma$ of 51 Myr above, to match the Bell et al. (2015) age for the AB Dor moving group of 149${}^{+51}_{-19}$ Myr. The symmetric gaussian, computed for comparison, is given a $\sigma$ of 20 Myr (13% age uncertainty). As expected, assuming younger ages result in more detectable planets compared to older ages, however the effect on the completeness plot is marginal. The 20% completeness contour, for example, moves at most 7% in planet mass, with a median shift of 1%, for the disjoint gaussian. For the symmetric gaussian, the effect is even smaller, maximum of 5% and median 0.2%. In the symmetric case, the contours move toward smaller planets (becoming more sensitive when using a symmetric gaussian age posterior compared to a single age), while the completeness plot is less sensitive for the disjoint gaussian. Similar results are found for the symmetric gaussian even if the value of $\sigma$ is doubled. We thus expect a minor effect on our results by including age posteriors for each target star, but this effect will be explored more in depth in a future paper on the statistical constraints from the full survey, once full age posteriors are available for all of our stars.

Utilizing mass posteriors is expected to have an even smaller effect on the completeness plots. Since simulated planets are generated as a function of semi-major axis, stellar mass has no effect on the completeness for a star observed at only a single epoch (273 out of our 300 stars), as planets are randomly assigned a mean anomaly, which does not require the period to be known. For stars observed at multiple epochs, however, period is needed to determine the amount of orbital motion between the two observations, but given the slow orbital speeds of these wide-separation planets, and that the observations described here span less than two years, this is a minor effect as well, especially as we typically reach $\lesssim$5% precision on stellar mass on our targets.

5.2 A correlation between stellar mass and planet occurrence rate

A shared feature of the planetary companions described in Section 4 is that all orbit the higher-mass stars in our sample: all host stars lie between 1.55 and 1.75 M⊙ (Figure 5). As shown in Figure 1, there is no clear mass bias in our sample toward higher masses; in fact, 177 out of 300 stars (59%) have masses below 1.5 M⊙. This trend is also the opposite from what we would expect for observational biases: lower mass stars are intrinsically fainter, and thus the same achievable contrast would correspond to a lower mass detectable planet.

We investigate the strength of this trend by comparing occurrence rates of wide-separation giant planets around higher-mass stars and lower-mass stars, here defined to have a mass determination above and below 1.5 M⊙. We consider planets with semi-major axis between 3 and 100 au, and planet mass between 2 and 13 MJup. These limits are chosen to encapsulate a region of planet parameter space to which our observations are most sensitive and encompass the detected planets in the GPIES sample, which are between 2.6–12.9 MJup and a projected separation from 4.84–53.65 au (as shown in Figure 4). While the depth of search is computed for semi-major axis, the most direct measurement we have of each detected planet is projected separation, so we have chosen a broad enough range of semi-major axes so that planets at these projected separations are unlikely to fall outside this range. Indeed, orbital monitoring of each of these planets (De Rosa et al., 2015; Wang et al., 2016; Rameau et al., 2016; Wang et al., 2018b) place them solidly in this range. We also begin by assuming planets are uniformly distributed over this range in log space in both semi-major axis and planet mass ($\frac{d^{2}N}{da\,dm}\propto a^{-1}m^{-1}$). This is not too dissimilar to the power law distributions for small-separation ($a<3.1$ au) giant planets found by Cumming et al. (2008) of $\frac{d^{2}N}{da\,dm}\propto a^{-0.61}m^{-1.31}$, obtained by converting the period distribution of $\frac{dN}{dP}\propto P^{-0.74}$ to semi-major axis for solar-type stars. For each star, then, we integrate the completeness to planets over this range with uniform weight given to each log bin in semi-major axis and mass.

Integrating over the entire range from 3–100 au and 2–13 MJup produces a completeness to massive giant planets of 17.5 stars among the higher-mass ($>1.5M_{\odot}$) sample, and 43.3 among the lower-mass sample. These two samples have 4 and 0 planetary systems detected, respectively, as HR 8799 counts as a single planetary system. We infer the frequency of planetary systems for each subsample using a Bayesian approach with a Poisson likelihood ($L=\frac{e^{-\lambda}\lambda^{k}}{k!}$). The expected number of planetary systems ($\lambda$) is the number of stars to which the subsample is complete multiplied by the frequency, and the measured number ($k$) is the number of detected planetary systems. The prior is taken to be the Jeffrey’s prior on the rate parameter of a Poisson distribution, $\mathrm{Prior}(\lambda)\propto\sqrt{\frac{1}{\lambda}}$. We compute the posterior over a well-sampled regular grid of frequency for both sub-samples, plotting the results in Figure 6.

Given no planets were detected around lower-mass stars ($<1.5$ M⊙), our measurement of the frequency for these stars represents an upper limit, $<$6.9% at 95% confidence. With four detected systems, the fraction of higher-mass stars with giant planet systems is 24${}^{+13}_{-10}$% (68% confidence interval). To evaluate the confidence with which we can conclude that the higher-mass (HM) star frequency is larger than the lower-mass (LM) star frequency, we compute the posterior of the difference, $\delta=f_{HM}-f_{LM}$, given by

[TABLE]

where $P_{HM}$ and $P_{LM}$ are the posteriors on the frequency ($f$) for higher- and lower-mass stars. We confirm this expression with a Monte Carlo method, where a set of two frequencies are randomly generated from the two posteriors and differenced. The two methods generated identical results, as shown in the right panel of Figure 6. For 99.92% of the samples, wide-separation giant planets around higher-mass stars are more common than around lower-mass stars, a $>$3$\sigma$ result.

We note that this result is not strongly dependent on our choice of limits in planet mass, semi-major axis, and stellar mass. When increasing the upper semi-major axis limit from 100 au to 300 au, the significance of the result that wide-separation giant planets are more common around higher-mass stars drops marginally from 99.92% to 99.85%. Since we are averaging over the depth of search plot in Figure 4, expanding the range to areas of lower sensitivity reduces the effective number of stars to which the survey is sensitive, and for lower-mass stars this number drops from 43.3 to 36.2, and higher-mass stars drops from 17.5 to 17.4, when expanding the semi-major axis range from 100 to 300 au. Increasing this range still further from 3–100 to 1-1000 au reduces the effective number of lower-mass stars to 24.9 and higher-mass stars to 12.0, and the significance in this case remains above 3$\sigma$ at 99.84%.

Similarly, dropping the planet mass range considered from 2–13 MJup to 1–13 MJup has a negligible effect, with the significance remaining at 99.92%. Dropping the 12 SBs from the sample is also a minor effect, removing mainly higher-mass stars, and raising the significance slightly to 99.95%. Stellar mass has a more significant effect, since moving the boundary between higher-mass and lower-mass stars results in fewer lower-mass stars in the sample. The significance drops to 99.70% when the boundary is moved to 1.35 M⊙ (an effective number of lower-mass stars of 37.9), and 96.99% at 1.10 M⊙ (24.3 stars). Thus the larger frequency of wide-separation giant planets around higher-mass stars is relatively robust to specific choices of which range of parameter space we choose.

5.3 Wide Separation Giant Planet Occurrence Rate around Higher-Mass Stars

Following the example of Cumming et al. (2008), we consider a model of planet distributions defined by power laws in mass and semi-major axis. We adopt a functional form of our model of:

[TABLE]

where $m$ and $a$ are planet mass and semi-major axis. We define this equation over a limited range, $[m_{1}\leq m\leq m_{2}]$ and $[a_{1}\leq a\leq a_{2}]$. The frequency ($f$) is the occurrence rate of planets in this range, and so the normalization constant $C_{1}$ is then given by

[TABLE]

Thus, if $\alpha\neq-1$ and $\beta\neq-1$:

[TABLE]

otherwise the two terms become $[\ln(a_{2})-\ln(a_{1})]^{-1}$ and $[\ln(m_{2})-\ln(m_{1})]^{-1}$.

In order to constrain the three free parameters of this model ($f$, $\alpha$, and $\beta$) we adopt a Bayesian approach utilizing Monte Carlo completeness calculations. Bayes’ Equation then becomes:

[TABLE]

and in standard Bayesian terminology, the first term is the posterior, the second the likelihood, and the final three terms are the priors. We follow a similar method to Kraus et al. (2008), dividing the two-dimensional observational space (mass vs. projected separation) into a series of bins. Similar Bayesian methods were also used by Biller et al. (2013), Wahhaj et al. (2013), and Brandt et al. (2014) to fit a power law model to imaged planets. 41 bins in mass are used, logarithmically spaced between 1 and 100 MJup, and 81 bins in separation, logarithmically spaced between 1 and 1000 au. In each bin we calculate the expected and actual number of planets in that bin; as a result, we choose a Poisson distribution for the likelihood in each bin, then find the product of probabilities across all the bins:

[TABLE]

where $i$ and $j$ are indices for bins in mass and projected separation, and $E_{i,j}$ and $O_{i,j}$ are the expected and observed number of planets in each bin, respectively.

Unlike Kraus et al. (2008), our model and likelihood use different parameters: Equation 2 is in terms of mass and semi-major axis, while Equation 6 is in terms of mass and projected separation, since directly-imaged planets are detected with a particular separation, but an (often unknown) semi-major axis. Using Monte Carlo simulations, however, we are able to project our modeled planets from semi-major axis into separation. As in Section 5.1, we inject simulated planets at a grid of mass and semi-major axis, and for each grid point save the fraction of simulated planets that are above the contrast curve (detectable with our observations). Additionally, we also sort the detectable planets into projected separation bins, tracking the projected separation of each detectable planet. The completeness map for each star, then, is a three-dimensional array across mass, projected separation, and semi-major axis. Marginalizing over the separation dimension, and summing over all stars, gives the depth of search shown in Figure 4. If we instead marginalize over semi-major axis, we derive a completeness map in mass and projected separation, the same parameters that define our observed planets.

For a given set of ($f$, $\alpha$, $\beta$), we follow this procedure for evaluating Equation 5:

For each star, generate a three-dimensional completeness map using Monte Carlo simulations 2. 2.

Sum the completeness maps over all stars 3. 3.

Weight each point in the map using Equation 2 and the values of ($f$, $\alpha$, $\beta$), giving the expected number of planets detected in each bin 4. 4.

Marginalize the map over semi-major axis, producing a map of expected planets vs. mass and separation 5. 5.

Produce a corresponding observed map of detected planets vs. mass and separation 6. 6.

Compute the likelihood from Equation 6. 7. 7.

Multiply by priors

We then utilize a Metropolis-Hastings MCMC procedure (as in Nielsen et al. 2014) to explore this parameter space and generate posteriors on frequency and power law indices. For priors we utilize uniform priors for $\alpha$ and $\beta$, and a Jeffrey’s prior for $f$ appropriate to a Poisson distribution, $\mathrm{Prior}(f)\propto\sqrt{\frac{1}{f}}$. Our results are not overly dependent on this choice of prior, we find similar posteriors with a prior that is uniform in $\ln{(f)}$. We first restrict ourselves to a region of parameter space where we have both good completeness and our detections, $[1~{}$M${}_{\rm Jup}\leq m\leq 13~{}$M${}_{\rm Jup}]$ and $[1~{}{\rm au}\leq a\leq 100~{}{\rm au}]$, and M${}_{*}>$ 1.5 M⊙. This range encompasses all six planets detected by the GPIES survey. Given our low sensitivity to low-mass, short period planets, there is a strong degeneracy between frequency and power law indices, since very negative indices (that place large numbers of planets at small separations and small masses) result in larger frequencies. Fits with large occurrence rates place most of these planets in the region of parameter space to which GPIES is least sensitive. In order to report an occurrence rate with high confidence, we define occurrence over a region of high completeness, and so we quote not the frequency over the entire range, but instead a modified frequency over the truncated range $[5~{}$M${}_{\rm Jup}\leq m\leq 13~{}$M${}_{\rm Jup}]$ and $[10~{}{\rm au}\leq a\leq 100~{}{\rm au}]$. Thus while 51 Eridani b, HD 95086 b, and $\beta$ Pic are excluded from this modified frequency, they fall within the broader range, and thus are all used to constrain the model. The resulting posteriors are shown in Figure 7, and Table 3.

We find a relatively high giant planet occurrence rate: 9${}_{-4}^{+5}$%, for stars with mass $>1.5M_{\odot}$, and planets in the range $5~{}M_{\rm Jup}\leq m\leq 13~{}M_{\rm Jup}$ and $10~{}{\rm au}\leq a\leq 100~{}{\rm au}$. Given the very negative values for $\alpha$ and $\beta$, if we were to extrapolate beyond the region of high completeness to lower planet masses and smaller semi-major axes, we would infer a much larger occurrence rate with significantly larger uncertainties. Vigan et al. (2012) combined the IDPS sample of 42 stars of AF spectral types with three literature stars, and assuming values of $\alpha=\beta=0$, found a fraction of high-mass stars with at least one 3–14 MJup planet between 5–320 au of 8.7${}^{+10.1}_{-2.8}$%. Nielsen et al. (2013) examined the 70 B and A stars observed by the Gemini NICI Planet-Finding Campaign, and with a null result for planets placed 2$\sigma$ upper limits on the fraction of high-mass stars with planets. These limits have $<10$% of high-mass stars having a 10–13 MJup planet with semi-major axis between 38–650 au, and $<20$% with a 4–13 MJup planet between 59–460 au. Rameau et al. (2013) considered 37 AF stars, including detections of planets around HR 8799 and $\beta$ Pic, finding a planet fraction of 7.4${}^{+3.6}_{-2.4}$% from 5–320 au, between 3–14 MJup, assuming $\alpha=\beta=0$. Bowler (2016) considered multiple direct imaging surveys, performing a meta-analysis on 384 unique target stars. This sample was then broken up by spectral type, into subsamples of M, FGK, and BA. For the BA sample (roughly equivalent to our $>$1.5 M⊙ sample), by assuming $\alpha=\beta=-1$, Bowler (2016) find a planet fraction of 5–13 MJup 30–300 au planets of 2.8${}^{+3.7}_{-2.3}$%.

We compare our results to large imaging surveys with detected planets that explicitly measure occurrence rate in the top left panel of Figure 8, where we plot each survey as a line showing the semi-major axis distribution. End points are taken to be the range of semi-major axes considered by each survey. We utilize the value of $\alpha$ from each survey to find the expected occurrence rate over the GPIES mass range (5-13 MJup). In most cases authors assumed a value of $\beta$ (and $\alpha$) of $-1$ or 0, as we did in Section 5.2, rather than fitting for the power law index directly. Vigan et al. (2012) and Rameau et al. (2013) find similar results, as expected given both had similar numbers of stars and detected planets. The Bowler (2016) meta-analysis included these surveys and others, for a total of 110 BA stars, and finds a lower occurrence rate, and assumes a distribution that is flat in log space, rather than one flat in linear space. In all cases, the GPIES result is similar in fraction of planets at $\sim$100 au, but rises more steeply toward lower semi-major axis. This is as we expect, given our finding of a more negative value of $\beta$ than $-1$ or 0, placing more planets close to the star.

We also plot possible extrapolations of RV occurrence rates for giant planets. Johnson et al. (2010a) reported a value of 20% planet fraction for retired A-stars, for semi-major axis less than 2.5 au, and K$<$20 m/s (which they note corresponds to a lower mass limit of 1.12 MJup for a stellar host of 1.6 M⊙). Johnson et al. (2010a) do not fit for $\alpha$, and do not give the upper mass for planets in their sample. We assume both are similar to the Cumming et al. (2008) sample for planets around FGK stars, and use the latter values for upper mass limit of 10 MJup and $\alpha=-1.31$, in order to convert the occurrence rate for all planets to 5–13 MJup. Similarly, while Johnson et al. (2010a) give an upper limit of 2.5 au, they do not state the minimum semi-major axis of planets in their sample, though we speculate that this minimum value is 0.08 au, for HD 102956 b (Johnson et al., 2010b). As there is no estimate for the value of $\beta$ for RV planets around higher-mass stars, we investigate three different values, flat in log space, the Cumming et al. (2008) value, and flat in linear space ($\beta$ of -1, -0.61, and 0, respectively), and plot these as gray lines in Figure 8.

For integrated occurrence rate, the Johnson et al. (2010b) value of 20% ($<$2.5 au) represents a rise from our value of 11.4% (10–100 au), though for different mass ranges. A more detailed study of the distributions of RV-detected giant planets around higher-mass stars is required to definitively determine the nature of the distribution of high mass planets within 10 au, and whether RV and imaged giant planets around higher-mass stars follow a continuous distribution, but with a break in the planet mass and semi-major axis power law indexes at $\sim$3 au.

5.4 Occurrence Rate with a Stellar Mass Dependence

We can instead fit the entire GPIES 300-star sample by including a term that allows the occurrence rate to vary as a function of stellar mass. We thus introduce an additional term into Equation 2:

[TABLE]

so that the overall occurrence rate $f$ now varies with stellar mass. We choose $1.75M_{\odot}$ as the normalization location given this is where most of our planet-hosting stars lie, and so occurrence rate is most constrained at this mass. Similarly, Equation 5 becomes:

[TABLE]

We solve this equation as before, but with one additional change: instead of summing the completeness map over all stars, we instead introduce an extra dimension in stellar mass, and bin the sample by stellar mass. This produces a four-dimensional completeness map, as a function of companion mass, companion separation, companion semi-major axis, and stellar host mass. Similarly, the map of expected planets is now over three dimensions, companion mass, companion separation, and stellar host mass.

The results from the MCMC fit are plotted in Figure 9. As expected from the prevelance of detected planets around higher-mass stars, the distribution for $\gamma$ is peaked at relatively high values, $\gamma=2.0\pm 1.0$. Interestingly, though likely coincidentally, this value of $\gamma$ is close to the negative of the Salpeter IMF power law index of -2.35 ($\frac{dN_{*}}{dM_{*}}\propto M_{*}^{-2.35}$, Salpeter 1955), which would mean that for a volume-limited sample of young stars, while the number of giant planets per star would increase with increasing stellar mass, the number of giant planets per stellar mass bin would only change slightly, decreasing as the cube root of the stellar mass ($\frac{dN_{p}}{dM_{*}}\propto M_{*}^{-0.35}$). When we compare this fit to our earlier fit to only higher-mass stars (Figure 7), the planet mass and semi-major axis power law indexes, $\alpha$ and $\beta$, are similar for both fits. Occurrence rate, $f$, is significantly lower for the fit utilizing a stellar mass dependence. While a more positive value of $\gamma$ ($\sim$4) better explains the additional 177 lower-mass stars without planets in our sample, such a distribution begins to overpredict the number of detectable planets around the highest mass stars in our sample ($\gtrsim$ 2 M⊙). Thus the best fit is found by combining a moderate value of $\gamma$ ($2.0\pm 1.0$) and a lower overall occurrence rate. While this fit incorporates the entire 300-star sample, we find the fit exclusively to the 123 higher-mass stars presented above to be more robust, as it centers around the detected planets. Multiple models could account for the stellar mass dependence of giant planet occurrence seen in Figure 5. Here we have introduced a scaling of overall occurrence rate with stellar mass, but such a result could also be explained by a more negative power law fit to either companion mass or semi-major axis for lower-mass stars compared to $>1.5$M⊙ stars, or some combination of varying power-law indices and occurrence rate as a function of stellar mass. With no detections around our lower-mass star sample, we are unable to determine which model best represents these wider-separation giant planets.

Galicher et al. (2016) described statistical results from the IDPS survey of 292 stars and two previous surveys, analyzing a sample of 356 stars with spectral types between B and M. Assuming values of $\alpha=\beta=0$, a planet fraction was inferred of 1.05${}^{+2.80}_{-0.70}$% between 0.5–14 MJup and 20–300 au, based on detections of planetary systems around HR 8799 and HIP 30034 (AB Pic). An analysis of the VLT/NACO large program by Vigan et al. (2017) found similar results to IDPS with similar assumptions and analysis methods. The survey observed 85 stars with NACO, which was then combined with previously published samples to reach a final sample of 199 FGK stars. This final sample contained two detections of brown dwarfs, GSC 08047-00232 B and PZ Tel B, along with AB Pic b. With $\alpha=\beta=0$, the occurrence rate of 0.5–14 MJup, 20–300 au planets was found to be between 0.85 and 3.65% at 1$\sigma$, with a best value of 1.15%. From the Bowler (2016) meta-analysis, assuming $\alpha=\beta=-1$, an occurrence rate of 0.6${}^{+0.7}_{-0.5}$% was derived for planets between 5–13 MJup and 30–300 au around stars between B and M spectral types.

The top right panel of Figure 8 compares direct measurements of the planet occurrence rates to our findings for GPIES, where again our steeper slope places more planets at smaller semi-major axis. We have used our median value of $\gamma$ to convert these occurrence rates from 1.75 to 1 M⊙, to better match literature values. As before, we use the values of $\alpha$ assumed from each reference to convert the given mass ranges to the GPIES range of 5–13 MJup. The occurrence rates estimates for different imaging surveys give about the same value at $\sim$100 au. Cumming et al. (2008) and Fernandes et al. (2018) analyze RV surveys and find values of $\beta=-0.61\pm 0.15$ and $\beta=-0.025^{+0.30}_{-0.22}$, respectively, at small separations. While Cumming et al. (2008) fit a single power law, Fernandes et al. (2018) find evidence for a broken power law, with a break at $\sim$3 au. At these wider separations the Fernandes et al. (2018) value of $\beta=-1.975^{+0.22}_{-0.30}$ (from their symmetric epos model) is a good match to our value of $\beta=-1.68^{+0.57}_{-0.50}$, though we note that while these RV surveys contain mainly FGK stars, all our planet detections are orbiting more massive stars, and so our GPIES results for 1 M⊙ come from assuming a continuous distribution jointly constrained by non-detections on these lower-mass stars and six detections around higher-mass stars in the sample.

A number of analyses have also considered the giant planet occurrence rate at these separations for specific categories of target stars. Biller et al. (2013) considered 80 stars that are members of moving groups that were observed by the Gemini NICI Planet-Finding Campaign, and used a Bayesian technique to fit a variant of Equation 2, where the upper semi-major axis cut-off was a free parameter. This analysis included $\beta$ Pic b (which was not detected in the first NICI campaign observation of the star), and AB Pic b, but the small number of detections meant that limits could not be placed on the values of $\alpha$ and $\beta$, though a maximum in probability was reached for a planet fraction of 4%, from 10–150 AU. Wahhaj et al. (2013) performed a similar analysis on a combination of the NICI debris disk subsample and the debris disk hosts in the IDPS high-mass sample (Vigan et al., 2012), with detections of $\beta$ Pic b and HR 8799 bcd. This analysis found 2$\sigma$ limits that placed $\beta$ approximately between -2.5 and -1, and $\alpha$ between -1 and 2.5, with a peak value for the fraction of stars with planets of $\sim$10%. Wahhaj et al. (2013) also fit a model with a stellar mass dependence to the average planet multiplicity ($F_{P}$), which is a separate term from the fraction of stars with planets. These results are suggestive of a rise of planet mulitiplicity with stellar mass, but with a large degree of uncertainty. For their relation $F_{P}\propto M_{*}^{\gamma}$, the 1$\sigma$ results place $\gamma$ between 0.3 and 3.2, and a wider 2$\sigma$ range between -0.5 and 4.7, which encompasses the value of 0, corresponding to no stellar mass dependence.

Bryan et al. (2016) examined the occurrence rate of wide-separation giant planets in systems with a previously-known inner RV planet. By combining long-period RV trends with (generally non-detections from) AO imaging, the possible mass and semi-major axis of the companions could be partially constrained. Their high occurrence rate (62.1${}^{+5.4}_{-5.7}$% for substellar companions between 1-20 MJup between 5-100 au) is noticably higher than what we find with GPIES, suggesting that the presence of a close-in giant planet substantially boosts the likelihood of a wider substellar companion. While they also fit values of $\alpha$ and $\beta$, Bryan et al. (2016) note that these posteriors are strongly dependent on the choice of mass and semi-major axis limits, as those limits truncate the probability distributions for each companion.

A recent paper by Stone et al. (2018) examined the constraints on occurrence rate from a sample of 98 stars observed in the LBT LEECH survey. While a number of planets were observed in the characterization sample, the statistical sample had no planet detections. For all stars, using the COND models Baraffe et al. (2003), they find a 2$\sigma$ upper limit of 25% on the occurrence rate of planets between 4–14 MJup and 5–100 au, consistent with our results.

5.4.1 Brown Dwarf Companions

In addition to detecting six planets, the GPIES campaign has also detected three brown dwarfs, and so we apply our Bayesian framework to inferring the underlying population of brown dwarfs. As before, we fit the four-parameter model of Equation 7 to brown dwarfs in the sample, restricting our fit to semi-major axis values between 1 and 100 au, and companion mass values between 13 and 100 MJup.

The depth of search, as shown by Figure 4, is much more uniform over this mass range, and represents a much deeper sensitivity to these companions: the three brown dwarfs all lie within or near the 160-star contour. Figure 10 displays the results from the brown dwarf analysis from the GPIES 300-star sample.

There is a tail to the posterior at large values of semi-major axis power law, which is truncated by the introduction of a tophat prior, uniform between values of -5 and 5, and 0 elsewhere. This tail results from the limited number of brown dwarfs, and our model being defined in semi-major axis while our detections are made in separation. A very positive power law index would place most brown dwarfs at the largest possible semi-major axis allowed, 100 au in our model. Even if all brown dwarfs have a semi-major axis of 100 au, there is a non-zero probability of detecting three companions at projected separations between 9.9–23.7 au.

As for planets, we again truncate the range for reporting the occurrence rate, 10–100 au and 13–80 MJup. As expected given the deeper sensitivity and fewer detected objects, the brown dwarf occurrence rate is noticeably lower than for giant planets, 0.8${}^{+0.8}_{-0.5}$%.

This occurrence rate is consistent with previous literature measurements of the brown dwarf companion frequency. Metchev & Hillenbrand (2009) conducted a deep survey of 100 stars (out of a larger 266-star sample), and determined an occurrence rate of brown dwarfs (13–75 MJup) orbiting between 28–1590 au of 3.2${}^{+3.1}_{-1.7}$%, assuming $\alpha=0$ and $\beta=-1$. The occurrence rate of a few percent for wide-separation brown dwarfs has been noted by a number of studies, including Oppenheimer et al. (2001), Brandt et al. (2014), Uyama et al. (2017), and Cheetham et al. (2015).

In the bottom left panel of Figure 8, we compare the semi-major axis distributions of GPIES and Metchev & Hillenbrand (2009). Our result is for a lower occurrence rate than that of Metchev & Hillenbrand (2009), with a more negative power-law slope. This decreasing of slope as we move closer to the star is consistent with results from RV surveys for an even smaller brown dwarf occurrence rate at closer separations, with a brown dwarf frequency $<$1% within 5 year periods (Grether & Lineweaver, 2006).

Additionally, we find no significant evidence for a stellar mass dependence of the brown dwarf occurrence rate. We find a value for $\gamma$, the stellar mass power-law, of $-0.9\pm 1.5$. Thus a value of 0, representing no stellar mass dependence, lies within the 1$\sigma$ confidence interval. A similar conclusion was reached by Lannier et al. (2016), who computed companion frequency for stars from the MASSIVE survey of 58 M-stars (which detected two substellar objects) and the survey of 59 mostly-dusty AF stars of Rameau et al. (2013). Focusing on companions with an intermediate mass ratio of 0.01-0.05 (corresponding to 10.5 to 52 MJup for a solar-type primary), Lannier et al. (2016) find no statistical difference between the higher mass and lower mass sample. Similarly, the Gemini-NICI Planet-Finding Campaign surveyed young 70 B and A stars, finding three brown dwarfs (Nielsen et al., 2013), while Bowler et al. (2015) found four brown dwarfs from a survey of 78 single M stars, again suggesting no strong stellar mass dependence in the wide-separation brown dwarf occurrence rate.

5.4.2 A Single Population of Substellar Companions at Wide Separations

Finally, we consider the entire range of substellar objects, and attempt to fit planets and brown dwarfs with the same distribution. GPIES has detected a total of 9 substellar objects, spanning a factor of 25 in companion mass, a factor of 11 in projected separation, and a factor of 1.6 in stellar mass. Thus if a single power-law distribution can describe all these objects, the GPIES 300-star sample provides excellent constraints on the parameters of this distribution.

Figure 11 displays the results from fitting the entire GPIES substellar object yield with Equation 7. Compared to the other subsamples, the power law indices are most constrained with this fit, as expected from the larger number of objects to fit. It is also not surprising that the medians of the posteriors of the three power law indices for the fit to all substellar objects are intermediate to posteriors for the fits to only planets and to only brown dwarfs.

A similar analysis was undertaken by Brandt et al. (2014) based on the SEEDS sample and previously published surveys. Over a similar range to ours, 5–70 MJup and 10–100 au, Brandt et al. (2014) found values of $\alpha=-0.65\pm 0.60$ and $\beta=-0.85\pm 0.39$, compared to the values we derive here of $\alpha=-2.1\pm 0.4$ and $\beta=-1.4\pm 0.5$. Thus we are finding significantly more negative power laws from our GPIES sample, placing more planets at smaller semi-major axis, and smaller masses. Brandt et al. (2014) also found an overall occurrence rate of substellar companions in this regime of $1.7^{+1.1}_{-0.5}$%. While we explicitly fit for $\gamma$, Brandt et al. (2014) do not, but instead adopted a fixed value of $\gamma=1$. The comparison is shown in the bottom right panel of Figure 8, with the two distributions consistent at $\sim$100 au, but the number of planets more steeply rising for GPIES when moving to smaller separations. Vigan et al. (2017) also computed the occurrence rates for substellar companions (0.5–75 MJup, 20–300 au) of 2.1${}^{+1.95}_{-0.60}$%, assuming $\alpha=\beta=0$. With the same assumed distribution, Galicher et al. (2016) derived an occurrence rate of 0.35%, with 95% confidence between 0.10-1.95%, for substellar companions between 5–70 MJup and 10–100 au.

The different values found of $\alpha$, $\beta$, and $f$ by GPIES and Brandt et al. (2014) can be accounted for by differences in the sensitivity and detected companions between the two samples. Brandt et al. (2014) found one planetary mass companion (GJ 504 b), five brown dwarfs, and one intermediate object crossing the planet/brown dwarf boundary ($\kappa$ And b). By contrast, GPIES detected more lower-mass companions, finding six planets and three brown dwarfs. This leads to a more negative companion mass power law, predicting more lower mass companions. In addition, the better IWA of GPI has led to three detections of substellar objects closer than 20 au (51 Eri b, $\beta$ Pic b, and HD 984 B), which in turn drives the semi-major axis power law more negative as well. Finally, the larger number of detections, as well as power laws placing more planets at smaller masses (with lower sensitivity), are responsible for the larger occurrence rate we infer.

5.5 Cold Start

In addition to the ‘hot start’ evolutionary models of giant planets we have considered up until this point, there is an alternate set of ‘cold start’ models. In the cold start models (Marley et al., 2007; Fortney et al., 2008), as the planet forms most of the gravitational potential energy of the infalling gas is radiated away very early on ($\lesssim 5$ Myr) in a shock, so that young planets have comparatively less internal energy to radiate away compared to their hot start counterparts, and so are cooler and less luminous at early ages. At older ages ($\gtrsim 300$ Myr) the cold start and hot start luminosity tracks converge.

The direct observational consequence of the cold start models is fainter, more difficult to detect planets than would be predicted from the hot start models. Additionally, these models predict a pile-up in luminosity at the youngest ages, where planets between 2 and 10 MJup maintain luminosities between about 1-3$\times 10^{-6}$ L⊙ for the first 30 Myr of their evolution (Fortney et al., 2008).

We use cold start evolutionary models (Fortney et al., 2008), and couple them with the AMES-COND atmospheric models (Allard et al., 2001; Baraffe et al., 2003) by matching values of luminosity and temperature predicted by the cold start evolutionary grid with gridpoints in the AMES-COND models to extract $H$ magnitude. We then produce a depth of search plot from the GPIES survey, given in the left panel of Figure 12. There are currently no cold start models for brown dwarf masses, and of the six planets detected by GPIES, five have luminosities too large to have formed according to the cold start models. Only 51 Eri b could potentially be a cold start planet. However, due to the pile-up in luminosity of this model grid, the estimated mass can only be constrained to be between 4 and 10 MJup (Macintosh et al., 2015).

As there is only one possible detection of a cold start planet, rather than attempt to fit a power law distribution, we instead solve for the occurrence rate of cold start giant planets with an assumed power law distribution. Following the method described in Section 5.2, we adopt a uniform distribution in log space, $\alpha=\beta=-1$. The average value over the depth of search between 2–13 MJup, between 3-100 au shown in Figure 12 is then 10.9 stars. From Bayes’ equation with a Poisson likelihood and a Jeffrey’s prior, and one detection over this range, we find an occurrence rate of 15${}^{+15}_{-9}$%. This is intermediate to the value we found for hot start giant planets around high-mass stars when assuming $\alpha=\beta=-1$, 24${}^{+14}_{-10}$% and the 2$\sigma$ upper limit for low-mass stars of $<6.9$% (we caution that 51 Eridani contributed to these occurrence rates for both sets of models). It is not yet certain which set of evolutionary models 51 Eri b follows, and so this cold start occurrence rate could be overestimated if it is indeed a hot start planet. Additionally, these hot start and cold start models represent extremes in evolutionary models, there also exist “warm start” models (Mordasini et al., 2017; Spiegel & Burrows, 2012) that give mass as a function of the initial entropy of the planet (e.g. Figure 17 in Rajan et al. 2017). Direct measurements of the mass of 51 Eri b from Gaia astrometry of the host star could break the degeneracy between the different models.

We also consider the new Sonora 2018 grid (Marley et al. 2019 in prep) with cloud-free atmospheres, and a new cold-start evolutionary grid. This new set of cold start models predicts that young planets are generally hotter by $\sim$200 K and $\sim$5x more luminous than the Fortney et al. (2008) models, as shown in Figure 13. Key differences between the two model sets include the assumed helium abundance (Y=0.24 and 0.28 for Marley et al. 2019 and Fortney et al. 2008, respectively), and the mass of the baryon in the equation of state, with Fortney et al. (2008) using $1\,\rm amu$ and Marley et al. (2019) using the mass of a hydrogen atom.

When using the Sonora 2018 grid, the mass of 51 Eri b is better defined, no longer lying in the overlap region of multiple mass tracks, but at a higher mass of 5.2 MJup. As before, no other detected GPIES planets are low luminosity enough to lie in the range of $H$ magnitudes predicted for 1–10 MJup planets. The depth of search with the Sonora 2018 grid is given in the right panel of Figure 12. The generally brighter planets results in deeper sensitivity to high-mass planets, though still not as deep as for the hot-start BT-Settl models. For the lower-mass giant planets ($\lesssim$3 MJup), however, the Fortney et al. (2008) models predict slightly brighter planets between $\sim$10–50 Myr. Due to a small number of nearby stars with ages $\lesssim$50 Myr, then, the 4-star contour of the depth of search reaches 1.5 MJup for the Fortney et al. (2008) models, compared to 3 MJup for Sonora 2018. Additionally, the Sonora 2018 grid has GPIES able to reach planets at larger semi-major axis, since a high luminosity planet means it can be detected around more distant target stars, so cold start planets can be detected around more of the GPIES sample. Repeating the occurrence rate calculation above gives a slightly lower value, given the greater sensitivity, with 11${}^{+11}_{-6}$% of stars hosting a Sonora 2018 cold start planet between 2–13 MJup and 3–100 au.

Compared to hot start planets, potential cold start planets are very rare among the directly imaged exoplanets. 51 Eri b and GJ 504 b are the only two imaged planets have low enough luminosity to be consistent with cold start models, and so how such a population depends on planet mass, semi-major axis, and stellar mass is not yet well understood. Unlike previous generations of direct imaging exoplanet surveys, GPIES is reaching higher contrasts within 1”, however surveys such as NICI (Nielsen et al., 2013; Wahhaj et al., 2013; Biller et al., 2013; Liu et al., 2013) and IDPS (Vigan et al., 2012; Galicher et al., 2016) reached higher contrasts than GPIES at $\gtrsim 2^{\prime\prime}$, further from the glare of the central star. For young, nearby moving group stars this corresponds to sensitivity to cold start planets to separations of $\gtrsim 50$ au. The relative dearth of detections could indicate a lower occurrence rate than hot start planets, or that cold start planets are found predominantly at a much smaller semi-major axis. If 51 Eri b is indeed a cold start planet, and the occurrence rate of $\sim$10% we infer here is accurate, then future surveys, either with an upgraded GPI (Chilcote et al., 2018) or an upgraded SPHERE, or an ELT high-contrast imager, could allow us to better study these planets.

6 Discussion

At 300 stars, and with detections of six planets and three brown dwarfs, the early GPIES Campaign data is allowing us to probe the distribution of directly-imaged substellar companions at lower masses and closer separations than has been previously possible.

6.1 Do Wide Separation Brown Dwarfs and Giant Planets Follow the Same Distribution?

The 13 MJup boundary that notionally separates giant planets from brown dwarfs represents the onset of deuterium burning, and a longstanding question has been whether there is a significant difference in formation mechanism between these two classes of objects (e.g., Schlaufman 2018). Brown dwarfs and perhaps also giant planets of the highest masses are thought to represent the low-mass outcome of the process by which stars form, i.e., gravitational instability (Kratter & Lodato, 2016; Forgan & Rice, 2013) or turbulent fragmentation (e.g., Hopkins, 2013). Objects of lower mass (Jupiters and below) may have formed from the “bottom up” via core accretion, i.e., by coagulation of solids (e.g., Goldreich et al. 2004; Ormel 2017) and subsequent accretion of gas (e.g., Harris 1978; Mizuno et al. 1978; Mizuno 1980; Pollack et al. 1996). Evidence for this dichotomy is seen for close-in companions observed by RV and transit surveys, where there exists a brown dwarf desert lacking companions near 30 $M_{\rm Jup}$ (e.g., Grether & Lineweaver 2006; Triaud et al. 2017). Such a deep minimum in the occurrence rate vs. companion mass distribution has not been seen at wider separations (Gizis et al., 2001; Metchev & Hillenbrand, 2009). Brandt et al. (2014) found that a single population can fit both brown dwarfs and giant planets at wide separations.

While turbulent fragmentation from a molecular cloud—a process integral to the formation of stars—has been invoked to explain isolated planet-mass companions in systems with primaries not too different in mass (e.g., 2MASSWJ 1207334-393254 and 2MASS J11193254-1137466; Chauvin et al. 2004; Best et al. 2017), such a mechanism is not naturally applied to other directly imaged planets, which show evidence for having formed in a disk. The orbit of $\beta$ Pic b is nearly coplanar with the star’s debris disk (Lagrange et al., 2010; Nielsen et al., 2014; Macintosh et al., 2014); HD 95086 b appears to exhibit this same coplanarity with its star’s disk (Rameau et al., 2016); and the HR 8799 planets are likely coplanar with each other and with the surrounding debris disk (Wang et al., 2018b). Coplanarity points to these planets forming in a disk, either by core accretion or gravitational instability, rather than by turbulent fragmentation from a cloud.

In Section 5.3 we fitted for the parameters to Equation 7 for three cases: using all substellar companions, only giant planets, and only brown dwarfs. In Figure 14 we overplot the posteriors on the four parameters for brown dwarfs alone and planets alone. Since the posteriors overlap, we cannot completely rule out the possibility that the two populations were drawn from the same underlying distribution. Nevertheless, when we follow the same procedure as in Section 5.2, we find the probabilities that the values of $\alpha$ and $\beta$ are larger for brown dwarfs at 87.3% and 74.3%, while the probabilities that $\gamma$ and $f$ are smaller for brown dwarfs are 94.4% and 95.1%, respectively. Thus, we find it likely that planets and brown dwarfs at 10-100 au separation follow different underlying distributions—planets obey a more bottom-heavy mass distribution, and are concentrated toward smaller SMA—at a confidence level between 1 and 2$\sigma$.

The greatest degree of overlap in the posteriors is for $\beta$, while the other three parameters show a more significant difference for planets and brown dwarfs. In particular, for wide separation giant planets this distribution is more likely to have lower companion masses and higher stellar host mass, with a higher overall occurrence rate compared to brown dwarfs.

6.2 Comparing to Demographics of RV-Detected Planets

The RV planet-search method represents one of the more robust ways to measure the underlying distributions of giant planets with periods $\lesssim$ 10 years. While Kepler’s transiting planets greatly outnumber RV-detected planets, the RV method does not suffer as much bias toward longer periods, especially with precision RV surveys operating for multiple decades (e.g., Fischer et al. 2014; Bonfils et al. 2013; Astudillo-Defru et al. 2017). Cumming et al. (2008, C08) presented an analysis of planets orbiting a uniform sample of Keck data of FGK target stars, and fitted a double power law model to mass $m$ and period $P$ (similar in form to our Equation 2), finding that $d^{2}N/(dP\,dm)\propto P^{-0.74}m^{-1.31}$, with an overall planet occurrence rate of 10.5% between 2–2000 days and 0.3–10 MJup. This period distribution is equivalent to a semi-major axis ($a$) distribution of $dN/da\propto a^{-0.61}$ over a range of 0.031–3.1 au for solar-mass stars. Past direct imaging surveys have sought to determine the extent to which this distribution holds for wider-separation planets (Nielsen & Close, 2010; Nielsen et al., 2013; Brandt et al., 2014). While our occurrence rate $f$ gives the number of planets per star, C08 only consider the planet with the largest Doppler amplitude in each system, and therefore calculate the number of planetary systems per star. However, C08 note that this is not a strong bias, as only $\sim$10% of their systems have multiple detected RV planets, and we proceed to use the C08 occurrence rate as an estimate of the number of planets per star.

Figure 7 is inconsistent with C08’s value of $\alpha$ (planet mass index) of -1.31 at the 1$\sigma$ level, and excludes their value of $\beta$ (semi-major axis index) of -0.61 at $\sim$2$\sigma$. If we insist on the C08 values of $\alpha$, $\beta$, and $f$, and extend the distribution out to wider separations, we would predict a planet fraction of 4.7% between 5–13 MJup and 10–100 au—this is also disfavored by GPIES at the 1$\sigma$ level. A more significant difference, though, is in host stellar mass, since the Cumming et al. (2008) sample consisted entirely of planets around FGK stars, while the GPIES planets are exclusively orbiting higher-mass stars.

We further investigate the extent to which the C08 RV distribution is consistent with the GPIES results by setting up a Monte Carlo simulation where planets are drawn from the C08 power laws and occurrence rate, but extrapolated out to an upper semi-major axis cut-off of 100 au, to account for planets seen out to $\sim$70 au (HD 95086 b, HR 8799 b). For now, we assume that the Cumming et al. (2008) distribution fitted to FGK stars applies to stars of all masses, so that planet frequency is not a function of stellar mass ($\gamma=0$), and $\alpha$ and $\beta$ are the same for lower-mass and higher-mass stars. We use the same method that generated the GPIES completeness maps to produce an ensemble of simulated planets drawn from a distribution using the C08 values of $\alpha$, $\beta$, and $f$ around each target star, and determine which planets are detectable given the contrast curve for that star (Nielsen et al., 2008; Nielsen & Close, 2010). We then record the mass, projected separation, and stellar host mass for each detectable planet from the ensemble. Finally, we conduct 10,000 mock GPIES surveys, randomly assigning each star a chance of having a planet in the given mass and semi-major axis range, then using the per-star completeness to determine whether that planet would have been detected.

Figure 15 compares these mock surveys to the results from GPIES. Shown in red are cumulative histograms of the six planets detected by GPIES in planet mass, projected separation, and stellar host mass. Overlaid are these same parameters as measured from the mock surveys. There is generally good agreement in planet mass, as expected from the aforementioned overlap between our measured value of $\alpha$ and that of Cumming et al. (2008). Separation also appears consistent, though the GPIES distribution tends to slightly smaller values of separation than the median of the mock survey, again as we would expect from the GPIES-only fit for $\beta$ finding a more negative value than for the RV planets. The largest departure, however, is clearly for stellar mass, where the median mock survey finds planets around stars over a range of stellar masses (by construction), whereas GPIES observes all planet hosts to be the higher-mass stars in the sample. So while the overall occurrence rate of planets appears similar to predictions of the C08 distribution extended to 100 au, the breakdown by stellar mass is significantly different.

To investigate the dependence of these results on the adopted upper cut-off in semi-major axis, which is currently set to 100 au, and to better account for the difference in stellar masses, we divide our sample into lower-mass stars ($<1.5$ M⊙) and higher-mass stars $\geq$1.5 M*⊙*), and conduct a second series of mock surveys for the two subsamples. This time we keep the value of $\alpha$ and the total occurrence rate within 3.1 au constant using the Cumming et al. (2008) values, and vary the value of the semi-major axis power law index ($\beta$) and the upper cut-off, which sets the maximum semi-major axis planets in the distribution may have. For each combination of these two parameters, we compute the expected number of planets GPIES should have detected using our Monte Carlo simulations. We then use the Poisson distribution to derive the probability that we achieved our actual number of detections in each subsample, displaying the results in Figure 16.

For lower-mass stars ($<1.5$ M⊙), with a null result for planets, we can only set upper limits for the two parameters, by taking the probability of detecting 0 planets given the expectation value. For the Cumming et al. (2008) value of $\beta$=-0.61, we find the upper cut-off must be less than 15 and 34 au at 1 and 2$\sigma$, respectively. Thus it is plausible that for solar-mass stars the Cumming et al. (2008) distribution extends into the giant planet region of our own solar system, as our GPIES sensitivity does not extend inside 10 au.

For higher-mass stars ($>1.5$ M⊙), since we have detected six planets, we can better constrain the values of $\beta$ and the upper cut-off within this modified Cumming et al. (2008) distribution. We again use the Poisson distribution to compute the probability of detecting six planets around the GPIES higher-mass stars at each point in the grid, then draw contours encircling 68%, 95%, and 99.7% of the total probability. We set a prior that imposes a lower limit on the semi-major axis cut-off of 70 au, corresponding to HR 8799 b (Wang et al., 2018a). Though we only detected HR 8799 cde given the field of view of GPI, we assert that a distribution fitted to the inner three planets of the system should not exclude the outer planet b. The Cumming et al. (2008) value of $\beta$=-0.61 intersects the 1$\sigma$ contour at 264 au, and intersects the 2$\sigma$ contour at a semi-major axis cut-off of 85 au. Such a large value for the upper cut-off is difficult to reconcile with the relative lack of detections of planets around higher-mass stars beyond $\sim$100 au by previous generations of imaging surveys. An extrapolation of the C08 distribution would predict an equal number of planets between 10-80 au (where all six GPIES planets reside) and between 80–264 au. The stated errors of C08 in the period distribution ($\frac{dN}{dlnP}\propto P^{0.26\pm 0.1}$), if they are Gaussian, would correspond to errors in $\beta$ for semi-major axis of $-0.61\pm 0.15$, and so a 1$\sigma$ value of $-0.46$ would be more consistent with our results, reaching the 1$\sigma$ value in the right panel of Figure 16 at 79 au. Adopting such a value for $\beta$ would make the upper limits for lower-mass stars more stringent, moving from 15 and 34 au at 1 and 2$\sigma$ to 12 and 23 au, respectively. Previously we had fitted for all three parameters in Equation 2 simultaneously—$\alpha$, $\beta$, $f$—while fixing the value of the upper cut-off to 100 au. Here, as we fix $\alpha$ and $f$ to the Cumming et al. (2008) values and let $\beta$ and the upper cut-off float, we reach a similar conclusion, namely that GPIES excludes C08’s value of $\beta$ at $\sim$1$\sigma$, and that adjusting the value of the upper cut-off is not a straightforward path to removing this discrepancy.

The upper cut-off analysis above considers only the occurrence rate when evaluating the consistency of the extended Cumming et al. (2008) distribution with our GPIES data. In order to consider not just occurrence rate but also the properties of the detected planets, we repeat the mock surveys on only the 123 stars above $1.5M_{\odot}$, assuming an upper semi-major axis cut-off of 100 au, and the C08 values of $\alpha$ and $\beta$. Given this modified distribution and our completeness to planets around our 123 higher-mass stars, the mean number of planets found is only 2.6. To more closely match our actual observed number of 6 planets, we increase the $<3.1$ au planet occurrence rate from the Cumming et al. (2008) value by a factor of 2 (which for a 1.5 M⊙ star would represent a value of $\gamma$ near our median value from Figure 9 of 2.03). This moves the mean number of detections up to 5.2, with 42% of the mock surveys now producing 6 or more detected planets. As before, we assume there is no stellar mass dependence in planet occurrence within this subsample of higher-mass stars ($\gamma=0$).

Results from this third set of mock surveys restricted to higher-mass stars ($>1.5$ M⊙) are plotted in Figure 17. While the planet mass distribution and overall occurrence rate appear consistent with the observed GPIES detections, there remain significant differences for the projected separation and stellar host mass distributions. Again, the Cumming et al. (2008) separation distributions place more planets at wider separation than are observed. We saw the same disagreement with the Cumming et al. (2008) distribution in our full power-law fit (Figure 7), where we found a more negative value of $\beta$ from the GPIES sample. Indeed, 80.2% of mock surveys have a detected planet at a projected separation beyond 53.65 au, the value for HD 95086 b. We could try to reduce this discrepancy by boosting $f$ still further and decreasing the upper cut-off in the semi-major axis distribution. But lowering the cut-off would exclude, without physical justification, observed planets like HD 95086 b in our sample at 55 au, as well as the widest-separation planet in the HR 8799 system, planet b at 70 au. Though HR 8799 b does not appear in our sample, the packed nature of the system (Wang et al., 2018b) leads us to reject a model fit to the inner three planets that does not fit the fourth planet.

Figure 17 also shows that despite our restriction to $>1.5$ M⊙ stars, a mismatch persists in the stellar mass distributions: the mock surveys predict planets to be detected around more 2 $M_{\odot}$ stars than are actually observed. For such high-mass stars in the GPIES sample, there is no obvious bias where more massive stars are systematically older or more distant. The persistent discrepancy in stellar mass distributions likely arises because the C08 distribution places about as many planets at large orbital distances as at small ones, and those at large distances are detectable in the mock surveys of higher-mass stars. In the actual GPI survey, as the absolute $H$ magnitude decreases by $\sim$1 mag from 1.5 to 2 M⊙, planets like 51 Eri b, $\beta$ Pic b, and HR 8799 e that cluster at small orbital separations become undetectable around higher-mass stars.

In sum, we find that the semi-major axis distribution from C08 of RV-discovered giant planets is a poor fit to directly imaged giant planets at wider separations. For lower-mass stars ($<1.5$ M⊙), while we cannot exclude the possibility that the Cumming et al. (2008) RV population extends to orbital distances of 15–34 au, we do not have any data to support this model, either. For higher-mass stars, the data from GPIES place more planets at smaller semi-major axes than at larger ones; GPIES prefers a semi-major axis power-law index of $\beta\simeq-2$, which is discrepant from the C08 value of -0.61 at 2$\sigma$.

Fernandes et al. (2018) analyze RV and Kepler giant planet occurrence rates, and observe a turnover in occurrence rate as a function of semi-major axis, with a peak at $\sim$3 au. Their symmetric epos fit gives a value of the planet mass power law $\alpha=-1.45\pm 0.05$ (similar to the C08 value of $\alpha=-1.31\pm 0.2$), and a broken power law fit to semi-major axis, with $\beta=-0.025^{+0.30}_{-0.22}$ for small periods and $\beta=-1.975^{+0.22}_{-0.30}$ for large periods, again converting a period power law index into semi-major axis. At small semi-major axis, then, Fernandes et al. (2018) find a more positive value of power law index compared to the C08 value of $\beta=-0.61\pm 0.15$, and a value of $\beta$ similar to our fit to planets around higher-mass stars ($\beta=-2.0^{+0.6}_{-0.5}$) and planets around stars with a stellar mass term ($\beta=-1.7^{+0.6}_{-0.5}$). For the mass power law index, we find a more negative value of $\alpha=-2.3^{+0.8}_{-0.6}$, and while the 1$\sigma$ confidence intervals overlap, the offset could be a sign of a different planet mass distribution at larger separations, or between higher-mass and lower-mass stars. Fernandes et al. (2018) find the location of the period break from their fit at $P=1581^{+894}_{-392}$ days, which corresponds to a semi-major axis for solar-type stars of $2.66^{+0.93}_{-0.46}$ au. The asymmetric epos fit gives similar values for $\alpha$, $\beta$ within the break, and overall occurrence rate, as well as consistent values for $\beta$ beyond the break and the location of the break, though with significantly larger error bars: $\beta=-2.8^{+1.4}_{-1.9}$ and a break at $3.2^{+1.1}_{-1.4}$ au. We find support for the conclusions of Fernandes et al. (2018), that the power law index $\beta$ measured by C08 cannot continue to large values of semi-major axis, but must turn over within $\lesssim$10 au (modulo the difference in stellar mass between the planet hosts in the C08 and Fernandes et al. 2018 samples and GPIES).

Finally, we combine RV results for small-separation planets and GPIES results at larger separations to estimate the total number of planets within 100 au, as shown in Figure 18. In the top two panels, we couple GPIES and the C08 occurrence rate to calculate the number of planets as a function of semi-major axis between 1–13 MJup around 1 M⊙ stars. Draws are taken from the MCMC chain of our fit to giant planets including a stellar mass term, and combined with Monte Carlo draws of $\alpha$, $\beta$, and occurrence rate from C08, assuming each is Gaussian distributed and that the three parameters are independent, thus producing two power law distributions. We attempt to make the combined distribution continuous by finding a point where the two power laws cross between 1–20 au; if the two do not intersect in this range, the handover point is taken to be 10 au, with a discontinuity in the final combined distribution. Finally, we integrate the distribution to find the posterior on the total number of giant planets. We find a total occurrence rate (number of planets per star) of 0.09$\pm$0.02 from $0.031-10$ au, and 0.04${}^{+0.09}_{-0.03}$ from $10-100$ au, producing a total occurrence rate within 100 au of 0.14${}^{+0.10}_{-0.05}$. Repeating this analysis with the Fernandes et al. (2018) epos symmetric distribution produces similar results, with occurrence rates of 0.07${}^{+0.03}_{-0.02}$, 0.04${}^{+0.09}_{-0.03}$, and 0.12${}^{+0.10}_{-0.05}$ respectively.

As GPIES only detected planets around higher-mass stars, the posterior for giant planet occurrence approaches zero between 10–100 au, though the GPIES data do set a clear upper limit given the lack of detections of giant planets around these lower-mass stars. By assuming this power-law model, we find a 2$\sigma$ upper limit on the occurrence rate of 1–13 MJup planets around 1 M⊙ stars of 39.0% and 38.1% for the C08 and Fernandes et al. (2018) distributions, respectively.

For higher-mass stars, we repeat this analysis, this time using the Johnson et al. (2010a) occurrence rate, stellar mass power law ($\gamma$) and errors, and setting the mass power law $\alpha$ equal to the C08 value and uncertainties, and allowing $\beta$ to follow a uniform distribution between -1 and 0. We then draw from our MCMC fits to the high-mass star GPIES sample and continue as before. Extrapolating the Johnson et al. (2010a) results gives an occurrence rate within 10 au for 1.75 M⊙ stars of 0.42${}^{+0.20}_{-0.12}$. Given the large occurrence rate determined by GPIES for 5–13 MJup, and the very negative value of $\gamma$, the total occurrence rate for 1–13 MJup planets between 10–100 au orbiting higher-mass stars is even larger, 1.0${}^{+1.7}_{-0.6}$, with a total occurrence rate within 100 au of 1.5${}^{+1.7}_{-0.7}$ for these higher-mass stars. Given the negative value of the planet mass power law index ($\alpha=-2.4^{+0.9}_{-0.8}$) found from GPIES, many of these predicted planets have small masses. With GPIES having limited sensitivity to $\lesssim 3$ MJup planets orbiting higher-mass stars, we cannot rule out a break in the power law distribution at lower planet masses, which would result in a smaller occurrence rate.

6.3 Connecting to Planet Formation Theory

Our finding that giant planets outnumber brown dwarfs at more than 10-to-1 around $>1.5M_{\odot}$ stars argues against these planets having formed by gravitational instability. Fragmentation of marginally Toomre-unstable disks (Toomre, 1964; Goldreich & Lynden-Bell, 1965) is thought to produce a mass distribution weighted toward larger and not smaller masses. When fragments initially condense, they may already be at least as massive as Jupiter (e.g., Figure 3 of Kratter & Lodato, 2016; see also Rafikov 2005). In a population synthesis study predicated on gravitational instability, Forgan & Rice (2013) found the overwhelming majority of clumps that survived tidal disruption to lie above the deuterium burning limit. These initial fragments are expected to accrete still more mass from their strongly self-gravitating, dynamically active parent disks, which themselves may still be accreting from infalling progenitor clouds (Kratter et al., 2010).

In addition, the requirement that disks rapidly cool in order to fragment (Gammie, 2001) predicts that collapse occurs at large stellocentric distances where material is less optically thick, e.g., outside at least $\sim$50 au in a model inspired by the HR 8799 system (Figure 3 of Kratter et al., 2010; see also Matzner & Levin 2005 and Rafikov 2009). The generic predictions of gravitational instability—a top-heavy mass distribution extending to brown dwarfs and stars, and large orbital separations (modulo orbital migration)—are not borne out by our 2–13 $M_{\rm Jup}$ planets located between 10–60 au. Instead we infer mass distributions (using the fit to planets around higher-mass stars) that are bottom-heavy ($dN/dm\propto m^{\alpha}$ with $\alpha=-2.4\pm 0.8$) and separation distributions that favor small semi-major axes ($dN/da\propto a^{\beta}$ with $\beta=-2.0\pm 0.5$).

Formation by gravitational instability remains an option for GPIES brown dwarfs, for which we infer mass and separation indices of $\alpha=-0.5\pm 1.1$ (more top-heavy) and $\beta=-0.7^{+1.0}_{-1.5}$ (weighted more toward large distances). The three brown dwarfs in our survey are located at orbital distances of 10–30 au, somewhat smaller than the values cited above for fragmentation. However, previous surveys with wider fields of view have found brown dwarfs at larger separations, e.g., HR 7329 B at 200 au (Lowrance et al., 2000a), $\zeta$ Del B at 910 au (De Rosa et al., 2014), and HIP 79797 Ba/Bb, a pair of brown dwarfs in a 3 au orbit separated by 370 au from the primary star (Huélamo et al., 2010; Nielsen et al., 2013). Moe & Kratter (2018) present mechanisms (e.g., Lidov-Kozai cycles) for stellar/brown dwarf binary systems to decrease their separation over time.

Another observed property of brown dwarf companions that we have highlighted is their lack of preference for either higher-mass or lower-mass stars as host. This, too, is consistent with formation within gravitationally unstable disks. If the disk gas surface density $\Sigma_{\rm g}\propto M_{*}$ and the sound speed $c_{s}\propto T^{1/2}\propto M_{*}^{1/4}$ (where we have taken the gas temperature $T\propto L_{*}^{1/4}$ and used the pre-main-sequence luminosity-mass relation $L_{*}\propto M_{*}^{2}$), then Toomre’s $Q\propto 1/M_{*}^{1/4}$. This weak dependence on stellar mass is made even weaker when we fold in the large scatter in observed gas disk masses for given $M_{*}$ (e.g., Andrews et al., 2018). Deeper analyses that combine the $Q$ criterion with the requirement of fast cooling reveal that disks gravitationally fragment at orbital locations that are remarkably insensitive to system parameters related to stellar mass (Matzner & Levin, 2005; Rafikov, 2009).

Our finding that giant planets at 10–100 au are more commonly hosted by A stars than by lower-mass stars extends the trend determined by prior demographic studies that closer-in gas giants are more common around FGK stars than around M dwarfs (Clanton & Gaudi, 2014; Mulders et al., 2015; Ghezzi et al., 2018). The correlation with stellar mass presents a major constraint on theory. Having argued above against planet formation by gravitational instability (a ‘top-down’ scenario), we ask here whether the trend with stellar mass is consistent with core accretion, a ‘bottom-up’ scenario whereby giants nucleate from sufficiently massive rocky cores. Note that the high luminosities of directly imaged planets, though consistent with hot start models and inconsistent with cold start models (Section 5.5), do not necessarily rule out core accretion. Cold and hot start models differ according to the planet’s assumed initial entropy, which can vary significantly within core accretion depending on the shock dynamics of runaway gas accretion (see, e.g., Figure 12 of Berardo et al., 2017).

In a prescient analysis, Laughlin et al. (2004) anticipated fewer giant planets around lower-mass stars, arguing that core formation times would be longer around M dwarfs, and potentially too long compared to gas disk dispersal times, because their disks have longer orbital periods and are expected to have fewer solids. We update this reasoning within the emerging paradigm of ‘pebble accretion,’ in which cores grow by accreting small solid particles aerodynamically entrained in disk gas (see Ormel, 2017 and Johansen & Lambrechts, 2017 for reviews). As particles drift radially inward through the gas disk (e.g., Weidenschilling, 1977), a fraction of them are accreted by the protocore. The final core mass is a monotonic increasing function of $\int^{t_{\rm drift}}\Sigma_{\rm s}dt$: the disk solid surface density $\Sigma_{\rm s}$ (local to the protocore) integrated over the time it takes solids to drift past the protocore’s orbit and drain out from the disk (see equations 36, A7, and A9 of Lin et al., 2018). There are dependencies of the final core mass on variables apart from this integral, but they are unlikely to change the qualitative conclusion that pebble accretion breeds more massive cores around more massive stars.222In equations (A7) and (A9) of Lin et al. (2018) pertaining to the ‘2D’ limit where the pebble accretion radius exceeds the solid disk scale height, the extra variable dependencies only strengthen the trend with stellar mass. In the case of their equation 36 (‘3D’), the extra variables mute the trend but do not reverse it, unless the particle Stokes number $St$ (aerodynamic stopping time normalized to the orbit time) decreases more rapidly than $1/M_{*}^{3/2}$. Pascucci et al. (2016) support $St\propto 1/M_{*}^{1/2}$ for fragmentation-limited particle growth (Birnstiel et al., 2012); this relation steepens to $1/M_{*}$ if particle sizes are fixed and gas disk masses scale linearly with stellar mass. The surface density $\Sigma_{\rm s}$ scales with the total solid disk mass, which should increase with $M_{*}$ (a linear dependence is commonly assumed). The drift time $t_{\rm drift}$ also increases with $M_{*}$, as inferred empirically from millimeter-wave dust continuum observations showing that disks are systematically larger around higher-mass stars (Pascucci et al., 2016, in particular their section 6.2; see also Andrews et al., 2018). Thus pebble/core accretion gets at least the sign right with respect to the correlation of giant planet occurrence with host star mass, and it remains to work out more quantitative models that can satisfy the constraints reported in our Table 3.333Which is not to say that there are not outstanding problems intrinsic to the theory. Lin et al. (2018) point out that pebble drift times in Jupiter-breeding scenarios do not match observations. There may also be insufficient mass in pebbles (Manara et al. 2018). Size distributions of accreting particles are crucial but uncertain; the origin of seed cores is tied to the perennial problem of planetesimal formation. See Johansen & Lambrechts (2017) for some current thinking.

What about semi-major axis trends in the context of pebble accretion? The preference reported here for giant planets to be located at the lower end of the 10-100 au range makes sense insofar as protocores at smaller orbital distances access larger reservoirs of solids lying exterior to their orbits and drifting past. Pebble accretion probabilities also tend to be higher closer in because disk gas is denser there and more efficiently drags pebbles onto cores: accretion at small distances is more likely to be in the so-called ‘settling’ regime, with particle trajectories more strongly gravitationally focussed, as opposed to the less efficient ‘hyperbolic’ regime which obtains farther out (Ormel & Klahr, 2010; Lin et al., 2018, in particular their Figure 15 which contrasts fast settling accretion at 5 au with slow hyperbolic accretion at 40 au).

We can also try connecting the observed preference of directly imaged giant planets for smaller semi-major axes—which we reiterate applies to distances of 10–100 au—to the observed preference of giant planets discovered by radial velocities (RV) for larger semi-major axes, which applies to distances of 0.03–3 au (e.g., Cumming et al., 2008; Fernandes et al., 2018). Transit surveys report a similar rise in giant planets out to $\sim$300 days (e.g., Dong & Zhu, 2013; Petigura et al., 2018). Assuming these semi-major axis trends can be stitched together—keeping in mind the caveat that they derive from disjoint samples of higher-mass and lower-mass stars—we infer the giant planet occurrence rate peaks between $\sim$3–10 au. Remarkably, a similar conclusion has been found by Fernandes et al. (2018) in a combined RV-transit analysis of FGK stars: they find evidence for a turn-over in the giant planet occurrence rate at $\sim$2–3 AU. A giant planet occurrence rate that peaks at intermediate distances aligns at least qualitatively with pebble/core accretion theory. The aforementioned gains in growth with smaller orbital distance are eventually reversed, at the shortest distances, by other competing physical effects. Cores stop accreting pebbles when they gravitationally perturb the surrounding disk gas so strongly that the inward drift of pebbles stalls (Lambrechts et al., 2014). Such ‘pebble isolation’ sets an upper mass limit that decreases as the orbital distance decreases (Bitsch et al., 2018; Fung & Lee, 2018); cores do not grow past sub-Earth masses inside $\sim$0.1 AU (see Figure 12 of Lin et al., 2018, who rule out in-situ formation of hot Jupiters by pebble accretion on this basis). Another effect stifling growth at the shortest distances is that disk gas becomes hotter and more opaque there and therefore cools and accretes onto cores more slowly (Lee & Chiang, 2015, 2016). All the above theoretical reasons, together with the many arguments pointing to inward migration of hot Jupiters by Kozai cycles/tidal friction and perhaps also disk torques (Dawson & Johnson, 2018), lead us to posit a “sweet spot” for gas giant formation at stellocentric distances that are neither too large nor too small. From Figures 8 and 18, this sweet spot appears to be at $\sim$2–10 au.

7 Conclusions

From the first 300 stars observed out of the planned 600-star survey, reaching contrasts of 106 within 1*′′* radius, GPIES is one of the largest and deepest direct imaging surveys for exoplanets conducted to date. Our analysis of the data shows that there is a clear stellar mass dependence on planet occurrence rate, with stars $>1.5M_{\odot}$ more likely to host giant planets (5–13 $M_{\rm Jup}$) at wide separations (semi-major axes 10–100 au) than lower-mass stars. Around higher-mass stars, the total occurrence rate of such planets is 9${}^{+5}_{-4}$%. When fitting our data for higher-mass stars ($>1.5$ M⊙) with a power-law model of planet populations where $d^{2}N/(dm\,da)\propto m^{\alpha}a^{\beta}$, we find negative power law indices of $\alpha\simeq-2$ and $\beta\simeq-2$, values that favor smaller masses and smaller semi-major axes. All the statistical trends that we have uncovered—a bottom-heavy companion mass distribution, orbital distances concentrated more toward 10 au than 100 au, and a strong preference for higher-mass host stars—appear consistent with a “bottom-up” formation scenario: core/pebble accretion.

By comparison to giant planets, brown dwarfs appear to exhibit different statistics: a more top-heavy mass distribution ($\alpha\simeq-0.47$), orbital distances that are somewhat more concentrated toward 100 au than 10 au ($\beta\simeq-0.65$), and no particular preference for host stellar mass. These are consistent with a “top-down” formation scenario: gravitational instability and fragmentation of circumstellar disks. Our finding of a brown dwarf occurrence rate on the order of 1%, independent of stellar mass, is consistent with previous literature analyses of brown dwarfs from direct imaging surveys.

While planets around GPIES stars are more common around higher-mass stars ($>$1.5 M⊙), it is interesting to note that the most massive planet host in the sample is $\beta$ Pictoris (1.73 M⊙), yet there are 76 stars more massive than 1.8 M⊙ in our 300-star sample. Repeating the analysis of Section 5.2, but dividing the higher-mass sample between 1.5–1.8 M⊙ and $>$1.8 M⊙, we find at 99.7% confidence that the intermediate stellar mass bin has more planets than the highest stellar mass bin. But this result is dependent on our choice of bin boundary—when we extend the dividing line from 1.8 to 2.1 to 2.4 M⊙, our confidence decreases from 99.7% to 97.1% to 87.1%, respectively. It is thus difficult to know if we are observing a real effect or whether we are fine-tuning our boundary based on a small number of detections. For comparison, the limit we used here to separate higher-mass and lower-mass stars, 1.5 M⊙, must be reduced to 1.04 M⊙ (with only 67 stars out of 300 below this limit) before the confidence with which we can conclude that higher-mass stars have more planetary systems than lower-mass stars drops below 95%. If the frequency of wide-separation giant planet systems does indeed drop off beyond $\sim$2 M⊙, this would suggest that the mechanism that formed the giant planets in the GPIES sample becomes less efficient at stellar masses above this value. Jones et al. (2016) observe reminiscent behavior for close-in giant planets around giant stars, with a planet occurrence rate that increases out to $\sim$2.3 M⊙, but then decreases at higher stellar masses. Jones et al. (2016) model this as a Gaussian, though like us they have fewer targets in their sample at the high stellar mass end. Additionally, their RV targets have evolved off the main sequence, while ours are predominantly on the main sequence, and so it is unclear how to disentangle formation effects from post-main sequence evolution effects. Nevertheless, further observations of additional higher mass stars ($\gtrsim 2$ M⊙) can test whether this trend is real.

There are a limited number of bright, young, nearby stars available to high-contrast direct imagers like GPI and SPHERE. As a result, at the current achievable contrasts there are unlikely to be a significant number of new detections of imaged giant planets following the completion of the GPIES and SPHERE/SHINE (Feldt et al., 2017) campaigns. An improvement in performance from instrument upgrades (Chilcote et al., 2018), however, could unlock mass/separation phase space around these target stars that is currently inaccessible. For example, 51 Eridani was observed multiple times by planet-hunting surveys (e.g., Heinze et al. 2010, Rameau et al. 2013, Biller et al. 2013), yet the planet 51 Eri b was not discovered until it was imaged with the more advanced GPI (Macintosh et al., 2015). Additional planets discovered with upgraded GPI and SPHERE would allow us to test the robustness of the trends with stellar mass and planetary mass uncovered in this paper. The Gaia mission will identify astrometric signatures of planets with $\gtrsim 10$ year orbits that may be confirmed with direct imaging, which can also boost the number of companions in this separation range.

Alternatively, high contrast imaging at longer wavelengths offers the opportunity to expand the list of potential imaging targets to older stars ($\sim$300 Myr–1 Gyr). While a 5 MJup planet is 100 times brighter at 30 Myr than 600 Myr at $H$-band, this ratio is closer to a factor of 10 at $M^{\prime}$ (BT-Settl, Baraffe et al. 2015). For a uniform star formation rate, there should be nine times more stars between 100 Myr and 1 Gyr than there are stars younger than 100 Myr. JWST, GMT, TMT, and ELT will be well-suited to probe wide-separation giant planet populations around these intermediate-age stars at 5–10 $\mu$m. Going forward, combining different detection methods to produce a full accounting of giant planet and brown dwarf populations from 0.01 to 1000 au will provide still more constraints on formation mechanisms.

We thank the referee for helpful comments that improved the quality of this manuscript. Based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), and Ministério da Ciência, Tecnologia e Inovação (Brazil). This research has made use of the SIMBAD and VizieR databases, operated at CDS, Strasbourg, France. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. J.R., R.D. and D.L. acknowledge support from the Fonds de Recherche du Québec. JRM’s work was performed in part under contract with the California Institute of Technology (Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. Support for MMB’s work was provided by NASA through Hubble Fellowship grant #51378.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. Supported by NSF grants AST-1411868 (E.L.N., K.B.F., B.M., and J.P.), AST-141378 (G.D.), and AST-1518332 (R.D.R., J.J.W., T.M.E., J.R.G., P.G.K.). Supported by NASA grants NNX14AJ80G (E.L.N., S.C.B., B.M., F.M., and M.P.), NNX15AC89G and NNX15AD95G (B.M., J.E.W., T.M.E., R.J.D.R., G.D., J.R.G., P.G.K.), NN15AB52l (D.S.), and NNX16AD44G (K.M.M). K.W.D. is supported by an NRAO Student Observing Support Award SOSPA3-007. J.J.W. is supported by the Heising-Simons Foundation 51 Pegasi b postdoctoral fellowship. E.C. thanks Courtney Dressing for helpful discussions. Portions of this work were performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. This work benefited from NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. Gemini:South (GPI)

Appendix A Properties of Target Stars