Population-level Eccentricity Distributions of Imaged Exoplanets and Brown Dwarf Companions: Dynamical Evidence for Distinct Formation Channels^*

We targeted HD 49197 using the High Contrast Instrument for the Subaru Next Generation Adaptive Optics (HiCIAO; Hodapp et al. 2008; Suzuki et al. 2010) near-infrared imager coupled with the AO188 AO system (Hayano et al. 2010) at Subaru Telescope on UT 2011 December 28 (see Table 1 for details). Conditions were photometric, but the seeing was poor and variable; UKIRT reported K-band natural seeing measurements between 15 and 2'' during our observations. We acquired a total of 60 frames in K_S band, each with an integration time of 30 s. Observations were taken using natural guide star AO in pupil-tracking mode (angular differential imaging; Liu 2004; Marois et al. 2006), which uses field rotation to distinguish speckles from real point sources. HD 49197 was placed behind the 300 mas diameter opaque Lyot coronagraph during the ADI sequence, which spanned 19° of sky rotation.

Data reduction and point-spread function (PSF) subtraction follow the steps detailed in Bowler et al. (2015a). Systematic bias stripes from the detector readout electronics are measured and subtracted, cosmic rays and bad pixels are removed, then images are divided by a normalized flat field. The K_S-band distortion solution from Bowler et al. (2015a) is applied to each image to correct for optical aberrations. The corresponding K_S-band plate scale of 9.67 ± 0.03 mas pixel⁻¹ is adopted for this data set. The typical residual rms on the distortion correction is 1.2 pix, or 11.6 mas. Celestial North was found to be aligned with the detector columns to within the measurement errors, so no rotation was applied and a value of 00 ± 01 is adopted. Images are registered using a 2D elliptical Gaussian fit to the PSF wings surrounding the coronagraph and then assembled into a data cube. Finally, PSF subtraction is carried out using "aggressive" and "conservative" implementations of PSF subtraction with the locally optimized combination of images (LOCI; Lafrenière et al. 2007). The aggressive subtraction using LOCI parameters W = 8, N_A = 300, g = 1, N_δ = 0.5, and dr = 2 produced a higher signal-to-noise ratio (S/N) for the modest-contrast (ΔK_S ≈ 8 mag) companion HD 49197 B, so we adopt this version of the reduction here. The final processed image is shown in Figure 1.

We observed 13 targets with substellar companions between 2014 and 2019 (Table 1) with the NIRC2 camera behind natural guide star AO at Keck Observatory (Wizinowich et al. 2000). All observations were acquired with the narrow camera mode, which provides a plate scale of ≈10 mas pix⁻¹ and a field of view of 102 × 102. Observations of HD 19467, κ And, GJ 504, and HD 49197 were acquired in pupil-tracking mode to facilitate standard post-processing PSF subtraction. Total on-source integration times ranged from 3 to 81 minutes and the field-of-view rotation angle ranged from 3° to 119°. The brown dwarf companions to HD 1160, 1RXS0342+1216, CD 35–2722, DH Tau, HD 23514, Ross 458, TWA 5, 2M1559+4403, and 1RXS2351+3127 have lower contrasts and were observed with shorter integration times. The partly transparent 600 mas diameter focal plane coronagraph was used for all systems except Ross 458. Details of the observations can be found in Table 1.

After removing bad pixels and cosmic rays, images are bias subtracted, flat fielded, and corrected for optical distortions using the distortion solution from Yelda et al. (2010) for observations taken before 2015 April and from Service et al. (2016) for observations after that date. The corresponding plate scales are 9.952 ± 0.002 mas pixel⁻¹ and 9.971 ± 0.004 mas pixel⁻¹ for these respective dates, and the P.A. offsets are 0252 ± 0009 and 0262 ± 0020 with respect to the detector columns. The typical rms residual after distortion correction is ≈1 mas. PSF subtraction for the ADI sequences is carried out following the description in Bowler et al. (2015a): images are first registered using the position of the star visible behind the focal plane mask, then PSF subtraction is performed using LOCI. The median-combined PSF-subtracted image is then north aligned and a noise map is created by measuring the rms in annuli centered on the host star with a width of 3 pixels. The final PSF-subtracted images and S/N maps are shown in Figures 1–5. Note that HD 1160 was observed in field-tracking mode; PSF subtraction for this target entailed derotating to a common pupil angle, implementing PSF subtraction with LOCI, then rerotating to a common sky position angle before coadding the individual frames.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

No PSF subtraction is necessary to recover the other substellar companions with modest flux ratios, which are readily visible in the raw frames. After basic image reduction, these frames are registered, north aligned, and then coadded to produce the final images shown in Figure 6.

Zoom In Zoom Out Reset image size

Relative astrometry is measured in each final processed image. Separations in pixels are converted into angular distances using the detector plate scale. Following Bowler et al. (2018), uncertainties take into account random positional measurement errors, estimated to be 0.5 pix for HiCIAO and 0.1 pix for NIRC2 based on end-to-end injection-recovery tests of companions and host star positions in Bowler et al. (2015a) and Bowler et al. (2018); uncertainties in the plate scale; and rms errors from the distortion correction. Similarly, the P.A. uncertainty incorporates random measurement errors, uncertainty in the absolute orientation of celestial north on the detector, rms errors from the distortion correction, and angular uncertainty associated with average azimuthal shearing of the PSF within each image for observations taken in pupil-tracking mode. Bowler et al. (2018) found that systematic errors can dominate the astrometric uncertainty budget for ADI observations with NIRC2 using the 600 mas coronagraph. We therefore adopt conservative errors of 5 mas and 03 for our ADI observations in this work.⁷ Our final measurements are reported in Table 1.

The S/N of each companion detection is calculated using aperture photometry. Aperture radii of ≈λ/D are chosen to encapsulate the central Airy disk and minimize potential effects of oversubtraction at larger radii. Here D is the telescope diameter and λ is the central wavelength of the filter. We adopt a 6 pix aperture for observations using HiCIAO in K_S band, a 5 pix aperture for NIRC2 in K_S band, and a 4 pix aperture for NIRC2 in H band. Noise levels are derived using counts in 100 circular apertures at the same separation but different position angles in each image. These S/N measurements range from 9 (for κ And B) to 816 for our 2014 epoch of HD 49197 (see Table 1).

Figures 7 and 8 compare our new astrometry (Table 1) to published values in the literature (Appendix A). Orbital motion is clearly detected in many systems and generally extends and refines linear evolution in separation and P.A. with time, with the exception of Ross 458, which has a well-determined orbit. The addition of recent observations is especially important for systems like CD–35 2722 B, for which no astrometry has been published since its discovery by Wahhaj et al. (2011).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

There are several instances in which our observations were taken close in time with other published epochs, which provides an opportunity to compare both measurements for mutual consistency. Our astrometry of HD 49197 from UT 2016 March 22 was taken four months after the observations by Bottom et al. (2017) on 2015 November 22 UT with the Stellar Double Coronagraph at Palomar Observatory. Our measured separation of 874 ± 5 mas and P.A. of 765 ± 03 agrees well with the values of 862 ± 25 mas and 766 ± 18 from Bottom et al. We targeted GJ 504 about one week before an observation taken with SPHERE that was reported in Bonnefoy et al. (2018). They find {ρ = 2495 ± 2 mas, θ = 32248 ± 005} on UT 2016 March 29; our measurements of {ρ = 2504 ± 5 mas, θ = 3227 ± 04} were obtained on UT 2016 March 22. These P.A.s are consistent at the <1σ level, and our separation measurement is consistent within 2σ (9.0 ± 5.4 mas). Our UT 2018 January 30 astrometry of HD 1160 B ({ρ = 790.0 ± 5 mas, θ = 2451 ± 03}) is within 1σ of the values also taken with NIRC2 by Currie et al. (2018) on UT 2017 December 9 ({ρ = 784 ± 6 mas, θ = 2449 ± 03}), about two months before our observations. Despite this good agreement with previous observations, below we attempt to address any potential systematic errors that could result from combining measurements taken with different instruments. This is especially important because systematic errors or underestimated uncertainties can mimic a strong acceleration, which is most readily accounted for in orbit fits through a high eccentricity and time of periastron close to the epoch of observations (see, e.g., results for HR 8799 e from Konopacky et al. 2016).

Linear evolution in astrometry is generally expected for most companions because their orbital periods are long relative to the time baseline of the observations (≈5–20 yr). We adopt reduced χ² values as a metric to assess the fidelity of the astrometric uncertainties, both for our own measurements and for those from the literature. These sample a wide range of instruments, PSF subtraction algorithms, and approaches to measuring astrometry, all of which have the potential to introduce systematic errors in the final astrometry. Linear fits generally provide good matches to the data (${\chi }_{\nu }^{2}$ ≈ 1), but in a few instances, the reduced χ² value is unreasonably high, indicating biased astrometry or underestimated errors. DH Tau is an especially acute example: ${\chi }_{\nu }^{2}$ = 6.8 for $\rho (t)$ and ${\chi }_{\nu }^{2}$ = 96 for θ(t). The source of this scatter is most likely instrument-to-instrument calibration errors in distortion correction, plate scale, and north alignment.

To mitigate these biases for our orbit fits, we introduce a "jitter" term added in quadrature with the astrometric uncertainties (σ_jit,ρ for separation and σ_jit,θ for P.A.). This effectively increases the astrometric errors in a systematic fashion until they are consistent with the expected level due to random scatter about linear evolution in separation and P.A. with time (Figures 7 and 8). For fits with ${\chi }_{\nu }^{2}$ > 1, jitter is derived by iteratively increasing σ_jit,ρ and σ_jit,θ until the linear fits result in ${\chi }_{\nu }^{2}$ = 1. Most systems do not require any additional additive error term. When needed, the typical jitter level for separation is ≈2–6 mas, and for the P.A. it is ≈01–04. For DH Tau, values of σ_jit,ρ (4.9 mas) and σ_jit,θ (074) imply that especially strong systematic errors are present in this data set. Linear fits of the separation and P.A. over time incorporating this added uncertainty are provided in Table 2. Residuals of these fits do not show any trends, suggesting the primary source of the high ${\chi }_{\nu }^{2}$ values is systematic astrometric errors or underestimated uncertainties, rather than our assumption of linear evolution for curved orbital motion (acceleration).

Table 2.  Linear Evolution of Separation and P.A.

