CONSTRAINTS ON THE ARCHITECTURE OF THE HD 95086 PLANETARY SYSTEM WITH THE GEMINI PLANET IMAGER

Unlike for the majority of extrasolar planets detected through indirect techniques, the orbital periods of planets detected through direct imaging range from decades (e.g., Lafrenière et al. 2009; Chauvin et al. 2012) to thousands of years (e.g., Bailey et al. 2014; Naud et al. 2014). So far, only five directly imaged systems have had the orbital parameters of their planets constrained (e.g., Beust et al. 2014; De Rosa et al. 2015; Millar-Blanchaer et al. 2015; Sallum et al. 2015; Zurlo et al. 2016). In the case of β Pictoris, the astrometric measurements cover almost half of a complete orbit, while for the others the astrometric measurements only cover a small fraction of the total orbit. Orbital constraints derived from small astrometric arcs can still be valuable. In the case of Fomalhaut, revealing the high eccentricity of planet b (Kalas et al. 2013; Beust et al. 2014) led to the need for additional planet(s) to explain the dynamical history of the system (Faramaz et al. 2015). Dynamical simulations of the HR 8799 system might set upper limits on the masses of the four planets, above which the system would be dynamically unstable at the age of the star (e.g., Fabrycky & Murray-Clay 2010; Goździewski & Migaszewski 2014).

Orbital configuration of planetary systems can also offer insight into the interactions between planets and circumstellar material. In the case of HR 8799, a three-component debris disk is required to model its spectral energy distribution (SED) and resolved images with a dust-free gap (Su et al. 2009; Matthews et al. 2013). The semimajor axis of the four planets are both consistent with the location and geometry of this gap, evidence suggesting that the imaged planets are responsible for gravitationally sculpting the disk. This dust distribution is analogous to that seen in the HD 95086 system. HD 95086 (1.7 M_⊙, 17 ± 4 Myr, 90.4 ± 3.4 pc) has a large infrared excess (Rizzuto et al. 2011; Chen et al. 2012) modeled by either two- or three-component disks (Moór et al. 2013; Su et al. 2015) and marginally resolved with Herschel (Moór et al. 2013). At a projected separation of 56 au, intermediate to the two populations of dust, a planetary companion, HD 95086 b, with a model-dependent mass of 4.4 ± 0.8 M_Jup (De Rosa et al. 2016), was discovered by direct imaging (Rameau et al. 2013b).

HD 95086 b was imaged during the commissioning of the Gemini Planet Imager (GPI; Macintosh et al. 2014) in late 2013 (Galicher et al. 2014). In this Letter, we present astrometric measurements of HD 95086 b obtained with GPI at Gemini South Observatory, Chile, between 2014 and 2016, in which significant orbital motion is measured. By combining our new astrometric measurements and a reanalysis of those presented in Galicher et al. (2014), with previous measurements from Rameau et al. (2013a), we present the first constraints on the orbital parameters of HD 95086 b and discuss implications of the system architecture and properties of any additional planet.

Observations of HD 95086 b were made with GPI at four different epochs between 2013 December and 2016 March with either the H (1.5–1.8 μm) or K₁ (1.9–2.2 μm) filters (program IDs GS-ENG-GPI-COM, GS-2015A-Q-501, GS-2015B-Q-500, and GS-2016A-Q-12). HD 95086 was also observed on 2014 March 24 and 26 (program ID GS-ENG-GPI-COM) but the image quality was not sufficient to detect the planet. The observing sequences consisted of an observation of an argon arc-lamp to calibrate the offset between the position of the microspectra at the target elevation and at zenith caused by instrument flexure (Wolff et al. 2014), followed by between 30 and 100 minutes of on-source integration, centered at meridian passage to maximize the field rotation for angular differential imaging (Marois et al. 2006). A log of the observations is given in Table 1.

Table 1.  Observations and Orbital Parameters of HD 95086 b

UT Date	Instrument/	Tint	Field	Plate Scale	Position Angle	ρ	θ	References
	Filter	(minutes)	Rot. (deg)	(mas px−1)	Offset (deg)	(mas)	(deg)
