EDP Sciences
Free Access
Issue
A&A
Volume 592, August 2016
Article Number A147
Number of page(s) 14
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201526309
Published online 18 August 2016

© ESO, 2016

1. Introduction

HR 8799 (also called HD 218396 and HIP 114189) is a nearby (39.4 ± 1.1 pc; van Leeuwen 2007), young ( Myr assuming member of the Columba association; Marois et al. 2010) A5V (up to F0V) star with mass 1.5 ± 0.3M (Gray & Kaye 1999; Gray et al. 2003). HR 8799 has four confirmed, directly imaged planets (Marois et al. 2008, 2010) and a debris disk (Su et al. 2009; Matthews et al. 2014). The planet masses are estimated to be 5 (b) and 7 (cde) MJup (Marois et al. 2010). The on-sky separations between planets and star are estimated to be 67.9 (b), 38.0 (c), 24.5 (d), and 14.5 (e) AU. Matthews et al. (2014) estimated the inclination of the system with respect to the plane of the sky to be 26 ± 3° by observing the outer part of the debris disk and assuming coplanarity with the planetary orbits. They estimate the position angle of the inclination to be 64 ± 3° (G. Kennedy, private communication). Figure 1 shows all observations of the planets published at the time of writing.

The estimates of the planetary masses in the HR 8799 system depend on the age of the system through cooling models, indicating higher planet masses with higher age (Baraffe et al. 2003). The age of the HR 8799 system is uncertain and has been estimated using different techniques to range between 30 Myr and 1 Gyr (see Marois et al. 2010, 2008, and references therein). The planet mass estimates reach the brown dwarf regime when an age of ≳ 100 Myr is assumed. In this paper we consider HR 8799 bcde as planets; the case of brown dwarfs is more closely investigated by Moro-Martín et al. (2010). We assume the age 30 Myr in this paper as it is estimated from the Columba association. An age of 30 Myr leads to multi-Jovian masses that, however, are still in the planetary regime. Whether it is true that HR 8799 is associated with the Columba association could be revealed by Gaia in the near future (Perryman et al. 2001). Other uncertainties in system properties include distance (39.4 ± 1.1 pc, van Leeuwen 2007) and stellar mass (1.5 ± 0.3 M, Gray & Kaye 1999; Gray et al. 2003).

The formation of the massive HR 8799 bcde planets can be explained by the recently developed models of pebble accretion (Lambrechts & Johansen 2012), but the dynamics of the system has been a puzzle unless resonances are considered. The orbits of the HR 8799 bcde planets are difficult to constrain as the planets have been observed during a small fraction of their orbits (see Fig. 1). In recent work by Pueyo et al. (2015), however, fits are made to constrain the orbital elements and a variety of solutions are found, favouring a slightly more eccentric orbit for planet d. It is reasonable to believe that the planet orbits are close to circular, as with high eccentricity the orbits are likely to cross. However, the HR 8799 planetary system appears to be unstable on a timescale much shorter than the estimated age of the system when integrating the orbits of the planets, assumed initially circular, seen pole-on and non-resonant (Fabrycky & Murray-Clay 2010; Reidemeister et al. 2009; Sudol & Haghighipour 2012). A possible and popular explanation to the survival of the HR 8799 planetary system is that the planets are locked deep in strong mean motion resonances that can lead to stability timescales up to 1 Gyr (e.g. Fabrycky & Murray-Clay 2010; Goździewski & Migaszewski 2009; Reidemeister et al. 2009; Moro-Martín et al. 2010; Currie et al. 2011; Soummer et al. 2011; Esposito et al. 2013; Goździewski & Migaszewski 2014). These resonant configurations are characterised by small-amplitude librations of one or more resonant arguments corresponding to two-, three-, or four-planet mean motion resonances, thus ensuring that systems remain protected for extended periods from close encounters between the planets. Moore & Quillen (2013) find stabilising effects from the outer debris disk, although this requires a very massive debris disk.

thumbnail Fig. 1

Observations of the HR 8799 planets. The black dashed line corresponds to the position angle of the inclination as estimated by Matthews et al. (2014) with the ± 1σ uncertainty shown by the grey region.

Open with DEXTER

In this paper we show simulations of planetary systems that look like the HR 8799 system, are stable for longer than the estimated age, and are not stabilised by strong resonant lock. We create these systems by starting the planetary orbits slightly more widely separated than the observed on-sky separation assuming circular orbits. A more widely separated configuration can also look like the HR 8799 system as the orbits are not exactly circular. Larger separation between the planetary orbits results in a much longer stability timescale, which makes systems remain stable for longer than the estimated age of HR 8799. This solution is different from previously published explanations, as we do not need protective resonant lock or any other stabilising effects to create a long-lived system.

In Sect. 2 we describe how we simulate planetary systems. Sections 3 and 4 present first results consistent with earlier work and later long-lived solutions in agreement with the observed system. In Sect. 5 we introduce a concept of evolutionary phases in the lives of wide-orbit planetary systems. These phases could be used to better understand the past and future of a wide-orbit planetary system. We summarise our findings and conclusions in Sect. 6.

2. Simulations

Table 1

Initial conditions of the simulations.

2.1. Orbital parameters

All planets in all our simulations are analogues of the HR 8799 bcde planets and they are initiated on circular, heliocentric orbits. The initial orbit inclination we find by drawing a random inclination between 0° and 5° and random longitude of ascending node between 0 and 360° (β = 5°; see Johansen et al. 2012). The initial planet position on the orbit (true anomaly) is randomised. We assign the mass 1.5M to the star and the planets are given the estimated values of the masses, 5MJup (b) and 7MJup (cde) if nothing else is stated (the masses are decreased in simulations 5 and 6). The planets orbit the barycentre of the system, which is not centred on the star as a result of the high planet mass. This means that the planetary orbits have a low initial eccentricity as they are initiated on orbits circular around the star. The initial orbital eccentricities vary between different systems as the inclinations, longitude of ascending nodes, and true anomalies are randomised.

