You need an eReader or compatible software to experiencethe benefits of the ePub3 file format.

Recent works on planetary migration show that the orbital structure of the Kuiper Belt can be very well reproduced if, before the onset of planetary instability, Neptune underwent a long-range planetesimal-driven migration up to ∼28 au. However, considering that all giant planets should have been captured in mean motion resonances among themselves during the gas-disk phase, it is not clear whether such a very specific evolution for Neptune is possible, or whether the instability could have happened at late times. Here, we first investigate which initial resonant configuration of the giant planets can be compatible with Neptune being extracted from the resonant chain and migrating to ∼28 au before planetary instability. We address the late instability issue by investigating the conditions where the planets can stay in resonance for about 400 Myr. Our results indicate that this can happen only in the case where the planetesimal disk is beyond a specific minimum distanceδ_stab from Neptune. Then, if there is a sufficient amount of dust produced in the planetesimal disk, which drifts inwards, Neptune can enter a slow dust-driven migration phase for hundreds of Myr until it reaches a critical distance from the disk. From that point, faster planetesimal-driven migration takes over and Neptune continues migrating outward until the instability happens. We conclude that although an early instability more easily reproduces the evolution of Neptune required to explain the structure of the Kuiper Belt, such evolution is also compatible with a late instability.

It is now well-accepted that the giant planets of the solar system did not have their current orbits when the gas was removed from the protoplanetary disk. The first work showing that the planets had to migrate while interacting with the remaining planetesimal disk was by Fernandez & Ip (1984). This work demonstrated that in the planetesimal-scattering process, Saturn, Uranus, and Neptune on average gain angular momentum and thus migrate outwards, while Jupiter loses angular momentum and moves inwards. Since then, the process of planetesimal-driven migration was considered to be fundamental in shaping the current structure of the solar system. Malhotra (1993,1995) showed that planetesimal-driven migration can explain the origin of the orbit of Pluto and the resonant populations of Kuiper Belt objects (KBOs). However, it soon became clear that planetesimal-driven migration alone is too smooth to explain several aspects of the Kuiper Belt’s orbital distribution, e.g., the orbital excitation of non-resonant objects. Thommes et al. (1999) was the first to show that dynamical instability of the giant planets could lead to global orbital excitation of the Kuiper Belt, while delivering the giant planets to final orbits roughly similar to their current ones.

Building on these results, the so-called Nice model (Gomes et al.2005; Morbidelli et al.2005; Tsiganis et al.2005) combined the processes of planetesimal-driven migration and dynamical instability into a new global scenario. It suggested that the planets Jupiter, Saturn, Uranus, and Neptune formed in a compact configuration with Jupiter around 5.45 au, with Saturn slightly interior to the 2:1 mean motion resonance (MMR) with Jupiter ( < 8.65 au) and the initial semimajor axes of the ice giants (Uranus and Neptune) in the ranges 11–13 au and 13.5–17 au. The planets were surrounded by a disk of leftover planetesimals, accounting for a total of about 35M_⊕, from ∼16 to 30 au (whereM_⊕ ∼ 5.97 × 10²⁴ kg is the mass of the Earth). In this model, the planets interacted with the planetesimal disk and, through planetesimal-driven migration, their mutual separations slowly increased. When Jupiter and Saturn crossed their mutual 2:1 MMR, the full system of giant planets became unstable and spread farther through mutual close encounters. In particular, Uranus and Neptune were scattered onto orbits with relatively large eccentricities and semimajor axes, which crossed the original trans-Neptunian disk. Then, following a phase of dynamical friction caused by their interaction with the disk, the orbits of Uranus and Neptune damped in eccentricity and inclination. A residual planetesimal-driven migration placed the planets on orbits fairly similar to their current ones (Tsiganis et al.2005).

The Nice model had great success in explaining many features of our solar system, such as the capture of Jupiter’s Trojans (Morbidelli et al.2005), the capture of the irregular satellites of the giant planets (Nesvorný et al.2007), the absence of regular moons beyond Oberon’s orbit at Uranus (Deienno et al.2011), to cite a few. Moreover, it was shown that, provided the inner edge of the trans-Neptunian disk was at an appropriate distance from the original orbit of Neptune, the giant-planet instability could have occurred after hundreds of Myr of evolution, potentially explaining the origin of the so-called Late Heavy Bombardment (LHB) of the terrestrial planets (Gomes et al.2005).

However, the original version of the Nice model presented a major problem: the initial configuration of the planets was chosen in an ad hoc manner. Morbidelli et al. (2007) performed hydrodynamical simulations of the evolution of the giant planets when they were still embedded in a disk of gas and found that the planets should have acquired a fully resonant configuration, with Jupiter and Saturn preferentially locked in their mutual 3:2 MMR with Jupiter at ∼5.4 au, Saturn and Uranus (ice I) in 3:2 MMR, and Uranus (ice I) and Neptune (ice II) in 4:3 MMR (although other resonant configurations are also possible). The authors also concluded that such a multiresonant system is compatible with a subsequent evolution similar to that described in Tsiganis et al. (2005), although the instability in this case happens when a pair of planets leave the original resonant configuration, rather than when Jupiter and Saturn cross their 2:1 MMR.

Levison et al. (2011) investigated whether, under some conditions, this initial multiresonant configuration could still lead to alate giant-planet instability, in order to explain the putative LHB as in the original Gomes et al. (2005) paper. They found that late instabilities are possible, provided that the inner edge of the disk is several astronomical units beyond the outermost planet and about 1000 Pluto-sized planetesimals are embedded in the disk, so as to produce an important dynamical self-stirring of the planetesimal population.

During the dynamical instability of the giant planets, several dynamical paths are possible, even when restricting the evolutions to those leading to final planetary orbits similar to the current ones. The actual dynamical evolution of the giant planets at the time of instability can be constrained by investigating the consequences on the other components of the solar system: the terrestrial planets and the small-body populations.