2012 Jan 12	NaCo/L '	52.0	24.5	27.1 ± 0.06	0.38 ± 0.02	624 ± 8	151.8 ± 0.8	(1)
2013 Mar 14	NaCo/L '	54.0	31.4	27.1 ± 0.04	0.45 ± 0.09	626 ± 13	150.7 ± 1.3	(1)
2013 Jun 27	NaCo/L '	151.7	29.3	27.1 ± 0.2	0.48 ± 0.1	600 ± 11	150.9 ± 1.2	(1)
2013 Dec 10	GPI/K1	33.4	11.7	14.166 ± 0.007	−0.1 ± 0.13	619 ± 5	150.9 ± 0.5	(2)a
2013 Dec 11	GPI/H	41.6	15.0	14.166 ± 0.007	−0.1 ± 0.13	618 ± 11	150.3 ± 1.1	(2)a
2014 May 13	GPI/K1	34.0	16.9	14.166 ± 0.007	−0.1 ± 0.13	618 ± 8	150.2 ± 0.7	(2)
2015 Apr 06	GPI/K1	41.7	37.7	14.166 ± 0.007	−0.1 ± 0.13	622 ± 7	148.8 ± 0.6	(2)
2015 Apr 08	GPI/K1	99.4	39.0	14.166 ± 0.007	−0.1 ± 0.13	622 ± 4	149.0 ± 0.4	(2)
2016 Feb 29	GPI/H	38.0	16.9	14.166 ± 0.007	−0.1 ± 0.13	621 ± 5	147.8 ± 0.5	(2)
2016 Mar 06	GPI/H	78.0	47.6	14.166 ± 0.007	−0.1 ± 0.13	620 ± 5	147.2 ± 0.5	(2)
Parameter			Unit	${{ \mathcal L }}_{{\rm{max}}}$	Median	68% CI	95% CI	i prior
Semimajor axis (a)			au	50.5	61.7	53.3–82.4	44.8–137	$\mathrm{cos}(i)$
				49.3	59.7	52.7–72.0	44.5–105.5	155 ± 5°
Eccentricity (e)			⋯	0.019	0.135	<0.205	<0.443	$\mathrm{cos}(i)$
				0.168	0.124	<0.187	<0.397	155 ± 5°
Inclination (i)			deg	146.0	153.0	139.5–162.7	126.1–174.7	$\mathrm{cos}(i)$
				149.2	155.7	150.9–160.6	146.2–165.5	155 ± 5°
Argument of periastron (ω)			deg	39.7	92.5	27.3–153.7	3.73–176.2	$\mathrm{cos}(i)$
				99.2	92.6	28.1–153.0	3.92–176.1	155 ± 5°
Position angle of nodes (Ω)			deg	329.7	83.8	38.2–146.7	5.48–174.7	$\mathrm{cos}(i)$
				2.5	89.8	32.4–149.4	4.76–175.3	155 ± 5°
Epoch of periastron (T0)			⋯	2172.37	2155.16	2040.26–2367.95	2005.35–2809.37	$\mathrm{cos}(i)$
				2202.44	2151.83	2038.27–2337.07	2005.12-2580.66	155 ± 5°
Period (P)			year	257.1	370.9	296.8–572.6	227.9–1230.9	$\mathrm{cos}(i)$
				260.0	352.8	291.2–469.2	225.5–834.1	155 ± 5°

Note.

^aThese data were published in Galicher et al. (2014) but re-analyzed in this work.

References. (1) Rameau et al. (2013a), (2) this work.

Download table as:  ASCIITypeset image

The GPI observations on 2013 December were first published by Galicher et al. (2014) but were reprocessed in this work for a uniform set of measurements, owing to an improvement of the spot registration and astrometric calibration. The raw GPI images were reduced with the Data Reduction Pipeline (DRP; v1.3.0²³ ) to produce (x, y, λ) calibrated data cubes. The DRP was also used to measure the satellite spot positions for precise image registration (Wang et al. 2014). The reduced images were filtered with an unsharp mask (15 × 15 pixels), the speckles were then subtracted using LOCI (Lafrenière et al. 2007), with dr = 5 pixels, N_A = 300 times the full-width-at-half-maximum (FWHM), g = 1, and N_δ = 0.75 or 1 FWHM depending on the amount of field rotation. The residual images were finally rotated so that north was aligned with the vertical axis, combined with a trimmed mean (10%), and collapsed to build a broadband final image.

