ASTROMETRIC CONFIRMATION AND PRELIMINARY ORBITAL PARAMETERS OF THE YOUNG EXOPLANET 51 ERIDANI b WITH THE GEMINI PLANET IMAGER

Monitoring the orbital motion of exoplanets, through direct imaging (e.g., Chauvin et al. 2012; Kalas et al. 2013; Bonnefoy et al. 2014; Nielsen et al. 2014), or through indirect techniques such as radial velocity and transit measurements (e.g., Cumming et al. 2008; Howard et al. 2012; Marcy et al. 2014; Moutou et al. 2015), can provide a wealth of information about their properties, the processes through which they form, and how they interact dynamically with other bodies in the system. Accurately determining the orbital parameters of exoplanets can constrain their masses and densities (e.g., Charbonneau et al. 2000) and lead either to the discovery of additional planets in the system (e.g., Nesvorný et al. 2012), or to the exclusion of additional planets within a range of periods by invoking dynamical stability arguments (e.g., Correia et al. 2005). Orbital parameters also provide insight as to how planetary companions dynamically interact with circumstellar material (e.g., Millar-Blanchaer et al. 2015). While the orbital periods are typically decades or longer for directly imaged planets, accurate astrometric monitoring of these systems can lead to preliminary constraints on their orbital parameters before a significant portion of the orbit is observed (e.g., Pueyo et al. 2015).

51 Eridani (51 Eri) is a nearby (29.43 ± 0.29 pc; van Leeuwen 2007) member of the young (24 ± 3 Myr; Bell et al. 2015) β Pictoris moving group (Zuckerman et al. 2001). Recently, Macintosh et al. (2015) reported the discovery of a low-mass (2–10 M_Jup) planet at a projected separation of 13.2 ± 0.2 AU based on observations with the Gemini Planet Imager (GPI; Macintosh et al. 2014). Based on J- and H-band spectroscopy, and L' photometry, 51 Eri b was shown to have a spectrum with strong methane and water absorption, with a temperature of 600–750 K (Macintosh et al. 2015). Due to the short baseline between discovery and follow-up, it was only possible to rule out a stationary background object from a 2003 non-detection; closer brown dwarf interlopers with non-zero proper motions could not be excluded. A statistical argument based on the space density of T-dwarfs (Reylé et al. 2010; Burningham et al. 2013), and the allowed range of distances of a foreground or background object based on the apparent magnitude of 51 Eri b combined with the luminosity of T-dwarfs (Kirkpatrick et al. 2012), was used to limit the possibility of an unassociated field brown dwarf to a probability of 2.4 × 10⁻⁶.

In addition to the newly resolved planetary companion, 51 Eri has an infrared excess indicative of a circumstellar debris disk (Patel et al. 2014; Riviere-Marichalar et al. 2014). The debris disk has yet to be spatially resolved, so its geometry is unconstrained. At a projected separation of 1960 AU lies GJ 3305—an M-dwarf binary with a semimajor axis of 9.80 ± 0.15 AU (Montet et al. 2015)—which is co-moving with 51 Eri, forming a bound hierarchical system (Feigelson et al. 2006). While the inclination of the GJ 3305 binary is well constrained (i = 921 ± 02; Montet et al. 2015), the period of the wide 51 Eri–GJ 3305 binary (∼10⁴ years; Feigelson et al. 2006) precludes any such estimation of its inclination.

51 Eri was initially observed with GPI at Gemini South as a part of the GPI Exoplanet Survey (GPIES) on 2014 December 18 UT (GS-2014B-Q-500). A faint companion candidate was identified, and subsequent observations demonstrated that its spectral energy distribution was consistent with that of a low-temperature, low-surface gravity giant planet (Macintosh et al. 2015). In total, 51 Eri b has been successfully observed four times since discovery in 2014, three times with GPI (GS-2014B-Q-501, GS-2015A-Q-501), and once with NIRC2 at the W. M. Keck 2 telescope using the facility adaptive optics system (Wizinowich et al. 2000). A summary of these observations is given in Table 1. For each GPI epoch, the observing strategy was the same. The target was acquired before transit to maximize field rotation (Marois et al. 2006), and the observations were taken using the spectral coronagraphic mode, with either the J- or H-band filters. In addition to these successful observations, 51 Eri was also observed with GPI in J-band on 2015 January 29, where 51 Eri b was not recovered due to poor image quality, and in H-band on 2015 August 30 and 31 where, although 51 Eri b was recovered, the signal-to-noise ratio was significantly worse.