We vary the separation between the initial planet orbits (regarding the planets as test particles) in terms of the mutual Hill radius, which is defined as follows: (1)where r signifies star-planet separation in the orbital plane and M mass (Chambers et al. 1996). The indices p1 and p2 refer to the two planets and refers to the star. We then write the separation between the initial, circular orbits in mutual Hill radii as (2)where in this case p2 refers to the outer planet and p1 to the inner planet. For a given planetary system there is a critical Δ at which in 50% of the cases the system leads to orbit crossing within a given time (Chambers et al. 1996; Chatterjee et al. 2008), although for any Δ there is considerable scatter in the times in which particular realisations become unstable.

2.2. MERCURY6 parameters

We perform our simulations using the MERCURY6 package (Chambers 1999). We use the Bulirsch-Stoer integrator and a Jacobian coordinate system in the output files. Ejection distance is 11 000 AU and fragmentation from collisions is not allowed. We have redefined a close encounter as when two planets are separated by less than one mutual Hill radius, defined as in Eq. (1).

2.3. Simulation parameters

In all simulations apart from 4f_long we stop the integration when a close encounter has occurred or 100 Myr is reached. In 4f_long we continue running until 100 Myr to investigate the future of multi-planet systems in Sect. 5. We use output rates of 10, 100 or 1000 yrs. Every simulation contains 100 runs with the difference in initial placement on orbits and the inclinations of the orbits. The simulations we carried out are summarised in Table 1 and described in detail in Sects. 3 and 4. We provide the initial conditions of our example system in Table A.1, but we note that the system is so chaotic that a small change in time-step completely alters the lifetime of the system (see Sect. 4.3). A more successful method is to create statistically similar, long-lived systems as explained in this section and described in Sect. 4 (see also Table 2).

Throughout the paper we compare our simulations to the observations of Currie et al. (2011) from 8 October 2009, as these observations have low measurement uncertainties and all four planets are observed simultaneously.

3. HR 8799 – unstable as observed

Below we describe in more detail the three simulations 1, 2, and 3 seen in Table 1. We use the term stability timescale as the timescale for which the systems in a simulation have not had any close encounters. This timescale we find is roughly the same as the time it takes for the systems to lose a planet.

In Fig. 2 the number of systems with 4, 3, 2, and 1 planets in simulations 1, 2, and 3 is shown as a function of time. The black line indicates the estimated age of HR 8799. The figure shows that none of the systems in simulations 1, 2, or 3 still have four planets when the age of HR 8799 is reached. The figure also shows that ~80% of the systems have two planets left and ~20% have one planet left after 100 Myr.

Table 2

Observed wide-orbit planets around stars.

thumbnail Fig. 2

Number of planets left in the systems as a function of time. Dotted line corresponds to the face-on assumption (simulation 1), solid line to the inclination inferred from observations (simulation 2), and dashed line to the largest separation achieved from inclination (simulation 3). See Table 1 and Sect. 3 for descriptions about the simulations. The age estimate of HR 8799 ( Myr) is indicated with a solid line with the grey region corresponding to the uncertainty.

Open with DEXTER

3.1. Simulation 1: face-on assumption

For simplicity, we begin by assuming the HR 8799 system is seen entirely face-on, i.e. inclination i = 0° compared to the plane of the sky. We use the observed on-sky star-planet separations of 67.9 (b), 38.0 (c), 24.5 (d), and 14.5 (e) AU as initial star-planet separations. The initial separations between the orbits are then Δ = 4.14 (b-c), 3.02 (c-d) and 3.55 (d-e) RH. The dynamics of the face-on assumption has been tested earlier by Fabrycky & Murray-Clay (2010) in the three-planet case and Marois et al. (2010) with all four planets. We confirm their results by finding a median time at which the planetary systems have the first close encounter of 0.10 Myr and there is no planetary system without close encounters for longer than 9.4 Myr. We refer to this simulation as simulation 1, see Table 1 and dotted lines in Fig. 2. Although none of our systems live for the age of HR 8799, we point out here the large range of lifetimes with the longest-lived system surviving 100 times longer than the median.

3.2. Simulation 2: inclined system according to observations

thumbnail Fig. 3

Smallest orbit separation, Δmin, between any planet pair changes with inclination and position angle. We indicate the inclination and position angle of the HR 8799 system measured by Matthews et al. (2014) with a white square (simulation 2), which is i = 26 ± 3° and position angle 64 ± 3°. The largest Δmin occurs at i = 25° and position angle 42°, when Δmin = 3.56RH. This value is indicated with a filled, white circle (simulation 3). This plot is based on the observations from Marois et al. (2008, 2010).

Open with DEXTER

The disk of the HR 8799 system is observed to be inclined by i = 26° compared to the plane of the sky and at a position angle of 64° E of N (64 ± 3°, best fit; G. Kennedy, priv. comm.; Matthews et al. 2014). Here we assume that the planets and disk are mutually coplanar. These inclination parameters we use also later in simulations 4, 5, and 6. The observed, on-sky separations between planets and star are then 67.9 (b), 41.9 (c), 25.8 (d), and 14.5 (e) AU, which means that planets c and d are located slightly further away from the star compared to the face-on case (simulation 1). The initial separations between the orbits are then Δ = 3.47 (b-c), 3.31 (c-d), and 3.90 (d-e) RH. We refer to this simulation as simulation 2 (see Table 1, solid lines in Fig. 2 and white square marker in Fig. 3).

The median time during which the systems have not had any close encounter is in this case 0.094 Myr, while there are no systems that still have four planets after 9.39 Myr. Compared to simulation 1 the stability timescale is similar and it remains too short to explain the current state of HR 8799.

3.3. Simulation 3: optimal separation from inclination

Assume we did not have any system inclination estimate. Which inclination and position angle would then give the largest planetary orbit spacing and thus the longest stability timescale? In Fig. 3 the minimum separation between any two planetary orbits in the system is given for sets of inclination and position angle assuming circular orbits centred on the observed mean (using the observations of Marois et al. 2008, 2010). The largest Δmin is found with an inclination of i = 25° and a position angle of 42° . The separation between the orbits are then Δ = 3.56RH for all planet pairs, which translates into star-planet separations of 68.8 (b), 41.9 (c), 24.8 (d), and 14.7 (e) AU. We use these optimal separation parameters in simulation 3 (see Table 1, dashed lines in Fig. 2 and white circle in Fig. 3). The median time the systems have not had any close encounter is 0.57 Myr and all systems have had close encounters after 6.85 Myr. Figure 2 shows that an increase of 0.25RH in Δmin results in a factor of ~2 longer stability timescale (compare medians of simulation 2 with simulation 3).