The astrometry of HD 95086 b was extracted at all epochs following Marois et al. (2010) and Lagrange et al. (2010). A model point-spread function was built from the spots and injected into the data at an initial position. The ameba-simplex optimization method (Nelder & Mead 1965) was used to iterate over the position and flux to minimize the residuals in a wedge of Δr × Δθ = 2 × 2 FWHM. Errors were calculated by injecting fake planets at the same separations but 20 position angles uniformly distributed between 90° and 270° away from planet b and repeating the same exercise. Final uncertainties on the relative companion position were computed by adding quadratically the errors from the measurements, the spot registration, and the plate scale and position angle offset (orientation of north with respect to the vertical axis). For the GPI measurements we used the same astrometric solution as described in De Rosa et al. (2015) to convert the image position (x, y) into an on-sky separation and position angle. The results of this analysis are presented in Table 1. Measurements from the literature were also compiled in this study (Rameau et al. 2013a, 2013b). They were obtained with VLT/NaCo at L ' (3.5–4.1 μm) at three different epochs between 2012 January and 2013 June.

The astrometry of HD 95086 b given in Table 1 showed a roughly constant separation and a decreasing position angle with time, consistent with face-on, circular orbital motion (see Figures 1 and 2). The orbital motion was still clear within the GPI data sets that covered more than two years. As in De Rosa et al. (2015), we fitted the orbit using a rejection sampling technique that is more computationally efficient than traditional Markov Chain Monte Carlo (MCMC) techniques at generating posterior parameter distribution for short orbital arcs. Briefly, this technique sampled the posterior probability density function of orbital parameters by generating orbits from prior distributions, scaling the semimajor axis, and rotating the position angle of nodes to fit one of the observational epochs, rotating through each epoch. Final orbits were rejected or accepted by comparing the probability ($P\propto {e}^{-{\chi }^{2}/2}$) of the remaining astrometric epochs to a uniform random variable. This method, introduced in Blunt et al. (2016), will be described more in depth in S. C. Blunt et al. (2016, in preparation). To fit HD 95086 b, we used priors that are uniform in $\mathrm{log}a$, $\mathrm{cos}i$, ω, Ω, and T₀, and a linearly decreasing prior in e fit to radial velocity planets (Nielsen et al. 2008). Based on the inclination of the outer disk of i = 25 ± 5° (without knowledge of the rotation direction, both 25° and 180 − 25 = 155° are possible; Moór et al. 2013; Su et al. 2015), we also ran a fit to the astrometry with a restricted prior on the inclination. A pair of Gaussians centered on i = 25° (counter-clockwise motion) and 180 − 25 = 155° (clockwise motion), both with σ = 5°, was used as the prior on the inclination.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 plots the posterior distributions of the orbital parameters of HD 95086 b. Overall, the posteriors distributions correspond to face-on, circular orbits. Table 1 gives, for the two inclination priors, the most likely (${{ \mathcal L }}_{{\rm{max}}}$), median, and 68% and 95% confidence intervals (CIs) for each orbital parameter. For the inclination-restricted fit, the posterior distributions are not significantly changed, although the median of the inclination distribution was shifted to slightly higher values due to the shape of the prior. Our finding of a low eccentricity was relatively independent of our choice of the prior on either the inclination and eccentricity, currently e < 0.44 at 95% confidence using the linear prior in e, but switching to a uniform prior still resulted in e < 0.56 at 95% confidence. Fitting the data points without the outlier from 2013 June 27 did not affect our conclusions about the semimajor axis, eccentricity, and inclination of the orbits.

Zoom In Zoom Out Reset image size

We confirmed the reliability of the rejection sampling fits with the same Metropolis–Hastings MCMC orbit fitting used in Nielsen et al. (2014), using the same priors on parameters as above. As in De Rosa et al. (2015) we found excellent agreement between the two methods, as shown in the diagonal elements of Figure 3.

