VLT/SPHERE robust astrometry of the HR8799 planets at milliarcsecondlevel accuracy
Orbital architecture analysis with PyAstrOFit^{⋆}
^{1} Space sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, 19c Allée du Six Août, 4000 Liège, Belgium
^{2} European Southern Observatory, Alonso de Cordova 3107, Vitacura, Casilla 19001, Santiago de Chile, Chile
^{3} Department of Astronomy, California Institute of Technology, 1200 E. California Blvd, MC 24917, Pasadena, CA 91125, USA
^{4} Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA
^{5} Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Received: 15 April 2016
Accepted: 7 October 2016
Context. HR8799 is orbited by at least four giant planets, making it a prime target for the recently commissioned SpectroPolarimetric Highcontrast Exoplanet REsearch (VLT/SPHERE). As such, it was observed on five consecutive nights during the SPHERE science verification in December 2014.
Aims. We aim to take full advantage of the SPHERE capabilities to derive accurate astrometric measurements based on Hband images acquired with the InfraRed Dualband Imaging and Spectroscopy (IRDIS) subsystem, and to explore the ultimate astrometric performance of SPHERE in this observing mode. We also aim to present a detailed analysis of the orbital parameters for the four planets.
Methods. We performed thorough postprocessing of the IRDIS images with the Vortex Imaging Processing (VIP) package to derive a robust astrometric measurement for the four planets. This includes the identification and careful evaluation of the different contributions to the error budget, including systematic errors. Combining our astrometric measurements with the ones previously published in the literature, we constrain the orbital parameters of the four planets using PyAstrOFit, our new opensource python package dedicated to orbital fitting using Bayesian inference with MonteCarlo Markov Chain sampling.
Results. We report the astrometric positions for epoch 2014.93 with an accuracy down to 2.0 mas, mainly limited by the astrometric calibration of IRDIS. For each planet, we derive the posterior probability density functions for the six Keplerian elements and identify sets of highly probable orbits. For planet d, there is clear evidence for nonzero eccentricity (e ~ 0.35), without completely excluding solutions with smaller eccentricities. The three other planets are consistent with circular orbits, although their probability distributions spread beyond e = 0.2, and show a peak at e ≃ 0.1 for planet e. The four planets have consistent inclinations of approximately 30° with respect to the sky plane, but the confidence intervals for the longitude of the ascending node are disjointed for planets b and c, and we find tentative evidence for noncoplanarity between planets b and c at the 2σ level.
Key words: planetary systems / stars: individual: HR8799 / methods: data analysis
© ESO, 2017
1. Introduction
Since its discovery by Marois et al. (2008), the HR8799 planetary system has been and still remains one of the most intriguing among the thousands of known planetary systems. Composed of at least four giant planets in a range of angular separations of approximately to (Marois et al. 2010b), and of two dusty debris belts (Su et al. 2009; Hughes et al. 2011; Matthews et al. 2014; Booth et al. 2016), it has been the focus of many different studies, including dynamical stability analyses to constrain the global orbital motion and estimate the masses of the four planets (see e.g. Goździewski & Migaszewski 2009, 2014; Reidemeister et al. 2009; Fabrycky & MurrayClay 2010; Soummer et al. 2011; Currie et al. 2012, 2014; Maire et al. 2015). This dynamical approach allows the orbits of the four planets to be simultaneously constrained, but requires strong assumptions, such as coplanar (but eccentric) or circular (but not necessarily coplanar) orbits. The individual analysis of each planet offers an alternative method to constraint the orbital architecture. To this aim, nonlinear leastsquares fits of Keplerian elements (semimajor axis a, eccentricity e, inclination i, longitude of ascending node Ω, argument of the periastron ω, and time of periastron passage t_{p}) have been performed (see e.g. Lafrenière et al. 2009; Bergfors et al. 2011; Esposito et al. 2013; Zurlo et al. 2016).
Recently, Pueyo et al. (2015) proposed an indepth analysis of the HR8799bcde orbital motion. The authors carried out a Bayesian analysis based on MonteCarlo Markov chain (MCMC) techniques adopting both a Metropolis Hastings algorithm (MacKay 2003; Ford 2005, 2006) and an affineinvariant ensemble sampler (ForemanMackey et al. 2013). This approach echoes the works published in Chauvin et al. (2012) for β Pictoris b, in Kalas et al. (2013) for Fomalhault b and more recently in Beust et al. (2016) for Fomalhault b and PZ Telescopii B. Among other things, Pueyo et al. (2015) discussed the coplanarity of the system, the orbital eccentricities of the planets, the possibility for mean motion resonances, and the role of HR8799d in possible dynamical interactions during the youth of this system. They also estimated the dynamical masses of HR8799bcde by computing the fraction of allowable orbits that pass the socalled closeencounter test. As pointed out in Pueyo et al. (2015), unaccounted biases and/or systematically underestimated error bars on the planets astrometry affect the MCMC results (see e.g. Givens & Hoeting 2012) and may lead to a biased estimation of the confidence intervals for the orbital parameters. Studying the astrometric history of HR8799 reveals indeed that the errors affecting some positions are most probably underestimated, as one can readily identify pairs or sets of positions that are not consistent with each other within their error bars, or cannot be modeled with a unique orbit. This was one of the incentives of the study presented by Konopacky et al. (2016), who very recently rereduced all of the Keck/NIRC2 observations of HR8799 to produce a selfconsistent data set free from variable instrumentrelated biases. This consistent data set was then used to derive updated probability distributions for the elements of the planetary orbits based on Monte Carlo simulations.
With the advent of secondgeneration highcontrast planet imagers such as the SpectroPolarimetric Highcontrast Exoplanet REsearch (SPHERE, Beuzit et al. 2008) at the Very Large Telescope (VLT), obtaining astrometric measurements of directly imaged planets is now becoming routine. It is therefore more important than ever that the methods used to derive such astrometric measurements include a careful estimation of all error sources, including systematic biases that are expected to affect even the most advanced planetimaging instruments. Here, we propose to derive the astrometry of the four HR8799 planets based on a data set obtained with SPHERE during its science verification phase in December 2014. While this data set was already analyzed and presented in Zurlo et al. (2016) and Apai et al. (2016), our aim here is to propose a detailed description of all individual contributions to the astrometric error budget, including systematic biases, and to derive general recommendations for future studies aiming at an accurate estimation of astrometric error bars. In Sect. 2, we start by describing the observations, data reduction and image processing steps that allow the four planets to be revealed with a high signaltonoise ratio (S/N). Then, Sect. 3 discusses our method to derive the astrometry of the four planets, gives a thorough description of all major sources of astrometric errors, and evaluates their respective contribution. We present in Sect. 4 the new opensource PyAstrOFit package, fully dedicated to orbital fitting based on Bayesian inference using the MCMC approach, which we use to perform an updated analysis of the individual orbits of the four planets. Some aspects of our results differ from previous analyses published in the literature. A short discussion of their implication on the orbital dynamics of the system is given before concluding in Sect. 5.
2. Observations and data reduction
2.1. Observations
SPHERE performs highcontrast imaging by combining an extreme adaptive optics system (Fusco et al. 2006), several coronographic masks, and three science subsystems including the InfraRed Dualband Imager and Spectrograph (IRDIS, Dohlen et al. 2008). The observations of HR8799 were performed during five consecutive nights from December 4 to 8, 2014, using IRDIS in the broadband H filter (1.48−1.77 μm) with an apodized Lyot mask (Soummer 2005; Carbillet et al. 2011; Guerri et al. 2011) of diameter 185 mas together with an undersized Lyot stop. A beam splitter located downstream from the coronagraphic masks produces two identical parallel beams (Beuzit et al. 2008), which results in two wellseparated images per acquisition, hereafter referred to as the left and right images. Each of the five observing sequences lasted approximately half an hour, and consisted of 218 frames with a detector integration time (DIT) of 8 s per frame. All observing sequences were obtained under fair seeing conditions (between and ), except on December 7 where the seeing was above . The sequences were acquired in pupilstabilized mode to take advantage of the angular differential imaging (ADI, Marois et al. 2006) technique. Due to the low elevation of HR8799 as seen from Cerro Paranal (maximum altitude of 44°), the amount of parallactic angle rotation was however quite small, amounting to , , , and for the five nights, respectively. Four elongated diffraction spots, the socalled satellite spots, were created during the whole observing sequences by injecting a waffle pattern on the deformable mirror (Langlois et al. 2012) to help with the starcentering procedure, as explained in the next section.
2.2. Data reduction
The IRDIS raw frames were preprocessed using the SPHERE EsoRex pipeline. As a first step, master dark and flat frames were created from calibration data obtained for each night of observations. Then, EsoRex identified the outlying pixels in the master dark frame by using a sigma clipping procedure and built a static badpixel map. Each frame was reduced by subtracting the corresponding master dark, dividing by the master flat and interpolating the pixels flagged in the badpixel map. At this stage we obtained two calibrated data cubes per night, one for each side of the IRDIS detector, resulting in ten data cubes. From each data cube we discarded bad frames by measuring the correlation of each frame with a reference frame that was tagged as good by visual inspection. 85% to 95% of the most correlated frames were kept for postprocessing, depending on the night. The night of December 7 was discarded due to its poor data quality, as already proposed by Apai et al. (2016).
We deliberately chose to skip the centering of the individual frames proposed by EsoRex. Instead, we used custom python routines, available in the VIP package (Gomez Gonzalez et al. 2016a,b), to precisely measure the position of the star and the related uncertainty for each individual frame of all data cubes by exploiting the four satellite spots. Indeed, since the satellite spots have a high S/N and are designed to be symmetric with respect to the star, one can use them to infer the position of the star. In practice, due to their wavelengthdependent elongation and to residual atmospheric dispersion, the satellite spots are not perfectly symmetric with respect to the star (Pathak et al. 2016). However, the symmetry is preserved at any given wavelength, and the spectrumweighted astrometric position of the four satellite spots remains symmetric with respect to spectrumweighted astrometric position of the star. To avoid the astrometric bias on the determination of the star position described by Pathak et al. (2016), the following strategy was adopted. For a given frame, we carefully fitted an asymmetric 2d Gaussian to each of the satellite spots to determine their respective centroid. Then, opposite centroids were connected by lines and the resulting intersection determined the estimated position (x,y) of the star in detector coordinates. This was performed for each frame to get the offset of the star from the center of the frame. For each data cube, a histogram of these offsets was built, and global offsets were obtained as the median of the vertical and horizontal offsets (see Fig. 1). All the frames were then shifted by the same amount for each cube to cancel the global offset, and cropped to a useful fieldofview of 511 × 511 pixels to reduce computation time during postprocessing. Our analysis suggests that a framebyframe recentering of the cube would not improve the final results, because the accuracy with which the stellar position can be determined in an individual frame is generally not smaller than the width of the histogram shown in Fig. 1. More details about the uncertainty on the position of the star are given in Sect. 3.4.
Fig. 1 Histogram of the horizontal and vertical offsets of the star with respect to the center of the frame in the December 5 (right) data cube. The vertical line represents the median of the histogram, and was used to globally recenter the data cube. The horizontal axis is in pixels, one pixel corresponding to 12.25 mas on sky. 

