# Orbital and dynamical analysis of the system around HR 8799. New astrometric epochs from VLT/SPHERE and LBT/LUCI

A. Zurlo, K. Goździewski, C. Lazzoni, D. Mesa, P. Nogueira, S. Desidera, R.

Gratton, F. Marzari, M. Langlois, E. Pinna, et al.

► To cite this version:

A. Zurlo, K. Goździewski, C. Lazzoni, D. Mesa, P. Nogueira, et al., Orbital and dynamical analysis of the system around HR 8799. New astrometric epochs from VLT/SPHERE and LBT/LUCI. Astronomy and Astrophysics - A&A, 2022, 666, 10.1051/0004-6361/202243862. insu-03854051

HAL Id: insu-03854051

https://insu.hal.science/insu-03854051

Submitted on 15 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Orbital and dynamical analysis of the system around HR 8799
## New astrometric epochs from VLT/SPHERE and LBT/LUCI

A. Zurlo<sup>1,2,3</sup>①, K. Goździewski<sup>4</sup>, C. Lazzoni<sup>5,6</sup>, D. Mesa<sup>6</sup>, P. Nogueira<sup>1</sup>, S. Desidera<sup>6</sup>, R. Gratton<sup>6</sup>, F. Marzari<sup>7</sup>, M. Langlois<sup>8,3</sup>, E. Pinna<sup>9</sup>, G. Chauvin<sup>10</sup>, P. Delorme<sup>10</sup>, J. H. Girard<sup>11</sup>, J. Hagelberg<sup>12</sup>, Th. Henning<sup>13</sup>, M. Janson<sup>14</sup>, E. Rickman<sup>15</sup>, P. Kervella<sup>16</sup>, H. Avenhaus<sup>13</sup>, T. Bhowmik<sup>1</sup>, B. Biller<sup>17</sup>, A. Boccaletti<sup>16</sup>, M. Bonaglia<sup>9</sup>, M. Bonavita<sup>18,6</sup>, M. Bonnefoy<sup>10</sup>, F. Cantalloube<sup>3</sup>, A. Cheetham<sup>12</sup>, R. Claudi<sup>6</sup>, V. D'Orazi<sup>6</sup>, M. Feldt<sup>13</sup>, R. Galicher<sup>16</sup>, E. Ghose<sup>9</sup>, A.-M. Lagrange<sup>10</sup>, H. le Coroller<sup>3</sup>, R. Ligi<sup>20</sup>, M. Kasper<sup>21</sup>, A.-L. Maire<sup>10</sup>, F. Medard<sup>10</sup>, M. Meyer<sup>19</sup>, S. Peretti<sup>12</sup>, C. Perrot<sup>16</sup>, A. T. Puglisi<sup>9</sup>, F. Rossi<sup>9</sup>, B. Rothberg<sup>22,23</sup>, T. Schmidt<sup>16</sup>, E. Sissa<sup>6</sup>, A. Vigan<sup>3</sup>, and Z. Wahhaj<sup>24</sup>

*(Affiliations can be found after the references)*

Received 25 April 2022 / Accepted 15 June 2022

## ABSTRACT

**Context.** HR 8799 is a young planetary system composed of four planets and a double debris belt. Being the first multi-planetary system discovered with the direct imaging technique, it has been observed extensively since 1998. This wide baseline of astrometric measurements, counting over 50 observations in 20 years, permits a detailed orbital and dynamical analysis of the system.

**Aims.** To explore the orbital parameters of the planets, their dynamical history, and the planet-to-disk interaction, we made follow-up observations of the system during the VLT/SPHERE guaranteed time observation program. We obtained 21 observations, most of them in favorable conditions. In addition, we observed HR 8799 with the instrument LUCI at the Large Binocular Telescope (LBT).

**Methods.** All the observations were reduced with state-of-the-art algorithms implemented to apply the spectral and angular differential imaging method. We re-reduced the SPHERE data obtained during the commissioning of the instrument and in three open-time programs to have homogeneous astrometry. The precise position of the four planets with respect to the host star was calculated by exploiting the fake negative companions method. We obtained an astrometric precision of the order of 6 mas in the worst case and 1 mas in the best case. To improve the orbital fitting, we also took into account all of the astrometric data available in the literature. From the photometric measurements obtained in different wavelengths, we estimated the masses of the planets following the evolutionary models.

**Results.** We obtained updated parameters for the orbits with the assumption of coplanarity, relatively small eccentricities, and periods very close to the 2:1 resonance. We also refined the dynamical mass of each planet and the parallax of the system (24.49 ± 0.07 mas), which overlap with the recent Gaia eDR3/DR3 estimate. Hydrodynamical simulations suggest that inward migration of the planets caused by the interaction with the disk might be responsible for the planets being locked in resonance. We also conducted detailed N-body simulations indicating possible positions of a putative fifth planet with a mass below the present detection limits of ≈3 $M_{\text{Jup}}$.

**Key words.** planets and satellites: dynamical evolution and stability – planet-disk interactions – stars: individual: HR8799 – instrumentation: adaptive optics – astrometry – techniques: image processing

## 1. Introduction

Among the exoplanets discovered with the high-contrast imaging technique (e.g., Chauvin et al. 2004; Lagrange et al. 2010; Rameau et al. 2013; Keppler et al. 2018; Bohn et al. 2021), the system around HR 8799 is undoubtedly one of the most interesting. This is mainly because HR 8799 hosts a greater number of planets detected with high-contrast imaging than any other system, only three systems that host two planets were detected: PDS 70 (Keppler et al. 2018; Haffert et al. 2019), TYC 8998-760-1 (Bohn et al. 2020), and β Pic (Lagrange et al. 2010; Nowak et al. 2020). HR 8799 is a perfect laboratory with which to study dynamical interaction in young planetary systems, with its four planets (HR 8799 bcde; Marois et al. 2008, 2010) and a double debris belt (see, e.g., Su et al. 2009; Hughes et al. 2011; Matthews et al. 2014; Booth et al. 2016; Faramaz et al. 2021). The host star HR 8799 is a young (~42 Myr,

Hinz et al. 2010; Zuckerman et al. 2011; Bell et al. 2015), γ Dor-type variable star (Gray & Kaye 1999) with $\lambda$ Boo-like abundance patterns. The mass of the star is $1.47_{-0.17}^{+0.12} M_{\odot}$ (Sepulveda & Bowler 2022) and its distance is 40.88 ± 0.08 pc from Gaia measurements (Gaia Collaboration 2020).

This system is an optimal target for high-contrast imaging observations, as its four planets are easily detectable with state-of-the-art imagers; their contrasts are about 2–8 × 10<sup>-6</sup>; and the separation of the closest planet HR 8799 e is ~390 mas, which is further than the inner working angle of most current instruments designed for direct imaging. For these reasons, the system has been observed dozens of times starting from 1998, with different instruments, wavelengths, and configurations. Thanks to this rich pool of archival data, HR 8799 can be studied in detail from an astrometric point of view, being the only multi-planetary system for which there are tens of astrometric data points on a baseline of more than 20 yr. These measurements began with
uncertainties of ~20 mas, but have since reached unprecedented precision of below a milliarcsecond (0.1–0.2 mas) with optical interferometry (GRAVITY Collaboration 2019).

Works regarding the analysis of the dynamical interaction of the four planets have been presented; most of these agree on the near coplanarity of the planets, and that they are locked in a 1:2:4:8 resonance, (see, e.g., Fabrycky & Murray-Clay 2010; Esposito et al. 2013; Konopacky et al. 2016; Zurlo et al. 2016; Wang et al. 2018; GRAVITY Collaboration 2019; Goździewski & Migaszewski 2020). From the age of the system and the luminosity of the planets, evolutionary hot-start models derive masses of around 5–7 $M_{Jup}$ (Marois et al. 2010; Currie et al. 2011; Sudol & Haghighipour 2012). These masses are compatible with the dynamical studies because small masses help the orbits to remain stable (e.g., Wang et al. 2018). On the other hand, Brandt et al. (2021) found a dynamical mass for the innermost planet HR 8799 e of $9.6_{-1.8}^{+1.9} M_{Jup}$ assuming that planets c, d, and e share the same mass within ~20%. This result is 2 $M_{Jup}$ higher than the prediction from the evolutionary models.

In this paper, we present all the astrometric measurements obtained during the whole guaranteed time observation (GTO) of VLT/SPHERE for HR 8799 bcde. Ad hoc observations were designed for the orbital follow-up of the four planets; HR 8799 was observed 21 times in total with SPHERE, and only four observations were rejected based on quality criteria. To complement the SPHERE measurement, we obtained one astrometric epoch with the instrument LUCI, installed at the LBT. This instrument observed HR 8799 as part of the commissioning phase of the new adaptive optics (AO) system. With a total of 18 epochs and the addition of all the astrometric measurements available in the literature (69 astrometric epochs in total), we performed the dynamical analysis of the system and the planet-to-disk interaction.

The outline of the paper is as follows: in Sect. 2, we present SPHERE and LUCI observations; in Sect. 3, we describe the reduction methods applied and the astrometric results that we obtained. In Sect. 4, we present the astrometric fitting for the four planets of HR 8799 and in Sect. 5 a possible interpretation of the history of the system. We also explored the possibility of the presence of a fifth planet, studying the planets–disk interaction (Sect. 6). We provide our conclusions in Sect. 7.

## 2. Observations

### 2.1. VLT/SPHERE

HR 8799 was observed several times during the SpHere INfrared survey for Exoplanets (SHINE, papers I, II, III; Desidera et al. 2021; Vigan et al. 2021; Langlois et al. 2021) during the GTO of VLT/SPHERE. The instrument SPHERE (Beuzit et al. 2019) is a planet finder equipped with an extreme AO system (SAXO; Fusco et al. 2006; Petit et al. 2014) to characterize substellar companions with high-contrast imaging. The near-infrared (NIR) arm includes the IR dual-band imager and spectrograph (IRDIS; Dohlen et al. 2008) and an integral field spectrograph (IFS; Claudi et al. 2008). During the observations, these two subsystems observed the target in parallel. HR 8799 was periodically observed for astrometric monitoring, the setup of the observations included three different filter pairs: IRDIS in H2H3 bands ($\lambda_{H2} = 1.593$ µm, $\lambda_{H3} = 1.667$ µm), in BB_H ($\lambda_H = 1.625$ µm), and K1K2 bands ($\lambda_{K1} = 2.102$ µm, $\lambda_{K2} = 2.255$ µm). The coronagraph used for the shortest wavelengths was an apodized Lyot with a mask diameter of 185 mas and an inner working angle (IWA) of 0.09 (see Boccaletti et al. 2008), while for K1K2 band

<table><caption>**Table 1.** Summary of the observations of HR 8799 with IRDIS and IFS during SHINE.</caption><thead><tr><th>UT date</th><th>Filter</th><th>Δ<sub>FoV</sub> (°)</th><th>S/N<sub>Min</sub></th><th>S/N<sub>Max</sub></th></tr></thead><tbody><tr><td>2014-07-12</td><td>DB_H23</td><td>17</td><td>9</td><td>44</td></tr><tr><td>2014-08-13 <sup>(a)</sup></td><td>BB_J</td><td>23</td><td>0</td><td>36</td></tr><tr><td>2014-12-04 <sup>(b)</sup></td><td>BB_H</td><td>9</td><td>14</td><td>24</td></tr><tr><td>2014-12-05 <sup>(b)</sup></td><td>BB_H</td><td>8</td><td>11</td><td>25</td></tr><tr><td>2014-12-06 <sup>(b)</sup></td><td>BB_H</td><td>8</td><td>15</td><td>25</td></tr><tr><td>2014-12-08 <sup>(a,b)</sup></td><td>BB_H</td><td>8</td><td>6</td><td>25</td></tr><tr><td>2015-07-03</td><td>DB_K12</td><td>18</td><td>12</td><td>75</td></tr><tr><td>2015-07-29 <sup>(a,b)</sup></td><td>DB_J23</td><td>82</td><td>4</td><td>51</td></tr><tr><td>2015-07-30 <sup>(b)</sup></td><td>DB_K12</td><td>81</td><td>21</td><td>123</td></tr><tr><td>2015-09-27</td><td>DB_K12</td><td>24</td><td>10</td><td>71</td></tr><tr><td>2016-11-17</td><td>DB_H23</td><td>17</td><td>15</td><td>52</td></tr><tr><td>2017-06-14</td><td>DB_H23</td><td>19</td><td>14</td><td>67</td></tr><tr><td>2017-10-07 <sup>(a)</sup></td><td>BB_H</td><td>76</td><td>0</td><td>0</td></tr><tr><td>2017-10-11</td><td>BB_H</td><td>74</td><td>45</td><td>90</td></tr><tr><td>2017-10-12</td><td>BB_H</td><td>78</td><td>24</td><td>100</td></tr><tr><td>2018-06-18</td><td>DB_H23</td><td>34</td><td>20</td><td>78</td></tr><tr><td>2018-08-17</td><td>BB_H</td><td>73</td><td>29</td><td>57</td></tr><tr><td>2018-08-18</td><td>BB_H</td><td>73</td><td>42</td><td>106</td></tr><tr><td>2019-10-31</td><td>DB_K12</td><td>24</td><td>22</td><td>62</td></tr><tr><td>2019-11-01</td><td>DB_K12</td><td>54</td><td>38</td><td>114</td></tr><tr><td>2021-08-20</td><td>DB_H23</td><td>20</td><td>21</td><td>59</td></tr></tbody></table>

Notes. <sup>(a)</sup>Discarded for bad weather. <sup>(b)</sup>IRDIS alone.

the IWA is ~0.12. For a detailed description of the observing sequence, we refer the reader to Zurlo et al. (2014, 2016). In general, the working sequence includes an image of the off-axis star point spread function (PSF) for the flux calibration, a long coronagraphic sequence with the satellite spots (Langlois et al. 2013) mode, a second image of the stellar PSF, and the sky images. The waffle mode was used on purpose to assure maximum astrometric precision. The long waffle sequence was taken in pupil stabilized mode in order to apply the angular differential imaging (ADI; Marois et al. 2006) method. While IRDIS has a field of view of ~11 × 11", IFS is smaller (1.7" × 1.7"), and only the two interior planets are visible. The SHINE observations are summarized in Table 1. For almost all the epochs, the observation was with IRDIS and IFS working in parallel. On a few occasions, IRDIS was used alone; these observations are marked in the table. We discarded data for time periods where the conditions did not permit a clear detection of the planets or there was a very poor signal-to-noise ratio (S/N). In particular, we rejected data from 2014 August 13 (only IFS data were considered), 2014 December 08, 2015 July 29, and 2017 October 07. All the other observations are taken into account in this analysis.

### 2.2. Large Binocular Telescope/SOUL-LUCI

