The following article is Free article

Dynamical Constraints on the HR 8799 Planets with GPI

Jason J. Wang, James R. Graham, Rebekah Dawson, Daniel Fabrycky, Robert J. De Rosa, Laurent Pueyo, Quinn Konopacky, Bruce Macintosh, Christian Marois, Eugene Chiang

Published 2018 October 12 • © 2018. The American Astronomical Society. All rights reserved.
The Astronomical Journal, Volume 156, Number 5Citation Jason J. Wang et al 2018 AJ 156 192DOI 10.3847/1538-3881/aae150

PDF Opens in a new tab.

ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

Article metrics

3968 Total downloads
0 Video abstract views

Permissions

Get permission to re-use this article

Dates

Received 2018 July 20
Revised 2018 September 5
Accepted 2018 September 11
Published 2018 October 12

Keywords

astrometry; planet–disk interactions; planets and satellites: dynamical evolution and stability; planets and satellites: gaseous planets; stars: individual (HR 8799); techniques: high angular resolution

Journal RSS

Create or edit your corridor alerts

What are corridors?

Abstract

The HR 8799 system uniquely harbors four young super-Jupiters whose orbits can provide insights into the system’s dynamical history and constrain the masses of the planets themselves. Using the Gemini Planet Imager, we obtained down to one milliarcsecond precision on the astrometry of these planets. We assessed four-planet orbit models with different levels of constraints and found that assuming the planets are near 1:2:4:8 period commensurabilities, or are coplanar, does not worsen the fit. We added the prior that the planets must have been stable for the age of the system (40 Myr) by running orbit configurations from our posteriors through N-body simulations and varying the masses of the planets. We found that only assuming the planets are both coplanar and near 1:2:4:8 period commensurabilities produces dynamically stable orbits in large quantities. Our posterior of stable coplanar orbits tightly constrains the planets’ orbits, and we discuss implications for the outermost planet b shaping the debris disk. A four-planet resonance lock is not necessary for stability up to now. However, planet pairs d and e, and c and d, are each likely locked in two-body resonances for stability if their component masses are above 6 M_Jup and 7 M_Jup, respectively. Combining the dynamical and luminosity constraints on the masses using hot-start evolutionary models and a system age of 42 ± 5 Myr, we found the mass of planet b to be 5.8 ± 0.5 M_Jup, and the masses of planets c, d, and e to be ${7.2}_{-0.7}^{+0.6}\,{M}_{\mathrm{Jup}}$ each.

Export citation and abstractBibTeX RIS

Previous article in issue

Next article in issue

1. Introduction

High-contrast imaging spatially separates the faint light of planets from the bright glare of their host star. By monitoring exoplanetary systems with high-contrast imaging, we are able to obtain footage of these exoplanets in motion and trace out their orbits. Orbit analysis has been a powerful tool in characterizing the dynamics of directly imaged systems. Through orbital monitoring of β Pic b, we now know that the planet is responsible for inducing the observed warp in the circumstellar debris disk (Dawson et al. 2011; Lagrange et al. 2012), although it may not be alone in clearing out the cavity of the disk (Millar-Blanchaer et al. 2015). Precise orbital determination also has timed the Hill sphere transit of the planet to between 2017 April and 2018 January (Wang et al. 2016), which offered a unique opportunity to probe the circumplanetary environment of a young exoplanet (de Mooij et al. 2017; Mékarnia et al. 2017; Stuik et al. 2017). For HD 95086 b, by combining orbit fits with constraints on the debris disk geometry, Rameau et al. (2016) showed that the planet alone cannot be clearing out the gap in the system, and that additional planets reside closer in to the star. The orbit of Fomalhaut b was shown to cross the debris disk in the system, revealing that the planet cannot be a massive Jupiter-like planet, but rather a dwarf planet shrouded by dust (Kalas et al. 2013). Finally, future orbital monitoring of 51 Eri b could shed light on the interactions between the planet and the wide-separation binary GJ 3305 (De Rosa et al. 2015).

Long-term orbital monitoring can also lead to dynamical mass measurements of the planets themselves, which will assess evolutionary models of young giant planets that all current mass estimates of directly imaged exoplanets are based on (Baraffe et al. 2003; Marley et al. 2007). In the coming years, Gaia will measure the astrometric reflex motion of stars hosting planets (Perryman et al. 2014). Gaia astrometry combined with long-term orbital monitoring from direct imaging will provide the tightest model-independent constraints on the masses of the planets (Sozzetti et al. 2016). Alternatively, multi-planet systems where planets mutually perturb their orbits provide another method of constraining the masses of the planets in the system. In resonant systems where the dynamical timescales are close to the orbital timescales, such mutual perturbations have been measured in short period planets as variations in the host star’s radial velocity signature (e.g., Marcy et al. 2001; Rivera et al. 2010) and as transit timing variations (e.g., Agol et al. 2005; Holman & Murray 2005), leading to direct measurements of the masses. Due to the long orbital periods of known directly imaged systems, such a direct measurement of the mutual perturbations on the orbits has been impossible with the current observational baselines, none of which span a full orbital period. Still, upper limits on the masses of the planets based on dynamical stability can be obtained. Stability mass constraints have been used to characterize exoplanets discovered in compact systems, such as TRAPPIST-1 (e.g., Quarles et al. 2017; Tamayo et al. 2017), Kepler-36 (Deck et al. 2012), and the HR 8799 system discussed in this paper.

HR 8799 is unique among directly imaged systems, as it is the only one known to harbor four planets (Marois et al. 2008, 2010). The planets orbit ∼15–70 au from the star between two rings of rocky bodies, similar to the configuration of the giant planets in our own solar system (Su et al. 2009; Reidemeister et al. 2009). The outer belt has been resolved with far-infrared and millimeter observations, although the exact orientation and inner edge of the disk are not entirely agreed upon (Hughes et al. 2011; Matthews et al. 2014; Booth et al. 2016; Wilner et al. 2018). Assuming “hot-start” evolutionary models and an age of 30 Myr, Marois et al. (2008, 2010) translated the planet luminosities into masses: planet b is ∼5 M_Jup and the inner three planets are ∼7 M_Jup (Marois et al. 2008, 2010). However, as the evolutionary models are uncertain at these early ages, so are the exact masses of the planets. Fortunately, dynamics can provide an additional constraint on the masses of the planets, even if their long orbital periods mean we cannot detect planet–planet interactions and fully constrain the masses this way.

Since the discovery of the HR 8799 planets, their orbits have been closely monitored. Keplerian orbits have been fit to the astrometry obtained from many instruments using least-squares techniques that look for families of orbits or Bayesian parameter estimation with Markov-chain Monte Carlo (MCMC) methods that explore the full posterior of orbital parameters (Soummer et al. 2011; Currie et al. 2012; Esposito et al. 2013; Maire et al. 2015; Pueyo et al. 2015; Konopacky et al. 2016; Zurlo et al. 2016; Wertz et al. 2017). Fitting the planets independently, some studies have reported planet d to be misaligned in its orbit relative to the other planets (Currie et al. 2012; Esposito et al. 2013; Pueyo et al. 2015) or one of the inner planets to have eccentricity above 0.2 (Maire et al. 2015; Wertz et al. 2017). However, several of the authors have noted that unaccounted astrometric calibration offsets between instruments may be inducing inclination and eccentricity biases (Maire et al. 2015; Pueyo et al. 2015; Konopacky et al. 2016). Recently, Konopacky et al. (2016) presented self-consistent astrometry using only measurements from Keck and found that coplanar and low-eccentricity solutions were consistent with the data. Despite the uniform analysis, the 7 years of Keck data still only cover a short arc of these orbits, which have periods between ∼40–400 years, leaving many possible orbital configurations.

The measured astrometry is not the only constraint on the orbit of these planets. HR 8799 is part of the Columba moving group (Zuckerman et al. 2011), a group of stars that formed together ${42}_{-4}^{+6}$ Myr ago (Bell et al. 2015). Thus the four planets need to be stable dynamically for almost the same amount of time. Studies using N-body simulations have explored the dynamical constraints on the orbital parameters and masses. These studies have found stable orbits using the nominal luminosity-derived masses from Marois et al. (2010) without invoking orbital resonances (Sudol & Haghighipour 2012; Götberg et al. 2016), or to even higher masses assuming long-term resonance lock of the planets (Fabrycky & Murray-Clay 2010; Marois et al. 2010; Currie et al. 2011; Goździewski & Migaszewski 2014; Gozdziewski & Migaszewski 2018). However, many of these studies initialize or fit the simulated orbits to one astrometric measurement, leaving a gap between orbit fits from the data and dynamical constraints from simulations (Fabrycky & Murray-Clay 2010; Marois et al. 2010; Sudol & Haghighipour 2012; Götberg et al. 2016). To connect simulations to the data more rigorously, Goździewski & Migaszewski (2014) developed a novel technique to lock the planets into resonance and then search for times and orientations that matched all of the available data. However, their orbit and mass constraints only apply to the family of orbits that slowly migrated into a four-planet resonance lock.

A few attempts have been made to include stability in the orbit fitting of this system. Analytical prescriptions have been used to remove the orbits that are most obviously not dynamically stable (Pueyo et al. 2015; Konopacky et al. 2016). Esposito et al. (2013) ran N-body simulations on their orbital fits from a least-squares algorithm and only found stable orbits up to 5 M_Jup. In general, finding stable orbits in the orbit fits has been impractical with short orbital arcs. Having only the 2D sky projection of an arc of an orbit, even with milliarcsecond-level precision, cannot break many degeneracies in the orbital parameters, resulting in too wide a variety of orbital solutions which are nearly all unstable.