Open with DEXTER 
The parallactic angles corresponding to the individual frames of each data cube were independently calculated frame by frame. The MJD time at the middle of each frame was derived from the information given by the MJDOBS and HIERARCH ESO DET FRAM UTC header cards, which give the time at the start and the end of the observing sequence, respectively, by dividing the total integration time equally into 218 parts. The parallactic angles were calculated using the spherical trigonometry formula given in Meeus (1998) based on the equatorial coordinates precessed to the epoch of the observations and corrected for nutation, aberrations, and refraction.
2.3. Angular differential image processing
Fig. 2 Illustration of a fullframe PCA ADI postprocessed SPHERE/IRDIS image of HR8799 acquired with broadband H filter (left part) during the night of December 4, 2014. The central part was masked with a disk of radius 20 pixels. 

Open with DEXTER 
We carried out the data postprocessing with the opensource Vortex Imaging Processing^{1} package (VIP, Gomez Gonzalez et al. 2016a,b) written in Python 2.7. Our postprocessing is based on ADI techniques, which aim to reduce the quasistatic speckle noise and reveal the presence of offaxis sources by constructing and subtracting a reference onaxis pointspread function (PSF) from the individual frames of a data cube obtained in pupil tracking mode, where the star corresponds to the field rotation center (Marois et al. 2006; Lafrenière et al. 2007). Recently, Soummer et al. (2012) and Amara & Quanz (2012) proposed to take advantage of PCA to make ADI postprocessing more efficient. The PCAADI algorithm implemented in VIP is based on the approach presented in Amara & Quanz (2012), which can be summarized as follows:

construct a set of orthogonal reference images, the socalledprincipal components (PCs), using singular valuedecomposition (SVD, see e.g. Presset al. 2007) of the data cube;

project all the frames of the cube onto a truncated set of n_{PC} (<n_{frame}) PCs, where n_{frame} represents the total number of frames in a data cube;

reconstruct the data cube using a linear combination of PCs, and subtract the result from the original data cube to obtain a cube of residual frames;

rotate and collapse this cube of residuals to obtain the final image.
To optimize the determination of the astrometry, the S/N for each planet must be maximized. The S/N calculation implemented into VIP is based on a Student ttest (Student 1908) and follows the recommendation of Mawet et al. (2014) on small sample statistics (see Gomez Gonzalez et al. 2016b, for more details). The S/N of the planets in the final, postprocessed image depends mainly on the number of PCs used when building the reconstructed cube. A small number of PCs leads to an incomplete representation of the speckle noise, while a large number of PCs tends to capture the signal of the planet in the reconstructed cube, which leads to a lower algorithmic throughput for the planetary signal after subtraction. An optimum number of PCs can generally be found to maximize the planet S/N (Meshkat et al. 2014). For each data cube, we thus performed a grid search on the number of PCs to maximize the mean S/N in a region of one resolution element in diameter around each companion. The optimal n_{PC} for each data cube is reported in Table 1. Let us note that the PCA implemented in VIP comes with several SVD libraries, such as the efficient randomized SVD (Halko et al. 2011) and the wellknown LAPACK (see e.g. Anderson et al. 1990). We refer to Gomez Gonzalez et al. (2016a) for a complete discussion of all the SVD libraries available in VIP.
Final HR8799bcde astrometric measurements with respect to the host star for nights of December 4–6, and 8, 2014.
Figure 2 illustrates a VIP postprocessed image using fullframe PCA, where all the pixels of each frame are used at once to construct the reference images through SVD. The close region surrounding the host star is most affected by residual speckle noise and was masked with a disk of radius 20 pixels to better reveal the planets in Fig. 2. Throughout the present analysis, we also performed annuluswise PCA, which consists of performing PCA only for a thin annulus passing through a companion, with a typical width of a few resolution elements. Although fullframe PCA and annuluswise PCA may lead to slightly different results, this choice does not significantly affect the final astrometry, which is dominated by other sources of error (see Sect. 3). Furthermore, annuluswise PCA is significantly faster when performed on a single annulus, which is useful when dealing with large data cubes and/or when PCA is performed a large number of times (see Sect. 3.2).
3. Robust astrometry
The astrometric position of the HR8799bcde planets based on the December 2014 SPHERE/IRDIS data set has already been determined by Zurlo et al. (2016) and Apai et al. (2016). The Zurlo et al. (2016) final astrometry was obtained from the combination of four independent imageprocessing pipelines, by the quadratic sum of the error bar from each data reduction pipeline plus the standard deviation associated to the individual positions. Apai et al. (2016) used their implementation of the KLIP algorithm to derive the planets positions by injecting artificial planets with negative count rates, and used a manual inspection of the image quality and of the subtraction residuals to estimate the error bars. Here, we propose to go beyond these approaches and to study in detail the various contributions to the astrometric error budget, in an attempt to derive more reliable error bars. Our study is also meant to explore the ultimate astrometric accuracy of a stateoftheart instrument such as SPHERE, and to identify possible ways to improve the astrometric accuracy in future studies.
What we call robust astrometry consists of performing a proper evaluation of the statistical errors and systematic biases affecting the final astrometric estimation. The whole procedure consists of four steps: (i) the description and estimation of the instrumental calibration errors; (ii) the determination of the planets position with respect to the star and the related statistical error through Bayesian inference with MCMC sampling; (iii) the determination of the systematic error due to residual speckles, and (iv) the calculation of the error on the star position. Throughout the remainder of this section, we provide details for each step of the process.
3.1. Instrumental calibration and related errors
To derive accurate astrometric measurements from IRDIS images, various astrometric calibrations must be performed, namely the determination of the plate scale, the orientation of the north, and the optical distortion. Firstly, the plate scale, given in arcsecs per pixel, depends on the characteristics of all the optical elements composing the instrument. It allows conversion of positions given in pixels into arcsecs. Secondly, when observing in pupilstabilized mode, the vertical axis of the detector does not necessarily point towards north. Two contributions need to be taken into account: (i) the pupil offset, which accounts for the zero point position of the derotator and is assumed to be constant between runs; and (ii) the socalled true north, which accounts for a variation in the detector orientation with respect to the sky due to thermal or mechanical stresses, and which must be estimated during each observing run. Thirdly, the distortion in SPHERE/IRDIS is mainly dominated by an anamorphic magnification between the horizontal and vertical axes of the detector. This effect is due to the presence of toric mirrors in the common path of the instrument (see e.g., Zurlo et al. 2016).
Details of the observations used to derive those astrometric calibrations for IRDIS are described in Zurlo et al. (2016). We refer to that paper for the details, but we still provide the reader with the practical information used in this study. The astrometric calibrations were obtained from IRDIS observations of the globular cluster 47 Tuc acquired on December 15, 2014, with the same instrument setup and filter, and were compared to the Hubble Space Telescope data of the same field, precessed to the same epoch and corrected for the differential proper motions of the individual stars. The values derived by Zurlo et al. (2016) for the plate scale and true north based on this data set have recently been revised by the SPHERE consortium, using their improved knowledge of the instrument. This revised estimation, described in Maire et al. (2016), leads to a plate scale of 12.251 ± 0.009 mas/pixels and a true north orientation of . These values are valid for both the left and right parts of the IRDIS detector. The pupil offset, based on commissioning and guaranteedtime data obtained on several astrometric fields, is equal to . Finally, the IRDIS distortion measured on sky is dominated by an anamorphism of 0.60% ± 0.02% between the horizontal and vertical directions (Maire et al. 2016). Although the SPHERE calibration plan includes the daily measurement of distortion maps based on pinhole grids, we found that the quality of the astrometric estimations is not improved by using these maps. Prior to any postprocessing, we thus simply rescaled each frame of each cube by a factor 1.006 along the y axis. To take into account the uncertainty on this correction, an additional error of 0.02% on the plate scale will be considered in the following analysis.
3.2. Planet position and statistical error
The next step in the robust astrometry process consists of determining, for each data cube, the position of the planets with respect to the host star and estimating the statistical error related purely to the photon noise of the underlying thermal background and speckles through Bayesian inference based on
MCMC simulations.
This step does not describe the effect of the speckles themselves on the measured planet position, which will be discussed separately in Sect. 3.3. Our astrometric measurements are based on the negative fake companion technique (NEGFC, see e.g. Marois et al. 2010a; Lagrange et al. 2010), which consists of injecting a negative PSF template into the data cube with the aim of canceling out the companion (as well as possible) in the final postprocessed image based on a wellchosen merit function. The NEGFC technique is an iterative process, for which a step can be described as follows. For the chosen position/flux combination, a negative fake companion is injected into each frame of the data cube, and annularwise PCAADI processing is performed on a single annulus passing through the considered companion. The intensities I_{j} of N pixels are then extracted within a circular region with a radius equal to a few resolution elements, centered on a first guess position defined at the start of the iterative process (which means that the position of the circular aperture is fixed and does not change during the process). Assuming that the noise affecting the jth pixel value is given by (pure photon noise), we define the merit function as follows: (1)The position/flux of the NEGFC is then optimized to minimize the merit function in a threestep approach, as described in the following paragraphs. The resulting postprocessed images, before and after injection of a NEGFC, are represented in Fig. 3. Because no offaxis PSF was acquired in December 2014 with the same observing setup as for the HR8799 observations, the adopted PSF template corresponds to unsaturated offaxis images of β Pictoris obtained with SPHERE/IRDIS during science verification on January 30, 2015 (PI: A.M. Lagrange) with the same observing mode as for HR8799 (same coronagraph, same broadband H filter, and similar seeing ~1′′). The influence of this choice is discussed at the end of Sect. 3.2, together with a discussion of the effect of PSF chromatic dispersion on the measured planet position.
Fig. 3 Illustration of the annuluswise PCA postprocessing and merit function evaluation used in the negative fake companion technique. Left: no NEGFC was injected before annuluswise PCA processing. Right: a NEGFC was injected at the position and flux minimizing the merit function. The white circle illustrates the fixed circular aperture from which the pixel values I_{j} have been extracted to evaluate the merit function. The same colorscale was used for both images. 