HR 8799 was observed with LBT/SOUL-LUCI1 during its commissioning on night 2020 September 29. LUCI1 (Seifert et al. 2010) is a NIR spectro-imager that can work in a diffraction-limited regime with a sampling of 15 mas pix$^{-1}$. The AO
<table><caption>**Table 2.** List of astrometric points from IRDIS observations.</caption><thead><tr><th>Date</th><th>Planet b<br>ΔRA ; ΔDec (mas)</th><th>Planet c<br>ΔRA ; ΔDec (mas)</th><th>Planet d<br>ΔRA ; ΔDec (mas)</th><th>Planet e<br>ΔRA ; ΔDec (mas)</th></tr></thead><tbody><tr><td>2014-07-12</td><td>1570 ± 3; 704 ± 3</td><td>-521 ± 3; 790 ± 9</td><td>-391 ± 2; -530 ± 2</td><td>-387 ± 2; -10 ± 3</td></tr><tr><td>2014-12-04</td><td>1574 ± 3; 701 ± 2</td><td>-514 ± 3; 798 ± 4</td><td>-399 ± 4; -525 ± 4</td><td>-389 ± 8; 11 ± 4</td></tr><tr><td>2014-12-05</td><td>1574 ± 4; 701 ± 3</td><td>-512 ± 3; 798 ± 4</td><td>-400 ± 4; -523 ± 4</td><td>-390 ± 7; 12 ± 4</td></tr><tr><td>2014-12-06</td><td>1573 ± 3; 701 ± 3</td><td>-512 ± 3; 797 ± 4</td><td>-403 ± 4; -524 ± 4</td><td>-383 ± 8; 11 ± 4</td></tr><tr><td>2015-07-03</td><td>1579 ± 1; 694 ± 1</td><td>-498 ± 1; 806 ± 1</td><td>-417 ± 1; -517 ± 1</td><td>-383 ± 9; 33 ± 5</td></tr><tr><td>2015-07-30</td><td>1580 ± 5; 689 ± 3</td><td>-495 ± 2; 806 ± 2</td><td>-419 ± 2; -516 ± 1</td><td>-386 ± 1; 36 ± 1</td></tr><tr><td>2015-09-27</td><td>1580 ± 1; 688 ± 1</td><td>-494 ± 1; 811 ± 1</td><td>-426 ± 1; -512 ± 1</td><td>-382 ± 9; 50 ± 5</td></tr><tr><td>2016-11-17</td><td>1589 ± 2; 666 ± 1</td><td>-464 ± 2; 824 ± 2</td><td>-454 ± 2; -490 ± 2</td><td>-378 ± 4; 90 ± 2</td></tr><tr><td>2017-06-14</td><td>1591 ± 1; 653 ± 1</td><td>-449 ± 1; 835 ± 1</td><td>-472 ± 2; -482 ± 2</td><td>-373 ± 3; 118 ± 2</td></tr><tr><td>2017-10-11</td><td>1595 ± 1; 647 ± 1</td><td>-441 ± 1; 839 ± 1</td><td>-480 ± 1; -477 ± 1</td><td>-369 ± 1; 128 ± 1</td></tr><tr><td>2017-10-12</td><td>1595 ± 1; 647 ± 1</td><td>-441 ± 1; 839 ± 1</td><td>-480 ± 1; -477 ± 1</td><td>-371 ± 2; 128 ± 2</td></tr><tr><td>2018-06-18</td><td>1601 ± 1; 635 ± 1</td><td>-424 ± 1; 848 ± 1</td><td>-497 ± 2; -463 ± 2</td><td>-358 ± 2; 156 ± 2</td></tr><tr><td>2018-08-17</td><td>1601 ± 2; 632 ± 3</td><td>-421 ± 1; 850 ± 1</td><td>-502 ± 1; -461 ± 1</td><td>-357 ± 1; 162 ± 2</td></tr><tr><td>2018-08-18</td><td>1600 ± 1; 632 ± 1</td><td>-421 ± 1; 851 ± 1</td><td>-502 ± 2; -458 ± 1</td><td>-358 ± 2; 163 ± 1</td></tr><tr><td>2019-10-31</td><td>1606 ± 2; 615 ± 2</td><td>-392 ± 2; 875 ± 2</td><td>-532 ± 3; -425 ± 2</td><td>-338 ± 2; 215 ± 2</td></tr><tr><td>2019-11-01</td><td>1611 ± 1; 611 ± 1</td><td>-388 ± 2; 870 ± 2</td><td>-530 ± 2; -430 ± 1</td><td>-335 ± 1; 210 ± 1</td></tr><tr><td>2021-08-20</td><td>1626 ± 1; 578 ± 2</td><td>-339 ± 2; 890 ± 2</td><td>-563 ± 2; -391 ± 2</td><td>-287 ± 4; 272 ± 2</td></tr></tbody></table>

correction is provided by SOUL (Pinna et al. 2019), the upgrade of the FLAO system (Esposito et al. 2010). During this observation, the SOUL system was correcting 500 modes with a rate of 1.7 kHz. No coronagraphic masks are available on LUCI, and therefore to avoid any saturation, HR 8799 was observed adopting a sub-windowing of 256 × 256 pixels (corresponding to a field of view (FoV) of 3.84 × 3.84′′). In this way, we were able to set a minimum exposure time of 0.34 s. We observed this target in pupil stabilization mode with two narrow-band filters: FeII ($\lambda_{FeII} = 1.646$ µm; $\Delta\lambda = 0.018$ µm) and H2 ($\lambda_{H2} = 2.124$ µm; $\Delta\lambda = 0.023$ µm). The observations with the two filters were alternated during the whole sequence.

## 3. Data reduction

### 3.1. SPHERE

The reduction of the IRDIS data was carried out entirely with the Data Center (DC; Delorme et al. 2017) which uses the standard SPHERE pipeline SpeCal (Galicher et al. 2018). For the astrometric calibration, we used a true north orientation of – 1.75 deg and a pixel scale value of 12.25 mas as reported in Maire et al. (2016). The reduction algorithm used by the DC for HR 8799 is the T-LOCI (Marois et al. 2014). To provide a homogeneous reduction of all the IRDIS data, we processed the observations presented in Zurlo et al. (2016) again, which were previously reduced with custom routines, which included seven epochs in 2014 (we excluded the sequence of 2014 August 13 for poor quality). The four epochs taken for variability monitoring – once per night (2014 December 4, 2014 December 5, 2014 December 6, 2014 December 8) – and presented in Apai et al. (2016) were reprocessed with the DC. As in Apai et al. (2016), we excluded the sequence of December 8 for poor conditions. In the same way, we used the DC to reduce the data from the variability monitoring of Biller et al. (2021), with observation dates:

2015 July 29, 2015 July 30, 2017 October 07, 2017 October 11, and 2017 October 12 (two epochs were excluded for poor quality), and 2018 August, 17–18. Finally, we re-processed the IRDIS data presented by Wahhaj et al. (2021) and taken on October 31, 2019, and November 1, 2019. In the DC, a “fake negative planets” (e.g., Zurlo et al. 2014, and references therein) algorithm is implemented to calculate the precise position of each planet with respect to the central star. A summary of all the SPHERE astrometric points, both from the SHINE survey and published observations presented in Zurlo et al. (2016), Apai et al. (2016), Biller et al. (2021), and Wahhaj et al. (2021) is listed in Tables 2 (IRDIS) and 3 (IFS). We refer the reader to these latter publications for further information about the observing conditions. The results of the updated astrometry are also listed in Table A.1, together with all the other data from different instruments presented in the literature.

From the IRDIS photometry, we also estimated the mass of each planet using evolutionary models. In particular, we applied evolutionary models from Baraffe et al. (2003, 2015) to the IRDIS filters photometry. The age of the system used in this estimation is 42 Myr and the parallax is 24.46 ± 0.05 mas, as in Gaia EDR3. The error reported is the standard deviation between different epochs taken with the same filter. Results are listed in Table 4.

### 3.2. SOUL+LUCI

The data were taken with a randomized offset of the position of the PSF on the detector. This was used to create a background obtained by calculating the median of all the frames. This background was then subtracted from each science frame. The final dataset was composed of 60 and 61 files for the FeII and H2 observations, respectively. In both cases, each file was composed of 112 frames for a total exposure time of 2284.80 s and 2322.88 s, respectively. During the observations, the derotator
<table><caption>**Table 3.** List of astrometric points from IFS observations.</caption><thead><tr><th rowspan="2">Date</th><th colspan="2">Planet d</th><th colspan="2">Planet e</th></tr><tr><th>ΔRA ; ΔDec (mas)</th><th></th><th>ΔRA ; ΔDec (mas)</th><th></th></tr></thead><tbody><tr><td>2014-07-12</td><td>-400 ± 4; -512 ± 4</td><td></td><td>-389 ± 1; -22 ± 2</td><td></td></tr><tr><td>2014-08-13</td><td>-396 ± 1; -524 ± 1</td><td></td><td>-389 ± 1; -17 ± 2</td><td></td></tr><tr><td>2015-07-03</td><td>-424 ± 4; -509 ± 3</td><td></td><td>-391 ± 1; 33 ± 2</td><td></td></tr><tr><td>2015-09-27</td><td>-420 ± 4; -513 ± 4</td><td></td><td>-392 ± 1; 39 ± 3</td><td></td></tr><tr><td>2016-11-17</td><td>-464 ± 1; -486 ± 2</td><td></td><td>-382 ± 2; 94 ± 6</td><td></td></tr><tr><td>2017-06-14</td><td>-473 ± 2; -476 ± 2</td><td></td><td>-377 ± 1; 115 ± 3</td><td></td></tr><tr><td>2017-10-11</td><td>-480 ± 2; -478 ± 2</td><td></td><td>-371 ± 1; 129 ± 2</td><td></td></tr><tr><td>2017-10-12</td><td>-492 ± 5; -463 ± 6</td><td></td><td>-370 ± 2; 135 ± 3</td><td></td></tr><tr><td>2018-06-18</td><td>-495 ± 1; -460 ± 2</td><td></td><td>-360 ± 2; 162 ± 2</td><td></td></tr><tr><td>2018-08-17</td><td>-509 ± 2; -452 ± 3</td><td></td><td>-361 ± 2; 166 ± 1</td><td></td></tr><tr><td>2018-08-18</td><td>-503 ± 1; -456 ± 2</td><td></td><td>-359 ± 1; 167 ± 2</td><td></td></tr><tr><td>2019-10-31</td><td>-527 ± 3; -432 ± 3</td><td></td><td>-337 ± 2; 210 ± 3</td><td></td></tr><tr><td>2019-11-01</td><td>-528 ± 1; -435 ± 2</td><td></td><td>-337 ± 2; 208 ± 1</td><td></td></tr><tr><td>2021-08-20</td><td>-569 ± 1; -390 ± 2</td><td></td><td>-292 ± 2; 276 ± 2</td><td></td></tr></tbody></table>

was switched off to allow the rotation of the FOV and be able to use the ADI method. Before performing the high-contrast imaging method, we obtained the position of the stellar PSF for each frame using the FIND procedure as implemented for IDL. Exploiting these positions, we then precisely registered the whole dataset, positioning the stellar PSFs at the center of each frame. The ADI was then implemented by exploiting the principal component analysis (PCA; Soummer et al. 2012) technique. For the reduction, we tested a different number of principal components, but found that the best solution was to use ten of them. We were able to recover all the known companions with S/N ranging between 7 and 10.

We obtained the precise position of each planet, introducing fake negative companions and changing its position to minimize the standard deviation in a small region around the planet itself. The results of this procedure are listed in Table 5.

# 4. Astrometric fit of the orbits

## 4.1. Dynamical constraints on the astrometry

Given the literature regarding astrometric fits of the HR 8799 system (e.g., Wang et al. 2018; Goździewski & Migaszewski 2020, and references therein), it is now widely recognized that the present observation time-window does not make it possible to determine long-term stable orbital solutions of HR HR 8799 without invoking particular dynamical constraints. Such constraints may arise from the coplanarity of the planets' orbits and their relatively small eccentricities, which implies a ratio of orbital periods close to 2:1 for successive pairs of planets, that is, 2:1 mean-motion resonances (MMRs). This assumption can be further supported by the likely origin of the system from planetary migration (e.g., Wang et al. 2018); see also Sect. 5 in this paper. Also, the recent analysis of high-contrast images and resulting orbital solutions in Wahhaj et al. (2021) indicate that planets HR 8799e and HR 8799c are most consistent with co-planar and resonant orbits. The dynamical multiple MMR scenario therefore seems to be well justified and documented

in the literature. Recently, Goździewski & Migaszewski (2020) (GM2020 hereafter) reported an exact-resonance configuration (or “periodic orbits model” (PO) hereafter). This model is designed to explain the astrometric measurements through constraints on geometry, planetary masses, and astrometric parallax of the system. Direct determination of planetary masses is crucial for calibrating astrophysical cooling models and for determining the origin and long-term orbital evolution of the entire system, including its debris disk components.

Here, we consider a less stringent condition on the stability of the system, assuming that it may be close to exact resonance but not necessarily fully periodic. This assumption may be natural in the sense that migration may result in a system that is not exactly in a configuration of PO (e.g. Ramos et al. 2017). We want to determine whether such near-resonance models can yield statistically different or perhaps even better best-fitting solutions than those that appear in the PO scenario. We define the following merit functions as in Goździewski & Migaszewski (2018, 2020):

$$ \ln \mathcal{L}(\mathbf{x}) = -\frac{1}{2} \chi^2(\mathbf{x}) - \frac{1}{2} \sum_{i=1}^{N_{\text{obs}}} \left[ \ln \theta_{i,\alpha}^2 + \ln \theta_{i,\delta}^2 \right] - N_{\text{obs}} \ln(2\pi), \\ \chi^2(\mathbf{x}) = \sum_{i=1}^{N_{\text{obs}}} \left[ \frac{[\alpha_i - \alpha(t_i, \mathbf{x})]^2}{\theta_{i,\alpha}^2} + \frac{[\delta_i - \delta(t_i, \mathbf{x})]^2}{\theta_{i,\delta}^2} \right], \quad (1) $$

where $(\alpha_i, \delta_i)$ are the right ascension (RA) and declination (Dec) measurements at time $t_i$; $\alpha(t_i, \mathbf{x})$, $\delta(t_i, \mathbf{x})$ are for their ephemeris (model) values at the same moments as implied by the adopted $\mathbf{x}$ parameters; and $\theta_{i,\alpha}^2$ and $\theta_{i,\delta}^2$ are the nominal measurement uncertainties in RA and Dec scaled in quadrature with the so-called error floor, $\theta_{i,\alpha}^2 = (\sigma_{i,\alpha}^2 + \sigma_{\alpha,\delta}^2)$ and $\theta_{i,\delta}^2 = (\sigma_{i,\delta}^2 + \sigma_{\alpha,\delta}^2)$, for each datum, respectively. Also, if $N_{\text{obs}}$ is the number of observations then $N = 2N_{\text{obs}}$, because RA and Dec are measured in a single detection. The error floor can be introduced for unmodeled measurement uncertainties, such that the resulting best-fitting model should yield the reduced $\chi_v^2 \approx 1$. However, in this work, as we rely on MCMC sampling and the error floor introduces little qualitative variability into the solutions, we omit this parameter, thereby also simplifying the astrometric model.

## 4.2. Mass and parallax priors

We conducted MCMC sampling from the posterior defined through the general merit function in Eq. (1) and astrophysical and geometrical priors. The most important astrophysical priors are the masses of the star and the planets, and the geometrical prior is the parallax.