In this paper, we present an analysis that better bridges the gap between orbit fits and dynamical constraints by incorporating N-body simulations as a rejection sampling step of our Bayesian orbit fit to enforce stability. In Section 2, we show we have obtained 1 milliarcsecond astrometry of all four planets using the Gemini Planet Imager (GPI; Macintosh et al. 2014) and the open-source pyKLIP data reduction package (Wang et al. 2015). In Section 3, we combine the precise GPI measurements with the uniformly reduced Keck astrometry measured by Konopacky et al. (2016) and fit multiple orbital models with different assumptions about coplanarity and resonance using MCMC techniques that sample the full posterior of possible orbital configurations. In Section 4, we take the posteriors of orbits from our Bayesian analysis and simulate them for 40 Myr using the REBOUND N-body integrator (Rein & Liu 2012) to find the posterior of stable orbits after applying a dynamical stability prior. We discuss the consequences of our results, such as planets shaping the cold debris disk, the necessity of orbital resonances for stability, dynamical limits on the masses of the planets, and the future stability of the system.

2. Observations and Data Reduction

To obtain astrometry of the planets, we used three epochs of observations of HR 8799 taken with the integral field spectroscopy (IFS) mode of GPI. Two epochs were from instrument commissioning (Gemini program GS-ENG-GPI-COM) and one epoch from the GPI Exoplanet Survey (Gemini program GS-2015B-Q-500; PI: Macintosh). Details of the three observations are listed in Table 1. While HR 8799 b is normally located outside of the field of view of GPI, we steered the field of view on the detector during the 2014 September 12 observations to see planet b, although the conditions in this data set were too poor to see planet e.

Table 1. GPI Observations of HR 8799

UT Date	Filter	Time (s)	Rotation (°)	Imaged
		Exposure	Field	Planets
2013 Nov 17	K1	2130	17	cde
2014 Sep 12	H	3107	19	bcd
2016 Sep 19	H	3579	21	cde

Download table as: ASCII Typeset image

Raw IFS data from each epoch were processed to create 3D spectral datacubes using the automated data reduction system for the GPI Exoplanet Survey (Wang et al. 2018). Briefly, the data were dark subtracted, individual micro-spectra on the detector were extracted to form spectral datacubes, bad pixels were corrected, distortion in the image was corrected, and satellite spots, fiducial diffraction spots centered about the location of the star, were located. The star center in each wavelength channel is estimated using the satellite spots to correct any remaining differential atmospheric refraction not removed by the atmospheric dispersion corrector. See Appendix A of Wang et al. (2018) for details.

We used the Karhunen–Loève Image Projection algorithm (KLIP; Soummer et al. 2012; Pueyo et al. 2015) to subtract the stellar glare and the Bayesian KLIP-FM Astrometry (BKA) technique (Wang et al. 2016) to measure the astrometry of each planet. BKA forward models the distortions to the planet point-spread function (PSF) induced by KLIP in subtracting the stellar PSF and fits for the planet position while also accounting for the correlated noise in the image as a Gaussian process. In Wang et al. (2016), we used this technique to obtain 1 milliarcsecond astrometry on β Pic b. We used the KLIP and BKA implementations available in the pyKLIP package (Wang et al. 2015) from commit 4f56e34. For all the reductions, we first ran a high-pass filter to suppress the low spatial frequency background, constructed the instrumental PSF from the satellite spots, selected an annulus containing each planet to run KLIP on, and averaged the data in time and wavelength. To optimize the detection of each planet, we varied the number of Karhunen–Loève (KL) modes to model the stellar PSF, and the minimum number of pixels the planet needed to move in the reference images due to angular differential imaging (Marois et al. 2006a) and spectral differential imaging (Marois et al. 2000). We listed these parameters in Table 2. To measure the planets’ astrometry, we used the emcee package (Foreman-Mackey et al. 2013) to sample the posterior distribution for the location of the planet while also fitting the noise as a Gaussian process with spatial correlation described by the same Matérn covariance function as used in Wang et al. (2016). For each planet, our Markov-chain Monte Carlo sampler used 100 walkers, and each walker was run for 800 steps, with a “burn-in” of 300 steps beforehand that corresponded to at least three autocorrelation times for any chain. We then added additional terms in our astrometric error budget in quadrature: a 0.05 pixel uncertainty in locating the central star (Wang et al. 2014), a plate scale of 14.166 ± 0.007 mas lenslet⁻¹, and a residual North offset of 0 fdg 10 ± 0 fdg 13 (De Rosa et al. 2015).

Table 2. Astrometric Measurements of the HR 8799 Planets

UT Date	Planet	KL Modes	Criterion (pixels)	Separation (mas)	(°)
			Exclusion	Radial	Position Angle
2013 Nov 17	c	10	3	949.5 ± 0.9	325.18 ± 0.14
	d	10	3	654.6 ± 0.9	214.15 ± 0.15
	e	20	1.5	382.6 ± 2.1	265.13 ± 0.24
2014 Sep 12	b	10	1.5	1721.2 ± 1.4	65.46 ± 0.14
	c	10	1.5	949.0 ± 1.1	326.53 ± 0.14
	d	10	1.5	662.5 ± 1.3	216.57 ± 0.17
2016 Sep 19	c	10	2	944.2 ± 1.0	330.01 ± 0.14
	d	10	2	674.5 ± 1.0	221.81 ± 0.15
	e	10	1	384.8 ± 1.7	281.68 ± 0.25

Download table as: ASCII Typeset image

Our final astrometric results are listed in Table 2. We achieved down to 1 mas precision on the astrometry of planets b, c, and d. For these three planets, which are further from the star, the dominant sources of uncertainty are from the location of the star and the astrometric calibration of GPI. We achieved 1–2 mas precision on planet e, which is limited by the signal to noise ratio of the planet. This is 1.5–2 times more precise than the SPHERE astrometry from Wertz et al. (2017) and at least three times more precise than the Keck astrometry from Konopacky et al. (2016).

3. Orbit Fitting

To investigate the possible orbital solutions for the HR 8799 planets, we combined our GPI measurements with those from Keck that were reported in Konopacky et al. (2016). We chose to consider only these two data sets to minimize unknown systematic errors in the astrometric calibration across instruments. Specifically, GPI is astrometrically calibrated against the NIRC2 instrument at Keck, the same instrument used for the Keck HR 8799 observations, so systematic offsets between the two data sets are minimized (Konopacky et al. 2014; De Rosa et al. 2015). While Hubble Space Telescope data from 1998 provides an additional 6 years of baseline, the 20–30 mas 1σ uncertainties are not particularly constraining, so we did not use them (Soummer et al. 2011).

In this section, we fit the four planet orbits to four orbital configurations with increasing constraints: first, fitting four Keplerian orbits that share the same parallax and stellar mass to the data (Section 3.1); second, forcing coplanarity of the four planets (Section 3.2); third, forcing the four planets to be near 1:2:4:8 period commensurabilities but with no coplanarity constraints (Section 3.3); and lastly, forcing both coplanarity and the periods to be near a 1:2:4:8 ratio (Section 3.4). The constraints are intended to tighten the parameter space around stable orbits, but we are not directly considering stable orbits in these orbit fits. Dynamical stability constraints will be added in Section 4.

3.1. Unconstrained Orbits

First, we fit four independent Keplerian orbits to the data. We employed the same Bayesian framework as Wang et al. (2016) that used MCMC to sample the posterior distribution of orbital elements. For each planet, we fit for the conventional Keplerian orbital elements: semimajor axis (a), epoch of periastron after MJD 50,000 in units of fractional orbital period (τ), argument of periastron (ω), longitude of the ascending node (Ω), inclination (i), and eccentricity (e). Our conventions follow those defined in Alzner & Argyle (2012) for binary stars. In this approach, each planet’s orbital properties are independent, except we require that the four planets’ orbits use the same parallax and total system mass, which we take to be the stellar mass. To account for the uncertainties in the parallax and stellar mass, we assumed a Gaussian prior for the system parallax of 24.76 ± 0.64 mas (Gaia Collaboration et al. 2016) and a Gaussian prior for the stellar mass of 1.52 ± 0.15 M_⊙, which is mass reported by Baines et al. (2012) but with an additional 10% uncertainty to account for systematic model errors, as was done in Konopacky et al. (2016). This case covers the full range of orbital parameters that are consistent with the data; the three following orbit fits will explore subsets of this parameter space. Due to the high dimensionality of the orbital parameters (26 in total), it will be incredibly difficult to find the dynamically stable orbits if they reside in a very small subspace. Regardless, this orbital fit is an important fiducial case to be used as a baseline model with minimal assumptions. We will refer to this orbital fit as the “Unconstrained” fit.

We generally used uniform priors on our orbital parameters. For each planet, the prior on a was uniform in log(a) between 1 and 100 au; the prior on τ was uniform between 0 and 1; the priors on ω and Ω were uniform from 0 to 2π; the prior on i is the geometric $\sin (i)$ prior between 0 and π; and the prior on e was uniform between 0.000001 to 0.999. We note that our choice of orbital parameters will result in dual peaks in the ω and Ω posteriors that reflect our ignorance of the planets’ radial velocities.