Open with DEXTER 
First guess estimation.
With the optimal number of PCs in hand (see Table 1), we derive a first guess of the position of each companion in each PCAADI postprocessed image by simply identifying the highest pixel value in the close vicinity of the companion. We then derive a first guess of the flux of the companion by injecting a NEGFC at that position and by evaluating the merit function for a grid of possible fluxes. Only the flux is optimized during this stage, while the companion position is fixed to our first guess.
NelderMead optimization.
Although the first guess estimation results in a rough determination of the position/flux, and would constitute a valid initial set of parameters to start an MCMCbased Bayesian inference process (as presented in the next paragraph), it may turn out to be very time consuming due to the large number of merit function evaluations required to reach convergence in the MCMC and properly sample the posterior distributions. Thus, we propose to refine the first guess of the position/flux of the companions for the purpose of initializing the MCMC sampling close to the highly probable solution. To this end, we use the first position/flux estimation as an initial guess for a NelderMead simplexbased optimization (Nelder & Mead 1965) implemented into the SciPy Python library^{2}. The adopted merit function is defined in Eq. (1), and the position (r,θ) of the NEGFC is now allowed to vary during the fit in addition to its flux. As expected, this leads to a significant improvement of the position/flux determination. The right panel of Fig. 3 illustrates the result of an annuluswise PCAADI postprocessing performed on a single annulus passing through HR8799b, after injecting a NEGFC characterized by a position and flux minimizing the merit function.
The PCAADI algorithm that we first used in VIP relied on a randomized SVD library, which approximates the SVD of the data cube by using random projections and thereby provides increased computational efficiency (for details see Gomez Gonzalez et al. 2016a). Although this randomized approach is very efficient, the random process induces random variations in the merit function that can be significant compared to the variations of the merit function between two steps, especially when approaching the minimum. This can prevent the optimization process from reaching the true minimum of the merit function, or even from converging. Therefore, we decided to use the more classical, yet slower, deterministic SVD approach proposed in the LAPACK library for our PCAADI processing in the simplex optimization, as well as for the rest of this study. This choice is all the more important when the companion is located in a region dominated by residual speckle noise.
MCMC and final positions.
Fig. 4 Illustration of a typical corner plot obtained from the MCMC simulations using the NEGFC technique. The target companion is HR8799b observed during the night of December 6, 2014. The radial distance r (in pixels) and azimuth θ (in degrees) are detector coordinates with respect to the host star. The diagonal panels illustrate the posterior PDFs while those offaxis illustrate the correlation between them. 

Open with DEXTER 
Because the merit function used in the NelderMead optimization is not strictly convex, it is not guaranteed that the optimization will converge at the exact position of the planet, as it could potentially get stuck in a local minimum. Although this behaviour was generally not observed (as shown in Morzinski et al. 2015), we decided to use the NEGFC technique coupled with an MCMC approach to obtain the final flux and position of the HR8799 planets, expressed in polar coordinates, with respect to the host star. Let us recall briefly that the MCMC approach aims to sample the posterior probability density function (PDF), that is the probability of the position/flux parameters given the data cube and the prior knowledge (see e.g. Hogg et al. 2010). The VIP module dedicated to the NEGFC technique embeds the emcee package (ForemanMackey et al. 2013), which implements an affineinvariant ensemble sampler for MCMC proposed by Goodman & Weare (2010). Such an ensemble is composed of walkers, which can be considered as MetropolisHastings chains. The main difference between walkers and MetropolisHastings chains lies in the fact that the proposal distribution for a given walker depends, at a given step, on the position of all other walkers in the ensemble. Conversely, the proposal distributions involved in the MetropolisHastings algorithm are independent. Besides being more efficient in terms of the number of calls to the cost function, one major advantage of emcee is that it relies on only two calibration parameters, in comparison to the ~N^{2} parameters required for a MetropolisHastings algorithm in an Ndimensional parameter space to properly sample the PDF and speed up the process (for more details, see ForemanMackey et al. 2013; Goodman & Weare 2010, and references therein).
Fig. 5 Astrometry for HR8799bcde observed during the nights of December 4–6, and 8, 2014. The positions obtained from the left (resp. right) data cubes are represented with downward (resp. upward) black triangles. The error bars on the individual data points take into account all the contributions discussed in Sects. 3.1–3.4. The red dots correspond to the final astrometric measurements for each planet, together with the final error bar discussed in Sect. 3.5. The dashed lines represent the best orbital solutions for each planet in terms of reduced χ^{2} reported in Table 5. 

Open with DEXTER 
For each data cube and each companion, we carried out MCMC simulations to sample posterior PDFs related to the planet polar coordinates (r,θ) with respect to the host star and the planet flux f. For each MCMC simulation, we used 200 walkers firstly initialized in a small ball around the solution obtained from the NelderMead optimization. The chain was sufficiently close to convergence to allow Bayesian inference after, typically, 200 steps. More details concerning convergence statistical tests are given in Sect. 4.1, where we describe the PyAstrOFit Python package. In Fig. 4, the socalled corner plot illustrates the posterior PDFs and the correlation between the parameters (r,θ,f) for HR8799b observed during the night of December 6, 2014. Similar results were obtained for other planets and observing nights. Although a flux estimation for each planet is obtained, we focus our analysis on only the astrometry in this paper.
Taking into account the plate scale, the true north and pupil offset orientation (see Sect. 3.1), we have projected the HR8799bcde highly probable sets of polar coordinates onto the north and east directions. As a result, the eight HR8799bcde final positions for the four nights (left and right parts) are reported in Table 1 and displayed in Fig. 5. These positions will be used in Sect. 3.5 to deduce the final HR8799 astrometry for epoch 2014.93. In addition to obtaining the highly probable position/flux for a given companion, the MCMC simulations give a robust estimation of the statistical error on the astrometry (i.e., related purely to photon noise). This error, reported in Cols. 7 and 8 of Table 1, generally constitutes a minor contribution to the error budget, as discussed in the following sections.
Fig. 6 Speckle noise estimation for HR8799b observed on December 6, 2014. The histograms illustrate the offsets between the true position/flux of a fake companion and its position/flux obtained from the NEGFC technique. The dashed lines correspond to the 1D Gaussian fit from which we determine the speckle noise. 