Preferred estimates of the HR 8799 system inclination from a dynamical point of view lie around 20−30° (Reidemeister et al. 2009; Sudol & Haghighipour 2012), while astroseismological measurements favour inclinations higher than 40° assuming that the planetary orbits are aligned with the stellar spin (Wright et al. 2011).

None of simulations 1, 2 or 3 manage to explain the existence of the HR 8799 system as they simply fall apart before the age of the system is reached.

4. A stable solution without resonant lock

We find long-lived, simulated systems that look like HR 8799 during some time in their evolution. These systems are created with initially wider separation between the orbits than calculated from observations assuming circular orbits. This means that we assume the planets are not moving on completely circular orbits.

Our simulated planetary orbits have a low eccentricity of e ≲ 0.05 from the beginning of the integration (for description on how this is set up see Sect. 2.1). The eccentricity enables planets on more widely separated orbits to look like HR 8799 during parts of the orbits. A wider separation between the orbits in a system can make the system survive longer than the estimated age of HR 8799 as the stability timescale is strongly dependent on orbital separation (Davies et al. 2014).

In this section we present our simulations 4, 5, and 6 (see Table 1), which all contain long-lived systems that look like HR 8799. In Sect. 4.1 we describe simulation 4 and show an example of a long-lived system that seems to fit the observations of HR 8799. In Sect. 4.2 we go into detail on how we determine if a simulated system fits the observations. Section 4.3 presents a way of finding even more well-fitting systems. In Sect. 4.4 we analyse the simulations in terms of resonant behaviour. Section 4.5 describes simulations 5 and 6, which consider lower planet masses that are still within the estimated range.

4.1. Simulation 4: wider initial orbital separation

We test a wider-orbit assumption by simulating systems with the following initial equal separations between the planet orbits: Δ = 3.6, 3.65, 3.7, 3.75, 3.8, 3.85, 3.9, 3.95, and 4 RH. We place planet e at an initial star-planet separation of 14.3 AU, which then fixes the initial semi-major axes of the other planets. The planet masses are set to the nominal masses M = 5 MJup (b) and 7 MJup (cde) (Marois et al. 2010). These simulations are called simulation 4a–i; see Table 1. Simulation 4 thus contains 900 runs compared to simulations 1, 2, and 3, which only contain 100 runs each. It makes sense to talk about simulations 4a–i for the initialisation of the runs, but not in terms of the architecture of the outcome. The measured minimum orbit separation Δmin for each system ranges between 3.6 and 4 RH in almost each of simulations 4a–i, depending on the orbital eccentricities. Therefore we bin the simulations 4 with Δmin, regardless of the initial set-up. We make nine bins with each 100 runs. These bins contain different ranges of orbit separations (Δmin) and thus represent slightly different architectures of the planetary system.

thumbnail Fig. 4

Number of systems without close encounters shown with time. The different coloured lines correspond to bins of each 100 systems in simulation 4, sorted in measured minimum orbit separations (Δmin). The legend indicates the bins and their systems’ mean value of Δmin. The black line and the grey zone show the age estimate with uncertainty of HR 8799.

Open with DEXTER

Figure 4 shows the time during which the systems in the simulation 4 bins have not had any close encounter. The figure shows the trend of stability timescale with minimum orbit separation, Δmin. Compared to the systems in Fig. 2 a substantial fraction of systems in Fig. 4 have not experienced any close encounter at the estimated age of HR 8799 (~30 Myr). Chatterjee et al. (2008) find that a larger number of mutual Hill radii are needed to keep their planetary systems stable compared to what we find. For various reasons (e.g. our higher planet masses, lower initial eccentricities, and different stability criterion), we consider both results valid and our results not contradicting the results of Chatterjee et al. (2008). We see that at any given separation there is a wide range of stability times, a fact to which we return later.

Do the surviving systems in simulation 4 look like HR 8799? In Fig. 5 we compare a system from simulation 2 (left panel, system 8) with a system from simulation 4e (right panel, system 94). Both runs seem to fit the observations, but the system from simulation 4 survives over 45 Myr because of the slightly more widely separated planetary orbits. We want to stress that system 94, simulation 4e is not a special system, but seems normal compared to the other systems in simulation 4. As we explain in Sect. 4.2 we find 12 systems fitting the observations in simulation 4. If we consider a distance to the HR 8799 system 1σ away from the mean of the estimate (39.4 + 1.1 pc = 40.5 pc) we almost double the number to 22 fitting systems in simulation 4. How we determine if a system fits the observations is discussed in detail in Sect. 4.2.

To demonstrate the similarity of our example system, system 94, to the observed HR 8799 system, we show a snapshot of the system at 37.6 Myr in Fig. 6.

thumbnail Fig. 5

Comparison between a system in simulation 2 and a system in simulation 4e. Both systems seem to fit the observed positions of the planets, but only the system in simulation 4e is stable for longer than the estimated age of HR 8799. This is an on-sky projection of the simulations assuming a distance of 39.4 pc (van Leeuwen 2007) and inclination 26° (position angle 64° , Matthews et al. 2014). The axis around which the system is inclined is indicated with a dashed line. The position of the star is indicated with an asterisk.

Open with DEXTER

4.2. Determining whether a system fits the observations

In this section we describe in detail how we determine how each simulated system fits the observations. We have chosen to compare the simulated systems with the observations taken on 8 October 2009 by Currie et al. (2011) because the planets are observed simultaneously during this observation and the measurement uncertainties are comparably small.