Target	Sep. a0	Sep. a1	Sep. rms	σjit,ρ	Sep.	P.A. b0	P.A. b1	P.A. rms	σjit,θ	P.A.
Name	(mas)	(mas yr−1)	(mas)	(mas)	${\chi }_{\nu }^{2}$	(°)	(° yr−1)	(°)	(°)	${\chi }_{\nu }^{2}$
HD 984 B	−24198.3242 ± 6652.9297	12.1124 ± 3.3012	13.1	3.2	0.99	17555.8809 ± 867.8156	−8.6683 ± 0.4306	0.85	0.00	0.49
HD 1160 B	−2723.1067 ± 1785.6223	1.7393 ± 0.8863	9.7	0.0	0.10	261.4736 ± 72.9082	−0.0084 ± 0.0362	0.71	0.24	0.99
HD 19467 B	12941.7510 ± 2004.5902	−5.6069 ± 0.9957	5.6	1.8	0.99	1304.3900 ± 107.5698	−0.5277 ± 0.0534	0.21	0.06	0.99
1RXS0342+1216 B	22939.0391 ± 327.6723	−10.9844 ± 0.1632	5.3	0.0	0.94	−378.7789 ± 63.8557	0.1975 ± 0.0317	0.45	0.23	0.96
51 Eri b	6720.5132 ± 2516.3901	−3.1086 ± 1.2477	3.8	0.0	0.31	11322.8174 ± 343.0927	−5.5345 ± 0.1701	0.39	0.00	0.26
CD-35 2722 B	57221.8008 ± 1387.3359	−26.9053 ± 0.6888	4.3	3.5	0.98	849.6935 ± 127.4335	−0.3016 ± 0.0633	0.35	0.40	1.00
HD 49197 B	13782.7002 ± 843.4277	−6.4062 ± 0.4196	14.4	4.6	0.98	327.0609 ± 44.3513	−0.1244 ± 0.0220	0.35	0.00	0.81
HR 2562 B	−50161.6797 ± 2145.8904	25.1873 ± 1.0638	1.2	0.0	0.28	−202.7896 ± 257.1518	0.2482 ± 0.1275	0.17	0.00	0.44
HR 3549 B	17045.1328 ± 9505.5195	−8.0338 ± 4.7161	0.9	0.0	0.02	1324.7207 ± 498.3914	−0.5798 ± 0.2473	0.29	0.00	0.20
HD 95086 b	−1867.0281 ± 1941.0756	1.2346 ± 0.9629	5.2	0.0	0.31	2351.5354 ± 193.0089	−1.0930 ± 0.0957	0.32	0.00	0.32
GJ 504 B	−1710.9081 ± 1831.3406	2.0867 ± 0.9085	10.4	0.6	0.99	2354.0681 ± 51.8511	−1.0076 ± 0.0257	0.39	0.03	1.00
HIP 65426 b	6849.5601 ± 2470.6326	−2.9850 ± 1.2246	8.0	0.0	0.82	511.8435 ± 213.0702	−0.1794 ± 0.1056	0.55	0.12	0.98
PDS 70 b	2706.3250 ± 3142.1650	−1.2465 ± 1.5588	6.2	0.0	0.45	5107.3291 ± 817.6153	−2.4576 ± 0.4055	1.35	0.00	0.28
PZ Tel Ba	−58956.1328 ± 835.7415	29.5029 ± 0.4153	6.0	3.8	0.98	326.9745 ± 62.2090	−0.1328 ± 0.0309	0.47	0.16	0.98
HD 206893 B	16907.2207 ± 2168.6890	−8.2524 ± 1.0749	2.4	0.0	0.81	18744.5137 ± 537.0878	−9.2639 ± 0.2663	0.64	0.29	0.98
κ And B	50347.6562 ± 3829.7366	−24.4935 ± 1.9010	10.6	5.2	0.99	2160.1741 ± 167.6283	−1.0457 ± 0.0832	0.38	0.00	0.63
HD 23514 B	1142.6167 ± 531.7578	0.7459 ± 0.2644	6.7	0.0	0.39	346.3043 ± 28.8196	−0.0591 ± 0.0143	0.33	0.05	0.88
DH Tau B	1396.2018 ± 757.6589	0.4703 ± 0.3766	5.7	4.9	0.98	136.6982 ± 81.7527	0.0012 ± 0.0407	0.69	0.74	0.99
2M1559+4403 B	12130.2617 ± 2763.0696	−3.2307 ± 1.3737	25.6	10.5	1.00	359.0865 ± 37.8433	−0.0371 ± 0.0188	0.24	0.00	0.82
TWA 5 B	12507.5928 ± 623.8889	−5.2820 ± 0.3101	9.2	2.3	0.99	1161.5425 ± 68.8053	−0.4009 ± 0.0342	0.63	0.58	1.00
1RXS2351+3127 B	788.8538 ± 609.4363	0.7955 ± 0.3026	1.6	0.3	0.99	352.0048 ± 30.3868	−0.1293 ± 0.0151	0.07	0.05	0.90

Notes. This table provides relations that describe the linear evolution of separation and P.A. as a function of calendar year (t): ρ(t) = a₀ + a₁t and θ(t) = b₀ + b₁t. These incorporate astrometric jitter terms (σ_jit,ρ and σ_jit,θ), which are added in quadrature with the quoted astrometric errors to ensure that the reduced χ² values are ≤1.0.

^aNote that PZ Tel B shows signs of slight curvature, which is not reflected in this linear fit.

Download table as:  ASCIITypeset image

Metchev & Hillenbrand (2004) discovered this companion as part of their Palomar AO search for substellar companions to young stars in the solar neighborhood. They measure a spectral type of L4 ± 1 from a K-band spectrum and infer a mass of ≈63 M_Jup for HD 49197 B based on its absolute magnitude. We assume a stellar mass of 1.11 ± 0.06 M_⊙ from Mints & Hekker (2017), a total system mass of 1.17 ± 0.06 M_⊙, and a parallax of 23.99 ± 0.05 mas from Gaia DR2 in our orbit analysis.

Two epochs of astrometry were obtained by Metchev & Hillenbrand (2004) in 2002 and 2003, at which point the separation was 095 and P.A. was ≈78°. Additional epochs were obtained in 2006 by Serabyn et al. (2009) and in 2015 by Bottom et al. (2017), who first noted orbital motion toward the star. We have been monitoring this companion at a low cadence since 2011 and confirm this motion to smaller separations at a rate of −6.4 mas yr⁻¹ based on four epochs with HiCIAO and NIRC2. We also detect significant evolution in P.A. for the first time at a rate of −012 yr⁻¹ in the clockwise direction. Significant astrometric jitter is required to lower ${\chi }_{\nu }^{2}$ to ≈1 for the separation measurements (σ_jit,ρ = 4.6 mas). Overall, HD 49197 B has moved inward by about 01 and moved by about 2° over the past 15 yr.

Our best-fitting orbital solution has a semimajor axis of ${29}_{-10}^{+7}$ au and an orbital period of ${150}_{-70}^{+50}$ yr (2σ credible interval: 80–450 yr). This implies that ≈10% of its orbit has currently been mapped. The eccentricity of HD 49197 B is poorly constrained. Continued monitoring will be needed to refine the orbital elements for this companion. Results from the orbit fit are shown in Figure 9 and Appendix B.

Zoom In Zoom Out Reset image size

Wahhaj et al. (2011) identified this young L4 ± 1 brown dwarf companion as part of the Gemini NICI Planet-Finding Campaign (Liu et al. 2010). The host is an M1 member of the ≈120 Myr AB Dor association, implying a mass of 31 ± 8 M_Jup for CD–35 2722 B based on its luminosity and young age. We adopt a host mass of 0.40 ± 0.05 M_⊙ from Wahhaj et al. (2011), a total system mass of 0.43 ± 0.05 M_⊙, and a parallax of 44.635 ± 0.026 mas from Gaia DR2 in our orbit analysis.

Wahhaj et al. presented two epochs from 2009 and 2010, but no additional astrometry has been reported since then. Our epoch from 2018 reveals significant orbital motion at a rate of −27 mas yr⁻¹ toward the host star in separation and −03 yr⁻¹ in P.A. in the clockwise direction (Figure 7). A modest level of additional jitter is needed—σ_jit,ρ = 3.5 mas and σ_jit,θ = 040—although we note that with only three epochs of data, it is difficult to precisely determine the magnitude of any systematic error.

Our orbit fit is shown in Figure 9 and Appendix B and summarized in Table 3. We find a semimajor axis of ${240}_{-130}^{+90}$ au, a high eccentricity of ≈0.94, and an orbital period of ${5400}_{-3900}^{+2900}$ yr (2σ credible interval: 1000–32,000 yr). These results are largely consistent with the orbit constraints from Blunt et al. (2017), which only made use of the 2009 and 2010 epochs, but our new epoch now excludes low-eccentricity solutions.

GJ 504 B was discovered by Kuzuhara et al. (2013) as part of the Strategic Exploration of Exoplanets and Disks with the Subaru Telescope direct imaging survey (Tamura 2016). Depending on the system age (21 ± 2 Myr or 4.0 ± 1.8 Gyr) and assuming hot-start evolutionary models, this T8–T9.5 companion has a mass that may be as low as ≈1 M_Jup or as high as ≈33 M_Jup (Bonnefoy et al. 2018). For this study we assume a stellar host mass of 1.10–1.25 M_⊙ from Bonnefoy et al. (2018), a companion mass of ≈23 M_Jup from Bonnefoy et al. (2018), a total system mass of 1.20 ± 0.04 M_⊙, and a parallax of 57.02 ± 0.25 mas from Gaia DR2.

Our observations continue the astrometric trends noted by Kuzuhara et al. (2013) and Bonnefoy et al. (2018) from their data spanning 2011 to 2017: GJ 504 B is slowly moving away from the host star at a rate of 2.1 mas yr⁻¹, while most of the motion is in the P.A. with a rate of −10 yr⁻¹ in the clockwise direction. Modest jitter levels of σ_jit,ρ = 0.6 mas and σ_jit,θ = 003 are added in quadrature with the available astrometry for this system. Our orbit fit for GJ 504 B indicates a semimajor axis of ${41}_{-10}^{+4}$ au, a modest eccentricity of ${0.22}_{-0.17}^{+0.11}$, and a period of ${240}_{-80}^{+40}$ yr (2σ credible interval: 140–470 yr). Our posteriors generally resemble those from Blunt et al. (2017) and Bonnefoy et al. (2018), with slightly improved constraints due to the extended astrometric coverage. Results from the orbit fit are shown in Figure 9 and Appendix B.

HD 19467 B was first imaged by Crepp et al. (2014) based on a long-term radial acceleration of its Sun-like host star. Photometry and spectroscopy from Crepp et al. (2015) imply a mid-T spectral type and a mass of at least 52 M_Jup. We assume a stellar mass of 0.95 ± 0.02 M_⊙ from Crepp et al. (2014), a total system mass of 1.00 ± 0.02 M_⊙, and a parallax of 31.226 ± 0.041 mas from Gaia DR2 in our orbit analysis.

Published astrometry from Crepp et al. (2014) and Crepp et al. (2015) between 2011 and 2014 showed clockwise orbital motion and motion toward the host star. Our new astrometry in 2018 confirms this trend, revealing evolution in P.A. at a rate of −05 yr⁻¹ and in separation at a rate of −5.6 mas yr⁻¹. A modest amount of additional jitter (σ_jit,ρ = 1.8 mas, σ_jit,θ = 006) is needed to bring the ${\chi }_{\nu }^{2}$ to unity for a linear evolution of the separation and P.A. HD 19467 B has an orbital period of ${420}_{-250}^{+170}$ yr (2σ credible interval: 160–1530 yr) with a semimajor axis of ${56}_{-25}^{+15}$ au. It does not appear to have a high eccentricity (e < 0.8). Results from the orbit fit to the astrometry for this system are shown in Figure 9 and Appendix B.

Nielsen et al. (2012) identified this low-mass companion to the A0 star HD 1160 as part of the Gemini-NICI Planet-Finding Campaign, along with the low-mass stellar companion HD 1160 C at a wider separation (52). Garcia et al. (2017) found a mid-M spectral type for HD 1160 B and a mass range between 35 and 90 M_Jup depending on the system age. Recently, Curtis et al. (2019) found that this system may belong to the proposed ≈120 Myr Psc-Eri stream, which stretches 120° across the sky and about 400 pc in space. At this age, the implied mass for the companion HD 1160 B would be about 0.12 M_⊙. However, details about the origin, metallicity, membership probabilities, and potential for age gradients in this new proposed stream have not yet been established. For this study we adopt the more uncertain age and mass constraints for HD 1160 B from Nielsen et al. (2012) and Garcia et al. (2017), but caution that this companion may reside above the hydrogen-burning limit if the age is indeed older. The mass we adopt for HD 1160 A is 1.95 ± 0.05 M_⊙ (Blunt et al. 2017), the total system mass we use is 2.01 ± 0.06 M_⊙, and the Gaia DR2 parallax is 7.942 ± 0.076 mas.