The terrestrial planets and the asteroid belt give similar constraints, as investigated in Brasser et al. (2009) and Morbidelli et al. (2010). If the orbits of Jupiter and Saturn migrate away from each other smoothly and slowly, as in planetesimal-driven migration, theg₅ secular frequency (dominant in the precession of Jupiter’s longitude of perihelion) decreases slowly in value and becomes very similar for some time to theg₁ secular frequency, dominant in the precession of Mercury’s perihelion. As a consequence, the orbit of Mercury is destabilized (Brasser et al.2009). Similarly, theg₆ secular frequency, dominant in the precession of Saturn’s perihelion, becomes equal to the frequencies characterizing the precession of the asteroids at low inclination in the inner asteroid belt as it slows down, so that this region becomes severely depleted (Morbidelli et al.2010). Both papers concluded that the ratio of the orbital periods of Jupiter and Saturn had to evolve almost discontinuously from <2.1 to >2.3 (the current value is). This happens if Jupiter and Saturn encounter another planet (Uranus or Neptune or a third planet of comparable mass). The sequence of encounters causes a divergent jump in the semimajor axes of the orbits of Jupiter and Saturn, and therefore a jump in the ratio. This scenario is commonly known as the jumping-Jupiter model. Notice that if the giant-planet instability occurred early on, the terrestrial planets were possibly not yet fully formed and therefore the constraint on the evolution ofg₅ relative tog₁ does not apply (Kaib & Chambers2016). However, the constraint set by the asteroid belt is still valid, supporting the need for a jumping-Jupiter evolution regardless of the time of the instability (Walsh & Morbidelli2011; Toliou et al.2016).

Given that in a jumping-Jupiter evolution the planet that encounters Jupiter is often ejected onto a hyperbolic orbit, Nesvorný (2011) and Batygin et al. (2012) suggested that the solar system originally had a third ice-giant planet with mass comparable to Uranus or Neptune. This scenario increases the probability that a jumping-Jupiter evolution ends with four planets near their current orbits, once the putative fifth planet is ejected after its encounter with Jupiter.

Then, Nesvorný & Morbidelli (2012) considered many different initial multiresonant configurations for Jupiter, Saturn, Uranus, Neptune, and the rogue fifth planet, and determined for each of them the probability that the ratio had a jump across the 2.1–2.3 range and the final planets reached orbits compatible with those of our current solar system. The successful simulations of Nesvorný & Morbidelli (2012) have then been used to investigate again the capture of the Trojans of Jupiter (Nesvorný et al.2013), the capture of irregular satellites at Jupiter (Nesvorný et al.2014a), the excitation of the inclination of the moon Iapetus at Saturn (Nesvorný et al.2014b), and the evolution of the asteroid main belt and of its collisional families (Roig & Nesvorný2015; Brasil et al.2016; Deienno et al.2016). They have also been shown to be compatible with the survival of the Galilean satellites (Deienno et al.2014), despite the encounters of Jupiter with the ejected planet.

However, when investigating the dynamical sculpting of the Kuiper Belt, Nesvorný (2015a,2015b) and Nesvorný & Vokrouhlický (2016) realized that, of all possible evolutions of Neptune consistent with the five-planet jumping-Jupiter model, only some are consistent with the current orbital structure of the trans-Neptunian population. In order to get the correct inclinations of the hot population of the KBOs, Nesvorný (2015a) found that Neptune should have migrated more than 5 au from its original resonance due to planetesimal-driven migration before the planets became unstable. This migration should have occurred on a timescaleτ ≥ 10 Myr, with the planet on a quasi-circular planar orbit ( and) to avoid excessive orbital excitation of the cold classical belt, in the ∼42–45 au region. The planetary instability should have happened when Neptune was already at 28 au (Nesvorný2015b) in order to explain the so-called Kernel of the cold Kuiper Belt (a clump of objects around 44 au; Petit et al.2011). In this case the Kernel would have formed from objects transported outwards in 2:1 MMR with Neptune and released from the resonance when Neptune’s orbit had a jump of about 0.5 au in semimajor axis, due to a close encounter with another planet. Finally, in order to reproduce the observed ratio between the resonant and non-resonant Kuiper Belt populations, Nesvorný & Vokrouhlický (2016) concluded that the planetesimal-driven migration of Neptune should have been characterized by several small amplitude jumps, such as those caused by a population of 1000–4000 Pluto-size objects in the disk.

The brief review reported above sets the context that motivated the study reported in this paper.

It is clear that the work by Nesvorný & Morbidelli (2012) needs to be revisited to test which initial resonant configurations are the most compatible with the new constraints on the dynamical evolution of Neptune. This is done in the first part of the paper (Section2). In the second part of the paper (Section3) we investigate whether the initial configurations selected in the first part of the work can lead to alate planetary instability, while still showing an evolution of Neptune consistent with Kuiper Belt constraints.

The possibility of alate instability is far from trivial. In fact, in order to exhibit long-range planetesimal-driven migration before planet instability, Neptune needs to be embedded in the planetesimal disk. But in this case, Neptune would start migrating away from the other planets from the very beginning of the simulation (corresponding to the gas-removal time). This is likely to lead to an instability after just a few tens of Myr. If instead the planetesimal disk is relatively far from Neptune, the planet can be extracted from its initial resonance at alate time, as in Levison et al. (2011), but the instability is likely to happen as soon as the initial resonant configuration is broken, i.e., before Neptune can migrate far away. For this problem we will introduce a new concept: that of dust-driven migration, i.e., migration due to the gravitational interaction of a planet with the dust produced in a distant planetesimal disk, spiraling inward due to Poynting–Robertson drag. We will study the evolution of the system as a function of the dust properties and the numerical resolution used in modeling the planetesimal disk. Finally, our results are summarized and discussed in Section4.

Nesvorný & Morbidelli (2012) characterized the evolution of the giant planets from several possible initial resonant configurations. They assumed that the solar system could have formed with four, five, or six giant planets. The additional planets had a total mass comparable to that of Uranus or Neptune. They analyzed the rate of success of each configuration by measuring the fraction of the planetary orbital evolutions that fulfilled a set of constraints, namely, the final planetary system is made up of four giant planets (criterion A), the orbits of the four remnant planets resemble the present ones (criterion B), the proper mode of Jupiter’s eccentricity () is greater than 0.022 in the final system, i.e., at least half of its current value (criterion C), and the ratio evolves from <2.1 to >2.3 in less than 1 Myr (criterion D). Criterion C is important because it is found that although encounters between Jupiter and an ice giant can excite, not all encounters can do it as much as needed. Criterion D is important because it is linked with the evolution of the terrestrial planets and of the asteroid belt (Brasser et al.2009; Morbidelli et al.2010; Toliou et al.2016). If Jupiter and Saturn migrate too slowly through the period ratio 2.1–2.3, theg = g₅ andg = g₆ resonances sweep through the terrestrial planet region and the inner asteroid belt, resulting in orbital distributions inconsistent with those observed. Then, from the success rates, Nesvorný & Morbidelli (2012) concluded that a planetary system initially having five planets should be preferred.