We start by checking whether the system has a longer stability timescale than the estimated age of HR 8799 (30 Myr, Marois et al. 2010). If that is the case, we proceed with checking whether the system looks like HR 8799. We check whether the architecture of the simulated system is similar to the architecture of HR 8799 in the following way:

  • 1.

    For every output, move each planet along its orbit to theobserved position angle. This is done because it is highly unlikelythat we find a single output where all four planets happen to belocated within the uncertainties of the observations. The spacethe observations take on compared to the orbit is about0.04 AU/200 AU =2 × 10-4. The chance of then catching a point in time when all four planets are within the observations would be (2 × 10-4)3 = 8 × 10-12 for a very well-fitting system (fixing for example b to always be at the observed position angle). To catch one output where all planets fit simultaneously would therefore require ~ 1011 outputs – a number about three orders of magnitude higher than our output rate. As during the simulation the planets are at some point located at exactly the observed position angle, but we do not catch that moment, it is validated to move the planets along their orbits in this way to really measure how close they come to the observations.

  • 2.

    Draw one million random star-planet separations from the distributions created in step number 1 for each planet. This is carried out to create a new distribution of star-planet separations with a method called bootstrapping with resampling.

  • 3.

    To simulate the uncertainty of the observation, draw a normally distributed error to each star-planet separation in the bootstrapped distribution using the standard deviation of the observations. We show these distributions for the simulated planets in system 94, simulation 4e in Fig. 7.

  • 4.

    Measure for each planet the quantile for the observed star-planet separation in the created distribution. The quantile is the fraction of outputs residing to the left of the observed star-planet separation. The star-planet separations are indicated with black, dark grey, and light grey lines in Fig. 7 for Earth-HR 8799 distances of 39.4 pc, then + 1 and + 2σ on the distance estimate.

The quantiles are the measures for how well the observed planet position and the simulated planet positions fit. Quantiles close to 0.5 would indicate a very good fit. We assess systems with quantiles between 10-4 and 1−10-4 as systems fitting to the observations. A low limit for the quantiles means that we allow systems that rarely are found with the architecture of the HR 8799 system to also count as fitting systems. Our example system, system 94, simulation 4e actually has a lower number for the quantile for planet b and does not count as a fitting system.

A visualisation of the assessment procedure can be seen in Fig. 7 where we show, for system 94, simulation 4e, the star-planet distributions plotted together with the star-planet distances inferred from the Earth-HR 8799 distances 39.4, 40.5, and 41.6 pc (corresponding to the mean of the estimate, then one and two standard deviations away from the mean). The quantiles for system 94, simulation 4e are qe = 0.89, qd = 0.89, qc = 0.02 and qb = 8.7 × 10-5 using the mean of the Earth-HR 8799 distance estimate. If we consider the distance one standard deviation away from the estimate (40.5 pc to HR 8799), however, we find the quantiles qe = 0.98, qd = 0.99, qc = 0.16 and qb = 0.01. At this distance we consider the system 94, simulation 4e to fit the observations well. To not miss systems such as system 94, simulation 4e in our analysis, we consider two values of the Earth-HR 8799 mean of the estimate distance (39.4 pc) and the distance one standard deviation away from the mean of the estimate (40.5 pc).

thumbnail Fig. 6

System 94, simulation 4e at 37.6 Myr (small black dots) compared to the observations of Currie et al. (2011) (coloured, larger dots). The axis around which the system is inclined is shown with a dashed line (position angle 64° and inclination 26° , Matthews et al. 2014). The position of the star is indicated with an asterisk.

Open with DEXTER

We summarise our findings of simulation 4 in the top two panels of Fig. 8. Again we have binned the simulation 4 with measured minimum orbit separation, Δmin, 100 systems in each bin. This makes in total 9 bins and 900 systems. We plot in light blue the fraction of the systems in each bin that have a stability timescale longer than the estimated age of HR 8799. In dark blue we show the fraction of systems in each bin that fit the observations of the HR 8799 system according to the above described method. The grey bars show the fraction of systems in each bin that fulfil both criteria and thus fit the HR 8799 system. We find for the nominal Earth-HR 8799 distance (39.4 pc), in the best case (where Δmin = 3.64−3.67) 5 systems of 100 that fit the observations. In the entire range of orbit separations, we find 12 of 900 systems that fit the observations (see Fig. 8a). For the Earth-HR 8799 distance of 40.5 pc we find, in the best case (where Δmin = 3.74−3.79) 8 systems of 100 fitting the observations and 22 of 900 when including all orbit separation ranges (see Fig. 8b).

In simulation 4 we create many varieties of architectures in the planetary systems that they would fit many ranges of inclinations and position angles. Therefore we are not concerned by the exact values of the inclination and position angle that we apply (from Matthews et al. 2014).

4.3. Simulations 2α and 4α: verification

thumbnail Fig. 7

Comparison of simulations to observed planet positions for planet e, d, c, and b subsequently. The histograms show the star-planet distance in our system 94, simulation 4e, after having moved all output to the observed position angle of the planet, bootstrapped a data set, and added a normally distributed error to the simulation data. We show, in black, dark grey, and light grey lines, the observed on-sky star-planet distances from Currie et al. (2011) assuming the mean Earth-HR 8799 distance (39.4 pc) and adding 1 and 2σ of this distance estimate (van Leeuwen 2007). This procedure is described in detail in Sect. 4.2.

Open with DEXTER

thumbnail Fig. 8

We show how many of our systems in simulation 4 (top and middle panels) and simulation 5 (bottom panel) fit the observations. We binned the simulations in ranges of minimum orbit separations, Δmin, 100 systems in each bin. The light blue shows the fraction of the systems in the bin that have a stability timescale longer than the estimated age of HR 8799. The dark blue curve shows the fraction of the systems in the bin that have architectures fitting to that of HR 8799. We determine whether a system fits through the procedure described in Sect. 4.2. The grey bars indicate the fraction of systems in each bin that fulfil both criteria. The top panel shows the results from simulation 4 using the mean of the Earth-HR 8799 distance estimate (39.4 pc), the middle panel shows the results from simulation 4 using a one standard deviation longer Earth-HR 8799 distance (40.5 pc), and the bottom panel shows the results from simulation 5, which has slightly lower planet masses, but is still within the estimated range.

Open with DEXTER

thumbnail Fig. 9

Comparison between simulations 4α (violet) and 2α (yellow). The lines show the number of systems in the simulation that have not experienced a close encounter, plotted with time. The slightly wider orbits in simulation 4α compared to simulation 2α makes the systems close encounter time shift about an order of magnitude. The black line and the grey zone show the age estimate and uncertainty of HR 8799.

Open with DEXTER