We used the parallel-tempered affine-invariant sampler (Goodman & Weare 2010) implemented in emcee (Foreman-Mackey et al. 2013) using 15 temperatures and 1500 walkers per temperature. To improve the speed of convergence of the orbit fit, we initialized the walkers by drawing from allowed orbital parameters of individual fits to each planet using the same process. We ran each walker for 125,000 steps, after an initial burn-in of 95,000 steps. Convergence was assessed using the autocorrelation time and confirming by eye that ω and Ω had symmetric peaks. On a 32 core machine with AMD Opteron 6378 processors clocked at 2.3 GHz, this took 7 days to complete, although we note that we did not make an attempt to optimize the code. We then thinned the chains by a factor of 75 to mitigate any correlation in the Markov chains. Taking only the lowest temperature walkers, we then were left with 2,499,000 samples of the posterior distribution. The posterior distributions are plotted in Figure 1 and reported in Table 3.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Posteriors of each planet’s orbital parameters for each of the four different models considered in Section 3. Each row contains the four planet’s posteriors (color coded by planet) for one model. For the coplanar models, the planets have the same Ω and i, so only one is plotted.
Download figure:
Standard image High-resolution image

Table 3. Orbital Parameters of HR 8799 bcde from Different Models

Body	Element	Unconstrained	Coplanar	Near 1:2:4:8	Coplanar Low-e
	Orbital				Near 1:2:4:8
b	a_b (au)	${69.5}_{-7.0}^{+9.3}$	${66.4}_{-3.6}^{+4.1}$	${69.4}_{-4.0}^{+3.1}$	${69.5}_{-2.8}^{+2.6}$
	τ_b	${0.54}_{-0.16}^{+0.14}$	${0.46}_{-0.06}^{+0.05}$	${0.40}_{-0.15}^{+0.11}$	${0.38}_{-0.12}^{+0.11}$
	ω_b (°)	${92}_{-34}^{+30}$	92 ± 15	${95}_{-41}^{+49}$	${102}_{-46}^{+39}$
	Ω_b (°)	${127}_{-93}^{+32}$	${126}_{-14}^{+12}$	${82}_{-16}^{+36}$	${78}_{-10}^{+13}$
	i_b (°)	${29}_{-8}^{+7}$	23 ± 5	${23}_{-6}^{+4}$	24 ± 3
	e_b	0.15 ± 0.05	0.15 ± 0.06	${0.07}_{-0.05}^{+0.06}$	${0.05}_{-0.03}^{+0.04}$
c	a_c (au)	${37.6}_{-1.7}^{+2.2}$	${40.5}_{-1.7}^{+2.7}$	${41.2}_{-1.6}^{+2.3}$	${43.3}_{-1.7}^{+1.9}$
	τ_c	${0.50}_{-0.18}^{+0.10}$	${0.14}_{-0.11}^{+0.18}$	${0.13}_{-0.09}^{+0.15}$	${0.09}_{-0.07}^{+0.02}$
	ω_c (°)	${65}_{-29}^{+59}$	${52}_{-35}^{+83}$	${48}_{-29}^{+62}$	${63}_{-29}^{+60}$
	Ω_c (°)	${110}_{-47}^{+38}$	${126}_{-14}^{+12}$	${112}_{-26}^{+17}$	${78}_{-10}^{+13}$
	i_c (°)	${20}_{-5}^{+4}$	23 ± 5	${21}_{-4}^{+3}$	24 ± 3
	e_c	0.09 ± 0.04	${0.05}_{-0.03}^{+0.05}$	${0.04}_{-0.03}^{+0.05}$	${0.03}_{-0.02}^{+0.04}$
d	a_d (au)	${27.7}_{-1.7}^{+2.2}$	${25.3}_{-1.1}^{+1.3}$	${25.6}_{-1.3}^{+1.2}$	${25.6}_{-0.9}^{+1.0}$
	τ_d	${0.79}_{-0.18}^{+0.07}$	${0.872}_{-0.016}^{+0.019}$	0.85 ± 0.03	0.839 ± 0.20
	ω_d (°)	${144}_{-23}^{+13}$	${133}_{-11}^{+15}$	${148}_{-137}^{+22}$	${165}_{-157}^{+11}$
	Ω_d (°)	${92}_{-15}^{+27}$	${126}_{-14}^{+12}$	${86}_{-16}^{+26}$	${78}_{-10}^{+13}$
	i_d (°)	33 ± 4	23 ± 5	${23}_{-6}^{+5}$	24 ± 3
	e_d	0.15 ± 0.11	0.28 ± 0.04	0.20 ± 0.05	${0.18}_{-0.03}^{+0.02}$
e	a_e (au)	${15.3}_{-1.1}^{+1.4}$	${14.0}_{-0.6}^{+0.7}$	${15.7}_{-0.07}^{+0.06}$	15.4 ± 0.06
	τ_e	${0.71}_{-0.33}^{+0.17}$	${0.91}_{-0.06}^{+0.05}$	${0.10}_{-0.06}^{+0.15}$	${0.07}_{-0.04}^{+0.05}$
	ω_e (°)	${100}_{-49}^{+27}$	${128}_{-18}^{+25}$	${86}_{-36}^{+44}$	${76}_{-21}^{+25}$
	Ω_e (°)	117 ± 17	${126}_{-14}^{+12}$	${90}_{-19}^{+9}$	${78}_{-10}^{+13}$
	i_e (°)	31 ± 5	23 ± 5	${28}_{-4}^{+3}$	24 ± 3
	e_e	${0.13}_{-0.05}^{+0.06}$	0.16 ± 0.05	${0.05}_{-0.04}^{+0.07}$	${0.08}_{-0.04}^{+0.03}$
A	Parallax (mas)	24.70 ± 0.16	24.60 ± 0.56	24.38 ± 0.62	${24.38}_{-0.53}^{+0.54}$
	M_⋆ (M_⊙ )	${1.48}_{-0.04}^{+0.05}$	${1.46}_{-0.11}^{+0.12}$	${1.42}_{-0.11}^{+0.12}$	1.40 ± 0.10

	${\chi }_{\nu }^{2}$	${1.01}_{-0.07}^{+0.10}$	0.88 ± 0.07	${0.95}_{-0.06}^{+0.07}$	${0.94}_{-0.05}^{+0.12}$
	ΔBIC	${0}_{-6}^{+7}$	−34 ± 0.06	$-{4}_{-6}^{+5}$	$-{29}_{-4}^{+6}$

Stable orbits (first 10⁶ draws)		0	0	1	441

Note. The quoted values for ω and Ω are wrapped to be between 0° and 180°, so posterior percentiles describe one of the two symmetric peaks. For each parameter, the median value is reported with the superscript and subscript corresponding to the 84th and 16th percentiles of the distribution, respectively. For a normal distribution, these values correspond to the mean and 1σ range.

Download table as: ASCII Typeset image

Following similar analyses from previous orbit fitting studies (e.g., Konopacky et al. 2016; Wertz et al. 2017), we investigate the mutual inclination of the planets’ orbits by plotting in Figure 2 Ω and i, the two orbital elements that describe the orientation of the orbital plane. We will assume the planets orbit in the same direction. A planet with Ω differing by 180^◦ would be in a retrograde orbit relative to the other planets, which we do not consider here. We see that the 1σ contours for the four planets do overlap near i ∼ 30° and Ω ∼ 100°, indicating coplanar orbital solutions exist. This result agrees with the assessment of coplanarity by Konopacky et al. (2016) using similar arguments, although they preferred a different Ω. However, we note that only 0.005% of our sampled orbits have all four planets mutually inclined by <10°. This result likely indicates that without any constraints on the orientation of the orbital planes, it is extremely inefficient to sample coplanar orbits in large quantities. This is not surprising since the near-coplanar solutions are just a small subset of an eight-dimensional space in which we have chosen uniform, uncorrelated priors on each parameter. To more rigorously test coplanar orbits, we will fit directly for them (Section 3.2) and assess the fits (Section 3.5).

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Posteriors of each planet’s angles Ω and i for the Unconstrained fit. Blue, magenta, green, and yellow correspond to planets b, c, d, and e, respectively. The 1σ contour is plotted on top of each planet’s histogram. Overlapping regions indicate where coplanar orbits reside. Note that Ω is wrapped to only consider angles between 0^◦ and 180^◦, as the posterior is identical between 180^◦ and 360^◦ by construction.
Download figure:
Standard image High-resolution image

3.2. Coplanar Orbits

As planets form from the circumstellar disk, it would not be surprising to find the planets residing in coplanar orbits. The posteriors from the fit without constraints are consistent with coplanarity, but do not strongly favor it. Here we will explicitly fit for coplanar orbits, and in Section 3.5, we will assess if this approach fits the data as well as the unconstrained one. We modify our fit so that all four planets share the same values of Ω and i, reducing the fit to 20 orbital parameters. We will refer to this orbital fit as the “Coplanar” fit.

We used a parallel-tempered sampler with 15 temperatures and 1500 walkers per temperature. We ran each walker for 87,500 steps, after an initial burn-in of 132,500 steps. Convergence was assessed in the same way as in Section 3.1. We again thinned the chains by a factor of 75, and formed our posterior distribution from the lowest temperature chains. Our posterior distribution has 1,749,000 samples. The posteriors are plotted in Figure 1 and reported in Table 3.

From the posteriors in Figure 1, we see that the angles Ω and i that define the orientation of the orbital plane are consistent with the orbital planes of the four planets of the Unconstrained fit. The Coplanar orbits favor inclinations between 20° and 30°, which is ∼10° more face-on than the solutions from Konopacky et al. (2016) with just the Keck data alone. We still find Ω > 90°, which is not preferred for coplanar orbits in Konopacky et al. (2016). The solutions where Ω ≈ 90°, which are in agreement with Konopacky et al. (2016), however, favor lower inclinations near i ≈ 20°. While there are some differences on the preferred values, we note that many of these values are not ruled out by Konopacky et al. (2016) in their analysis.

