Issue 
A&A
Volume 666, October 2022



Article Number  A133  
Number of page(s)  25  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202243862  
Published online  19 October 2022 
Orbital and dynamical analysis of the system around HR 8799
New astrometric epochs from VLT/SPHERE and LBT/LUCI
^{1}
Núcleo de Astronomía, Facultad de Ingeniería, Universidad Diego Portales,
Av. Ejercito 441,
Santiago, Chile
email: alice.zurlo@mail.udp.cl
^{2}
Millennium Nucleus on Young Exoplanets and their Moons (YEMS), Universidad Diego Portales,
Av. Ejercito 441,
Santiago, Chile
^{3}
AixMarseille Université, CNRS, LAM (Laboratoire d’Astrophysique de Marseille) UMR 7326,
13388
Marseille, France
^{4}
Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University,
Grudziadzka 5,
87–100
Toruń, Poland
^{5}
University of Exeter, Astrophysics Group, Physics Building,
Stocker Road,
Exeter
EX4 4QL, UK
^{6}
INAF – Osservatorio Astronomico di Padova,
Vicolo dell’Osservatorio 5,
35122
Padova, Italy
^{7}
Department of Physics and Astronomy, University of Padova,
Via Marzolo 8,
35131
Padova, Italy
^{8}
CRAL, UMR 5574, CNRS, Université de Lyon, ENS,
9 avenue Charles André,
69561
SaintGenisLaval Cedex, France
^{9}
INAF – Osservatorio Astrofisico di Arcetri,
L.go E. Fermi 5,
50125
Firenze, Italy
^{10}
IPAG, Univ. Grenoble Alpes, CNRS, IPAG,
38000
Grenoble, France
^{11}
Space Telescope Science Institute (STScI),
3700 San Martin Dr,
Baltimore, MD
21218, USA
^{12}
Geneva Observatory, University of Geneva,
51 ch. Pegasi,
1290
Versoix, Switzerland
^{13}
Max Planck Institute for Astronomy,
Königstuhl 17,
69117
Heidelberg, Germany
^{14}
Department of Astronomy, Stockholm University,
106 91
Stockholm, Sweden
^{15}
European Space Agency (ESA), ESA Office, Space Telescope Science Institute,
3700 San Martin Drive,
Baltimore, MD
21218, USA
^{16}
LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris,
5 place Jules Janssen,
92195
Meudon, France
^{17}
SUPA, Institute for Astronomy, University of Edinburgh,
Blackford Hill,
Edinburgh
EH9 3HJ, UK
^{18}
School of Physical Sciences, Faculty of Science, Technology, Engineering and Mathematics, The Open University,
Walton Hall,
Milton Keynes
MK7 6AA, UK
^{19}
Institute for Particle Physics and Astrophysics, ETH Zurich,
WolfgangPauliStrasse 27,
8093
Zurich, Switzerland
^{20}
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange,
06304
Nice, France
^{21}
European Southern Observatory,
KarlSchwarzschildStrasse 2,
85748
Garching, Germany
^{22}
LBT Observatory, University of Arizona,
933 N.Cherry Ave,
Tucson AZ
85721, USA
^{23}
George Mason University, Department of Physics & Astronomy,
MS 3F3,
4400 University Drive,
Fairfax, VA
22030, USA
^{24}
European Southern Observatory,
Alonso de Cordova 3107, Casilla
19001
Vitacura, Santiago 19, Chile
Received:
25
April
2022
Accepted:
15
June
2022
Context. HR 8799 is a young planetary system composed of four planets and a double debris belt. Being the first multiplanetary 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 planettodisk interaction, we made followup 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 stateoftheart algorithms implemented to apply the spectral and angular differential imaging method. We rereduced the SPHERE data obtained during the commissioning of the instrument and in three opentime 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 Nbody simulations indicating possible positions of a putative fifth planet with a mass below the present detection limits of ≃3 M_{Jup}.
Key words: planets and satellites: dynamical evolution and stability / planetdisk interactions / stars: individual: HR8799 / instrumentation: adaptive optics / astrometry / techniques: image processing
© A. Zurlo et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
Among the exoplanets discovered with the highcontrast 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 HR8799 is undoubtedly one of the most interesting. This is mainly because HR 8799 hosts a greater number of planets detected with highcontrast 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 89987601 (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 (HR8799bcde; 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), γ Dortype variable star (Gray & Kaye 1999) with λ Boolike abundance patterns. The mass of the star is ${1.47}_{0.17}^{+0.12}$ M_{⊙} (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 highcontrast imaging observations, as its four planets are easily detectable with stateoftheart imagers; their contrasts are about 2–8 × 10^{−6}; 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 multiplanetary 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 & MurrayClay 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 hotstart 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 HR8799bcde. Ad hoc observations were designed for the orbital followup 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 planettodisk 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 planetsdisk 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 highcontrast imaging. The nearinfrared (NIR) arm includes the IR dualband 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 (λ_{H2} = 1.593 µm, λ_{H3} = 1.667 µm), in BB_H (λ_{H} = 1.625 µm), and K1K2 bands (λ_{K1} =2.102 µm, λ_{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 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 offaxis 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 signaltonoise 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.
Summary of the observations of HR 8799 with IRDIS and IFS during SHINE.
2.2 Large Binocular Telescope/SOULLUCI
HR 8799 was observed with LBT/SOULLUCI1 during its commissioning on night 2020 September 29. LUCI1 (Seifert et al. 2010) is a NIR spectroimager that can work in a diffractionlimited regime with a sampling of 15 mas pix^{−1}. The AO 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 subwindowing 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 narrowband filters: FeII (λ_{FeII} = 1.646 µm; Δλ = 0.018 µm) and H2 (λ_{H2} = 2.124 µm; Δλ = 0.023 µm). The observations with the two filters were alternated during the whole sequence.
List of astrometric points from IRDIS observations.
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 TLOCI (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 reprocessed 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 was switched off to allow the rotation of the FOV and be able to use the ADI method. Before performing the highcontrast 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.
List of astrometric points from IFS observations.
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; GoZdziewski & Migaszewski 2020, and references therein), it is now widely recognized that the present observation timewindow does not make it possible to determine longterm 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 meanmotion 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 highcontrast images and resulting orbital solutions in Wahhaj et al. (2021) indicate that planets HR8799e and HR8799c are most consistent with coplanar 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 exactresonance 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 longterm 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 nearresonance models can yield statistically different or perhaps even better bestfitting solutions than those that appear in the PO scenario. We define the following merit functions as in Goździewski & Migaszewski (2018, 2020): $$\begin{array}{c}\mathrm{ln}\mathcal{L}(x)=\frac{1}{2}{\chi}^{2}(x)\frac{1}{2}{\displaystyle \sum _{i=1}^{{N}_{\text{obs}}}\left[\mathrm{ln}{\theta}_{i,\alpha}^{2}+\mathrm{ln}{\theta}_{i,\delta}^{2}\right]}{N}_{\text{obs}}\mathrm{ln}(2\pi ),\\ {\chi}^{2}(x)={\displaystyle \sum _{i=1}^{{N}_{\text{obs}}}\left[\frac{{[{\alpha}_{i}\alpha ({t}_{i},x)]}^{2}}{{\theta}_{i,\alpha}^{2}}+\frac{{[{\delta}_{i}\delta ({t}_{i},x)]}^{2}}{{\theta}_{i,\delta}^{2}}\right]},\end{array}$$(1)
where (α_{i}, δ_{i}) are the right ascension (RA) and declination (Dec) measurements at time t_{i}; α(t_{i}, x), δ(t_{i}, x) are for their ephemeris (model) values at the same moments as implied by the adopted 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 socalled 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_{obs} is the number of observations then N = 2N_{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 bestfitting model should yield the reduced ${\chi}_{v}^{2}\simeq 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_{⋆} = 1.52 M_{⊙} determined by Baines et al. (2012) as ${m}_{\star}={1.516}_{0.024}^{+0.038}$ M_{⊙} for the age of ≃33 Myr, and ${m}_{\star}={1.47}_{0.17}^{+0.12}$ M_{⊙}, 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_{⋆} = 1.50 ± 0.08 M_{⊙}, which overlaps with the two estimates. We note here that the stellar mass is fixed in our astrometric Nbody model because it also constrains the parallax Π of the whole system. These two parameters are strongly correlated through the nearIII Keplerian law. Moreover, the parallax tends to be systematically tightly bounded in subsequent Gaia catalogs. The most recent Gaia eDR3 estimate of Π = 24.460 ± 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 Π = 24.462 ± 0.046 mas.
We also implied Gaussian priors N_{m} = [8.7 ± 1.7, 8.7 ± 1.0, 8.7 ± 1.0, 6 ± 0.7] M_{Jup} for the masses of the innermost to the outermost planet, respectively, based on the hotstart 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 HR8799e and HR 8799c appear strongly anticorrelated. Masses for D_K2 seem to be too far into the highend 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 ≃7 M_{Jup}, given the results of dynamical Nbody simulations. In particular, Wang et al. (2018) report difficulty in finding dynamically stable solutions for planet masses larger than ≃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, ≃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 ≃9 M_{Jup} range. Brandt et al. (2021) estimated the mass of HR8799e as ${9.6}_{1.8}^{+1.9}$ 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 ± 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 Π = 24.46 ± 0.05 mas. This value is significantly shifted from Π ≃ 24.22 ± 0.08 mas in the Gaia DR2 and 24.8 ± 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 ± 0.05 mas (Kervella et al. 2022). This apparently very small shift of the order of 1σ uncertainty may still be meaningful when compared to the dynamical estimates, as is found below.
Values for the mass of each planet calculated using the evolutionary models in the different IRDIS filters.
List of astrometric points from the LBT/LUCI observation.
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 , …, 4, the osculating period ratio κ for the two innermost planets – which selects the particular fourbody Laplace resonance chain, as well as the reference epoch τ, three Euler angles (I, ω_{rot}, Ω) rotating the coplanar resonant configuration to the sky (observer) plane, and the system parallax Π. This means 11 free parameters. We note that all remaining orbital parameters are constrained through the Nbody 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 coplanar Nplanet system, as described in GM2020 and implemented by Cezary Migaszewski (priv. comm.). To conduct the MCMC sampling, we used the affine sampler in emcee package (ForemanMackey et al. 2013). We initiated 1536 walkers in a small hyperball in the parameter space around the initial condition in Table 6 determined with a search with the Powell nongradient 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 2dim projections and 1dim 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} anticorrelation, there is also m_{e}−κ correlation and m_{d} − κ anticorrelation. 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 bestfitting 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_{PO} for the two stellar masses of m_{⋆} = 1.52 M_{⊙} and m_{⋆} = 1.47 M_{⊙}, respectively.
For a stellar mass of m_{⋆} = 1.47M_{⊙}, the HR8799d mass seems to be constrained to ≃9.2 M_{Jup} within a relatively very small uncertainty of ≃0.1 M_{Jup}. This mass estimate is slightly larger than that found by GM2020. The mass of the outermost planet HR 8799b may be determined as ≃5.8 M_{Jup} which is similar to ≃5.5 M_{Jup} found by these latter authors, yet with also smaller uncertainty of ≃0.3 M_{Jup}. Moreover, for the larger stellar mass of 1.52 M_{⊙}, we derived slightly larger masses of HR 8799d ≃9.6 ± 0.1 M_{Jup} and HR 8799b ≃6.0 ± 0.4 M_{Jup}.
The astrometric model for the updated data window yields Π ≃ 24.26 ± 0.05 mas, which is consistent with the earlier estimate by GM2020. The parallax estimate of Π ≃ 24.53 ± 0.04 mas derived for 1.47 M_{⊙} overlaps to 1σ with Π ≃ 24.460 ± 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 Π = 24.50 ± 0.05 mas which is even closer to the value derived from the PO model. This relation confirms the predicted strong stellar massparallax correlation and indirectly favors m_{⋆} = 1.47 M_{⊙}. 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 anticorrelation, 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 ≃20 yr around the bestfitting periodic model with deviations and uncertainties of ~1 mas. It appears that all correlations except those between the masses and κ 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.
Osculating, heliocentric elements of the bestfitting solutions at the epoch of 1998.830.
Fig. 1 MCMC posterior for the bestfitting, strictly resonant model derived for m_{⋆} = 1.47 M_{Sun}. Parameters included in the diagram are for the planet masses and orbital period ratio κ = P_{d}/P_{e} = κ, and parallax Π, inclination I, and the nodal angle Ω (the orbital scale ρ, PO epoch shift t_{0} and the common rotation angle ω_{rot} in the orbital planet are skipped). The crossed lines mark the bestfitting PO solution in terms of the smallest χ^{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. 
Fig. 2 Temporal evolution of the osculating orbital elements in the astrocentric (blue curves) and Jacobian (yellow curves) frames, respectively, for the same Nbody 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. 
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 Nbody dynamics. We resolve these frequencies with the refined Fourier frequency analysis (e.g., Laskar 1993; Nesvorny & FerrazMello 1997, Numerical Analysis of Fundamental Frequencies, NAFF). To perform the NAFF, we consider the time series $$\begin{array}{cc}{\mathcal{S}}_{\lambda ,i}=\{{a}_{i}({t}_{j})\mathrm{exp}\text{i}{\lambda}_{j}({t}_{j})\},& {\mathcal{S}}_{\mathcal{M},i}=\{{a}_{i}({t}_{j})\mathrm{exp}\text{i}{\mathcal{M}}_{j}({t}_{j})\},\end{array}$$
where j = 0,1,…, t_{j} = jΔt, and a_{i}(t), λ_{i}(t) and M_{i} are the canonical osculating semimajor axis, the mean longitude, and the mean anomaly for planet i, i = 1, 2, 3, 4, respectively, and Δ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 semimajor 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 Nbody 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 Nbody initial conditions, and the canonical elements are more suitable for the qualitative characterization of the orbits.
Considering a particular fourbody MMR chain that can explain astrometric observations of the HR 8799 system, the zerothorder meanmotion resonance implies one of the possible linear combinations of the fundamental frequencies (the proper meanmotions), $$\Delta f(k)\equiv {\displaystyle \sum _{i=1}^{N}{k}_{i}{f}_{i},}$$(2)
where k = [k_{1}, …, k_{N}] is a vector of nonzero integers, and ${\sum}_{i=1}^{N}\left{k}_{i}\right\ne 0$, which yields small Δf. Simultaneously, we may determine the critical angle corresponding to η = arg min Δf(n): $${\theta}_{\eta}(t)\equiv {\displaystyle \sum _{i=1}^{N}{\eta}_{i}{\lambda}_{i},}$$(3)
if Δf(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(k)={n}_{1}2{n}_{2}+{n}_{3}2{n}_{4}.$$(4)
Similarly to the classic Laplace resonance in the Galilean moons system, we consider the generalized MMR as a chain of twobody MMRs, $$\begin{array}{ccc}2{n}_{i+1}{n}_{i}+{g}_{j}\simeq 2{n}_{i+1}{n}_{i}\simeq 0,& i=1,2,3,& j=i,i+1\end{array},$$
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 {\displaystyle \sum _{i=0}^{3}\left(\frac{{n}_{i}}{{n}_{i+1}}2\right)}=0,$$(5)
given the proper meanmotions resolved from the time series S_{M,i}. Numerically, we find that the exact resonance yields Δn ≃ 10^{−13} rad d^{−1} and Δf ≃ 10^{−15} rad d^{−1} for 3 × 8192 samples of Δ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 σ_{Δn} to be determined later.
We determined the frequency variance σ_{Δn} based on extensive numerical experiments that relied on calibrating the Δn range with the Lagrange (geometric) stability interval, which should not be shorter than the stellar age of 40–0 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 Δn prior is crucial to bound (resolve) a longterm stable system in a very short integration time that spans merely a few hundred of the outermost orbits, ≃2 × 10^{5} 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 longterm stable if Δn is roughly less than 10^{−5}–10^{−6} rad d^{−1}, as described above.
4.5 Astrometric fits to nearresonance 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 nearresonant Nbody solutions, besides uniform priors for orbital parameters in sufficiently wide ranges. In this sense, the assumed prior set is minimal. Moreover, the Nbody model is parameterized with (m_{i}, a_{i}, x_{i} ≡ e_{i} cos ϖ_{i}, y_{i} = e_{i} sin ϖ_{i}, M_{i}), i = 1, …, 4, that is, mass, semimajor 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, Ω) determining the orbital plane of the system, and the astrometric parallax Π. 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 256000 iterations, we could obtain an acceptance ratio as low as 0.1. Therefore, we understand the MCMC sampling with the Δn prior as dynamically constrained optimization (rejection sampling) which helps to determine parameter ranges for the bestfitting, 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 nearresonance models consistent with the astrometric observations.
We performed many MCMC experiments varying cγ_{δb} within the range of 10^{−7}–10^{−9} rad d^{−1} and parameters of the emcee samplers. The example posterior for σ_{Δn} = 10^{−8} rad d^{−1} 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 semiamplitude of the critical angles of the Laplace resonance for the integration time spanning ≃200 outermost periods. Fig. 3 shows a collection of samples with Δn < 10^{−5}rad d^{−1}, and the semiamplitude of Laplace resonance argument Δθ < 60°. 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 quasiperiodic 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 Δ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 ℒ or χ^{2} may be somewhat improved with the quasiperiodic model, the RMS of the bestfitting models remains almost the same in the two approaches. It may be concluded that the quasiperiodic model does not lead to any qualitative change in the dynamical status of the system, apart from making it possible to systematically explore longterm stable solutions, which are not necessarily regular. We discuss this further in Sect. 4.6.
Figure 4 illustrates the bestfitting orbital solutions (red curves) overplotted for models within 1σ (grey curves) randomly selected from the MCMC samples such that their RMS varies between 7 and 9 mas. The bestfitting 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 ( topmiddle panel) and HR8799d (bottomright 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, shortterm 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 topright 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 bestfitting model in Table 6 in the ΔRA–ΔDecplane 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 bestfitting models marked in Fig. 4. We selected MCMC samples in Fig. 3 for m_{⋆} = 1.47 M_{⊙} and from other MCMC sampling for m_{⋆} = 1.52 M_{⊙}. 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 longterm 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 twobody 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 Δn ≃ 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 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 〈Y〉 is a numerical technique designed to efficiently characterize the stability of the Nbody 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 CPUparallelized µ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 Δn ≃ 10^{−5} 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 semiamplitudes of other critical angles are similar to those for quasiperiodic 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 longterm stable models characterised by Δn ≃ 10^{−6}rad d^{−1}. All these models selfdestruct 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 ≃10 M_{Jup} range, and Model 4 is especially interesting given its low χ^{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 HR8799e is consistent with the dynamical estimate in Brandt et al. (2021).
These examples are to justify that the system may be dynamically longterm stable in the planet mass range of ≃10 M_{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 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, quasiperiodic, 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 χ^{2} and RMS. The models are difficult to distinguish statistically and visually from each other. This may mean that for a coplanar, nearresonance, 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.
Fig. 3 Posterior samples for stable solutions characterized with Δn < 10^{−5} rad d^{−1}, and the semiamplitude of the libration of the Laplace resonance argument Δθ < 60°. 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. 
Fig. 4 Bestfitting solutions to the fourplanet model in Table 6 for m_{⋆} = 1.47 M_{⊙}, illustrated on the sky plane as randomly selected MCMC samples from the Δn posterior. The yaxis corresponds to the north (N), and the xaxis corresponds to the east (E) direction, respectively (we note that the numerical values of are opposite in sign with respect to the formal lefthand direction of the RA α). Redfilled circles are for the new or rereduced IRDIS measurements in this paper, yellow hexagons are for the IFS measurements, green hexagons are for the LUCI points, lightblue filled circles are for measurements in Konopacky et al. (2016), and darkblue 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 ≃ 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 ≃ 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 closeups and interesting regions. 
Fig. 5 Residuals to the selected bestfitting, nearresonant fourplanet Model 2 in Table 7 for planets HR 8799b and HR 8799c; stellar mass is m_{⋆} = 1.47 M_{⊙}. The yaxis corresponds to the north (N), and the xaxis corresponds to the east (E) direction, respectively (we note that the numerical values of ΔRA are signopposite to regarding the formal lefthand 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 darkblue and lightblue (lightgrey) 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. 
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 socalled 〈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 semimajor 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 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_{Jup}.
We considered two types of probe objects: Cereslike asteroids with a mass of 10^{−6} M_{Jup} and Jupiterlike planets with a mass of 1 M_{Jup}. To calculate the 〈Y〉 values for the synthetic systems, we integrated the Nbody equations of motion and their variational equations with the GraggBulirsch–Stoer integrator (Hairer et al. 2000) for 3 Myr, which covers 6 × 10^{4} orbital periods of HR8799e and 7 × 10^{3} orbital periods of HR 8799b, respectively. That integration time is consistent with the typical 10^{4} outermost orbits required to achieve 〈Y〉 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 Cereslike asteroids are illustrated in Fig. 10 and Fig. 11. In this experiment, we sampled the semimajor axis a_{0} ∈ [4, 18] au and eccentricity e_{0} ∈ [0, 0.8] of these objects. We collected 3 × 10^{5} Nbody initial conditions with 〈Y〉 − 2 < 0.05 that represent the structure of longtermstable 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 Cereslike objects defined at the epoch 1998.830 up to the given final epoch t. For instance, t = 2009.575 (middleright 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 overplotted 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 HR8799e 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 HR8799e; 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 orbital phases, (a_{0}, λ_{e} − λ_{0}) (bottom panel). In both graphs, we mark the test particles with the same color scheme as in the previous plot, and the loworder 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 Δa_{0} = 4 au towards the star and clearly marks a boundary of stable particles.
The bottom graph shows objects involved in loworder MMRs with HR 8799e that are apparently widely spread on the sky plane, but they appear in narrow and welllocalized islands on the planes of orbital elements; this is particularly clear on the (a_{0}, λ_{e} − λ_{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 loweccentricity 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 counterintuitive. 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 longtermstable but mutually noninteracting asteroids in the system, on the same orbit as the innermost planet, whose orbital evolution is governed by the gravitational interactions with the massive Jovian companions. Such objects should not be significantly affected by nongravitational forces, such as PoyntingRobertson 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 indepth analysis of the dynamical structure of the warm dust disk to a future paper.
Fig. 6 Residuals to the bestfitting, nearresonant fourplanet Model 2 in Table 7 for m_{⋆} = 1.47 M_{⊙}, a continuation of Fig. 8 for planets HR8799d and HR 8799e. See the caption of Fig. 5 for labels. The yellow star marks the GRAVITY measurement for HR 8799e. 
Fig. 7 Orbital evolution of selected longtermstable 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 topleft corner of each column. 
Fig. 8 Orbital evolution of selected longtermstable 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. 
4.8 Putative fifth Jupiterlike 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 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 lowmass 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 Δ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 loworder 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 ≃7 au, they should be involved in a 2:1e, 3:1e, or 5:2e MMR with the inner planet, respectively. These loworder 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_{Jup}, respectively. Neither did these authors detect any new candidate companions at the smallest observable separation, of namely 0.1″ or ≃4.1 au, overlapping with the semimajor axes range in Figs. 1012. Wahhaj et al. (2021) conclude that the planet may still exist with a mass of 23.6 MJup at 7.5 au (3:1e MMR) or 1.5–2.8 M_{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 ≃1 M_{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 ≃10 au under the current observational configuration of four planets (Fig. 11), and related to the abovementioned 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.
Fig. 9 Dynamical maps in the Keplerian, osculating astrocentric semimajoraxiseccentricity plane for the quasiperiodic, nearresonant configuration described by Model 2. All masses and elements but the map coordinates are fixed at their bestfitting values in Table 7. The osculating epoch is 1998.830. Stable solutions are determined with 〈Y〉 − 2 ≃ 0 and marked with blue color. 
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 multipleplanet 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 insideout depletion of gas in the disk due to the presence of the planet itself and photoevaporation lead to different migration speeds for the planets, which may end up in resonance. Numerical modeling by Hands et al. (2014) and Szuszkiewicz & PodlewskaGaca (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 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_{⊙}. 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 Σ = Σ_{0}r^{−1/2} with Σ_{0} = 50 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 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 semimajor 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 Σ_{0} to 100 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 semimajor 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 farther out. It is noteworthy that the wider oscillations observed in the critical arguments of the resonances and the semimajor axis of the planets compared to a pure Nbody problem are related to the presence of the perturbations of the disk on the planets.
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_{Jup} at the epoch labeled in the lowerleft corner of each graph. We collected 3.3 × 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 quasihomogeneous disk (gray points). The 2024.83 epoch covers roughly one orbital period in a 2:1e MMR resonance (dark blue circles) with the innermost planet. 
Fig. 11 Structure of the inner debris disk composed of 3 × 10^{5} Cereslike asteroids of mass 10^{−6} M_{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, quasihomogeneous 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 Δa_{0} ≃ 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. 
6 Planetsdisk 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 Kuiperlike 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 dustfree 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 Nbody 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 semimajor 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_{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 Nbody 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 planetdisk 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 Nbody 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_{Jup}. These simulations are based on a quasiperiodic, nearresonant orbital model 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_{⋆} = 1.52M_{⊙} and a larger parallax Π = 25.4 mas, the whole system appears more compact than predicted at present – these authors found the osculating astrocentric semimajor axis of HR 8799b ≃67.1 au, compared to ≃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 semimajor axis >90 au. These planets could be involved in loworder 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.
Fig. 12 Possible instant positions (small gray and colored dots) of a hypothetical and still undetected planet of mass 1 M_{Jup} on the sky plane, at the epoch labeled in the lowerleft 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 lightgray curve marks the collision zone with HR 8799e, the thinner curve marks the stability zone for m_{0} = 10^{−6} M_{Jup} (see Fig. 10), and the red curve is the image of the geometric collision curve shifted by Δa_{0} ≃ 6 au toward the star. 
Fig. 13 HR 8799 architecture: two belts (blue lines), four detected planets (pink circles), and a putative fifth planet (lightblue 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. 
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 stateoftheart algorithms that apply the ADI technique, and for consistency, we also repeated the reduction of published SPHERE data from opentime 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 rederived 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 longterm stability of the orbits (over a few hundred 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 ≃6 M_{Jup}. Moreover, the dynamical parallax is consistent with 1σ uncertainty with the corrected, independently determined Gaia eDR3 value of 24.50 ± 0.05 mas, reinforcing both the adopted mass of the host star m_{⋆} = 1.47 M_{⊙} and the assumed resonant or closetoresonant configuration of the system.
Regarding the quality of the Nbody astrometric models considered in this work, we did not find any significant or qualitative improvement of multiparameter nearresonant configurations 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 nearresonant, 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 welldefined PO Nbody 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 quasiexact resonance.
We looked for the dynamical space where a putative fifth planet could sit using analytical formulae and estimates, and suggest that only a lowmass 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 highcontrast 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 Nbody simulations. The results may be useful for predicting the positions of smallmass objects below the present detection limits of ≃3 M_{Jup}.
Fig. 14 Resonance trapping of the HR 8799 system during the migration phase. The upper panel shows the evolution of the semimajor axis of the outer planet at the moment of capture in resonance with the inner three bodies. After the resonance trapping, the semimajor 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. 
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 semimajor axes of the planets rescaled with constant values for ease of visualization in the same plot. 
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 Poznan Supercomputer and Network Centre (PCSS, Poland) for the generous longterm support and computing resources (grant No. 529). K.G. also thanks the staff of the Tricity Supercomputer Centre, Gdansk (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 RII3Ct2004001566 for FP6 (2004–2008), grant number 226604 for FP7 (2009–2012) and grant number 312430 for FP7 (2013–2016).
Appendix A Astrometric points used in the analysis
List of the astrometric points.
References
 Apai, D., Kasper, M., Skemer, A., et al. 2016, ApJ, 820, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Baines, E. K., White, R. J., Huber, D., et al. 2012, ApJ, 761, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Baraffe, I., Chabrier, G., Barman, T. S., Allard, F., & Hauschildt, P. H. 2003, A&A, 402, 701 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beaugé, C., Michtchenko, T. A., & FerrazMello, S. 2006, MNRAS, 365, 1160 [Google Scholar]
 Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593 [Google Scholar]
 Bergfors, C., Brandner, W., Janson, M., Köhler, R., & Henning, T. 2011, A&A, 528, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Biller, B. A., Apai, D., Bonnefoy, M., et al. 2021, MNRAS, 503, 743 [NASA ADS] [CrossRef] [Google Scholar]
 Boccaletti, A., Abe, L., Baudrand, J., et al. 2008, SPIE Conf. Ser., 7015, 70151B [NASA ADS] [Google Scholar]
 Bohn, A. J., Kenworthy, M. A., Ginski, C., et al. 2020, ApJ, 898, L16 [Google Scholar]
 Bohn, A. J., Ginski, C., Kenworthy, M. A., et al. 2021, A&A, 648, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Booth, M., Jordan, A., Casassus, S., et al. 2016, MNRAS, 476, 2845 [Google Scholar]
 Boss, A. P. 1997, Science, 276, 1836 [Google Scholar]
 Brandt, G. M., Brandt, T. D., Dupuy, T. J., Michalik, D., & Marleau, G.D. 2021, ApJ, 915, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, A. G. W. 1978, Moon and Planets, 18, 5 [Google Scholar]
 Chauvin, G., Lagrange, A.M., Dumas, C., et al. 2004, A&A, 425, L29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chen, C. H., Mittal, T., Kuchner, M., et al. 2014, ApJ, 211, 25 [Google Scholar]
 Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Physica D Nonlinear Phenomena, 182, 151 [NASA ADS] [CrossRef] [Google Scholar]
 Claudi, R. U., Turatto, M., Gratton, R. G., et al. 2008, SPIE Conf. Ser., 7014, 70143E [Google Scholar]
 Crida, A., Sándor, Z., & Kley, W. 2008, A&A, 483, 325 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Currie, T., Burrows, A., Itoh, Y., et al. 2011, ApJ, 729, 128 [Google Scholar]
 Currie, T., Fukagawa, M., Thalmann, C., Matsumura, S., & Plavchan, P. 2012, ApJ, 755, L34 [NASA ADS] [CrossRef] [Google Scholar]
 Currie, T., Burrows, A., Girard, J. H., et al. 2014, ApJ, 795, 133 [Google Scholar]
 D’Angelo, G., & Marzari, F. 2012, ApJ, 757, 50 [CrossRef] [Google Scholar]
 De Rosa, R. J., Esposito, T. M., Gibbs, A., et al. 2020, SPIE Conf. Ser., 11447, 114475A [NASA ADS] [Google Scholar]
 Delisle, J.B. 2017, A&A, 605, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delorme, P., Meunier, N., Albert, D., et al. 2017, in SF2A2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. C. Reylé, P. Di Matteo, F. Herpin, et al. [Google Scholar]
 Desidera, S., Chauvin, G., Bonavita, M., et al. 2021, A&A, 651, A70 [EDP Sciences] [Google Scholar]
 Dohlen, K., Langlois, M., Saisse, M., et al. 2008, SPIE Conf. Ser., 7014, 70143L [Google Scholar]
 Esposito, S., Riccardi, A., Fini, L., et al. 2010, SPIE Conf. Ser., 7736, 107 [Google Scholar]
 Esposito, S., Mesa, D., Skemer, A., et al. 2013, A&A, 549, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabrycky, D. C., & MurrayClay, R. A. 2010, ApJ, 710, 1408 [NASA ADS] [CrossRef] [Google Scholar]
 Faramaz, V., Marino, S., Booth, M., et al. 2021, AJ, 161, 271 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Fukagawa, M., Itoh, Y., Tamura, M., et al. 2009, ApJ, 696, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Fusco, T., Rousset, G., Sauvage, J.F., et al. 2006, Opt. Express, 14, 7515 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration 2020, VizieR Online Data Catalog: I/350 [Google Scholar]
 Galicher, R., Marois, C., Macintosh, B., Barman, T., & Konopacky, Q. 2011, ApJ, 739, L41 [Google Scholar]
 Galicher, R., Boccaletti, A., Mesa, D., et al. 2018, A&A, 615, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goździewski, K., & Migaszewski, C. 2018, ApJS, 238, 6 [Google Scholar]
 Goździewski, K., & Migaszewski, C. 2020, ApJ, 902, L40 [CrossRef] [Google Scholar]
 Goździewski, K., Bois, E., Maciejewski, A. J., & KiselevaEggleton, L. 2001, A&A, 378, 569 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 GRAVITY Collaboration (Lacour, S., et al.) 2019, A&A, 623, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gray, R. O., & Kaye, A. B. 1999, AJ, 118, 2993 [NASA ADS] [CrossRef] [Google Scholar]
 Hadjidemetriou, J. D. 1977, Celest. Mech., 16, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
 Hairer, E., Nørsett, S., & Wanner, G. 2000, Solving Ordinary Differential Equations I Nonstiff Problems, 2nd edn. (Berlin: Springer) [Google Scholar]
 Hands, T. O., Alexander, R. D., & Dehnen, W. 2014, MNRAS, 445, 749 [NASA ADS] [CrossRef] [Google Scholar]
 Hinz, P. M., Rodigas, T. J., Kenworthy, M. A., et al. 2010, ApJ, 716, 417 [NASA ADS] [CrossRef] [Google Scholar]
 Hughes, A. M., Wilner, D. J., Andrews, S. M., et al. 2011, ApJ, 740, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kervella, P., Arenou, F., & Thévenin, F. 2022, A&A, 657, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Konopacky, Q. M., Marois, C., Macintosh, B. A., et al. 2016, AJ, 152, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Lafrenière, D., Marois, C., Doyon, R., & Barman, T. 2009, ApJ, 694, L148 [CrossRef] [Google Scholar]
 Lagrange, A. M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57 [Google Scholar]
 Langlois, M., Vigan, A., Moutou, C., et al. 2013, in Proceedings of the Third AO4ELT Conference, ed. S. Esposito & L. Fini, 63 [Google Scholar]
 Langlois, M., Gratton, R., Lagrange, A. M., et al. 2021, A&A, 651, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Laskar, J. 1993, Celest. Mech. Dyn. Astron., 56, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Lazzoni, C., Desidera, S., Marzari, F., et al. 2018, A&A, 611, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596 [Google Scholar]
 Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4 [EDP Sciences] [Google Scholar]
 Maire, A. L., Skemer, A. J., Hinz, P. M., et al. 2015, A&A, 579, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maire, A.L., Langlois, M., Dohlen, K., et al. 2016, SPIE Conf. Ser., 9908, 990834 [Google Scholar]
 Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556 [Google Scholar]
 Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348 [Google Scholar]
 Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 2010, Nature, 468, 1080 [NASA ADS] [CrossRef] [Google Scholar]
 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 [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masset, F., & Snellgrove, M. 2001, MNRAS, 320, L55 [Google Scholar]
 Matthews, B., Kennedy, G., Sibthorpe, B., et al. 2014, ApJ, 780, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Mizuno, H. 1980, Progr. Theor. Phys., 64, 544 [CrossRef] [Google Scholar]
 Moorhead, A. V., & Adams, F. C. 2005, Icarus, 178, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Nesvorny, D., & FerrazMello, S. 1997, A&A, 320, 672 [NASA ADS] [Google Scholar]
 Nowak, M., Lacour, S., Lagrange, A. M., et al. 2020, A&A, 642, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Papaloizou, J. C. B. 2015, Int. J. Astrobiol., 14, 291 [Google Scholar]
 Petit, C., Sauvage, J.F., Fusco, T., et al. 2014, SPIE Conf. Ser., 9148, 0 [Google Scholar]
 Pinna, E., Rossi, F., Puglisi, A., et al. 2019, in Adaptive Optics for Extremely Large Telescopes 6 – Conference Proceedings [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Rameau, J., Chauvin, G., Lagrange, A.M., et al. 2013, ApJ, 772, L15 [Google Scholar]
 Ramos, X. S., Charalambous, C., BenítezLlambay, P., & Beaugé, C. 2017, A&A, 602, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424 [Google Scholar]
 Seifert, W., Ageorges, N., Lehmitz, M., et al. 2010, in Proc. SPIE, Vol. 7735, Groundbased and Airborne Instrumentation for Astronomy III, 77357W [Google Scholar]
 Sepulveda, A. G., & Bowler, B. P. 2022, AJ, 163, 52 [NASA ADS] [CrossRef] [Google Scholar]
 Soummer, R., Hagan, J. B., Pueyo, L., et al. 2011, ApJ, 741, 55 [Google Scholar]
 Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28 [Google Scholar]
 Su, K. Y. L., Rieke, G. H., Stapelfeldt, K. R., et al. 2009, ApJ, 705, 314 [Google Scholar]
 Sudol, J. J., & Haghighipour, N. 2012, ApJ, 755, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Szuszkiewicz, E., & PodlewskaGaca, E. 2012, Origins Life Evol. Biosphere, 42, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Thommes, E. W. 2005, ApJ, 626, 1033 [NASA ADS] [CrossRef] [Google Scholar]
 Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72 [EDP Sciences] [Google Scholar]
 Wahhaj, Z., Milli, J., Romero, C., et al. 2021, A&A, 648, A26 [EDP Sciences] [Google Scholar]
 Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, AJ, 156, 192 [Google Scholar]
 Wertz, O., Absil, O., Gómez González, C. A., et al. 2017, A&A, 598, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wilner, D. J., MacGregor, M. A., Andrews, S. M., et al. 2018, ApJ, 855, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, J. T., Veras, D., Ford, E. B., et al. 2011, ApJ, 730, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Zuckerman, B., Rhee, J. H., Song, I., & Bessell, M. S. 2011, ApJ, 732, 61 [Google Scholar]
 Zurlo, A., Vigan, A., Mesa, D., et al. 2014, A&A, 572, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zurlo, A., Vigan, A., Galicher, R., et al. 2016, A&A, 587, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
Values for the mass of each planet calculated using the evolutionary models in the different IRDIS filters.
Osculating, heliocentric elements of the bestfitting solutions at the epoch of 1998.830.
All Figures
Fig. 1 MCMC posterior for the bestfitting, strictly resonant model derived for m_{⋆} = 1.47 M_{Sun}. Parameters included in the diagram are for the planet masses and orbital period ratio κ = P_{d}/P_{e} = κ, and parallax Π, inclination I, and the nodal angle Ω (the orbital scale ρ, PO epoch shift t_{0} and the common rotation angle ω_{rot} in the orbital planet are skipped). The crossed lines mark the bestfitting PO solution in terms of the smallest χ^{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. 

In the text 
Fig. 2 Temporal evolution of the osculating orbital elements in the astrocentric (blue curves) and Jacobian (yellow curves) frames, respectively, for the same Nbody 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. 

In the text 
Fig. 3 Posterior samples for stable solutions characterized with Δn < 10^{−5} rad d^{−1}, and the semiamplitude of the libration of the Laplace resonance argument Δθ < 60°. 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. 

In the text 
Fig. 4 Bestfitting solutions to the fourplanet model in Table 6 for m_{⋆} = 1.47 M_{⊙}, illustrated on the sky plane as randomly selected MCMC samples from the Δn posterior. The yaxis corresponds to the north (N), and the xaxis corresponds to the east (E) direction, respectively (we note that the numerical values of are opposite in sign with respect to the formal lefthand direction of the RA α). Redfilled circles are for the new or rereduced IRDIS measurements in this paper, yellow hexagons are for the IFS measurements, green hexagons are for the LUCI points, lightblue filled circles are for measurements in Konopacky et al. (2016), and darkblue 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 ≃ 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 ≃ 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 closeups and interesting regions. 

In the text 
Fig. 5 Residuals to the selected bestfitting, nearresonant fourplanet Model 2 in Table 7 for planets HR 8799b and HR 8799c; stellar mass is m_{⋆} = 1.47 M_{⊙}. The yaxis corresponds to the north (N), and the xaxis corresponds to the east (E) direction, respectively (we note that the numerical values of ΔRA are signopposite to regarding the formal lefthand 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 darkblue and lightblue (lightgrey) 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. 

In the text 
Fig. 6 Residuals to the bestfitting, nearresonant fourplanet Model 2 in Table 7 for m_{⋆} = 1.47 M_{⊙}, a continuation of Fig. 8 for planets HR8799d and HR 8799e. See the caption of Fig. 5 for labels. The yellow star marks the GRAVITY measurement for HR 8799e. 

In the text 
Fig. 7 Orbital evolution of selected longtermstable 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 topleft corner of each column. 

In the text 
Fig. 8 Orbital evolution of selected longtermstable 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. 

In the text 
Fig. 9 Dynamical maps in the Keplerian, osculating astrocentric semimajoraxiseccentricity plane for the quasiperiodic, nearresonant configuration described by Model 2. All masses and elements but the map coordinates are fixed at their bestfitting values in Table 7. The osculating epoch is 1998.830. Stable solutions are determined with 〈Y〉 − 2 ≃ 0 and marked with blue color. 

In the text 
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_{Jup} at the epoch labeled in the lowerleft corner of each graph. We collected 3.3 × 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 quasihomogeneous disk (gray points). The 2024.83 epoch covers roughly one orbital period in a 2:1e MMR resonance (dark blue circles) with the innermost planet. 

In the text 
Fig. 11 Structure of the inner debris disk composed of 3 × 10^{5} Cereslike asteroids of mass 10^{−6} M_{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, quasihomogeneous 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 Δa_{0} ≃ 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. 

In the text 
Fig. 12 Possible instant positions (small gray and colored dots) of a hypothetical and still undetected planet of mass 1 M_{Jup} on the sky plane, at the epoch labeled in the lowerleft 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 lightgray curve marks the collision zone with HR 8799e, the thinner curve marks the stability zone for m_{0} = 10^{−6} M_{Jup} (see Fig. 10), and the red curve is the image of the geometric collision curve shifted by Δa_{0} ≃ 6 au toward the star. 

In the text 
Fig. 13 HR 8799 architecture: two belts (blue lines), four detected planets (pink circles), and a putative fifth planet (lightblue 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. 

In the text 
Fig. 14 Resonance trapping of the HR 8799 system during the migration phase. The upper panel shows the evolution of the semimajor axis of the outer planet at the moment of capture in resonance with the inner three bodies. After the resonance trapping, the semimajor 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. 

In the text 
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 semimajor axes of the planets rescaled with constant values for ease of visualization in the same plot. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.