Regarding the stellar mass, we considered two fixed values of $m_\star = 1.52 M_\odot$ determined by Baines et al. (2012) as $m_\star = 1.516^{+0.038}_{-0.024} M_\odot$ for the age of $\approx 33$ Myr, and $m_\star = 1.47^{+0.12}_{-0.17} M_\odot$, following the most recent dynamical estimate by Sepulveda & Bowler (2022). The same value was used by Wang et al. (2018) in their earlier work. Kervella et al. (2022) determined $m_\star = 1.50 \pm 0.08 M_\odot$, which overlaps with the two estimates. We note here that the stellar mass is fixed in our astrometric N-body model because it also constrains the parallax $\Pi$ of the whole system. These two parameters are strongly correlated through the near-III Keplerian law. Moreover, the parallax tends to be systematically tightly bounded in subsequent Gaia catalogs. The most recent Gaia eDR3 estimate of $\Pi = 24.460 \pm 0.045$ mas is accurate to 0.1%, and it is a meaningful prior for the MCMC sampling. We note that the parallax did not change in the final Gaia DR3 catalog; now it is listed as $\Pi = 24.462 \pm 0.046$ mas.
<table><caption>**Table 4.** Values for the mass of each planet calculated using the evolutionary models in the different IRDIS filters.</caption><thead><tr><th>Filter</th><th>Planet e (M<sub>Jup</sub>)</th><th>Planet d (M<sub>Jup</sub>)</th><th>Planet c (M<sub>Jup</sub>)</th><th>Planet b (M<sub>Jup</sub>)</th></tr></thead><tbody><tr><td>BB_H</td><td>8.9 ± 1.1</td><td>8.4 ± 0.9</td><td>8.5 ± 0.9</td><td>6.4 ± 0.6</td></tr><tr><td>D_H2</td><td>7.0 ± 0.3</td><td>6.6 ± 0.3</td><td>6.7 ± 0.3</td><td>4.6 ± 0.3</td></tr><tr><td>D_H3</td><td>8.5 ± 0.3</td><td>8.1 ± 0.2</td><td>8.3 ± 0.3</td><td>6.8 ± 0.3</td></tr><tr><td>D_K1</td><td>9.6 ± 0.3</td><td>9.3 ± 0.3</td><td>9.3 ± 0.2</td><td>6.7 ± 0.3</td></tr><tr><td>D_K2</td><td>10.5 ± 0.1</td><td>10.5 ± 0.1</td><td>10.4 ± 0.1</td><td>9.0 ± 0.3</td></tr></tbody></table>

Notes. The age of the system is assumed to be 42 Myr. The error is the standard deviation between different observations with the same filter.

<table><caption>**Table 5.** List of astrometric points from the LBT/LUCI observation.</caption><thead><tr><th>Date</th><th>Planet b<br>ΔRA; ΔDec (mas)</th><th>Planet c<br>ΔRA; ΔDec (mas)</th><th>Planet d<br>ΔRA; ΔDec (mas)</th><th>Planet e<br>ΔRA; ΔDec (mas)</th></tr></thead><tbody><tr><td>2020-09-29</td><td>1620 ± 1; 591 ± 3</td><td>-364 ± 2; 883 ± 1</td><td>-551 ± 1; -415 ± 4</td><td>-315 ± 3; 242 ± 5</td></tr></tbody></table>

We also implied Gaussian priors $N_m = [8.7 \pm 1.7, 8.7 \pm 1.0, 8.7 \pm 1.0, 6 \pm 0.7]$ $M_{Jup}$ for the masses of the innermost to the outermost planet, respectively, based on the hot-start cooling models by Baraffe et al. (2003), following discussion in Wang et al. (2018) and confirmed by our estimates collected in Table 4. The mass ranges for BB_H are consistent with the priors in Wang et al. (2018), and D_H3, D_K1 overlap with the determinations from the earlier astrometric model in GM2020. As also demonstrated by GM2020 and confirmed below, there are many discrepancies regarding the mass hierarchy; according to the resonant model, HR 8799d is the most massive one, and the masses of HR 8799e and HR 8799c appear strongly anti-correlated. Masses for D_K2 seem to be too far into the high-end range regarding the dynamical stability, especially that of planet HR 8799b.

It is common in the literature to assume smaller planet masses, that is, of $\approx 7 M_{Jup}$, given the results of dynamical N-body simulations. In particular, Wang et al. (2018) report difficulty in finding dynamically stable solutions for planet masses larger than $\approx 7 M_{Jup}$. Simultaneously, Goździewski & Migaszewski (2018, 2020) found rigorously stable systems locked in the Laplace resonance for higher mass ranges as well that is, of, $\approx 8-9 M_{Jup}$. Given systematic observationally constrained shifts to higher masses, as predicted by the exact Laplace MMR model, we decided to apply the mass priors in the $\approx 9 M_{Jup}$ range. Brandt et al. (2021) estimated the mass of HR 8799e as $9.6^{+1.9}_{-1.8} M_{Jup}$ based on the secular variation of the proper motion of the parent star in Gaia and HIPPARCOS catalogs. Moreover, in a very recent work, Kervella et al. (2022) similarly determined the dynamical mass of $12 \pm 3.5 M_{Jup}$ for HR 8799e. This latter is less precise than in the prior work of Brandt et al. (2021), because Kervella et al. (2022) did not account for the actual orbital geometry. These high mass ranges for HR 8799e are consistent with the larger mass priors. Therefore, in some experiments we also used the mass priors $N_m$ with the innermost mass changed to the value computed in Brandt et al. (2021).

As the parallax prior, we defined the most recent Gaia eDR3/DR3 estimate of $\Pi = 24.46 \pm 0.05$ mas. This value is

significantly shifted from $\Pi \approx 24.22 \pm 0.08$ mas in the Gaia DR2 and $24.8 \pm 0.7$ mas in Gaia DR1 catalog, respectively. The parallax determinations in the Gaia catalogs appear subtly biased, depending on the luminosity and spectral type (Lindegren et al. 2021). However, compared to the uncertainties of the dynamical estimates here, the predicted parallax correction for Gaia eDR3/DR3 of $<0.1$ mas (a few tens of µas) would be insignificant. Indeed, the parallax correction predicted by Lindegren et al. (2021) yields $24.50 \pm 0.05$ mas (Kervella et al. 2022). This apparently very small shift of the order of $1\sigma$ uncertainty may still be meaningful when compared to the dynamical estimates, as is found below.

### 4.3. The exact Laplace resonance revisited

In the first step, we verified whether or not the resonant model in GM2020 based on measurements in Konopacky et al. (2016) fits the recent measurements listed in Tables 2, 3, and 5. We found that this model visually and significantly deviates from the updated data set, especially for planet HR 8799b. Therefore, we refined the PO model using the same parametrization as for the merit function; see Eq. (1) in GM2020. In the formulation presented by these latter authors, the primary parameters of the astrometric model consist of all planet masses $m_i$, $i = 1, \dots, 4$, the osculating period ratio $\kappa$ for the two innermost planets – which selects the particular four-body Laplace resonance chain–, as well as the reference epoch $\tau$, three Euler angles $(I, \omega_{\text{rot}}, \Omega)$ rotating the coplanar resonant configuration to the sky (observer) plane, and the system parallax $\Pi$. This means 11 free parameters. We note that all remaining orbital parameters are constrained through the N-body dynamics confined to the manifold of strictly periodic solutions.

We performed the MCMC sampling from the posterior defined with the merit function in Eq. (1). This merit function is based on a numerical procedure for computing periodic solutions in a co-planar N-planet system, as described in GM2020 and implemented by Cezary Migaszewski (priv. comm.). To conduct the MCMC sampling, we used the affine sampler in emcee
<table><caption>**Table 6.** Osculating, heliocentric elements of the best-fitting solutions at the epoch of 1998.830.</caption><thead><tr><th colspan="8">Model IV<sub>PO</sub>: m<sub>★</sub> = 1.52 M<sub>☉</sub>, Π = 24.2585686 mas (24.26 ± 0.05) mas, χ<sup>2</sup><sub>v</sub> = 4.32, v = 441, dim p = 11</th></tr><tr><th></th><th>m (M<sub>Jup</sub>)</th><th>a (au)</th><th>e</th><th>I (deg)</th><th>Ω (deg)</th><th>ω (deg)</th><th>M (deg)</th></tr></thead><tbody><tr><td rowspan="2">HR 8799e</td><td>8.1 ± 0.7</td><td>16.28 ± 0.03</td><td>0.147 ± 0.001</td><td></td><td></td><td>110.8 ± 0.5</td><td>336.8 ± 0.5</td></tr><tr><td>8.24465258</td><td>16.2817072</td><td>0.147508399</td><td></td><td></td><td>110.765660</td><td>336.793683</td></tr><tr><td rowspan="2">HR 8799d</td><td>9.5 ± 0.1</td><td>26.78 ± 0.08</td><td>0.115 ± 0.001</td><td></td><td></td><td>29.4 ± 0.9</td><td>59.8 ± 0.8</td></tr><tr><td>9.48438158</td><td>26.7888442</td><td>0.01154150</td><td></td><td></td><td>29.4676779</td><td>59.6448288</td></tr><tr><td rowspan="2">HR 8799c</td><td>7.7 ± 0.5</td><td>41.40 ± 0.10</td><td>0.055 ± 0.002</td><td>26.9 ± 0.2</td><td>62.2 ± 0.4</td><td>92.1 ± 0.5</td><td>145.6 ± 0.8</td></tr><tr><td>7.64047381</td><td>41.3945527</td><td>0.05485344</td><td>26.8923597</td><td>62.2416412</td><td>92.1126122</td><td>145.7203489</td></tr><tr><td rowspan="2">HR 8799b</td><td>6.0 ± 0.4</td><td>71.95 ± 0.13</td><td>0.017 ± 0.001</td><td></td><td></td><td>44.4 ± 3.0</td><td>309.2 ± 2.6</td></tr><tr><td>6.01476193</td><td>71.9477733</td><td>0.017597398</td><td></td><td></td><td>44.6384383</td><td>309.291519</td></tr></tbody></table><table><thead><tr><th colspan="8">Model IV<sub>PO</sub>: m<sub>★</sub> = 1.47 M<sub>☉</sub>, Π = 24.5256617 mas (24.53 ± 0.04) mas, χ<sup>2</sup><sub>v</sub> = 4.32, v = 441, dim p = 11</th></tr><tr><th></th><th>m (M<sub>Jup</sub>)</th><th>a (au)</th><th>e</th><th>I (deg)</th><th>Ω (deg)</th><th>ω (deg)</th><th>M (deg)</th></tr></thead><tbody><tr><td rowspan="2">HR 8799e</td><td>7.6 ± 0.9</td><td>16.10 ± 0.03</td><td>0.148 ± 0.001</td><td></td><td></td><td>110.8 ± 0.5</td><td>336.8 ± 0.3</td></tr><tr><td>7.4060918</td><td>16.1038912</td><td>0.147675953</td><td></td><td></td><td>110.816707</td><td>336.828050</td></tr><tr><td rowspan="2">HR 8799d</td><td>9.2 ± 0.1</td><td>26.46 ± 0.09</td><td>0.115 ± 0.002</td><td></td><td></td><td>29.1 ± 1.0</td><td>60.0 ± 1.0</td></tr><tr><td>9.18999484</td><td>26.45172667</td><td>0.114600334</td><td></td><td></td><td>29.0594738</td><td>60.1752106</td></tr><tr><td rowspan="2">HR 8799c</td><td>7.7 ± 0.7</td><td>40.97 ± 0.10</td><td>0.054 ± 0.002</td><td>26.9 ± 0.2</td><td>62.2 ± 0.4</td><td>92.3 ± 0.6</td><td>145.5 ± 0.7</td></tr><tr><td>7.79611845</td><td>40.9861801</td><td>0.05380902</td><td>26.8734271</td><td>62.1852189</td><td>92.3879462</td><td>145.5347082</td></tr><tr><td rowspan="2">HR 8799b</td><td>5.8 ± 0.4</td><td>71.14 ± 0.15</td><td>0.017 ± 0.002</td><td></td><td></td><td>43.50 ± 3.0</td><td>310.35 ± 3.0</td></tr><tr><td>5.80235941</td><td>71.1427851</td><td>0.01667073</td><td></td><td></td><td>42.9704380</td><td>310.895065</td></tr></tbody></table><table><thead><tr><th colspan="8">Model IV<sub>Δn</sub>: m<sub>★</sub> = 1.47 M<sub>☉</sub>, Π = 24.5256617 mas (24.49 ± 0.07) mas, χ<sup>2</sup><sub>v</sub> = 3.9, v = 429, dim p = 23</th></tr><tr><th></th><th>m (M<sub>Jup</sub>)</th><th>a (au)</th><th>e</th><th>I (deg)</th><th>Ω (deg)</th><th>ω (deg)</th><th>M (deg)</th></tr></thead><tbody><tr><td rowspan="2">HR 8799e</td><td>7.5 ± 0.6</td><td>16.00 ± 0.10</td><td>0.150 ± 0.004</td><td></td><td></td><td>108.8 ± 1.3</td><td>336.00 ± 1.5</td></tr><tr><td>7.40609179</td><td>16.1038912</td><td>0.147675954</td><td></td><td></td><td>110.816707</td><td>336.828050</td></tr><tr><td rowspan="2">HR 8799d</td><td>9.3 ± 0.5</td><td>26.55 ± 0.12</td><td>0.114 ± 0.004</td><td></td><td></td><td>28.66 ± 1.2</td><td>60.0 ± 1.0</td></tr><tr><td>9.18999484</td><td>26.45172667</td><td>0.11460033</td><td></td><td></td><td>29.05947383</td><td>60.1752106</td></tr><tr><td rowspan="2">HR 8799c</td><td>7.8 ± 0.6</td><td>41.00 ± 0.15</td><td>0.050 ± 0.003</td><td>26.5 ± 0.5</td><td>63.1 ± 0.5</td><td>91.8 ± 1.0</td><td>145.0 ± 1.0</td></tr><tr><td>7.79611845</td><td>40.9861801</td><td>0.05380902</td><td>26.8734271</td><td>62.1852189</td><td>92.3879462</td><td>145.5347082</td></tr><tr><td rowspan="2">HR 8799b</td><td>5.8 ± 0.8</td><td>71.3 ± 0.2</td><td>0.017 ± 0.002</td><td></td><td></td><td>42.3 ± 1.6</td><td>310.8 ± 2.0</td></tr><tr><td>5.80235941</td><td>71.1427851</td><td>0.01667073</td><td></td><td></td><td>42.9704380</td><td>310.895065</td></tr></tbody></table>

**Notes.** The stellar mass is $m_{\star} = 1.52 M_{\odot}$ or $m_{\star} = 1.47 M_{\odot}$. Uncertainties for parameters in models IV<sub>PO</sub> and IV<sub>Δn</sub> are determined as the 16th and 86th percentiles of the samples. Parameter values with 8–9 significant digits are provided in order to reproduce the results of particular integrations. Parameter $v = N_{\text{obs}} - \text{dim } p$ is for the degrees of freedom, where $N_{\text{obs}}$ denotes the number of $(\alpha, \delta)$ measurements, and $\text{dim } p$ is the number of free parameters.

package (Foreman-Mackey et al. 2013). We initiated 1536 walkers in a small hyper-ball in the parameter space around the initial condition in Table 6 determined with a search with the Powell non-gradient minimization algorithm. Since the autocorrelation time appears relatively short, of ~100 iterations only, we ended up with 1000 iterations that make it possible to derive the posterior distribution, given the large number of walkers.

The derived posterior is illustrated in 2-dim projections and 1-dim histograms of the MCMC samples for the primary parameters (Fig. 1). This experiment reveals significant correlations between different parameters already reported in GM2020 that still cannot be reduced with the new data. Besides a strong $m_e - m_c$ anti-correlation, there is also $m_e - \kappa$ correlation and $m_d - \kappa$ anti-correlation. The remaining two masses are free from significant mutual correlations, however, the mass of HR 8799b

appears strongly correlated with geometric parameters of the system, i.e., the Euler angles and the parallax. These parameters show mutually weaker but still significant correlations.

Parameters of the best-fitting solutions with their formal uncertainties and particular solutions from the posterior that may be useful for future numerical studies are collected in Table 6, referred to as Model IV<sub>PO</sub> for the two stellar masses of $m_{\star} = 1.52 M_{\odot}$ and $m_{\star} = 1.47 M_{\odot}$, respectively.