We also find that forcing the system to be coplanar causes the eccentricity of planet d to be much higher, with <2% of the allowed orbits having e < 0.2. This was due to nearly all of planet d’s low-eccentricity orbits from the Unconstrained fit lying outside of the range of allowed orbital planes from the Coplanar fit. As Ω and i are constrained by the other three planets, raising e_d provided a way to obtain the best fits to the data. Note that the systems with Ω near 90° did have the lowest eccentricities for planet d.

3.3. Near 1:2:4:8 Period Ratio Orbits

We then investigated resonant orbits, focusing in particular on the 1:2:4:8 resonance, where consecutive pairs of planets are in 2:1 period resonance. We will first choose to be agnostic about the four planets’ mutual inclinations. Because these planets are not massless, even if they are in resonance, they do not necessarily reside at the exact period commensurabilities. Additionally, precession of the planets’ longitude of periastrons can further offset the observed period ratios from exact integer values. We note that previous orbit fitting work has assumed exact period commensurabilities when assessing if the fits were consistent with certain resonances.

At high planet masses like the HR 8799 planets, stable period ratios for the 2:1 two-body resonance tend to be larger than 2 due to instability caused by resonance overlaps at smaller period ratios (Morrison & Kratter 2016). Thus, instead of fixing the period ratio of the planet pairs, we use a parameter that gives each period ratio room to float. We picked our priors empirically from our own preliminary analysis of where the stable orbits existed. Our prior on the period ratio between b:c and c:d is a uniform distribution between 1.8 and 2.4. For the d:e period ratio, we chose a narrower uniform prior between 1.8 and 2.2, because we found all of the dynamically stable orbits were in this more narrow range and limiting it as such improved the efficiency of finding dynamically stable orbits (Section 4). We will show in Section 4.2 that our choices for our priors did not exclude stable orbits. Note that because we effectively replaced the parameters for the semimajor axes of the outer planets with their period ratios, we did not reduce the number of parameters in our MCMC fit, even though the parameter space shrank in size. We will refer to this orbital fit as the “Near 1:2:4:8” fit, which as the naming implies, only places the period ratios near resonance and does not guarantee the planets are indeed in resonance at all.

We initialized the walkers using coplanar solutions, which delayed convergence and caused the walkers to take a considerable amount of time to fully explore all of the allowed parameter space. We ran our parallel-tempered sampler with 15 temperatures and 1500 walkers per temperature for 75,000 steps, after a burn-in of 495,000 steps that was chosen using the same metric for convergence as Section 3.1. We performed the same thinning of the chains by a factor of 75. The resulting posterior was taken from the lowest temperature walkers and has 1,500,000 samples. The posteriors are plotted in Figure 1 and reported in Table 3.

3.4. Near 1:2:4:8 Period Ratio Coplanar Low-e Orbits

Lastly, we looked at coplanar resonant orbits. We applied both the coplanarity and period ratio constraints from Sections 3.2 and 3.3. We also applied an additional constraint that the eccentricity of all of the orbits have to be less than 0.2. In Section 3.2, we found that <2% of the coplanar orbits have e_d < 0.2. From preliminary analysis done concurrently with the orbit fits, we could only find stable orbits when all planets had e < 0.2. This fact will be further reinforced by the analysis in Section 4.2. Thus, we do not believe we lost stable orbits by applying this additional constraint, and merely improved the efficiency of finding stable orbits. We will refer to this orbital fit as the “Near 1:2:4:8 Coplanar” fit, which, like the Near 1:2:4:8 fit, does not guarantee the planets are actually in resonance.

For this 20-parameter orbit fit, we used a parallel-tempered sampler with 15 temperatures and 1500 walkers per temperature. We ran each walker for 125,000 steps, after an initial burn-in of 95,000 steps. Convergence was confirmed using the metrics defined in Section 3.1. We again thinned the chains by a factor of 75, and formed our posterior distribution from the lowest temperature chains. This resulted in 2,499,000 samples of the posterior. The posteriors are plotted in Figure 1 and reported in Table 3.

As we anticipated, the posterior for the eccentricity of planet d runs right up against our prior bounds. Without stability constraints, higher eccentricity orbits are favored. Just like in the Coplanar orbit fit, we find an orbital inclination for the system in the 20° to 30° range. However, Ω is now in agreement with that found in Konopacky et al. (2016) for coplanar orbits, unlike our previous orbit fits. It is likely this was a small family of orbits that was not represented in the 1σ range of our previous analyses.

3.5. Goodness of Fit

We used the reduced chi-squared ( ${\chi }_{\nu }^{2}$ ) statistic to measure the goodness of fit of a model. Since the highest likelihood model often does not represent the whole posterior of possible orbital configurations, we compute ${\chi }_{\nu }^{2}$ on 1000 randomly drawn allowed orbits for each model. We list the 16th, 50th, and 84th percentiles in Table 3. For the Unconstrained model, we found ${\chi }_{\nu }^{2}\approx 1$ , indicating the unconstrained Keplerian orbits can suitably describe the data as one might expect if the uncertainties are estimated accurately, given it is a physical model. The other three models have ${\chi }_{\nu }^{2}$ similarly close to unity, showing they also fit the data well.

We also investigated if more-restrictive models with additional, dynamically motivated constraints better describe the data than the fiducial Unconstrained case. We calculated the Bayesian Information Criteria (BIC; Schwarz 1978; Liddle 2007) as a simplified alternative to full Bayesian model comparison. The BIC assesses how well a model fits the data and penalizes models that have more free parameters. Models with lower BIC are preferred. We define the ΔBIC as the difference between the BICs of a more-restrictive model and the median BIC of the Unconstrained model. We also calculate ΔBIC using the same 1000 randomly drawn orbits for each model, and list the 16th, 50th, and 84th percentiles for this value in Table 3.

We find the ΔBIC is negative for the other three models relative to the Unconstrained fit. This indicates that adding constraints that tend the fits toward what we believe are stable orbits improves agreement with the data, as we discard some parameter space containing likely unstable orbits that do not reflect reality. We also note that ${\chi }_{\nu }^{2}$ and ΔBIC are not perfect metrics, as they only consider the number of free parameters in the models, and not the total parameter space being considered. In particular, when we limit the period ratios, this does not decrease the number of free parameters, while significantly limiting the space of possible orbits. Thus, we see these goodness of fit metrics favor the coplanar solutions as they explicitly reduce the number of parameters in the model. It would be better to have computed the Bayes factor between each pair of models to more rigorously compare models, but the Bayes factor is computationally difficult to calculate with a high-dimensional problem like this and our MCMC samplers were only set up to perform parameter estimation. Because of this, we do not think it is valid to conclude from solely these two metrics that coplanar orbits are favored. However, we can assert that adding constraints to the orbit fit does not worsen the fit from the fiducial case, and thus the constraints are reasonable given the current astrometric data. This conclusion agrees with the analysis from Konopacky et al. (2016), who found that coplanar orbits and orbits near 1:2:4:8 period ratios were fully consistent with the Keck astrometry. While the metrics we have employed cannot decide which orbit model should be favored, the stability constraints to the system that are investigated in the following section will clearly show what the realistic orbits are.

4. Dynamical Constraints

Keplerian motion is not the only constraint on the orbits of the planets. We also know that these four planets must also have survived from their formation up to this point. The latest estimate for the age of the star is ${42}_{-4}^{+6}$ Myr old (Bell et al. 2015), based on its membership in the Columba moving group (Torres et al. 2008; Zuckerman et al. 2011). This stellar age is further supported by interferometric measurements of the stellar radius (Baines et al. 2012). These orbits must have been stable for roughly the lifetime of the star, since giant planets likely formed quickly before the gas disk dispersed in the first few Myr (Williams & Cieza 2011). As the gas disk is difficult to model and exists for only a short period of the system’s lifetime, we do not simulate the time between planet formation and gas disk dispersal. There could be additional bodies in the system, but it is impractical to consider them without making assumptions about their nature. Instead, our analysis will focus on eliminating unstable orbital configurations based on the four planets alone. If additional bodies are detected, they could further constrain these orbits.

Thus, we investigated which orbital configurations allowed by astrometric measurements are also stable if we simulate the four planets’ orbits backwards in time for 40 Myr. In this section, we will apply this dynamical constraint on each of the four orbit fits (Unconstrained, Coplanar, Near 1:2:4:8, Near 1:2:4:8 Coplanar) and investigate the family of stable orbits that arise.

4.1. Stability of Orbital Models

We used the REBOUND N-body simulation package (Rein & Liu 2012) with the WHFast integrator (Rein & Tamayo 2015). To set up a simulation, we added particles for planets e, d, c, and b, in that order, using a chosen set of orbital parameters from our fits and placing them at the predicted location on MJD 56609, the date of the first GPI epoch. We drew masses for each of the planets in a process described in the following paragraph, and set the primary mass to be the stellar mass from our orbit fits. We then reversed the present velocities of the planets and integrated the system for 40 Myr to simulate the past dynamical history of the system, using fixed timesteps equal to 1% of planet e’s initial orbital period. We considered a configuration unstable if two planets passed too close, or if one planet was ejected from the system. We considered an encounter too close if any two planets passed with a distance less than the initial mutual Hill radius of planets d and e, which we approximated as

where M_e, M_d, and M_⋆ refer to the masses of planet e, planet d, and the primary, respectively. We considered a planet to be ejected if it moved further than 500 au from the star. Any orbit that survived for 40 Myr without encountering either condition is considered stable.