Table 1.  Observations and Orbital Parameters of 51 Eridani b

UT Date	MJD	Instrument	Filter	tint	Ncoadd	Nexp	Plate Scale	Position Angle	ρ	θ
				(s)			(mas px−1)a	Offset (deg)	(mas)	(deg)
2014 Dec 18	57009.1292	GPI	H	59.6	1	38	14.166 ± 0.007	−1.10 ± 0.13	449.8 ± 6.7	171.0 ± 0.9
2015 Jan 30	57052.0572	GPI	J	59.6	1	70	14.166 ± 0.007	−1.10 ± 0.13	454.0 ± 6.4	170.6 ± 1.0
2015 Jan 31	57052.9753	GPI	H	59.6	1	64	14.166 ± 0.007	−1.10 ± 0.13	461.8 ± 7.1	170.5 ± 0.9
2015 Feb 01	57054.0364	NIRC2	L'	0.9	60	62	9.952 ± 0.002b	0.252 ± 0.009b	461.5 ± 23.9	170.4 ± 3.0
2015 Sep 01	57266.4052	GPI	H	59.6	1	93	14.166 ± 0.007	−1.10 ± 0.13	454.7 ± 5.7	166.5 ± 0.6
Preliminary Orbital Parameters of 51 Eridani b
Parameter		Unit	${\chi }_{{\rm{min}}}^{2}$		Median		Lower	Upper
Semimajor axis (a)		AU	18.61		14		11	21
Eccentricity (e)		...	0.2804		0.21		0.06	0.40
Inclination (i)		deg	126.7		138		125	153
Argument of periastron (ω)		deg	4.233		90		32	148
Position angle of nodes (Ω)		deg	166.6		75		47	100
Epoch of periastron (T0)		...	2016.05		2014.79		2003.31	2026.11
Period (P)		year	60.68		41		29	76

Notes.

^aYelda et al. (2010). ^bIn reduced GPI data cubes, one pixel is equivalent to one lenslet.

Download table as:  ASCIITypeset image

The GPI observations obtained on 2015 September 1 were reduced using the GPI Data Reduction Pipeline (DRP; Perrin et al. 2014).²⁸ The dark current was subtracted and bad pixels were identified and fixed. The shift in the position of the microspectra on the detector due to mechanical flexure was measured by comparing reference argon arcs taken monthly, to arcs obtained after target acquisition (Wolff et al. 2014). The microspectra were extracted, converting the 2D image into a 3D (x, y, λ) data cube. These were then divided by a flat field to correct for lenslet throughput, and were interpolated along the wavelength axis to a common wavelength vector across the bandpass. Finally, the optical distortion was corrected for using measurements obtained with a pinhole mask (Konopacky et al. 2014).

To minimize potential biases between the astrometry presented in Macintosh et al. (2015) and the new measurements presented here, we used Pipelines 1 and 3 from Macintosh et al. (2015) to both perform the point-spread function (PSF) subtraction, and to extract the astrometry of 51 Eri b. Pipeline 2 was switched to a Python implementation (Wang et al. 2015) of the Karhunen–Loève Image Projection algorithm (Soummer et al. 2012), and uses a forward-modeled PSF to perform a Markov Chain Monte Carlo (MCMC) analysis to determine the posterior distributions of the separation and position angle of 51 Eri b.

The plate scale and position angle of GPI have been monitored by continually observing a set of astrometric calibrators with well-determined orbital solutions or contemporaneous NIRC2 measurements, which has an accurate astrometric solution (Yelda et al. 2010). In addition to observations listed in Konopacky et al. (2014), we have observed the θ¹ Ori B quadruple system an additional four times, the HD 157516 binary twice, and the HIP 80628 binary once. These observations were reduced as above, and PSF subtraction was not required. The pixel positions of each component were measured as in Konopacky et al. (2014).