For a stellar mass of $m_{\star} = 1.47 M_{\odot}$, the HR 8799d mass seems to be constrained to ≈9.2 $M_{\text{Jup}}$ within a relatively very small uncertainty of ≈0.1 $M_{\text{Jup}}$. This mass estimate is slightly larger than that found by GM2020. The mass of the outer-most planet HR 8799b may be determined as ≈5.8 $M_{\text{Jup}}$ which is similar to ≈5.5 $M_{\text{Jup}}$ found by these latter authors, yet with also smaller uncertainty of ≈0.3 $M_{\text{Jup}}$. Moreover, for the larger stellar
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** MCMC posterior for the best-fitting, strictly resonant model derived for $m_* = 1.47 M_{\text{Sun}}$. Parameters included in the diagram are for the planet masses and orbital period ratio $\kappa = P_d/P_e = \kappa$, and parallax $\Pi$, inclination $I$, and the nodal angle $\Omega$ (the orbital scale $\rho$, PO epoch shift $t_0$ and the common rotation angle $\omega_{\text{rot}}$ in the orbital planet are skipped). The crossed lines mark the best-fitting PO solution in terms of the smallest $\chi^2$ that was used as a starting point to initiate the emcee walkers; see Tables 6 and 7 (Model 1) for parameters of this solution. Uncertainties for the parameters are determined as the 16th and 86th percentile of the samples around the median values.</figcaption></figure>

mass of $1.52 M_{\odot}$, we derived slightly larger masses of HR 8799d $\simeq 9.6 \pm 0.1 M_{\text{Jup}}$ and HR 8799b $\simeq 6.0 \pm 0.4 M_{\text{Jup}}$.

The astrometric model for the updated data window yields $\Pi \simeq 24.26 \pm 0.05$ mas, which is consistent with the earlier estimate by GM2020. The parallax estimate of $\Pi \simeq 24.53 \pm 0.04$ mas derived for $1.47 M_{\odot}$ overlaps to $1\sigma$ with $\Pi \simeq 24.460 \pm 0.045$ mas in Gaia eDR3/DR3. Moreover, as noted above, when applying the parallax correction of the Lindegren et al. (2021) eDR3 catalog, Kervella et al. (2022) obtained $\Pi = 24.50 \pm 0.05$ mas which is even closer to the value derived from the PO model.

This relation confirms the predicted strong stellar mass–parallax correlation and indirectly favors $m_* = 1.47 M_{\odot}$. The dynamical parallax determinations appear to be very accurate and meaningful, overlapping closely with the independent geometric Gaia measurements.

Regarding the HR 8799e–HR 8799d mass anti-correlation, we conducted additional MCMC sampling for the primary parameters in the PO model for fake data series. We prepared two or three synthetic observations per year, extending the real observations window by $\simeq 20$ yr around the best-fitting periodic
<figure><img src="image_3.png" /><figcaption>Fig. 2. Temporal evolution of the osculating orbital elements in the astrocentric (blue curves) and Jacobian (yellow curves) frames, respectively, for the same N-body initial condition implying a stable, near 8:4:2:1 MMR system. We note that the elements overlap for the first planet according to the construction of the Jacobian reference frame.</figcaption></figure>

model with deviations and uncertainties of ~1 mas. It appears that all correlations except those between the masses and $\kappa$ would be almost eliminated. A probable cause of that effect may be the almost perfectly aligned gravitational tugs of these planets in the present particular time window of the observations. If the PO model is correct, then likely only highly precise data similar to the most accurate GRAVITY measurement in GRAVITY Collaboration (2019) could break this degeneracy. This was indicated by GM2020. We note that the GRAVITY measurement is depicted with a star in all figures illustrating orbital solutions.

### 4.4. Near 8:4:2:1 MMR model

To characterize the system near the 2:1 MMR chain, we rely on the notion of the proper mean motions as the fundamental frequencies $f_i$ in the framework of conservative N-body dynamics. We resolve these frequencies with the refined Fourier frequency analysis (e.g., Laskar 1993; Nesvorny & Ferraz-Mello 1997, Numerical Analysis of Fundamental Frequencies, NAFF). To perform the NAFF, we consider the time series

$$S_{\lambda,i} = \{a_i(t_j) \exp i\lambda_i(t_j)\}, \quad S_{M,i} = \{a_i(t_j) \exp iM_i(t_j)\},$$

where $j = 0, 1, \dots, t_j = j\Delta t$, and $a_i(t), \lambda_i(t)$ and $M_i$ are the canonical osculating semi-major axis, the mean longitude, and the mean anomaly for planet $i$, $i = 1, 2, 3, 4$, respectively, and $\Delta t$ is the sampling interval ($i$ is the imaginary unit). These canonical elements inferred in the Jacobi or Poincaré reference frame make it possible to account for mutual interactions between the planets to the first order in the masses.

We note that the canonical orbital elements evolve differently from the common astrocentric elements. This is illustrated in Fig. 2. Subsequent panels are for the semi-major axis and eccentricity of HR 8799d,c,b, respectively, computed for the interval of a few tens of outermost orbits, and marked with different colors. It turns out that these elements span much wider ranges in the common astrocentric Keplerian frame. This effect is particularly noticeable for the two outermost planets. A small change in the N-body initial condition may result in substantial displacement of the orbit in the Keplerian $(a, e)$-plane. Therefore, the geometric elements should be understood as a formal representation of the Cartesian N-body initial conditions, and the canonical elements are more suitable for the qualitative characterization of the orbits.

Considering a particular four-body MMR chain that can explain astrometric observations of the HR 8799 system, the

zeroth-order mean-motion resonance implies one of the possible linear combinations of the fundamental frequencies (the proper mean-motions),

$$\Delta f(\mathbf{k}) \equiv \sum_{i=1}^{N} k_i f_i, \qquad (2)$$

where $\mathbf{k} \equiv [k_1, \dots, k_N]$ is a vector of non-zero integers, and $\sum_{i=1}^N |k_i| \neq 0$, which yields small $|\Delta f|$. Simultaneously, we may determine the critical angle corresponding to $\eta = \arg \min |\Delta f(\mathbf{n})|$:

$$\theta_{\eta}(t) \equiv \sum_{i=1}^{N} \eta_i \lambda_i, \qquad (3)$$

if $\Delta f(\mathbf{k})$ accounts only for the proper mean motions according to the d'Alembert rule. The resonant configuration may therefore be called the generalized Laplace resonance (Papaloizou 2015). This type of MMR, which possibly drives the orbital evolution of HR 8799, is characterized with the proper mean motions fulfilling the following relation (e.g., GM2020):

$$\Delta f(\mathbf{k}) = n_1 - 2n_2 + n_3 - 2n_4. \qquad (4)$$

Similarly to the classic Laplace resonance in the Galilean moons system, we consider the generalized MMR as a chain of two-body MMRs,

$$2n_{i+1} - n_i + g_j \approx 2n_{i+1} - n_i \approx 0, \quad i = 1, 2, 3, \quad j = i, i+1,$$

for subsequent pairs of planets. Here, $g_j$ are for the proper frequencies (mode) associated with the pericenter rotation. In the exact resonance, which we consider as the PO in the reference frame rotating with a selected planet (Hadjidemetriou 1977), the apsidal lines of all planets rotate synchronously, with the same angular velocity relative to the inertial frame. The condition in Eq. (4) for the exact resonance (periodic motion in the rotating frame) may be expressed through the proper mean motions of the planets (e.g., Delisle 2017):

$$\Delta n \equiv \sum_{i=1}^{3} \left( \frac{n_i}{n_{i+1}} - 2 \right) = 0, \qquad (5)$$

given the proper mean-motions resolved from the time series $S_{M,i}$. Numerically, we find that the exact resonance yields $\Delta n \approx 10^{-13}$ rad d$^{-1}$ and $\Delta f \approx 10^{-15}$ rad d$^{-1}$ for 3 × 8192 samples of $\Delta t = 2048$ days.
We used the resonance constraint in Eq. (5) as one of the priors for the MCMC sampling of the posterior defined through the maximum likelihood function in Eq. (1). This indirect frequency prior is Gaussian with a small variance $\sigma_{\Delta n}$ to be determined later.

We determined the frequency variance $\sigma_{\Delta n}$ based on extensive numerical experiments that relied on calibrating the $\Delta n$ range with the Lagrange (geometric) stability interval, which should not be shorter than the stellar age of 40–50 Myr, as predicted in the most recent work (Brandt et al. 2021). These authors also rule out ages for HR 8799 of greater than ≈300 Myr. The proper choice of the $\Delta n$ prior is crucial to bound (resolve) a long-term stable system in a very short integration time that spans merely a few hundred of the outermost orbits, ≈2 × 10⁵ yr, which is also the instability timescale outside the Laplace resonance. We find, as described below in Sect. 4.6, that the astrometric HR 8799 solutions are long-term stable if $|\Delta n|$ is roughly less than 10⁻⁵–10⁻⁶ rad d⁻¹, as described above.

### 4.5. Astrometric fits to near-resonance configurations

To sample the posterior defined with Eq. (1), we adopted priors described in Sect. 4.2. We did not imply any other particular limits on the anticipated near-resonant *N*-body solutions, besides uniform priors for orbital parameters in sufficiently wide ranges. In this sense, the assumed prior set is minimal. Moreover, the *N*-body model is parameterized with ($m_i, a_i, x_i \equiv e_i \cos \bar{\omega}_i, y_i \equiv e_i \sin \bar{\omega}_i, M_i$), $i = 1, \dots, 4$, that is, mass, semi-major axis, Poincaré elements composed of eccentricity and argument of pericenter, and the mean anomaly at the osculating epoch for each orbit, as well as two angles ($I, \Omega$) determining the orbital plane of the system, and the astrometric parallax $\Pi$. This means 23 free parameters.

It should be noted that the implicitly defined frequency prior (Eq. (5)) makes it difficult to derive the posterior distribution. In many experiments performed with a different number of MCMC walkers of up to 256 and up to 256 000 iterations, we could obtain an acceptance ratio as low as 0.1. Therefore, we understand the MCMC sampling with the $\Delta n$ prior as dynamically constrained optimization (rejection sampling) which helps to determine parameter ranges for the best-fitting, stable solutions. Also, due to implicit frequency prior definition, which has to be calibrated consistently with the expected dynamical stability of the system, the method has a heuristic character. Nevertheless, we may gain much insight into the parameter ranges of stable near-resonance models consistent with the astrometric observations.

We performed many MCMC experiments varying $\sigma_{\Delta n}$ within the range of 10⁻⁷–10⁻⁹ rad d⁻¹ and parameters of the *emcee* samplers. The example posterior for $\sigma_{\Delta n} = 10^{-8}$ rad d⁻¹ is illustrated in Fig. 3 for parameters selected consistently with the posterior for PO shown in Fig. 1. During the sampling, we recorded all elements of the tested models, and the semi-amplitude of the critical angles of the Laplace resonance for the integration time spanning ≈200 outermost periods. Fig. 3 shows a collection of samples with $|\Delta n| < 10^{-5}$ rad d⁻¹, and the semi-amplitude of Laplace resonance argument $|\Delta\theta| < 60^\circ$. This choice is explained in the following Sect. 4.6, in which we analyze the orbital evolution of particular selected solutions.

Overall, when comparing the posterior distributions in Fig. 1 derived for the PO model, and in Fig. 4 for the quasi-periodic model, we may notice similar ranges for the parameters. A basic distinction relies on the different character of the solutions: while the PO posterior covers strictly stable, resonant solutions, the $\Delta n$ posterior also involves weakly chaotic solutions around the

Laplace resonance center, which is selected as the starting PO model and that may explain missing mass correlations found for the PO model. While the fit quality in terms of $\ln L$ or $\chi^2$ may be somewhat improved with the quasi-periodic model, the RMS of the best-fitting models remains almost the same in the two approaches. It may be concluded that the quasi-periodic model does not lead to any qualitative change in the dynamical status of the system, apart from making it possible to systematically explore long-term stable solutions, which are not necessarily regular. We discuss this further in Sect. 4.6.

Figure 4 illustrates the best-fitting orbital solutions (red curves) over-plotted for models within $1\sigma$ (grey curves) randomly selected from the MCMC samples such that their RMS varies between 7 and 9 mas. The best-fitting models with RMS ≈7.6 mas are marked as red curves. All of those solutions are depicted for roughly one osculating period for each planet. Clearly, the orbits do not close for all companions. This effect is best visible for the outermost (top-middle panel) and HR 8799d (bottom-right panel) planet, respectively. In the latter case, the range of the orbit splitting can be compared to a relatively large uncertainty on the HST measurements; for more accurate measurements, the splitting would be even more significant. This effect illustrates strong, short-term mutual interactions. Here, we follow Goździewski & Migaszewski (2018), who indicated that Keplerian orbits are already inadequate for representing the proper astrometric solutions despite limited observational orbital arcs.

In the top-right panel Fig. 4, for HR 8799c, blue circles representing measurements in Konopacky et al. (2016) and red circles and yellow hexagons for SPHERE measurements are slightly shifted with respect to each other, and relative to the orbital model arcs. This is further visible in Figs. 5 and 6, which illustrate residuals of the best-fitting model in Table 6 in the ΔRA–ΔDec-plane as well as individual ΔRA(t) and ΔDec(t) panels. While the SPHERE data are roughly uniformly clustered around the origin (0,0), the blue circles exhibit systematic patterns, especially for planets HR 8799b and HR 8799c, besides a much larger spread. This effect is curious once we recall the two data sets consisting of uniformly reduced, essentially homogeneous measurements in Konopacky et al. (2016) and this work, respectively. It is likely that the instruments and/or reduction pipelines are not fully compatible. The effect is quite subtle but still noticeable.

Another trend of residuals is revealed in Fig. 6 for planet HR 8799e. The SPHERE data seem to be clearly arranged along a curve in the RA coordinate, revealing a systematic deviation in time with the ephemeris. This effect is especially clear for IRDIS data and may be noted in the lower right ΔRA–ΔDec panel and the ΔRA(t) panel, as the yellow hexagons spread along the ΔRA = 0 axis. This may indicate a systematic effect in IRDIS data reduction because it does not seem to appear for IFS data. However, the time–ΔRA(t) panel is suggestive of this trend for both sets of measurements. If it has no instrumental origin, then it may indicate an unmodeled component of the astrometric model, such as the presence of a yet unknown object (e.g., a distant moon) or a different curvature of the orbit due to different parameters and geometry (e.g., noncoplanarity). We note here the perfect position of the GRAVITY measurement, but this one point is insufficient to dismiss the trend.

### 4.6. Frequency priors versus the dynamical stability

Regarding the dynamical character of solutions derived in Sect. 4.5, we analyzed representative examples of the best-fitting models marked in Fig. 4. We selected MCMC samples in Fig. 3
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Posterior samples for stable solutions characterized with $|\Delta n| < 10^{-5}$ rad d$^{-1}$, and the semi-amplitude of the libration of the Laplace resonance argument $|\Delta\theta| < 60^{\circ}$. The starting solution that initiated nearby MCMC walkers in the emcee samplers is the periodic configuration in Tables 6 and 7 (Model 1). Uncertainties for the parameters are determined as the 16th and 86th percentiles of the samples around the median values.</figcaption></figure>

for $m_{\star} = 1.47 M_{\odot}$ and from other MCMC sampling for $m_{\star} = 1.52 M_{\odot}$. These initial conditions were integrated with the IAS15 integrator from the REBOUND package (Rein & Spiegel 2015) for 1 Gyr or up to disruption of the system when a collision or ejection of a planet occurs. The results are shown in a sequence of plots in Figs. 7 and 8.