Here, we revisit the Nesvorný & Morbidelli (2012) work, adding a new constraint to the orbital evolution of Neptune from the work of Nesvorný (2015b), namely that Neptune should have migrated to ∼28 au before the onset of planetary instability. The works by Nesvorný (2015a,2015b) and Nesvorný & Vokrouhlický (2016) also give other constraints on the migration timescale of Neptune and its eccentricity and inclination before the instability. The eccentricity and inclination constraints are met fairly easily because in this scenario, Neptune evolves in a planetesimal-driven migration without planetary close encounters until the instability phase. As for the migration timescale, Nesvorný (2015a) claims that before the planetary instability, once out of resonance, Neptune can migrate according toτ = 10 Myr for disks with 20M_⊕ and according toτ = 50 Myr for disks with 15M_⊕. Considering that those disks have the same extension, this exemplifies how sensitive the migration timescale can be with respect to the density of the disk. Therefore, we decide not to consider the migration timescale as an important parameter at this point, because getting a correct timescale of the migration is probably only a matter of tuning the density and profile of the planetesimal disk. Thus, we add just one new criterion (criterion E): at the time of the instability, has to be in the interval 27–29 au. Although criteria A to D indicate that the initial solar system most probably hosted five giant planets (Nesvorný & Morbidelli2012), criterion E is a very important criterion that can be used to constrain in which resonant chain these five planets likely were at the end of the gas-disk phase.

For simplicity, we consider only the cases of Nesvorný & Morbidelli (2012) that involved five planets. We also always considered Jupiter to be initially at ∼5.4 au, and the distance between the outermost ice giant and the inner edge of the planetesimal diskδ = 1 and 2 au and extending to 30 au (Gomes et al.2004), and the disk composed of 1000 planetesimals (which is the resolution considered by Nesvorný & Morbidelli2012), with two values of surface density () equal to 0.02 and 0.03. Such configuration implies that the mass of the planetesimal disk ranges from ∼13M_⊕ (in the more relaxed configuration—3:2, 3:2, 2:1, 2:1,⁵ with 2 au and) to ∼65M_⊕ (in the more compact configuration—3:2, 3:2, 4:3, 5:4, withδ = 1 au and).

To perform our simulations, we used the Hybrid integrator in theMercury package (Chambers1999), with a time step of half a year.

Table1 gives in column 6 the rate of success for our new criterion E, after 20 runs of each initial configuration from Nesvorný & Morbidelli (2012, Table 6 in their paper). Those 20 simulations were composed of 5 simulations for each pairδ = [1, 2] au, and. Although we have used two different values of surface density and two values ofδ for the planetesimal disk, we decided to compress the results of all 20 simulations for each configuration in Table1, column 6. This is because, concerning our criterion E, only minor differences were observed with these changes. Thus, when important, those minor differences will be reported. For completeness and to increase statistics, Table1 also reports the average of the rate of success for the other criteria taken from Nesvorný & Morbidelli (2012, Table 6 in their paper). Because the simulations performed in this work are essentially the same as in Nesvorný & Morbidelli (2012), the success rate for criteria A to D are consistent (within the small number statistics) with those reported in Nesvorný & Morbidelli (2012). Thus, we did not report our own results for the criteria A, B, C, and D. Rather, we decided to check only the success rate relative to the new criterion E, and then combine our findings with those of Nesvorný & Morbidelli (2012), which have much higher statistics, in order to determine which of the initial planetary configurations are the most compatible with the constraint on the dynamical evolution of Neptune, unveiled in Nesvorný (2015a,2015b).

Our results can be interpreted as follows:

2:1 Jupiter–Saturn resonance: In these configurations, as it also appears in Nesvorný & Morbidelli (2012, Figure 16), it is possible that Neptune leaves the resonant configuration and migrates ahead of the other planets, but not always in the way required in Nesvorný (2015a,2015b). Indeed, a configuration with initial resonances 2:1, 3:2, 3:2, 3:2 never leads to an instability that fulfilled criterion E in our 20 simulations. The evolutions from the 2:1, 4:3, 3:2, 3:2 configuration, on the other hand, had one of the largest rates of success for criterion E. However, such configurations present other major problems, mostly concerning the final orbital separation between Jupiter and Saturn, which is typically too big. The point is that because the ratio has to change by less than 0.5 if Jupiter and Saturn are initially in 2:1 MMR, the mass of the fifth planet should be ∼1/2 of the mass of Neptune. Instead, we assumed for this planet a mass equal to that of Uranus–Neptune, as in Nesvorný & Morbidelli (2012). A smaller mass planet, however, would not excite Jupiter’s eccentricity enough, so that criterion C would not be fulfilled. Moreover, hydrodynamical simulations (Morbidelli et al.2007; Pierens & Nelson2008) show that in a massive gas disk, Jupiter and Saturn would be locked in their mutual 3:2 MMR, rather than 2:1. The Grand Tack model (Walsh et al.2011) also implies that Jupiter and Saturn were in their 3:2 MMR at the disappearance of the disk of gas. So, we do not consider that configurations with Jupiter and Saturn in 2:1 MMR are very realistic.

3:2, 3:2, 4:3, 4:3(5:4): These two configurations are too compact. Only 1 out of 20 simulations for each configuration fulfilled criterion E. Usually, the evolution of the planetary system is as in Nesvorný & Morbidelli (2012, Figure 7), where the instability takes place as soon as one planet goes out of resonance. As a result, no large-range planetesimal-driven migration occurs before the planetary instability. Instead, the instability is violent and the evolution of Neptune is largely different from that constrained in Nesvorný (2015a,2015b) (Figure1 top left, in this work). Nevertheless, although they are low, none of the success rates reported in Table1 for the configuration 3:2, 3:2, 4:3, 4:3 are null. Thus, even being very unlikely, due to the reasons presented, this configuration cannot be ruled out.

Zoom InZoom OutReset image size

3:2, 3:2, 2:1, 2:1: This configuration is too loose. Nesvorný & Morbidelli (2012) already pointed out that Neptune stops farther than 30 au for whatever reasonable starting disk parameters for this configuration. Here, we found that even with the planetesimal disk of lowest mass (m_disk ∼ [13–24] M_⊕), ranging froma_disk ∼ [25(26)–30] au, Neptune goes out of resonance and usually migrates all the way up to 30 au or more, before the planetary instability happens (Figure1 top right). However, sometimes the instability can also occur before Neptune reaches 30 au, say, at ∼28 au. Therefore, although this configuration may also present some problems with the extension of the planetesimal disk, as will be discussed in the following sections, it does have a good rate of success in fitting criterion E (25% success rate).