Combining these measurements results in a plate scale of 14.166 ± 0.007 mas lenslet⁻¹, and a position angle offset of −110 ± 013. The position angle offset is defined as the angle between the lenslet y-axis, and the elevation axis of the telescope, measured east from north. Since version 1.2, the DRP takes into account an offset of −1° while processing the data, and as such the difference between the true position angle and the measured position angle in a DRP-reduced image is ${\theta }_{{\rm{true}}}-{\theta }_{{\rm{measured}}}=-0\buildrel{\circ}\over{.} 10\pm 0\buildrel{\circ}\over{.} 13$. There is no evidence of variations of either the plate scale or position angle offset between observing runs within the measurement uncertainties, so a single value is adopted for all epochs. The revised astrometric calibration was used to recalculate the astrometry from Macintosh et al. (2015), which are consistent with the previous values, and are shown in Table 1.

Typically, confirmation of common proper motion is achieved by comparing the motion of a candidate with respect to a background track for a stationary background object with negligible parallax, so that the only relevant movement is the parallax and proper motion of the primary star (e.g., Nielsen et al. 2013). Such an analysis for 51 Eri b is shown in Figure 1, with the astrometry diverging from the stationary background object track (2.8-σ in ρ, 8.8-σ in θ). The spectrum of 51 Eri b excluded the distant background star hypothesis, leaving two plausible fits to the spectrum: an unbound field brown dwarf, or a bound planet, as described by Macintosh et al. (2015), who found a probability of 2.4 × 10⁻⁶ that the object was an unbound field brown dwarf. By including the latest epoch, this probability can be further reduced by solving for the allowable parallax and proper motion of an unbound object, and reducing the volume in which an unbound T-dwarf could exist.

Zoom In Zoom Out Reset image size

Using the measured parallax and proper motion of 51 Eri (van Leeuwen 2007) the relative astrometry of 51 Eri b was converted into absolute astrometry. A Metropolis–Hastings MCMC technique (e.g., Ford 2006) was then used to fit the proper motion and parallax of 51 Eri b to both the absolute astrometry, and the non-detection from 2003 discussed in Macintosh et al. (2015). A uniform prior in the two proper motion directions was used, and a p(d) ∝ d² prior for the distance, with a maximum distance of 200 pc. Good convergence was achieved in the MCMC chains, with a Gelman–Rubin statistic <1.00007. The posteriors from the MCMC fit are shown in Figure 2. The median distance from the posteriors is 31 pc, with a 68% confidence interval between 25 and 39 pc, consistent with the distance to 51 Eri (29.4 ± 0.3 pc). The offset in proper motion with respect to the star is consistent with orbital motion occurring over the span of the observations.

Zoom In Zoom Out Reset image size

The calculation of the unbound brown dwarf probability from Macintosh et al. (2015) is updated using this new distance constraint. The previous calculation found the product of the number density of T-dwarfs and the volume of a cone with the angular width of the GPI detector out to the largest distance a T-dwarf could be seen in the reduced GPI image of a given target, repeated for each of the 44 stars observed as a part of GPIES at the time of the discovery. There was a further correction factor based on the proper motion constraints from the 2003 non-detection, as only 34% of the expected distribution of field brown dwarfs had the motion required to place the object behind the star in 2003. This correction is not used here as the revised proper motion estimate for the unbound brown dwarf scenario is 3σ discrepant from the proper motion of the star due to orbital motion (Figure 2). Instead, the volume calculation of the cones was convolved with the distance posterior distribution, lowering the posterior probability that 51 Eri b is an unbound brown dwarf to 2 × 10⁻⁷.