Systems illustrated in Fig. 7 are for long-term stable solutions that survived for at least 1 Gyr. The left column is for exactly resonant Model 1 and its elements displays Table 7 and also Table 6. This solution is stable forever “by design” in the framework of the *N* body dynamics. All but one critical angle of the two-body MMRs between subsequent pairs of planets librate around centers in the [-90°, 90°]-range. The exceptional critical

argument is for the outermost pair of planets and librates around the center at ≈180°. As it comes to the critical argument of the Laplace resonance, it librates with a very small amplitude of just ≈5° around the resonance center at ≈15°. The periodic, perfectly regular character of this solution manifests as time evolution of the eccentricities.

The middle column (Model 2) in Fig. 7 is for a model slightly displaced from the exact resonance, as measured by $|\Delta n| \approx 10^{-7}$ rad d$^{-1}$, yet it yields slightly better fit quality. The critical arguments evolve in the same manner, as for the exact resonance, but their amplitudes and eccentricities become noticeably larger. This solution also yields a fully stable configuration of the planets. It is regular in the sense of the Lyapunov exponent, as it
<figure><img src="image_5.png" /></figure>

**Fig. 4.** Best-fitting solutions to the four-planet model in Table 6 for $m_* = 1.47 M_⊙$, illustrated on the sky plane as randomly selected MCMC samples from the $\Delta n$ posterior. The y-axis corresponds to the north ($N$), and the x-axis corresponds to the east ($E$) direction, respectively (we note that the numerical values of $\Delta\alpha$ are opposite in sign with respect to the formal left-hand direction of the RA $\alpha$). Red-filled circles are for the new or re-reduced IRDIS measurements in this paper, yellow hexagons are for the IFS measurements, green hexagons are for the LUCI points, light-blue filled circles are for measurements in Konopacky et al. (2016), and dark-blue circles are for data in other references collected in Tables 2, 3, and 5 and summarised in Table A.1. Diamonds are for the GPI data, and a star symbol is for the most accurate GRAVITY point. Red curves mark stable solutions with the lowest RMS $\approx 7.6-8$ mas, randomly selected from the MCMC samples, and darker grey curves are for other orbital arcs derived for stable models up to RMS $\approx 9-10$ mas. All model orbits have been derived for all measurements available to date, and cover roughly one osculating orbital period for every planet, respectively. The osculating epoch is 1998.830. We also labeled the observational epochs encompassing the orbital arcs with data. The *upper left panel* is for the global view of the system on the sky plane, and subsequent panels are for its close-ups and interesting regions.

can be seen in Fig. 9 showing dynamical maps in terms of the MEGNO fast indicator (Cincotta et al. 2003). The mean exponential growth factor of nearby orbits (MEGNO; also known as $\langle Y \rangle$) is a numerical technique designed to efficiently characterize the stability of the $N$-body solutions in terms of the maximal Lyapunov exponent (MLE). We implemented MEGNO for the planetary problem in Goździewski et al. (2001) and in our CPU-parallelized µFARM numerical package. It is also clear that this solution lies close to the edge of the resonance island (dark blue color).

Finally, the right column (Model 3) in Fig. 7 is for a marginally stable solution, in the sense of nonzero MLE, with $|\Delta n| \approx 10^{-5} \text{ rad d}^{-1}$. Although it survived for the 1 Gyr integration interval, one of the critical angles in the outermost 2:1 MMR rotates. This solution is peculiar in the sense that, despite rotations of one of the 2:1 MMR critical arguments, the semi-amplitudes of other critical angles are similar to those for quasi-periodic Model 2.

Model 3 in Fig. 7 is also remarkable when we compare it with a sequence of solutions in Fig. 8 illustrating chaotic, yet still long-term stable models characterised by $|\Delta n| \approx 10^{-6} \text{ rad d}^{-1}$. All these models self-destruct between 800 Myr (Model 4, the

left column in Fig. 8), 600 Myr (Model 5, the middle column), and just 200 Myr (Model 6, the right column). Despite this, for the initial few hundred million years, extending safely beyond most approximations of the stellar age, which range between 30–160 Myr (Baines et al. 2012; Sepulveda & Bowler 2022), the systems are bounded and locked in the resonance. Simultaneously, Models 3 and 4 yield masses of the inner planet in the $\approx 10 M_{\text{Jup}}$ range, and Model 4 is especially interesting given its low $\chi^2$ compared to the initial starting PO configuration. These models yield a declining mass hierarchy that resembles that of the outer Solar System, and a mass of HR 8799e is consistent with the dynamical estimate in Brandt et al. (2021).

These examples are to justify that the system may be dynamically long-term stable in the planet mass range of $\approx 10 M_{\text{Jup}}$, even if detuned from the exact resonance and mildly chaotic. In all unstable cases, one of the critical arguments of the 2:1 MMR of the two outermost planets progressively increases its libration amplitude and eventually begins to rotate. In this sense, the outermost pair HR 8799b–c is the weakest link in the resonance chain, provoking instability of the whole system displaced from the resonance. The time for the onset of instability can be relatively
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Residuals to the selected best-fitting, near-resonant four-planet Model 2 in Table 7 for planets HR 8799b and HR 8799c; stellar mass is $m_{\star} = 1.47 M_{\oplus}$. The y-axis corresponds to the north (N), and the x-axis corresponds to the east (E) direction, respectively (we note that the numerical values of ΔRA are sign-opposite to regarding the formal left-hand direction of the RA axis). Red filled circles, and yellow and green hexagons are for the IRDIS, IFS, and LUCI measurements reported here, respectively, and dark-blue and light-blue (light-grey) filled circles are for measurements in previous papers and in Konopacky et al. (2016), respectively. Yellow diamonds are for GPI, the grey diamonds are for the early HST data.</figcaption></figure>

very long, as shown by the evolution of Models 4 and 5. In any case, the libration amplitudes of the critical angles seem to be a weak indicator of instability, because there is no clear relationship between these amplitudes and the time of instability. Also, Model 5 illustrates the difficulty in predicting the system behavior based on the variation of critical angles, especially if the system stability is tested for a limited period of time. The apparently regular, quasi-periodic, and bounded evolution of the critical angles for ≈200 Myr does not prevent the system from eventually becoming unstable after ≈400–500 Myr. A similar effect may be observed for Model 6, although in this case, the critical argument of the Laplace resonance varies irregularly at the beginning of the integration.

It is worth noting that all of the example solutions yield astrometric fits that differ little in quality in terms of $\chi^2$ and RMS. The models are difficult to distinguish statistically and visually from each other. This may mean that for a coplanar, near-resonance, or resonance model we can hardly differentiate between perfectly regular and chaotic evolution of the system, as long as the instability time is sufficiently long relative to the age of the star. However, these two types of configurations occur near the exact Laplace resonance.

## 4.7. Resonant structure of the inner debris disk

Based on the updated orbital solutions in Table 6, we revised and extended simulations of the dynamical structure of the inner debris disk in Goździewski & Migaszewski (2018). These experiments rely on the concept of the so-called 〈*Y*〉-model. We assume that the planets form a system of primaries in safely stable orbits robust to small perturbations. We then inject bodies with masses significantly smaller than those of the primaries and on orbits with different semi-major axes and eccentricities spanning the interesting region, and randomly selected orbital angles. Next, we integrate the synthetic configuration and determine its stability with the MEGNO (〈*Y*〉) fast indicator. Calibration experiments spanning orbital evolution of debris disks in HR 8799 for up to 70 Myr are described in detail in Goździewski & Migaszewski (2018). A comparison of the results of direct numerical integrations with the outcomes of the 〈*Y*〉-model confirms that these two approaches are consistent with one another, yet the 〈*Y*〉-based method is much more CPU efficient and therefore makes it possible to obtain a clear representation of the structure of stable solutions. This algorithm is especially effective for strongly interacting systems. Also, these simulations
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Residuals to the best-fitting, near-resonant four-planet Model 2 in Table 7 for $m_* = 1.47 M_\odot$, a continuation of Fig. 8 for planets HR 8799d and HR 8799e. See the caption of Fig. 5 for labels. The yellow star marks the GRAVITY measurement for HR 8799e.</figcaption></figure>

revealed that the close proximity of the four primaries to the exact Laplace resonance makes the system stability robust even to apparently significant perturbations caused by probe masses as large as ≈2 $M_{\text{Jup}}$.

We considered two types of probe objects: Ceres-like asteroids with a mass of $10^{-6} M_{\text{Jup}}$ and Jupiter-like planets with a mass of $1 M_{\text{Jup}}$. To calculate the $\langle Y \rangle$ values for the synthetic systems, we integrated the N-body equations of motion and their variational equations with the GraggBulirsch–Stoer integrator (Hairer et al. 2000) for 3 Myr, which covers $6 \times 10^4$ orbital periods of HR 8799e and $7 \times 10^3$ orbital periods of HR 8799b, respectively. That integration time is consistent with the typical $10^4$ outermost orbits required to achieve $\langle Y \rangle$ convergence. We also note that the most significant interactions are exerted by the inner planets. We chose the integration time that is optimal from the CPU overhead point of view.

The results for the less massive Ceres-like asteroids are illustrated in Fig. 10 and Fig. 11. In this experiment, we sampled the semi-major axis $a_0 \in [4, 18]$ au and eccentricity $e_0 \in [0, 0.8]$ of these objects. We collected $3 \times 10^5$ N-body initial conditions with $|\langle Y \rangle - 2| < 0.05$ that represent the structure of long-term-stable orbits in the inner debris disk. Subsequent panels in Fig. 10 illustrate snapshots of the disk at other epochs following the first observation ($t_0 = 1998.830$) as seen on the sky plane. To obtain the snapshots, we numerically integrated the whole set of initial conditions for planets and the Ceres-like objects defined at

the epoch 1998.830 up to the given final epoch *t*. For instance, *t* = 2009.575 (middle-right panel) is for the first detection of the innermost planet in Marois et al. (2010).

We mark the model orbits of the two inner planets with red curves, and the astrometric measurements are over-plotted on them. Instant positions of the planets are marked with shaded large filled circles (for the first HST epoch), and large filled circles for positions of the planets at the particular epoch labeled in the lower right corner of each panel. The time range of about 30 yr is more than half the orbital period of HR 8799e and roughly one orbital period of an asteroid involved in the 2:1e MMR with this planet.

Positions of the test particles are marked with different colors depending on their dynamical status: yellow and orange dots are for objects involved in 1:1e MMR with the innermost planet HR 8799e; green dots are for the 3:2e MMR, blue dots are for the 2:1e MMR, and gray dots are for the other stable asteroids. As can be seen in Fig. 10, the outer edge of the disk is highly nonsymmetric and quickly evolves in time. Its temporal structure strongly changes even for a relatively very short interval of 25 years spanned by the astrometric observations of the system.

The structure of the disk has a clear resonant structure that was also noted by Goździewski & Migaszewski (2018). This is demonstrated in Fig. 11, which shows the distribution of canonical osculating elements inferred in the Jacobi reference frame, in the $(a_0, e_0)$-plane (top panel), and on the plane of relative
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Orbital evolution of selected long-term-stable solutions with similar astrometric fit quality, but exhibiting qualitatively different stability signatures. The left column is for the exact Laplace resonance, the middle column is for a model displaced from the PO but rigorously stable, and the right column is for a >1 Gyr stable solution characterized by rotation of one critical angle of the 2:1 MMR in the outermost pair. The rows from the top to bottom are for all critical angles of the 2:1 MMR, the critical angle of the generalized Laplace MMR, and eccentricities for subsequent pairs of planets. See Table 7 for the initial conditions labeled from 1 to 3 at the top-left corner of each column.</figcaption></figure>

orbital phases, ($a_0$, $\lambda_e$ – $\lambda_0$) (bottom panel). In both graphs, we mark the test particles with the same color scheme as in the previous plot, and the low-order resonances with HR 8799e are labeled. We also plot the geometrical collision curve (in gray) with planet HR 8799e constrained through the apocenter distance of the inner orbit equal to the pericenter distance of the planet, $a_0(1 + e_0) = a_e(1 - e_e)$. Also, the curve depicted in red is an image of the collision curve shifted by $\Delta a_0 = 4$ au towards the star and clearly marks a boundary of stable particles.

The bottom graph shows objects involved in low-order MMRs with HR 8799e that are apparently widely spread on the sky plane, but they appear in narrow and well-localized islands on the planes of orbital elements; this is particularly clear on

the ($a_0$, $\lambda_e$ – $\lambda_0$)-plane. Both graphs mark two families of asteroids in 1:1e MMR – depicted in gold colour ($e_0 < 0.3$) and in orange ($e_0 > 0.3$). The low-eccentricity objects resemble classic Trojan asteroids in the Solar System, relative by ±60° to the planet. The second unusual family of highly eccentric particles in 1:1e MMR can be found on the sky plane far beyond the orbit of HR 8799e. At some epochs (e.g., 2011.83), these appear closer to the orbit of HR 8799d, which may be counter-intuitive. Given the 1:1e MMR dynamical classification, they would be expected to share the same orbit as the planet.

This simulation reveals the dynamical structure composed of long-term-stable but mutually noninteracting asteroids in the system, on the same orbit as the innermost planet, whose orbital
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Orbital evolution of selected long-term-stable chaotic solutions with similar astrometric fit quality. The rows from the top to bottom are for all critical angles of the 2:1 MMR, the critical angle of the generalized Laplace MMR, and eccentricities for subsequent pairs of planets. See Table 7 for the initial conditions labeled from 5 to 8 at the top of each column.</figcaption></figure>

evolution is governed by the gravitational interactions with the massive Jovian companions. Such objects should not be significantly affected by nongravitational forces, such as Poynting-Robertson drag or the Yarkovsky effect. However, such forces may be crucial for the investigation and modeling of the actual emission profile of dust produced by collisional dynamics. The results indicate that the asteroid breakup events may be frequent and violent. Given the large eccentricities and the instant mixing of different resonant fractions of these objects in the wide inner region between 4 au (and below) and 18 au, including wide Lagrange clumps of planet HR 8799e, their orbital velocity dispersion may be significant. The interplay of all these factors – and especially the clearing rate of the dust due to radiation of the young and bright host star – determines the emission profile. The present observational evidence is limited, and the position of

the warm disk is estimated down to 610 au (e.g., Su et al. 2009; Matthews et al. 2014; Chen et al. 2014), which coincides with the 2:1e MMR (Fig. 11); see also Sect 5. Given the topic is complex and somewhat out of the scope of the present work, we postpone in-depth analysis of the dynamical structure of the warm dust disk to a future paper.

## 4.8. Putative fifth Jupiter-like planet

We performed a very similar experiment for test particles of mass $1 M_{Jup}$, a hypothetical planet below the current detection level that may exist in the inner part of the HR 8799 system. The results are illustrated in Fig. 12 in a similar way to in Figs. 10–11. We collected $10^5$ stable solutions. Again, the positions of the hypothetical planet are marked at characteristic epochs; for
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Dynamical maps in the Keplerian, osculating astrocentric semi-major-axis–eccentricity plane for the quasi-periodic, near-resonant configuration described by Model 2. All masses and elements but the map coordinates are fixed at their best-fitting values in Table 7. The osculating epoch is 1998.830. Stable solutions are determined with $|\langle Y \rangle - 2| \approx 0$ and marked with blue color.</figcaption></figure>

example, near the first epoch of the HST detection, up to the last epoch of measurements in this work, and one epoch ahead of that time. Compared to the previous case, the distribution of stable objects on the plane of the sky is much narrower and more restricted. As in the case of low-mass asteroids, the positions of the putative planets can be seen to change rapidly relative to the background model orbits and astrometric measurements.