Astrometry from 2002 to 2018 is presented by Nielsen et al. (2012), Maire et al. (2016), and Currie et al. (2018). Our new epoch in 2018 is consistent with that of Currie et al. (2018), which was taken about two months earlier. We find slow motion away from HD 1160 A at a rate of 1.7 mas yr⁻¹, but no significant change in P.A. over time. There is no evidence that the separation measurement errors are underestimated, but we find that a jitter level of 024 is needed to inflate the P.A. measurements. Our orbit fit (Figure 9 and Appendix B) implies a semimajor axis of ${80}_{-30}^{+20}$ au and an orbital period of ${520}_{-270}^{+200}$ yr (2σ credible interval: 230–3700 yr). The eccentricity is poorly constrained. These results are consistent with those from Blunt et al. (2017).

κ And B (Carson et al. 2013) is a substellar companion orbiting a young B9 star. The companion mass falls near the brown dwarf/planetary-mass boundary at 22 ± 9 M_Jup. Follow-up photometry and spectroscopy from Hinkley et al. (2013), Bonnefoy et al. (2014b), and Currie et al. (2018) indicate an early-L dwarf with red colors similar to low-gravity planets and brown dwarfs. We adopt a host star mass of 2.8 ± 0.1 M_⊙ (Jones et al. 2016), a total system mass of 2.82 ± 0.10 M_⊙, and a parallax of 19.975 ± 0.342 mas from Gaia DR2.

Our astrometry from 2016 combined with measurements from the literature obtained between 2012 and 2018 reveal rapid motion toward the star at a rate of −24 mas yr⁻¹ and P.A. evolution at a rate of −10 yr⁻¹ in the clockwise direction. Significant jitter is required for the separation measurements (σ_jit,ρ = 5.2 mas), but no jitter is needed for the P.A. uncertainties. Our orbit fit is shown in Figure 9 and the corner plot is displayed in Appendix B; we find a semimajor axis of ${80}_{-30}^{+20}$ au, an orbital period of ${420}_{-210}^{+150}$ yr (2σ credible interval: 170–1300 yr), and a high eccentricity of ${0.74}_{-0.08}^{+0.10}$. These results are consistent with the orbit fit from Currie et al. (2018).

This 35 ± 8 M_Jup companion to the M4 dwarf 1RXS J034231.8+121622 (2MASS J03423180+1216225) was discovered by Bowler et al. (2015a) as part of the Planets Around Low-Mass Stars (PALMS) survey. Bowler et al. (2015b) determine a spectral type of L0 from near-infrared spectroscopy and identified significant orbital motion. For our orbit analysis we adopt a host mass of 0.20 ± 0.05 M_⊙ (Bowler et al. 2015b), a total system mass of 0.23 ± 0.05 M_⊙, and parallax of 30.308 ± 0.067 mas from Gaia DR2.

Our new epoch was obtained in 2018, significantly extending the orbital coverage since the last published epoch acquired in 2013. 1RXS0342+1216 B continues to approach its host star at a rate of −11.0 mas yr⁻¹ in a counterclockwise direction at 02 yr⁻¹. We find that a jitter term of σ_jit,θ = 023 is needed to inflate the P.A. uncertainties. Our orbit fit is shown in Figure 10 and Appendix B; we find a semimajor axis of ${37}_{-14}^{+8}$ au and an orbital period of ${470}_{-240}^{+150}$ yr (2σ credible interval: 160–1040 yr). At present, the orbital eccentricity is effectively unconstrained.

Zoom In Zoom Out Reset image size

This substellar companion to the unusually dusty Pleiad HD 23514 was first identified by Rodriguez et al. (2012). Bowler et al. (2015b) determined a spectral type of M8 for HD 23514 B from H- and K-band spectroscopy. We adopt a host star mass of 1.42 ± 0.15 from Huber et al. (2016), a total system mass of 1.48 ± 0.15 M_⊙, and parallax of 7.210 ± 0.054 mas from Gaia DR2.

When combined with published astrometry by Rodriguez et al. (2012) and Yamamoto et al. (2013), our new astrometry from 2018 establishes the first signs of orbital motion for HD 23514 B: the separation is increasing by 0.7 mas yr⁻¹ and the P.A. is evolving by −006 yr⁻¹ in the clockwise direction. Over the ≈11 year baseline since the first detection of this system in 2006, this amounts to a change of about 8 mas and 07. No jitter is necessary for the separation measurements, and only a slight amount (005) is needed for the P.A.s. Results from the orbit fit are shown in Figure 10 and Appendix B. Despite the limited orbit coverage, eccentricities above ≈0.8 can be ruled out, and current astrometry favors modest values near 0.3. The semimajor axis is ${420}_{-160}^{+100}$ au with a long and highly uncertain orbital period of ${7000}_{-3900}^{+2500}$ yr (2σ credible interval: 2100–27,400 yr).

Itoh et al. (2005) discovered this brown dwarf companion to DH Tau, a young (≈2 Myr) M1 star in Taurus. Bonnefoy et al. (2014a) determine a spectral type of M9.25 and infer a mass of 8–21 M_Jup. Ginski et al. (2014) examined astrometry for this system, but did not find signs of orbital motion at that time. For our orbit fit we adopt a host mass of 0.64 ± 0.04 M_⊙ from Kraus & Hillenbrand (2009), a total system mass of 0.65 ± 0.04 M_⊙, and a parallax of 7.388 ± 0.069 mas from Gaia DR2.

Our new observations from 2018 add to extensive published astrometry for this system, dating back to 1999. However, these published values are highly discrepant—especially among the P.A. measurements—most likely indicating significant instrument-to-instrument calibration errors or underestimated uncertainties. This is evidenced by the excess astrometric jitter needed to bring the ${\chi }_{\nu }^{2}$ value for a linear fit to unity: 4.9 mas for the separations and 074 for the P.A.s. We find modest evidence for slight motion away from the host star at a rate of 0.5 mas yr⁻¹, but no signs of any change in P.A. over time. Nevertheless, these data offer some broad constraints on the orbital properties of DH Tau B (Figure 10 and Appendix B). The semimajor axis is ${330}_{-160}^{+90}$ au and the orbital period is ${7500}_{-4800}^{+3200}$ yr (2σ credible interval: 2400–60,000 yr). The eccentricity is unconstrained.

2MASS J15594729+4403595 B is a widely separated (56; 250 au) brown dwarf companion to a young M1.5 star and was first identified by Janson et al. (2012). Bowler et al. (2015b) measured strong lithium absorption in the optical spectrum of the companion and found a spectral type of M7.5. The system age (≲200 Myr) and luminosity of 2M1559+4403 B imply a mass of 43 ± 9 M_Jup. The host is likely a single-lined spectroscopic binary (Bowler et al. 2015b), which is bolstered by the strong astrometric excess noise (1.3 mas) reported from the Gaia DR2 astrometric fit. For this work we adopt a host mass of 0.54 ± 0.10 M_⊙ from Muirhead et al. (2018), a total system mass of 0.58 ± 0.10 M_⊙, and the Gaia DR2 parallax of the host (22.340 ± 0.218 mas) for our orbit analysis.

Janson et al. (2012), Bowler et al. (2015b), Janson et al. (2014), and Bowler et al. (2015b) presented astrometry for this companion spanning 2008 to 2014, which did not reveal signs of orbital motion. We find evidence for orbital motion toward the host star in separation with our new epoch from 2018 at a rate of −3.2 mas yr⁻¹. 2M1559+4403 B also shows changes in P.A. at the 2σ level with a rate of −004 yr⁻¹. Astrometric jitter is required for the separation measurement uncertainties (σ_jit,ρ = 10.5 mas), but not for the P.A. errors. Our orbit constraints for this system are presented in Figure 10 and Appendix B. 2M1559+4403 B has a semimajor axis of ${290}_{-140}^{+80}$ au and an orbital period of ${6500}_{-4300}^{+2800}$ yr (2σ credible interval: 1800–33500 yr). Its eccentricity is poorly constrained, but values above ∼0.9 are disfavored.

TWA 5 is a young (≈10 Myr) triple system comprising a close (≈3 au) stellar binary, TWA 5 Aab (Macintosh et al. 2001), and a wider brown dwarf companion located at ≈80 au (Lowrance et al. 1999; Webb et al. 1999). The 6 yr orbit of TWA 5 Aab has been well mapped over the past two decades (Konopacky et al. 2007; Köhler et al. 2013), but no orbit determination has been made for TWA 5 B. TWA 5 Aab is unresolved in our observations in 2018, but these data combined with published astrometry dating back to 1998 provide a 20-year baseline to assess the orbit of TWA 5 B.¹⁰ For our orbit analysis we adopt the dynamical mass of 0.9 ± 0.1 M_⊙ from Köhler et al. (2013) for the inner binary TWA 5 Aab, a total system mass of 0.92 ± 0.10 M_⊙, and a parallax of 20.252 ± 0.059 mas from Gaia DR2.

TWA 5 B is approaching the host pair at a rate of −5.3 mas yr⁻¹ and we measure evolution in P.A. at a rate of −04 yr⁻¹ in the clockwise direction, which takes into account astrometric jitter levels of σ_jit,ρ = 2.3 mas and σ_jit,θ = 058. The orbit of TWA 5 B about Aab is shown in Figure 10 and Appendix B, and the orbital elements are summarized in Table 3. We find a semimajor axis of ${150}_{-50}^{+30}$ au, a moderate eccentricity of ${0.43}_{-0.10}^{+0.10}$, and an orbital period of ${1810}_{-840}^{+540}$ yr (2σ credible interval: 740–4400 yr). The orbital inclination of TWA 5 B about Aab is ${151}_{-13}^{+12}^\circ$, which is misaligned with the measured orbital inclination of TWA 5 Aab (975 ± 01) from Köhler et al. (2013) at the 4σ level. It therefore appears that the orbits of TWA 5 Aab and TWA 5 B about Aab are not coplanar.

Ross 458 is a triple system comprising a close binary, Ross 458 AB (Heintz 1990), and a comoving late-T dwarf planetary-mass companion at a projected separation of about 1200 au (Goldman et al. 2010; Scholz 2010). The system age is estimated to be between 150 and 800 Myr (Burgasser et al. 2010). From its absolute magnitude of M_K ≈ 9.7 (Beuzit et al. 2004), Ross 458 B resides below the hydrogen-burning limit based on the Baraffe et al. (2015) evolutionary models if its age is younger than about 300 Myr. If the age is closer to 800 Myr, then its implied mass is about 0.1 M_⊙.

Heintz (1990) first detected unresolved astrometric perturbations of Ross 458 with an orbital period of 13.5 yr. This system has since been monitored with AO for more than a complete orbital cycle. We adopt a host star mass of 0.61 ± 0.03 M_⊙ from Neves et al. (2013), a total system mass of 0.68 ± 0.03 M_⊙, and a parallax of 86.857 ± 0.152 mas from Gaia DR2. Our observations in 2016 and 2018 supplement astrometry from Mann et al. (2019) and Ward-Duong et al. (2015), which together extend the astrometric coverage back to 2000. The orbit for Ross 458 B is therefore very well determined (see Figure 10 and Appendix B); we measure an orbital period of ${13.528}_{-0.014}^{+0.02}$ yr (2σ credible interval: 13.49–13.56 yr), a semimajor axis of 4.95 ± 0.01 au, and an orbital eccentricity of 0.242 ± 0.001.

1RXS J235133.3+312720 B is an early-L type substellar companion discovered as part of the PALMS survey (Bowler et al. 2012). The M2 host is a member of the ≈120 Myr AB Dor moving group (Shkolnik et al. 2012; Schlieder et al. 2012; Malo et al. 2013); at this age, the inferred mass of 1RXS2351+3127 B is 32 ± 6 M_Jup. We adopt a host mass of 0.45 ± 0.05 M_⊙ from Bowler et al. (2012), a total system mass of 0.48 ± 0.05 M_⊙, and parallax of 23.218 ± 0.052 mas from Gaia DR2.