To assess the dynamical stability of allowed orbital configurations for each model, we performed rejection sampling to assess which orbit models contained significant amounts of stable orbits. We drew 1 million random orbital configurations from each of the four model posteriors. The orbit fits do not specify the mass of the planets, so we needed to add additional parameters for them. For simplicity, we set planets c, d, and e to be equal in mass, as we would expect due to their similar luminosities (Marois et al. 2010). We drew the mass of these planets, M_cde, from a uniform prior between 4 and 11 M_Jup to encompass the uncertainty on the luminosity-derived masses from Marois et al. (2010). We drew the mass of planet b, M_b, from a uniform prior between 3 M_Jup and M_cde to account for its lower luminosity. We will discuss using a more informative prior based on the planets’ luminosities in Section 4.4. We ran each of the configurations through the REBOUND setup described previously. Our dynamical stability prior sets the probability of an orbit to be 0 if the system is not stable, and 1 if it is stable, discarding the unstable orbits in our rejection sampling. In Table 3, we record the number of stable orbits from each of the configurations.

We found that the Unconstrained and Coplanar orbit solutions did not yield any stable orbits after 1 million draws. Especially for the Unconstrained case, the lack of stable draws does not mean that these models are inconsistent with stable orbits, but rather that the islands of stability in this high-dimensional space are small and were not sampled even after millions of MCMC draws. Simply, these models do not currently allow for a practical search of stable orbits.

Both models that assume that the planets’ orbital periods are near the 1:2:4:8 period ratio do yield stable orbits, with the model not assuming coplanarity, the Near 1:2:4:8 model, resulting in just a single stable orbit after 1 million draws. This model encompasses all of the parameter space explored by the Near 1:2:4:8 Coplanar model, so it is the more general model. We explored this model further by running 20 million REBOUND simulations in total, leading to 50 stable orbits.

We find that for masses of the inner three planets greater than 5 M_Jup, the maximum mutual inclination between any pair of planets in a stable system is <8°, although we only have few samples in this regime (14 stable orbits spanning 8 unique present-day orbital configurations). We also do not find stable orbits above 6 M_Jup, which likely reiterates the difficulty of finding stable noncoplanar orbits due to the high dimensionality of the problem. Thus, with the limited orbital arcs we have so far, looking for noncoplanar stable orbits is impractical. Since in Section 3 we found that our astrometry is consistent with the system being coplanar, we will focus on those orbits since we can find many stable orbits with this assumption (hundreds per million tries). We will leave the thorough exploration of orbits with mutual inclinations for future work with longer astrometric baselines and more computation time. However, in our preliminary analysis, it seems that the mutual inclinations are probably small in order for the system to be stable.

4.2. Stable Coplanar Orbital Solutions

For the rest of the analysis, we focus on the orbital parameters from the Near 1:2:4:8 Coplanar fit. We increase the number of N-body simulations from 1 million to 22 million, obtaining 9792 stable orbital configurations. We plot the initial osculating orbital elements (i.e., those on MJD 56609) of these stable orbits in Figure 3 and list them in Table 4.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Posterior of stable orbital elements for coplanar configurations of the four planets. These posteriors show all stable orbits with M_cde > 4 M_Jup and 3 M_Jup < M_b < M_cde. As discussed in Section 4.2, solutions with higher planet masses lie within a smaller region of this space.
Download figure:
Standard image High-resolution image

Table 4. Stable Coplanar Orbital Parameters of HR 8799 bcde

Body	Element	Coplanar
	Orbital	Stable
b	a_b (au)	${70.8}_{-0.18}^{+0.19}$
	τ_b	${0.46}_{-0.26}^{+0.31}$
	ω_b (°)	87 ± 58
	Ω_b (°)	${67.9}_{-5.2}^{+5.9}$
	i_b (°)	26.8 ± 2.3
	e_b	${0.018}_{-0.013}^{+0.018}$
c	a_c (au)	${43.1}_{-1.4}^{+1.3}$
	τ_c	${0.43}_{-0.24}^{+0.15}$
	ω_c (°)	${67}_{-39}^{+59}$
	Ω_c (°)	${67.9}_{-5.2}^{+5.9}$
	i_c (°)	26.8 ± 2.3
	e_c	${0.022}_{-0.017}^{+0.023}$
d	a_d (au)	${26.2}_{-0.7}^{+0.9}$
	τ_d	${0.839}_{-0.017}^{+0.020}$
	ω_d (°)	${17}_{-11}^{+12}$
	Ω_d (°)	${67.9}_{-5.2}^{+5.9}$
	i_d (°)	26.8 ± 2.3
	e_d	${0.129}_{-0.025}^{+0.022}$
e	a_e (au)	16.2 ± 0.5
	τ_e	${0.124}_{-0.013}^{+0.019}$
	ω_e (°)	110 ± 9
	Ω_e (°)	${67.9}_{-5.2}^{+5.9}$
	i_e (°)	26.8 ± 2.3
	e_e	${0.118}_{-0.028}^{+0.019}$
A	Parallax (mas)	${24.30}_{-0.69}^{+0.49}$
	M_⋆ (M_⊙)	${1.47}_{-0.08}^{+0.11}$

	${\chi }_{\nu }^{2}$	${1.01}_{-0.05}^{+0.06}$
	ΔBIC	$-{22}_{-4}^{+6}$

Note. Values are reported in the same way as Table 3.

Download table as: ASCII Typeset image

We find that the posteriors have tightened significantly after applying the dynamical stability constraint. Figure 4 visually compares the spread of possible orbits on the 2D sky plane for the orbit fits with increasing constraints placed on them. The stable coplanar orbits appear as a well-defined ellipse with minimal uncertainty for each planet’s orbit. This is also reflected visually and numerically in the posterior percentiles. The middle 68%, the difference between the 84th and 16th percentiles, of the semimajor axes of the planets decreased by 1.5–4.5 times when compared to the Unconstrained case, and by a factor of 1.17 to 1.50 when compared to the Near 1:2:4:8 Coplanar fit that the stable orbits were drawn from. Similarly, the middle 68% of the eccentricities also decreased by a factor between 2.2 and 4.7 compared to the Unconstrained fits. In fact, the fractional uncertainty on the semimajor axes is about the same as the fractional uncertainty of the Gaia DR1 parallax of the system (≈2%). The inclusion of the parallax from Gaia Data Release 2, released after this analysis was completed, should reduce its contribution to the semimajor axis and total system mass uncertainties by a factor of 7 (Gaia Collaboration et al. 2018). We have chosen not to rerun our analysis since the conclusions in this paper do not strongly depend on the exact semimajor axes of the orbits, and we will leave this for a future work.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Comparison of 200 allowed orbits from the Unconstrained (Section 3.1), Near 1:2:4:8 Coplanar (Section 3.4), and dynamically stable coplanar solutions (Section 4.2) projected onto the sky plane. The black star in the middle represents the location of the star, the black circles are the measured astrometry (uncertainties too small to show on this scale), and the current orbit for each planet is colored in the same way as in Figures 1 and 3 (i.e., planet b is blue, c is red, d is green, and e is yellow).
Download figure:
Standard image High-resolution image

The stable orbits, despite being much more restrictive, are good fits to the data. Aggregating 1000 random orbits, we find a ${\chi }_{\nu }^{2}={1.01}_{-0.05}^{+0.06}$ that is just as good as the fiducial Unconstrained model. The ΔBIC is similarly comparable to the other models. Although, once again, we note that BIC does not account for the narrower parameter space due to the additional stability constraint. We conclude that these stable orbits are a small, but allowed part of a much larger space that we have explored through our Bayesian analysis.

The masses of the stable configurations are plotted in Figure 5. We will discuss mass constraints in Section 4.4 in detail. Briefly here, we can see that stable orbits exist with the mass of the inner three planets at almost 9 M_Jup, and separately with the mass of planet b to be nearly 7 M_Jup. Also, the majority of stable orbits we found are low mass. A total of 95.6% of the orbits have M_cde < 6 M_Jup, and 73.0% of the orbits have M_cde < 5 M_Jup. This highlights the difficulty in finding stable high-mass solutions when starting with our current orbit fits.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Distribution of masses of the stable orbits from Section 4.2 (blue) and comparison to the priors from which the masses were drawn (gray). The main plot in the bottom left shows the 2D distribution of masses. The contour lines represent 15th, 35th, 55th, 75th, and 95th percentiles of the distribution, with everything outside the 95th percentile plotted individually as points. The top and right panels show 1D histograms for M_cde and M_b, respectively, with the frequency in each bin plotted on a logarithmic scale to highlight the high-mass bins. The gray priors are plotted in the same fashion as the blue posteriors.
Download figure:
Standard image High-resolution image

With these stable orbits, we can look at the mass dependence on the orbital parameters to justify our choices of prior constraints in doing the Near 1:2:4:8 Coplanar orbit fit. In Figure 6, we plot the range of period ratios and eccentricities of stable orbits as a function of M_cde. We see that for M_cde > 6 M_Jup, none of the period ratios or the eccentricities are close to the bounds set by our priors. Below 6 M_Jup, the period ratio of planet d to e as well as the eccentricities of d and e are near the upper bound in the extreme case, indicating our prior may be excluding some low-mass stable orbits. Since the interquartile range of these parameters is far away from these bounds, only a few extreme low-mass cases have been excluded, so the effect should be minimal. As the masses increase, we see the range in the allowed parameter space decreases, indicating that the highest mass stable orbits reside in a subspace of the parameters we are exploring. Thus, we conclude that we are not unnecessarily excluding stable orbit configurations with our choice of priors that were designed to improve the efficiency of finding stable orbits.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Distribution of current period ratios and eccentricities as a function of mass of the inner three planets for stable orbits. For M_cde < 7 M_Jup, the data is binned into one box plot per M_Jup. Each box shows the 25th, 50th, and 75th percentiles of the given distribution, while the whiskers show the extrema. Above 7 M_Jup, points are plotted individually, as they are sparse enough. The bounds of the priors are plotted as gray dot-dashed lines, except for the upper bound of the d:e period ratio, which is the yellow dotted line. These plots show how the range of allowed period ratios and eccentricities decrease as planet mass increases. Above M_cde > 6 M_Jup, the full range of stable orbits are not near the prior bounds.
Download figure:
Standard image High-resolution image

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Demonstration of the algorithm to identify librating and circulating segments of each critical angle. The top, middle, and bottom plots show a purely librating, transitioning, and purely circulating angle, respectively. The points colored red show regions identified as circulating, and the points colored blue show regions identified as librating. In this example, θ_d:e,e is always librating with a librating center of 0° and a libration amplitude of 36°. θ_c:d:e is librating only 33% of the time and θ_b:c,b is circulating 98% of the time, so librating amplitudes and centers are not well defined.
Download figure:
Standard image High-resolution image