We note that the stability limit (red curve) in the bottom plot for the $(a_0, e_0)$-plane is offset by $\Delta a_0 = 6$ au with respect to the collision curve with planet HR 8799e shown in gray. In addition to the clear dependence of the structure of stable solutions on the probe mass, this simulation indicates that the hypothetical planet may be located only on very narrow islands in the orbital parameter space, limited to low-order resonances with the innermost planet. The temporal evolution of these islands in the sky can be useful for analyzing AO images and provides clues as to where an additional, Jovian planet might be expected. If such objects exist beyond $\approx 7$ au, they should be involved in a 2:1e, 3:1e, or 5:2e MMR with the inner planet, respectively. These low-order resonances may be preferred if the system has undergone a migration in the past. We discuss this further in Sect. 5.

Attempts to find the hypothetical innermost fifth planet have so far been unsuccessful. Very recently, Wahhaj et al. (2021), using SPHERE measurements also reported here, did not detect planet HR 8799f at the most plausible locations, namely 7.5 and 9.7 au, down to mass limits of 3.6 and 2.8 $M_{\text{Jup}}$, respectively. Neither did these authors detect any new candidate companions at the smallest observable separation, of namely 0.1" or $\approx 4.1$ au, overlapping with the semi-major axes range in Figs. 10–12. Wahhaj et al. (2021) conclude that the planet may still exist with a mass of 2–3.6 $M_{\text{Jup}}$ at 7.5 au (3:1e MMR) or 1.5–2.8 $M_{\text{Jup}}$ at 10 au, which is in the region we closely investigated in this section. The contrast curve in the interesting zone between roughly

7 and 16 au lies above roughly three Jupiter masses, as also concluded here; see Fig. 13 and discussion in Sect. 6. Therefore, we believe that steadily improved imaging techniques, gaining better contrast and lower detection limits, combined with dynamical simulations similar to those conducted in this section may eventually reveal the fifth planet. However, if a relatively massive planet of mass $\approx 1 M_{\text{Jup}}$ exists in the inner part of the system, the outer edge of the debris disk carved out by this planet would have a much smaller radius than predicted at $\approx 10$ au under the current observational configuration of four planets (Fig. 11), and related to the above-mentioned detections of warm dust at 6–10 au (e.g., Su et al. 2009; Hughes et al. 2011; Matthews et al. 2014; Chen et al. 2014). We defer the analysis of this situation to a future work as well, because of its complexity arising from different MMR scenarios.

## 5. Possible history of the system

Convergent migration due to tidal interaction with the nesting circumstellar disk is thought to be the most reliable mechanism for producing resonance trapping in multiple-planet systems (Masset & Snellgrove 2001; Lee & Peale 2002; Moorhead & Adams 2005; Thommes 2005; Beaugé et al. 2006; Crida et al. 2008; D'Angelo & Marzari 2012). The different masses of the planets and progressive inside-out depletion of gas in the disk due to the presence of the planet itself and photo-evaporation lead to different migration speeds for the planets, which may end up in resonance. Numerical modeling by Hands et al. (2014) and Szuszkiewicz & Podlewska-Gaca (2012) shows that during inward migration, planets are often trapped in resonances like the 2:1 and 3:2 (the most frequent).

While in most resonant cases, such as Gliese 876, HD 82943, and HD 37124 (Wright et al. 2011), the planets are close to the
<table><caption>**Table 7.** Osculating, heliocentric elements of the best-fitting solutions at the epoch of 1998.83 for models illustrated in Fig. 4 and in Figs. 7–8.</caption><thead><tr><th></th><th>m<sub>*</sub> = 1.47 M<sub>☉</sub>,<br>m(M<sub>Jup</sub>)</th><th>Π = 24.525662 mas,<br>a (au)</th><th>e</th><th>χ<sup>2</sup> = 1593.22,<br>I (deg)</th><th>RMS = 7.41 mas,<br>Ω (deg)</th><th>log |Δn[rad d<sup>-1</sup>]| = -12.0<br>ω (deg)</th><th>M (deg)</th></tr></thead><tbody><tr><td>Model 1:</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>HR8799e</td><td>7.4060918</td><td>16.1038912</td><td>0.1476760</td><td>26.8734271</td><td>62.1852189</td><td>110.8167067</td><td>-23.1719500</td></tr><tr><td>HR8799d</td><td>9.1899948</td><td>26.4517267</td><td>0.1146003</td><td>26.8734271</td><td>62.1852189</td><td>29.0594738</td><td>60.1752106</td></tr><tr><td>HR8799c</td><td>7.7961185</td><td>40.9861801</td><td>0.0538090</td><td>26.8734271</td><td>62.1852189</td><td>92.3879462</td><td>145.5347082</td></tr><tr><td>HR8799b</td><td>5.8023594</td><td>71.1427851</td><td>0.0166707</td><td>26.8734271</td><td>62.1852189</td><td>42.9704380</td><td>-49.1049354</td></tr><tr><td>Model 2:</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>HR8799e</td><td>7.8680000</td><td>15.9753973</td><td>0.1524708</td><td>26.2394739</td><td>63.7770393</td><td>108.3445443</td><td>335.6210695</td></tr><tr><td>HR8799d</td><td>9.6490000</td><td>26.6396632</td><td>0.1145654</td><td>26.2394739</td><td>63.7770393</td><td>28.3126706</td><td>59.7394429</td></tr><tr><td>HR8799c</td><td>7.3081300</td><td>40.9406369</td><td>0.0517422</td><td>26.2394739</td><td>63.7770393</td><td>91.3551907</td><td>145.0637721</td></tr><tr><td>HR8799b</td><td>5.6590700</td><td>71.3289844</td><td>0.0173411</td><td>26.2394739</td><td>63.7770393</td><td>41.6627885</td><td>310.5199311</td></tr><tr><td>Model 3:</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>HR8799e</td><td>8.9817200</td><td>16.1163303</td><td>0.1548634</td><td>27.3470617</td><td>62.9005413</td><td>108.0955330</td><td>336.2997163</td></tr><tr><td>HR8799d</td><td>9.0680400</td><td>26.8768634</td><td>0.1068892</td><td>27.3470617</td><td>62.9005413</td><td>29.2629241</td><td>60.0277429</td></tr><tr><td>HR8799c</td><td>7.2212900</td><td>41.0832702</td><td>0.0627005</td><td>27.3470617</td><td>62.9005413</td><td>91.9708727</td><td>144.7587008</td></tr><tr><td>HR8799b</td><td>6.1950900</td><td>71.7822044</td><td>0.0192741</td><td>27.3470617</td><td>62.9005413</td><td>44.4483815</td><td>308.9107410</td></tr><tr><td>Model 4:</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>HR8799e</td><td>10.6577500</td><td>16.1274692</td><td>0.1496852</td><td>26.9747476</td><td>63.6149718</td><td>107.6133926</td><td>336.6560551</td></tr><tr><td>HR8799d</td><td>9.5739000</td><td>26.8852716</td><td>0.1169394</td><td>26.9747476</td><td>63.6149718</td><td>27.5720349</td><td>60.0200173</td></tr><tr><td>HR8799c</td><td>6.8413200</td><td>41.2618982</td><td>0.0555455</td><td>26.9747476</td><td>63.6149718</td><td>93.1349542</td><td>143.0235158</td></tr><tr><td>HR8799b</td><td>5.7085500</td><td>71.9446080</td><td>0.0217920</td><td>26.9747476</td><td>63.6149718</td><td>45.1149621</td><td>307.6876430</td></tr><tr><td>Model 5:</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>HR8799e</td><td>7.1130900</td><td>15.9976012</td><td>0.1545297</td><td>26.5969851</td><td>62.8989533</td><td>108.8687608</td><td>336.1919760</td></tr><tr><td>HR8799d</td><td>9.8864900</td><td>26.5290674</td><td>0.1122973</td><td>26.5969851</td><td>62.8989533</td><td>27.8297522</td><td>61.2151462</td></tr><tr><td>HR8799c</td><td>7.7520500</td><td>41.1188833</td><td>0.0497536</td><td>26.5969851</td><td>62.8989533</td><td>92.4070918</td><td>144.9763772</td></tr><tr><td>HR8799b</td><td>5.5380200</td><td>71.3013150</td><td>0.0172797</td><td>26.5969851</td><td>62.8989533</td><td>44.6821435</td><td>308.5634968</td></tr><tr><td>Model 6:</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>HR8799e</td><td>7.3387300</td><td>15.9246991</td><td>0.1553777</td><td>26.5990408</td><td>62.7197878</td><td>108.3820177</td><td>336.0139713</td></tr><tr><td>HR8799d</td><td>10.0516400</td><td>26.5092120</td><td>0.1117405</td><td>26.5990408</td><td>62.7197878</td><td>29.2671129</td><td>59.8339479</td></tr><tr><td>HR8799c</td><td>8.0973100</td><td>40.9565104</td><td>0.0509428</td><td>26.5990408</td><td>62.7197878</td><td>91.8212168</td><td>145.6257159</td></tr><tr><td>HR8799b</td><td>5.8932800</td><td>71.1073028</td><td>0.0167537</td><td>26.5990408</td><td>62.7197878</td><td>41.9648749</td><td>311.2692413</td></tr></tbody></table>

star, in the case of HR 8799 the planets are significantly far away. This implies that the primordial disk from which they formed should have extended beyond 100 au which is compatible with the mass of the star being ~1.5 M<sub>☉</sub>. The most robust scenario is that in which the planets were fully formed further out in the disk and then migrated inwards until they became progressively trapped in the multiple resonances. This scenario raises the question of whether these planets were formed by core accretion (Mizuno 1980; Pollack et al. 1996) or by gravitational instability (Cameron 1978; Boss 1997).

The resonance trapping in this system at that distance from the star is confirmed by hydrodynamical simulations performed with the FARGO code (Masset 2000). Fully radiative models have been adopted where the energy equation contains viscous heating and radiative cooling through the disk surface. A polar grid with 682 × 512 elements is used to cover the disk, extending from 1 to 120 au. The initial surface gas density is given by $\Sigma = \Sigma_0 r^{-1/2}$ with $\Sigma_0 = 50 \text{ g cm}^{-2}$. This low density is motivated by the evolved state of the disk when the planets are fully formed. To test the sequential resonance trapping, we start the
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Evolution of the inner debris disk as seen on the plane of the sky. The subsequent panels show the instantaneous positions of two inner planets and small asteroids in stable orbits (small colored dots) with a mass of $10^{-6} M_{\text{Jup}}$ at the epoch labeled in the lower-left corner of each graph. We collected $3.3 \times 10^5$ particles. Their distribution in the $(a_0, e_0)$-plane of osculating, canonical elements at the initial epoch 1998.830 is shown in Fig. 11. The largest filled circles correspond to the positions of the planets HR 8799e and HR 8799d at the snapshot epoch (large filled circles) and the initial osculating epoch (shaded filled circles), respectively. The asteroid colors are for the lowest order MMR types, as indicated in the *bottom two plots*, or those belonging to the innermost quasi-homogeneous disk (gray points). The 2024.83 epoch covers roughly one orbital period in a 2:1 e MMR resonance (dark blue circles) with the innermost planet.</figcaption></figure>

inner planets on already resonant orbits; they are not affected by the disk, and so do not migrate. The outer planet instead feels the disk perturbations and migrates inward until it is trapped in the 2:1 resonance and stops migrating. This behavior is illustrated in Fig. 14 where in the top panel the semi-major axis of the fourth planet is shown while in the middle and bottom panels the critical argument of the 2:1 resonance between the third and fourth planet and the Laplace resonance argument is illustrated.

It is noteworthy that, once all the planets are trapped in resonance, they share a common gap and migrate inwards. We performed an additional simulation where all four planets are affected by the disk perturbations and can freely migrate. In this

simulation, we increased $\Sigma_0$ to $100 \text{ g/cm}^{-2}$ in order to speed up the migration. Figure 15 shows the gas density distribution after 8 Kyr from the beginning of the simulation. The planets create a common gap where some overdensities are still present in the corotation regions of each planet. The semi-major axes of the planets are illustrated in the bottom panel, rescaled to be shown in a single plot. They migrate inwards while trapped in resonance, suggesting that the present position of the planets is not necessarily the one where they became trapped in resonance. The capture in the mutual resonance may have occurred farther out, followed by an inward migration until the disk dissipates. This may push the location where the planet initially grew even
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Structure of the inner debris disk composed of $3 \times 10^5$ Ceres-like asteroids of mass $10^{-6} M_{\text{Jup}}$ on stable orbits, seen in the $(a_0, e_0)$-plane of Jacobian (canonical) osculating elements at the initial 1998.830 epoch. The colors of the asteroids indicate the lowest order MMR in which they are involved with HR 8799e, as indicated in the diagrams, or belong to the innermost, quasi-homogeneous disk (gray points) according to Fig. 10. The light gray curve indicates the collision of the orbits with HR 8799e and the red curve is an image of the collision curve shifted by $\Delta a_0 \approx 4$ au toward the star. We note that loose points above the collision curve in the top panel represent immediately scattered asteroids on regular but hyperbolic (open) orbits.</figcaption></figure>

farther out. It is noteworthy that the wider oscillations observed in the critical arguments of the resonances and the semi-major axis of the planets compared to a pure *N*-body problem are related to the presence of the perturbations of the disk on the planets.

**6. Planets–disk interaction**

Along with the four giant planets detected so far, the architecture of HR 8799 is enriched by the presence of an extended debris disk with two components. The cold Kuiper-like component was resolved in the FIR (Matthews et al. 2014) and at millimeter wavelengths (Booth et al. 2016; Wilner et al. 2018) and was detected up to 450 au from the star, with a halo extending to thousands of astronomical units (au). The warm component was instead never fully resolved and its inferred position of 9.3 au is given by modeling the IR excess at 155 K (Chen et al. 2014) in the SED of the star under the assumption of dust particles behaving like black bodies. The two components are separated by a huge dust-free gap that extends from 109 au (Wilner et al. 2018) to regions interior to the orbit of HR 8799e.

As the presence of the gap is tightly correlated with the planets in it, we can assess the planets–disk interactions and infer whether or not the planets detected are responsible for carving the entire gap or whether or not there is further dynamical space for hosting other undetected planets. For this analytical study, the only planets taken into account are the two closest to the edges of the gap, HR 8799 b and e. The orbital parameters and masses adopted are the ones listed in Tables 6 and 7 for Model 1. To calculate the extension of the region from which dust particles are scattered away from the orbit of the planet (chaotic zone), we used Eqs. (9) and (10) of Lazzoni et al. (2018) for HR 8799b and HR 8799e, respectively. As a result, we obtained that the chaotic zone of the inner planet extends down to 9.8 au, which is consistent with the direct *N*-body simulations in Sect. 4.7. Given the uncertainties on the position of the inner belt and on its width, we can safely state that HR 8799e is likely responsible for shaping

the inner edge of the gap. On the other hand, the chaotic zone for HR 8799b can clear the orbit of the planet up to 91.4 au. Given the position of the outer edge of the gap at 109 au, we can try to infer the characteristics of a fifth planet able to extend the chaotic zone up to the detected position of the edge. Following a very simple approach, we can use the same analytical tools to estimate the mass and semi-major axis of a planet able to carve a gap extending from 91.4 au to 109 au, that is, the free dynamical region left beyond the orbit of HR 8799b. As a result, we obtain a $0.12 M_{\text{Jup}}$ planet orbiting at 100 au.

Figure 13 shows a summary of the results discussed above. The four detected planets are represented as pink circles together with the extension of their chaotic zone (pink shaded area). As it emerges from the image, the inner extension of the clearing zone gets very close to the inner belt (blue vertical line) whereas it stops roughly 20 au from the outer belt. Even a very low mass planet could carve the remaining part of the gap and, consistently with our observations, it would be too small to be detected with an instrument such as SPHERE, as proven by the detection limits (red curve; IRDIS data taken on 2017 October 12) shown in the figure. Simulations of the inner debris disk in the framework of the *N*-body resonant model of the system in Sect. 4.7 indicate possible localizations of such a planet. We want to stress that this is only a preliminary analysis of the planet–disk interaction. More detailed studies will have to consider the dynamical interaction of the fifth planet with all four of the detected ones as well as the possibility of having a larger object locked in MMR with at least one of the other known planets, as discussed in the following.

Indeed, the results of numerical *N*-body simulations in Goździewski & Migaszewski (2018) based on the ⟨*Y*⟩-model described in Sect. 4.7 partially confirm the analytical predictions and address the likely real resonant configuration of the planets. These latter authors simulated the inner edge of the outer debris disk, also accounting for the presence of a hypothetical fifth planet with a mass between 0.1 and $1.66 M_{\text{Jup}}$. These simulations are based on a quasi-periodic, near-resonant orbital model
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Possible instant positions (small gray and colored dots) of a hypothetical and still undetected planet of mass $1 M_{\text{Jup}}$ on the sky plane, at the epoch labeled in the lower-left corner of each graph. We collected $10^5$ initial conditions. The bottom plot shows their distribution in the $(a_0, e_0)$-plane of canonical osculating elements at the initial epoch of 1998.830, similarly to Fig. 11. As in Fig. 10, the largest filled circles correspond to the positions of the planets HR 8799e and HR 8799d at the snapshot epoch (large filled circles) and the initial osculating epoch (shaded filled circles), respectively. Colors mark the lowest order MMR types, as labeled in the bottom plot. The light-gray curve marks the collision zone with HR 8799e, the thinner curve marks the stability zone for $m_0 = 10^{-6} M_{\text{Jup}}$ (see Fig. 10), and the red curve is the image of the geometric collision curve shifted by $\Delta a_0 \approx 6$ au toward the star.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 13.** HR 8799 architecture: two belts (blue lines), four detected planets (pink circles), and a putative fifth planet (light-blue circle). The pink shaded area represents the extension of the chaotic zone whereas the gray shaded zone represents the extension of the outer disk. The contrast curve for the system is shown in red.</figcaption></figure>

of the four known planets derived through migration simulations. However, given that Goździewski & Migaszewski (2018) adopted the larger mass for the star of $m_\star = 1.52 M_\odot$ and a larger parallax $\Pi = 25.4$ mas, the whole system appears more compact than predicted at present – these authors found the osculating astrocentric semi-major axis of HR 8799b $\approx 67.1$ au, compared to $\approx 71$ au in the present models. The updated values of the parallax and the stellar mass translate to orbits expanded by a few au. According to the simulations, the inner edge of the outer disk is globally highly nonsymmetric, and planets with a mass of between 0.1 and $1.66 M_{Jup}$ may be present beyond the clearing zone, with semi-major axis $> 90$ au. These planets could be involved in low-order resonances, such as 3:2b, 5:3b, 2:1b, or 5:2b with the outermost planet HR 8799b and in low and moderate eccentricity orbits; see Fig. 9 in Goździewski & Migaszewski (2018). The fifth planet, even if it had a small mass, would be strongly affecting the shape of the inner parts of the outer debris disk.

## 7. Conclusions

The system around HR 8799 is a unique laboratory with which to study the mutual gravitational interactions between the planets and their relation with the circumstellar disk in the early evolutionary stages of the system. To better understand the dynamics of this system, we followed up HR 8799 with SPHERE at the Very Large Telescope and with LUCI at the Large Binocular Telescope to refine the orbital parameters of the system. We reduced the new data with state-of-the-art algorithms that apply the ADI technique, and for consistency, we also repeated the reduction of published SPHERE data from open-time programs. Precise astrometry for the four planets was obtained for 21 epochs from SPHERE and one epoch from LUCI.

We performed a detailed exploration of the orbital parameters for the four planets with dynamical constraints imposed by the 8:4:2:1 Laplace resonance. This model was updated using 68 epochs from our reduction and the literature. As a result, we derived the dynamical masses of the planets and the parallax of the system with minimal prior information. The derived masses are consistent with the prediction from the evolutionary models and allow long-term stability of the orbits (over a few hundred

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Resonance trapping of the HR 8799 system during the migration phase. The upper panel shows the evolution of the semi-major axis of the outer planet at the moment of capture in resonance with the inner three bodies. After the resonance trapping, the semi-major axis remains constant. The *middle panel* shows the critical argument of the resonance between the two outer planets. After the capture in resonance, the critical argument librates. The same occurs for the Laplace critical argument shown in the *bottom panel*.</figcaption></figure>

million years), even if they are ~2 $M_{Jup}$ bigger than the masses proposed in the previous dynamical studies. We find masses of 8–9 $M_{Jup}$ for planets HR 8799 e, d, and c, while for the exterior planet HR 8799b we estimate a smaller mass of $\approx 6 M_{Jup}$. Moreover, the dynamical parallax is consistent with $1\sigma$ uncertainty with the corrected, independently determined Gaia eDR3 value of $24.50 \pm 0.05$ mas, reinforcing both the adopted mass of the host star $m_\star = 1.47 M_\odot$ and the assumed resonant or close-to-resonant configuration of the system.

Regarding the quality of the *N*-body astrometric models considered in this work, we did not find any significant or qualitative improvement of multi-parameter near-resonant configurations
<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Orbital evolution of the HR 8799 system during the migration phase. The *upper panel* shows the gas density of the disk. The *bottom panel* shows the semi-major axes of the planets rescaled with constant values for ease of visualization in the same plot.</figcaption></figure>

parametrized with 23 free orbital elements and masses compared to the previously proposed, strictly resonant periodic configuration described by 11 free parameters by GM2020. However, given that near-resonant, chaotic configurations may survive for hundreds of millions of years, it is still very difficult to differentiate such solutions from a strictly resonant configuration of the system. On the other hand, the well-defined PO N-body model with a minimal number of free parameters may be used as a template solution in numerical experiments and simulations requiring sufficiently long stable orbital evolution of the planets. These may include simulations of debris disks and yet undetected objects.

As the system shows planets locked in MMRs, we ran numerical simulations to explore the possible migration histories which could have led to the present architecture of the system. Hydrodynamical modeling suggests that the planets migrated inwards and that a different migration speed or different formation times drove the system into the present quasi-exact resonance.

We looked for the dynamical space where a putative fifth planet could sit using analytical formulae and estimates, and suggest that only a low-mass planet in the outer part of the system would agree with the currently accepted scenario. Its location would be in between planet HR 8799b and the external debris belt, and its mass would be too small ($<0.1 M_{Jup}$) to be detected by current instrumentation for high-contrast imaging. On the other hand, the interior belt is heavily shaped by

planet HR 8799e, and there is no need for an additional inner planet to explain it. However, we also determined the dynamical structure of this region with direct N-body simulations. The results may be useful for predicting the positions of small-mass objects below the present detection limits of $\approx 3 M_{Jup}$.

**Acknowledgements.** We are thankful to the referee for the constructive report. A.Z. acknowledges support from the FONDECYT Iniciación en investigación project number 11190837 and ANID - Millennium Science Initiative Program - Center Code NCN2021_080. K.G. is very grateful to Dr Cezary Migaszewski for sharing a template C code for computing periodic orbits in a few planet systems and explanations regarding the PO approach. K.G. thanks the staff of the Poznań Supercomputer and Network Centre (PCSS, Poland) for the generous long-term support and computing resources (grant No. 529). K.G. also thanks the staff of the Tricity Supercomputer Centre, Gdańsk (Poland) for computing resources on the Tryton cluster that were used to conduct some numerical experiments and for their great and professional support. A.V. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme, grant agreements no. 757561 (HiRISE). The LBTO AO group would like to acknowledge the assistance of R.T. Gatto with night observations. We would also like to thank J. Power for assisting in preparation of the observations using ADI mode. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any author accepted manuscript version arising from this submission. SPHERE is an instrument designed and built by a consortium consisting of IPAG (Grenoble, France), MPIA (Heidelberg, Germany), LAM (Marseille, France), LESIA (Paris, France), Laboratoire Lagrange (Nice, France), INAF - Osservatorio di Padova (Italy), Observatoire de Genève (Switzerland), ETH Zurich (Switzerland), NOVA (Netherlands), ONERA (France), and ASTRON (Netherlands), in collaboration with ESO. SPHERE was funded by ESO, with additional contributions from CNRS (France), MPIA (Germany), INAF (Italy), FINES (Switzerland), and NOVA (Netherlands). SPHERE also received funding from the European Commission Sixth and Seventh Framework Programmes as part of the Optical Infrared Coordination Network for Astronomy (OPTICON) under grant number RII3-Ct-2004-001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016).