3:2, 3:2, 3:2, 3:2: This configuration allows Neptune to migrate away from the other planets before the instability, but usually not far enough. The instability commonly happens when Neptune is around 20 au (see Nesvorný & Morbidelli2012, Figures 8 and 9, and Figure1 bottom left in this work). Only for high-density disks can Neptune migrate fast enough and reach ∼28 au before the planetary instability. If the disk’s density is low, Neptune migrates slowly and cannot go too far into the planetesimal disk. So, this configuration has a low, but non-negligible, success rate of 15% and cannot be ruled out.

3:2, 3:2, 2:1, 3:2: This is the configuration that presents the highest rate of success, considering all criteria together, thus it is highlighted in bold in Table1. This is also the most successful configuration relative to criterion E, with a success rate of 65% independent of the disk setup. The evolution of Neptune and the rest of the planetary system is usually similar to that shown in Nesvorný & Morbidelli (2012, Figures 12, 13, and 14), and satisfies the constraint identified in Nesvorný (2015a,2015b) (Figure1 bottom right in this work).

Figure1 shows examples of the evolution of the systems corresponding to the descriptions given for each of the resonant chains discussed above (with the exception of the 2:1 Jupiter–Saturn resonance, due to the considerations previously made). Almost all evolutions depicted in Figure1 fulfill criteria A, B, C, and D. Criterion C was not fulfilled only for configuration 3:2, 3:2, 4:3, 4:3 (top left), where, although the eccentricity of Jupiter was excited during the instability, it was quickly damped afterwards (e₅₅ ∼ 0.0075), due to the strong dynamical friction provided by the planetesimal disk. Only the configuration 3:2, 3:2, 2:1, 3:2 (bottom right) fulfilled all constraints together, including our new criterion E.

Summarizing, when the new constraints on the evolution of Neptune from Nesvorný (2015a,2015b); Nesvorný & Vokrouhlický (2016) are considered, the most probable initial configuration of the giant planets was 3:2, 3:2, 2:1, 3:2. The configurations 3:2, 3:2, 2:1, 2:1; 3:2, 3:2, 3:2, 3:2; or 3:2, 3:2, 4:3, 4;3 also appear possible, but less probable. This seems to be a consequence of the more relaxed orbital separation between the first and second ice giants. Therefore, one would expect that any configuration having the first and the second ice giant in distance from the Sun in 2:1 MMR or similarly relaxed could perhaps also work. However, for simplicity, we decide to restrict ourselves only to the cases of Nesvorný & Morbidelli (2012) for which we know the rate of success in fulfilling criteria A, B, C, and D.

The analysis performed in this section does not take into account the timing of the planetary instability. With the location of the inner edge of the disk that we adopted, all evolutions are relatively fast and lead to instabilities after a few 10s of Myr. In the next section we analyze whether it is possible to get a late instability out of the 3:2, 3:2, 2:1, 3:2 configuration, under some assumption on the geometry of the disk, still satisfying criterion E on the evolution of Neptune.

Since Tsiganis et al. (2005) presented the Nice model, the existence of LHB of the terrestrial planets has been associated with the possibility that the giant-planet instability occurred after hundreds of Myr of evolution (Gomes et al.2005; Bottke et al.2012). Actually, the Nice model, both in its initial version and from generic initial resonant configurations (Levison et al.2011), is compatible with both early and late instabilities, depending on the distance of the inner edge of the original trans-Neptunian disk and other disk properties (mass, excitation, etc.). Thus, the Nice model does not prove the existence of a late instability. The timing of the instability can only be determined by looking at other constraints, such as the lunar impact chronology. Unfortunately, the latter is still under debate.

However, it is not obvious that the Nice model can still be compatible with a late instability, if we fix a specific initial resonant configuration and impose a specific evolution of Neptune. As explained in the introduction, if the disk is sufficiently close to Neptune at the beginning, it is relatively easy for Neptune to migrate several au away before the instability happens, but in this case the full dynamical evolution is rapid, as in the previous section. If the initial location of the disk is sufficiently far from Neptune, the instability may happen late, but Neptune is unlikely to migrate significantly before the instability happens. Moreover, in the case of the most relaxed resonant configurations (e.g., 3:2, 3:2, 2:1, 2:1, which has Neptune initially around ∼24.5 au), the requirement that the inner edge of the disk is several au beyond Neptune while the outer edge is at ∼30 au reduces the disk to a ring. This seems implausible, particularly when considering the constraint that the disk should have hosted 1000–4000 Pluto-size bodies (Nesvorný & Vokrouhlický2016), which would probably be problematic for the late instability as well, because the spreading of the narrow disk would not permit a delay.

If one could show that the most developed and constrained version of the Nice model is incompatible with a late instability, this would be very important in clarifying the chronology of events in the solar system. Thus, in this section, we try hard to figure out whether the instability could have happened late and under which conditions. We let the reader judge whether these conditions are realistic or plausible, with our considerations summarized in Section4.

For this goal, we proceed in steps. First, in Section3.1, we determine two critical values for the distance of the inner edge of the disk from Neptune. One,δ_mig, is the maximum distance that allows Neptune to have long-range planetesimal-driven migration away from the other planets, as in the previous section. The other one,δ_stab, is the minimum distance of the disk for which the five-planet system remains stable in resonance at all times. Clearlyδ_stab > δ_mig andδ_mig is larger than or equal to the Neptune–disk distances considered in the previous section. The problem is that bothδ_stab andδ_mig depend on disk resolution (number of planetesimals and individual masses that are used to model the disk) and that realistic values cannot be used in practical simulations. So, we need to estimate trends to evaluate the values of the critical distances.

In a second step, in Section3.2, we study how a flux of dust from the planetesimal disk, due to collisional planetesimal grinding and Poynting–Robertson drag, influences planet migration. This depends on total mass flux and individual dust size. The idea is that, under an appropriate dust flux, Neptune may migrate away from its initial location, so that its distance from the planetesimal disk decreases from >δ_stab to <δ_mig. At this point, planetesimal-driven migration starts, as in the previous section. If the dust-driven migration is slow enough, the planet instability can occur late, even though it occurs only a few millions of years after the beginning of planetesimal-driven migration.

This concept is then proven in Section3.3, where we combine in a single simulation dust-driven migration and planetesimal-driven migration, and we obtain an instability after 300 Myr, with Neptune’s evolution satisfying criterion E. This simulation is of course performed with an unrealistically low disk resolution, but is a proof of concept that, for whatever disk resolution, an appropriate dust flux can be found.

We finally discuss what a realistic dust flux could be in Section3.4.