Bowler et al. (2015b) detected hints of orbital motion based on observations spanning 2011 to 2013. Our epoch from 2019 clearly establishes evolution in the astrometry (Figure 10). 1RXS2351+3127 B is moving away from its host star at a rate of 0.8 mas yr⁻¹ in a clockwise direction at −013 yr⁻¹. We find jitter terms of σ_jit,ρ = 0.3 mas and σ_jit,θ = 005. 1RXS2351+3127 B has a semimajor axis of ${100}_{-40}^{+30}$ au and an orbital period of ${1490}_{-870}^{+630}$ yr (2σ credible interval: 570–10100 yr). Eccentricities above ∼0.8 are disfavored (Appendix B).

We have constructed a large (and to our knowledge complete) sample of 125 imaged substellar companions to stars. This list draws from previous catalogs of low-mass companions from Zuckerman & Song (2009), Faherty et al. (2010), Deacon et al. (2014), and Bowler (2016), together with additional discoveries compiled from the literature. Orbital motion is detected for 36 companions based on multi-epoch astrometry, both from the literature and from our new observations in this work (Figure 11).

The following criteria are used to define our samples of brown dwarfs and giant planets, with the goal of identifying targets that are most representative of the initial dynamical conditions of this population and least likely to have been influenced by significant orbital migration as a result of dynamical encounters with a third body. Companions must have projected separations between 5 and 100 au at the time of discovery and have measured orbital motion from multi-epoch astrometry. The hosts must be stars (>75 M_Jup), rather than brown dwarfs, to prevent biasing the samples toward binary star-like mass ratio distributions. Because the eccentricities of substellar companions can also be influenced by post-formation interactions in hierarchical triple systems (e.g., Fabrycky & Tremaine 2007; Allen et al. 2012; Reipurth & Mikkola 2015), we require the host stars not be known binaries. We similarly limit the sample to companions that are themselves single. Companions whose existence or characteristics are a matter of ongoing debate—in particular, LkCa15 bcd (Kraus & Ireland 2012; Ireland & Kraus 2014; Sallum et al. 2015; Thalmann et al. 2016; Currie et al. 2019) and HD 100546 bc (e.g., Quanz et al. 2013; Currie et al. 2015; Rameau et al. 2017)—are excluded from this analysis. This amounts to 27 systems that represent our full sample of directly imaged substellar companions undergoing orbital motion. Finally, we isolate two main subsamples within 5–100 au based on their masses as inferred from hot-start evolutionary models: 18 brown dwarfs spanning 15–75 M_Jup , and 9 giant planets between 2 and 15 M_Jup. Characteristics of the host stars and companions are summarized Table 4.

Table 4.  Sample of Substellar Companions between 2 and 75 M_Jup and 5–100 au Exhibiting Orbital Motion

Target	α	δ	Host	Age	πa	M*	Mcompb	M2/M1	Comp.	Sep.	Sep.	IWA/ac	Δt/Pd	Ref.
Name	(J2000.0)	(J2000.0)	SpT	(Myr)	(mas)	(M⊙)	(MJup)	(×10−2)	SpT	('')	(au)
HD 984 B	00 14 10.25	−07 11 56.8	F7	30–200	21.7806 ± 0.0564	1.15 ± 0.06	34–94	${5.3}_{-1.4}^{+2.0}$	M6.5 ± 1.5	0.2	10	0.26	0.05	1, 2, 3
HD 1160 B	00 15 57.30	+04 15 04.0	A0	80–125	7.9417 ± 0.0764	1.95 ± 0.10	35–90	${3.1}_{-0.8}^{+1.0}$	$M{5.5}_{-0.5}^{+1.0}$	0.8	100	1.18	0.03	4, 5, 6, 7
HD 4113 C	00 43 12.60	−37 58 57.47	G5	${5000}_{-1700}^{+1300}$	23.8531 ± 0.0536	1.05 ± 0.03	66 ± 5	6.0 ± 0.5	T9 ± 1	0.5	22	0.23	0.008	8
HD 4747 B	00 49 26.76	−23 12 44.9	G9	${3300}_{-1900}^{+2300}$	53.1836 ± 0.1264	0.82 ± 0.08	66 ± 3	7.7 ± 0.9	L9–T1	0.6	11	0.56	0.27	9, 10, 11
HD 19467 B	03 07 18.57	−13 45 42.4	G3	4600–10000	31.2255 ± 0.0410	0.95 ± 0.02	${57}_{-7}^{+5}$	5.7 ± 0.6	T5.5 ± 1	1.6	51	0.47	0.02	12, 13
1RXS0342+1216 B	03 42 31.80	+12 16 22.5	M4	60–300	30.3075 ± 0.0668	0.20 ± 0.05	35 ± 8	17 ± 5	L0 ± 1	0.8	19	0.12	0.02	14, 15
51 Eri b	04 37 36.13	−02 28 24.7	F0	23 ± 3	33.5770 ± 0.1354	1.75 ± 0.05	2 ± 1	0.11 ± 0.05	T4–T8	0.45	13	0.71	0.11	16, 17, 18, 19
β Pic b	05 47 17.09	−51 03 59.4	A6	23 ± 3	50.6231 ± 0.3339	1.84 ± 0.05	13 ± 3	0.7 ± 0.16	L1 ± 1	0.4	9	0.33	0.53	19, 20, 21, 22
CD–35 2722 B	06 09 19.21	−35 49 31.1	M1	${133}_{-20}^{+15}$	44.6346 ± 0.0262	0.40 ± 0.05	31 ± 8	7 ± 2	L4 ± 1	3.1	67	0.08	0.002	23, 24
Gl 229 B	06 10 34.62	−21 51 52.7	M1	2000–6000	173.6955 ± 0.0457	0.54 ± 0.04	72 ± 5	12.7 ± 1.3	T7p ± 0.5	6.8	39	0.79	0.05	25, 26, 27, 28
HD 49197 B	06 49 21.33	+43 45 32.8	F5	260–790	23.9903 ± 0.0486	1.11 ± 0.06	${63}_{-21}^{+13}$	5.4 ± 1.4	L4 ± 1	1	40	1.6	0.11	29, 30, 3
HR 2562 B	06 50 01.01	−60 14 56.9	F5	200–750	29.3767 ± 0.0411	1.37 ± 0.02	32 ± 14	2.2 ± 1.0	T2–T3	0.6	20	0.37	0.01	31, 32
HR 3549 B	08 53 03.78	−56 38 58.1	A0	100–150	10.4864 ± 0.0898	${2.3}_{-0.1}^{+0.2}$	40–50	1.9 ± 0.2	M9–L0	0.9	80	0.46	0.005	33, 34
HD 95086 b	10 57 03.01	−68 40 02.4	A8	17 ± 4	11.5684 ± 0.0325	1.7 ± 0.1e	4.4 ± 0.8	0.25 ± 0.05	L1–T3	0.6	52	0.79	0.02	35, 36, 37
GJ 504 Bf	13 16 46.52	+09 25 26.9	G0	4000 ± 1800	57.0186 ± 0.2524	1.10–1.25	${23}_{-9}^{+10}$	1.8 ± 0.8	T8–T9.5	2.5	44	0.79	0.03	38, 39
HIP 65426 b	13 24 36.09	−51 30 16.0	A2	14 ± 4	9.1566 ± 0.0626	1.96 ± 0.04	8 ± 1	0.39 ± 0.05	mid-L	0.8	90	0.36	0.003	40, 41
PDS 70 b	14 08 10.15	−41 23 52.57	K7	5 ± 1	8.8159 ± 0.0405	0.76 ± 0.02	4–10	0.9 ± 0.3	L:	0.2	22	0.76	0.05	42, 43
PZ Tel B	18 53 05.88	−50 10 49.9	K0	23 ± 3	21.2186 ± 0.0602	${1.25}_{-0.20}^{+0.05}$	38–72	${4.2}_{-0.8}^{+1.0}$	M7 ± 1	0.3	17	0.52	0.07	44, 45, 19, 5
Gl 758 B	19 23 34.01	+33 13 19.1	G8	6000–10000	64.0623 ± 0.0218	0.83 ± 0.11	${38.1}_{-1.5}^{+1.7}$	4.4 ± 0.6	T7–T8	1.8	25	0.26	0.04	46, 47, 48, 49, 11
HR 7672 B	20 04 06.22	+17 04 12.6	G0	3000–5000	56.4256 ± 0.0690	0.96 ± 0.05	72.7 ± 0.8	7.2 ± 0.4	L4.5 ± 1.5	0.8	14	0.63	0.11	50, 51, 11
HR 8799 b	23 07 28.72	+21 08 03.3	A5	40 ± 5	24.2175 ± 0.0881	1.48 ± 0.05	5 ± 1	0.32 ± 0.07	L/Tpec	1.7	68	0.48	0.03	52, 53, 54, 55
HR 8799c	23 07 28.72	+21 08 03.3	A5	40 ± 5	24.2175 ± 0.0881	1.48 ± 0.05	7 ± 1	0.45 ± 0.07	L/Tpec	0.95	38	0.90	0.08	52, 53, 54, 55
HR 8799 d	23 07 28.72	+21 08 03.3	A5	40 ± 5	24.2175 ± 0.0881	1.48 ± 0.05	7 ± 1	0.45 ± 0.07	L6–L8pec	0.6	24	0.75	0.16	52, 53, 54, 55
HR 8799 e	23 07 28.72	+21 08 03.3	A5	40 ± 5	24.2175 ± 0.0881	1.48 ± 0.05	7 ± 1	0.45 ± 0.07	L6–L8pec	0.4	14	0.54	0.13	52, 53, 54, 55, 56
HD 206893 B	21 45 21.90	−12 47 00.1	F5	${250}_{-200}^{+450}$	24.5062 ± 0.0639	1.32 ± 0.02	15–50	2.3 ± 0.9	L3–L5pec	0.3	10	0.72	0.10	57, 58
κ And B	23 40 24.51	+44 20 02.2	B9	${50}_{-40}^{+30}$	19.9751 ± 0.3418	2.8 ± 0.1	22 ± 9	0.8 ± 0.3	L0–L1	1.1	55	0.41	0.01	59, 60, 61
1RXS2351+3127 B	23 51 33.66	+31 27 22.9	M2	${133}_{-20}^{+15}$	23.2183 ± 0.0524	0.45 ± 0.05	32 ± 6	6.8 ± 1.5	L0 ± 1	2.4	100	0.26	0.005	62, 24

Notes.

^aParallaxes from Gaia DR2. ^bAssumes hot-start cooling history. ^cIWA/a refers to the inner working angle at the contrast of the companion and at the time of discovery, divided by the maximum a posteriori of the semimajor axis distribution. The IWA was visually estimated based on the discovery papers. This value—IWA/a—represents a useful metric to assess the potential impact of discovery bias on the inferred population-level eccentricity distribution (see Dupuy & Liu 2011 and Section 4.4). Values above 1.0, between 0.5 and 1.0, and below 0.5 indicate strong, moderate, and minimal discovery bias, respectively. ^dFractional orbital coverage, computed from the time baseline over which this companion has been imaged (Δt) and the best-fit orbital period (P). ^eUncertainty in the host mass is approximated from multiple assessments in the literature (Meshkat et al. 2013; Chen et al. 2014). ^fBonnefoy et al. (2018) find two degenerate solutions for the host age (21 ± 2 Myr or 4.0 ± 1.8 Gyr), which results in two solutions for the companion mass (${1.3}_{-0.3}^{+0.6}$ or ${23}_{-9}^{+10}$ M_Jup) and mass ratio (0.11 ± 0.04 × 10⁻² or 1.8 ± 0.8 × 10⁻²). For this study we adopt the older age and correspondingly higher companion mass.