## References

Apai, D., Kasper, M., Skemer, A., et al. 2016, ApJ, 820, 40

Baines, E. K., White, R. J., Huber, D., et al. 2012, ApJ, 761, 57

Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701

Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42

Beaugé, C., Michtchenko, T. A., & Ferraz-Mello, S. 2006, MNRAS, 365, 1160

Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593

Bergfors, C., Brandner, W., Janson, M., Köhler, R., & Henning, T. 2011, A&A, 528, A134

Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Biller, B. A., Apai, D., Bonnefoy, M., et al. 2021, MNRAS, 503, 743

Boccaletti, A., Abe, L., Baudrand, J., et al. 2008, SPIE Conf. Ser., 7015, 70151B

Bohn, A. J., Kenworthy, M. A., Ginski, C., et al. 2020, ApJ, 898, L16

Bohn, A. J., Ginski, C., Kenworthy, M. A., et al. 2021, A&A, 648, A73

Booth, M., Jordan, A., Casassus, S., et al. 2016, MNRAS, 476, 2845

Boss, A. P. 1997, Science, 276, 1836

Brandt, G. M., Brandt, T. D., Dupuy, T. J., Michalik, D., & Marleau, G.-D. 2021, ApJ, 915, L16

Cameron, A. G. W. 1978, Moon and Planets, 18, 5

Chauvin, G., Lagrange, A.-M., Dumas, C., et al. 2004, A&A, 425, L29

Chen, C. H., Mittal, T., Kuchner, M., et al. 2014, ApJ, 211, 25

Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Physica D Nonlinear Phenomena, 182, 151

Claudi, R. U., Turatto, M., Gratton, R. G., et al. 2008, SPIE Conf. Ser., 7014, 70143E

Crida, A., Sándor, Z., & Kley, W. 2008, A&A, 483, 325

Currie, T., Burrows, A., Itoh, Y., et al. 2011, ApJ, 729, 128

Currie, T., Fukagawa, M., Thalmann, C., Matsumura, S., & Plavchan, P. 2012, ApJ, 755, L34

Currie, T., Burrows, A., Girard, J. H., et al. 2014, ApJ, 795, 133

D'Angelo, G., & Marzari, F. 2012, ApJ, 757, 50

De Rosa, R. J., Esposito, T. M., Gibbs, A., et al. 2020, SPIE Conf. Ser., 11447, 114475A

Delisle, J.-B. 2017, A&A, 605, A96

Delorme, P., Meunier, N., Albert, D., et al. 2017, in SF2A-2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. C. Reylé, P. Di Matteo, F. Herpin, et al.
Desidera, S., Chauvin, G., Bonavita, M., et al. 2021, A&A, 651, A70

Dohlen, K., Langlois, M., Saisse, M., et al. 2008, SPIE Conf. Ser., 7014, 70143L

Esposito, S., Riccardi, A., Fini, L., et al. 2010, SPIE Conf. Ser., 7736, 107

Esposito, S., Mesa, D., Skemer, A., et al. 2013, A&A, 549, A52

Fabrycky, D. C., & Murray-Clay, R. A. 2010, ApJ, 710, 1408

Faramaz, V., Marino, S., Booth, M., et al. 2021, AJ, 161, 271

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Fukagawa, M., Itoh, Y., Tamura, M., et al. 2009, ApJ, 696, L1

Fusco, T., Rousset, G., Sauvage, J.-F., et al. 2006, Opt. Express, 14, 7515

Gaia Collaboration 2020, VizieR Online Data Catalog: I/350

Galicher, R., Marois, C., Macintosh, B., Barman, T., & Konopacky, Q. 2011, ApJ, 739, L41

Galicher, R., Boccaletti, A., Mesa, D., et al. 2018, A&A, 615, A92

Goździewski, K., & Migaszewski, C. 2018, ApJS, 238, 6

Goździewski, K., & Migaszewski, C. 2020, ApJ, 902, L40

Goździewski, K., Bois, E., Maciejewski, A. J., & Kiseleva-Eggleton, L. 2001, A&A, 378, 569

GRAVITY Collaboration (Lacour, S., et al.) 2019, A&A, 623, A11

Gray, R. O., & Kaye, A. B. 1999, AJ, 118, 2993

Hadjidemetriou, J. D. 1977, Celest. Mech., 16, 61

Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749

Hairer, E., Nørsett, S., & Wanner, G. 2000, Solving Ordinary Differential Equations I Nonstiff Problems, 2nd edn. (Berlin: Springer)

Hands, T. O., Alexander, R. D., & Dehnen, W. 2014, MNRAS, 445, 749

Hinz, P. M., Rodigas, T. J., Kenworthy, M. A., et al. 2010, ApJ, 716, 417

Hughes, A. M., Wilner, D. J., Andrews, S. M., et al. 2011, ApJ, 740, 38

Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44

Kervella, P., Arenou, F., & Thévenin, F. 2022, A&A, 657, A7

Konopacky, Q. M., Marois, C., Macintosh, B. A., et al. 2016, AJ, 152, 28

Lafrenière, D., Marois, C., Doyon, R., & Barman, T. 2009, ApJ, 694, L148

Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57

Langlois, M., Vigan, A., Moutou, C., et al. 2013, in Proceedings of the Third AO4ELT Conference, ed. S. Esposito & L. Fini, 63

Langlois, M., Gratton, R., Lagrange, A. M., et al. 2021, A&A, 651, A71

Laskar, J. 1993, Celest. Mech. Dyn. Astron., 56, 191

Lazzoni, C., Desidera, S., Marzari, F., et al. 2018, A&A, 611, A43

Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596

Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4

Maire, A. L., Skemer, A. J., Hinz, P. M., et al. 2015, A&A, 579, A2

Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, SPIE Conf. Ser., 9908, 990834

Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556

Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348

Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 2010, Nature, 468, 1080

Marois, C., Correia, C., Véran, J.-P., & Currie, T. 2014, in IAU Symposium, 299, IAU Symposium, eds. M. Booth, B. C. Matthews, & J. R. Graham, 48

Masset, F. 2000, A&AS, 141, 165

Masset, F., & Snellgrove, M. 2001, MNRAS, 320, L55

Matthews, B., Kennedy, G., Sibthorpe, B., et al. 2014, ApJ, 780, 97

Mizuno, H. 1980, Progr. Theor. Phys., 64, 544

Moorhead, A. V., & Adams, F. C. 2005, Icarus, 178, 517

Nesvorny, D., & Ferraz-Mello, S. 1997, A&A, 320, 672

Nowak, M., Lacour, S., Lagrange, A. M., et al. 2020, A&A, 642, A2

Papaloizou, J. C. B. 2015, Int. J. Astrobiol., 14, 291

