An exact, generalised Laplace resonance in the HR 8799 planetary system
DD RAFT VERSION S EPTEMBER
16, 2020Typeset using L A TEX twocolumn style in AASTeX63
An exact, generalised Laplace resonance in the HR 8799 planetary system K RZYSZTOF G O ´ZDZIEWSKI AND C EZARY M IGASZEWSKI Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus Univ., Grudzia¸dzka 5, Toru´n, Poland
ABSTRACTA system of four super-Jupiter planets around HR 8799 is the first multi-planet configuration discovered viathe direct imaging technique. Despite over decade of research, the system’s architecture remains not fully re-solved. The main difficulty comes from still narrow observing window of ∼ years that covers small arcs oforbits with periods from roughly to years. Soon after the discovery it became clear that unconstrainedbest-fitting astrometric configurations self-disrupt rapidly, due to strong mutual gravitational interactions be-tween (cid:39) -Jupiter-mass companions. Recently, we showed that the HR 8799 system may be long term stablewhen locked in a generalized Laplace 8:4:2:1 mean motion resonance (MMR) chain, and we constrained itsorbits through the planetary migration. Here we qualitatively improve this approach by considering the MMRin terms of an exactly periodic configuration. This assumption enables us to construct for the first time theself-consistent N -body model of the long-term stable orbital architecture, using only available astrometric po-sitions of the planets relative to the star. We independently determine planetary masses, which are consistentwith thermodynamic evolution, and the parallax overlapping to σ with the most recent GAIA DR2 value. Wealso determine the global structure of the inner and outer debris discs in the [8, 600] au range, consistent withthe updated orbital solution. Keywords: planets and satellites: dynamical evolution and stability — stars:individual (HR 8799) — astrometry— methods: numerical — celestial mechanics INTRODUCTIONSeveral approaches are being used to detect extrasolarplanets. Indirect methods, such as the radial velocity (Mayor& Queloz 1995), transits (Henry et al. 2000), timing (Wol-szczan & Frail 1992), classic astrometry (Muterspaugh et al.2010) rely on studying radiation of the central star, whilethe planets themselves are not observed. The imaging tech-nique detects the planets directly, given their own infra-red(IR) radiation. This method, limited by the contrast, stabil-ity and resolution of the images is sensitive for massive andyoung planets in wide orbits. Therefore, even for a nearbystar HR 8799 located ∼ pc from the Sun (Brown et al.2018) it is only possible to detect the planets with long peri-ods of – years, as discovered by Marois et al. (2008,2010). This makes the orbits determination a very difficulttask. The measured astrometric positions of the planets rela-tive to the star are typically uncertain to a few . (cid:48)(cid:48) (mas),but this is still not sufficient to uniquely constrain the orbits.Qualitatively different architectures are consistent with thepresent observations (e.g. Wertz et al. 2017). Astrometric(purely geometric) orbital models are strongly unstable (e.g.,Konopacky et al. 2016; Wang et al. 2018), however, thereare also reported dynamically tuned configurations which,though hardly chaotic in the Lyapunov exponent sense, can survive for hundreds of Myrs (Götberg et al. 2016), compa-rable with the age of HR 8799 of (cid:39) – Myr (Maroiset al. 2010; Wilner et al. 2018).On the other hand, a resonant or near resonant system re-sulting from the convergent migration was shown to explainthe observations as well (Go´zdziewski & Migaszewski 2014,2018, furthermore, GM14 and GM18, respectively). We jus-tified a rigorously stable 8:4:2:1 MMR as the most likely ar-chitecture on the dynamical grounds, consistently with therecent studies in (Konopacky et al. 2016; Wang et al. 2018).They independently found, imposing MCMC dynamical pri-ors, that coplanar orbits near the 8:4:2:1 MMR result in or-ders of magnitude more stable orbits than any other scenario,and provide adequate fits to the measurements. In the presentwork, we extend the MMR hypothesis by linking the putativeresonance chain with the planetary N -body periodic solu-tions (Hadjidemetriou 1976; Hadjidemetriou & Michalodim-itrakis 1981, see also citations therein). In the later paper,they computed families of periodic orbits (POs) for the 4-body planetary system, and applied the results to the Galileanmoons of Jupiter. This follows de Sitter’s mathematical the-ory of the Laplace resonance in a Newtonian framework. DeSitter found a family of stable POs as the PoincarÃl’ orbitsof the second kind (Broer & Hançmann 2016). These au-thors also propose that librations (quasi-periodic solutions) a r X i v : . [ a s t r o - ph . E P ] S e p G O ´ZDZIEWSKI & M
IGASZEWSKI near these POs may provide realistic explanation of the ob-servations. We apply a similar reasoning to the HR 8799system.Because the 8:4:2:1 MMR chain in the HR 8799 systemgeneralizes the Laplace resonance (Papaloizou 2015), thereis a fairly obvious link of this prior resonant model found inGM14 and GM18, with a PO interpreted as the MMR cen-ter. Here, similarly to (Hadjidemetriou & Michalodimitrakis1981), we consider a planetary system as a PO when the os-culating orbital elements of the orbits are periodic in timew.r.t. non-uniformly rotating reference frame tied to the os-culating apsidal line of a selected planet. The present study isinspired by our finding that the 8:4:2:1 MMR configurationsfitting the observations of HR 8799 are in fact very close toan exact PO of the 5-body system.Periodic configurations are known to result from smoothconvergent migration in systems of two and more planets(e.g. Beaugé et al. 2006; Migaszewski 2015). In our nu-merical simulations of convergent migration (see GM14 andGM18), three– and four–planet MMR chains of the Laplaceresonance appear naturally, in wide ranges of the migrationtime scales and planet masses. The migration quickly drivesplanets to long-term stable systems, typically in a few Myrtime scale. This indicates that the 8:4:2:1 MMR capture maybe an efficient process, weakly dependent on physical con-ditions in the protoplanetary disk. Simultaneously, stablesystems are confined to tiny, isolated islands in the orbitalparameter space, as narrow as (cid:39) . au in semi-major axesand (cid:39) . in eccentricity (GM14 and GM18), and that re-flects deterministic character of the migration. Also the in-creased eccentricities of the innermost planets observed inthe best-fitting, stable solutions and in our simulations is con-sistent with an early evolutionary period of convergent in-ward migration of all four planets, trapping them pairwise in2:1 MMRs and pumping of the orbital eccentricities while inresonance lock (Wang et al. 2018; Yu & Tremaine 2001).Our new method improving the migration constrained op-timisation (MCOA) in GM14 and GM18 relies on two at-tributes of periodic configurations. When modelling the datawith stable PO, the long-term dynamical stability is guaran-teed per-se. Also, similarly to MCOA, instead of exploringlarge, n -dimensional parameter space, where the number offree parameters n = 23 for a co-planar system of four plan-ets, including their masses and the parallax, we may limit theoptimisation to a p -dimensional manifold embedded in thisspace; here, as explained below, p = 11 , or – if the massesand the parallax are fixed (given a’priori) – n = 18 vs. p = 6 .Therefore the MMR (periodic) constraint makes it possibleto substantially reduce the number of free parameters char-acterising the orbital configuration, and to avoid degeneracycaused by a small ratio of the data points to the degrees offreedom. But the advantage of the PO–constrained method over thestandard orbit fitting lies not only in the reduction of the pa-rameter space. In modeling a generic planetary configura-tion, the orbital elements as well as the planetary masses mustbe all free parameters, independent one of another. There-fore, the masses cannot be determined from the astromet-ric observations, unless they sufficiently map the orbits orare sufficiently precise to make it possible to detect mutualgravitational perturbations. In turn, the orbital elements ofa periodic (exactly resonant) solution strongly depend on theplanet masses and the total angular momentum of the system,assuming the linear scale of the system and the central star’smass to be given (Hadjidemetriou 1976; Hadjidemetriou &Michalodimitrakis 1981). Therefore these critical parame-ters may be derived with relatively short orbital arcs, and in-dependently of the planets’ cooling theory. As we also showbelow, because the PO impose tight timing on the orbital evo-lution, it is already possible to indirectly measure the systemparallax . The planet masses and parallax determined fromthe relative astrometry which are self-consistent with the as-trophysically fixed stellar mass, establish a test-bed and abenchmark for our hypothesis.We describe the results of the PO model of the HR 8799system in Sect. 2, the global structure of debris discs inSect. 3, and main conclusions in Sect. 4. The details andsupplementary material are given in Appendix. FITTING THE EXACT LAPLACE RESONANCEIn order to test the PO hypothesis, we used the mostearly HST observations in (Lafrenière et al. 2009; Soum-mer et al. 2011), a homogeneous, uniformly reduced datain (Konopacky et al. 2016), as well as the most recent re-fined Gemini Planet Imager (GPI) observations in (De Rosaet al. 2020), and the most accurate detection of HR 8799ein (Lacour et al. 2019) with the GRAVITY instrument. Thisprimary set does not contain all observations available in theliterature, and we limited the data in order to reduce pos-sible observational biases due to different instruments andpipelines, but to extend the observational window as much aspossible. This approach follows our earlier work GM18 andalso Wang et al. (2018). The data set D consists of N obs = 65 astrometric planet positions (the right ascension RA i ≡ α i ,and declination DEC i ≡ δ i , i = 1 , . . . , N obs ) relative to thestar with the mean uncertainty (cid:39) mas. However, there isa particular datum from the GRAVITY with ( α, δ ) errors assmall as 0.07 mas and 0.2 mas, respectively. This precisiondetection with the optical interferometry seems to be criticalfor constraining the best-fitting solutions. As the epoch t of the osculating initial condition we chose the date of thefirst HST observation t = 1998 . in (Lafrenière et al.2009; Soummer et al. 2011). We also tested the PO modelagainst all N obs = 127 measurements available in the litera- N EXACT L APLACE RESONANCE IN THE
HR 8799
PLANETARY SYSTEM Table 1.
The best-fitting, strictly periodic model of the HR 8799 planetary system. Median values and the σ uncertainties resulting fromthe DE-MC sampling are given in top row of each parameter, while values in the bottom row with more significant digits, in order to closelyreproduce the PO solution, correspond to the best-fitting parameters in terms of the χ statistics. Gaussian planet mass priors from the hot-startevolutionary models are (5 . ± . m Jup for HR 8799b, and (7 . ± . m Jup for all other planets (Wang et al. 2018), and the parallax prior is (24 . ± . mas (Brown et al. 2018). The first part of the Table is for the primary fit parameters x , and the last four rows are for the inferredosculating, astrocentric Keplerian elements at the epoch of t = 1998 . , the first measurement in (Lafrenière et al. 2009; Soummer et al.2011). The mass of the parent star (1 . ± . M (cid:12) (Konopacky et al. 2016) is fixed at its nominal value, in order to avoid the mass–orbitalscale correlation. For the model with p ≡ dim x = 11 free parameters, χ = 142 . and an RMS (cid:39) . mas.parameter/planet HR 8799e HR 8799d HR 8799c HR 8799bPlanet mass, m ( m Jup ) 7 . ± . . ± . . ± . . ± . . . . . Longitude of ascending node, Ω (deg) . ± . . Inclination, I (deg) . ± . . Parallax, Π (mas) . ± . . Period ratio, P d /P e . ± . . Scale factor, ρ . ± . . Relative Phase, t phase (yr) . ± . . Rotation angle, ω rot (deg) ± . Semi-major axis, a ( au ) 16 . ± .
04 26 . ± .
08 41 . ± .
11 71 . ± . . . . . Eccentricity, e . ± . . ± . . ± . . ± . . . . . Argument of pericentre, ω (deg) . ± . ± . ± . . ± . . . . . Mean anomaly, M (deg) − . ± . . ± . . ± . − . ± . − . . . − . ture, as listed in (Wertz et al. 2017), and updated with neweror re-reduced KECK, GPI and GRAVITY points (Appendix,Sect. AI.2). In both cases, the fitting results and conclusionsclosely overlap.Fitting a PO to the astrometric data is similar to ourapproach in GM14 and GM18, in which a migration-constrained co-planar solution is appropriately transformedto be consistent with the observations. Here, the optimisationprocess is essentially deterministic (fully reproducible), bet-ter constrained, regarding the masses and parallax as free pa-rameters of the dynamical model, and CPU-efficient — com-putations may be performed on a single workstation. Insteadof simulating the migration, for given masses m e , m d , m c , m b and C or, equivalently, the period ratio κ = P d /P e ofthe inner pair of planets, we find a strictly periodic, co-planar resonant solution (see Appendix for details), howeverwith some arbitrary relative phases of the planets. Then, using the N -body dynamics scale invariance, the inferred“raw” semi-major axes are linearly re-scaled with the fac-tor ρ , the orbital plane rotated to the sky plane by threeEuler angles ( I, Ω , ω rot ) , and planets propagated along thePO with the N -body integrator to epoch t phase equivalent tothe first observation epoch t . For given observation epochs t i , i = 1 , . . . , N obs , Cartesian coordinates of the planetsare re-scaled by the parallax Π in order to obtain the an-gular positions ( α i , δ i ) . Then the ephemeris may be com-pared with the observations. To quantify that, we constructthe merit function, e.g., the common χ or other goodness-of-fit measure, such as the maximum likelihood function L or the Bayesian posterior P . In the most general set-tings, the merit function depends on free parameters, x ≡ ( m e , m d , m c , m b , κ, ρ, I, Ω , ω, t phase , Π ), i.e., masses ofthe planets, period ratio of two inner planets, linear scale fac-tor, Euler angles, the epoch and the parallax. At this point G O ´ZDZIEWSKI & M
IGASZEWSKI
Figure 1.
Astrometric observations (red, green, blue and magenta points for planets HR 8799e,d,c,b, respectively) in the set used for theanalysis. The grey thick curves illustrate the best-fitting orbits (Tab. 1). The light grey curves in the left-hand panel mark referencing circularorbits of radii , , , . . . , au in the orbital plane of the system. The red, blue, green and magenta arrows point to the periastron ofeach orbit. The close-up of the GRAVITY datum is illustrated in a small panel in the top-right corner of the left plot. Grey curves represent randomly chosen orbits from the DE-MC sampling, the black dots mark positions at the orbits in the epoch of the GRAVITY observation( . ). The graph is centred at the datum, the axes are expressed in mili-arcseconds (mas). Black curves are for σ , σ and σ confidenceintervals for the model position at the right ascension–declination, (RA, DEC)–plane at the epoch of . , derived from the DE-MCsampling. The right-hand panel illustrates the observations and the model orbits as the time-functions of RA and DEC. Figure 2.
Astrocentric positions of the HR 8799 planets over Gyr N -body integration presented in the orbital plane co-rotating withHR 8799b. The red, green and blue curves illustrate the orbits ofHR 8799e, d c, respectively. Big red, green, blue and magenta sym-bols mark the positions of planets e, d, c and b during the conjunc-tion of planets b and c, while smaller symbols denote the positionsduring their opposition. The yellow symbol at the origin marks theparent star. the data modelling becomes almost the standard optimization— almost, and not really trivial, since we now must seek thebest-fitting x constrained to a manifold of a stable PO family,representing the particular MMR chain. One needs to find anextremum of the merit function, as well as to estimate uncer-tainties of the best-fitting parameters (e.g., Gregory 2010).Details and variants of our experiments, regarding optimi-sation on the PO manifold, are described in Appendix. Here,we quote the final best-fitting parameters in Table 1. Thefirst part of this Table shows the primary fit parameters x ,and the bottom part is for the derived, osculating astrocentricKeplerian elements at the epoch t = randomly chosen synthetic orbitsare shown with grey curves, while black points denote thepositions of the synthetic solutions in the epoch of the obser-vation. Grey oval contours mark σ , σ and σ confidenceintervals stemming from the DE-MC sampling. The right-hand panel illustrates the RA and DEC of the model and ob- N EXACT L APLACE RESONANCE IN THE
HR 8799
PLANETARY SYSTEM Figure 3.
Posterior probability distributions of the planets’ masses, inclination, the longitude of the ascending node, and the parallax. Theshades of grey indicate the σ, σ and σ confidence ranges of the parameters. Red symbols with circle/ellipse contours in the top panels showthe astrophysical mass constraints (Wang et al. 2018), which are the priors in the DE-MC sampling. The red curve in the bottom-right panel isthe prior put on the parallax, according to the GAIA DR2 catalogue (Brown et al. 2018). See also Table 1 and its caption. servations as functions of the epoch. The PO described inTable 1 yields the reduced χ ν (cid:39) . for p = 11 free param-eters, N obs = 65 , ν = 119 and the RMS (cid:39) . mas comparesto the mean uncertainty of the measurements (cid:39) mas. It ad-equately explains the data in a statistic sense. In particular,the time– and sky-plane– synchronisation of the model withthe GRAVITY datum (left panel in Fig. 1) is apparently per-fect.The orbital evolution of this best-fitting system, integratedfor Gyr, is presented in Fig. 2. This figure shows orbitsof HR 8799e, HR 8799d and HR 8799c in a reference frameco-rotating with HR 8799b. All trajectories are closed, con-sistent with the periodic evolution of the system. The po-sitions of the planets are shown only in epochs of conjunc-tions between HR 8799b and HR 8799c (big filled circles)as well as their oppositions (small circles). Both the con-junctions and the oppositions repeat in the same pattern. Thesystem is then an exact 8:4:2:1 MMR chain, consisting oftriple two-body 2:1 MMRs of subsequent pairs of planets,with librations of the critical angle of the zero-th order 4-body MMR φ = λ e − λ d − λ c + 2 λ b (where λ b,c,d,e are the mean longitudes of the planets) with a small ampli-tude (cid:39) degrees (Fig. A2). It is worth noting that while themean orbital osculating period ratios are (cid:39) . , (cid:39) . and (cid:39) . for the innermost to outermost pairs of planets, re- spectively, the canonical (proper) mean motion frequencies(Morbidelli 2002) ratios are equal to / , indicating exact 2-body 2:1 MMRs. Therefore the MMR chain is understoodas the generalized Laplace resonance. In order to illustratethe long term stability of the model, we computed dynami-cal maps in terms of the Mean Exponential Growth factor ofNearby Orbits (MEGNO aka (cid:104) Y (cid:105) , Cincotta et al. 2003) foreach planet. The integration interval of 10 Myrs translatesto (cid:39) , outermost orbits, sufficient to detect short-term,MMR-induced instability. Remarkably, the maps (Fig. A3)are similar to our earlier Fig. 9 in GM14, illustrating theMCOA model build upon much narrower data window, andstill consistent with the updated periodic model of the sys-tem.Uncertainties of the parameters are illustrated in Fig. 3(also Figs. A8 and A9 in Appendix). Two top panels are forthe mass–mass diagrams. Red points with shaded ellipsesindicate Gaussian priors imposed on the masses consistentwith the hot-start cooling theory (Wang et al. 2018), whilegrey filled contours denote σ , σ and σ confidence inter-vals of the posterior probability distributions. Apart of theHR 8799d mass, the posterior closely fits with the astrophys-ical constraints. The bottom-left panel shows the posteriordistribution of the orbital inclination and the longitude of as-cending node. These parameters exhibit substantial correla- G O ´ZDZIEWSKI & M
IGASZEWSKI
Figure 4.
The inner part of the outer debris disk revealed by (cid:39) . × (cid:104) Y (cid:105) -stable orbits found in the whole debris discs sim-ulation, illustrated as a snapshot of astrocentric coordinates ( x, y ) and osculating orbital eccentricities e of these orbits at the initialepoch, color-coded and labeled in the top bar. Initial positions ofplanets are marked with filled circles. Gray rings illustrate their or-bits integrated in a separate run for 10 Myr. tions, yet much reduced thanks to the priors. The bottom-right panel is for the parallax, nominally agreeing to (cid:39) . with the GAIA DR2 value. The Gaussian prior as the GAIAparallax (Brown et al. 2018, the red curve) closely overlapswith the DE-MC posterior. RESONANT STRUCTURE OF DEBRIS DISCSThe orbits of the planets likely share the common planewith the outer debris disk (Matthews et al. 2014; Booth et al.2016; Read et al. 2018; Wilner et al. 2018). Determinationof the debris disk structure with the infra-red and millime-tre observations is still not fully conclusive, both in termsof the orientation as well as the inner edge r inner of the disk(Booth et al. 2016). They argue that the structure of thedisk might be a footprint of a fifth, yet unseen planet be-yond HR 8799b. Read et al. (2018) proposed such an addi-tional planet HR 8799f with the mass and semi-major axisof . m Jup and 138 au that could predict the outer belt’sedge and explain the ALMA observations. Later, Wilneret al. (2018), with observations at the Submillimeter Ar-ray at 1340 µ , detected the inner edge of the debris diskat , − au, and the disk extending to (cid:39) au.They also constrained the mass of outer planet HR 8799bto (cid:39) m Jup . Remarkably, it is close to our best-fitting value.Furthermore, Geiler et al. (2019) found that a single, wide planetesimal disk does not reproduce the observed emissionsand proposed a two-population model, comprising a Kuiper-Belt-like structure of low-eccentricity planetesimals and ascattered disk comprising of high-eccentricity population ofcomets.With the new, strictly resonant configuration of the fourplanets, including their updated masses and the parallax,we conducted preliminary N -body simulations resulting in . × small-mass asteroids, that reveal the global dynam-ical structure of the debris discs (see Appendix, Sect. AIII fordetails). The inner border of the outer disk (Figs. 4 and A15)is significantly non-symmetric, with non-uniform density ofasteroids, which may bias the disk orientation angles derivedfrom simple models assuming the axial symmetry. The in-ner edge from our simulations agrees with the observationalmodel of Wilner et al. (2018). Moreover, we found a ringof high-eccentricity asteroids at (cid:39) − au (Fig. A15),close to the inner edge reported in (Booth et al. 2016; Readet al. 2018), which results in locally increased velocity dis-persion. The velocity dispersion could impose higher dustproduction rate and stronger emission, making the disk ra-dial intensity profile no longer consistent with a simple powerlaw. DISCUSSION AND FUTURE RESEARCHUnder the PO hypothesis, which is justified on the dynami-cal and the system formation grounds, the present astrometricdata of the HR 8799 planets make it possible to determine notonly the parallax but also their masses, independently fromthe cooling theory. In order to illustrate this prediction, wesimulated new synthetic observations around the best-fittingmodel in Table 1 with fixed m e = 7 m Jup , and Gaussiannoise equal to the GRAVITY datum uncertainty. We per-formed the χ minimisation without the planets’ masses pri-ors, adding new synthetic measurements after the last epochof each planet. The resulting time-series of the best-fittingm e and its σ range indicate (Fig. A1) that with merely onemore epoch (cid:39) . , all masses become meaningfully con-strained without prior information. If the HR 8799 systemis indeed represented by a PO, or a nearby stable resonantconfiguration, then it may be possible to determine the plan-ets’ masses basing solely on the relative astrometry. Thiscould be a test-bed for the cooling theory of HR 8799–like,massive planets, and possibly other multiple planetary sys-tems discovered via the direct imaging. The deterministic POmodel may serve as a reference configuration useful for theastrometric and physical characterisation of such resonant orclose to resonant systems.The PO hypothesis may be naturally confronted with com-pact multiple Kepler and super-Earth systems which are pre-dominately close to, but not actually inside of MMRs (e.g.Fabrycky et al. 2014). The planetary migration might eas- N EXACT L APLACE RESONANCE IN THE
HR 8799
PLANETARY SYSTEM
Beaugé, C., Michtchenko, T. A., & Ferraz-Mello, S. 2006,MNRAS, 365, 1160, doi: 10.1111/j.1365-2966.2005.09779.xBergfors, C., Brandner, W., Janson, M., Köhler, R., & Henning, T.2011, A&A, 528, A134, doi: 10.1051/0004-6361/201116493Booth, M., Jordán, A., Casassus, S., et al. 2016, MNRAS, 460,L10, doi: 10.1093/mnrasl/slw040Broer, H., & Hançmann, H. 2016, Indagationes Mathematicae,27, doi: 10.1016/j.indag.2016.09.002Brown et al. 2018, A&A, 616, A1,doi: 10.1051/0004-6361/201833051Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Physica DNonlinear Phenomena, 182, 151,doi: 10.1016/S0167-2789(03)00103-9Currie, T., Fukagawa, M., Thalmann, C., Matsumura, S., &Plavchan, P. 2012, ApJL, 755, L34,doi: 10.1088/2041-8205/755/2/L34Currie, T., Burrows, A., Itoh, Y., et al. 2011, ApJ, 729, 128,doi: 10.1088/0004-637X/729/2/128Currie, T., Burrows, A., Girard, J. H., et al. 2014, ApJ, 795, 133,doi: 10.1088/0004-637X/795/2/133De Rosa, R. J., Nguyen, M. M., Chilcote, J., et al. 2020, Journal ofAstronomical Telescopes, Instruments, and Systems, 6, 015006,doi: 10.1117/1.JATIS.6.1.015006Delisle, J. B. 2017, A&A, 605, A96,doi: 10.1051/0004-6361/201730857Esposito, S., Mesa, D., Skemer, A., et al. 2013, A&A, 549, A52,doi: 10.1051/0004-6361/201219212Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ,790, 146, doi: 10.1088/0004-637X/790/2/146Galicher, R., Marois, C., Macintosh, B., Barman, T., & Konopacky,Q. 2011, ApJL, 739, L41, doi: 10.1088/2041-8205/739/2/L41 Geiler, F., Krivov, A. V., Booth, M., & Löhne, T. 2019, MNRAS,483, 332, doi: 10.1093/mnras/sty3160Götberg, Y., Davies, M. B., Mustill, A. J., Johansen, A., & Church,R. P. 2016, A&A, 592, A147,doi: 10.1051/0004-6361/201526309Go´zdziewski, K., & Migaszewski, C. 2014, MNRAS, 440, 3140,doi: 10.1093/mnras/stu455—. 2018, ApJS, 238, 6, doi: 10.3847/1538-4365/aad3d3Gregory, P. 2010, Bayesian Logical Data Analysis for the PhysicalSciences (Cambridge University Press)Hadjidemetriou, J. D. 1976, Ap&SS, 40, 201,doi: 10.1007/BF00651199Hadjidemetriou, J. D., & Michalodimitrakis, M. 1981, A&A, 93,204Henry, G. W., Marcy, G. W., Butler, R. P., & Vogt, S. S. 2000,ApJL, 529, L41, doi: 10.1086/312458Hinz, P. M., Rodigas, T. J., Kenworthy, M. A., et al. 2010, ApJ,716, 417, doi: 10.1088/0004-637X/716/1/417Ida, S., Tanaka, H., Johansen, A., Kanagawa, K. D., & Tanigawa,T. 2018, ApJ, 864, 77, doi: 10.3847/1538-4357/aad69cIzzo, D., Ruci´nski, M., & Biscani, F. 2012, The Generalized IslandModel, Vol. 415 (Berlin, Heidelberg: Springer-Verlag), 151–169Johansen, A., & Lambrechts, M. 2017, Annual Review of Earthand Planetary Sciences, 45, 359,doi: 10.1146/annurev-earth-063016-020226Konopacky, Q. M., Marois, C., Macintosh, B. A., et al. 2016, AJ,152, 28, doi: 10.3847/0004-6256/152/2/28Lacour et al. 2019, A&A, 623, L11,doi: 10.1051/0004-6361/201935253Lafrenière, D., Marois, C., Doyon, R., & Barman, T. 2009, ApJL,694, L148, doi: 10.1088/0004-637X/694/2/L148 G O ´ZDZIEWSKI & M
IGASZEWSKI
Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122,doi: 10.1088/0004-637X/761/2/122Maire et al., A. L. 2015, A&A, 579,doi: 10.1051/0004-6361/201425185eMarois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322,1348, doi: 10.1126/science.1166585Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., &Barman, T. 2010, Nature, 468, 1080, doi: 10.1038/nature09684Matthews, B., Kennedy, G., Sibthorpe, B., et al. 2014, ApJ, 780,97, doi: 10.1088/0004-637X/780/1/97Mayor, M., & Queloz, D. 1995, Nature, 378, 355,doi: 10.1038/378355a0Metchev, S., Marois, C., & Zuckerman, B. 2009, ApJL, 705, L204,doi: 10.1088/0004-637X/705/2/L204Migaszewski, C. 2015, MNRAS, 453, 1632,doi: 10.1093/mnras/stv1739Morbidelli, A. 2002, Modern celestial mechanics : aspects of solarsystem dynamics (CRC Press, Advances in Astronomy andAstrophysics)Muterspaugh, M. W., Lane, B. F., Kulkarni, S. R., et al. 2010, AJ,140, 1657, doi: 10.1088/0004-6256/140/6/1657Papaloizou, J. C. B. 2015, International Journal of Astrobiology,14, 291, doi: 10.1017/S1473550414000147Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P.2002, Numerical recipes in C++ : the art of scientific computing(Cambridge University Press) Price, K., Storn, R. M., & Lampinen, J. A. 2005, DifferentialEvolution: A Practical Approach to Global Optimization(Natural Computing Series) (Berlin, Heidelberg:Springer-Verlag)Pueyo et al., L. 2015, ApJ, 803, doi: 10.1088/0004-637X/803/1/31Ramos, X. S., Charalambous, C., Benítez-Llambay, P., & Beaugé,C. 2017, A&A, 602, A101, doi: 10.1051/0004-6361/201629642Read, M. J., Wyatt, M. C., Marino, S., & Kennedy, G. M. 2018,MNRAS, 475, 4953, doi: 10.1093/mnras/sty141Šidlichovský, M., & Nesvorný, D. 1996, Celestial Mechanics andDynamical Astronomy, 65, 137, doi: 10.1007/BF00048443Soummer, R., Brendan Hagan, J., Pueyo, L., et al. 2011, ApJ, 741,55, doi: 10.1088/0004-637X/741/1/55Ter Braak, C. J. F. 2006, Statistics and Computing, 16, 239,doi: 10.1007/s11222-006-8769-1Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, AJ, 156, 192,doi: 10.3847/1538-3881/aae150Wertz, O., Absil, O., Gómez González, C. A., et al. 2017, A&A,598, A83, doi: 10.1051/0004-6361/201628730Wilner, D. J., MacGregor, M. A., Andrews, S. M., et al. 2018, ApJ,855, 56, doi: 10.3847/1538-4357/aaacd7Wolszczan, A., & Frail, D. A. 1992, Nature, 355, 145,doi: 10.1038/355145a0Yu, Q., & Tremaine, S. 2001, AJ, 121, 1736, doi: 10.1086/319401Zurlo et al., A. 2016, A&A, 587,doi: 10.1051/0004-6361/201526835
N EXACT L APLACE RESONANCE IN THE
HR 8799
PLANETARY SYSTEM Figure A1.
The best-fitting value of the innermost planet with σ uncertainty for the data set with additional synthetic measurements given insubsequent epochs up to the year of . See the main text (Sect. 4) for details. Figure A2.
Temporal evolution, for the first 10 Kyr, of the canonical, osculating orbital elements, expressed in the Jacobian reference frame,for the PO configuration in Tab. 1. Top-left panel: osculating period ratios for subsequent pair of planets, and their mean values (horizontal redlines) from top to bottom, P b /P c (cid:39) . , P c /P d (cid:39) . , and P d /P e (cid:39) . . We note that the proper orbital periods, in the sense of themean motions as fundamental frequencies (Morbidelli 2002), expressed in Julian years of 365.25 days, are . , . , . ,and . for planets HR 8799e,d,c,b, respectively, forming an exact 2:1, 2:1, 2:1 MMR chain. Top-right panel: eccentricities of theplanets HR 8799e,d,c,b from top to bottom, respectively. Bottom-left panel: one of the elements x i ≡ e i cos (cid:36) i , ( i = HR 8799b,c,d,e), usedto compute the secular frequency of the apsides rotation. The second component of the quasi-periodic signal (not shown) is y i ≡ e i sin (cid:36) i .Bottom-right panel: the critical argument of the zero-th order, four-body generalized Laplace resonance for the same initial condition thatlibrates around (cid:39) degrees. APPENDIX0 G
O ´ZDZIEWSKI & M
IGASZEWSKI
Figure A3.
Dynamical maps in the semi major-axis–eccentricity plane, in terms of the MEGNO fast indicator. Each point in the mapswas integrated for 10 Myrs, equivalent to more than 20,000 orbits of planet HR 8799b. Stable systems are confined to (cid:104) Y (cid:105) (cid:39) (blue) andyellow colour marks strongly unstable configurations typically self-disrupting in less than 1 Myr. The star symbol marks the nominal osculatingelements in the initial condition in Table 1, for each planet in subsequent panels. The resolution of each map is × pixels.AI. NUMERICAL SETUP AND ALGORITHMSAI.1.
Searching for periodic configurations
A co-planar orbital configuration of a planetary system is determined by a vector containing positions and velocities of the plan-ets { x i , y i , u i , v i } or, equivalently, by a vector consisting of astrocentric Keplerian elements of the orbits { a i , e i , (cid:36) i , M i } , i.e.,the semi-major axis, eccentricity, pericenter longitude, and the mean anomaly, respectively, where i = e, d, c, b or, equivalently, i = 1 , , , . Both the state vectors are given at a particular epoch and the state of the system at another epoch may be obtainedthrough propagating the initial condition with the numerical integration of the N -body equations of motion. A configuration ofthe planetary N -body problem is called periodic if after time interval T (called the period of PO) it returns to its initial state inthe reference frame rotating non-uniformly with the temporal (osculating) pericenter of a selected planet (Hadjidemetriou 1976).For the Keplerian representation of orbits we have the boundary conditions: a i ( T ) = a i (0) , e i ( T ) = e i (0) , M i ( T ) = M i (0) forall orbits. Since the angular momentum C of the system must be conserved, the pericenter longitudes (cid:36) i ( T ) (cid:54) = (cid:36) i (0) . However, ∆ (cid:36) i,j ≡ (cid:36) i − (cid:36) j ( i (cid:54) = j ) must be the same after T , ∆ (cid:36) i,j ( T ) = ∆ (cid:36) i,j (0) . This means that the relative configuration of theplanets remains fixed after the period, while the whole system rotates by a certain angle around its angular momentum vectoraligned with the z -axis (e.g. Lithwick et al. 2012).Here, we consider the evolution of the Keplerian elements periodic in a reference frame co-rotating with the apsidal line of theinnermost orbit — the reference orbit can be chosen freely. Equivalently, when considering the Cartesian coordinates, one needsto search for configurations whose positions and velocities expressed in a reference frame co-rotating with one of the planets ( ¯ x , ¯ y , ¯ u , ¯ v ) fulfil the periodicity conditions ¯ x i ( T ) = ¯ x i (0) , ¯ y i ( T ) = ¯ v x,i (0) , ¯ v y,i ( T ) = ¯ v y,i (0) (Hadjidemetriou 1976). We used N EXACT L APLACE RESONANCE IN THE
HR 8799
PLANETARY SYSTEM κ = ( P d /P e ) | M e =0 . Formally, for a chain ofthe 8:4:2:1 MMR, there are different epochs in which the innermost planet is in the pericentre. We select one of them and keepthis choice when continuing a given family w.r.t. other parameters of the solution, denoted as a generic parameters vector x .After testing various parameterisations of the PO in terms of numerical efficiency and reliability, we decided to use the follow-ing set of components of the state vector X , each of which being a function of the astrocentric, osculating Keplerian elements: X = log e e , X = log e d , X = log e c , X = log e b , X = P c /P d − ( P c /P d ) nom , X = P b /P c − ( P b /P c ) nom , X = (cid:36) e − (cid:36) d , X = (cid:36) d − (cid:36) c , X = (cid:36) c − (cid:36) b , X = M d , X = M c , X = M b , where ( P c /P d ) nom = (cid:104) C { − ( P d /P e ) } (cid:105) − , ( P b /P c ) nom = (cid:104) C { − ( P c /P d ) nom } (cid:105) − , C = j qi ( p + j ) , C = k pj ( r + k ) . The nominal values of the period ratios ( P c /P d ) nom and ( P b /P b ) nom corresponds to a chain of exact MMRs (Delisle 2017), P d /P e ≈ ( q + i ) /q , P c /P d ≈ ( p + j ) /p , P b /P c ≈ ( r + k ) /r . Therefore, for the case of the 2:1, 2:1, 2:1 MMRs chain, both thefactors C = C = 1 / . Although the relations given above were designed for weakly interacting systems, whose evolution iswell described with the averaging approach (Delisle 2017), we found that such a representation enables appropriately controllingthe period ratios.For a given (fixed) period ratio κ and masses m e , m d , m c , and m b , parameterising a given family of PO, its member is beingsearched for with the Newton method for nonlinear equations (Press et al. 2002), in 12-dimensional X space. Since ∆ X i ≡ X i ( T ) − X i (0) , where T is fixed so the innermost planet completes exactly full revolutions (counted from its pericentre topericentre), the set of non-linear equations to be solved reads as ∆ X i ( X , X , ..., X ) = 0 , where i = 1 , . . . , . At first, thestarting point is drawn randomly, yet around X , X (cid:39) (or close to an approximate solution, which we already know, such asthe 8:4:2:1 MMR fits found in GM14 and GM18) until the algorithm finds a solution with ∆ X i ≈ . Next, this solution can becontinued for changed κ and the planet masses. In general settings, the continuation of PO is a complex problem, since manystable and unstable families may exist, in different parameter ranges (Hadjidemetriou & Michalodimitrakis 1981, and referencestherein). AI.2. Data fitting on the parametric grid of PO
At the first attempt, we performed the optimisation similarly to the MCOA variant in GM18. Here, instead of CPU-demandingmigration simulations, which result in resonant but not necessarily periodic systems, we continued POs in a -dimensional Z grid ( Z = m e , Z = m d , Z = m c , Z = m b , Z = κ ), and we found ∼ configurations covering the interesting region of the8:4:2:1 MMR chain. In this sense, we obtained the exactly resonant (periodic) configurations which might fit the observationsas well, each being the 8:4:2:1 MMR center for different masses and the inner orbits period ratio. From this point, the analysisis essentially similar to GM18. In order to fit a given PO to the measurements, one needs to find a minimum of χ in a few-dimensional space of model parameters. We performed the optimisation experiments using the same parameters, as in GM18: thescale parameter ρ , the phase of a periodic configuration corresponding to the reference epoch t phase , and the 3-1-3 Euler angles ( I, Ω , ω rot ) , fixing orientation of the orbital plane w.r.t. the sky (observer) frame. In such settings, χ (equivalently, the maximumlikelihood L ) depends on Y , whose components are Y = ρ , Y = t phase , Y = Π , Y = I , Y = Ω , Y = ω rot . We also updatedthe parameter vector Y by the system parallax Π , and the so called error floor σ α,δ , re-scaling the nominal uncertainties, in orderto account for possibly underestimated errors and biases of the observations.The merit function L ( Y ) is defined for this variant of parametrisation as follows: ln L ( Y ) = − χ ( Y ) − N obs (cid:88) i =1 (cid:2) ln θ i,α + ln θ i,δ (cid:3) − N obs ln(2 π ) , (A1) χ ( Y ) = N obs (cid:88) i =1 (cid:34) [ α i − α ( t i , Y )] θ i,α + [ δ i − δ ( t i , Y )] θ i,δ (cid:35) , (A2)where ( α i , δ i ) are the measurements at time t i , α ( t i , Y ) , δ ( t i , Y ) are the ephemeris values, θ i,α and θ i,δ are the nominal mea-surements uncertainties in RA and DEC scaled in quadrature with the error floor, θ i,α = ( σ i,α + σ α,δ ) or θ i,δ = ( σ i,δ + σ α,δ ) ,2 G O ´ZDZIEWSKI & M
IGASZEWSKI for each datum, respectively, and N obs is the number of observations. Also, N = 2 N obs , since RA and DEC are measured in asingle detection. The ln L function in Eq. A1 is defined in the way that assuming the uncertainties as Gaussian and uncorrelated,the resulting best-fitting models should yield the reduced χ ν (cid:39) . The merit function was optimized with the help of genetic andevolutionary algorithms (Izzo et al. 2012).In order to illustrate the results of the PO-grid approach, we invoke a particular experiment in which we fitted the ln L functionin Eq. A1 to all measurements available at the moment in the literature: the early HST data in (Lafrenière et al. 2009) and(Soummer et al. 2011), a homogeneous data set in (Konopacky et al. 2016), GPI data in (Wang et al. 2018), and the GRAVITYmeasurement in (De Rosa et al. 2020) updated with mostly VLT/SPHERE, SUBARU and LBT data, collected in (Wertz et al.2017), from Metchev et al. (2009); Hinz et al. (2010); Currie et al. (2011, 2012, 2014); Bergfors et al. (2011); Pueyo et al. (2015);Galicher et al. (2011); Esposito et al. (2013); Maire et al. (2015); Zurlo et al. (2016). This set consists of N obs = 127 (RA, DEC)observations, some of them clearly deviating from any astrometric model. Because preliminary PO and MCOA fits indicated thereduced χ ν (cid:39) . for the best-fitting models, we introduced the error floor σ α,δ in order to account for possible data biases andun-modeled errors.The results are illustrated in Fig. A4, which shows selected fitted parameters, which are gathered on the pre-computed PO gridand plotted vs. ∆ ln L ≡ ln L max − ln L , relative to the best-fitting value ln L max found in the search. Most of the primary Y parameters, such as the mass of HR 8799d (top-left panel) and HR 8799b (not shown), the error floor σ α,δ (the bottom-rightpanel), and the system parallax (the bottom-left panel) exhibit clear extrema. Also the PO–constrained eccentricities (such as forHR 8799e in top-right panel) are quasi-parabolically bounded. We found particularly surprising that the masses of HR 8799dand HR 8799b could be potentially constrained, although the ln L extremum is apparently shallow. That also regards the parallax Π , which, as the free parameter of the astrometric model , overlaps with the GAIA DR2 trigonometric parallax Π (cid:39) (24 . ± . mas within its formal σ uncertainty – our best-fitting value Π (cid:39) . mas is accurate to the fourth significant digit( (cid:39) . , in other fits up to (cid:39) . ). We consider this as a meaningful benchmark of the self-consistency of the astrometricmodel, the parallax and the physical characteristics of the system: the derived masses of the planets and the adopted stellar massof 1.52M (cid:12) (Konopacky et al. 2016). We note that the stellar mass must be fixed or tigthly constrained, since otherwise it wouldintroduce a strong mass–period (linear scale) correlation through the Keplerian law (actually, the N -body dynamics scaling).The orbital geometry of the best-fitting models is illustrated in Fig. A5. Panels are for the close-ups of the sky-plane forsubsequent planets, with all data points marked with symbols for different sub-sets, gray curves for random solutions drawnwithin (cid:39) mas range around the best-fitting models which are illustrated with the red curves, and yielding ln L marginally worsethan the extremum value. The orbits are plotted for the time interval between t = 1998 . and the last GRAVITY epoch. The ∆ ln L –range in Fig. A4 translates to a substantial spread of the models, which is best visible for HR 8799d. Although we didnot compute formal confidence intervals, the spread indicates sensitivity of ∆ ln L to a variation of the parameters.The orbits plotted globally (Fig. A6), for the osculating orbital period of each planet separately, appear well bounded, and thisis particularly apparent for the innermost, fast moving planets. The model might predict their geometric positions close to thebest-fitting PO motion for a long time, in spite of the parameter uncertainties.The bad message received from the PO–grid fitting is a strong anti-correlation between the masses of HR 8799c,e (m c , m e );moreover, the best-fitting configurations exhibit very small innermost mass m e (cid:39) . m Jup . The grid method also depends onthe resolution, which should be individually tuned for each parameter in 5-dim space, as illustrated in the HR 8799d mass scan(top-left panel in Fig. A4). We found that the error floor does not help in eliminating nor even reducing the mass correlation.These preliminary PO–grid experiments provide interesting and useful hints for the final approach, described below (Sect.AI.3), yet we could not consider them as fully conclusive. Unfortunately, the parametric grid approach is tedious and introduceslarge CPU-overhead. The PO continuation and sampling must be multi-dimensional and that implies not only the need ofcomputing huge sets of the solutions but also optimising the initial conditions one by one — although we rather performed theMonte Carlo search on the grid. Determining the best-fitting model to the present observations is also difficult due to the (m c ,m e ) anti-correlation. Getting rid of that degeneracy needs additional, prior information, such as the planet masses estimated onthe grounds of the thermodynamical evolution and cooling tracks.In order to address these issues, we improved the grid MMR-constrained optimisation in ways making it CPU-efficient andindependent on the grid resolution, described below. In the final experiments, we examined the reduced data set ( D from hereafter)comprising of HST data in (Lafrenière et al. 2009; Soummer et al. 2011), a homogeneously reduced data set in (Konopacky et al.2016), and GPI data in (Wang et al. 2018) derived in a similar instrument, as well as the GRAVITY datum in (De Rosa et al.2020). We used the HST and GRAVITY data to extend the time–window of the observations. For this set, we did not considerthe error–floor, since the PO models with the basic p = 11 free parameters yield χ ν (cid:39) . and an RMS (cid:39) mas, consistentwith the mean uncertainty of the measurements in the reduced set, (cid:39) mas. N EXACT L APLACE RESONANCE IN THE
HR 8799
PLANETARY SYSTEM Figure A4.
Plots illustrating the PO–grid search for the best fitting parameters to astrometric measurements of the HR 8799 system. In thisexperiment, we fitted seven free parameters, including the orbit scale ρ , three Euler angles, the initial epoch, system parallax Π , and the errorfloor σ α,β , to all measurements gathered in the literature, N obs = 127 data points. We optimised the likelihood function ln L with the errorfloor term, and without any priors. The extremum of ln L (cid:39) − . , the error floor (cid:39) . mas. Fitted systems in two Monte-Carlo samplingruns on the PO grid, are marked with filled circles and different colors. Top-left panel: mass of HR 8799d with a clear extremum. Top-rightpanel: inferred eccentricity of HR 8799e from the PO models. Bottom-left panel: system parallax, with the GAIA DR2 nominal value (redline) marked with its σ (gold rectangle) and σ (grey rectangle) confidence intervals. Bottom-right plot: the error floor σ α,δ , as the globalcorrection factor of the measurements uncertainty, yielding the reduced χ ν (cid:39) . AI.3.
Optimisation on a manifold of periodic configurations
From the mathematical point of view, our goal is to find the best-fitting PO in the space of vectors Z = ( m e , m d , m c , m b , κ ) .We note that the single κ period ratio of the innermost pair of planets is sufficient to identify the required PO, since we seekthe PO in a small range around the nominal, remaining period ratios, X and X with fixed MMR factors C , C . As explainedbefore, in the MCOA-like optimisation the space of vectors Z is explored in a grid of pre-computed PO. Unfortunately, even inthe 5-dimensional space of Z , the number of solutions to be data–fitted becomes huge. Moreover, the mass anti-correlation withm e tending toward very small and non-realistic values implies a difficulty in estimating parameters ranges of the grid as well asits resolution. Even if we reach the neighborhood of the merit function’s extremum, it is difficult to find and tune the proper gridresolution for all parameters of the model, and the 5-dim PO grid must be updated in subsequent iterations. Still, the method isuseful in investigating the parameter space in wide ranges, and provides a good starting point (solutions) for a more refined andaccurate method.Clearly, a mathematically correct algorithm must explore the model parameter sub-space (a manifold) fixed by the requirementof PO. In order to implement this manifold fitting, for a given (prescribed) Z , a co-planar PO is searched for, resulting in theCartesian state vector X . Next, we select the best-fitting parameters as a vector Y , in the orbital and geometric elements space.In this way, for the given Z , the objective function, such as χ = χ ( Z , Z , Z , Z , Z ) is optimised under assumptions thati) the model orbits are periodic, ii) the POs are optimised through the linear scaling, time-phasing, its distance (parallax), andspatial orientation ( Y ). We closed the whole algorithm in a single procedure, being numerical implementation of the merit4 G O ´ZDZIEWSKI & M
IGASZEWSKI
Figure A5.
A PO–model to all HR 8799 observations in the literature (see the text for details). Big rotated rectangles are for the earliestHST data in (Lafrenière et al. 2009; Soummer et al. 2011), circles are for the uniformly reduced data set in (Konopacky et al. 2016), the starsymbol is for the GRAVITY measurement in (De Rosa et al. 2020), pentagons are for GPI data in (Wang et al. 2018), and diamonds are for theVLT/SPHERE, SUBARU and LBT data collected in (Wertz et al. 2017). Black curves illustrate synthetic orbits in the sky-plane, derived fromthe PO grid search, and plotted between epochs t = 1998 . and t = 2018 . of the GRAVITY datum, for ln L < − . , and the redcurves are for solutions providing ln L < − . , see Fig. A4. Subsequent panels are for close-ups of the data and orbital arcs for each planet,respectively. Note that ∆ RA is labeled negative w.r.t. the RA direction. function for the optimisation of Z . Both steps that consist of continuing the PO and its final fitting to the data, in the Y and Z spaces equivalent to the -dim vector of sampled parameters x = ( Z , Z , Z , Z , Z , Y , Y , Y , Y , Y , Y ) , explicitly x = ( m e , m d , m c , m b , κ, ρ, I, Ω , ω rot , t phase , Π) , are being done with a help of the Levenberg–Marquardt (LM) algorithm (Presset al. 2002). The iterative scheme enables us to find the best-fitting PO in a small number of steps, typically a few tens iterations.AI.4. DE-MC sampling and uncertainties of the best-fitting parameters
In order to assess realistic uncertainties, and to investigate possible parameter correlations, we performed the DifferentialEvolution Markov Chain (DE-MC) sampling (Ter Braak 2006). Recalling some well known elements of the Bayesian statisticsand the Markov Chain Monte Carlo sampling (e.g. Gregory 2010), we consider the posterior probability distribution P ( x |D ) ∼L ( D| x ) P ( x ) , where D denotes the data set, L ( D| x ) represents a probability that parameters x explain the data set D , and P ( x ) is the prior information imposed on x . We define the ln L function the same, as in Eq. A1, ln L ( D| x ) = − χ ( x ) − N obs (cid:88) i =1 [ln σ i,α + ln σ i,δ ] − N obs ln(2 π ) ≡ − χ ( x ) + const , (A3)but skipping the error floor, since we performed the MCMC sampling on the reduced data set D described in Sect. 2, and forthese measurements the best-fitting models yield χ ν ∼ – there is no need to account for the uncertainties correction.The DE-MC sampling, which is a variant of the canonical Metropolis-Hastings algorithm, occurs according to the probabilityof moving from a starting point x ( i )1 to a new point x ( i )2 in the parameter space, p ( x ( i )2 | x ( i )1 ) , which is a product of q ( x ( i )2 | x ( i )1 ) and α ( x ( i )1 , x ( i )2 ) , where q ( x ( i )2 | x ( i )1 ) is a probability of choosing a candidate point x ( i )2 when starting from x ( i )1 . The superscriptdenotes the i -th chain from a population of n = 100 chains which are evolved in parallel. The candidate point of the i -th chain is N EXACT L APLACE RESONANCE IN THE
HR 8799
PLANETARY SYSTEM − − ΔRA [arcsec] − − Δ D E C [ a r c s e c ] Figure A6.
Global view of the HR 8799 system geometry in the plane of the sky for all data in literature, and the astrometric model derived inthe PO–grid search, see caption to Fig. A5. chosen according to: x ( i )2 = x ( i )1 + γ (cid:16) x ( j )1 − x ( k )2 (cid:17) + Uniform (cid:16) − ∆ ( i ) , +∆ ( i ) (cid:17) , (A4)where the chains j and k ( j (cid:54) = k and j, k (cid:54) = i ) are chosen randomly, while ∆ ( i ) is chosen individually for each parameter. In orderto obtain (cid:39) acceptance rate we chose γ = 0 . and ∆ ( i ) were (10 − m Jup , . × − m Jup , . × − m Jup , × − m Jup , × − , × − yr, . × − , − mas, . × − deg , × − deg , × − deg) , for subsequent components of x (see the previous subsection). The second term in Eq. A4 is the Metropolis-Hastings ratio: α ( x ( i )1 , x ( i )2 ) = min (cid:20) , P ( x ( i )2 ) P ( x ( i )1 ) (cid:21) , (A5)which denotes the acceptance probability of x ( i )2 when starting from x ( i )1 . Importantly, the DE-MC algorithm propagates anumber of Markov chains in parallel, starting from different initial positions in the parameters space, and introduces mixingof the solutions in the chains through the Differential Evolution (Price et al. 2005). That makes this algorithm both simpleand computationally efficient. We also note that the DE-MC approach is crucial for our optimisation problem, given the need ofcomputationally complex PO continuation w.r.t. model parameters, since the PO cannot be updated sufficiently freely, as requiredby the Markov chain propagation.Priors P ( x ) for the masses were set as Gaussian with mean values and standard deviations according to Wang et al. (2018),from the hot-start evolutionary models, to (5 . ± . m Jup for HR 8799b, and (7 . ± . m Jup for all other planets. Similarly, theparallax Gaussian prior is
Π = (24 . ± . mas (Brown et al. 2018). For the remaining parameters, the prior distributionswere uniform, in sufficiently wide ranges.We initiated the DE-MC sampling by choosing solutions from the vicinity of the best-fitting model in Table 1. Theevolution of the whole population of Markov chains is illustrated in Fig. A7 with black curves, while one selected, examplechain is depicted with the red colour. At the beginning (first ∼ iteration steps) all the chains evolve closely to the initialcondition. Since the differences x ( j )1 − x ( k )2 increase, the sampling begins to occur over wider part of the parameter space. After ∼ steps the chain is already burnt-out. Those first steps were not included in the final statistics of solutions obtainedafter iterations. In this DE-MC experiment we did not estimate the auto-correlation time for the Markov chains, sinceclearly the relatively small number of iterations already leads to a smooth approximation of the posterior. Also, as illustrated6 G O ´ZDZIEWSKI & M
IGASZEWSKI
Figure A7.
The DE-MC sampling for the mass of HR 8799e. Black curves indicate the evolution of all individual chains, while the redcurve illustrates the behaviour of one chosen chain. First steps are treated as burn-in steps. The total number of iterations is , in thetop panel first steps were shown for clarity. in Fig. A7, each of the chains quickly reach the random-walk state and explore the whole parameter space. Remarkably, thisbehaviour is much different from the MCMC sampling with the full Keplerian or even the N -body models (e.g. Konopacky et al.2016; Wertz et al. 2017; Wang et al. 2018, GM18, and references therein), that notoriously exhibit parameter correlations andlong auto-correlation times ∼ , due to multi-modal posteriors and ill-constrained optimisation problem implied by a smallratio of the measurements to the number of free parameters and narrow time window of the data.According to the final results of the DE-MC sampling, as well as to the scans of χ function, the best-fitting configuration ismeaningfully constrained w.r.t. all parameters. In particular, the bottom-right panel in Fig. 3 (also Fig. A8) illustrates the Gaussianprior as GAIA parallax (Brown et al. 2018) (the red curve) over-plotted with the DE-MC posterior. The distributions closelyoverlap. We also recall the grid-based experiments indicating that the best-fitting parallax may be determined independently ofthe GAIA measurements. The -dim posterior probability distributions of all the free parameters determined with the DE-MCsampling are shown in Fig. A8, while -dimensional contour plots of the posteriors for the Keplerian elements are illustrated inFig. A9. The parameters uncertainties derived from the sampling are listed in Tab 1. AII.
MASSES W.R.T. DATA BIASESWith the improved PO algorithm described in Sect. AI.3, we systematically explored the χ minimum in 2-dimensional planetmass planes, without (Fig. A10) and with independently determined astrophysical priors (Fig. A11), regarding the measurementsset D with N obs = 65 observations, and also reducing it by particular GPI points (data set D from hereafter), as explained below.For a given, fixed point in selected 2-dim masses plane, the remaining two masses and the κ period ratio are optimised in termsof the best-fitting χ . In the χ optimisation, the mass priors may be included as additional terms in the χ function (Eq. A2): ( m i − m i, cooling ) / ( σ i, cooling ) , where i = e, d, c, b and m i, cooling denotes the planet i mass constraint (prior), while σ i, cooling isits σ uncertainty. We set the mass priors after (Wang et al. 2018), the same as in the DE-MC experiments. The χ -scans in themass planes for this enhanced model are shown in Fig. A11. We note that in the χ experiments the parallax was treated as a freeparameter with no prior.In the first experiment for data set D and without considering mass priors, we found the best-fitting m e (cid:46) . m Jup . Also, thebest-fitting and astrophysical masses are significantly different, as marked in the top row of Fig. A10, particularly for HR 8799d.In this figure, the mass estimates from the cooling theory are shown for a reference. Masses of HR 8799e and HR 8799c arestrongly anti-correlated and not bounded at all, since the best-fitting m e converges towards very small and non-realistic values.When fixing the inner planet’s mass at m Jup , the anti-correlation disappears, and masses m d and m b become much betterconstrained, but their values are still significantly shifted with respect to the astrophysical priors (Fig. A11, top row). N EXACT L APLACE RESONANCE IN THE
HR 8799
PLANETARY SYSTEM D set – rows from top to bottom are for subsequent planets.While the most precise GRAVITY datum is modelled apparently perfectly, there are precision GPI observations (Wang et al.2018) significantly deviating from the astrometric model, compared to the uncertainties. In the next experiment, we temporarilyremoved these points from the data set and we found a new best-fitting model for this modified set D , with residuals shown inthe right panel of Fig. A12. The GPI points, over-plotted with bigger grey symbols, reveal systematic shifts w.r.t. this best-fittingmodel.Figure A13 illustrates the residuals in the (RA, DEC)-plane. All the GPI points exhibit a systematic positive RA shift withrespect to the model, and apart from one point, all of them have negative DEC deviations (the top-left panel). Moreover, most ofthe data points deviate from the model by more than σ (the bottom-left panel in A13). Observations D , without the GPI data,are distributed uniformly in the (RA, DEC)–residuals plane, as expected for a statistically valid solution (the right column). Thismay suggest a bias in the GPI data w.r.t. the other measurements, yet of an unknown origin.The obtained χ minima with mass priors for data set D , presented in top row of Fig. A11 overlap with the results of the DE-MCsampling around the best-fitting model illustrated in Figs. A8 and A9. As noted above, the best-fitting mass of planet HR 8799d isthe only one significantly inconsistent with the priors from thermodynamical tracks by (cid:39) m Jup (top row of Fig. A11). However,when the GPI measurements are excluded from the data set, then the difference reduces by factor of (cid:39) , making the astrometricmodel results marginally consistent with the HR 8799d mass determined from the cooling theory (bottom row of Fig. A11). Allmasses become constrained much better for the reduced data set D than D and are marginally consistent with the astrophysicalvalues, although their uncertainties are still significant. This experiment demonstrates the sensitivity of the astrometric model tothe most accurate data points. We also recall, that with added just one GRAVITY-like measurement for each planet close to thepresent epoch, the astrometric data alone might fully constrain the masses (see the main part, and Fig. A1). AIII.
THE DEBRIS DISCS SIMULATIONAIII.1.
The N -body model of the debris discs Given the ongoing discussion in the literature, as summarized in the main part, we aim to resolve the dynamical structureof the debris disks composed of small, Kuiper-belt like objects. Such a structure may reflect unique characteristics implied bythe strictly resonant motion of the planets. Here, we essentially follow the approach in GM18. The numerical model relies indetermining the orbital stability of small-mass particles in the HR 8799 system through resolving the chaotic or regular characterof their motion with the MEGNO (cid:104) Y (cid:105) fast indicator (Cincotta et al. 2003). We dubbed it the (cid:104) Y (cid:105) -model. As we found in GM18with the long-term, direct N -body integrations, the (cid:104) Y (cid:105) -model reproduces closely the dynamical structure of the debris discsfound with the direct integrations, yet in much shorter computation time.Here, we conducted an extensive (cid:104) Y (cid:105) -model simulation of the debris discs co-planar with the planets involved in the exactLaplace resonance (Tab. 1). We considered three mixed fractions of asteroids with masses of − m Jup , similarly to GM18, aswell as − m Jup and − m Jup . As the initial Keplerian osculating elements, we randomly draw the semi-major axis a ∈ [10 , au, the pericenter longitude and the mean anomaly (cid:36) , M ∈ [0 ◦ , ◦ ) . For the inner part of the disk ( a < au),we selected e ∈ [0 , . , and for the outer part, beyond planet HR 8799b, e < − ( a /a b )(1 + e b ) , i.e., under the collisioncurve of asteroids with HR 8799b in the ( a , e ) -plane. We integrated the equations of motion and the variational equations forthe whole N -body system of the observed planets (Table 1, primaries), updated by a test particle, with the Bulirsh-Stoer-Gragg(BGS) integrator. The local and absolute accuracy of the integrator set to (cid:15) = 10 − provided the relative energy error as small as − for the total integration time of 10 Myrs. The BGS algorithm has been proved reliable for collisional and chaotic dynamics,which may be anticipated on the basis of previous works (GM18).Concerning the appropriate integration time required to reliably characterize the orbits of asteroids, we note that the four,massive planets are locked deeply in the Laplace resonance (Fig. A2), and each planet is located in the center of the stability zone(Fig. A3). The system stability is robust to perturbations of quite massive additional companions (see also simulations in GM18for “asteroid” masses as large as – m Jup ). Therefore, the (cid:104) Y (cid:105) integrations of the best-fitting initial condition extended by theelements of a test asteroid reveal the dynamical character of its motion, and the orbits of the primaries are not affected.The geometric structure of the debris discs is illustrated in Figs. 4, A14 and A15. In the numerical experiment, we collected (cid:39) . × (cid:104) Y (cid:105) -stable orbits. Astrocentric positions of asteroids are marked at at the end of integration time (top-left panelof Fig. A15) and at the initial epoch (top-right panel in Fig. A15), and colour-coded, according to their osculating eccentricity.Such snapshots represents a population of quasi-periodic and resonant orbits of the asteroids with various orbital phases andeccentricity, while their semi-major axes may overlap. We note, following GM18, that the orbits might be potentially present in8 G O ´ZDZIEWSKI & M
IGASZEWSKI e (m Jup ) 0 5 10 6 7 8 9 10 0 5 10 6 7 8 9 10m d (m Jup ) 0 5 10 5 6 7 8 9 10 0 5 10 5 6 7 8 9 10m c (m Jup ) 0 5 10 4 5 6 7 8 0 5 10 4 5 6 7 8m b (m Jup ) 0 5 10 1.97 1.98 1.99 2 0 5 10 1.97 1.98 1.99 2 κ phase (yr) 0 5 10 15 1.05 1.06 0 5 10 15 1.05 1.06 ρ Π (mas) 0 5 10 25 26 27 28 0 5 10 25 26 27 28I (deg) 0 5 10 60 61 62 63 64 0 5 10 60 61 62 63 64 Ω (deg) 0 5 10 154 156 158 160 0 5 10 154 156 158 160 ω rot (deg) posteriorGaussianprior 1 σ σ σ σ σ σ po s t e r i o r s a m p l e s ( pe r c en t ) Figure A8.
Posterior probability distributions for the free parameters of the PO model of the HR 8799 system (grey curves). The shades ofgrey are for σ, σ and σ confidence intervals. Red curves mark the Gaussian priors set for out of model parameters. The red areasunder the curves indicate σ ranges of each parameter. The remaining parameters of the model have the uniform (non-informative) priors,not shown in the plot. See Table 1 and the text in Appendix. the real system, but the actual population of asteroids may depend on the formation history of the whole system, its migrationhistory, as well as locally variable density of asteroids.In regions interior to, and beyond orbit of planet HR 8799b, the test orbits are extremely chaotic besides particular resonantsolutions. Such (cid:104) Y (cid:105) -unstable orbits are also strongly unstable in the Lagrangian (geometric) sense — particles are ejected orcollide with the primaries in the time scale of a few Myrs only. We found this after testing the semi-major axis–eccentricityevolution in time for orbits selected in a strip of 1000 initial conditions marked with red filled points in the right panel ofFig. A14. It shows the proper (canonical) elements (Morbidelli 2002) of dynamically stable asteroids in the semi-major axis–eccentricity plane ( a , e ) , see the right panel. In order to study unstable motions, test particles were randomly placed underthe collision curve with planet HR 8799b. The initial eccentricity of their orbits is slightly larger than the respective limit of (cid:104) Y (cid:105) -stable motions, and the initial semi-major axes a ∈ [100 , au as well as initial phases are also random. We investigatedclosely orbits of all these test asteroids by integrating them for 10 Myrs. In this set, 466 asteroids collided with planet HR 8799b,128 objects collided with the star, and 351 asteroids were ejected from the system beyond 5000 au, leaving the radius of autypically in a few Myrs, and less than the maximum interval of 10 Myrs. Only (cid:39) objects located in stable, resonant regionssurvived for the maximum integration interval.Moreover, with the Modified Fourier Transform or the Fundamental Frequency Analysis (Šidlichovský & Nesvorný 1996) ofthe canonical Jacobi elements z i = e i (cos (cid:36) i + √− (cid:36) i ) , ( i = b , c , d , e), illustrated in the bottom-left panel in Fig. A2, wecomputed the frequency spectrum of planet pericenters rotation ˙ (cid:36) ( t ) . Since the motion of the planets is strictly periodic, the N EXACT L APLACE RESONANCE IN THE
HR 8799
PLANETARY SYSTEM Figure A9.
Posterior probability distributions presented in diagrams of semi-major axis vs. eccentricity (left four panels), as well as theargument of pericentre vs. the mean anomaly (right four panels). The shades of grey, from darkest to lightest, represent σ, σ and σ confidence levels, respectively. z i ( t ) signals involve common leading frequency f (cid:36) = − . (cid:48)(cid:48) yr − equivalent to the retrograde rotation of the system withthe period P (cid:36) (cid:39) . years, i.e., only (cid:39) orbits of planet HR 8799b. Besides the leading frequency there are a few evenlarger, with periods smaller than years.Since the dynamics is governed by short-term MMRs and, possibly, secular resonances, we could fix the same integrationtime of Myr across the whole disk. That integration time corresponds to (cid:39) (cid:39) au which is sufficient to resolve the dynamical character of asteroid orbits. Particles markedas (cid:104) Y (cid:105) -stable for that interval of time should persist for more than 10 times longer interval in Lagrange-stable orbits, roughly – Myr (see GM14, GM18), which is much longer than typical estimates of the parent star lifetime, (cid:39) Myrs (Wang et al.2018), in the 30–60 Myr range earlier adopted in (Marois et al. 2010). At the outer edge of the disk (cid:39) au, as determined byBooth et al. (2016), the integration interval translates to a few thousands of orbital periods, which is still meaningful to determinethe stability border in the ( a , e ) -plane, as we justified above. Moreover, given strong instability generated by the short terminteractions, the (cid:104) Y (cid:105) integrations may be stopped as soon as (cid:104) Y (cid:105) > , sufficiently different from (cid:104) Y (cid:105) (cid:39) for stable systems.That makes it possible to examine large sets of a few test orbits, orders of magnitude larger than they could be sampled withthe direct N -body integrations. The most complex and interesting parts of the debris discs may be then mapped in detail with the (cid:104) Y (cid:105) -model. AIII.2. Dynamical structure and features of the debris discs
Regarding the inner part of the system, we found the same irregular inner boundary of the outer disk, similarly to simulations inGM18. In order to understand this feature, we analyse the ( a , e ) –diagram, shown in Fig. A14. Comparing the left-hand panelof Fig. A14 with the disk structure illustrated in Fig. 4, we find that the inner edge of the outer disk is significantly asymmetric dueto low and moderate eccentricity orbits in the 1:1 and 3:2 MMR with planet HR 8799b. Low density of asteroids around (cid:39) auappears due to unstable 2:1 MMR and higher order resonances forming a kind of thickening “comb” with increasing semi-majoraxis. It forms a border of stable orbits shifted below the collision curve with planet HR 8799b by a substantial value of ∼ . . Wecan now understand and interpret the strongly unstable orbital evolution of the tested asteroids in this zone (Fig. A14, the rightpanel). The strong instability is caused by overlapping two-body MMRs, multi-body MMRs and (possibly) mixed secular–meanmotion resonances. The pericenter frequency of the system is commensurate with the mean motion of asteroids n in this zone,for instance, regarding absolute values of the frequencies, f (cid:36) : 1 n at (cid:39) au, f (cid:36) : 2 n at (cid:39) au, f (cid:36) : 1 n at (cid:39) au, f (cid:36) : 1 n at (cid:39) au. However, the resonances are retrograde for the disk rotating with the same spin direction, as the planets,therefore we did not observe their direct or clear dynamical influence on the asteroids. A streaking feature of stable zone beyondHR 8799b is the presence of low density rings, which could be identified with higher-order resonances with this outermost planet,such as 2:1, 3:2, 3:1 and 5:2, extending up to (cid:39) au (Fig. A14 and Fig. 4 in the main part).0 G O ´ZDZIEWSKI & M
IGASZEWSKI
Figure A10. A χ -scan of PO models in the plane of masses HR 8799e and HR 8799c (left panel), and HR 8799d and HR 8799b (rightpanel) for the data set with (top row, D ) and without the GPI measurements (bottom row, D ). The red filled symbols and the circle/ellipsemark reference astrophysical masses m e = (7 . ± . m Jup and m c = (7 . ± . m Jup (left panel) and m d = (7 . ± . m Jup andm b = (5 . ± . m Jup (right panel), following (Wang et al. 2018). The green filled point denotes the position of the minimum of χ function,while the green contour denotes the level of min χ + 1 . The white and black curves denote the levels of min χ + 2 , . . . , +6 , apart from thebottom-left panel for which the white contours denote the confidence levels of min χ + 0 . , +0 . , +0 . . A stable 1:1 MMR with planet HR 8799b forms huge, symmetric Lagrangian areas of low eccentricity objects extending for70–80 au and (cid:39) au across. The Lagrangian 1:1 MMRs governed by inner planets are non-symmetric in respective pairs. Thereare also islands of the 2:1 MMR and 3:2 MMR with HR 8799d and HR 8799e. In these islands, eccentricity of the asteroidsreaches e (cid:39) . (yellow colour in Fig. 4). The outer, continuous edge of the inner debris disk appears at (cid:39) au. (The dynamicalstructure of the inner disk was investigated in more detail in GM18).In the top panels of Fig. A15, we present the global view of the debris discs revealed by (cid:39) . × (cid:104) Y (cid:105) -stable orbits in thewhole simulation. Similarly to Fig. 4, the panels represent snapshots of astrocentric coordinates ( x, y ) of asteroids and theirosculating orbital eccentricities e (color-coded and labeled in the top bar), at the initial epoch (the right panel) and at the end ofthe integration interval of 10 Myrs (the left panel). We selected the end epoch in order to illustrate a saturation of asteroids after asubstantial interval of thousands of orbital periods. Initial positions of the planets are marked with filled circles. For a reference,gray rings illustrate their orbits integrated for the same interval of 10 Myr, with the initial conditions in Table 1, independently ofthe disk integrations.The top-left panel of Fig. A15 shows the debris disks in the orbital plane at the final epoch t = 10 Myr, and the top-right panelis for the sky-view of the disks at the initial epoch t , rotated by the inclination and nodal angle in the initial condition (Tab. 1), N EXACT L APLACE RESONANCE IN THE
HR 8799
PLANETARY SYSTEM Figure A11. A χ -scan of PO models in the masses plane of HR 8799e and HR 8799c (left panel), and HR 8799d and HR 8799b (rightpanel) for the data set with (top row, D ) and without the GPI measurements (bottom row, D ), as well as with the the hot-start cooling theorypriors. The red filled symbols and the circle/ellipse correspond to m e = (7 . ± . m Jup and m c = (7 . ± . m Jup (left panel) andm d = (7 . ± . m Jup and m b = (5 . ± . m Jup (right panel), following Wang et al. (2018). The green filled point denotes the position ofthe minimum of χ function, while the green contour denotes the level of min χ + 1 . The white and black curves denote the confidence levelsof min χ + 2 , . . . , +6 . respectively. These global representations for the outer disk reveal a ring of high-eccentric orbits between (cid:39) au and (cid:39) auand broad outer ring forming a diffuse outer edge of the disk. We note that the inner ring is substantially shifted with respect tothe inner edge of the disk found at (cid:39) au. The dynamical structure of the whole disk is also illustrated in the ( a , e ) -planeof the canonical Poincaré elements in the right plot of Fig. A14. In the top-right panel of Fig. A15 we also marked the innerand outer boundary of the disk model in Booth et al. (2016), according with their estimate of the inclination I = 41 ◦ and nodalangle Ω = 50 ◦ . These values appear substantially different from the inclination and nodal angle of the best-fitting elements ofthe planetary system orbital plane (Tab. 1).In order to interpret the ring structure around (cid:39) au we plotted (not shown here) the canonical, osculating eccentricity e vs. the astrocentric radius of particles at the epoch t , and also at t + 10 Myr. They reveal that the excess of particles with higheccentricity seems to be a real feature, unlikely due to a particular sampling or plotting order of the particles. It is also clear that e is a very steep function of the radius r at the innermost part of the outer disk.Given the variation of eccentricity across the disk, we computed the Keplerian velocity dispersion of the particles. We binnedasteroids in the region covering the whole disk, x, y ∈ [ − , au in square bins of au × au. In each box, with non-zeronumber of particles, we computed σ v = (cid:80) ni ( v i − v ) /n, where v i is the velocity module of a particle i in the given bin, n is the2 G O ´ZDZIEWSKI & M
IGASZEWSKI
Figure A12.
Residuals of the best-fitting model to the data set with (left panel, D ) and without the GPI observations (right panel, D ). In theright-hand panel the GPI data, not included in the fitted data set, are marked with big grey symbols without error bars. The GRAVITY datumfor HR 8799 e is additionally enlarged in the top row (grey rectangles). counted number of particles in this bin and v is the mean velocity module. The results are illustrated in the bottom-left panel ofFig. A15. The ring structure associated with high eccentricity asteroids and the gradient of e ( r ) implies the velocity dispersion σ v a few times larger than in the inner parts of the disk. It could imply more intense dust production due to both locally largerdensity of objects and higher velocity during their collisions. We may note that the inner disk boundary fitted by Booth et al.(2016) seems to overlap with the eccentricity ring edge, which could suggest a systematic shift of the detected emission w.r.t.the actual dynamical border of the disk at (cid:39) au. It might actually confirm the results of Wilner et al. (2018), in their morerecent model of the disk predicting the inner edge also at (cid:39) au. Such border is better consistent with our updated orbitalmodel of the HR 8799 system, regarding the present parallax estimate in the GAIA DR2 catalogue, and the resulting, true lineardimensions of the system.Finally, we simulated the relative intensity image of the disk. The relative intensity is defined the same as in (Read et al. 2018), I ν ( r ) ∼ K Σ( r ) r − / , where Σ( r ) is the surface density and K is the scaling factor. In order to estimate Σ( r ) , we used the countsof asteroids in the same × au bins used for computing the velocity dispersion. The results are illustrated in the bottom-rightpanel in Fig. A15. The bright rings are associated with fractions of stable asteroids in the 3:2 MMR and 2:1 MMR with planetHR 8799b.While interpretation of the results needs more work, we might briefly conclude that the disk simulation reveals features relatedto the resonant character of the system. They consits of asymmetry of the inner edge of the outer debris disk, and highly variabledensity of asteroids in its inner part, due to low-order MMRs with planet HR 8799b, including large Lagrangian clouds. Thereare also possible two rings of high-eccentricity asteroids around – au and at the outer edge (cid:39) au. These features mayinfluence the intensity images used for modeling the emission in different wavelengths, and likely they should be accounted forin order to avoid biases in the emission models. N EXACT L APLACE RESONANCE IN THE
HR 8799
PLANETARY SYSTEM Figure A13.
Residuals of the best-fitting model without the mass priors and without the GPI data. The top panels show the residuals togetherwith their uncertainties, while the bottom panels show the residuals weighted with the uncertainties. Circles of radii equal to and marks σ and σ deviations of a given data point from the model. In the right-hand column, the fitted astrometric measurements are illustrated.Deviations of the GPI data from the model are presented in the left-hand panels. Figure A14.
Canonical, Poincaré elements ( a , e ) of (cid:104) Y (cid:105) -stable solutions at the end of the integration interval of 10 Myr. Grey lines are forthe collision curve of orbits with planet HR 8799b. Approximate positions of a few low-order MMRs with planets HR 8799b and HR 8799c arelabeled. The left-hand panel is for the inner part of the disk, and the right panel is for the whole simulation. Red filled circles in the right-handpanel illustrate 1,000 initial condition of test orbits, analysed in order to explain wide instability zone below the collision curve. O ´ZDZIEWSKI & M
IGASZEWSKI
Figure A15.
Top-left panel: the global view of debris disk revealed by (cid:39) . × (cid:104) Y (cid:105) -stable orbits of the whole simulation illustrated asa snapshot of astrocentric coordinates ( x, y ) of asteroids at the end of integration interval of 10 Myrs. Osculating orbital eccentricities e ofthese orbits are color-coded and labeled in the top bar. Initial positions of planets are marked with filled circles. Gray rings illustrate their orbitsintegrated for 10 Myr. Top-right panel: similar to Fig. 4 and the top-left panel, but the initial astrocentric ( x, y ) -coordinates of asteroids in (cid:104) Y (cid:105) -stable orbits are rotated by the inclination and nodal angles of the initial condition in Tab. 1. Gray ellipses illustrate the disk boundaries r in = 145 au and r out = 429 au, respectively, fitted in (Booth et al. 2016), and rotated by inclination I = 41 ◦ and nodal angle Ω = 50 ◦ derivedin this paper. Bottom-left panel: velocity dispersion in the debris discs evaluated at au × au bins and color-coded over initial astrocentric ( x, y ) -coordinates of asteroids in (cid:104) Y (cid:105) -stable orbits. Bottom-right: The relative Planck intensity of the outer disk, I ν ∼ K Σ( r ) r − / , as inRead et al. (2018), where Σ( r ) is proportional to the number of asteroids in au ×2