References. (1) Meshkat et al. (2015); (2) Johnson-Groh et al. (2017); (3) Mints & Hekker (2017); (4) Nielsen et al. (2012); (5) Maire et al. (2016); (6) Garcia et al. (2017); (7) Blunt et al. (2017);(8) Cheetham et al. (2018); (9) Crepp et al. (2016); (10) Crepp et al. (2018); (11) Brandt et al. (2019a); (12) Crepp et al. (2014); (13) Crepp et al. (2015); (14) Bowler et al. (2015b); (15) Bowler et al. (2015a); (16) Macintosh et al. (2015); (17) Simon & Schaefer (2011); (18) Rajan et al. (2017); (19) Mamajek & Bell (2014); (20) Lagrange et al. (2010); (21) Dupuy et al. (2019); (22) Chilcote et al. (2017); (23) Wahhaj et al. (2011); (24) Gagné et al. (2018); (25) Nakajima et al. (1995); (26) Oppenheimer et al. (1995); (27) Burgasser et al. (2006); (28) Brandt et al. (2019b); (29) Metchev & Hillenbrand (2004); (30) Metchev & Hillenbrand (2009); (31) Konopacky et al. (2016); (32) Mesa et al. (2018); (33) Mawet et al. (2015); (34) Mesa et al. (2016); (35) Rameau et al. (2013b); (36) Meshkat et al. (2013); (37) De Rosa et al. (2014); (38) Kuzuhara et al. (2013); (39) Bonnefoy et al. (2018); (40) Chauvin et al. (2017); (41) Cheetham et al. (2019); (42) Keppler et al. (2018); (43) Müller et al. (2018); (44) Biller et al. (2010); (45) Mugrauer et al. (2010); (46) Thalmann et al. (2009); (47) Brewer et al. (2016); (48) Nilsson et al. (2017); (49) Bowler et al. (2018); (50) Liu et al. (2002); (51) Crepp et al. (2012); (52) Marois et al. (2008); (53) Wang et al. (2018); (54) Barman et al. (2015); (55) Bonnefoy et al. (2016); (56) Marois et al. (2010); (57) Milli et al. (2017); (58) Delorme et al. (2017); (59) Carson et al. (2013); (60) Jones et al. (2016); (61) Currie et al. (2018); (62) Bowler et al. (2012).

Download table as:  ASCIITypeset image

We compiled all available astrometry of targets in our full sample with small fractional orbit coverage (see Appendix A) to reassess their orbits in a uniform manner instead of relying on orbit determinations in the literature, which can be influenced by different priors and approaches to fitting short orbital arcs. For the remaining targets with well-constrained (or moderately well-constrained) orbits we adopt eccentricity posteriors from the literature. This includes the HR 8799 planets (bcde), for which we use the "unconstrained" eccentricity posteriors from Wang et al. (2018), Gl 229 B from Brandt et al. (2019b), and approximately Gaussian-shaped eccentricity distributions for the following five systems: β Pic b (e = 0.24 ± 0.06; Dupuy et al. 2019), HD 4113 C (e = 0.38 ± 0.06; Cheetham et al. 2018), HD 4747 B (e = 0.735 ± 0.003; Brandt et al. 2019a), Gl 758 B (e = 0.40 ± 0.09; Brandt et al. 2019a), and HR 7672 (e = 0.542 ± 0.018; Brandt et al. 2019a). Note that because the orbits of these systems are well constrained, the choice of priors does not meaningfully influence the posteriors. Figure 12 shows the distribution of fractional orbital coverage for our sample. Most systems have traced out <10% of their orbits. Only 11 have been imaged for over 5% of their orbital periods.

Zoom In Zoom Out Reset image size

For the remaining nine systems that were not already discussed in Section 3, we apply the same Bayesian priors and fitting approach as previously discussed. Astrometric jitter for each system is assessed and implemented by identifying the amplitude of excess noise, which, when added in quadrature to quoted errors in separation and separately for P.A., results in a linear fit with a reduced χ² value of 1.0 (Table 2). No additional uncertainty is added if ${\chi }_{\nu }^{2}$ is lower than 1.0 using the raw unadjusted errors.

Results from the orbit fits are shown in Figures 13 and 14 and summarized in Table 3. Our constraints are generally consistent with published orbits that make use of the same astrometry. Below we compare our fits with those in the literature with a particular focus on the eccentricity posteriors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

HD 984 B. This substellar companion was found by Meshkat et al. (2015) and was further characterized by Johnson-Groh et al. (2017) with Gemini/GPI as part of the GPIES survey. We adopt a stellar mass of 1.15 ± 0.06 M_⊙ from Mints & Hekker (2017), a total system mass of 1.21 ± 0.06 M_⊙, and a parallax of 21.781 ± 0.056 mas from Gaia DR2 for our orbit fit. Our orbital constraints are similar to those of Johnson-Groh et al. using the same five epochs from 2012 to 2016 (Figure 13 and Appdendix B). We include a jitter term of 3.2 mas in separation, but this does not significantly alter the results. Our eccentricity posterior peaks near 0.0 with a long tail out to e ≈ 0.8. Although this constraint is broad, circular orbits are clearly favored for this system.

51 Eri b. This ≈2 M_Jup planet was discovered by Macintosh et al. (2015) and has been reimaged several times with GPI and SPHERE since then. The M+M binary GJ 3305 AB orbits 51 Eri Ab at ≈2000 au (Feigelson et al. 2006; Montet et al. 2015). A total system mass of 1.75 ± 0.05 M_⊙ (Macintosh et al. 2015) and parallax of 33.577 ± 0.135 mas from Gaia DR2 are fixed for our orbit fit. We make use of 11 epochs for our analysis—five from De Rosa et al. (2015), four of which were originally presented in Macintosh et al. (2015), and six from Maire et al. (2019). No additional jitter is included for this system.

Maire et al. find hints of curvature in the orbit and note a small (≈1σ) offset in P.A. between GPI and SPHERE for observations taken at similar epochs. They proceeded to add a 10 offset to the GPI data to reduce this difference and also removed the epoch from 2015 January from consideration because the separation is somewhat larger than for other epochs. For this study we take the astrometry at face value; no recalibrations are applied, and we consider all published astrometry. Any slight differences in astrometry for individual epochs are at the ≈1σ level and are well accounted for with larger uncertainties. Indeed, linear fits to separation and P.A. over time yield ${\chi }_{\nu }^{2}$ values below unity (Table 2), suggesting that reported uncertainties may even be slightly overestimated for this system. Nevertheless, our results are very similar to those of Maire et al. (Figure 13 and Appendix B); we find a significant eccentricity for 51 Eri b of e = ${0.50}_{-0.08}^{+0.11}$, in good agreement with their median value of e = 0.45 and a 68% credible interval of e = 0.30–0.55.

HR 2562 B. Konopacky et al. (2016) presented four epochs of HR 2562 B with GPI in their discovery paper. They confirmed that the companion shares a common proper motion with the host star and found that its P.A. and the orientation of the disk are consistent with a nested orbit inside the disk. Maire et al. (2018) added six additional measurements over three epochs in 2016 and 2017 with SPHERE. They found orbital motion and confirmed that the orbital plane of HR 2562 B may be aligned with the debris disk. They also identified a wide range of eccentricities when all the observations were used, but values of e ≲ 0.3 were preferred when a prior imposed by the debris disk was included.

Our orbit fit using a stellar mass of 1.37 ± 0.02 M_⊙ (Mesa et al. 2018), a total system mass of 1.40 ± 0.02 M_⊙, and Gaia DR2 parallax of 29.377 ± 0.041 mas is broadly consistent with the results from Maire et al. All eccentricities are allowed (the 2σ credible interval spans e = 0.037–1.0; see Figure 13 and Appendix B). We also find that HR 2562 B is on a nearly edge-on orbit, which is apparent from the large radial motion of HR 2562 B with respect to its host star, whereas the P.A. is nearly unchanging. Note that no additional jitter is included in our orbit fit for this system. We do not incorporate additional constraints on orbital solutions from the debris disk to avoid having to make assumptions about disk-planet coplanarity for HR 2562 and other systems with disks in this study.

HR 3549 B. Limited astrometry has been published for the substellar companion HR 3549 B. Mawet et al. (2015) obtained two epochs with VLT/NaCo from 2013 and 2015, and Mesa et al. (2016) acquired two astrometric measurements taken at the same epoch in 2015 with SPHERE. Mesa et al. examined orbital constraints for HR 3549 B with and without constraints from the disk. They found that high eccentricities (e ≳ 0.3) are preferred for uninformed priors, and orbits informed by the presence of the inner disk suggest even higher values (e ≳ 0.5).

For this system we use a stellar mass of 2.3 ± 0.15 M_⊙ from Mawet et al. (2015), a total system mass of 2.34 ± 0.15 M_⊙, and a parallax of 10.486 ± 0.090 mas from Gaia DR2. No additional jitter is added to the published astrometry. Our orbit fit is shown in Figure 13 and Appendix B. Our results disagree with the eccentricity distribution from Mesa et al.; we find that all eccentricities are possible, and that the highest eccentricities above e ≈ 0.9 are least preferred (the 2σ credible interval is 0.0–0.88). The origin of this discrepancy is not immediately obvious, but may arise from the difference between the Bayesian orbit fitting of orbitize! compared to the least-squares Monte Carlo technique.

HD 95086 b. HD 95086 b was identified by Rameau et al. (2013b) and confirmed by Rameau et al. (2013a). This companion has been continuously monitored with GPI (Rameau et al. 2016) and SPHERE (Chauvin et al. 2018) since its discovery. Both of these studies found that the companion's modest orbital motion points to lower eccentricities (e ≲ 0.5), with the most likely values being near zero.

We adopt a stellar (and total system) mass of 1.7 ± 0.1 M_⊙ (see Table 4) and a parallax of 11.568 ± 0.033 mas from Gaia DR2. No additional jitter is needed for this companion. Our results using all 20 epochs of published astrometry (Figure 13 and Appendix B) are in good agreement with those of Rameau et al. and Chauvin et al.; our 2σ credible interval for the eccentricity is 0.0–0.48, with low values near zero being preferred.

HIP 65426 b. Chauvin et al. (2017) discovered this long-period planet with SPHERE and demonstrated that this object shares a common proper motion with its host star, but no orbital motion was detected. Cheetham et al. (2019) presented follow-up multi-epoch astrometry with VLT/NaCo and SPHERE that enabled them to carry out the first orbital analysis for this system. They found that the eccentricity is essentially unconstrained and is dominated by the choice in prior rather than the likelihood (the data).

Our results are shown in Figure 14 and Appendix B. We adopt a stellar mass of 1.96 ± 0.04 M_⊙ (Chauvin et al. 2017), a total system mass of 1.97 ± 0.04 M_⊙, and parallax of 9.157 ± 0.063 mas from Gaia DR2. We find similar results as Cheetham et al. using a modest addition of 012 of jitter in P.A. (but no adjustment to the uncertainties in separation): the eccentricity distribution is essentially unconstrained, and our 2σ credible spans 0.03 to 0.97 as a result of the sparse orbital coverage for this system.

PDS 70 b. This young planet was found by Keppler et al. (2018) nested in the transition disk of its host star using NICI, NaCo, and SPHERE as part of the SHINE survey. Müller et al. (2018) presented additional observations with SPHERE and carried out preliminary orbit constraints, finding an eccentricity posterior that peaks near e = 0.0 with a tail to a value of about 0.6. Additional astrometry was published by Wagner et al. (2018) using Magellan/MagAO, Christiaens et al. (2019) using VLT/SINFONI, and Haffert et al. (2019) using VLT/MUSE.