System 94, simulation 4e (see e.g. Fig. 5b) is not a special system, but there are many evolutions leading to systems that fit the observations of HR 8799. To show this, we simulate 100 copies of system 94, simulation 4e with the only difference being the simulation output rate, which changes the timestep slightly. This very small timestep variation changes the fate of the systems in ~ 106 years because they are chaotic. The systems have similar architecture, but with different evolution. We call them simulation 4α.

Because our simulated systems are very chaotic it is likely that an integration made with our exact initial conditions would evolve differently from that which we present here. A small change in timestep alters the evolution on a timecale of 106 yrs and our interesting systems have a stability timescale much longer than that. Therefore, we note that our initial conditions in Table A.1 are likely to not reproduce our system 94, simulation 4e. Instead, we present our recipes on how to create similar systems.

We also create simulation 2α from system 8, simulation 2, using the same procedure as when creating simulation 4α. Figure 9 shows with time the number of systems in simulations 2α and 4α that have not experienced any close encounter. The slopes of the curves show that the systems indeed have different evolutions despite their identical initial conditions. The location of the curves in the diagram shows that the stability timescale of simulation 4α is about one order of magnitude longer than the stability timescale of simulation 2α. This leads to 26 systems in simulation 4α surviving longer than the estimated age of HR 8799 and 10 of these also matching the observations of HR 8799.

4.4. Mean motion resonance analysis

In this section we investigate whether resonances affect our systems in simulation 4. Stabilisation because the system is deep in a mean-motion resonance is a common explanation to the survival of the HR 8799 system (see e.g. Fabrycky & Murray-Clay 2010; Goździewski & Migaszewski 2009, 2014; Reidemeister et al. 2009; Moro-Martín et al. 2010; Currie et al. 2011; Soummer et al. 2011; Esposito et al. 2013). In the following we show that we simulate systems that fit the HR 8799 system without being stabilised by strong resonant lock.

We consider resonant angles correlated to the Laplace resonances, 2:1, 4:2:1, and 8:4:2:1 (Murray & Dermott 1999). We also consider angles related to the secular motion. All these angles are calculated using the mean longitude (λ) and longitude of periastron (ϖ), which are calculated from the longitude of ascending node (Ω), argument of periastron (ω), and mean anomaly (M) for each planet in the following way: The resonant angles we consider of the 2:1 resonance are which for the three planet pairs (b-c, c-d and d-e) result in six angles to consider.

For the 4:2:1 resonance we consider the angle (7)which results in two angles as there are two sets of three-planet systems; b-c-d and c-d-e.

For the 8:4:2:1 resonance we take into account the four angles as follows: which we now label with the planet names.

We also look at the difference in longitude of periastron, which is an angle related to secular motion. It is defined as (12)which leads to three angles from the three planet pairs (b-c, c-d and d-e).

In Fig. 10 we show the evolution of four of the above angles in system 94, simulation 4e. We want to stress that system 94, simulation 4e is a system similar to many of our other systems in the resonant behaviour seen in Fig. 10 also. The evolution of the angles in our panels 2, 3, and 4 in Fig. 10 are also shown in the mean motion resonance analysis of Goździewski & Migaszewski (2009, their Figs. 4, middle panel, and 6, top panel) and Goździewski & Migaszewski (2014, their Fig. 1, bottom right panel). The previously presented solutions (see e.g. Goździewski & Migaszewski 2009, 2014; Fabrycky & Murray-Clay 2010) all show small-amplitude librations of one or more resonant angles. In contrast, our solution shows at most transitions between large-amplitude libration and circulation. This configuration does not guarantee the indefinite protection from close encounters that small-amplitude librations provide.

To show the overall behaviour of our simulation 4 we bin the systems by stability time into 9 bins, each containing 100 systems. For each system we assess whether the resonant angles presented in Eqs. (5)(12) are circulating, librating, or falling in and out of resonance. We visualise the resonant behaviour of our simulation 4 in Fig. 11. Figure 11 shows no trend of longer stability time with librating angles (right panel) and no clear trend of stability time angles transitioning between librational and circulatory behaviour. We infer from this that our systems are not being stabilised by tight resonant lock. However, the separatrix-crossing behaviour we observe in some systems is characteristic of chaotic evolution driven by multiple resonances (Chirikov 1979). Indeed, this may offer an explanation as to why there is a long-lived tail of surviving systems with no consistent resonant protection: Shikita et al. (2010) describe how some trajectories stick for extended periods at the edge of the chaotic region of phase space.

thumbnail Fig. 10

Evolution of (from top to bottom) the resonant angles φ2,1c:2d, φ1,1b:2c:4d, φ1,1b:2c:4d:8e and ϖcϖd for system 94, simulation 4e. The black dots show the output from our simulation, sampled every 100 years. In the top panel (evolution of the φ2,1c:2d) the output clusters around 0 rad, which is a sign of weak resonant behaviour. The system alternates between libration and circulation in this resonant angle, which is a different behaviour than what is seen in systems affected by strong resonance in which all output would be narrowly confined about 0 rad. None of the bottom three panels show signs of resonant libration, unlike in the cases considered by Goździewski & Migaszewski (2009, 2014). The red line indicates the estimated age of HR 8799.

Open with DEXTER

thumbnail Fig. 11

Resonant behaviour of our simulation 4, binned with stability time (x-axis) and 100 systems per bin. Colour shows for different resonant angles (y-axis) how many systems in the bin are circulating (left panel), librating (right panel) or transitioning between a resonant and non-resonant state (middle panel). We consider the angles in Eqs. (5)(12) (see Sect. 4.4). Our stability times in simulation 4 stretch between 0.04 Myr up to stable systems after 100 Myr. We see a resonant behaviour in the 2:1 resonant angles where it is common that the systems transition between libration and circulation. However, there seems to be no dependence on stability time with resonances.

Open with DEXTER

4.5. Simulations 5 and 6: lower planet masses

Another way to increase the stability timescale of a planetary system is to decrease the planet masses. In terms of Δ (Eqs. (1) and (2)), lowering planet masses widens the orbit separations and thus increases the stability timescale.

We use the exact same initial conditions for simulation 4 and lower the planet masses to create simulation 5 and 6 (see Table 1). Simulation 5 has planet masses 4 (b) and 6 (cde) MJup and simulation 6 has planet masses 3 (b) and 5 (cde) MJup. We want to point out that these planet masses are still comfortably within the estimated masses from cooling models (5 ± 2MJup (b) and (cde); Marois et al. 2008, 2010).