There are several notable features in our posteriors of stable orbital configurations. The bimodality of the eccentricity posteriors is clear. The outer planets b and c have e ∼ 0, while the inner planets d and e have e ∼ 0.1. These eccentricities agree well with what was found by Goździewski & Migaszewski (2014), who migrated planets into resonance lock, rotated the orientations to match the astrometry, and selected orbital configurations with a ${\chi }_{\nu }^{2}$ cutoff. Given that this conclusion was reached by two completely different analysis methods, the fact the inner two planets have slightly eccentric orbits while the outer two planets are in near-circular orbits is a notable result that seems to be required for most stable orbital configurations that are consistent with the measured astrometry. The increased eccentricities of planets d and e, and the proximity of all four planets to 1:2:4:8 period commensurabilities, are consistent with an early evolutionary period of convergent inward migration of all four planets, trapping of planet pairs d & e and c & d into 2:1 resonances, and pumping of the orbital eccentricities of d and e by continued migration while in resonance lock (e.g., Yu & Tremaine 2001; see also Section 4.3).

Also, comparing our stable orbits with those of Goździewski & Migaszewski (2014) and Gozdziewski & Migaszewski (2018), we note that most orbital parameters agree fairly well except for the semimajor axes of the planets, which we find to be significantly larger. For example, only 0.11% of our orbital solutions have a_c ≤ 39.4 au, the best fit solution of Goździewski & Migaszewski (2014). Our uncertainties in parallax and stellar mass are consistent with the fixed parallax and stellar mass they used, so the difference in a is not just a result of different system parameters. We also generally have larger uncertainties on our values, which can be due to a combination of allowing lower mass orbits, not strictly enforcing 1:2:4:8 resonance lock, and a more systematic exploration of parameter space.

As the orbits are coplanar, each planet’s argument of periastron can be used to measure the relative orientation of the planets’ orbits. As both planet b and c have extremely low eccentricities, ω is basically unconstrained for these planets and is not notable since their orbits are near circular. The significant nonzero eccentricity of planets d and e, however, corresponds to sharp peaks in ω for both planets. While a broader peak was already seen in ω_d in the Near 1:2:4:8 Coplanar fits before enforcing stability constraints, the addition of dynamical stability has disallowed circular orbits of planet e, giving rise to a sharp peak in ω_e. Interestingly, the orientation of the orbits of planets d and e are not aligned, with ${\omega }_{e}-{\omega }_{d}={94}_{-9}^{+11}$ degrees, essentially perpendicular to being aligned. We note that Goździewski & Migaszewski (2014) also found a similar result.

The period ratios of the planets shown in Figure 6 heavily favor period ratios above the nominal 2:1. We find period ratios P_b/P_c = 2.11 ± 0.06, ${P}_{c}/{P}_{d}={2.10}_{-0.05}^{+0.04}$ , and ${P}_{d}/{P}_{e}\,={2.06}_{-0.04}^{+0.03}$ . For P_c/P_d and P_d/P_e, the data favors period ratios above 2. This can been seen by computing the period ratios using the median a for the Near 1:2:4:8 Coplanar fits from Table 3 that do not have a dynamical prior applied. It seems that these period ratios are at these high values to satisfy the astrometry. For P_b/P_c, the data allows both period ratios above and below 2, so having it strongly favor values above 2 (only 2% of stable configurations have P_b/P_c < 2) indicates that spacing the planets slightly farther apart enhances stability. The period ratios driven by the astrometry could be indicative of a primordial period ratio. In particular, the planets could have experienced eccentricity damping while in resonance and were hence repelled to period ratios greater than 2 while still maintaining resonance lock (Lithwick & Wu 2012; Batygin & Morbidelli 2013). As disk gas is the primary mechanism for eccentricity dissipation of Jupiter-mass planets at large-separations,³⁵ this may indicate that the planets were in or near their current location during the gas disk stage (Dong & Dawson 2016).

We found the system has an inclination of i = 26 fdg 8 ± 2 fdg 3 and longitude of ascending node of ${\rm{\Omega }}=68\buildrel{\circ}\over{.} {0}_{-5.3}^{+5.9}$ , consistent with the work by Konopacky et al. (2016) fitting coplanar orbits, but a bit more precise. Both i and Ω match the debris disk inclination of 26° ± 3° and position angle of 62° ± 3° derived from far-infrared Herschel observations (Matthews et al. 2014). While the inclination is also consistent with the millimeter observations of the debris disk by the Submillimeter Array (SMA) and the Atacama Large Millimeter/Submillimeter Array (ALMA), Ω is higher than the position angle of the disk of $35\buildrel{\circ}\over{.} {6}_{-10.1}^{+9.4}$ in the millimeter (Wilner et al. 2018). This implies the disk in the millimeter is mutually inclined from the planetary orbital plane by ${16}_{-11}^{+22}$ degrees. If this offset is real, we would be observing a process that decouples the millimeter planetesimals from the planets and smaller dust grains probed by Herschel. However, it is not clear what could cause this. As reported in Section 4.1, we do not find high-mass solutions that are mutually inclined by more than 8°, so it seems unlikely that planet b is torquing the disk, although we cannot definitely exclude planet b being inclined from the rest of the planets. Still, it would not explain why the smaller dust seen by Herschel are indeed coplanar with the planets. Deeper observations of the debris disk are needed to determine if the debris disk is coplanar to the planets.

Previous works have considered the need for an additional planet to carve the inner edge of the outer belt of debris (Booth et al. 2016; Read et al. 2018). Here we investigate whether planet b is consistent with sculpting the inner edge of the outer belt, assuming the planets are coplanar or nearly so to the disk; the possible millimeter-wave offset of ${16}_{-11}^{+22}$ deg cited previously is assumed negligible in this regard. Since planet b’s orbit is likely near circular, with 95% of the allowed stable orbits having e_b < 0.05, we can compute the clearing zone of planet b using the following equation from Morrison & Malhotra (2015) for the outer edge of a planet’s chaotic zone, validated for high-mass planets like HR 8799 b:

Here R_in is the inner radius of the disk and corresponds to the outer edge of a planet’s clearing zone, a_p is the semimajor axis of the planet, M_p is the mass of the planet, and M_⋆ is the mass of the star. Plugging in the numbers from our dynamically stable orbits, we find ${R}_{\mathrm{in}}={89}_{-2}^{+3}$ au when considering all stable orbits, and ${R}_{\mathrm{in}}={93}_{-2}^{+3}$ au when considering only stable orbits with M_b > 5 M_Jup. When compared to the inner edge of ${104}_{-12}^{+8}$ au derived by Wilner et al. (2018), our median value of the inner edge when only considering M_b > 5 M_Jup is consistent with their middle 68% credible interval, while our median value when considering all of our stable solutions equally is slightly below this credible interval. However, given that our quoted numbers on the inner disk edge depend on our priors on the mass of the planet, and given the uncertainty in the inner disk edge (Booth et al. 2016; Wilner et al. 2018), a 1σ disagreement is not significant. Our orbit fits place planet b at a location consistent with sculpting the inner edge of the debris disk, although finer studies of the dynamical interactions of system and more refined system parameters will help clarify the picture.

4.3. Orbital Resonances

Having stable orbits near integer period ratios does not guarantee resonance. To explore possible orbital resonances in our stable configurations, we saved the state of each stable system every 200 years using the SimulationArchive feature of REBOUND (Rein & Tamayo 2017). We looked at the resonant angles of each system as a function of time to infer the resonant nature of the system: planets in resonance will have a corresponding resonant angle that is librating, but the angle will circulate if the planets are not in resonance. For this four-planet system, we considered nine resonant angles. The first six angles are two-body resonant angles that look at whether consecutive planets are in 2:1 resonance. These 2:1 resonant angles are defined by

Here, ϖ = Ω + ω is the longitude of periastron, and λ = ϖ + M is the mean longitude, where M is the mean anomaly. The labels i and j refer to the labels of a consecutive pair of planets, with the inner planet being j, and k refers to either i or j, resulting in two resonant angles per pair of planets and thus six two-body resonant angles in total for the three consecutive pairs of planets. For example, the 2:1 resonant angle for planets c and d using planet d’s ϖ would be written as θ_c:d,d. In a similar notation, the three-body 1:2:4 Laplace resonance can be written as

where the innermost planet is k and the outermost planet is i. Lastly, we looked at the same four-body 1:2:4:8 resonant angle as Goździewski & Migaszewski (2014):