In addition to these Monte Carlo techniques, the method for constraining orbital parameters over short orbital arcs presented in Pearce et al. (2015) was applied to the astrometry of HD 95086 b. The angle between the projected separation and velocity vectors was calculated as $\varphi ={96.7}_{-9.2}^{+9.4}\;$ deg and a value of the dimensionless parameter B of ${0.52}_{-0.16}^{+0.20}$. Comparing these to the minimum inclination and eccentricity contours of Pearce et al. (2015), the orbital parameters of HD 95086 b can only be constrained to $e\geqslant {0.15}_{-0.15}^{+0.31}$ and $i\leqslant {59.1}_{-13.8}^{+9.6}$ deg (corresponding to $i\geqslant {120.9}_{-9.6}^{+13.8}$ deg). In each case, uncertainties on the measured positional offset of HD 95086 b were propagated in a Monte Carlo fashion. While these limits are consistent with the values in Table 1, they are significantly less constraining.

The point-source sensitivity of the GPI data was estimated by measuring the noise in concentric annuli in the residual LOCI images. The throughput was computed by injecting fake planets in the raw data that were reduced with the same coefficients as the science images. The most sensitive GPI observations were obtained on 2015 April 8 at K₁ band (see Figure 1), and the deepest NaCo L ' detection limit was taken from Rameau et al. (2013a). The planet-to-star contrast was converted into predicted mass with the AMES-COND (Baraffe et al. 2003) model. An optimized version of the Monte-Carlo-based MESS tool (Bonavita et al. 2012) was used to generated random on-sky positions of planets in order to compute detection probabilities over a separation range of 1–1000 au, with a 2 au step size, and a planet mass range of 0.5–20 M_Jup, with a 0.5 M_Jup step size. The distributions of the orbital parameters were the same as for the orbit fitting. For each point in the mass–semimajor axis grid, 10,000 orbits were randomly generated, and the fraction of planets that would have been detected in either the GPI or NaCo observations was used as the completeness at that point. The final completeness map built is shown in Figure 4 (left panel).

Zoom In Zoom Out Reset image size

The architecture of the HD 95086 system can now be constrained based on the first estimates of the orbital parameters of planet b and on the detection limits reached with the current ensemble of observations. Su et al. (2015) suggested that the belts that produce the dust properties inferred from the SED and the Herschel images are separated by a dust-free gap from ∼8 to ∼80 au. They proposed four non-exhaustive scenarios to explain the large size of the gap. While there are currently large uncertainties in the exact location of the dust rings, and a large number of alternative multiple planet configurations which may explain the disk gap, the current observational constraints can be used to explore the four specific scenarios presented in Su et al. (2015). By assuming a system in which the disk and planet(s) are coplanar, the orbital parameters of HD 95086 b, and the associated completeness for additional companion (Figure 4, right panel), allows us to:

• 1.  
  Rule out scenario A in which the planet b would be responsible for clearing the entire gap with an eccentricity of ∼0.7. Based on the orbit fit presented in Table 1, an eccentricity of e > 0.40 for planet b can be excluded at the 95% confidence level, and therefore planet b is unlikely to alone account for the gap. Moreover, Su et al. (2015) argued that the outer disk might have e < 0.25. This morphology is consistent with the probable orbit of the planet.
• 2.  
  Neither validate nor eliminate scenario B in which planet b resides on a moderately eccentric orbit (∼0.3), within the allowed range of orbital parameters presented in Table 1 and with a semimajor axis of ∼40 au at the edge of the 95% CI. Another inner planet with moderate eccentricity at ∼16 au is required but is not within the scope of our observations. Other configurations proposed by Su et al. (2015) with a more massive inner planet on a higher eccentric orbit and planet b with reduced eccentricity are unconstrained by our observations.
• 3.  
  Reconsider scenario C in which two additional 7 M_Jup planets are needed at ∼12 and ∼26 au with low eccentricities. An eccentricity of ∼0.1 and a semimajor axis of ∼56 au for planet b as suggested by this model are permitted by the orbital fitting. However, the second planet in this scenario at 26 au with a mass of 7 M_Jup and e < 0.3 would have been detected in our observations. A conservative 75% threshold was reached at this separation at 4 M_Jup, or at 25 au with the expected mass, thus very close to the predicted properties of the second planet. The inner planet remained undetectable.
• 4.  
  Rule out scenario D in which three equal-mass (5 M_Jup) planets on circular orbits at ∼11, ∼19, and ∼34 au populate the gap in addition to planet b. The third 5 M_Jup planet would have been detected at 100% confidence at the predicted position with a null eccentricity.