We adopt a stellar host mass of 0.76 ± 0.02 M_⊙ from Müller et al. (2018), a total system mass of 0.77 ± 0.02 M_⊙, and a parallax of 8.816 ± 0.041 mas from Gaia DR2. Our orbit fit using all 14 published epochs is shown in Figure 14 and Appendix B, and is summarized in Table 3. No jitter is required. Our 2σ credible interval spans e = 0.0–0.59, with lower values being preferred.

PZ Tel B. This brown dwarf companion was independently found by Biller et al. (2010) with NICI and Mugrauer et al. (2010) with NaCo. Despite having only two epochs available, Biller et al. were able to establish a high eccentricity of >0.6 for this system. This has been bolstered by additional astrometry and progressively more refined orbit constraints over the past decade (Mugrauer et al. 2012; Ginski et al. 2014; Beust et al. 2016; Maire et al. 2016).

We adopt a stellar host mass of 1.25 ± 0.10 M_⊙ from Biller et al. (2010), a total system mass of 1.3 ± 0.1 M_⊙, and a parallax of 21.219 ± 0.060 mas from Gaia DR2 for our orbit analysis. All 26 published epochs are used for our orbit fit. There is some indication that on average, the reported uncertainties are underestimated; 3.8 mas and 016 of astrometric jitter are required to bring the separation and P.A. measurements in statistical agreement with a linear fit. Note that we find significant curvature in both separation and P.A. (Figure 14), suggesting we may be slightly overestimating the value of the jitter for this star. This curvature was first noted by Ginski et al. (2014), Mugrauer et al. (2012), and Maire et al. (2016). With this additional adjustment, we find that high values of eccentricity above 0.6 are strongly preferred, with the 2σ credible interval spanning 0.74–1.0. Results for PZ Tel B can be found in Figure 14 and Appendix B.

HD 206893 B. Five astrometric epochs are available for this companion: two with SPHERE and NaCo as part of the discovery observations by Milli et al. (2017), and three additional observations with SPHERE by Delorme et al. (2017) and Grandjean et al. (2019). Delorme et al. determined initial constraints on the orbital elements for this companion and found that all eccentricities are allowed with a preference against the highest values above about 0.9. When we consider only orbits that are coplanar with the debris disk, this collapses to a range of about e = 0.0–0.4. More recently, Grandjean et al. (2019) combine relative astrometry, a radial acceleration, and astrometric acceleration measured between Hipparcos and Gaia to constrain the orbit and dynamical mass of HD 206893 B. They find a low eccentricity (e ≲ 0.4) and potential evidence of a second companion based on the larger-than-expected radial acceleration.

For this study we fix the host mass at 1.32 ± 0.02 M_⊙ from Delorme et al. (2017), adopt a total system mass of 1.35 ± 0.02 M_⊙, and use a parallax of 24.506 ± 0.064 mas from Gaia DR2. No jitter is required in separation, but modest jitter of 029 is inferred for the P.A. measurements. Results from our orbit fits are presented in Figure 14 and Appendix B. We determine an eccentricity of e = ${0.25}_{-0.14}^{+0.17}$ with a 2σ credible interval of 0.0–0.44. Our constraints are narrower than the general unrestricted fit from Delorme et al. and are similar to the recent results from Grandjean et al. (2019), who also included relative astromety, radial velocities, and absolute astrometry of the host star from Hipparcos and Gaia in their orbit fit.

Our goal is to determine the underlying behavior of the substellar eccentricity distribution at the population level given a sample of measured eccentricity posteriors for individual systems—some precisely determined, others broadly constrained, and most of which are asymmetric and non-Gaussian in shape (Figure 15 and Table 5). Incorporating the structure of the posterior distribution in these constraints is especially important in this study because the number of systems under consideration is relatively modest: 9 long-period planets and 18 brown dwarfs, totaling 27 substellar companions altogether.

Zoom In Zoom Out Reset image size

Hierarchical Bayesian modeling offers a natural framework to incorporate this type of probabilistic information at multiple (individual and population) levels. This tool is gaining popularity in astronomy (see, e.g., Loredo 2013) and particularly within the field of exoplanets; for example, it has been used to determine exoplanet eccentricity distributions (Hogg et al. 2010; Shabram et al. 2016; Van Eylen et al. 2019), the value of η_⊕ (Foreman-Mackey et al. 2014), the planet mass–radius relationship (Rogers 2015; Wolfgang et al. 2016), and host star obliquities (Morton & Winn 2014; Campante et al. 2016).

Hierarchical Bayesian modeling enables simultaneous inference for parameters of individual systems (${\boldsymbol{\theta }}$) and hyperparameters governing the underlying behavior of the population (${\boldsymbol{\Lambda }}$), given the data ${\boldsymbol{d}}$. Here variables in bold denote vectors of multiple values, parameters, or data sets. Bayes' theorem for this multi-level modeling becomes

$$\begin{eqnarray}&&p({\boldsymbol{\theta }},{\boldsymbol{\Lambda }}| {\boldsymbol{d}})=\displaystyle \frac{p({\boldsymbol{d}}| {\boldsymbol{\theta }})p({\boldsymbol{\theta }}| {\boldsymbol{\Lambda }})p({\boldsymbol{\Lambda }})}{p({\boldsymbol{d}})},\end{eqnarray} \tag{ 1 }$$

where $p({\boldsymbol{\theta }},{\boldsymbol{\Lambda }}| {\boldsymbol{d}})$ is the joint posterior distribution for the individual and population parameters, $p({\boldsymbol{d}}| {\boldsymbol{\theta }})$ is the likelihood function of the data, $p({\boldsymbol{\theta }}| {\boldsymbol{\Lambda }})$ is the prior on the individual systems conditioned on the set of population-level hyperparameters, $p({\boldsymbol{\Lambda }})$ is the hyper-prior—the prior distribution on the set of population-level parameters ${\boldsymbol{\Lambda }}$—and $p({\boldsymbol{d}})$ is the marginalized likelihood.

For this problem, the orbit fits are carried out separately, and the individual posteriors for the eccentricity distributions (and all other orbital elements) are available. We therefore seek to constrain the hyperparameters of a parameterized model for the underlying eccentricity distribution based on observations of N systems. The posterior probability distribution of ${\boldsymbol{\Lambda }}$ is simply

$$\begin{eqnarray}&&p({\boldsymbol{\Lambda }}| {\boldsymbol{d}})\propto { \mathcal L }({\boldsymbol{d}}| {\boldsymbol{\Lambda }})\pi ({\boldsymbol{\Lambda }}).\end{eqnarray} \tag{ 2 }$$

Here ${ \mathcal L }({\boldsymbol{d}}| {\boldsymbol{\Lambda }})$ is the likelihood function and $\pi ({\boldsymbol{\Lambda }})$ is the set of priors on the hyperparameters.

Hogg et al. (2010) describe an importance sampling approach to hierarchical Bayesian modeling with a specific application to exoplanet eccentricities. We follow their method by making use of the eccentricity posterior distributions of individual systems from our orbit fits to inform the population-level likelihood function for parameters in the underlying eccentricity distribution, $f(e| {\boldsymbol{\Lambda }})$. In this case, the sampling approximation to the likelihood function in Equation (2) is

$$\begin{eqnarray}&&{ \mathcal L }({\boldsymbol{d}}| {\boldsymbol{\Lambda }})\approx \displaystyle \prod _{n=1}^{N}\displaystyle \frac{1}{K}\displaystyle \sum _{k=1}^{K}\displaystyle \frac{f({e}_{{nk}}| {\boldsymbol{\Lambda }})}{\pi ({e}_{{nk}})},\end{eqnarray} \tag{ 3 }$$

where N is the number of substellar companions being considered, K is the number of samples from the posterior eccentricity distribution, e_nk is the kth random draw for each eccentricity distribution n, $f({e}_{{nk}}| {\boldsymbol{\Lambda }})$ is the probability density of the population-level eccentricity distribution evaluated at e_nk and conditioned on the hyperparameters ${\boldsymbol{\Lambda }}$, and π(e_nk) is the probability density of the prior probability distribution evaluated at e_nk.

We follow the approach of Hogg et al. (2010), Kipping (2013), and Van Eylen et al. (2019), by adopting a standard Beta distribution for our population-level eccentricity distribution, $f(e| {\boldsymbol{\Lambda }})$. The advantage of the standard Beta distribution functional form is that it is flexible, spans a range of [0, 1], and is described by only two shape parameters, ${\boldsymbol{\Lambda }}$ ≡ (α, β), both of which take on values >0:

$$\begin{eqnarray}&&f(e| \alpha ,\beta )=\displaystyle \frac{{\rm{\Gamma }}(\alpha +\beta )}{{\rm{\Gamma }}(\alpha ){\rm{\Gamma }}(\beta )}{e}^{\alpha -1}{\left(1-e\right)}^{\beta -1}.\end{eqnarray} \tag{ 4 }$$

Here Γ represents the Gamma function. The Beta distribution can capture a wide range of shapes, including uniform (α = 1, β = 1); ${ \mathcal U }$-shaped with a single anti-mode; unimodal with positive, negative, or no skew; ${ \mathcal J }$-shaped or reverse ${ \mathcal J }$-shaped; bell-shaped; and triangular. α governs the function's behavior at low eccentricities and β influences its shape at high eccentricities. Low values of α and β correspond to high probability densities near e = 0 and e = 1, while high values of these parameters correspond to low probability densities near zero.¹¹

This simple parameterization is especially convenient for this study because it can qualitatively reproduce a wide range of potential physical outcomes of the planet formation and migration process: outward scattering, which is expected to excite eccentricities (Scharf & Menou 2009; Veras et al. 2009); cloud fragmentation, which should result in a broad range of eccentricities; dynamically relaxed systems that follow the thermal eccentricity distribution (f(e) ∼ 2e; Ambartsumian 1937); and formation in a disk, in which case companions should retain nearly circular orbits if they are dynamically undisturbed.

Our primary goal is therefore to constrain the hyperparameters α and β of the Beta distribution for our sample of imaged substellar companions undergoing measured orbital motion. We also examine the eccentricity distribution of other subsamples: a division of giant planets and brown dwarfs based on companion mass, a subdivision by mass ratio, imaged planets including and excluding the HR 8799 system, and divisions based on orbital separation and age. These are discussed separately in more detail below.

For each of these cases we use the Metropolis–Hastings MCMC algorithm (Metropolis et al. 1953; Hastings 1970) to sample the posterior distributions of the model parameters α and β. Linearly uniform hyperpriors are chosen for α from 0 to 1000, β from 0 to 1000, and π(e) from 0 to 1. With only two parameters for the parent model and uncomplicated covariance, we use a single chain typically comprising 10⁶ links and K = 1000 samples from each individual system-level posterior to explore the population-level eccentricity posterior.¹² Convergence is monitored using the Gelman–Rubin (GR) statistic (Gelman & Rubin 1992), which compares the variance within chains to the variance between chains (here our single chain divided into subchains). In all cases the GR statistic is lower than 1.1 within 10⁵ links, and in most cases it is lower than 1.01, indicating that the chains are well mixed. No burn-in is warranted because of the low dimensionality of the model, and the starting points are all near the equilibrium point of the posterior distributions. We adopt normal proposal distributions with standard deviations set to avoid acceptance rates that are too high (near 1), in which case step sizes become too small to efficiently map posterior space, and too low (near 0) where large jumps mean that only few proposed values are accepted and convergence is slow. Most of our final acceptance rates fall between 0.3 and 0.8.

4.3.1. Example: Recovering the Short-period Exoplanet Eccentricity Distribution