Open with DEXTER 
Influence of the template PSF.
Since a nonsaturated, offaxis PSF was not obtained in the same observing mode during the nights where HR8799 was observed, we chose, as a PSF template for our NEGFC analysis, the closest offaxis PSF in time obtained with the same observing mode under similar weather conditions, which turned out to be an offaxis PSF of beta Pictoris obtained on January 30, 2015. The fact that both the instrument and the atmospheric conditions may have changed within the interval leads to a possible bias in our measurement of the planets’ position, which could vary from night to night. To evaluate this bias, we took a series of twelve offaxis PSFs observed in the same mode under good atmospheric conditions, obtained in 2015 in the context of the SHARDDS survey (J. Milli, priv. comm.). For each planet and each observing night in our HR8799 data set, we successively used the twelve offaxis PSFs as templates for the NEGFC technique, and derived the planets’ astrometry using the method described above. The dispersion of the astrometric measurements gives us an estimation of the bias that can be introduced by using a noncontemporaneous PSF. The observed dispersion does not depend significantly on the planet nor on the observing night, and has an overall standard deviation of 0.6 mas. This error bar will be added quadratically to the other error sources in Sect. 3.5.
Influence of residual dispersion.
Another source of imperfection in the recovery of the planets astrometry for broadband observations is the residual atmospheric dispersion after correction by the atmospheric dispersion correctors (ADC) included in the SPHERE optical path. While the small angular separation between the star and planets ensures the residual dispersion to be almost perfectly equal for all of them, their different spectra can result in a chromatic offset between their measured positions. The residual dispersion after correction by the SPHERE ADC has been shown to be smaller than 1.2 mas rms for zenith angles as large as the maximum of 54° encountered in the present data set (Hibon et al. 2016). Taking into account the Hband spectrum of the star and of the four planets (Bonnefoy et al. 2016), we estimate that the maximum astrometric offset between the star and planets due to residual dispersion cannot be larger than 0.25 mas in the worst case where residual dispersion shows a linear trend across the H band. This contribution is negligible in our final astrometric error budget.
3.3. Systematic error due to residual speckles
Performing PCAADI removes a large fraction of the quasistatic speckle noise and significantly improves the S/N of the companions. Although highly effective, this process is not perfect and some level of residual speckle noise remains in the postprocessed images. Such noise has a major impact on photometric and astrometric measurements (Guyon et al. 2012), and needs to be taken into account in the error budget. Since speckle noise is known to have a radial dependence, we propose to estimate its impact by injecting fake companions into the data cube at the radial distance of the real planets but for a wide range of angular positions, and by testing the ability of the NelderMead optimization to find their position and flux through the NEGFC technique. The first step in this process is to create an “empty” data cube by injecting four NEGFCs characterized by the highly probable positions/fluxes derived from the previous MCMC simulations. In the empty cube, we inject a fake companion characterized by a flux f_{true} and a radial distance r_{true}, both corresponding to the highly probable solution, but at an arbitrarily chosen angular coordinate θ_{true,i}. Using the NEGFC technique coupled with the NelderMead optimization, we determine the position/flux (r_{i},θ_{i},f_{i}) of the fake companion. We then compute the offsets Δr_{i} = r_{true}−r_{i}, Δθ_{i} = θ_{true,i}−θ_{i}, Δf_{i} = f_{true}−f_{i} between the known position/flux characterizing the fake companion and the solution obtained from the optimization process. The same process is repeated for a series of 360 azimuths equally spaced between 0° and 360°. These 360 realizations are used to build three normalized histograms, for Δr, Δθ and Δf, respectively. The histograms for Δr and Δθ are then fitted with a Gaussian function, and the standard deviations σ_{spec,r} and σ_{spec,θ} of the Gaussian functions are used as an estimation of the speckle noise affecting the radial and azimuthal coordinates. A similar approach was already used by Maire et al. (2015), for example. We illustrate in Fig. 6 the three histograms for HR8799b observed on December 6, 2014. The results obtained for all planets and data cubes are reported in Cols. 9 and 10 in Table 1. It appears clear that the error induced by speckle noise increases for decreasing angular separations of the companion with respect to the host star. Indeed, the brightness of the residual speckles increases closer to the star. We also note that speckle noise is consistently larger than statistical noise, except for HR8799b.
Another possible way to evaluate speckle noise is to measure the influence of the number of PCs used in the PCA postprocessing on the position/flux determination, as proposed by Pueyo et al. (2015), for example. Indeed, the residual speckle pattern changes as a function of the number of PCs. To verify the consistency of this method with the one proposed above, we determined the position/flux of each companion in each data cube using the NEGFC technique with the NelderMead optimization using a number of PCs ranging from 5 to 90 (for a number of PCs > 90, the companion selfsubtraction becomes too significant to get a high S/N). We then constructed three normalized histograms, for r, θ and f. As expected, the standard deviations of these histograms are similar to those deduced above.
Finally, we note that the residual speckle noise estimated here is in good agreement with the semiempirical estimation of the astrometric accuracy based on the planet S/N proposed in the case of pure photon noise by Guyon et al. (2012, Eq. (A1)), provided that we extrapolate this relation to the speckledominated regime in the following way, as already proposed by Mawet et al. (2015): σ_{1D} [λ/D] = 1/(πS/N). Using such a semiempirical formula therefore seems to be a possible method to acquire a rapid estimation of the astrometric error bar related to speckle noise, although we recommend going through the analysis presented in this section to obtain a robust estimation.
3.4. Error on the star position
Inside SPHERE, a dedicated differential tiptilt sensor is used to obtain an image of the PSF immediately upstream of the coronagraph, and is used as an input for closedloop control of the star position with respect to the coronagraph, thereby ensuring a stable star centering (Fusco et al. 2006; Baudoz et al. 2010). Based on laboratory measurements, the expected accuracy of the star centering is considered to be approximately 0.5 mas on sky (Baudoz et al. 2010). As mentioned in Sect. 2.2, no individual frame centering was applied to the data cubes in our analysis, but rather a global centering of all frames in each individual cube using the same x,y offsets.
Here, we independently estimate the uncertainty on the mean star position for each data cube. The evaluation of this uncertainty is based on the histogram of the x and y offsets measured for all individual frames by the centroid plus intersection method described in Sect. 2.2. The mean position of the star in a given data cube can be obtained by a Gaussian fit of the two histograms, as illustrated in Fig. 1. Based on this figure, in the following discussion we assume that the histograms follow a Gaussian distribution, so that the accuracy on the determination of the mean stellar position in a given cube is given by the standard deviation of the bestfit Gaussian divided by the square root of the number of realizations. The standard deviations of the bestfit Gaussians are given in Table 2 in terms of right ascension (RA) and declination (Dec), by projecting the (σ_{⋆,x},σ_{⋆,y})error ellipses expressed in detector coordinates onto the north and east directions. We note that the derived stellar jitter estimation is slightly larger than predicted in Baudoz et al. (2010), with values varying from 0.76 mas to 1.44 mas depending on the night (i.e., around 0.1 pixel in detector coordinates). Based on these values, and taking into account the ~200 frames present in each data cube, the error bar on the mean stellar position in any given cube amounts to less than 0.1 mas, and is therefore completely negligible in our final noise budget.
Estimation of the stellar jitter in the eight data cubes.
However, this contribution represents only the purely statistical error on the determination of the star position. We also need to take into account possible systematic biases on the determination of the star position based on the satellite spots. To this end, we obtained a data set on a relatively bright star, using the waffle mode of the DM, but without coronagraph. The star was mildly saturated at its center to increase the S/N on the satellite spots. We determined the center of the star based on a truncated Moffat profile to reject the saturated part of the PSF, and compared this estimation with the prediction based on the satellite spots. We verified that the two estimations match with an accuracy better than 0.1 pixels, which represents our best estimation of an upper limit on a possible bias. This also confirms that the method proposed in Sect. 2.2, to determine the stellar position from the satellite spots, does not lead to major astrometric bias, even in the presence of residual atmospheric dispersion. Here, we conservatively assume that a bias of 0.1 pixels (1.2 mas) affects our determination of the mean star position in all cubes.
3.5. Final astrometry
Particular care must be taken when combining the results and error bars of several astrometric measurements, especially in the presence of correlated errors. How the various error bars add up requires specific discussion. Firstly, we note that our experimental determination of the error bar related to residual speckles inherently takes into account the contribution of photon noise. Indeed, the empirical intensity of the speckles includes the contribution of the photon noise associated to all sources of signal at any given location (stellar residuals, planet, sky emission and thermal background). This is supported by the fact that the error bar associated to speckle noise generally dominates the error bar associated to photon noise. Only in the case of planet b are they of the same order of magnitude, which reflects the fact that residual speckles are very faint compared to residual background noise at that angular distance from the star.
Secondly, we make the conservative assumption that the errors related to speckle noise are fully correlated, not only between the left and right data cubes obtained on the same night, but also between all nights. The assumption of full correlation between the left and right data cubes is justified by the fact that the signals recorded by the two parts of the detector are almost identical (to within photon noise and some minor differential aberrations that amount to a few nm rms at most), and is backed up by the fact that the estimated error bars are almost identical for the left and right sides for most of the nights and planets (see Table 1). The assumption that speckle noise is fully correlated from night to night is more debatable. It is indeed expected that speckle noise will be partly correlated between successive nights, because residual speckles are often associated to noncommon path aberrations in the instrument that can vary on very long timescales. To remain on the conservative side, we will assume a full correlation of speckle noise in all data sets. The error bar on the final astrometry regarding speckle noise should then be computed as the median of all speckle noiserelated error bars. We note however that the estimations of the speckle noiserelated error bars significantly vary from one night to another (see Table 1), which suggests that this noise is at least partly uncorrelated, and that our final error bars will be pessimistic.
The final HR8799bcde astrometric measurements with respect to the star for epoch 2014.93.
Thirdly, we proposed in the previous section that the final error bar related to the determination of the star position is dominated by a systematic bias that can amount to 1.2 mas, and that the variability of the PSF shape can induce a bias of up to 0.6 mas. These biases will be added quadratically to our final astrometric error bar for all planets. The same applies to instrumental calibration errors, which are supposed to affect all data cubes in the exact same way. Indeed, appropriate observations of astrometric fields were not performed on each of the five HR8799 observing nights. We therefore had to rely on an astrometric calibration carried out by the SPHERE consortium one week later (see Sect. 3.1), which was used as a reference for all five nights. Although we could not check the stability of the calibration over several nights, we note that the latest IRDIS astrometric calibrations by the SPHERE consortium show that the time variations of plate scale and true north are mostly within their estimated error bars, based on two years of astrometric field observations, while the pupil offset and anamorphic factor are mostly constant (Maire et al. 2016). This suggests that our final estimation of the astrometric error bar should not include any unaccounted bias related to the variability of the IRDIS astrometric calibration. That being said, we still recommend that, in future observing programs dedicated to precise astrometric measurements, observations of standard astrometric fields be obtained during each individual night to ensure a high astrometric robustness.
Based on these assumptions, the computation of the final astrometry and related error bars proceeds as follows for each planet:

define the final astrometry of the four planets as the weightedmean of the eight individual positions (left and right parts of thedetector for the four nights), using the inverse of the variance ofspeckle noise as a weight;

estimate the final error bar related to speckle noise as the median of the individual error bars on the eight astrometric measurements;