Therefore, in the context of a ∼8–80 au gap with coplanar planets with masses derived from the AMES-COND models, the scenario with two planets (B) cannot be ruled out, whereas the scenario with three 7 M_Jup planets (C) is inconsistent with our observations; instead, the second planet would require a higher eccentricity with either a lower mass or smaller semimajor axis in order to be consistent with our observational constraints. The tempting four equal-mass (5 M_Jup) planet scenario (D), analogous to the HR 8799 system, was rejected based on the non-detection of the third planet. All the scenarios might, however, still be accepted if the predicted masses are wrong by a factor of a few, assuming other evolutionary models (e.g., Spiegel & Burrows 2012; Marleau & Cumming 2014).

As noted previously, the architecture of the debris belts was inferred from the SED and marginally resolved image data of the system and provided dynamical constraints on the location of planets. However, the geometry of a debris disk based solely on SED modeling can be wrong by a factor of a few (Booth et al. 2013; Morales et al. 2013; Pawellek et al. 2014). Grain properties and dust position are degenerate. Resolved images of the different disk components are necessary to break this degeneracy (e.g., Lebreton et al. 2012). This behavior was mentioned in Su et al. (2015), in which the warm and cold belts might be located between 7 and 8 au, and between 60 and 125 au, respectively. The outer edge of the inner belt might also be poorly determined if the emission is dominated by its innermost region, artificially making a more extended inner disk into a narrow and closer-in ring. These uncertainties affect the predicted semimajor axes, eccentricities, and masses of the planets in the aforementioned scenarios. Even if a conservative 75% detection threshold was applied to assess the different hypotheses, a change by a factor of a few in mass or semimajor axis might result in a drastic change in detection probability at the predicted properties of the planets.

We have presented four epochs of astrometric measurements of the young giant planet HD 95086 b with the GPI. We were able to detect a significant change in the position angle of the planet between the end of 2013 and early 2016, whereas the separation remained similar. By using Monte Carlo methods to sample probable orbits, we placed the first constraints on the orbital parameters of HD 95086 b. The semimajor axis distribution peaks near its median value of 61.7 au, while the eccentricity and inclination were limited to e < 0.21 and $i=153\fdg {0}_{-13.5}^{+9.7}\;$, respectively, at 68% confidence.

We also discussed the system architecture in the context of a coplanar planet(s)–debris disk configuration. Based on the preliminary constraints on the orbit of planet b, and on the detection probabilities, we tested the scenarios from Su et al. (2015) to explain the sustainability of the ∼8–80 au dust-free gap inferred from SED fitting and resolved submillimeter imaging. Of the four scenarios, only the two-planet scenario (B) cannot be excluded based on the observational constraints presented here. The upper limit on the eccentricity of planet b confidently ruled out a single-planet scenario (A), and the GPI sensitivity excluded the three and four equal-mass planet scenarios (C and D) by a non-detection of the second and third planet on low eccentricity orbits at 26 and 34 au, respectively. The four scenarios presented in Su et al. (2015) are not exhaustive, and modifications to masses and orbits of additional companions are still possible to sustain the gap, while also being dynamically stable for 17 Myr.

Due to the degeneracy between the properties of additional planets, extensive exploration of the parameter space is necessary to further assess the architecture of the system beyond the analytical statements made above. Such an analysis will also require the precise geometry of the debris. Resolved images of the outer disk with ALMA, or by polarimetry with GPI or SPHERE (Beuzit et al. 2008), combined with the detection of the inner disk by interferometry will not only constrain the locations of the inner and outer edges of the gap, but also the eccentricities of the warm and cold belts.

We thank Kate Y. L. Su for interesting discussions about the debris system and Virginie Faramaz for her fruitful advices on handling chaotic zones. 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, Tecnología e Innovación Productiva (Argentina). J.R., R.D., and D.L. acknowledge support from the Fonds de Recherche du Québec. Supported by NSF grants AST-1518332 (R.J.D.R., J.R.G., J.J.W., T.M.E., P.K.), AST-1411868 (B.M., A.R., K.W.D.), and AST-141378 (G.D.). Supported by NASA grants NNX15AD95G/NEXSS and NNX15AC89G (R.J.D.R., J.R.G., P.K., J.J.W., T.M.E.), and NNX14AJ80G (E.L.N., S.C.B., B.M., F.M., M.P.). 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 (S.M.A.).

Facility: Gemini:South (GPI) - Gemini South Telescope.