In our simulations, we found that these resonant angles varied in behavior, with some continuously librating (i.e., locked in resonance for 40 Myr), some continuously circulating, and some transitioning between the two over the 40 Myr orbit integration. To analyze all of the simulations uniformly, we developed an algorithm to identify libration and compute the fraction of time a resonant angle is librating or circulating over the course of a simulation. The algorithm takes advantage of the fact that librations oscillate around a fixed value while circulating angles are monotonically changing. Briefly, the algorithm uses a Fourier transform to identify the periodicity of the data, smooths it on that scale, and computes the time derivative of the smoothed angle over the time series. Any sections of the time derivative with significant deviations from zero are deemed circulating, and the rest are deemed librating. Figure 7 shows an example of this algorithm classifying librating and circulating sections of a few resonant angles. We note that this method is not perfect and requires a subjective threshold to determine when a deviation is significant. However, inspecting the results from several resonant angles by eye, the algorithm seems comparable to by-eye identification, and seems to accurately identify resonant angles that are librating continuously (i.e., librating 100% of the time). For systems that by eye are transitioning quickly or that are always circulating, we estimate that we misidentified ∼5% of the time series. This error is small, and the gain in having an automated algorithm to uniformly analyze all of these time series is large.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Distribution of systems in the space of libration fraction versus planet mass. Contours are 25th, 50th, 75th, and 97th percentiles (e.g., 97% of systems lie within the 97th percentile contour). Above that, individual points that correspond to particular stable orbital configurations are plotted. Libration fractions of 100% indicate resonance lock, while libration fractions <5% indicate the planets are probably never in resonance.
Download figure:
Standard image High-resolution image

We apply this algorithm to all nine resonant angles for each of our simulations. We plot the libration fraction, the fraction of time in the last 40 Myr during which that angle is librating, for each angle in Figure 8. At low masses (M_cde < 6 M_Jup), the scatter in the libration fraction is high for all angles, indicating resonance lock for any subset of the planets is not necessary for stable orbits at low masses. Several of the angles never reach 100% libration fraction at any mass, indicating that all of the stable orbital configurations we found do not have all four planets in resonant lock. While we do not see any four-planet resonant chains, it might be possible they reside in a small island of parameter space that our MCMC did not sample. All resonance angles involving planet b never reach 100% libration fraction. Planet b may occasionally come into resonance with the inner planets, but does not remain there. It is not too surprising that planet b does not favor resonance, insofar as the magnitude of the resonant potential associated with θ_b:j,b is proportional to the planet’s orbital eccentricity e_b (e.g., Murray & Dermott 1999), and its median eccentricity is the lowest among the four planets.

The inner three planets do favor resonance more than planet b. Above ∼6 M_Jup, a large majority of stable orbits have θ_d:e,e librating 100% of the time. Similarly, above ∼7 M_Jup, θ_c:d,d is always librating for most stable orbits. For all masses, these two angles are always librating some of the time. However, the other two resonant angles, θ_d:e,d and θ_c:d,c, appear primarily transitioning between libration and circulation, with θ_d:e,d trending to libration at high masses. This behavior is also reflected in the three-body resonant angle between the inner three planets, with only 1.2% of the stable orbits having this three-body angle librating for at least 90% of the time. Still, this behavior indicates that the inner planets being in a 1:2:4 three-body resonance is both consistent with the data and dynamically stable for masses up to 8 M_Jup. In these cases, the libration center of θ_c:d:e often jumps between ∼90° and ∼−90°, but typically keeping $| {\varpi }_{d}-{\varpi }_{e}| \sim 90^\circ$ . However, the 1:2:4 three-body resonance is not required for stability, even at high masses.

When the three planets are not locked in resonance together, pairs of planets can be in resonance. We find, averaged across the ensemble of simulations, these two-body resonant angles librate around 0°. However, in a single simulation, the libration center can be offset from 0°, a phenomenon know as asymmetric libration (e.g., Murray-Clay & Chiang 2005) that is observed for two-body angles in other resonant chain systems (e.g., Kepler-80, MacDonald et al. 2016) and is caused by the gravitational effect of a third planet. For the situation where θ_d:e,e librates but θ_d:e,d does not, the conjunction of planets d and e always occurs at the periastron of planet e’s orbit, but is completely uncorrelated with planet d’s orbit. In the case where both θ_d:e,e and θ_c:d,d librate but the three-body angle θ_c:d:e does not, planet e’s orbit orients itself so that it is lined up to the conjunction with planet d, while planet d’s longitude of periastron is driven by the conjunction with planet c, which is not locked in with planet e. In this case, consecutive planet pairs appear to be locked in resonance for 40 Myr, but three-body resonance lock does not exist.

Having planets d and e and planets c and d locked in two-body resonances fits well with the picture that they were locked in resonance quickly after formation, before the gaseous protoplanetary disk disappeared. After the planets migrated into resonance lock, the eccentricities of e and d were both amplified by their resonant migration in the gas disk and damped by the gas, pushing the planets to larger period ratios. After the gas dispersed, the planets maintained their primordial eccentricities and period ratios, with these parameters only oscillating as the planets exchange energy and angular momentum in resonance.

4.4. Dynamical Mass Limits

As we only have a short orbital arc of data, we are limited on the mass constraints we can place based solely on dynamical considerations. We have not yet measured the perturbations of the planets’ orbits by each other. Without seeing a significant effect, we cannot place a lower bound on the masses of the planets dynamically. Impractically, we may need to measure the change in the orbital elements of the planets after many orbits, akin to the masses derived from transit timing variations (Agol et al. 2005; Holman & Murray 2005). Thus, dynamical constraints based on short orbital arcs cannot fully constrain the masses alone.

Looking at Figure 5 again, our stable orbits heavily favor low masses, since we cannot place a lower bound on the masses and we used a uniform prior. For high masses, M_cde > 8 M_Jup, we have only a few stable orbits, all with the mass of the planet b near the lower bound of what we would expect, indicating a possible upper limit to the masses of the planets. However, we cannot verify if this is because of a lack of sampling of stable orbits. That is, the probability of drawing a stable orbital configuration with M_cde > 8 M_Jup and M_b > 5 M_Jup might be so small that we do not expect to find one with our current sampling. We are potentially limited by the fact that 2.5 million samples of the posterior are not sufficient when the posterior has 20 dimensions and the islands of stability at extremely high masses are extremely small. Thus, there is no indication of a sharp drop-off that would point to a firm upper limit on the mass based on dynamical considerations. In fact, Goździewski & Migaszewski (2014; 2018) found stable orbits above 9 M_Jup for the inner three planets, indicating higher-mass stable orbits exist if one forces the entire system to be in a four-planet resonance lock. It might be that we missed those systems because we did not enforce such a global lock and therefore had a larger parameter space of allowed stable orbits. If we acknowledge that we are dependent on our choice of priors for the mass of the planets, we can say that 99.9% of the orbits that are dynamically stable for the last 40 Myr have M_cde < 7.6 M_Jup and M_b < 6.3 M_Jup for this particular choice of prior. These mass upper limits are consistent with the luminosity-derived masses of M_cde = 7 M_Jup and M_b = 5 M_Jup based on hot-start evolutionary models (Marois et al. 2008, 2010).

We also investigated using a mass prior based on the measured luminosity and hot-start evolutionary models. Since it is computationally intensive to rerun all 22 million N-body simulations, we instead weighted each sample drawn from the previous prior distribution with the relative change in probability due to switching to a luminosity-based prior, resulting in a down-weighting of the lowest mass configurations, which are disfavored by the measured luminosities. To accomplish this, we reevaluated the luminosity-derived masses with newest age estimates for the Columba moving group from Bell et al. (2015), which we approximate as a Gaussian distribution of 42 ± 5 Myr. This age has better stated uncertainties than the 30–60 Myr range given in Marois et al. (2010). We still used the same luminosities measured by Marois et al. (2008, 2010), as there has not been an update to them with stated uncertainties. Using the Baraffe et al. (2003) hot-start cooling tracks, we get model-dependent masses of 6.0 ± 0.7, 8.7 ± 1.0, 8.7 ± 1.0, and 8.7 ± 1.7 M_Jup for planets b, c, d, and e, respectively. To stay self-consistent with our simulations that fix the masses of the inner three planets to be the same, we choose to use the 8.7 ± 1.0 M_Jup of planets c and d for the mass of the planet e also. Otherwise, planet e will tend toward lower masses and be inconsistent with our simulation assumption of equal masses for the inner three planets. Indeed, spectrophotometric measurements from the latest generation of high-contrast imagers confirm that planet e has similar near-infrared fluxes to planet c and d, and not planet b (Zurlo et al. 2016; Greenbaum et al. 2018), so our assumption should be robust. Altogether, when combining the luminosity and dynamical constraints for the mass assuming hot-start evolutionary tracks, we get a mass of 5.8 ± 0.5 M_Jup for planet b and ${7.2}_{-0.7}^{+0.6}\,{M}_{\mathrm{Jup}}$ for planets c, d, and e. These mass distributions are plotted in Figure 9. Note that our mass estimates depend on the resonances we found the planets to be in, and that higher masses can be achieved by assuming a four-planet resonance lock (Goździewski & Migaszewski 2014; Gozdziewski & Migaszewski 2018).

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** Histogram of dynamically stable masses after being weighted by the luminosity prior on the planet masses assuming hot-start evolutionary tracks. Planet b is in blue, and planets c, d, and e are assumed to follow the same red histogram.
Download figure:
Standard image High-resolution image