add quadratically; the contribution of speckle noise, the upper limit on the stellar centering bias, and the contribution of instrumental calibration errors to obtain the final astrometric error bars.
All these calculations are performed in polar coordinates, reflecting the fact that error bars generally have different behaviors along the radial and azimuthal directions. The last step is based on the following formulae: where r is the radial distance in pixels, σ_{r,spec} and σ_{θ,spec} the final radial (pixels) and azimuthal (degrees) error bars related to speckle noise, σ_{r,⋆} and σ_{θ,⋆} the radial (pixels) and azimuthal (degrees) stellar centering biases, σ_{r,PSF} and σ_{θ,PSF} the radial (pixels) and azimuthal (degrees) error bars related to the imperfection of the PSF template in the NEGFC analysis, σ_{r,AF} and σ_{θ,AF} the radial and azimuthal errors on the anamorphic factor expressed as percentages, and where PLSC refers to the plate scale in ″/pixel, PO to the pupil offset and TN to the true north, both in degrees. The final astrometries and related error bars are given for the four planets in Table 3 and are illustrated in Fig. 5. Table 3 includes a projection of the error bars onto the RA and Dec directions, to comply with the usage. However, we suggest that expressing the error bars in polar coordinates is more appropriate, because polar coordinates usually correspond to the major and minor axes of the error ellipse. Another, even more appropriate way to proceed would be to specify the error ellipse by its three parameters (two axes and position angle). In the present case, the error bars are sufficiently symmetric to proceed with RA/Dec error bars, even though we note that the HR8799b error bars are significantly asymmetric, the angular error bar being twice as large as the radial one. This is mostly due to the large uncertainty on the pupil offset (, see Sect. 3.1), which severely affects planets located far from the star.
To verify the consistency of our error bars, we compared the statistical distribution of the eight individual data points obtained for each planet to the individual error bars on the eight data points. Table 4 shows that the final error bars are generally approximately twice larger than the dispersion of the individual data points. This is related to the fact that the major error sources (speckle noise, stellar position bias, and instrumental calibration) are supposed to be fully correlated between individual measurements, so that the final error bar has a similar size as the individual ones. This suggests that an improvement by up to a factor two in astrometric accuracy could be achieved by improving the astrometric calibration. That said, the individual error bars are in relatively good adequacy with the dispersion of the data points (see Fig. 5), although we note a significant asymmetry in the distribution of the data points towards the NESW direction. This asymmetry appears relatively consistent between the four planets, and we therefore suggest that it comes from a time variability in the bias on the stellar position measurement (the only error source that is naturally expressed in RA/Dec), which could be related to variations in the PSF shape and/or in the diffraction pattern created by the DM on a nighttonight timescale. This variation remains within the expected amplitude of approximately 0.1 pixels for the star position bias.
For planet b, the main contribution to the error budget comes from the imperfect astrometric calibration and from the uncertainty on the star position, while speckle noise is negligible. This is consistent with the fact that HR8799b lies in a region that is not significantly affected by residual speckles (see Fig. 2). For planet c, although speckle noise significantly increases, the noise budget remains dominated by the astrometric calibration and stellar position uncertainties. The dominance of stellar centering noise in the astrometric error budget of these two planets is supported by the fact that the dispersion in the individual astrometric measurements for planets b and c has a similar amplitude and shape (see Fig. 5), as is expected for a global centering error. For the two innermost planets (d and e), speckle noise progressively becomes the dominant contributor to the error budget, and once again this is consistent with Fig. 5, where the dispersion of the astrometric data points increases significantly, especially for planet e. We finally note that our astrometric measurements are in general agreement with the astrometric measurements derived in Zurlo et al. (2016) and Apai et al. (2016) to within error bars, but that our error bars are two to three times smaller, thanks to a careful evaluation of all systematic error sources. For the orbital architecture analysis presented in the next section, we will thus only use our data reduction for the IRDIS data set of December 2014.
4. Orbital fitting analysis
4.1. The PyAstrOFit Python package
To perform our analysis of the HR8799bcde orbital architecture, we adopted the Bayesian framework. With the aim of making our results reproducible as well as allowing anyone to straightforwardly perform similar analyses, we introduce the PyAstrOFit package^{3} implemented in Python 2.7, which is fully dedicated to orbital fitting using the MCMC approach. The code is open source and has been used to carry out our analysis and to produce all the figures presented in this section. PyAstrOFit is composed of several modules, but the core of the package relies on three main modules, referred to as Orbit, Sampler and Inference:

The Orbit module is used to instantiate an orbit object, whichincludes the required data to model or represent any bound orbit.Unbound orbits are not considered because they would requirethe use of universal Keplerian variables and Stumpff functions(Beust et al. 2016).

The Sampler module constitutes the core of the Markov chains construction. It embeds the emcee package (ForemanMackey et al. 2013), which implements the affineinvariant ensemble sampler for MCMC proposed by Goodman & Weare (2010).

The Inference module is dedicated to Bayesian inference. Its main purpose is to represent both the posterior PDFs and the correlations between parameters from the Markov chain, to determine the confidence intervals, and to derive a set of allowable orbits or the best solution in terms of reduced χ^{2}.
The PyAstrOFit sampler comes with several convergence diagnostic tools. In practice, one can never be sure that a chain has actually converged, but there exists several tests to evaluate whether the chain appears to be close to convergence (or more precisely, far from nonconvergence^{4}):

The acceptance rate (MacKay 2003), whichcorresponds to the fraction of accepted to proposed candidates,can be monitored: if the acceptance rate is too high, the chain isprobably not mixing well, while a low acceptance rate indicatesthat too many proposed candidates are rejected (which issymptomatic of a walker stuck in a given position).

The GelmanRubin statistical test (Gelman & Rubin 1992; Ford 2006; Gelman et al. 2014) compares, for each parameter, the variance estimated from nonoverlapping parts of the chain to the variance of their estimates of the mean. A large value may arise from slow chain mixing or multimodality (Cowles & Carlin 1996). Conversely, a value close to 1 indicates that the Markov chain is close to convergence.

The lag ρ_{k} autocorrelation corresponds to the correlation between every draw and its kth lag. A relatively high ρ_{k = K} value for a given K indicates a high degree of correlation between the draws, a slow mixing and a chain far from convergence.

The integrated autocorrelation time τ (see e.g. ForemanMackey et al. 2013; Christen & Fox 2010; Goodman & Weare 2010), also called inefficiency factor, aims to give an estimate of the number of posterior PDF evaluations required to draw an independent sample. The smaller τ, the better.
Following Ford (2006) and Chauvin et al. (2012), the sampling can be done on three different state vectors noted x = (a,e,i,ω,Ω,t_{p}), x′ = (log P,e,cosi,Ω + ω,Ω−ω,t_{p}) where P represents the orbital period, and u(x) defined at Eq. (A.1) in Chauvin et al. (2012). As suggested by Ford (2006), adopting a uniform prior distribution for x′ may help to improve the convergence of the chain. Indeed, a uniform prior distribution for cosi in the interval [− 1,1] implies a prior distribution proportional to sini in the interval [− 90°,90°] for the inclination. As a consequence, it implies that orbits characterized by i ≃ 0° (faceon) are considered intrinsically less probable than those characterized by i ≃ 90° (edgeon).
Which statistical test one should adopt as a convergence criterion is a question with no trivial solution. Nonetheless, some recommendations can be found in the literature (see e.g. Cowles & Carlin 1996). For instance, Gelman et al. (2014) theoretically demonstrated that the optimal Markov chain mixing to sample normally distributed posterior PDFs is characterized by an acceptance rate equal to ~0.44 when adopting the MetropolisHastings algorithm. It is widely agreed that an acceptance rate between 0.2 and 0.5 constitutes an appropriate value to ensure a good Markov chain mixing. The criteria adopted in this analysis are defined in the following section where we address the HR8799 orbit fitting, and are generally similar to those used by Pueyo et al. (2015) in their study of the HR8799 orbital parameters. When all the statistical tests meet the criteria, we consider that the part of the chain satisfying them has converged. The Inference module can then be used to draw independent samples from the chain to construct the final posterior PDFs for each Keplerian parameter. The confidence intervals are defined in terms of highest density regions (HDR, Hyndman 1996), also referred to as highest posterior density intervals. For a given confidence level 1−α (0 ≤ α ≤ 1), the idea is to take a horizontal line and shift it up until the area under the regions of the PDF located above this line represents a fraction 1−α of the total area under the PDF. The projection to the x axis of this area defines the 100(1−α)% HDR. To infer our confidence intervals, we chose to use the 68.3% HDRs. In the ideal case where the PDF f is unimodal, the HDR corresponds to the smallest of all intervals [a,b] that satisfy Pr(a ≤ x ≤ b) = 1−α, which happens such that f(a) = f(b).
The Inference module also comes with various tools to display the results, such as corner plots to illustrate the PDFs and their corresponding correlation, walk plots to illustrate the mixing of the chain, and the illustration of allowable orbits together with the data. More information about all the PyAstrOFit possibilities and tutorials dedicated to each module can be found in the GitHub repository.
4.2. HR8799 orbital fitting with PyAstrOFit
Here, we revisit the MCMCbased Bayesian analysis described in Pueyo et al. (2015) using more robust convergence criteria before using our chains for inference, and using an extended data set by adding not only the SPHERE astrometric data presented in this work but also the latest astrometric measurements from the literature (see below). We update those results to what is described in this section. As a consequence, this analysis supersedes that presented in Pueyo et al. (2015).
For our orbital analysis, we assumed the distance of the system and the mass of the host star to be 39.4 pc (van Leeuwen 2007) and 1.51 M_{⊙} (Baines et al. 2012), respectively. The possibility of inferring these parameters from the orbital fitting modules has not been implemented into PyAstrOFit yet, but could be the subject of a future update. The astrometric positions used for the orbit fitting come from the works of Marois et al. (2008), Lafrenière et al. (2009), Fukagawa et al. (2009), Metchev et al. (2009), Hinz et al. (2010), Currie et al. (2011, 2012, 2014), Bergfors et al. (2011), Galicher et al. (2011), Soummer et al. (2011), Esposito et al. (2013), Maire et al. (2015), Zurlo et al. (2016), and Konopacky et al. (2016) along with the positions derived in the present work. A compilation of all the astrometric measurements for HR8799bcde is provided in Table A.1. All these data are included in the PyAstrOFit source code, and can easily be queried by an interested user.
Fig. 7 Results of the MCMC simulations for HR8799b, displayed as a corner plot for the Keplerian elements a, e, i, Ω and ω. The diagonal panels illustrate the posterior PDFs while the offaxis panels illustrate the correlation between the parameters. The yellow lines and crosses correspond to the best solution in terms of reduced χ^{2}. 

Open with DEXTER 
Fig. 8 Same as Fig. 7, but for HR8799c. 