3.1. Disk Resolution versus Stability

The interaction of a planetary system with a planetesimal disk can be very resolution dependent, i.e., sensitive to the number of planetesimals used to describe the disk in the simulation and to their individual masses.

In some cases, as in those of Section2, this problem of sensitivity to disk resolution is not large because the distance between the outermost planet and the inner edge of the disk is small, so that there are from the very beginning many planetesimals on unstable orbits that encounter the planet, force its migration and, consequently, bring many new planetesimals to planet-crossing orbit. This is the reason why Nesvorný & Morbidelli (2012) found similar results for test simulations conducted with 100, 300, 1000, 3000, and 10,000 planetesimals in the disk (same total disk masses), hence they decided to use only 1000 planetesimals in their final simulations. Of course, even in these cases, the disk resolution has an influence on how grainy or smooth the planet migration will be, but it does not have a major influence on the amplitude and the timescale of the migration.

On the other hand, the resolution of the planetesimal disk becomes a real issue when the inner edge of the disk is far enough from the outermost planet. In this situation, the fraction of the disk that is unstable is small, mostly limited to the volume of the strongest MMRs with the planet. The resonant planetesimals develop planet-crossing eccentricities and, with close encounters, force the planet to move a bit outwards. This moves the resonances outwards as well, destabilizing new planetesimals (those that were not in resonance at the beginning but are swept by the resonant motion). Depending on whether the mass of the new unstable planetesimals is larger or smaller than the mass of the original unstable planetesimals, the planet migration will accelerate or damp (Gomes et al.2004; Levison et al.2007). In the damped migration mode, the planet eventually stops migrating while the resonances become empty.

If the number of planetesimals is small, the planetesimals are individually more massive in order to compose a disk with the same total mass. Thus, the migration of the planet is grainier because each encounter between the planet and a planetesimal causes a larger jump in the planet’s semimajor axis. Of course, there are fewer jumps, so that the planet migrates, on average, by the same amount. But the jumpier migration of the planet causes the resonances to jump around as well, thus destabilizing more mass from the disk than in the case of a smoother planet migration with the same average range. Because more mass is transferred from the disk to the planet-crossing orbit, the migration of the planet can change, in some cases, from the damped to the accelerated mode.

With this issue in mind, we proceed in our study looking for two critical disk distances, for a series of disk resolutions. First, we look for the minimum distance (δ_stab) the inner edge of the planetesimal disk has to have from Neptune (the outermost planet for our case) in order to keep the planets stable for at least 400 Myr. Second, we search for the maximum distance (δ_mig) that Neptune can have from the inner edge of the planetesimal disk to evolve as in Section2. From these results, we will extrapolate the values ofδ_stab andδ_mig to realistic, almost infinite resolution.⁶

Once again we consider disks with fixed surface densities (scaling as), because in this case we always guarantee the same amount of mass in the planets’ MMRs, independently of the resolution of the disk (apart from stochastic sampling). The initial distribution of eccentricity and inclination of the disk was such thate_disk = 0 and (the same initial distribution for all disks considered in this work). This is the simplest configuration that represents a cold planetesimal disk emerging from the end of the dissipation of the gas nebula. Even if the disk were initially more dynamically excited, not many differences in the overall evolution are expected because the MMRs inside the disk are very efficient in exciting the disk’s eccentricities to planet-crosser values. So, in a very short timescale, an initially cold disk would appear like an initially excited one. Finally, we consider disks with resolution of 1000, 5000, 10,000, 15,000, and 20,000 planetesimals. The outer edge of the disk is fixed at 30 au. The distanceδ between the inner edge of the disk and the planet ranges from 1 to 7 au, with increments of 1 au.

Given the great demand of CPU time required to simulate high-resolution disks, we only considered the most successful initial configuration in Section2, with planets in the 3:2, 3:2, 2:1, and 3:2 MMRs. For this configuration, Neptune is initially ∼20 au. For each value ofδ, we considered six different values of Σ₀.

We observed that, with decreasing resolution or with increasing disk surface density, the inner edge of the disk has to be farther and farther away from the planet to ensure the stability of the planetary system. For instance, Table2 showsδ_stab (with a 1 au accuracy) for different values of Σ₀ in the case of the less resolved disk (1000 planetesimals). Indeed,δ_stab ranges from 4 au for the disk with the lowest density, reaching up to 7 au in the case of Σ₀ = 0.06. Unfortunately, we cannot run disks for 400 Myr with about 15,000–20,000 planetesimals. So, we can only judge from a few experiments that in increasing the disk resolution,δ_stab tends to be between 3–5 au depending on Σ₀. For example, when Σ₀ = 0.02, experiments withδ = 4 au (unstable for disks with 1000 and 5000 planetesimals) remained stable for a resolution of 10,000 planetesimals. The same resolution makes the system unstable forδ = 3 au, with Neptune starting a fast migration at around 40 Myr. But, with the same configuration, resolutions of 5000 and 1000 planetesimals started Neptune’s fast migration in about 25 Myr and 15 Myr, respectively. This trend in instability time versus disk resolution suggests that disks with resolutions of 15,000 and 20,000 planetesimals could perhaps also be stable forδ = 3 au, but they were not performed due to excessive CPU cost. However, it is clear that a system withδ = 0 has to be unstable for whatever disk resolution, so the value ofδ_stab in the limit of infinite resolution has to exist and be positive. Our experiments suggest that the limit value is about 1–2 au smaller than the values ofδ_stab in Table2, and never significantly smaller than 3 au.

Table 2. δ_stab as a Function of Σ₀ for the Less Resolved Planetesimal Disk (1000 Planetesimals)

Σ0 ()	0.01	0.02	0.03	0.04	0.05	0.06
δstab (au)	4	5	5	6	6	7
mdisk (M⊕)	10.2	17.2	25.8	28.0	35.2	32.2

Note.δ_stab: minimum distance that the inner edge of the planetesimal disk has to have from Neptune (the outermost planet for our case) in order to keep the planets stable for at least 400 Myr. Σ₀: density of the planetesimal disk, scaled as Σ(r) = Σ₀/r. The mass of the planetesimal disk,m_disk, in Earth masses, for each pair Σ₀–δ_stab for the initial configuration 3:2, 3:2, 2:1, 3:2, with Neptune at ∼20 au and the outer disk edge at 30 au also shown for reference

Download table as: ASCIITypeset image