It is uncertain exactly how bright planets are during the first 100 Myr, as this depends on uncertain formation mechanisms. Planet cooling tracks are instead paramterized by a quantity like the initial entropy of the material that formed the planet (Spiegel & Burrows 2012; Marleau & Cumming 2014). Since the hot-start models really are the high-entropy upper limit with regard to planet formation models, the masses estimated assuming these tracks are the lowest masses for the planets. Thus, we can use the hot-start model to quote a lower limit on the mass. Combining the dynamical constraints with the luminosity prior on the masses from the hot-start model, we find 95% of the stable systems have M_b > 4.9 M_Jup and M_cde > 6.1 M_Jup. Alternatively, we can use the upper limits on the masses from dynamical stability alone to constrain the initial conditions of the cooling tracks. Using the same stated luminosities and age of the system, and also using the Spiegel & Burrows (2012) warm-start models, our upper limits of M_cde < 7.6 M_Jup and M_b < 6.3 M_Jup would then correspond to lower limits on the initial entropy of 9.5 k_B per baryon for the inner three planets and 9.2 k_B per baryon for planet b. This excludes the most-extreme cold-start formation models, but is consistent with a range of higher entropy models, as warm- and hot-start models have similar luminosities at this age.

4.5. Long-term Dynamical Stability

Even though the system has been stable for ∼40 Myr, we investigate whether the system we see today is reflective of the final state of system, or whether this configuration is a transient one. We used the saved SimulationArchive of each stable orbit, reversed the velocities again, and now integrated them forward in time for up to 1 Gyr. We used the same criteria to assess stability as before. We plot the amount of time in the future each system is stable as a function of the mass of the inner three planets in Figure 10.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Stability time of the dynamically stable orbits integrating forward in time for up to 1 Gyr. Contours are 25th, 50th, 75th, and 97th percentiles. The horizontal blue dotted line indicates 40 Myr, the time we integrated the systems backward as a reference.
Download figure:
Standard image High-resolution image

We find orbits that are stable for 1 Gyr up to 6 M_Jup. Above 7 M_Jup, none of our orbits are stable for more than 100 Myr. Thus, if these planets have masses greater than 7 M_Jup, as is favored by our combined dynamical and luminosity constraint, then the system is not likely stable. As there is evidence in our own solar system of dynamical upheaval of our less-tightly packed gas giants early on (Tsiganis et al. 2005), it is not surprising to find that the HR 8799 system will become unstable.

Certain resonances do seem to improve stability. In particular, in systems where the three-body angle θ_c:d:e is librating at least 50% of the time, systems are 10 times less likely to go unstable in 1 Myr and 4 times less likely to become unstable in the 10 Myr than systems with θ_c:d:e librating less than 50% of the time. There does not appear to be strong correlation between the two-body angles and stability except for θ_b:c,c. Systems for which θ_b:c,c is libration more than 30% of the time are 5 and 15 times more likely to be stable for at least 1 Myr and 10 Myr, respectively, than systems where this angle librates less than 30% of the time. While certain resonances seem to prevent systems from short-term instabilities, there is no indication that spending more time in resonance or achieving resonance lock improves stability at the 1 Gyr level. This is perhaps due to the fact that we did not find systems with the four planets locked in resonance, which could have improved longer term stability.

5. Conclusion

This paper has aimed to explore the dynamically stable orbits of the HR 8799 system. In the first part of the paper, we continued to demonstrate the precise astrometry that can be achieved by GPI and explored various assumptions on the orbits of the planets.

1.
Using GPI IFS data from 2014 to 2016, we measured the astrometry of the HR 8799 with 1 milliarcsecond precision using the open-source pyKLIP package.
2.
We utilized MCMC methods to fully explore the 20+ dimensional space of orbital configurations using Bayesian parameter estimation and different assumptions on the coplanarity and resonant nature of the system.
3.
We found that assuming the system is coplanar or the system is near the 1:2:4:8 period ratio resonance does not significantly worsen the fit to the data, and in fact might make it better. We find including both assumptions provides adequate fits to the data, agreeing well with the conclusions of Konopacky et al. (2016).

In the second half of this paper, we presented the first attempt to rigorously fold dynamical constraints into orbit fits of a directly imaged system, and demonstrated the power of including a dynamical prior.

1.
We performed rejection sampling on our posteriors of orbit fits to apply our prior of dynamical stability. Using the REBOUND N-body integrator, we ran orbits backward in time for 40 Myr, the age of the system, varying the masses of the planets, and looked for the stable orbital configurations.
2.
We find that coplanar orbits near the 1:2:4:8 resonance produce orders of magnitude more stable orbits than any other scenario. We find a few orbits near 1:2:4:8 resonance with some mutual inclinations that are stable, but the inefficiency of finding them makes studying that family of orbits impractical with current astrometric data.
3.
As demonstrated by Figure 4, the stable coplanar orbits lie within a small fraction of the allowed orbital space. In this subspace, we find the outer two planets have near zero eccentricity, while the inner two planets have e ∼ 0.1.
4.
Our orbits are consistent with being coplanar with the Herschel-derived debris disk plane, but misaligned with the plane derived from SMA and ALMA by ${16}_{-11}^{+22}$ degrees. Our fitted orbit for planet b is consistent with that planet sculpting the inner edge of the debris disk in the millimeter, assuming the orbits are close enough to coplanar.
5.
If M_cde ≳ 6 M_Jup, planet e needs to be locked in resonance with planet d in order for the system to be stable. Likewise, if M_cde ≳ 7 M_Jup, planet d is likely in resonance lock with planet c. Although we find stable configurations where the inner three planets are in a 1:2:4 Laplace resonance, such a three-body resonance is not required, as we found many stable configurations where only pairs of planets are in resonance. Planet b does not need to be in resonance for this system to be stable so far.
6.
Using a uniform prior on M_cde and a slightly low-mass-favored prior on M_b, we find 99.9% of our stable orbits have M_b < 6.3 M_Jup and M_cde < 7.6 M_Jup. Folding in the mass constraints from the planet luminosities, hot-start evolutionary models, and a system age of 42 ± 5 Myr, we find M_b = 5.8 ± 0.5 M_Jup and ${M}_{{cde}}={7.2}_{-0.7}^{+0.6}\,{M}_{\mathrm{Jup}}$ . Either way, our mass constraints are consistent with hot-start evolutionary tracks.
7.
We do not find systems with the inner planets above 7 M_Jup that are stable for the next 1 Gyr.

In the future, as more of the orbital arcs are traced out with precise astrometry, it will become clearer where in the 20+ dimensional space the planets’ true orbital configuration lie. In the meantime, using more computational power, we can attempt search for stable orbits that are not forced to be exactly coplanar, are stable at higher masses, or are stable at other resonances. Additionally, the new parallax value from Gaia Data Release 2 (Gaia Collaboration et al. 2018) should help tighten the constraints on the semimajor axes of the planets and total mass of the system for the stable coplanar solutions. Even a decade since its discovery, the HR 8799 planetary system is one of the most unique and interesting systems that we know, and combining both detailed dynamical studies with atmospheric characterization will help us understand how these planets formed and how they will interact.

Our technique of performing rejection sampling to apply dynamical constraints after MCMC sampling of the orbital parameters can also be applied to other directly imaged multi-planet systems to better constrain the orbits with just a short orbital arc. Here we have shown that even with <15% of the full orbit covered, we can constrain orbital parameters to a few percent. This can remove the orbital uncertainty that comes with exoplanets discovered through imaging alone, where we typically need to wait for long-period planets to complete a orbital revolution before fully constraining its orbit. This is also potentially valuable for future space-based imaging missions that search for exo-Earths in multi-planet systems, since this method can estimate the mass and orbit with a small orbital arc. This can allow for the mission to better prioritize which exo-Earth candidates are observed with expensive spectroscopic observations by determining which are most likely Earth-mass and orbiting at a favorable distance from the star.

Simulations in this paper made use of the REBOUND code, which can be downloaded freely at http://github.com/hannorein/rebound. We thank Dan Tamayo for offering REBOUND tutorials and help with setting up the simulations. We also thank Eve Lee for helpful discussions on dynamics. J.J.W., J.R.G., P.K., and R.J.D.R. were supported in part by NSF AST-1518332, NASA NNX15AC89G, and NNX15AD95G. This work benefited from NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate. R.I.D. gratefully acknowledges support from NASA XRP 80NSSC18K0355. The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium. 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.

The GPI project has been supported by Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), 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). This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

Software: Gemini Planet Imager Data Reduction Pipeline (Perrin et al. 2014, 2016), pyKLIP (Wang et al. 2015), emcee (Foreman-Mackey et al. 2013), REBOUND (Rein & Liu 2012), Astropy (The Astropy Collaboration et al. 2018), matplotlib (Hunter 2007).

Facility: Gemini: South (GPI). -

Footnotes

35
Note, however, that damping of planetary eccentricity by gas dynamical friction does not conserve the planet’s orbital angular momentum, contrary to damping of eccentricity by tidal dissipation in the planet; the latter, not the former, is considered by Lithwick & Wu (2012) and Batygin & Morbidelli (2013). However, the general mechanism of resonant repulsion also occurs for eccentricity damping by gas dynamical friction (e.g., as simulated for giant planets at wide separations by Dong & Dawson 2016).

Please wait… references are loading.

Dynamical Constraints on the HR 8799 Planets with GPI

Authors

Author notes

Notes

Article metrics

Permissions

Share this article

Article information

Dates

Abstract

Related links

1. Introduction

2. Observations and Data Reduction

3. Orbit Fitting

3.1. Unconstrained Orbits

3.2. Coplanar Orbits

3.3. Near 1:2:4:8 Period Ratio Orbits

3.4. Near 1:2:4:8 Period Ratio Coplanar Low-e Orbits

3.5. Goodness of Fit

4. Dynamical Constraints

4.1. Stability of Orbital Models

4.2. Stable Coplanar Orbital Solutions

4.3. Orbital Resonances

4.4. Dynamical Mass Limits

4.5. Long-term Dynamical Stability

5. Conclusion

Footnotes

Show References