Open with DEXTER 
Fig. 9 Same as Fig. 7, but for HR8799d. 

Open with DEXTER 
Fig. 10 Same as Fig. 7, but for HR8799e. 

Open with DEXTER 
The presence of systematic errors, whose careful calibration is described in Sect. 3, turned out to be a significant nuisance when trying to reconcile contemporaneous astrometric measurements from various instruments. For instance, the difference in astrometry for HR8799b between Currie et al. (2014) and Pueyo et al. (2015) could be explained by an offset in the true north position between Palomar and Keck (this offset has a lesser impact for the planets at smaller separations). Given the relative paucity of data in the astrometric calibrator (based on a single astrometric binary) presented in Pueyo et al. (2015), compared to the long history of high precision Keck astrometry (see e.g., Yelda et al. 2010; Service et al. 2016), we chose to include only the Currie et al. (2014) points in our analysis. Another example would be the discrepancy for the position of HR8799d between these two papers, which can be easily traced back to the presence of a bright residual speckle near planet d in the Palomar/P1640 data. Because there is no CPUefficient method to carry out the negative injection in IFS data without ADI, Pueyo et al. (2015) could not use the method described in Sect. 3.1, and only used a method similar to that presented in Sect. 3.3 (at other azimuth angles, where there was no bright speckle). As a consequence the Pueyo et al. (2015) uncertainties on HR8799d are most likely underreported, and we instead use the contemporaneous estimate of Currie et al. (2014). These two examples illustrate the complexity of precision astrometry in high contrast imaging and the importance of carrying out all the steps of robust astrometry as described in the present paper. Because the HR8799 system has been observed by multiple instruments since 2012, and because we had access to the P1640 data, we could conduct these instrument to instrument sanity checks and choose the most robust published astrometry. Before 2012, the measurements are more sparse and we thus decided to include all of them in the orbit fitting in the absence of further information.
The affine invariant sampler implemented in emcee comes with only two hyperparameters to be tuned: the number of walkers and an adjustable scale parameter a> 1, which has a direct impact on the acceptance rate of each walker (see Goodman & Weare 2010). All our simulations were performed with 1200 walkers. For a given Keplerian parameter, the set composed of the first element of each walker constitutes the initial distribution, which depends on how we decide to initialize the walkers. The equilibrium distribution, which we expect to be close to the posterior distribution when the chain has converged, should not depend on this initial distribution (see e.g. Meyn & Tweedie 2009). As discussed in ForemanMackey et al. (2013), starting the simulation with an initial distribution close to the expected posterior distribution speeds up the convergence. This can be done by initializing the walkers in a small Ndimensional ball in the parameter space around a highly probable solution. However, such an approach not only requires a priori knowledge of the main posterior distribution peak, but can also jeopardize the chain convergence if the posterior distribution is multimodal. Various alternatives are proposed in ForemanMackey et al. (2013), and we have tested some of them. In particular, we can start the walkers uniformly over a given range in the parameter space. All our tests have led to identical results for all the Keplerian elements of each planet.
We have started the chain construction with a minimum of 1000 steps per walker before beginning any convergence tests. During the MCMC run, the acceptance rate was monitored and the hyperparameter a (initially set to 2) was dynamically tuned in order to ensure an acceptance rate between 0.2 and 0.5 for at least 75% of the walkers. The walkers for which acceptance rate was outside [0.2,0.5] were discarded and not used for Bayesian inference. We considered that a chain has converged when the GelmanRubin statistical test is satisfied three times in a row for all Keplerian parameters. The convergence was reached after typically 40 000 steps per walker, for a total computing time equal to three hours using a computer equipped with 28 processing units.
Confidence intervals (1σ) and best solutions in terms of for all the Keplerian elements, as well as for the orbital period P.
For each planet, the corner plots are depicted in Figs. 7 to 10. They illustrate the resulting posterior PDFs and the correlation between the Keplerian parameters a, e, i, Ω and ω. The yellow lines represent the best solution in terms of reduced χ^{2}. These solutions do not necessarily coincide with the peaks of the posterior PDFs. Following the procedure described in Sect. 4.1, the confidence intervals were inferred from the posterior PDFs and are summarized in Table 5 together with the best solution in terms of reduced χ^{2}. In addition, Figs. B.1 and B.2 illustrate a set of 1000 allowable orbits characterized by and for which all the Keplerian parameters are in the confidence intervals reported in Table 5.
4.3. Discussion
One of the most striking results of our MCMC analysis concerns the eccentricity of planet d, which shows a clear peak around e_{d,peak} ~ 0.35, and rejects the circular orbit hypothesis outside its 1σ confidence interval (although a significant set of solutions characterized by e_{d}< 0.2 cannot be ruled out). Actually, all four planets show possible signs of noncircular orbits at various levels, as the eccentricity posterior PDFs are generally rather broad, or in some cases (planets d and e) not monotonically decreasing.
Fig. 11 Coplanarity test for HR8799bcde. The scatter plot illustrates the vector coordinates on the sky plane (see text) for the HR8799b (blue), HR8799c (red), HR8799d (orange) and HR8799e (green) allowable orbits. All the points on a given circle refer to orbits characterized by the same inclination. Similarly, all points on a given radial branch refer to orbits characterized by the same longitude of the ascending node. The angle between two spokes corresponds to 20°. The solid lines represent the 68% confidence interval around the most probable vector coordinates, which are depicted by thick points. 

Open with DEXTER 
Another striking result concerns the orientation of the orbits. The inclination of all the planets is similar and lies between approximately 20° and 38°. It clearly rules out a family of solutions previously proposed in the literature (see e.g. Marois et al. 2010b; Currie et al. 2011), for which the planets have a faceon orbit (i = 0°). The longitude of the ascending node Ω generally shows a large confidence interval for all planets, but we can readily note a significant difference (at the >2σ level) between the Ω derived for planets b and c. Comparing the individual inclinations and longitudes of ascending nodes has, however, a limited usefulness, and we therefore propose to compare the threedimensional relative orientations of the orbits for all four planets, in order to test the coplanarity of the system. This can be done by projecting onto the sky plane the normalized vector orthogonal to the orbital plane, defined by: (4)The scatter plot represented in Fig. 11 illustrates the vector coordinates on the sky plane for the allowable orbits of the four planets. The pole of the polar grid locates the projected vector that points towards Earth. All the points on a given arc of a circle refer to orbits characterized by the same inclination. Similarly, all points on a given spoke refer to orbits characterized by the same ascending node longitude. Planets d and e have very wide distributions of orientations, which are compatible with any other individual planet. We can even note that the 68% confidence interval is disjointed for planet d, which echoes the bimodal PDFs seen in Fig. 9. However, there is a clear discrepancy between the orbital planes of planets b and c, for which the orbital planes show a mutual inclination of 35°, and the 68% confidence intervals are largely disjointed. This suggests that the system might not be coplanar, at a significance level of approximately 2σ. Taken together with the evidence for nonzero eccentricity for planet d, this represents new empirical constraints that may be hard to reconcile with the meanmotion resonance scenarios currently proposed in the literature (Fabrycky & MurrayClay 2010; Goździewski & Migaszewski 2014). Longlived, nonresonant orbital architectures do not seem to predict these peculiar features either (Gotberg et al. 2016).
A thorough analysis of the system dynamics and stability, taking into account the most recent positions, would be of great interest but is beyond the scope of this paper. Yet, an interesting clue in this context is the compatibility of the planet periods with meanmotion resonances (Fabrycky & MurrayClay 2010; Goździewski & Migaszewski 2014). We therefore propose to simply identify meanmotion resonances compatible with our results. To this aim, we illustrate in Fig. 12 the distribution of the ratios of periods, P_{b}/P_{c}, P_{c}/P_{d} and P_{d}/P_{e}, obtained by dividing the respective PDFs. These PDFs were obtained by applying the Kepler’s third law to the semimajor axis PDFs illustrated in Figs. 7 to 10, assuming a mass of 1.51 M_{⊙} (Baines et al. 2012) for the starplanets system. The highly probable period ratios and the associated confidence intervals are , and , respectively. It appears clear that the most compatible meanmotion resonance between planets b and c is 1c:2b. For the other two planets, the most probable meanmotion resonances are 2d:5c and 2e:3d, while the 1c:2d and 1d:2e solutions also have a high probability. In Table 6, we derive the probability of various meanmotion resonances by determining the fraction of orbits with ratios of periods in the range of 10% around the selected resonance, that is, comprised in the [1.9,2.1] interval for the 1:2 resonance, as proposed in Currie et al. (2012). The 1e:2d:4c:8b configuration has been identified in the literature as a stable resonant configuration for the system (see e.g. Soummer et al. 2011; Currie et al. 2012; Goździewski & Migaszewski 2014). The probabilities reported in Table 6 show that this configuration is likely, with probabilities equal to 16.25%, 12.07% and 13.52%, respectively, for the b–c, c–d and d–e pairs.
Comparing our results with other works is not an easy task, due to various assumptions made in the literature. We can still note that the confidence intervals reported in Table 5 include a significant fraction of the solutions derived in Goździewski & Migaszewski (2014), Currie et al. (2012), Maire et al. (2015) and Zurlo et al. (2016). In the latter two papers, only circular orbits characterized by welldefined meanmotion resonances are considered for the orbit fitting. This may explain the slight discrepancies between the results obtained from the two different approaches. In particular, the semimajor axis a derived in Maire et al. (2015) and Zurlo et al. (2016) corresponds systematically to the upper boundary of our confidence intervals, which is linked to almost circular orbits. The only works that can be directly compared with ours are the studies of Pueyo et al. (2015) and Konopacky et al. (2016), which include PDFs for all the orbital elements resulting from an MCMC posterior sampling or from Monte Carlo simulations. We refrain from giving a direct comparison of our study with the work of Pueyo et al. (2015) due to the discrepancies identified in Sect. 4.2, and briefly discuss the compatibility of our results with those of Konopacky et al. (2016).
Comparing our PDFs with those presented by Konopacky et al. (2016) shows a broad consistency between the two analyses, although their constraints on the orbital parameters are generally broader than ours, especially regarding the semimajor axis and inclination. This is most probably due to the shorter time baseline used in their analysis. Some discrepancies can be noted, however. One of them concerns the argument of the periastron for the four planets, which are generally not consistent. This is not unexpected since this orbital parameter is particularly difficult to constrain based on such small phase coverage for orbits with low eccentricity. Besides the argument of periastron, we can note some intriguing differences. The main one relates to coplanarity: our MCMC analysis does not support their conclusion that the system is most probably coplanar. More precisely, Konopacky et al. (2016) do not favour the solutions for planet c with a longitude of ascending node Ω_{c} in the range [125, 153] deg derived in our study, while a broad peak can also be noticed at approximately Ω_{c} = 130 deg in their PDF. They argue that a higher number of loweccentricity solutions does favor an Ω_{c} near ~50 deg. This is at odds with our analysis, as we do not find any peak at approximately 50 deg in the PDF of Ω_{c}, while the vast majority of our solutions for planet c have an eccentricity lower than 0.2. Another difference concerns the inclination found for the orbital planes of planets d and e, which peaks around 45 deg in their analysis, while our results suggest approximately 30 deg for both planets. The discrepancy between these results reaches a significance level of approximately 2σ. It is difficult to assess whether these discrepancies could be (partly) related to instrumentspecific biases and/or to the different approaches used to sample the PDFs, or if they only result from the increased time baseline and enlarged data set in our study.
5. Conclusions
In the first part of this work, we present an independent data reduction of SPHERE/IRDIS images of the planetary system HR8799, acquired in the broadband H filter using coronagraphic imaging during the December 2014 science verification run. To achieve a robust determination of the astrometric position of the planets with respect to the star, we took advantage of the angular differential imaging postprocessing algorithm implemented in the VIP pipeline (Gomez Gonzalez et al. 2016a,b), based on principal component analysis, and performed a detailed analysis of the various contributions to the astrometric error budget. The resulting astrometric positions agree within 1σ with previous estimations based on the same data set (Zurlo et al. 2016; Apai et al. 2016), with error bars two to three times smaller thanks to a careful estimation of systematic errors. The main contribution to the astrometric error depends on the angular distance from the star: the error budget is dominated by the uncertainty on the stellar position (~1 mas) and instrumental calibration errors for planet b, while residual speckle noise increases for smaller angular separations and becomes dominant for planet e. We note that these revised error bars match the early expectations of SPHERE in terms of astrometric accuracy (~2 mas), and suggest that the astrometric accuracy could even be further improved (especially for planets located outside the specklenoise dominated regime) by a more careful IRDIS astrometric calibration and by improving upon our estimation of the bias on the star center determination using dedicated calibration programs. In practice, nothing seems to prevent SPHERE/IRDIS from reaching a 1 mas astrometric accuracy in the future, based on a careful calibration plan.
Fig. 12 Histograms of the period ratios P_{b}/P_{c} (left panel), P_{c}/P_{d} (middle panel) and P_{d}/P_{e} (right panel) derived from the semimajor axis posterior PDFs given in Figs. 7 to 10 and adopting 1.51 M_{⊙} for the host star mass. These period ratios allow to identify which meanmotion resonances are compatible with the allowable orbits derived from our MCMC analysis. It appears clear that the most compatible meanmotion resonance is 1e:2d:4c:8b. 