MCMC orbit fitting is slow to converge for sparsely sampled astrometry or short orbital arcs, and so we implemented a more computationally efficient Monte Carlo method to generate plausible distributions of orbital parameters based on astrometry covering only a small fraction of an orbital period. For four of the orbital parameters, a large number of samples were drawn from appropriate probability distributions: uniform for argument of periastron (ω) and epoch of periastron passage (T₀), uniform in cos(i) for inclination angle (i), and eccentricity (e) following the linear fit to radial velocity planets of Nielsen et al. (2008). Initial orbits were generated using these four parameter distributions and fixed values of semimajor axis (a) and position angle of nodes (Ω), the values of a was scaled and the value of Ω rotated to reproduce the astrometry at the first epoch. Period (P) was not fit and was calculated assuming a stellar mass of 1.75 M_⊙ (Simon & Schaefer 2011). Astrometric errors were incorporated into the generated parameters by adding random offsets to the separation and position angle of the first epoch before each orbit is shifted and rotated. These uncertainties were randomly drawn from Gaussian distributions equivalent to the first epoch astrometric uncertainties.

With multiple epochs, we proceeded to iteratively reject sets of orbital parameters that do not match later measurements, with acceptance probability given by a two-dimensional Gaussian with the astrometric errors as the standard deviations. Generated orbits that were closest to the observed separation and position angle at the corresponding observational epoch were more likely to be accepted. We then obtained distributions of fitting orbital parameters given the input astrometry. Varying the epoch chosen to initialize the procedure had little effect on the distribution of accepted orbital parameters.

We validated this method by generating 10³ orbits with one fixed orbital parameter and randomly sampling the other parameters, and then creating five epochs from one year of simulated astrometry for each orbit, each with an observational uncertainty. The spacing of the epochs and astrometric errors were chosen to be the same as the measurements of 51 Eri b. We then applied our method to the artificial astrometry and examined the returned distributions of orbital parameters. Of the simulated orbits with semimajor axis of 13 AU, the median was 13.3 AU, and in 74% of orbits 13 AU was within the 68% confidence interval of each individual orbital fit for semimajor axis. This suggests the generated distributions are reasonable representations of the posterior probability distributions, and comparison with MCMC (described below) further supports this. This procedure will be described further in S. C. Blunt et al. (2015, in preparation).

This technique was applied to the astrometry of 51 Eri b. One hundred orbits from the fit are shown in Figure 3, and the distributions of accepted orbital parameters are shown in Figure 4. While the marginal distribution for inclination using the discovery epoch follows the shape of the prior, the distribution after incorporating data from 2015 peaks at ∼138°, showing the new observations have provided significant orbital constraints. More moderate changes are seen in eccentricity and semimajor axis, with less eccentric orbits, and orbits closer to 14 AU, becoming more favored with increasing number of measurements. The orbital elements corresponding to the minimum χ², representing the best-fitting orbit generated, are given in Table 1, along with the median and limits corresponding to where 68% of the Monte Carlo generated orbits were found. In order to compute the median and 68% intervals, ω was wrapped to be within 0°–180°, Ω to within 30°–120°, and T₀ to within 1995–1995+P.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The shapes of the distributions of accepted orbital parameters were confirmed using an affine invariant MCMC ensemble sampler (Foreman-Mackey et al. 2013)²⁹ , with the same prior on each parameter. The shapes of the posterior distributions are similar (Figure 4); the median and 1σ confidence intervals for a and i were found to be ${14.1}_{-3.0}^{+8.2}$ AU and ${135}_{-13}^{+15}$ deg, consistent with the values in Table 1, but taking two orders of magnitude longer to compute. In addition to these Monte Carlo techniques, the astrometry of 51 Eri b was applied to the technique for constraining orbital parameters over short orbital arcs presented in Pearce et al. (2015). The angle between the projected separation and velocity vectors was calculated as $\phi ={84.5}_{-14.3}^{+14.2}$ deg, and a value for their dimensionless parameter B of ${0.23}_{-0.09}^{+0.12}.$ Comparing these to the minimum inclination and eccentricity contours of Pearce et al. (2015), the eccentricity of 51 Eri b is unconstrained, and the inclination is restricted to $i\lt {77}_{-7}^{+5}$ deg (or $i\gt {103}_{-5}^{+7}$ deg). While consistent with the values in Table 1, these limits are significantly less constraining.