Another consideration that we can retrieve from Table2 is that values ofδ_stab of 4–5 au preclude giant-planet configurations with large separations among the planets, like 3:2, 3:2, 2:1, 2:1 (where Neptune is initially ∼24.5 au). In fact, in this case the “disk” extension would be less than 1.5 au, which is implausible, particularly taking into account that 1000–4000 Pluto mass bodies had to exist in this disk (Nesvorný & Vokrouhlický2016). Given that in the 3:2, 3:2, 3:2, 3:2 and 3:2, 3:2, 4:3, 4:3 configurations Neptune is unlikely to migrate far enough before the instability happens (if the migration rate is slow; Section2), the 3:2, 3:2, 2:1, 3:2 configuration is left as the only one with the chance to be compatible with alate instability. This is the configuration we are considering for the experiments in this section, and is also the one with the higher rate of success in fulfilling criterion E (Table1 in Section2).

From another similar set of simulations, where we restricted ourselves to Σ₀ = 0.02 and 0.03M_⊕ au⁻², we ran cases with 1000, 5000, 10,000, 15,000, and 20,000 planetesimals, observing Neptune’s fast migration (as in Section2) for all resolutions whenδ ≤ 2 au. Forδ ≥ 3 au, the evolution of the planet becomes resolution dependent, i.e., it suffers long-range migration if the disk is under-resolved, but not if the disk is modeled with tens of thousands of planetesimals (as previously described). Thus, we defineδ_mig = 2 au as the distance that Neptune needs to be from the inner edge of the disk to behave as in Section2 simulations, independently of the resolution of the planetesimal disk.

In order to achieve a late instability of the planetary system, we now need to find a mechanism that can slowly migrate Neptune from a distanceδ from the disk that is larger thanδ_stab to withinδ_mig. For this purpose, we investigate below the effects on planet migration of a flow of dust through Neptune’s orbit, produced by the slow collisional grinding of the disk population.

3.2. Dust Flux and Its Effect

In this section, we simulate different fluxes of dust, with different mass flows and different velocities of inward migration. The latter is characterized by the parameterβ, defined by the ratio between the radiation and gravitational forces, and directly depends on the size of the dust grains (Burns et al.1979). As usual in dust simulations, we will use super-particles, i.e., bodies that have the value ofβ of a small dust grain but carry an important amount of massm_dust, like a collection of many real dust grains.

We started our analyses by considering fixed values ofβ = 0.1, 0.05, and 0.01, corresponding with dust grain radii of about 2.3, 4.5, and 17.0μm, respectively (Moro-Martin2013, p. 431). Values ofβ > 0.1 will not be considered because those values are close enough to the limit of ejectionβ = 0.5. So, we expect their effects on the planets to be negligible. On the other hand, particles withβ < 0.01 are also not considered because those grains move very slowly. Thus, simulating a significant mass flux of these particles would require integrating the evolution of a huge number of super-particles or increasing their mass too much. So, these simulations become prohibitive.

For these experiments, we modeled the Poynting–Robertson drag in themfo_user subroutine from theMercurypackage (Chambers1999), so the dust grains could feel this drag.

The dust grains were created in the code by randomly selecting massive parent planetesimals, one at a time, in the disk and turning part of its mass into a collection of dust super-particles. In other words, the mass of each selected planetesimal is decreased as, wherek represents the number of dust super-particles created andm_dust is their individual mass. The same planetesimal can be selected more than once, provided that. The parent planetesimal keeps its original position and velocity, and the dust grains are generated with the velocity of the parent body and just a smallx–y random displacement to avoid them having subsequent identical evolutions. In fact, chaotic dynamics ensures that the super-particles will have independent evolutions for virtually any small relative displacement. This procedure is adopted in order to avoid discontinuities in the time distribution of the masses, positions, and velocities in the disk, which would occur if we created dust at random locations and/or without changing the planetesimal’s mass.

Following Table2, the planetesimal disk was considered far enough from Neptune (), with density and mass such that, without the generation of dust, the planets would not migrate and would be stable. In this way we can appreciate the effects of the dust on the dynamics of the planets.

We carried the experiments for three different values ofβ (0.1, 0.05, and 0.01) also considering different dust fluxes (ϕ_dust = 1, 10, and 20M_⊕ per 400 Myr, hereafter,, and, respectively). A collection ofk dust grains was created every 10,000 years of integration, from a different planetesimal. Thus, 2.5 × 10⁻⁵, 2.5 × 10⁻⁴, and 5.0 × 10⁻⁴M_⊕ of dust were created every 10,000 years, respectively, for the three values ofϕ_dust we investigate. The mass of a single dust grain, then, is defined by the value ofk. After several preliminary short time simulations (t < 50 Myr) consideringk = 1, ..., 40 (which means 100, ..., 4000 particles per Myr), we saw virtually the same evolution of Neptune. This is expected because even the largest “dust grain” (5.0 × 10⁻⁴ M_⊕ fork = 1 and) is still very small compared to Neptune’s mass (∼15M_⊕). What is important for the planet dynamics is not the size of each dust grain, but instead, itsβ and the considered mass flux. So, for the main simulations we adoptedk = 1.

What was observed (Figure2) is a first phase of a few Myr during which the planet does not migrate; a second phase of fast migration, independent ofβ and mostly related to the distant encounters between the fifth and fourth planets after they have gone out of resonance, which brings the planet a few tenths of au away from its initial position in ∼20 Myr; then a third phase of slow, constant migration. This last phase is the one driven by the interaction with the dust flow. We observe that forβ = 0.1, the planet basically does not migrate (dark gray in Figure2). Forβ = 0.05 (black) and 0.01 (light gray), we clearly see outward migration, with a rate linearly anticorrelated withβ.

Zoom InZoom OutReset image size

Next, we ran several additional cases for dust fluxes,, and, and assigned a random value ofβ between 0.1 and 0.01 with a uniform distribution for the dust grains created. We think that this approach is valid and somewhat more realistic to represent a size distribution of dust grains.

These last runs resulted in an evolution very similar to those usingβ = 0.05 (Figure2 dashed black). This is because the existence of dust grains withβ close to 0.1 (ineffective in driving migration) is compensated by the existence of dust grains withβ ∼ 0.01 (more effective).

Adopting this prescription forβ, Figure3 shows the evolution of Neptune’s semimajor axis under the effects of the dust flux (ϕ_dust) as a function of the time. As in Figure2, one can immediately notice three changes in the migration regime: during the first ∼10 Myr Neptune does not migrate; then it migrates very rapidly up to ∼25 My; finally it migrates more slowly and linearly with time, with a rate proportional to the dust-mass flux as expected.

Zoom InZoom OutReset image size