An example of a system that fits from simulation 5 is system 24, simulation 5d seen in Figs. 12 and 13. The system does not experience any close encounter during the 100 Myr integration time and we measure the quantiles to be qe = 0.93, qd = 0.98, qc = 0.15 and qb = 0.07 at the nominal Earth-HR 8799 distance 39.4 pc. We consider this system to be a good fit to the observations of HR 8799. The system is not stabilised by low-amplitude librational behaviour.

Using the technique described in Sect. 4.2 we find that for simulation 5, in the best case (where Δmin = 3.81−3.85) 18 systems of 100 fit the observations. Looking at a broader range of orbit separations in total 82 out of 900 systems fit the observations (see Fig. 8c). Considering these 1 MJup lower planet masses in simulation 5 the architecture of HR 8799 is becoming common in our systems also after 30 Myr. The boost of fitting systems in simulation 5 is clearly seen in Fig. 8 when comparing the bottom panel (simulation 5) with the top two panels (simulation 4).

We include simulation 6 to prove the point that lower planet masses implies much longer stability timescale. In simulation 6 we find 64 fitting systems of 100 in the best case (when Δmin = 3.98−4.02) and 241 of 900 when considering all orbit separation ranges.

Considering a slightly larger Earth-HR 8799 distance increased the number of fitting systems with a factor of two (Sect. 4.1). Considering slightly lower planet masses boosts the number of fitting systems by a factor of ten.

5. Evolutionary phases

In this section we use simulation 4f_long (see Table 1) to investigate the evolution of wide-orbit planetary systems. Simulation 4f_long is created in the same way as simulation 4f, but the integrations are not stopped after a close encounter has occurred and therefore always continue up to 100 Myr.

We find three distinct phases of dynamical evolution in our wide-orbit planetary systems. We call the phases A, B, and C and they correspond to the time before the first close encounter between planets (A), the time between the first and last close encounter (B), and the time after the last close encounter (C). We define a close encounter between planets to be whenever two planets enter their mutual Hill radius, see Eq. (1). We can use these evolutionary phases to predict what will happen and what has happened to the HR 8799 system and systems similar to it.

thumbnail Fig. 12

System 24 in simulation 5d (using lower planet masses of 4 MJup for planet e and 6 MJup for planets b, c, and d). Black dots are our output data points assuming the inclination and position angle of 26° and 64° from Matthews et al. (2014). The coloured dots are the observations of the HR 8799 planetary system. The system matches the data and does not experience any close encounter during our 100 Myr integration time.

Open with DEXTER

thumbnail Fig. 13

Same as Fig. 7 but for the system 24, simulation 5d (planet masses are lower than nominal, M = 4 MJup(b), and 6 MJup (cde)).

Open with DEXTER

thumbnail Fig. 14

Eccentricity and semi-major axes of the planets in each system of simulation 4f_long at a randomly picked time for each system in each phase. From top to bottom we show the cases of 1-, 2-, 3-, and 4-planet systems and from left to right we show phase A, B, and C. All phase A systems have low eccentricity and four planets. In phase B, the systems have three or more planets. Typically, systems in phase C have two planets, one at a ~ 10 AU and one at 30 ≲ a ≲ 1000 AU. Both phase B and C show generally higher eccentricities.

Open with DEXTER

Figure 14 shows how the eccentricity and semi-major axes of the planetary orbits change between phase A, B, and C. In the figure we picked one random time in phase A, B, and C in each system in simulation 4f_long, respectively.

Phase A is characterised by low-eccentricity orbits, not deviating much from the initial set-up. Phase B and C show an even distribution of eccentricities (0 <e< 1) and planets pile-up around a ~ 10 AU. All systems in phase B have three planets or more, while the most common number of planets in phase C is two planets. In the two-planet case of phase C, one planet typically has a semi-major axis of a ~ 10 AU and the other a ~ 30−1000 AU. The timescale of phase A is what we referred to as the stability timescale, which is the time it takes before the first close encounter occurs. This timescale is strongly dependent on orbital separation (Davies et al. 2014). The timescale of phase B is independent of initial orbit separation and we measure it to vary between 104 to 107 years, peaking at 106 years. These timescales for the clearing of wide-orbit planets are consistent with previous studies (Veras et al. 2009). This timescale is typically shorter than the duration of phase A and the lifetime of the observed systems (see Table 2). Therefore statistically it is unlikely to observe a system in phase B unless the star is very young.

5.1. Application to observed wide-orbit planetary systems

In this section we predict the evolutionary phase of directly imaged, massive, wide-orbit planetary systems assuming their evolution is similar to that which we find for HR 8799 in simulation 4f_long. Table 2 gives some observed properties of the systems we include.

1RXS J1609

1RXS J1609 b is found at a projected distance of 330 AU and is probably scattered to this distance as a protoplanetary disk of this size is unlikely. Its age is not sufficient to distinguish between phase B or C, but if there are no other planets present we can predict phase C.

AB Pic

AB Pic b is also found at a very large distance from the star (~ 260 AU). It is highly likely the planet has been scattered there and also taking the age of ~ 30 Myr into account, we predict phase C as phase B generally lasts for shorter time.

β Pic

β Pic b is separated from β Pic by 8−15 AU. The fact that no outer planets have been found indicates either that the system is a one-planet system in phase C or that the outer planets are not massive enough to be visible and the system is in phase A or C. However, the low eccentricity of the orbit (e ~ 0.06, Macintosh et al. 2014) favours phase A.

Fomalhaut

Fomalhaut b is found at a projected distance of ~ 119 AU, which is far out and makes in situ formation less probable. The predicted eccentric orbit necessary to fit the observations would indicate phase C considering the advanced age of the system (~ 440 Myr).

GJ 504

The recently imaged planet around the solar-type star GJ 504 (Kuzuhara et al. 2013) is found at a projected star-planet distance of 43.5 AU. GJ 504 b could be an outer planet in a phase C system (consider the age of ~ 160 Myr), which means that there might be another planet on a closer orbit.

HD 106906