We carried out a series of experiments using mock data sets drawn from the actual eccentricity distribution of warm Jupiters to establish how the individual-level eccentricity measurement precision (σ_e) and number of systems under consideration (N) influence the ability to constrain the shape parameters α and β. We adopt Beta parameters α = 0.867 and β = 3.03 from Kipping (2013) as underlying "truth" for this exercise and test four samples comprising 5, 10, 20, and 50 systems. For each sample we draw random eccentricity mock measurements from the population-level distribution and assign three sets of Gaussian uncertainties to these data sets, σ_e = 0.2, 0.1, and 0.05. These individual synthetic measurements are meant to mimic random sampling from the parent distribution and are truncated below e = 0 and above e = 1.

Results of these 12 experiments are summarized in Table 6, and one example with N = 20 and σ_e = 0.05 is shown in Figure 16. As expected, the fidelity with which the parameters of the true distribution are recovered strongly depends on the number of measurements as well as on the precision of each measurement. As more planets are "observed" and the constraints on the orbital eccentricity improve, the median value of the posteriors for α and β approach the true values and the uncertainties tend to decrease.

Zoom In Zoom Out Reset image size

Table 6.  Experiments Recovering the Warm Jupiter Eccentricity Distribution

Na	σeb	αc	βc
5	0.20	${15.94}_{-12.46}^{+7.86}$	${31.75}_{-8.08}^{+18.24}$
5	0.05	${2.03}_{-1.35}^{+0.83}$	${3.65}_{-2.33}^{+1.53}$
5	0.01	${2.23}_{-1.27}^{+0.85}$	${14.01}_{-8.35}^{+5.91}$
10	0.20	${8.13}_{-5.39}^{+3.60}$	${35.91}_{-6.59}^{+14.05}$
10	0.05	${4.72}_{-2.97}^{+1.83}$	${21.24}_{-12.87}^{+8.26}$
10	0.01	${1.97}_{-0.86}^{+0.73}$	${9.43}_{-4.73}^{+3.34}$
20	0.20	${8.69}_{-5.55}^{+3.99}$	${31.35}_{-9.05}^{+16.97}$
20	0.05	${1.44}_{-0.54}^{+0.46}$	${4.73}_{-1.86}^{+1.51}$
20	0.01	${0.80}_{-0.36}^{+0.25}$	${4.89}_{-2.32}^{+1.48}$
50	0.20	${8.60}_{-3.16}^{+3.04}$	${37.70}_{-6.23}^{+12.29}$
50	0.05	${1.25}_{-0.32}^{+0.23}$	${3.84}_{-0.97}^{+0.78}$
50	0.01	${0.66}_{-0.15}^{+0.14}$	${2.22}_{-0.59}^{+0.46}$

Notes.

^aNumber of draws from the underlying Beta distribution (α = 0.867, β = 3.03) from Kipping (2013). ^bMeasurement errors for each mock realization. σ_e denotes the standard deviation of the Gaussian probability distribution for the orbital eccentricity of each system. ^cRecovered α and β shape parameters of the standard Beta distribution from hierarchical Bayesian modeling.

Download table as:  ASCIITypeset image

However, it is also clear from this exercise that randomness in the draws from the underlying eccentricity distribution can occasionally produce results that formally agree with the true parameters at the 2–3σ level, but which have the potential to be mis- or overinterpreted. This is especially important for small samples, where stochasticity in the draws have a higher probability of producing results that qualitatively differ from the true underlying distribution. For example, in the first experiment with N = 5 and σ_e = 0.2, we find values of α = ${15.9}_{-12.5}^{+7.9}$ and β = ${31.8}_{-8.1}^{+18.2}$. The resulting Beta distribution peaks at higher eccentricities than the actual distribution, which would overestimate how dynamically hot the warm Jupiter exoplanet population really is.

Although the actual values of α and β are within 3σ of the joint constraints for all 12 experiments, we conclude that the results for small sample sizes (N ≲ 10) should be interpreted as an indication of the qualitative behavior of the population. Larger sample sizes are needed to more precisely uncover the detailed shape and quantitative constraints of the underlying distribution. Even these rely on the assumption that our input model f(e) is correct, however. It is certainly possible that another functional form better emulates the true population-level eccentricity distribution, but we avoid this type of model comparison in our orbit study because of the limited sample and the broad constraints on the eccentricity of each individual system. More targets and longer-term orbit monitoring will enable this type of model selection in the future.

To assess how readily different distributions can be distinguished from each other, we performed the same experiment for a uniform distribution as well as the warm Jupiter distribution mirrored at high eccentricities (Beta shape parameters α = 3.03 and β = 0.867). Because we are searching for population-level differences between giant planets and brown dwarfs, these results better reflect the goals of this study. Results are summarized in Figure 17. It is clear from these tests that larger samples and better measurement precision more reliably reproduce the underlying distribution, but differences between distributions can readily be inferred from small samples. For example, even in the most pessimistic case (N = 5, σ_e = 0.2), the warm Jupiter and high-e samples are noticeably distinct, although in this case the uniform distribution happens to resemble the warm Jupiter sample. We conclude that it is easier to determine whether two very different parent distributions are distinct from each other than it is to establish the exact shape of the underlying distribution.

Zoom In Zoom Out Reset image size

4.3.2. Eccentricity Distribution for the Full Sample

Our results for the underlying eccentricity distribution of the full sample of 27 substellar companions spanning 5–100 au and 2–75 M_Jup are shown in Figure 18 and summarized in Table 7. The best-fitting values of α and β are ${0.95}_{-0.43}^{+0.41}$ and ${1.30}_{-0.46}^{+0.61}$, respectively, with positive covariance between the two parameters. This corresponds to an approximately flat distribution across the entire range of eccentricities. Uncertainties in the posterior are larger at the lowest and highest eccentricities and allow for some flexibility at both ends of the distribution, especially near e = 0.

Zoom In Zoom Out Reset image size

Table 7.  Population-level Eccentricity Distributions

	Separation	Mass	Other	Sample	α	β
Sample	Range	Range	Constraints	Size	(1σ C.I.)	(1σ C.I.)
Full Sample	5–100 au	2 MJup ≤ M2 < 75 MJup	⋯	27	0.950 (0.520–1.36)	1.30 (0.840–1.91)
Giant Planets	5–100 au	2 MJup ≤ M2 < 15 MJup	⋯	9	30.0 (58.2–156)	200 (542–1000)
Brown Dwarfs	5–100 au	15 MJup ≤ M2 < 75 MJup	⋯	18	2.30 (1.02–3.68)	1.65 (0.860–2.55)
High Mass Ratio	5–100 au	0.01 < M2/M1 < 0.2	⋯	17	1.85 (0.930–3.34)	1.25 (0.800–2.44)
Low Mass Ratio	5–100 au	0.001 < M2/M1 < 0.01	⋯	10	1.50 (0.500–3.68)	4.50 (1.52–12.8)
Close Separations	5–30 au	2 MJup ≤ M2 < 75 MJup	⋯	14	1.70 (0.830–2.87)	2.75 (1.52–4.45)
Wide Separations	30–100 au	2 MJup ≤ M2 < 75 MJup	⋯	13	0.650 (0.260–1.28)	0.850 (0.510–1.56)
Planets Excluding HR 8799	5–100 au	2 MJup ≤ M2 < 15 MJup	No HR 8799 bcde	5	90.0 (86.9–311)	300 (555–1000)
Only HR 8799	5–100 au	2 MJup ≤ M2 < 15 MJup	Only HR 8799bcde	4	120 (43.8–127)	960 (570–1000)
Short-Period Brown Dwarfs	5–30 au	15 MJup ≤ M2 < 75 MJup	⋯	9	3.25 (1.18–8.05)	3.25 (1.37–7.22)
Long-Period Brown Dwarfs	30–100 au	15 MJup ≤ M2 < 75 MJup	⋯	9	3.12 (1.11–7.69)	1.75 (0.680–3.39)
Young Brown Dwarfs	5–100 au	15 MJup ≤ M2 < 75 MJup	<1 Gyr	11	1.50 (0.210–5.61)	1.05 (0.410–2.29)
Old Brown Dwarfs	5–100 au	15 MJup ≤ M2 < 75 MJup	>1 Gyr	7	4.70 (2.37–9.60)	4.35 (2.33–7.87)
Nearby Brown Dwarfs	5–100 au	15 MJup ≤ M2 < 75 MJup	<40 pc	8	4.50 (1.69–9.33)	2.25 (0.910–4.22)
Distant Brown Dwarfs	5–100 au	15 MJup ≤ M2 < 75 MJup	>40 pc	10	1.20 (0.500–7.06)	1.20 (0.680–6.91)
Small IWA	5–100 au	2 MJup ≤ M2 < 75 MJup	IWA/a < 0.5	12	1.00 (0.650–2.37)	1.20 (0.880–2.97)
Large IWA	5–100 au	2 MJup ≤ M2 < 75 MJup	IWA/a > 0.5	15	0.600 (0.430–1.46)	0.800 (0.770–2.25)

Download table as:  ASCIITypeset image

Four systems have eccentricities that are peaked near 1.0 with substantial power spanning all values: HD 49197 B, HD 1160 B, HR 2562 B, and 1RXS0342+1216 B. The origin of this shape is unclear, although it may be a result of an incorrect stellar mass or perhaps underestimated astrometric uncertainties. To test whether these objects bias the results, we ran the same analysis except with uniform eccentricity distributions for these systems. The results for the full sample are nearly indistinguishable from the original case. This also holds true for all of the additional experiments we carry out in the sections below.

One of the practical implications of this result is that this flat distribution can be reliably used as a Bayesian prior for orbit fits of new substellar companions that are discovered in the future. A uniform eccentricity distribution has generally been adopted for orbit fits in the past as an uninformative prior, but this can now justifiably be used as an informed prior: in the absence of discovery bias (see Section 4.4), a newly identified substellar companion is effectively equally likely to have a low, moderate, or high eccentricity. A somewhat more precise prior would make use of the actual best-fitting Beta distribution we identify, which differs slightly from a uniform distribution.

The clearest implication of a flat posterior is that regardless of the formation or migration processes that produce this population of low-mass companions, they do not appear to imprint a strong preference for high, intermediate, or low eccentricities, at least when marginalized over other parameters such as stellar host mass, substellar mass, and system age. This analysis for the full population also implicitly assumes that the underlying eccentricity distribution of this population does not vary as a function of companion mass or separation. That is, there is no strong gradient in the orbital properties of the sample within our sample spanning 5–100 au and 2–75 M_Jup. Below we test these assumptions by subdividing this full sample based on companion mass, mass ratio, orbital separation, and system age to determine whether there are signs of population-level changes in the eccentricities of these companions.

4.3.3. Giant Planets and Brown Dwarf Companions

We further subdivide our full sample of 27 substellar companions by mass to assess whether there is evidence for a difference between the population-level eccentricities of giant planets and brown dwarf companions. We adopt a dividing line of 15 M_Jup in this study, which is motivated by the approximate transition at smaller separations between the planet mass function and that of stellar or substellar companions from radial velocity samples (e.g., Udry & Santos 2007; Schneider et al. 2011). These subsamples consist of 9 directly imaged long-period giant planets (2–15 M_Jup) and 18 brown dwarf companions (15–75 M_Jup).

Fits to the giant planet subsample are presented in Figure 19 and Table 7. The best-fitting values are α = 30 and β = 200, with 1σ credible intervals spanning 58–156 and 540–1000, respectively. Note that the peak of this joint distribution does not fall within the marginalized 1σ intervals because of the strong positive covariance between these parameters. Nevertheless, despite this broad range of values, the two components in each {α, β} pair tend to balance each other to produce a fairly narrow range for the resulting eccentricity distribution function between e = 0.05–0.25, and with a peak value at $\bar{e}$ = 0.13.