Open with DEXTER 
In the second part of the paper, we presented the opensource PyAstrOFit package written in Python 2.7, fully dedicated to orbital fitting within the Bayesian framework. Thanks to PyAstrOFit, we performed the orbital motion analysis and derived posterior PDFs for the six Keplerian elements of each planet. While planets b, c and e are characterized by small eccentricities, the eccentricity of planet d clearly peaks at e_{d} = 0.35, yet without ruling out solutions with smaller eccentricities. The combination of the posterior PDFs for the inclination and the longitude of the ascending node allowed us to evaluate the distribution of the projection on the sky of the vector perpendicular to the highly probable orbits. Our results seem to rule out the coplanarity between the four planets, especially for planets b and c, whose orbital planes are shown to be inclined by 35° with respect to each other. If confirmed by future astrometric monitoring, the significant eccentricity of planet d and the noncoplanarity of planets b and c would become stringent constraints for dynamical models of the HR8799 planetary system. We note in particular that an eccentricity larger than 0.3 is generally not compatible with the predictions based on systems locked in multiple meanmotion resonances (Goździewski & Migaszewski 2014), nor with longlived nonresonant systems (Gotberg et al. 2016). New dynamical simulations based on the latest astrometric measurements will be particularly useful in constraining the origin and fate of this system. We also recommend future studies to give up on the circular and coplanar assumptions, which are not backed up by our detailed MCMC simulations, although we should not overlook the fact that orbital fitting can be strongly affected by unaccounted astrometric biases or underestimation of astrometric errors.
Relative probabilities (in %) associated with different meanmotion resonances between consecutive pairs of planet.
Acknowledgments
The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (ERC Grant Agreement No. 337569), and from the French Community of Belgium through an ARC grant for Concerted Research Action.
References
 Amara, A., & Quanz, S. P. 2012, MNRAS, 427, 948 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, E., Bai, Z., Bischof, C., et al. 1990, Comput. Sci. Dept. Technical Report, CS90105 [Google Scholar]
 Apai, D., Kasper, M., Skemer, A., et al. 2016, ApJ, 820, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Baines, E. K., White, R. J., Huber, D., et al. 2012, ApJ, 761, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Baudoz, P., Dorn, R. J., Lizon, J.L., et al. 2010, in Groundbased and Airborne Instrumentation for Astronomy III, Proc. SPIE, 7735, 77355 [CrossRef] [Google Scholar]
 Bergfors, C., Brandner, W., Janson, M., Köhler, R., & Henning, T. 2011, A&A, 528, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beust, H., Bonnefoy, M., Maire, A.L., et al. 2016, A&A, 587, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beuzit, J.L., Feldt, M., Dohlen, K., et al. 2008, in Groundbased and Airborne Instrumentation for Astronomy II, Proc. SPIE, 7014, 701418 [CrossRef] [Google Scholar]
 Bonnefoy, M., Zurlo, A., Baudino, J. L., et al. 2016, A&A, 587, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Booth, M., Jordán, A., Casassus, S., et al. 2016, MNRAS, 460, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Carbillet, M., Bendjoya, P., Abe, L., et al. 2011, Exp. Astron., 30, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Chauvin, G., Lagrange, A.M., Beust, H., et al. 2012, A&A, 542, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Christen, J. A., & Fox, C. 2010, Bayesian Anal., 5, 263 [CrossRef] [Google Scholar]
 Cowles, M. K., & Carlin, B. P. 1996, J. Am. Statist. Assoc., 91, 883 [CrossRef] [Google Scholar]
 Currie, T., Burrows, A., Itoh, Y., et al. 2011, ApJ, 729, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Currie, T., Fukagawa, M., Thalmann, C., Matsumura, S., & Plavchan, P. 2012, ApJ, 755, L34 [NASA ADS] [CrossRef] [Google Scholar]
 Currie, T., Burrows, A., Girard, J. H., et al. 2014, ApJ, 795, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Dohlen, K., Langlois, M., Saisse, M., et al. 2008, in Groundbased and Airborne Instrumentation for Astronomy II, Proc. SPIE, 7014, 70143 [CrossRef] [Google Scholar]
 Esposito, S., Mesa, D., Skemer, A., et al. 2013, A&A, 549, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabrycky, D. C., & MurrayClay, R. A. 2010, ApJ, 710, 1408 [NASA ADS] [CrossRef] [Google Scholar]
 Ford, E. B. 2005, AJ, 129, 1706 [NASA ADS] [CrossRef] [Google Scholar]
 Ford, E. B. 2006, ApJ, 642, 505 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Fukagawa, M., Itoh, Y., Tamura, M., et al. 2009, ApJ, 696, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Fusco, T., Rousset, G., Sauvage, J.F., et al. 2006, Opt. Exp., 14, 7515 [NASA ADS] [CrossRef] [Google Scholar]
 Galicher, R., Marois, C., Macintosh, B., Barman, T., & Konopacky, Q. 2011, ApJ, 739, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Gelman, A., & Rubin, D. B. 1992, Statist. Sci., 7, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. 2014, Bayesian data analysis, Vol. 2 (Chapman & Hall/CRC) [Google Scholar]
 Givens, G. H., & Hoeting, J. A. 2012, Comput. statist., Vol. 710 (John Wiley & Sons) [Google Scholar]
 Gomez Gonzalez, C., Wertz, O., Christiaens, V., Absil, O., & Mawet, D. 2016a, Astrophysics Source Code Library [record ascl: 1603.003] [Google Scholar]
 Gomez Gonzalez, C. A., Absil, O., Absil, P.A., et al. 2016b, A&A, 589, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [CrossRef] [MathSciNet] [Google Scholar]
 Gotberg, Y., Davies, M. B., Mustill, A., Johansen, A., & Church, R. P. 2016, A&A, 592, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goździewski, K., & Migaszewski, C. 2009, MNRAS, 397, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Goździewski, K., & Migaszewski, C. 2014, MNRAS, 440, 3140 [NASA ADS] [CrossRef] [Google Scholar]
 Guerri, G., Daban, J.B., RobbeDubois, S., et al. 2011, Exp. Astron., 30, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Guyon, O., Bendek, E. A., Eisner, J. A., et al. 2012, ApJS, 200, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Halko, N., Martinsson, P.G., & Tropp, J. A. 2011, SIAM Rev., 53, 217 [CrossRef] [Google Scholar]
 Hibon, P., Vigan, A., Dohlen, K., et al. 2016, in Groundbased and Airborne Instrumentation for Astronomy VI, Proc. SPIE, 9908, 99080 [Google Scholar]
 Hinz, P. M., Rodigas, T. J., Kenworthy, M. A., et al. 2010, ApJ, 716, 417 [NASA ADS] [CrossRef] [Google Scholar]
 Hogg, D. W., Bovy, J., & Lang, D. 2010, ArXiv eprints [arXiv:1008.4686] [Google Scholar]
 Hughes, A. M., Wilner, D. J., Andrews, S. M., et al. 2011, ApJ, 740, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Hyndman, R. J. 1996, Am. Statist., 50, 120 [Google Scholar]
 Kalas, P., Graham, J. R., Fitzgerald, M. P., & Clampin, M. 2013, ApJ, 775, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Konopacky, Q. M., Marois, C., Macintosh, B. A., et al. 2016, AJ, 152, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, ApJ, 660, 770 [NASA ADS] [CrossRef] [Google Scholar]
 Lafrenière, D., Marois, C., Doyon, R., & Barman, T. 2009, ApJ, 694, L148 [NASA ADS] [CrossRef] [Google Scholar]
 Lagrange, A.M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Langlois, M., Vigan, A., Moutou, C., et al. 2012, in Adaptive Optics Systems III, Proc. SPIE, 8447, 84473 [CrossRef] [Google Scholar]
 MacKay, D. J. 2003, Information theory, inference and learning algorithms (Cambridge University Press) [Google Scholar]
 Maire, A.L., Skemer, A. J., Hinz, P. M., et al. 2015, A&A, 576, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maire, A.L., Langlois, M., Dohlen, K., et al. 2016, in Groundbased and Airborne Instrumentation for Astronomy VI, Proc. SPIE, 9908, 990834 [CrossRef] [Google Scholar]
 Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556 [NASA ADS] [CrossRef] [Google Scholar]
 Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Marois, C., Macintosh, B., & Véran, J.P. 2010a, in Adaptive Optics Systems II, Proc. SPIE, 7736, 77361 [Google Scholar]
 Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 2010b, Nature, 468, 1080 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Matthews, B., Kennedy, G., Sibthorpe, B., et al. 2014, ApJ, 780, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Mawet, D., David, T., Bottom, M., et al. 2015, ApJ, 811, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Meeus, J. 1998, Astronomical algorithms (Richmond: WillmannBell) [Google Scholar]
 Meshkat, T., Kenworthy, M. A., Quanz, S. P., & Amara, A. 2014, ApJ, 780, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Metchev, S., Marois, C., & Zuckerman, B. 2009, ApJ, 705, L204 [NASA ADS] [CrossRef] [Google Scholar]
 Meyn, S., & Tweedie, R. L. 2009, Markov Chains and Stochastic Stability, 2nd edn. (New York: Cambridge University Press) [Google Scholar]
 Morzinski, K. M., Males, J. R., Skemer, A. J., et al. 2015, ApJ, 815, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Nelder, J. A., & Mead, R. 1965, Comput. J., 7, 308 [CrossRef] [Google Scholar]
 Pathak, P., Guyon, O., Jovanovic, N., et al. 2016, PASP, 128, 124404 [NASA ADS] [CrossRef] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes: The Art of Scientific Computing, 3rd edn. (New York: Cambridge University Press) [Google Scholar]
 Pueyo, L., Soummer, R., Hoffmann, J., et al. 2015, ApJ, 803, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Reidemeister, M., Krivov, A. V., Schmidt, T. O. B., et al. 2009, A&A, 503, 247 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Service, M., Lu, J. R., Campbell, R., et al. 2016, PASP, 128, 095004 [NASA ADS] [CrossRef] [Google Scholar]
 Soummer, R. 2005, ApJ, 618, L161 [NASA ADS] [CrossRef] [Google Scholar]
 Soummer, R., Brendan Hagan, J., Pueyo, L., et al. 2011, ApJ, 741, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Student, W. S. G. 1908, Biometrika, 6, 1 [CrossRef] [Google Scholar]
 Su, K. Y. L., Rieke, G. H., Stapelfeldt, K. R., et al. 2009, ApJ, 705, 314 [NASA ADS] [CrossRef] [Google Scholar]
 van Leeuwen, F. 2007, A&A, 474, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yelda, S., Lu, J. R., Ghez, A. M., et al. 2010, ApJ, 725, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Zurlo, A., Vigan, A., Galicher, R., et al. 2016, A&A, 587, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Archival astrometric data