HD 106906 b has an on-sky projected separation of ~ 650 AU, which makes it clear that the planet is scattered and the system is in phase B or C as such a large protoplanetary disk is unlikely. It may be that there is also an inner planet in the HD 106906 system as the inner part cannot yet be resolved.

HD 95086

HD 95086 b has a projected distance to the star of 56.4 ± 0.7 AU and could be an outer planet in a phase B or C configuration, which then indicates possible inner companions. It could also be a phase A system with an unresolved inner region because of the distance to the star (90.4 ± 3.4, van Leeuwen 2007). The debris disk of HD 95086 is found to be similar to that of HR 8799 (Su et al. 2015).

HR 8799

The number of planets in HR 8799 and their low eccentricities make us classify the system as still being in phase A.

κ And

The κ And system seems like a higher mass version of the systems HD 95086 and GJ 504. The inner regions might hide a phase A system whilst the outer regions might reveal an outer planet indicating phase C. The system is old enough for us to predict phase B to be unlikely.

6. Summary and conclusions

The star HR 8799 is observed to be accompanied by four massive planets on wide orbits. Puzzlingly, the planetary system has been shown to be unstable in a timescale much shorter than its estimated age. A common explanation of its existence has been stabilisation by mean motion resonances (Goździewski & Migaszewski 2014).

We have simulated long-lived systems similar to HR 8799 without forcing them into resonance. The initial orbits of our planets are slightly more widely separated than circular orbits that intersect the observed star-planet on-sky separations. Slightly larger separations between the orbits result in a significantly longer stability timescale. The low eccentricity of the orbits can in some cases make a system have the observed architecture of HR 8799. We place the planets with random azimuth angle, assuming initial low eccentricities.

When using the estimated planetary masses (5, 7, 7 and 7 MJup) together with the estimated Earth-HR 8799 distance we rarely find a fitting system; in the best case, Δmin = 3.64−3.67, we find 5 of 100 fitting systems, and in the broader range of orbit separations we find 12 out of 900 fitting systems. We almost double the number of fitting systems by considering the 1σ longer distance to HR 8799 of 40.5 pc; the best case now considers orbit separation range Δmin = 3.74−3.79 and yields 8 of 100 fitting systems: in total we find 22 of 900 fitting systems. It is easier to find a fitting system by considering slightly lower planet masses (4, 6, 6, and 6 MJup) that are still comfortably within the planet mass estimate. That results in the best case (Δmin = 3.98−4.02) 18 of 100 systems fit the observations and in the broader range of orbit separations we find 82 out of 900 systems fitting. The above described numbers are visualised in Fig. 8. Examples of a system that fits the observations can be seen for the nominal planet masses in Figs. 5b, 6, and 7 and for the slightly lower planet masses in Figs. 12 and 13.

We test our systems for mean motion resonances (see Sect. 4.4) and find that most of them fall sporadically in and out of the 2:1 resonance between various planet pairs, however, we see no trend in long stability timescale in the systems with more resonant behaviour. This means that these resonances do not work in a stabilising way, unlike the resonances treated earlier in the literature for stability of the HR 8799 system.

We conclude that the HR 8799 system does not necessarily have to be caught in a strong resonance. We simulated many systems that at some point in their evolution look like HR 8799 and we want to stress that our example system (system 94, simulation 4e) is not a lucky coincidence. We present a different solution than before as our fitting systems are stable for longer than the lifetime of HR 8799 and do not show signs that they are deep in resonance.

A bonus from our simulations is a prediction of the future evolution of the HR 8799 system. We see typical architectures in different evolutionary phases of the system. Before close encounters between the planets, the orbits have low eccentricity (e< 0.2). After the first close encounter has occurred, the first planet ejection follows within ~ 106 years. The close encounter and ejection change the architecture of the system and the eccentricities are now more evenly distributed (0 <e< 1). In ~ 80% of the cases, the systems eventually lose another planet and end up with two planets on eccentric orbits: one with a semimajor axis of ainner ~ 10 AU and the other at aouter ~ 30−1000 AU (see Fig. 14). This configuration is comparably more similar to other massive, wide-orbit planets, such as AB Pic b or HD 106906 b, than the current HR 8799 bcde.

Our HR 8799 analogues naturally evolve into architectures that resemble other wide-orbit planetary systems. Therefore we find it most plausible that the HR 8799 system is not of a different kind than other observed wide-orbit planetary systems, but simply in an earlier evolutionary phase. This phase seems to be relatively short for massive, wide-orbit planetary systems and it thus makes sense that only one system has been observed in this phase, that is HR 8799.

Acknowledgments

The authors thank the anonymous referees for their comments which improved the manuscript. The authors thank Phil Uttley for expertise in statistics. Y.G. acknowledges the support of a Ph.D studentship at the University of Amsterdam under the supervision of Selma E. de Mink. A.J. and Y.G. were partially funded by the European Research Council under ERC Starting Grant agreement 278675–PEBBLE2PLANET and by the Swedish Research Council (grant 20103710). MBD was supported by the Swedish Research Council (grant 2011-3991). R.P.C. was supported by the Swedish Research Council (grants 20122254 and 20125807). Calculations presented in this paper were carried out at LUNARC using computer hardware funded in part by the Royal Fysiographic Society of Lund. A.M. was supported by the Swedish Research Council (grant 20113991) and KAW 2012.0150 from the Knut and Alice Wallenberg foundation.

References

Appendix A: Initial conditions of system 94, simulation 4e

Table A.1

Initial conditions of system 94, simulation 4e.

All Tables

Table 1

Initial conditions of the simulations.

Table 2

Observed wide-orbit planets around stars.

Table A.1

Initial conditions of system 94, simulation 4e.

All Figures

thumbnail Fig. 1

Observations of the HR 8799 planets. The black dashed line corresponds to the position angle of the inclination as estimated by Matthews et al. (2014) with the ± 1σ uncertainty shown by the grey region.

Open with DEXTER
In the text
thumbnail Fig. 2

Number of planets left in the systems as a function of time. Dotted line corresponds to the face-on assumption (simulation 1), solid line to the inclination inferred from observations (simulation 2), and dashed line to the largest separation achieved from inclination (simulation 3). See Table 1 and Sect. 3 for descriptions about the simulations. The age estimate of HR 8799 ( Myr) is indicated with a solid line with the grey region corresponding to the uncertainty.