In fact, in this third stage, Neptune’s dust-driven evolution enters a steady-state regime. Each dust particle is first captured in MMR with Neptune, pushing the planet inwards, then is released and, through close encounters, moves the planet outwards. The net effect of the dust is a slow outward migration of the planet, much slower than that obtained from an equivalent mass of planetesimals encountering the planet until their ultimate dynamical removal. This is because, at any given time, while some dust particles are having close encounters with Neptune pushing it outwards, some other particles are still in MMR with Neptune pushing it inwards (recall that planetesimals do not push the planet inwards, even if they are captured in MMR because they do not feel the Poynting–Robertson drag). Moreover, the lifetime of dust in the planet-encountering regime is shorter than for planetesimals (dust can drift across the planet’s orbit without close encounters because of its radial migration).

In this phase, the migration is linear in time and the migration rate depends linearly on the dust-mass flux. We fit the averaged migration of Neptune with the formula, with and, 0.002, and 0.004 au Myr⁻¹ for (dashed black line), (white line), and (solid black line), respectively. From the linear relationship of on, one can easily estimate the migration of Neptune for any dust-mass flux.

3.3. Delayed Evolution

In this section we combine the results of Sections3.1 and3.2 to obtain a late planet instability starting from the configuration (3:2, 3:2, 2:1, 3:2), while fulfilling all constraints on the planetary evolutions (particularly criterion E presented in Section2).

The idea is as follows. Under the effect of dust flux, Neptune leaves its MMR with the adjacent planet in about 10 Myr, and its orbital semimajor axis rapidly moves outwards by ∼0.4 au. If at this point the planet is still more thanδ_mig away from the inner edge of the disk, the evolution of the planet will be dominated by dust-driven migration (with one caveat discussed below). If the dust flux has the appropriate magnitude, Neptune can reach the distanceδ_mig from the disk after hundreds of Myr. When this point is reached, rapid planetesimal-driven migration will bring the planet to its final position again in ∼10 Myr, and during this rapid migration the planetary instability will occur, as in Section2.

The caveat that we anticipated above is that when the distance of the planet from the disk is larger thanδ_mig but less thanδ_stab, some planetesimals start to interact closely with the planet. Thus, the migration of the planet is no longer driven only by dust, but by a combination of dust and planetesimal effects. In the absence of dust, the planetesimal-driven migration would rapidly damp out, because the distance of the planet from the disk is larger thanδ_mig and because of the very definition of this parameter (the largest distance from which planetesimal-driven migration is not in a damped mode). But, because the dust imposes linear migration on the planet, planetesimals keep escaping from the disk through the moving MMRs with Neptune and add their contribution to the migration rate.

To evaluate this contribution, we proceeded with a series of simulations. We considered a disk with an inner edge at 25 au, Σ₀ = 0.03M_⊕ au⁻², andm_disk = 25.8M_⊕, modeled with 1000, 5000, and 10,000 planetesimals. We imposed a linear outward migration to Neptune with the rates found in the previous section, and we measured the actual migration rate of the planet, resulting from the combination of the migration rate we impose and of close encounters with planetesimals escaping from the disk.

With these series of nine simulations (three resolutions and three values of), we found that the actual migration rate is about 2.8 times the imposed rate, whatever the value of and the disk’s resolution. Obviously, the value of 2.8 is likely linear on the surface density of the disk, but should also depend on the actual distance of the planet from the disk (i.e., which MMRs sweep the disk). But because the disk parameters used in these simulations are close to those of the disk that we will use below, we will adopt the factor 2.8 to set up our ultimate simulation.

Our ultimate simulation will be a proof of concept that it is possible to obtain an instability delayed by about 300–400 Myr from the (3:2, 3:2, 2:1, 3:2) initial planetary configuration. The simulation accounts for the combined effects of dust and planetesimal scattering, with no imposed migration. We consider this simulation a “proof of concept” because a constant dust-mass flux is applied (a simplistic assumption; see next section) and the disk is quite under-resolved (1000 planetesimals). But, given that a non-zero value ofδ_stab exists for whatever disk resolution, it is clear that with enough computational resources, a similar evolution could be obtained for a more realistic disk resolution, rescaling the dust-mass flux in an appropriate way (see Table3). Also, according to Izidoro et al. (2015), if Uranus and Neptune formed from a set of proto-cores migrating inwards and stuck in resonance with Jupiter and Saturn, some gap has to be expected.

Table 3. Dust Flux (ϕ_dust) as a Function of the Initial Separationδ_stab Between Neptune and the Disk, for the Initial Configuration (3:2, 3:2, 2:1, 3:2), Necessary to Bring Neptune from 20.5 au toδ_mig ≤ 2 au in ∼380 Myr when Considering Disks with Enough Resolution

δstab (au)	4	5	6	7
ϕdust (M⊕/400 Myr)	7.05	11.5	16.45	21.15

Using the information obtained in the previous sections, one can compute the dust flux needed in order to achieve an instability at a desired (late) time. For instance, let us suppose a disk (arbitrarily) withM_⊕ au⁻² and an inner edge 5 au beyond the initial position of the planet, i.e., beyond the valueδ_stab for a high-resolution disk. We know that, after about 10 Myr of evolution, Neptune will be at about 20.4–20.6 au, namely, ∼4.5 au away from the disk. We also know from Section3.1 that, for the most resolved diskδ_mig ≤ 2 au, and that after Neptune reaches this distance, it takes about 10 Myr to the instability point. Thus, we need the planet to migrate 2.5 au (i.e., from 4.5 au to 2 au from the disk), in say, 380 Myr. This implies Myr⁻¹ and, from the simulations above and the results of Section3.2, Myr⁻¹ (where/2.8) andM_⊕ over 400 Myr. The values ofϕ_dust for other initial planet–disk separations, and considering the value ofδ_mig for the large disk resolution limit, are reported in Table3. Only separations ≥4 au are considered because those ≤3 au would lead to an early migration.

In our proof of concept simulation, however, we can only deal with a low resolution disk due to computational limitations. We use a disk resolution of 1000 planetesimals and we assume. According to Table2, in this caseδ_stab = 7 au. So, we place the inner disk’s edge at this distance from Neptune. Thus, repeating the calculations above for this distance between Neptune and the disk, we findM_⊕/400 Myr. The values ofk andβ were chosen ask = 1 andβ uniformly distributed and randomly selected in the range between 0.1 and 0.01. Finally, to reduce the effects of the coarse resolution of the disk on Neptune’s migration, every time a planetesimal becomes Neptune’s crosser, it is replaced by 10 new planetesimals each with one-tenth of the mass of the parent planetesimal. This results in a disk that when in contact with Neptune could have a maximum of 10,000 planetesimals.