Petit, C., Sauvage, J.-F., Fusco, T., et al. 2014, SPIE Conf. Ser., 9148, 0

Pinna, E., Rossi, F., Puglisi, A., et al. 2019, in Adaptive Optics for Extremely Large Telescopes 6 – Conference Proceedings

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62

Rameau, J., Chauvin, G., Lagrange, A.-M., et al. 2013, ApJ, 772, L15

Ramos, X. S., Charalambous, C., Benítez-Llambay, P., & Beaugé, C. 2017, A&A, 602, A101

Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424

Seifert, W., Ageorges, N., Lehmitz, M., et al. 2010, in Proc. SPIE, Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III, 77357W

Sepulveda, A. G., & Bowler, B. P. 2022, AJ, 163, 52

Soummer, R., Hagan, J. B., Pueyo, L., et al. 2011, ApJ, 741, 55

Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28

Su, K. Y. L., Rieke, G. H., Stapelfeldt, K. R., et al. 2009, ApJ, 705, 314

Sudol, J. J., & Haghighipour, N. 2012, ApJ, 755, 38

Szuszkiewicz, E., & Podlew ska-Gaca, E. 2012, *Origins Life Evol. Biosphere*, 42, 113

Thommes, E. W. 2005, ApJ, 626, 1033

Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72

Wahhaj, Z., Milli, J., Romero, C., et al. 2021, A&A, 648, A26

Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, AJ, 156, 192

Wertz, O., Absil, O., Gómez González, C. A., et al. 2017, A&A, 598, A83

Wilner, D. J., MacGregor, M. A., Andrews, S. M., et al. 2018, ApJ, 855, 56

Wright, J. T., Veras, D., Ford, E. B., et al. 2011, ApJ, 730, 93

Zuckerman, B., Rhee, J. H., Song, I., & Bessell, M. S. 2011, ApJ, 732, 61

Zurlo, A., Vigan, A., Mesa, D., et al. 2014, A&A, 572, A85

Zurlo, A., Vigan, A., Galicher, R., et al. 2016, A&A, 587, A57

<footnote n="1">¹ Núcleo de Astronomía, Facultad de Ingeniería, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile
e-mail: alice.zurlo@mail.udp.cl</footnote>

<footnote n="2">² Millennium Cluster on Young Exoplanets and their Moons (YEMS), Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile</footnote>

<footnote n="3">³ Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, 13388 Marseille, France</footnote>

<footnote n="4">⁴ Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Toruń, Poland</footnote>

<footnote n="5">⁵ University of Exeter, Astrophysics Group, Physics Building, Stocker Road, Exeter EX4 4QL, UK</footnote>

<footnote n="6">⁶ INAF – Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy</footnote>

<footnote n="7">⁷ Department of Physics and Astronomy, University of Padova, Via Marzolo 8, 35131 Padova, Italy</footnote>

<footnote n="8">⁸ CRAL, UMR 5574, CNRS, Université de Lyon, ENS, 9 avenue Charles André, 69561 Saint-Genis-Laval Cedex, France</footnote>

<footnote n="9">⁹ INAF – Osservatorio Astrofisico di Arcetri, L.go E. Fermi 5, 50125 Firenze, Italy</footnote>

<footnote n="10">¹⁰ IPAG, Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="11">¹¹ Space Telescope Science Institute (STScI), 3700 San Martin Dr, Baltimore, MD 21218, USA</footnote>

<footnote n="12">¹² Geneva Observatory, University of Geneva, 51 ch. Pegasi, 1290 Versiox, Switzerland</footnote>

<footnote n="13">¹³ Max Planck Institute for Astronomy, Königstuhl 17, 69117, Heidelberg, Germany</footnote>

<footnote n="14">¹⁴ Department of Astronomy, Stockholm University, 106 91 Stockholm, Sweden</footnote>

<footnote n="15">¹⁵ European Space Agency (ESA), ESA Office, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="16">¹⁶ LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195 Meudon, France</footnote>

<footnote n="17">¹⁷ SUPA, Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="18">¹⁸ School of Physical Sciences, Faculty of Science, Technology, Engineering and Mathematics, The Open University, Walton Hall, Milton Keynes MK7 6AA, UK</footnote>

<footnote n="19">¹⁹ Institute for Particle Physics and Astrophysics, ETH Zurich, Wolfgang-Pauli-Strasse 27, 8093 Zurich, Switzerland</footnote>

<footnote n="20">²⁰ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, 06304 Nice, France</footnote>

<footnote n="21">²¹ European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany</footnote>

<footnote n="22">²² LBT Observatory, University of Arizona, 933 N.Cherry Ave, Tucson AZ 85721, USA</footnote>

<footnote n="23">²³ George Mason University, Department of Physics & Astronomy, MS 3F3, 4400 University Drive, Fairfax, VA 22030, USA</footnote>

<footnote n="24">²⁴ European Southern Observatory, Alonso de Cordova 3107, Casilla 19001 Vitacura, Santiago 19, Chile</footnote>
<table><caption>Table A.1: List of the astrometric points.</caption><thead><tr><th rowspan="2">Date</th><th colspan="3">Planet b</th><th colspan="3">Planet c</th><th colspan="3">Planet d</th><th colspan="3">Planet e</th><th rowspan="2">Ref.</th></tr><tr><th>ΔRA</th><th>σ<sub>ΔRA</sub></th><th>ΔDec</th><th>σ<sub>ΔDec</sub></th><th>ΔRA</th><th>σ<sub>ΔRA</sub></th><th>ΔDec</th><th>σ<sub>ΔDec</sub></th><th>ΔRA</th><th>σ<sub>ΔRA</sub></th><th>ΔDec</th><th>σ<sub>ΔDec</sub></th><th>ΔRA</th><th>σ<sub>ΔRA</sub></th><th>ΔDec</th><th>σ<sub>ΔDec</sub></th></tr></thead><tbody><tr><td>1998.83</td><td>1411</td><td>9</td><td>986</td><td>9</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Latfremière et al. (2009)</td></tr><tr><td>1998.83</td><td>1418</td><td>22</td><td>1004</td><td>20</td><td>-837</td><td>-</td><td>26</td><td>483</td><td>23</td><td>133</td><td>35</td><td>-533</td><td>34</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Soummer et al. (2011)</td></tr><tr><td>2002.54</td><td>1481</td><td>23</td><td>919</td><td>17</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Fukagawa et al. (2009)</td></tr><tr><td>2004.53</td><td>1471</td><td>6</td><td>884</td><td>6</td><td>-739</td><td>6</td><td>612</td><td>6</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2005.54</td><td>1496</td><td>5</td><td>856</td><td>5</td><td>-713</td><td>5</td><td>630</td><td>5</td><td>-87</td><td>10</td><td>-578</td><td>10</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Currie et al. (2012)</td></tr><tr><td>2007.58</td><td>1504</td><td>3</td><td>837</td><td>3</td><td>-683</td><td>4</td><td>671</td><td>4</td><td>-179</td><td>5</td><td>-588</td><td>5</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2007.81</td><td>1500</td><td>7</td><td>836</td><td>7</td><td>-678</td><td>7</td><td>676</td><td>7</td><td>-175</td><td>10</td><td>-589</td><td>10</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2008.52</td><td>1527</td><td>4</td><td>799</td><td>4</td><td>-658</td><td>4</td><td>701</td><td>4</td><td>-208</td><td>4</td><td>-582</td><td>4</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Marois et al. (2008)</td></tr><tr><td>2008.61</td><td>1527</td><td>2</td><td>801</td><td>2</td><td>-657</td><td>2</td><td>706</td><td>2</td><td>-216</td><td>2</td><td>-582</td><td>2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Marois et al. (2008)</td></tr><tr><td>2008.71</td><td>1516</td><td>4</td><td>818</td><td>4</td><td>-663</td><td>3</td><td>693</td><td>3</td><td>-202</td><td>4</td><td>-588</td><td>4</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2008.89</td><td>1532</td><td>20</td><td>796</td><td>20</td><td>-634</td><td>20</td><td>700</td><td>20</td><td>-217</td><td>20</td><td>-608</td><td>20</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Hinz et al. (2010)</td></tr><tr><td>2009.02</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-612</td><td>30</td><td>665</td><td>30</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Hinz et al. (2010)</td></tr><tr><td>2009.58</td><td>1526</td><td>4</td><td>797</td><td>4</td><td>-639</td><td>4</td><td>712</td><td>4</td><td>-237</td><td>3</td><td>-577</td><td>3</td><td>-306</td><td>7</td><td>-211</td><td>7</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2009.62</td><td>1536</td><td>10</td><td>785</td><td>10</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Currie et al. (2011)</td></tr><tr><td>2009.70</td><td>1538</td><td>30</td><td>777</td><td>30</td><td>-634</td><td>30</td><td>697</td><td>30</td><td>-282</td><td>30</td><td>-590</td><td>30</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Hinz et al. (2010)</td></tr><tr><td>2009.76</td><td>1535</td><td>20</td><td>816</td><td>20</td><td>-636</td><td>40</td><td>692</td><td>40</td><td>-270</td><td>70</td><td>-600</td><td>70</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Bergfors et al. (2011)</td></tr><tr><td>2009.77</td><td>1532</td><td>7</td><td>783</td><td>7</td><td>-627</td><td>7</td><td>716</td><td>7</td><td>-241</td><td>7</td><td>-586</td><td>7</td><td>-306</td><td>7</td><td>-217</td><td>7</td><td>-</td><td>Currie et al. (2011)</td></tr><tr><td>2009.83</td><td>1524</td><td>10</td><td>795</td><td>10</td><td>-636</td><td>9</td><td>720</td><td>9</td><td>-251</td><td>7</td><td>-573</td><td>7</td><td>-310</td><td>9</td><td>-187</td><td>9</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2009.84</td><td>1540</td><td>19</td><td>800</td><td>19</td><td>-630</td><td>13</td><td>720</td><td>13</td><td>-240</td><td>14</td><td>-580</td><td>14</td><td>-304</td><td>10</td><td>-196</td><td>10</td><td>-</td><td>Galicher et al. (2011)</td></tr><tr><td>2010.53</td><td>1532</td><td>5</td><td>783</td><td>5</td><td>-619</td><td>4</td><td>728</td><td>4</td><td>-265</td><td>4</td><td>-576</td><td>4</td><td>-323</td><td>6</td><td>-166</td><td>6</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2010.55</td><td>1547</td><td>6</td><td>757</td><td>9</td><td>-606</td><td>6</td><td>725</td><td>6</td><td>-269</td><td>6</td><td>-580</td><td>6</td><td>-329</td><td>6</td><td>-178</td><td>6</td><td>-</td><td>Currie et al. (2014)</td></tr><tr><td>2010.83</td><td>1535</td><td>15</td><td>766</td><td>15</td><td>-607</td><td>12</td><td>744</td><td>12</td><td>-296</td><td>13</td><td>-561</td><td>13</td><td>-341</td><td>16</td><td>-143</td><td>16</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2011.55</td><td>1541</td><td>5</td><td>762</td><td>5</td><td>-595</td><td>4</td><td>747</td><td>4</td><td>-303</td><td>5</td><td>-562</td><td>5</td><td>-352</td><td>8</td><td>-130</td><td>8</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2011.79</td><td>1579</td><td>11</td><td>734</td><td>11</td><td>-561</td><td>13</td><td>752</td><td>13</td><td>-299</td><td>13</td><td>-563</td><td>13</td><td>-326</td><td>13</td><td>-119</td><td>13</td><td>-</td><td>Esposito et al. (2013)</td></tr><tr><td>2011.86</td><td>1546</td><td>11</td><td>725</td><td>11</td><td>-578</td><td>13</td><td>767</td><td>13</td><td>-320</td><td>13</td><td>-549</td><td>13</td><td>-382</td><td>16</td><td>-127</td><td>16</td><td>-</td><td>Esposito et al. (2013)</td></tr><tr><td>2012.55</td><td>1545</td><td>5</td><td>747</td><td>5</td><td>-578</td><td>5</td><td>761</td><td>5</td><td>-339</td><td>5</td><td>-555</td><td>5</td><td>-373</td><td>8</td><td>-84</td><td>8</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2012.82</td><td>1549</td><td>4</td><td>743</td><td>4</td><td>-572</td><td>3</td><td>768</td><td>3</td><td>-346</td><td>4</td><td>-548</td><td>4</td><td>-370</td><td>9</td><td>-76</td><td>9</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2012.83</td><td>1558</td><td>6</td><td>729</td><td>9</td><td>-557</td><td>6</td><td>763</td><td>6</td><td>-343</td><td>6</td><td>-555</td><td>6</td><td>-371</td><td>6</td><td>-80</td><td>6</td><td>-</td><td>Currie et al. (2014)</td></tr><tr><td>2013.79</td><td>1545</td><td>22</td><td>724</td><td>22</td><td>-542</td><td>22</td><td>784</td><td>22</td><td>-382</td><td>16</td><td>-522</td><td>16</td><td>-373</td><td>13</td><td>-17</td><td>13</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2013.81</td><td>1562</td><td>8</td><td>713</td><td>13</td><td>-538</td><td>6</td><td>784</td><td>13</td><td>-377</td><td>7</td><td>-538</td><td>11</td><td>-394</td><td>11</td><td>-36</td><td>17</td><td>-</td><td>Maire et al. (2015)</td></tr><tr><td>2013.88</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-537</td><td>1</td><td>782</td><td>2</td><td>-370</td><td>1</td><td>-539</td><td>1</td><td>-381</td><td>2</td><td>-30</td><td>0</td><td>-</td><td>Wang et al. (2018)<sup>a</sup></td></tr><tr><td>2014.53</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-400</td><td>4</td><td>-512</td><td>4</td><td>-389</td><td>1</td><td>-22</td><td>IFS, This work</td></tr><tr><td>2014.53</td><td>1570</td><td>3</td><td>704</td><td>3</td><td>-521</td><td>3</td><td>790</td><td>9</td><td>-391</td><td>2</td><td>-530</td><td>2</td><td>-387</td><td>2</td><td>-10</td><td>3</td><td>-</td><td>IRDIS, This work</td></tr><tr><td>2014.54</td><td>1560</td><td>13</td><td>725</td><td>13</td><td>-540</td><td>13</td><td>799</td><td>13</td><td>-400</td><td>11</td><td>-534</td><td>11</td><td>-387</td><td>11</td><td>3</td><td>11</td><td>-</td><td>Konopacky et al. (2016)</td></tr><tr><td>2014.62</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-396</td><td>1</td><td>-524</td><td>2</td><td>-389</td><td>1</td><td>-17</td><td>IFS, This work</td></tr><tr><td>2014.70</td><td>1569</td><td>4</td><td>707</td><td>2</td><td>-519</td><td>1</td><td>794</td><td>2</td><td>-397</td><td>1</td><td>-530</td><td>2</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>Wang et al. (2018)<sup>1</sup></td></tr><tr><td>2014.93</td><td>1575</td><td>2</td><td>702</td><td>4</td><td>-511</td><td>2</td><td>799</td><td>2</td><td>-400</td><td>2</td><td>-523</td><td>2</td><td>-385</td><td>3</td><td>12</td><td>2</td><td>-</td><td>Wertz et al. (2017)</td></tr><tr><td>2014.93</td><td>1574</td><td>3</td><td>701</td><td>2</td><td>-514</td><td>3</td><td>798</td><td>4</td><td>-399</td><td>4</td><td>-525</td><td>4</td><td>-389</td><td>8</td><td>11</td><td>4</td><td>-</td><td>IRDIS, This work</td></tr></tbody></table>

<footnote>a Revised following the updated astrometric calibration in De Rosa et al. (2020)</footnote>