Open with DEXTER
In the text
thumbnail Fig. 3

Smallest orbit separation, Δmin, between any planet pair changes with inclination and position angle. We indicate the inclination and position angle of the HR 8799 system measured by Matthews et al. (2014) with a white square (simulation 2), which is i = 26 ± 3° and position angle 64 ± 3°. The largest Δmin occurs at i = 25° and position angle 42°, when Δmin = 3.56RH. This value is indicated with a filled, white circle (simulation 3). This plot is based on the observations from Marois et al. (2008, 2010).

Open with DEXTER
In the text
thumbnail Fig. 4

Number of systems without close encounters shown with time. The different coloured lines correspond to bins of each 100 systems in simulation 4, sorted in measured minimum orbit separations (Δmin). The legend indicates the bins and their systems’ mean value of Δmin. The black line and the grey zone show the age estimate with uncertainty of HR 8799.

Open with DEXTER
In the text
thumbnail Fig. 5

Comparison between a system in simulation 2 and a system in simulation 4e. Both systems seem to fit the observed positions of the planets, but only the system in simulation 4e is stable for longer than the estimated age of HR 8799. This is an on-sky projection of the simulations assuming a distance of 39.4 pc (van Leeuwen 2007) and inclination 26° (position angle 64° , Matthews et al. 2014). The axis around which the system is inclined is indicated with a dashed line. The position of the star is indicated with an asterisk.

Open with DEXTER
In the text
thumbnail Fig. 6

System 94, simulation 4e at 37.6 Myr (small black dots) compared to the observations of Currie et al. (2011) (coloured, larger dots). The axis around which the system is inclined is shown with a dashed line (position angle 64° and inclination 26° , Matthews et al. 2014). The position of the star is indicated with an asterisk.

Open with DEXTER
In the text
thumbnail Fig. 7

Comparison of simulations to observed planet positions for planet e, d, c, and b subsequently. The histograms show the star-planet distance in our system 94, simulation 4e, after having moved all output to the observed position angle of the planet, bootstrapped a data set, and added a normally distributed error to the simulation data. We show, in black, dark grey, and light grey lines, the observed on-sky star-planet distances from Currie et al. (2011) assuming the mean Earth-HR 8799 distance (39.4 pc) and adding 1 and 2σ of this distance estimate (van Leeuwen 2007). This procedure is described in detail in Sect. 4.2.

Open with DEXTER
In the text
thumbnail Fig. 8

We show how many of our systems in simulation 4 (top and middle panels) and simulation 5 (bottom panel) fit the observations. We binned the simulations in ranges of minimum orbit separations, Δmin, 100 systems in each bin. The light blue shows the fraction of the systems in the bin that have a stability timescale longer than the estimated age of HR 8799. The dark blue curve shows the fraction of the systems in the bin that have architectures fitting to that of HR 8799. We determine whether a system fits through the procedure described in Sect. 4.2. The grey bars indicate the fraction of systems in each bin that fulfil both criteria. The top panel shows the results from simulation 4 using the mean of the Earth-HR 8799 distance estimate (39.4 pc), the middle panel shows the results from simulation 4 using a one standard deviation longer Earth-HR 8799 distance (40.5 pc), and the bottom panel shows the results from simulation 5, which has slightly lower planet masses, but is still within the estimated range.

Open with DEXTER
In the text
thumbnail Fig. 9

Comparison between simulations 4α (violet) and 2α (yellow). The lines show the number of systems in the simulation that have not experienced a close encounter, plotted with time. The slightly wider orbits in simulation 4α compared to simulation 2α makes the systems close encounter time shift about an order of magnitude. The black line and the grey zone show the age estimate and uncertainty of HR 8799.

Open with DEXTER
In the text
thumbnail Fig. 10

Evolution of (from top to bottom) the resonant angles φ2,1c:2d, φ1,1b:2c:4d, φ1,1b:2c:4d:8e and ϖcϖd for system 94, simulation 4e. The black dots show the output from our simulation, sampled every 100 years. In the top panel (evolution of the φ2,1c:2d) the output clusters around 0 rad, which is a sign of weak resonant behaviour. The system alternates between libration and circulation in this resonant angle, which is a different behaviour than what is seen in systems affected by strong resonance in which all output would be narrowly confined about 0 rad. None of the bottom three panels show signs of resonant libration, unlike in the cases considered by Goździewski & Migaszewski (2009, 2014). The red line indicates the estimated age of HR 8799.

Open with DEXTER
In the text
thumbnail Fig. 11

Resonant behaviour of our simulation 4, binned with stability time (x-axis) and 100 systems per bin. Colour shows for different resonant angles (y-axis) how many systems in the bin are circulating (left panel), librating (right panel) or transitioning between a resonant and non-resonant state (middle panel). We consider the angles in Eqs. (5)(12) (see Sect. 4.4). Our stability times in simulation 4 stretch between 0.04 Myr up to stable systems after 100 Myr. We see a resonant behaviour in the 2:1 resonant angles where it is common that the systems transition between libration and circulation. However, there seems to be no dependence on stability time with resonances.

Open with DEXTER
In the text
thumbnail Fig. 12

System 24 in simulation 5d (using lower planet masses of 4 MJup for planet e and 6 MJup for planets b, c, and d). Black dots are our output data points assuming the inclination and position angle of 26° and 64° from Matthews et al. (2014). The coloured dots are the observations of the HR 8799 planetary system. The system matches the data and does not experience any close encounter during our 100 Myr integration time.

Open with DEXTER
In the text
thumbnail Fig. 13

Same as Fig. 7 but for the system 24, simulation 5d (planet masses are lower than nominal, M = 4 MJup(b), and 6 MJup (cde)).

Open with DEXTER
In the text
thumbnail Fig. 14

Eccentricity and semi-major axes of the planets in each system of simulation 4f_long at a randomly picked time for each system in each phase. From top to bottom we show the cases of 1-, 2-, 3-, and 4-planet systems and from left to right we show phase A, B, and C. All phase A systems have low eccentricity and four planets. In phase B, the systems have three or more planets. Typically, systems in phase C have two planets, one at a ~ 10 AU and one at 30 ≲ a ≲ 1000 AU. Both phase B and C show generally higher eccentricities.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.