The results of the “proof of concept” simulation are shown in Figures4 and5. Figure4 shows the evolution of the semimajor axis (a), the perihelion distance (q), and the aphelion distance (Q) for all planets as a function of time. Figure5 shows snapshots of (a versuse) and (a versusi) of the whole system for different times of the evolution.⁷

Zoom InZoom OutReset image size

Zoom InZoom OutReset image size

The evolution shown in these figures has an instability after about 300 Myr, and the evolution of the planets is in considerably good agreement with that deduced in Nesvorný (2015a,2015b) from Kuiper Belt constraints. In particular, the instability happens when Neptune is at 28 au as required. Moreover, as we can see from these figures, the evolution of Neptune before the onset of instability follows closely the assumptions/predictions of Tables2 and3, and the dust mechanism of Section3.2. The instability did not happen precisely at the time we targeted, due to the approximations inherent in our linear combination of different processes, which in reality have a nonlinear interplay. Nevertheless, the instability happens late enough to suggest that the association between a planetary instability and the trigger of the LHB is indeed possible.

3.4. Dust Generation and Lifetime

It is very likely that at the disappearance of the protoplanetary disk of gas, the giant planets of the solar system were all in MMRs with each other, because this is the typical outcome of gas-driven planet migration (Morbidelli et al.2007). Nesvorný & Morbidelli (2012) investigated which resonant configurations are more likely to reproduce the current orbital configuration of the giant planets via a phase of dynamical instability, such as that described in the Nice model. They considered four criteria of success and found a number of initial configurations that meet all of them in a non-negligible fraction of the simulations. All these configurations require the existence of a fifth giant planet, with a mass comparable to those of Uranus or Neptune, which was eventually ejected during the instability.

More recently, Nesvorný (2015a,2015b) showed that explaining and reproducing the current orbital structure of the Kuiper Belt requires that Neptune had a specific evolution, reaching a heliocentric distance of ∼28 au by planetesimal-driven migration before the instability happened.

Thus, in Section2, of this paper, we have revisited the work of Nesvorný & Morbidelli (2012), adding a new criterion of success from the results of Nesvorný (2015a,2015b). Of the configurations selected by Nesvorný & Morbidelli (2012), the one denoted 3:2, 3:2, 2:1, 3:2 (each integer ratio referring to the MMR between a pair of adjacent planets) stands out as the most likely to fulfill the new criterion of success. Other configurations (the 3:2, 3:2, 2:1, 2:1, the 3:2, 3:2, 3:2, 3:2, and the 3:2, 3:2, 4:3, 4:3) cannot be excluded, but their success rate relative to the new criterion is significantly lower.

An open issue with the Nice model is the timing of the giant-planet instability. If one wants to explain the LHB of the Moon as a spike in the impact crater rate (Tera et al.1974), the instability has to have happened late (Gomes et al.2005; Levison et al.2011; Bottke et al.2012), unlike in the simulations of Section2, in which the instability occurred after a few 10s of Myr. It is not obvious that a five-planet resonant configuration can remain stable for a long time and become unstable then, particularly with the specific evolution of Neptune such as that required to sculpt the Kuiper Belt. Thus, in Section3 we investigated whether a late instability from the 3:2, 3:2, 2:1, 3:2 configuration is possible and under which conditions. First, we looked for characteristics of the planetesimal disk that can maintain the planets in stable orbits for 400 Myr. We found that the evolution and stability of the planets are very sensitive to the resolution of the planetesimal disk. Thus, we determined the minimum distance that the inner edge of the planetesimal disk has to be away from Neptune to keep the planets in stable orbit for 400 Myr as a function of the surface density of the disk and its numerical resolution. We named this distanceδ_stab (Table2). In a second step, we found a characteristic Neptune–disk distance, which we namedδ_mig, from where Neptune starts long-range planetesimal-driven migration. Thus, for the system to exhibit a late instability together with an evolution of Neptune like the one that is needed to sculpt the Kuiper Belt, Neptune needs to be extracted from the 3:2, 3:2, 2:1, 3:2 resonance chain and to migrate outwards over hundreds of Myr until it reaches the distanceδ_mig. To understand how this could be possible, we introduced the concept of dust-driven migration, i.e., migration due to the gravitational interaction of a planet with the dust produced in a distant planetesimal disk, spiraling inward by Poynting–Robertson drag. We have shown that a late instability can be achieved when considering a flux of dust (ϕ_dust) ranging from 7 to 22 M_⊕ per 400 Myr, depending on the value ofδ_stab (Table3).

Continuing in Section3.4, we evaluated how much dust could be produced from the disk-like configurations that we assumed and found good agreement with the dust flux described above. Indeed, we found that the total amount of dust produced in 400 Myr may vary from ∼7M_⊕ to a maximum of about 32M_⊕ depending on the SFD and catastrophic disruption energies assumed for the planetesimal disk and its constituents. Finally, we also found that the dust production rate expected for a self-grinding disk is not constant, but it decays as 1/t. The abundance of dust, which is also limited by the dust collisional lifetime, should decay as. Thus, because the migration of the planet is linear in the amount of dust that it encounters, the evolution of its semimajor axis should be proportional to, which implies that small changes in the flux of dust can result in Neptune reaching theδ_mig distance either early or very late. Therefore, obtaining a late instability around 400 Myr requires a more delicate tuning of the dust production.

Summarizing, based on the work of Nesvorný & Morbidelli (2012) and considering the evolution of Neptune required to explain the structure of the Kuiper Belt (Nesvorný2015a,2015b), our results indicate that the preferred initial configuration of the giant planets at the time of the dissipation of the gas nebula is the 3:2, 3:2, 2:1, 3:2 configuration, and that, although an early instability reproduces more easily the evolution of Neptune required to explain the structure of the Kuiper Belt, such evolution is also compatible with a late instability.

We thank an anonymous reviewer for useful comments and suggestions on the submitted manuscript. R.D. acknowledges support provided by grants #2015/18682-6 and #2014/02013-5, São Paulo Research Foundation (FAPESP) and CAPES. A.M. acknowledges support by the French ANR, project number ANR-13-13-BS05-0003-01 project MOJO (Modeling the Origin of JOvian planets). R.S.G. acknowledges his grant No. 307009/2014-9 from CNPq, Conselho Nacional de Desenvolvimento Científico e Tecnológico, Brazil. The work of D.N. was supported by the NASA Emerging Worlds program.