Table A.1 compiles all astrometric measurements available in the literature for the HR8799 system.
Compilation of astrometric measurements for HR8799bcde available in the literature and used in our orbital analysis.
Appendix B: Illustration of the orbital solutions
Figures B.1 and B.2 give an illustration of the highly probable orbits resulting from our MCMC simulations.
Fig. B.1 Illustration of the onsky projection of 1000 allowable orbits for HR8799bcde, all characterized by and by Keplerian parameters in the confidence intervals reported in Table 5. The yellow orbit corresponds to the solution reported in Table 5. 

Open with DEXTER 
Fig. B.2 Same as Fig. B.1, zooming on the part of the orbit where astrometric measurements are available. 

Open with DEXTER 
All Tables
Final HR8799bcde astrometric measurements with respect to the host star for nights of December 4–6, and 8, 2014.
The final HR8799bcde astrometric measurements with respect to the star for epoch 2014.93.
Confidence intervals (1σ) and best solutions in terms of for all the Keplerian elements, as well as for the orbital period P.
Relative probabilities (in %) associated with different meanmotion resonances between consecutive pairs of planet.
Compilation of astrometric measurements for HR8799bcde available in the literature and used in our orbital analysis.
All Figures
Fig. 1 Histogram of the horizontal and vertical offsets of the star with respect to the center of the frame in the December 5 (right) data cube. The vertical line represents the median of the histogram, and was used to globally recenter the data cube. The horizontal axis is in pixels, one pixel corresponding to 12.25 mas on sky. 

Open with DEXTER  
In the text 
Fig. 2 Illustration of a fullframe PCA ADI postprocessed SPHERE/IRDIS image of HR8799 acquired with broadband H filter (left part) during the night of December 4, 2014. The central part was masked with a disk of radius 20 pixels. 

Open with DEXTER  
In the text 
Fig. 3 Illustration of the annuluswise PCA postprocessing and merit function evaluation used in the negative fake companion technique. Left: no NEGFC was injected before annuluswise PCA processing. Right: a NEGFC was injected at the position and flux minimizing the merit function. The white circle illustrates the fixed circular aperture from which the pixel values I_{j} have been extracted to evaluate the merit function. The same colorscale was used for both images. 

Open with DEXTER  
In the text 
Fig. 4 Illustration of a typical corner plot obtained from the MCMC simulations using the NEGFC technique. The target companion is HR8799b observed during the night of December 6, 2014. The radial distance r (in pixels) and azimuth θ (in degrees) are detector coordinates with respect to the host star. The diagonal panels illustrate the posterior PDFs while those offaxis illustrate the correlation between them. 

Open with DEXTER  
In the text 
Fig. 5 Astrometry for HR8799bcde observed during the nights of December 4–6, and 8, 2014. The positions obtained from the left (resp. right) data cubes are represented with downward (resp. upward) black triangles. The error bars on the individual data points take into account all the contributions discussed in Sects. 3.1–3.4. The red dots correspond to the final astrometric measurements for each planet, together with the final error bar discussed in Sect. 3.5. The dashed lines represent the best orbital solutions for each planet in terms of reduced χ^{2} reported in Table 5. 

Open with DEXTER  
In the text 
Fig. 6 Speckle noise estimation for HR8799b observed on December 6, 2014. The histograms illustrate the offsets between the true position/flux of a fake companion and its position/flux obtained from the NEGFC technique. The dashed lines correspond to the 1D Gaussian fit from which we determine the speckle noise. 

Open with DEXTER  
In the text 
Fig. 7 Results of the MCMC simulations for HR8799b, displayed as a corner plot for the Keplerian elements a, e, i, Ω and ω. The diagonal panels illustrate the posterior PDFs while the offaxis panels illustrate the correlation between the parameters. The yellow lines and crosses correspond to the best solution in terms of reduced χ^{2}. 

Open with DEXTER  
In the text 
Fig. 8 Same as Fig. 7, but for HR8799c. 

Open with DEXTER  
In the text 
Fig. 9 Same as Fig. 7, but for HR8799d. 

Open with DEXTER  
In the text 
Fig. 10 Same as Fig. 7, but for HR8799e. 

Open with DEXTER  
In the text 
Fig. 11 Coplanarity test for HR8799bcde. The scatter plot illustrates the vector coordinates on the sky plane (see text) for the HR8799b (blue), HR8799c (red), HR8799d (orange) and HR8799e (green) allowable orbits. All the points on a given circle refer to orbits characterized by the same inclination. Similarly, all points on a given radial branch refer to orbits characterized by the same longitude of the ascending node. The angle between two spokes corresponds to 20°. The solid lines represent the 68% confidence interval around the most probable vector coordinates, which are depicted by thick points. 

Open with DEXTER  
In the text 
Fig. 12 Histograms of the period ratios P_{b}/P_{c} (left panel), P_{c}/P_{d} (middle panel) and P_{d}/P_{e} (right panel) derived from the semimajor axis posterior PDFs given in Figs. 7 to 10 and adopting 1.51 M_{⊙} for the host star mass. These period ratios allow to identify which meanmotion resonances are compatible with the allowable orbits derived from our MCMC analysis. It appears clear that the most compatible meanmotion resonance is 1e:2d:4c:8b. 

Open with DEXTER  
In the text 
Fig. B.1 Illustration of the onsky projection of 1000 allowable orbits for HR8799bcde, all characterized by and by Keplerian parameters in the confidence intervals reported in Table 5. The yellow orbit corresponds to the solution reported in Table 5. 

Open with DEXTER  
In the text 
Fig. B.2 Same as Fig. B.1, zooming on the part of the orbit where astrometric measurements are available. 

Open with DEXTER  
In the text 