We have presented a new astrometric measurement of the young exoplanet 51 Eridani b, which provides further evidence supporting the bound companion hypothesis. Using this new measurement, we observe significant deviation (2.8 σ in ρ, 8.8 σ in θ) from the stationary background object track, with a reduced ${\chi }_{{\rm{bkg}}}^{2}$ of 11.19 (Nielsen et al. 2012). We compute a probability that it is an interloping field brown dwarf of 2 × 10⁻⁷, decreased by a factor of twelve relative to Macintosh et al. (2015). By implementing a computationally efficient Monte Carlo method to sample probable orbits, we place the first constraints on the orbital parameters of 51 Eri b. Based on the present astrometry, the median of the semimajor axis distribution is 14 AU, corresponding to a period of 41 years, assuming a mass of 1.75 M_⊙ for 51 Eri. While the additional astrometric epoch did not significantly change the distribution of accepted semimajor axes and periods, the range of allowed inclinations was significantly constrained, with a median value of 138°. The eccentricity of the orbit remains unconstrained, with circular orbits only marginally preferred relative to the prior distribution.

Based on these preliminary constraints, the orbit of 51 Eri b does not appear to be co-planar with the orbit of the M-dwarf binary GJ 3305 (i = 921 ± 02; Montet et al. 2015). The inclination of the orbital plane of the outer binary relative to 51 Eri is not known. The large separation between 51 Eri and GJ 3305, and the young age of the system, would suggest that secular Lidov–Kozai oscillations would not have had sufficient time to significantly alter the semimajor axis of the planet (Montet et al. 2015), although moderate changes in the inclination and eccentricity are not excluded (Fabrycky & Tremaine 2007). Given the long timescale of these oscillations (200 Myr for a perturber on a circular orbit, Montet et al. 2015), this effect could be completely suppressed by quicker secular precession due to, e.g., a relatively low-mass, undetected planet (Wu & Murray 2003). Alternatively, the planet may remain on its current orbit, despite the presence of the wide binary (Holman & Wiegert 1999). Continued astrometric monitoring of 51 Eri b over the next few years should be sufficient to detect curvature in the orbit, further constraining the semimajor axis and inclination of the orbit, and placing the first constraints on the eccentricity. Absolute astrometric measurements of 51 Eri with GAIA (e.g., Perryman et al. 2014), in conjunction with monitoring of the relative astrometry of 51 Eri b, will enable a direct measurement of the mass of the planet. Combined with the well-constrained age of 51 Eri b, such a determination would provide insight into the evolutionary history of low-mass directly imaged extrasolar planets, and help distinguish between a hot-start or core accretion formation process for this planet.

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 National Science Foundation (NSF) on behalf of the Gemini partnership: the NSF (United States), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), Ministério da Ciência, Tecnologia e Inovação (Brazil) and Ministerio de Ciencia, Tecnologi´a e Innovación Productiva (Argentina). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. Supported by NSF grants AST-0909188 and AST-1313718 (R.J.D.R., J.R.G., J.J.W., T.M.E., P.G.K.), AST-1411868 (B.M., K.F., J.L. P., A.R., K.W.D.), AST-141378 (P.A., G.D., M.P.F.), NNX11AF74G (A.Z.G., A.S.), and DGE-1232825 (A.Z.G.). Supported by NASA grants NNX15AD95G/NEXSS and NNX11AD21G (R.J.D.R., J.R.G., J.J.W., T.M.E., P.G.K.), and NNX14AJ80G (E.L.N., S.C.B., B.M., F.M., M.P.). J.R., R.D. and D.L. acknowledge support from the Fonds de Recherche du Québec. 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 (SMA). B.G. and M.J.G. acknowledge support from the National Sciences and Engineering Research Council of Canada. G.V. acknowledges a JPL Research and Technology Grant for improvements to the GPI CAL system.

Facility:  Gemini:South (GPI) - .