Zoom In Zoom Out Reset image size

However, as discussed in Section 4.3.1, we caution that interpretations of the resulting eccentricity posterior for small samples should be made with care. Here our sample size comprises nine systems, but several of the individual eccentricities are so broad that they provide little meaningful constraint at the population level. Moreover, the additional discovery of a single moderate- to high-eccentricity planet has the potential to substantially inflate this distribution. So while the resulting posterior is expected to capture the overall qualitative trend of the underlying eccentricities, our experiments in Section 4.3.1 indicate that small sample sizes can influence the detailed behavior of the reconstructed distribution, especially near the endpoints at e = 0 or e = 1. It is therefore unclear whether the lack of power at the lowest eccentricities reflects something genuine about this population or is perhaps just a reflection of small number statistics. New discoveries and continued orbit monitoring of these systems are needed to address this question in the future.

Results for the sample of brown dwarf companions are summarized in Figure 20. The best-fitting hyperparameters are α = ${2.3}_{-1.3}^{+1.4}$ and β = ${1.7}_{-0.8}^{+0.9}$, with significant positive covariance between the two parameters. The eccentricity distribution for brown dwarfs peaks between e ≈ 0.6–0.9 with a broad range between e ≈ 0.3–1.0. If a single Beta distribution adequately describes this population, it appears that near-circular orbits are rare, and there is an indication that the highest values above 0.9 are disfavored.

Zoom In Zoom Out Reset image size

The most striking aspect of this underlying distribution is its dissimilarity with results for the giant planets. As a population, widely separated brown dwarf companions are significantly more eccentric and span a wider range of eccentricities than their lower mass counterparts. This tendency is evident in Figure 21, which shows eccentricities for individual systems as a function of companion mass. There is a noticeable dearth of planets at high eccentricities; most of the cumulative power is focused at low values. On the other hand, brown dwarfs span a wide range of eccentricities, with fewer companions on near-circular orbits. We revisit the implications of these distinct distributions in Section 5.

Zoom In Zoom Out Reset image size

To quantify the difference between these two inferred eccentricity distributions, we calculate the probability that random variables p drawn from the brown dwarf distribution (BD), are greater than the giant planet (GP) distribution (e.g., Cook 2003; Raineri et al. 2014),

$$\begin{eqnarray}&&\begin{array}{l}P({p}_{\mathrm{BD}}\gt {p}_{\mathrm{GP}})\\ \quad \ =\ {\displaystyle \int }_{0}^{1}{f}_{\mathrm{BD}}(e| {\alpha }_{\mathrm{BD}},{\beta }_{\mathrm{BD}}){F}_{\mathrm{GP}}(e| {\alpha }_{\mathrm{GP}},{\beta }_{\mathrm{GP}}){de},\end{array}\end{eqnarray} \tag{ 5 }$$

where ${f}_{\mathrm{BD}}(e| {\alpha }_{\mathrm{BD}},{\beta }_{\mathrm{BD}})$ is the probability distribution function for brown dwarfs and ${F}_{\mathrm{GP}}(e| {\alpha }_{\mathrm{GP}},{\beta }_{\mathrm{GP}})$ is the cumulative distribution function (CDF) for giant planets. The parameterized model we adopt for the underlying population-level eccentricities is the Beta distribution. The CDF of the Beta distribution is the regularized incomplete beta function,

$$\begin{eqnarray}&&{F}_{\mathrm{GP}}(e| {\alpha }_{\mathrm{GP}},{\beta }_{\mathrm{GP}})=\displaystyle \frac{{\int }_{0}^{e}{t}^{{\alpha }_{\mathrm{GP}}-1}{\left(1-t\right)}^{{\beta }_{\mathrm{GP}}-1}{dt}}{{\int }_{0}^{1}{t}^{{\alpha }_{\mathrm{GP}}-1}{\left(1-t\right)}^{{\beta }_{\mathrm{GP}}-1}{dt}}.\end{eqnarray} \tag{ 6 }$$

Pairs of α and β are randomly drawn from the joint posteriors for both the giant planet and brown dwarf MCMC results, and Equation (5) is numerically integrated for each trial. This procedure is repeated 10⁵ times to produce a distribution of probabilities. We find that P(p_BD > p_GP) = 0.979 with a 2σ credible interval of 0.85–1.0. Two equivalent distributions will produce probabilities of 50%, so this value of ≈98% points to a significant difference between these two populations.

4.3.4. A Mass Ratio Threshold

4.3.5. A Separation Threshold

We also examine whether there are differences in the eccentricity distribution as a function of separation, which might be expected if, for example, the inner population of substellar companions predominantly originates within a disk while the outer population largely represents the product of cloud fragmentation. For this experiment, we adopt a threshold of 30 au, which is chosen so as to divide the full sample of substellar companions into two approximately equal bins.

Results for the subsample of 14 companions between 5 and 30 au and 13 companions between 30 and 100 au are shown in Figure 22 and are summarized in Table 7. We find values of α = ${1.70}_{-0.9}^{+1.2}$ and β = ${2.8}_{-1.2}^{+1.7}$ for the sample at close separations, which corresponds to a broad peak between e = 0.1–0.4 with significant power from e = 0.0–0.8. At wide separations, we find α = ${0.7}_{-0.4}^{+0.6}$ and β = ${0.9}_{-0.3}^{+0.7}$. This corresponds to a roughly flat distribution with somewhat higher power at the bounded endpoints. There is some evidence that the more closely separated population lacks companions at the highest eccentricities, but these two distributions are otherwise broadly similar and are not nearly as distinct as the giant planet and brown dwarf subsamples. The probability that wide companions have higher eccentricities than close companions is 52.3% (2σ credible interval: 24.4%–81.7%).

Zoom In Zoom Out Reset image size

4.3.6. Exploring the Influence of HR 8799

4.3.7. Short- and Long-period Brown Dwarfs

4.3.8. Young and Old Brown Dwarfs

A direct imaging survey will preferentially find close-in companions with higher eccentricities compared to companions that have the same semimajor axes but are on circular orbits (e.g., Dupuy & Liu 2011; Kane 2013). This is because more eccentric orbits will reach wider apastron distances and therefore companions will spend more time at large separations compared to those with lower eccentricities. This preference produces a bias that can skew the apparent eccentricity distribution of discoveries toward high eccentricities. Discovery bias is strongest when the semimajor axis is much smaller than the IWA (at the same contrast as the companion), and it asymptotically disappears for companions with semimajor axes well beyond the IWA. This metric—IWA/a—is therefore a useful tool to assess whether a given sample is skewed to higher eccentricities from discovery bias.

Our sample of substellar companions draws from an assortment of AO imaging surveys carried out over the past two decades that had a wide range of sensitivities, inner working angles, and resulting contrast curves. Moreover, targets in our sample span a broad range of distances from 5.8 pc (Gl 229) to 126 pc (HD 1160); for a given contrast curve, more distant systems may be more heavily biased in favor of higher eccentricities for the same reasons shorter period companions are preferentially selected against.

If our sample is strongly influenced by discovery bias, then we would expect the eccentricities of our targets to exhibit several trends that are characteristic of this preferential selection. Below we discuss five such correlations in detail to assess the potential impact of this bias on our results.

• 1.  
  The IWA at the time of discovery should be comparable to or larger than the semimajor axes of companions in our sample. High values of IWA/a would suggest that a strong bias is likely present, while low values below unity would indicate minimal bias. The IWA here should be at the same contrast as the companion and must be determined at the time the companion was first imaged. We therefore revisited the original discovery papers for each system in our sample and visually estimated the IWA from the discovery images. These are then converted into physical units using the distance to the system and divided by the maximum a posteriori of the semimajor axis posterior from Table 3 (or values from the literature for systems that we did not refit in this study). The distribution of IWA/a for our sample is shown in Figure 23, and individual values are listed in Table 4. Thirteen systems from our full sample have IWA/a < 0.5, 12 systems have 0.5 ≤ IWA/a < 1.0, and 2 systems (HD 1160 B and HD 49197 B) have IWA/a ≥ 1.0. Dupuy & Liu (2011) simulated the impact of discovery bias as a function of both IWA/a and eccentricity (see their Figure 1) and found that there is minimal impact on the relative fraction of detected systems across all eccentricities for values of IWA/a ≲ 0.5. For IWA/a values between ≈0.5–1, there is a modest preferential suppression of low-eccentricity orbits. For values above IWA/a ≈ 1.0, this suppression becomes severe and only the highest eccentricities are detectable. We expect that about half of our sample suffers from minimal discovery bias and about half is moderately influenced by this effect.
• 2.  
  There should be a correlation between IWA/a and eccentricity. Objects with semimajor axes well outside the IWA should be minimally biased, whereas those at or inside the IWA should preferentially be more eccentric.The lower panel of Figure 23 shows the eccentricity of each system as a function of IWA/a. There is a broad range of eccentricities for IWA/a values between 0.0–1.0. The two systems with the highest values of IWA/a (HD 1160 B and HD 49197 B) have poorly constrained eccentricities, but these distributions favor higher values of e. This suggests that these two companions are probably influenced by this bias, but there is no obvious trend for the rest of the sample. This is reflected in the reconstructed eccentricity distributions for systems with IWA/a < 0.5 and IWA/a > 0.5, which both have approximately flat shapes spanning e = 0–1.
• 3.  
  More distant systems should have higher eccentricities, on average, as IWA/a becomes larger and this bias becomes stronger. A positive correlation between the population-level eccentricity distribution and the distance to the system would be another indicator that discovery bias may be important. To test whether our results for brown dwarf companions in Section 4 may be skewed to high values because of discovery bias, we reassess the population-level eccentricities for nearby (<40 pc) and distant (>40 pc) subsamples. Results from this exercise are shown in Figure 24 and Table 7. The recovered eccentricity distribution for nearby brown dwarfs peaks at high values ($\bar{e}$ ≈ 0.7). The best-fitting distribution for distant systems is flat across all eccentricities, with an overdensity of posterior values at modest eccentricities near e = 0.5. This apparent tendency for closer systems to have higher eccentricities is opposite to the trend we would expect if discovery bias played a strong role in shaping the population-level eccentricity distribution. With only eight objects in the nearby sample and nine in the distant sample, we expect that the slight differences we observe are caused by small number statistics.
• 4.  
  Lower mass companions with higher contrasts relative to their host stars should be more strongly biased toward high eccentricities compared to higher mass (lower contrast) companions. This is because contrast curves are not constant, but typically curve to larger angular separations at higher contrasts. That is, the effective IWA increases at lower companion masses. In Section 4 we showed that the inferred eccentricity distribution of giant planets in our sample is more circular than the distribution for brown dwarf companions. This is opposite of what we would expect if discovery bias played a significant role in shaping these observed distributions.
• 5.  
  Companions with shorter orbital periods should be more eccentric than their counterparts at wider separations. In Section 4 we derived the underlying eccentricity distributions for substellar companions at separations between 5–30 au and 30–100 au. The close-separation subsample exhibited a broad peak at $\bar{e}\approx 0.3$ with a paucity of power at the highest eccentricities. The wide-separation subsample was essentially flat across all eccentricities. Once again, these trends are opposite of what we would expect if discovery bias played a dominant role in shaping these distributions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Altogether, these series of tests argue against discovery bias playing a major role in shaping the population-level eccentricity distribution of substellar companions in our sample. It is possible (and likely) that some bias is present in this sample, but the eccentricity distributions we inferred in Section 4 appear to predominantly reflect the intrinsic properties of substellar companions. We conclude that our primary finding that brown dwarf companions are more eccentric than giant planets is robust against discovery bias.