Dynamical Evolution of Closely Packed Multiple Planetary Systems Subject to Atmospheric Mass Loss

A gap in exoplanets ’ radius distribution has been widely attributed to the photoevaporation threshold of their progenitors ’ gaseous envelope. Giant impacts can also lead to substantial mass loss. The out ﬂ owing gas endures tidal torque from the planets and their host stars. Alongside the planet – star tidal and magnetic interaction, this effect leads to planets ’ orbital evolution. In multiple super-Earth systems, especially in those that are closely spaced and / or contain planets locked in mean motion resonances, modest mass loss can lead to dynamical instabilities. In order to place some constraints on the extent of planets ’ mass loss, we study the evolution of a series of idealized systems of multiple planets with equal masses and a general scaled separation. We consider mass loss from one or more planets either in the conservative limit or with angular momentum loss from the system. We show that the stable preservation of idealized multiple planetary systems requires either a wide initial separation or a modest upper limit in the amount of mass loss. This constraint is stringent for the multiple planetary systems in compact and resonant chains. Perturbation due to either impulsive giant impacts between super-Earths or greater than a few percent mass loss can lead to dynamical instabilities.

The locations of innermost planets in multiple super-Earth systems are often in the proximity of their host stars.These short-period planets are exposed to the intense flux of X-ray, FUV, and UV photons, especially during their host stars' infancy (Baraffe et al. 2005;Murray-Clay et al. 2009;Owen & Jackson 2012;Erkaev et al. 2016;Kubyshkina et al. 2018;Lalitha et al. 2018;Ionov et al. 2018;King & Wheatley 2021).Photoionization of H/He-dominant planetary envelopes can lead to a substantial amount of mass-loss from low-mass close-in planets, while only a small fraction of the gaseous envelope may be removed from the more massive and distant planets (Ionov et al. 2018).Theoretical analyses suggest the evaporation rate is a sensitive function of planets' core and envelope masses.An estimation on the evaporation threshold Wang & Lin leads to the prediction of a mass-loss dichotomy and a bimodal distribution in the apparent sizes among planets over a small range of masses (Wu & Lithwick 2013;Owen & Wu 2013, 2017;Lopez & Fortney 2013).
The success of the theoretical prediction introduces a new conundrum concerning the progenitors of the super-Earths on either side of the period gap.If the main cause of the gap is due to a critical condition for photo-evaporation, the progenitors of all planets on either side of the gap must have formed with nearly identical core masses (M c ∼ 10M ⊕ ) and a modest fraction (f ∼ a few %) of envelope masses (Madhusudhan 2019;Owen 2019;Wu 2019;Zhang 2020;Rogers & Owen 2021;Rogers et al. 2023).
The near uniformity of M c provides both clues and constraints on the core accretion scenario (Ida & Lin 2004) which assumes an initial emergence of planetesimals through runaway growth (Dullemond & Dominik 2005;Kenyon & Bromley 2009;Garaud et al. 2013), gravitational instability (Goldreich & Ward 1973;Weidenschilling & Cuzzi 1993;Youdin & Shu 2002;Garaud & Lin 2007), streaming instability (SI) (Youdin & Goodman 2005) or vortices trapping (Johansen & Youdin 2007) of small dust in protostellar disks.Thereafter, these building blocks rapidly coagulate until a few oligarchic embryos emerge (with masses from a fraction to a few M ⊕ ) which consumes all of the residual planetesimals in their feeding zones (Kokubo & Ida 1998).However, embryos' growth can be further supplied by fast migrating optimum-sized grains (i.e., approximately millimeterto meter-sized pebbles; (Ormel & Klahr 2010)), until they attain the pebble-isolation mass, typically a few to 10M ⊕ and evolve into protoplanetary cores.Subsequently, cores induce local maxima in the surface density and pressure distribution for the disk gas which suppresses their growth (Lambrechts & Johansen 2012;Lambrechts et al. 2014;Ormel 2017;Bitsch et al. 2018).This effect quenches the cores'growth via further pebble accretion.The gas accretion rate is regulated by both opacity and entropy advection (Ikoma et al. 2000;Rafikov 2006; Lee et al. 2014;Piso & Youdin 2014;Piso et al. 2015;Lee & Chiang 2015, 2016).The accumulation of pebbles near their tidal barrier also increases their collision frequency with fragmentary byproducts in the form of sub-mm grains which are well coupled to and dragged by the gas flow.This process elevates the effective dust opacity, reduces the radiative flux in the envelope, limits the envelope's asymptotic masses (M e ) through gas accretion, and preserves super-Earths in the stellar proximity at ∼ 0.1 AU (Chen et al. 2020).
These considerations provide a natural explanation for the nearly uniformity in M c of close-in super-Earths in the context of in-situ formation scenario (Chatterjee & Tan 2014;Liu & Ormel 2017).However, it is not clear whether there may be a dispersion in the metallicity and mass of their envelopes (M e ), especially in multiple planetary systems where the gas and grain supplies may vary among members with a range of semi-major axis.Moreover, some multiple systems are locked in mean motion resonances (MMRs), including chain configurations which suggest the likelihood of their orbital migration (Gillon et al. 2017;Leleu et al. 2021).Similar to hot Jupiters (Lin et al. 1996), some super-Earths may have formed at a few AU from their host stars and migrated over large distances.At these distances, the cores' critical masses for inducing pebble isolation and suppressing gas accretion may differ from those close to their host stars (Chen et al. 2020).Such diversity may also introduce dispersion in the asymptotic values of both M c and M e among infant super-Earths.
Rapid convergent migration can also lead to orbit crossing and giant impacts.Even if super-Earths' progenitors had a substantial gaseous envelope, such a population could be the byproducts of giant impacts (Liu et al. 2015) during the depletion of their natal disks.These events can lead to increases in their M c and losses in their envelope masses.
In this paper, we consider a potential implication of substantial mass-loss from super-Earths' gaseous envelope.This possibility would arise naturally if some super-Earths formed with predominantly gaseous envelopes or contained mostly volatile ices and subsequently relocated to the proximity of their host stars to Figure 1.The distribution of planetary systems which contain more than three planets observed by TESS mission.The grey circles mean the innermost planet in the system, and the inner one in a planet pair is labeled in red if the planet pair is in or near 2:1 MMR, red sold dots represent the period ratio deviations of 2.0 by no more than 2%, while the red circles represent the period ratio deviations of 2.0 by no more than 10%.The green sold dots and circles show the period ratio results of 3:2 MMR, and the blue sold dots show the results of 4:3 MMR.The grey sold dots mean no obvious sign of first order MMR relationship between two adjacent planets.
endure intense stellar irradiation.In principle, photoevaporation would remove a large fraction of their initial mass in their envelope if their core mass was below some threshold value for atmospheric retention.But, the outflowing gas is subjected to the tidal torque and the wind from their host stars (Carroll-Nellenback et al. 2017) which induce angular momentum transfer between the evaporated gas and the orbits of the planets.This interaction leads to fractional changes in the semi-major axes of the planets through an amount of mass they have lost.With a set of idealized models, we first demonstrate (in §2) that this effect can lead to dynamical instabilities in closely packed multiple planetary systems.For computational simplicity, we construct these idealized models with a set uniform planetary masses, scaled locations and separation (in units of Roche radii).
Although idealized models are useful to illustrate the destabilization tendency induced by the mass-loss, they need to be generalized to more realistic multiple planetary systems with a wide range of observed kinematic configurations.The list of known multiple planetary systems is likely to grow rapidly with the ongoing confirmation campaign for hundreds of planetary candidates in the Kepler dataset (Batalha et al. 2013;Mazeh et al. 2013;Fabrycky et al. 2014;Dressing et al. 2017).In addition, TESS mission has already found several multi-planet systems and many more are anticipated (Quinn et al. 2019), as shown in Figure 1 which includs candidate and confirmed systems.In some of the closely packed systems, modest orbital migration induced by the loss of envelope's mass or impulsive energy and angular momentum changes due to giant impacts may introduce non-negligible perturbation and disturb the stability of these systems especially if they are closely packed, similar to the TRAPPIST-1 system (Gillon et al. 2017), or TOI-178 system (Leleu et al. 2021).

Wang & Lin
Some compact multiple planetary systems reside in nearly commensurable orbits with their members' period ratios close to some low-order integer fractions.Presumably, these systems have maintained their present-day configurations over the age of their host stars.Nevertheless, the resonant planets exchange orbital energy and angular momentum which lead to variations in their semi-major axes, eccentricities, and libration in their relative orbital phases.In §3, we scout the island of stability of compact systems with multiple MMRs.Based on a migration model, we robustly reproduce a series of multiple planetary systems with MMRs.With these initial orbits, we show that these MMR configurations occupy narrowly confined islands of stability.We impose small impulsive changes in the spacing between adjacent planets in these systems and show that they can lead to rapid evolution in their orbital elements, including transitions to dynamical instability.We then introduce protracted mass-loss and orbital evolution to explore some quantitative constraints on the amount and rate of photo-evaporation based on its consequential change to the stability of the systems.Finally, we summarize our results and discuss their implications in §4.

STABILITY OF SYSTEMS WITH MASS-LOSS OF THE INNERMOST PLANETS
In this section, we illustrate the consequence of massloss on an idealized system of multiple planets with equal initial mass and scaled spacing, in units of their local Roche radius.This set of initial conditions are designed to highlight the destabilizing effect of mass-loss on compact multiple planetary systems.We assume envelope is lost after the depletion of these planets' natal disk and introduce a simple prescription to approximate its influence on the orbit of the mass-losing planet.

The timescales of mass-loss and orbital migration
Herein, we approximate the loss of the planet's envelope proceeds on a constant, characteristic time scale τ m such that where ṁp is the mass-loss rate of planet.The planet's orbital angular momentum L is such that If the evaporation of the planet proceeds after they are tidally synchronized, the planet's internal tidal dissipation may also be efficiently intense to render a non negligible ė.In this case, photo-ionization occurs mostly on the day side facing the host star.Around quiet mature stars (those without an intense stellar wind or/and strong magnetic field), gas leaves the planets' envelope via their inner Lagrangian point and mostly falls toward the star (Carroll-Nellenback et al. 2017).In this case, angular momentum of the systems is essentially conserved (i.e.L = 0) and the planets' orbits expand on a time scale τ a = a p / ȧp = τ m /2.
Planet-host stars' UV flux is most intense during their infancy when they also blow powerful wind which may exert strong ram pressure on the evaporated gas from the planets' surface.Planetary magnetic field may also modify the outflow direction of the evaporated gas (Bisikalo et al. 2013;Zhilkin & Bisikalo 2019).Under some circumstances, angular momentum may be lost from the systems (Carroll-Nellenback et al. 2017), while the planets' orbital decay on a similar time scale τ a (Tang et al in preparation).
In order to take into account of these diverse possible outcomes, we introduce a simple idealized prescription to characterize planets' mass-loss with Main model parameters are: 1) the planets' initial mass (m 0 ), 2) the fractional mass-loss amplitude A, and 3) the mass-loss time scale τ m .The model parameters used here are chosen for illustrative purposes and they can easily be scaled to the appropriate values applicable to real observed systems.We approximate the planets' orbital evolution with where we adopt f p = −1 for orbital expansion associated with conservative evaporation and f p > 0 for the possibility of orbital decay induced by mass and angular momentum loss from the systems.We incorporate the effect of planets' mass-loss by introducing a torque in their equation of motion where F θ = F θ θ is being applied at the radius r i in the azimuthal (θ) direction with In the above equation, J = GM * a p is the planets' specific angular momentum at the guiding center of their epicycles.Note that F θ vanishes for t >> τ m .Two additional forces F damp + F migI are included in the equation of motion to take into account the planets' interaction with their natal disks during their infancy in section 3 (see §3.2), in this section we only consider the effect of F θ .
For the numerical simulations to be presented below, Equation ( 6) is implemented into the Rebound Code (Rein & Tamayo 2015).The orbits of all planets are assumed to be coplanar around a M * = 1M ⊙ host star.At the onset of the calculation, we assume the planets' orbits to be circular with randomly distributed mean anomaly and argument of pericenter.Each model is run for a maximum time scale of 10 7 yrs or until the planets' orbits cross each other.

Orbit crossing in multiple planetary systems
For illustrative purpose, we consider an idealized system of nine equal-mass (M p ) planets.At the onset of the calculation, they are equally separated in logarithmic interval.The normalization is set, such that the semi-major axis of the third innermost planet is designated to be at 1AU and that of the i + 1 th planet is scaled to be at where R H = (2µ/3) 1/3 (a i + a i+1 )/2 is the Hill radius, µ = M p /M * , and k is a constant model parameter which represents the planetary separation in units of mutual Hill radii.One of the planets in the system is designated to lose a fraction A/(1 + A) of its initial mass due to the photo-evaporation and to undergo either inward or outward migration over a fraction ∼ −2f p A/(1 + A) distance.
In these equal-mass systems, the innermost planet is the most likely planet to endure photo-evaporation.For the first set of simulations, we set A = 0.05 and assume conservative photo-evaporation with f p = −1 (i.e.planets undergo outward migration as they loss mass).We consider three groups of simulations with M p = 1, 3, and 10 M ⊕ , respectively.For each group, 12 models are computed with different values of k (initial separation) and mass-loss time scale (τ m ).
We list the initial parameters in Table 1.Stable models (planets avoided orbit crossing for at least 10 7 yr) are denoted in red color whereas models which led to orbit crossing for any adjacent planets within 10 7 yr are marked in black color.The stabilities of the systems in Table 1 are also shown in the upper panels of Figure 3.The low-mass (1M ⊕ ) group have relatively small R H such that some compact systems (with modest k) may have overlapping MMRs.Only systems with large initial k(≥ 16) and long mass-loss time scale τ m ≥ 10 6 remain stable for more than 10 7 yrs.
In comparison with the M p = 1M ⊕ planets, the magnitude of R H is larger for the M p = 3M ⊕ planets.Adjacent planets with k ≥ 16 are stable because they avoid most MMR's, regardless the mass-loss time scale τ m .In Figure 2, we show the evolution of a typical stable (for > 10 7 yrs) system (with M p = 3M ⊕ , A = 0.05 and τ m = 10 6 yrs).As a consequence of innermost planet's outward migration, it is captured into a 4:3 MMR with the next planet.For the high mass (M p = 10M ⊕ ), R H is sufficiently large, all systems with k ≥ 10 are stable for more than 10 7 yrs.In general, systems with larger M p , longer τ m , larger initial k values are more stable.
These results are generalized to cases with nonconservative mass-loss.We consider the possibility that it may lead to inward migration with f p = 1.Under this assumption, planetary migration would be divergent and the systems would remain stable if the masslosing planet has the smallest semi-major axis.But, if the mass-losing planet is one of the more distant member of a closely-spaced multiple planetary system, nonconservative mass-loss would lead to convergent migration and the onset of dynamical instability.We have also examined the possibility that the mass-losing planet is initially located in the middle of the pack.In this case, either inward or outward migration can destabilize the system in accordance with the results presented in Table 1.
The gaseous envelope contributes ∼ 10% of the total mass in Uranus and Neptune.We consider a second series of models with larger fractional mass-loss (A = 0.1 and 0.2) from the planets with M p = 10M ⊕ , f p = −1 (outward migration), a range of values for k and τ m (Table 2).These models indicate that greater fractional mass-loss (ie larger A) leads to more extended migration and forces the planetary system to evolve into more compact configuration and become unstable, as shown in the bottom panels of Figure 3. Modest mass-loss (A = 0.2) can destabilize the systems even in the limit that they started with a wide initial separation (k = 16) regardless the masses of the planets and the mass-loss time scale.
Table 1.The initial parameters of the cases in section 2.2.A is the mass-loss fraction of the total planet (A=0.05 for the cases in this table), Mass represents the initial mass of planet in the system, k is the initial planetary separation in units of mutual Hill radii, τm is the mass-loss timescale.The cases shown in red color are stable for at least 10 7 yrs, while the cases in black represent the systems with orbital crossing happened before the end of simulations.The results in the previous section clearly show the influence of mass-loss on the stability of multiple planetary systems.However, the equal-mass systems we have adopted as initial conditions are highly idealized.Emerging multiple planets around the same host stars are likely to adjust their initial separation due to their interaction with their natal disk and each other.Their differential migration may lead them to be captured into each other's MMRs.

Case
In this section, we consider the impact of migration induced by either impulsive giant impacts or gradual photo-evaporation on some compact multiple planetary systems with MMRs.For initial conditions, we first construct a generic formation model for these systems in §3.1.Applying modest amount of uniform compression or expansion, we carry out numerical experiments to show that although equilibrium configuration may emerge after the disk depletion, their islands of stability are narrowly confined.Due to this effect, minor impulses and/or protracted loss of a modest amount of mass can indeed lead to dynamical instabilities for these systems.Figure 2. The evolution of a typical stable system.In this run, the masses of planets in the system are 3 M⊕, the initial separation between planets is 14 R hill , and the timescale of mass-loss is 10 6 yrs.The system is stable for more than 10 Myrs.The innermost planet loses 5% of its total mass in this process.Panel (a) shows the evolution of semi-major axis and mass of the innermost planet, he red line shows the evolution of mass and the black line represents the evolution of semi-major axes.Panel (b) represents the semi-major axis evolution of nine planets in the system, the black lines mean the evolution of the semi-major axis and the grey lines display the evolution of the pericenter and apocenter.Panel (c) displays the evolution of the relative space kji between two adjacent planets Pi and Pj (here, i represents the serial number of the planet from inner to outer, 1 i 8, j = i + 1), different colors display different kji as shown in the legend in panel (c).
Five out of the six planets in the TESS mission found system TOI-178 are in near chain resonances as shown in Figure 1 (Leleu et al. 2021).
Through exhaustive numerical integration, the present-day masses and orbital elements of these systems have been determined, under the constraint that they are stable for at least 10 7 yrs (Tamayo et al. 2017).These calculations do not, however, take into account the possibility of any mass-loss from one or more of their members, even though some of these planets have observed radii smaller than the lower bound of the gap in the super-Earths' radius distribution (presumably they do not have any significant atmosphere, possibly as a result of previous photo-evaporative mass-loss).Due to the vast range of potential initial orbital configurations and mass-loss efficiency, it would not be practical to separately examine the effects of mass-loss on all the individual observed systems.Instead, we adopt some generic formation scenarios for MMR systems and utilize the results generated from these models as the initial conditions for the systems' subsequent impulsive or prolonged orbital evolution.
Planets in observed chain MMRs or near MMRs systems have masses several times that of the Earth.Their orbital semi-major axes are typically an order of magnitude larger than the radii of their host stars.Insitu formation of these super-Earths requires disk surface density in planetesimals to be many orders of magnitude larger than that extrapolated from the minimum mass nebula scenario (Ida & Lin 2004;Wang et al. 2012;Chiang & Laughlin 2013;Tamayo et al. 2017).Although high-eccentricity migration due to the planet-planet scattering (Rasio & Ford 1996;Lin & Ida 1997;Beauge & Nesvorny 2012), Kozai resonance (Wu & Murray 2003;Anderson et al. 2016), and the secular chaos (Nagasawa et al. 2008;Wu & Lithwick 2011;Hamers et al. 2017) can lead to the relocation of close-in Jupiter-mass gas giants, the configurations of MMRs in multiple super-Earth systems are difficult to accomplish (Giacalone et al. 2017;Xue et al. 2017) in these circumstances.

Formation of multiple planets with MMRs
The most plausible scenario for forming super Earth planetary system with MMRs is orbital migration due to their tidal interaction with their natal disks (Goldreich & Tremaine 1980;Lin et al. 1996;Bryden et al. 2000;Masset & Snellgrove 2001;Migaszewski 2016;Izidoro, et al 2017;Liu et al. 2017;Liu & Ormel 2017;Ramos et al. 2017).These planets' masses are too small to clear deep gaps in the gas surface density distribution of their natal disks.Nevertheless, The cases with A=0.2 and A=0.1 are corresponding to that listed in table 2. The results shown in the upper and lower panels of the rightmost column are the same cases.
they exert a tidal torque on the disk which leads to planets' type I migration (Ward 1997;Tanaka et al. 2002).
The direction and pace of embedded planets' migration is determined by the total torque where γ = 1.4 is the adiabatic index, h is the aspect ratio, and Ω p is the Keplerian frequency at the location of the planet.We adopt the empirical minimum mass solar nebular model (MMSN; Hayashi 1981) for Σ p = Σ g (r p ) with where f g is the enhancement factor.Herein, we consider the gas depletion time scale τ dep in the range of [0.1, 3] Myrs (Haisch et al. 2001).
The total-torque coefficient f tot = f CR + f LB is the sum of corotation (f CR ) and Lindblad (f LB ) torque on the planets (Paardekooper et al. 2010(Paardekooper et al. , 2011)).The magnitude and sign of f CR and f LB are functions of the evolving distribution of s = ∂lnΣ g /∂lnr, β = ∂lnT /∂lnr, viscosity ν and thermal diffusivity ξ.For the mass range of super-Earths considered here, the corotation torque is unsaturated and maximally effective.It can lead to outward type I migration in the viscously heated inner disk region and inward migration in the irradiated outer disk region (Kretke & Lin 2012).Consequently, the migrating super-Earths have a tendency to converge toward the transition radius r t which separates these two domains.
In order to take into account these possibilities, we introduce a prescription to characterize this transition in f CR across r t such that where we consider a range for the torque coefficient coef (2, 5, or 10).While r t ∼ a few AU in a MMSN, it can reduce to a fraction of 1 AU during the disk depletion phase (Garaud & Lin 2007;Wang et al. 2021).For our simulations, we choose two sets of r t = 0.5 and 3 AU to represent the range of possibilities.
In most regions of the disk, the Lindblad torque is negative (f LB < 0).But, near their inner edge r m , disks are sharply truncated by the magnetosphere of their host stars and a steep positive surface density gradient leads to super Keplerian azimuthal velocity.In this region the Lindblad torque becomes positive (f LB > 0) such that inward migration of super-Earths may be stalled.In order to take into account this possibility, we introduce another prescription for Tidal torque due to disk-planet interaction leads to type I migration on a time scale where τ 0 = M * (h 2 γ/2µΩ p )/(Σ p r 2 p ) (Tanaka et al. 2002;Kley & Nelson 2012).The associated time scale for eccentricity damping is τ e = e/ ė = h 2 τ a .These contributions are incorporated into the Equation of motion (Eq 6) through Γ tot θ M p r and ( 14) Aside from F migI and F damp , we do not take into account the self gravity of the disk gas nor its effect on the procession of the planets' orbits.While the former contribution is negligible, the latter effect may lead to secular resonances among well separated planets in sparsely populated systems (Zheng et al. 2017).

Emergence of stable multiple planetary systems with MMRs
In order to obtain some reasonable realistic initial conditions, we adopt a generic formation scenario for multiple planetary systems with MMRs based on the assumption that their kinematic orders were established through differential and convergent type I migration (Liu et al. 2017;Liu & Ormel 2017).Similar to the previous section, we adopt systems of nine equal mass (M p = 10 −5 M ⊙ ) planets designated to be P1 to P9 in order to increasing a around a M * = 1M ⊙ (Zhou et al. 2007).At the onset of our calculation, we place these planets with equal fractional spacing kR H around the third (Earth-like) planet at 1AU.
In order to examine how planets' asymptotic semimajor axis and eccentricity distribution may depend on the evolution of disk properties (Wang & Ji 2017), we consider three groups of models (G1, G2, G3) with several different values of f g (= 0.1, 0.2, 0.5, 1) and a range of τ dep (= 0.1, 1/3, 1, 3 Myr).For the planets, we adopt the same value for coef(= 10) and r m (= 0.05 AU) for two initial sets of k(= 8, 16) and r t (= 0.5, 3 AU) (Table 3).These model parameters lead to three representative groups of multi-planet systems with diverse MMRs configurations: Group 1: All the planets start with a i > r t and initial separation of 8 R H .In the absent of gas, their orbits would be stable for at least 10 7 yrs (Zhou et al. 2007).In a low-mass (f g = 0.1) thin (h = 0.05(r/1AU) 0.25 ) disk with r t (= 0.5 AU), all planets undergo type I inward migration (panel (a), Figure 4).The light grey lines represent the planets' apo and peri center distances.The innermost planet first stalls at r t .Others follow and they force it to a smaller radius.Convergent migration leads to asymptotic separation in the range of 8-14.2R H . Several planets capture each other into a chain of discrete (3:2, 4:3, and 5:4) MMRs (panel b) with integer period ratios and small eccentricity (e ≤ 0.035).Planets P5, P6, and P7 form a 5:4:3 resonant chain and P7, P8, and P9 form a 4:3:2 MMRs.Therefore, according to the estimation of crossing time (Zhou et al. 2007), such configuration formed under inward migration with relative separation larger than 8 R H will be stable for more than 10 7 yrs.
Group 2: Widely separated (16R H ) planets are installed on either sides of r t (= 3 AU) in a low-mass (f g = 0.1) disk.In a gas-free environment, their orbit crossing time is considerably > 10 7 yrs (Zhou et al. 2007).In their natal disks, planets 7, 8, and 9 undergo inward type I migration, while other planets migrate outward.Panels (a) and (b) in Figure 5 show the evolution of a i and the separation between adjacent planets.After 10 5 yrs, the inner six planet pairs are captured into 3:2 MMRs, and planet 7 and 8 are in 4:3 MMR.P9 is also captured into a 3:2 MMR at about 2 × 10 5 yrs with P8 firstly, and into a 4:3 MMR at ∼ 5 × 10 5 yrs.Planets P8, P7, and P6 are also trapped in a 4:3:2 chain MMRs.All the resonant angles librate with small amplitude after the planets are trapped into MMRs (Panel c) with Table 3. Parameters used in the simulations, Group 1, Group 2, and Group 3. k0 = (ai+1 − ai)/RH represents the relative separation between two adjacent planets initially.rt is the transition radius, τ dep is the depletion timescale of the gas disk, coef is the coefficient in the corotation torque, h=H/r represents the disk's aspect ratio, and fg is the enhancement factor of the gas density.the exception of that between P7 and P8 (∼ 83 • ).Their eccentricities are ∼ 0.1 (see Table 4).Different from the results in Group 1, the relative separation between two adjacent decreases with their evolution.The multiple planetary system may become unstable after the inward and outward migration.
Group 3: The initial conditions are similar to that in Group 2, but with different f g (= 0.2 − 1) and depletion timescale (0.1-3 Myr) of the gas disk.In a minimum mass nebula (Model G3-1, panel (a) of Figure 6), the inner three planet pairs enter into 3:2 MMRs, P5 and P6 into 4:3 MMR, two pairs (P6 and P7, P7 and P8) into 5:4 MMRs.P8 and P9 into 6:5 MMR in less than 0.1 Myr.This system becomes unstable after ∼ 5 Myrs.In a lessmassive disk (Model G3-2 with f = 0.5 in panel b), the planets are captured into similar MMR's in ∼ 0.1 Myrs to that in G3-1.In the least-mass disk (Model G3-3 with f g = 0.2 in panel c), planet pairs enter into 3:2 and 4:3 MMRs with wider spacing between them.Thus, with a less dense gas disk, the configuration formed through orbital migration can be stable for longer time with larger relative separations (Zhou et al. 2007).In comparison with Model G3-3, we consider τ dep = 0.1, 1, 3 Myrs in models G3-4, G3-5, and G3-6.With a longer depletion timescale, the asymptotic configuration is more stable with smaller libration amplitude of their resonant angles (Figure 7).
The above results show the robust emergence of multiple planetary systems with MMRs with different disk properties, and reveal the relationship between the stability of configuration formed through orbital migration and the disk properties.Here, we utilize these results as initial conditions to study the consequence of mass-loss and its associated orbital evolution.In this section, we adopt an idealized prescription to demonstrate the limited extent of stable region in the proximity of a set of representative MMRs from the results of Model G2 (Figure 5).First, we uniformly redistribute the planets' azimuthal position without any change to their semi-major axes.Since the eccentricities of all the planets are very small, these changes correspond to the randomization of their resonant angles.For some pairs, this modification breaks their MMR since the modified resonant angles are outside their original range of libration.Nevertheless, these systems remain stable without orbit crossing over a time scale of 10 Myr.

Island of stability near the MMRs
In order to test the extent of these islands of stability, we introduce an artificial, impulsive, uniform, fractional, expansion or compression for all the planets in these systems.For computational simplicity, we multiply the separation between adjacent pairs with a con-stant factor of f k < 1 or f k > 1 respectively.We fix the semi-major axis of the middle planet (P5) and study the implications of artificial incremental changes with f k in the range between 0.8 and 1.25.This approach, though idealized, provides suggestive indication on the limited extent of islands of global stability in closely packed multiple systems.
The evolution of these systems is examined based on the numerical solution of Equation ( 6) with F θ , F damp , and F migI are set to be zero.The time scale for the first set of planet to cross each other's orbit is plotted in Figure 8.In general, a relatively small global contraction (with f k slightly less than unity) can reduce the crossing time scale to less than 1 Myr.But, the stability of this system is preserved (on time scale longer than 10 Myrs) with either a very limited changes (1 ≤ f k ≤ 1.01) or a large (f k ≥ 1.2) expansion factor, although some intermediate (1.01 < f k < 1.1) expansion factor can lead to relative short (< 2 Myrs) orbit crossing time.
We generalize these calculations for systems which contain diverse MMR chains.These configurations are formed in disks with a range of τ dep (= 0.1, 1/3, 1, and 10 Myr) for the same value of f g = 0.2 (models G3-3, G3-4, G3-5 and G3-6 in Table 3).Although their libration amplitudes for the adjacent resonant pairs differ considerably, their evolution subsequent to the incremental changes is very similar to that shown in Figure 8.For Model G3-3, the libration amplitude is relatively large (Fig. 7) such that very small changes in the planets' spacing can lead to the breaking of the MMRs.Subsequent secular interaction between these multiple planets leads them to orbital crossing on time scales less than 1 Myr.With a relatively large τ dep , Model G3-6 produces a more closely packed system.Its small libration amplitudes suggest the system contains deep and stable MMR chains.Nevertheless, relative small incremental changes in their spacing would shorten the crossing time to less than 1 Myr.These results imply that although multiple MMRs can be established in compact planetary systems while they are embedded in their natal disks, the extent of these islands of stability becomes narrowly confined after the depletion of the disk gas, 10% compression (f k ≤ 0.9) of its original stable configuration can lead to quick orbital crossing, the crossing time shorter than 10 4 yrs.

Preservation of MMRs during an extended epoch of mass-loss
In §3.3 we show that as a consequence of tidal interaction with their natal disks and type I migration, planets in multiple planetary systems with MMR's robustly emerge.In disks which deplete rapidly, these Panel (a) shows the evolution of semi-major axis, pericenter and apocenter of nine planets.The black lines display the evolution of semi-major axis, and the grey lines show the evolution of pericenter and apocenter of nine planets.Panel (b) shows the evolution of relative space between adjacent planet kji, the meanings of kji are the same as in Figure 4. Panel (c) shows the evolution of resonant angles θ = (p + q)λ ′ − pλ − q̟ ′ , where the unprimed and primed parameters refer to the inner and outer planets respectively, two planets are in or near (p+q):p resonance, λ and ̟ are the mean longitude and longitude of pericentre.Amji represents the resonant angle amplitude of the planet pair Pi and Pj (i represents the serial number of the planet from inner to outer, 1 i 8, j = i + 1).
systems have large libration amplitude whereas in more massive disks, they tend to be more compact.Under most circumstances, their MMR configuration can be preserved for at least 10 Myrs.But their islands of stability are surrounded by nearby unstable regions where planets' mutual perturbation can trigger dynamical instability and orbit crossing within a few Myrs (see §3.4).Small compressive or expansive impulsive adjustments in the distance between the adjacent resonant pairs are more prone to dynamical instability than those between the spaciously separated multiple planetary systems (see §2.2).We now consider the stability and dynamical evolution of these systems subjected to individual planet's atmospheric mass-loss and orbital adjustment (see §2.1).
Because of their confined islands of stability, the destiny of multiple planetary systems with MMRs may provide more stringent constraints on amount of mass-loss induced by the photo-evaporation of their atmosphere.
This process is expected to be protracted and persistent.In order to simulate such process, we adopt the asymptotic (at 10 Myrs) kinematic structure of an emerging multiple planetary system (from model G2) as the initial conditions for later dynamical evolution.Six typical cases are shown in Table 5.The parameters are adopted for the fractional amount (A) and time scale (τ m ) of mass-loss from the innermost planet.In all six cases, conservative mass-loss is applied to the same initial orbital elements of the nine planets (Table 4).Am21 Am32 Am43 Am54 Am65 Am76 Am87 Am98 .The amplitude of resonant angles in the system with different gas disk depletion timescales.Amji (1 i 8, j = i + 1) has the same definitions as in Figure 5. Red circles, blue asterisks and black squares represent depletion timescales of 10 5 , 10 6 and 3 × 10 6 yrs respectively.The distribution of crossing time vs the separation change of planet pairs.The original system which is obtained from G2 can be stable for more than 10 7 yrs as shown in Figure 5 .In order to find out the stable region for such multiple planetary system, we introduce an artificial factor f k acting on the original stable system to change the separation between adjacent planets.We fix the semi-major axis of the middle planet P5 and change the separation between adjacent by f k , 0.8 ≤ f k ≤ 1.25.f k < 1 refers to a global contraction comparing to the original system configuration and f k > 1 represents an expansion to the original stable system.
Table 5.The parameters of mass-loss in the cases in section 3.2.A is the mass-loss fraction, τm is the mass-loss timescale, θ/2π is the average libration amplitude of all the planet pairs, and ē is the average eccentricity of all planets in the system.In case 1-3, we consider the innermost planet lost 5% of its mass (A = 0.05) on timescales ranging from 10 6 − 10 4 yrs.Figures 9 and 10 show the evolution of Case 1 and 3, respectively.Case 1 remains stable for 10 Myrs despite the large eccentricity (0.24) for planet 2. The resonant angles of all the MMR pairs remain close to π.Their small average libration amplitudes (Figure 9) suggest that this system is locked into a deep resonant chain.The mass-loss time scale is sufficiently long to enable the planets to adiabatically adjust their orbits and to share among themselves the angular momentum The black line means the evolution of the semi-major axis and the red one represents the evolution of mass.In this case, the innermost planet will lose about 0.14 Earth-mass which is about 5% of the initial mass.Panel (b) shows the evolution of semi-major axis, percenter and apocenter.The black lines mean the evolution of the semi-major axis and the grey lines display the evolution of the pericenter and apocenter.At the end of the simulation, the eccentricities of planets will be excited.But the system is still stable due to the relative angles between them.Panel (c) shows the evolution of the resonant angles θji of two adjacent planet Pi and Pj, 1 i 8, j = i + 1. θji = (p + q)λ ′ − pλ − q̟ ′ , where the unprimed and primed parameters refer to the inner and outer planets respectively, two planets are in or near (p+q):p resonance, λ and ̟ are the mean longitude and longitude of pericentre.Panel (d) shows the evolution of eccentricities of planets in the system epi, pi represents the serial number of the planet from inner to outer, 1 pi 9, the eccentricities of different planets are labeled in different colors, as shown in the legend in panel (d).The evolution with A=0.1 and τm = 10 6 yrs in Case 2 is similar to that in Case 1.

Case
change due to the mass-loss from the innermost planet.Their MMR configuration is maintained.Although the magnitude of τ m in Case 2 is smaller than that in Case 1, the kinematic structure of the system is also maintained for more than 10 Myrs.In Case 3, (with τ dep = 10 4 yrs) the planets' eccentricities e and libration amplitude θ of the resonant angles between adjacent MMR pairs grow substantially (Figure 10).The planets' average eccentricities ē and average libration amplitude θ in Case 3 are much larger than those in Case 1 (Table 5).Although this system remains marginally stable after 10 Myrs, chaotic relaxation is likely to trigger dynamical instability and orbit crossing over a longer time scale.
In case 4-6, the fractional mass-loss is set to be 10% (i.e.A = 0.1, higher than that in case 1-3) on time scales τ m = 10 6 − 10 4 yrs.The growth of planets' eccentricities and resonant angles between adjacent MMR pairs in Case 4 and 5, shown in Table 5, are similar to that in Case 3. In Case 6, τ m = 10 4 yrs is comparable to the libration time scales for the resonant angles between the adjacent pairs (in the range of a few 10 4 yrs).The rapid change of the innermost planet's orbit can no longer be accommodated by adiabatic adjustment of other planets such that the libration amplitude of resonant angles ( θ) grow to ±π for some relevant pairs to break their MMRs (Liu et al. 2015), as shown in Figure 11.This transition from libration to circulation is followed by orbit crossing.Later close encounters excite large eccentricities.Thus, with a relatively small (∼ 5% − 10%) fractional massloss, compact systems with multiple MMRs have a ten- shows the evolution of a, a(1-e), and a(1+e) of all the planets.At the end of the simulation, the results are similar to that in Case 3. Panel (c) shows the evolution of the resonant angles.The resonant angles tend to be librate in a large range than in Case 1. Panel (d) shows the evolution of eccentricities of planets in the system.The meanings of different colors are the same as in Figure 9 .
dency to become unstable and undergo orbit crossing, especially for the system with higher mass-loss fraction or shorter mass-loss timescales.
For the systems with nine planets around the central star, the system tends to be destroyed by an orbital crossing event.However, if the number of planets in the system decreases, the results change.Up to now, the observed multiple planetary systems with more than three planets in the system are mainly focused on no more than five.Herein, we test the influence of mass-loss process in the system with three or five planets.The initial configurations are chosen from the innermost part of the results of G2 as shown in Table 4. Figure 12 and 13 show two typical results.The trends of these groups with fewer planets in the systems are similar to that in nine-planet systems.With shorter mass-loss timescale or higher mass-loss fraction, the orbital crossing events are easier to happen in the system.As shown in Figure 12, when the fraction of mass-loss reaches 15% of the innermost planet's total mass, the amplitude of the resonant angles increase with the evolution, planet pairs in the system try to escape from the MMRs which is similar to that happened in Case 3 as shown in Figure 10.But the critical fraction of mass-loss (A c ) that leads to the orbital crossing varies.With the decreasing of the number of planets in the system, A c increases.If the system contains five planets, A c increases to 25%, while when the innermost planet loses almost 35% of its total mass, orbital crossing happened in three-planet system.Additionally, after the orbital crossing happened, planetary systems which contains less than five planets are able to survive.Figure 13 shows the results with A =0.25 and τ m = 10 4 yrs in a five-planet system.After the orbital crossing, the inner three planets collide to be a larger one, and the outer two planets also collide and merge to a larger planet.At the end of the simulation, there are two super-Earths (with mass ratio m outer /m inner = 1.4) formed in the system with widely separated configura- .
tion (with period ratio of ∼ 3.6), like the configuration of system Kepler-1090 which contains two widely separated planets (Morton et al. 2016;Valizadegan et al. 2022), the ratio of planet radius of these two planets is ∼ 1.46 and the period ratio is ∼ 4.68.Therefore, there are no clear indications of significant orbital adjustment induced by atmospheric losses in observed closely packed multiple planetary systems.Planets in such systems probably emerged from their natal disks without much primary atmosphere.While for the widely separated planetary systems, the planets may preserve a certain amount of atmosphere once before they suffered from strong photo-evaporation rising from the central star.

SUMMARY AND DISCUSSIONS
In this work, we investigate the dynamical evolution of multiple terrestrial planetary systems under the influence of orbital migration and mass-loss of planets caused by photo-evaporation.
We showed that sparsely populated multiple systems (not initially in MMRs) can maintain their stability despite some modest fractional (∼ 10%) mass-loss especially if it proceeds on a relatively long time scale or high planetary mass or larger relative space between the nearby planets.Nevertheless, it is difficult for a system to preserve its stability if one or more planets lose 20% of their total mass.Therefore, for the systems with planets not in MMRs, the stability of the system determined by the intensity of mass-loss process on the planets, the lowest mass of the planets, and the distance between the adjacent planets.
We adopt the assumption that super-Earths undergo extensive type I migration while they are embedded in their natal disks.They emerge with diverse kinematic properties which are determined by the surface density distribution and the depletion rate of the gas disk.Some of these systems are spatially extended while others contain MMRs.The depletion timescale of gas disk will determine whether the system can maintain stability af- ter they are captured into MMRs.According to the observation data, the depletion timescale of gas disk is few million years (Haisch et al. 2001;Williams & Cieza 2017).Based on our results, the planetary system can maintain long-lived stability for τ dep > 1 Myr which results in deep MMRs of planet pairs.The surface density of gas disk determines how compact the final configuration which also related to the stability of the systems.In general, such systems are dynamically stable with relative small eccentricities and libration amplitudes.
Another factor influence the stability of the system during the orbital migration process is the location of the transition radius from the optically thick region to the optically thin region in a disk.If the planets formed in a disk with the transition radius very close to the central star, planets tend to undergo only inward migration and the system may be easy to be stable if there is a stopping mechanism.If the transition radius is far away from the central star, planets have large opportunity to form a compact configuration, just like the configuration shown in Group 2 after both inward and outward migration.However, the island of stability of systems after the orbital migration process is very small.Tiny perturbations will destroy the system.Densely packed systems with resonant chains can break their mean motion resonances by sufficiently large (≥ 10%) impulsive semi-major axis changes or persistent atmosphere loss by one or more of the planets.
However, if the system formed around a star with strong irradiation, planets which are in a MMR configuration may lose significant fractions of their mass.Such a process may destroy the stable configuration they have already formed.With the effect of mass-loss on the innermost planet, the final configuration is related to the timescale competition between the massloss process and the libration of the resonant angles.When the libration timescale is larger than the massloss timescale, the amplitude of the resonance angles shows the evolution of eccentricities of planets in the system.The meanings of different colors are the same as in Figure 9 .become large enough to break the resonance making orbital crossing happened in the system.With the timescale of mass-loss less than 2 × 10 4 years, planet pairs may escape from MMR.The system can be destroyed if the innermost planet lost ∼ 10 − 20% of its total mass for the system with nine planets which is consistent with the results obtained from the mass-loss process in all first-order resonances multiple planetary systems (Matsumoto & Ogihara 2020).The system, however, can survive in a widely separated configuration if they are formed from the system with less than five embryos initially.
Based on our work, the destiny of the closely packed multiple planetary systems can be classified into three possible types.(I) The system is far away from the central star with MMRs configuration.Multiple planets in such systems undergo orbital migration which make them in or near MMRs and find their stable configuration finally.(II)The innermost planet is close to the central star and formed with little atmosphere.The system can keep in stable and chain resonances configuration after little mass-loss ( less than ∼5-10% of its total mass).We find that systems similar in configuration to TRAPPIST-1 (Gillon et al. 2017;Tamayo et al. 2017) or K2-138 (Lopez et al. 2019) can be stable for more than 10 7 yrs if they are formed through orbital migration and without primary atmospheres.(III) The innermost planet is close to the central star and formed with a certain atmosphere.A stable configuration can be destroyed, leaving a single planet system or system with wide separated planets.From TESS mission results, there are more than 20 multiple planetary system candidates with at least three planets in the system have been observed as shown in Figure 1.35/59 planet pairs are in near first order MMRs, among them 20/59 are near 2:1 MMRs as labelled in red dots and circles, 17/59 are near 3:2 MMRs labelled in green dots and circles, and 3/59 are near 4:3 MMRs.A large proportion of systems contain MMRs configurations.Such systems may be stable for longer than 10 7 yr, if they are formed through orbital migration, and the innermost planets have a very small amount of atmosphere, such as the systems TOI-1136 (Dai et al. 2022), TOI-1749, TOI-175 (Fukui et al. 2021), andTOI-178 (Leleu et al. 2021), they can be classified as type II.The systems with wide separated planets not in MMRs, like TOI-880 (Sun et al. 2022), TOI-451 (Newton et al. 2021), TOI-4010 (Kunimoto et al. 2022), TOI-431 (Osborn et al. 2021), andTOI-561 (Lacedelli et al. 2020) may have gone through a similar formation process with type III planetary systems.
In this work, we mainly focus on the influence of massloss process on the dynamical evolution and stability of multiple planetary systems.More than 800 multiple planetary system have been confirmed including 67 four-planet systems, 27 five-planet systems, 10 sixplanet systems, 1 seven-planet system and 1 eight-planet system.For the systems with five or more planets in closely packed configurations, the innermost planet must not lose too much atmosphere, less than 10% of its total planetary mass, to keep the system stable.Planets in multiple planetary systems undergoing rapid photoevaporation, especially those that lose a large fraction of their mass, are prone to orbital crossing which may destroy the system or leave fewer planets in widely separated configurations after fierce collisions and merger processes.According to our results, long-term stable multiple planetary systems which formed through orbital migration from the outer disk favor a slower massloss process as concluded in the section §3.If the planets in multiple planetary system sustain mass-loss pro-cess for shorter than ∼ 2 × 10 4 yrs, the system will be faced with the challenge of being unstable or being destroyed.Therefore, for planets undergoing a corepowered mass-loss process, which acts on Gyrs instead of Myrs (Ginzburg et al. 2018;Gupta & Schlichting 2019), the multiple planetary systems may remain in stable configurations for more than 10 7 yrs.On the other hand, planets formed in-situ with minimal migration are likely to go through a "boil-off" process, with rapid atmospheric escape on timescales ∼ 10 5 yrs (Ikoma & Hori 2012;Owen & Wu 2016;Ginzburg et al. 2016).Such multiple planetary systems are likely to be unstable, leaving fewer planets in widely-separated configurations.If, however, the planets lose less than ∼ 5% of their total mass, the system may remain stable, even with a "boil-off" phase.

Figure 3 .
Figure3.The stabilities of multiple planetary systems are shown with varied mass-loss timescales and fractions.k represents the relative Hill Radius between two adjacent planets initially.The y-axis represents the logarithm of mass-loss timescale.The systems which are stable for at least 10 7 yrs are labelled in dark blue, while the systems with orbital crossing before the end of simulations are labelled in light blue.Upper three panels show the stabilities of systems with different planetary masses, but with same mass-loss fraction (A=0.05).The cases in the upper three panels are corresponding to that listed in table 1.The lower three panels show the stabilities of systems with different mass-loss fraction, but with same planetary masses (m = 10M⊕).The cases with A=0.2 and A=0.1 are corresponding to that listed in table 2. The results shown in the upper and lower panels of the rightmost column are the same cases. k0

Figure 4 .
Figure4.The results of G1.In this Group, all planets undergo inward type I migration.The transition radius is 0.5 AU.The planets are labelled in 1 to 9 from the innermost to the outermost.Panel (a) shows the evolution of semi-major axes of nine planets with black lines, the evolution of the pericenter and apocenter of nine planets are shown in grey lines.Panel (b) displays the evolution of the relative separation between planets.kji means the the relative separation between two adjacent planets Pi and Pj (here, i represents the serial number of the planet from inner to outer, 1 i 8, j = i + 1), different colors display different kji as shown in the legend in panel (b).

Figure 5 .
Figure5.The results of G2.In this Group, the transition radius locates at 3 AU.Thus planets 7, 8 and 9 undergo inward migration, while the other planets migrate outward.Planet pairs can be trapped into 3:2 (k=14.2) and 4:3 MMRs (k=10.1).Panel (a) shows the evolution of semi-major axis, pericenter and apocenter of nine planets.The black lines display the evolution of semi-major axis, and the grey lines show the evolution of pericenter and apocenter of nine planets.Panel (b) shows the evolution of relative space between adjacent planet kji, the meanings of kji are the same as in Figure4.Panel (c) shows the evolution of resonant angles θ = (p + q)λ ′ − pλ − q̟ ′ , where the unprimed and primed parameters refer to the inner and outer planets respectively, two planets are in or near (p+q):p resonance, λ and ̟ are the mean longitude and longitude of pericentre.Amji represents the resonant angle amplitude of the planet pair Pi and Pj (i represents the serial number of the planet from inner to outer, 1 i 8, j = i + 1).

Figure 6 .
Figure 6.The evolution of relative separation between adjacent planets in different cases with different gas density.Panel (a), (b), and (c) are corresponding to the cases with fg=1, 0.5, and 0.2, respectively.fg is the enhancement factor of the gas density as shown in equation (10).fg = 1 refers to the gas density profile of MMSN model.The colors of different lines represent the relative separation between different planet pairs as the legend shown in panel (a).
Figure7.The amplitude of resonant angles in the system with different gas disk depletion timescales.Amji (1 i 8, j = i + 1) has the same definitions as in Figure5.Red circles, blue asterisks and black squares represent depletion timescales of 10 5 , 10 6 and 3 × 10 6 yrs respectively.
Figure8.The distribution of crossing time vs the separation change of planet pairs.The original system which is obtained from G2 can be stable for more than 10 7 yrs as shown in Figure5.In order to find out the stable region for such multiple planetary system, we introduce an artificial factor f k acting on the original stable system to change the separation between adjacent planets.We fix the semi-major axis of the middle planet P5 and change the separation between adjacent by f k , 0.8 ≤ f k ≤ 1.25.f k < 1 refers to a global contraction comparing to the original system configuration and f k > 1 represents an expansion to the original stable system.

Figure 9 .
Figure9.The evolution of the system in Case 1, A=0.05 and τm = 10 6 yrs.Panel (a) shows the evolution of the innermost planet.The black line means the evolution of the semi-major axis and the red one represents the evolution of mass.In this case, the innermost planet will lose about 0.14 Earth-mass which is about 5% of the initial mass.Panel (b) shows the evolution of semi-major axis, percenter and apocenter.The black lines mean the evolution of the semi-major axis and the grey lines display the evolution of the pericenter and apocenter.At the end of the simulation, the eccentricities of planets will be excited.But the system is still stable due to the relative angles between them.Panel (c) shows the evolution of the resonant angles θji of two adjacent planet Pi and Pj, 1 i 8, j = i + 1. θji = (p + q)λ ′ − pλ − q̟ ′ , where the unprimed and primed parameters refer to the inner and outer planets respectively, two planets are in or near (p+q):p resonance, λ and ̟ are the mean longitude and longitude of pericentre.Panel (d) shows the evolution of eccentricities of planets in the system epi, pi represents the serial number of the planet from inner to outer, 1 pi 9, the eccentricities of different planets are labeled in different colors, as shown in the legend in panel (d).The evolution with A=0.1 and τm = 10 6 yrs in Case 2 is similar to that in Case 1.

Figure 10 .
Figure10.The evolution of the system in Case 3, A=0.05 and τm = 10 4 yrs.Panel (a) shows the evolution of the innermost planet.In this case, the innermost planet will lose about 0.14 Earth-mass which is about 5% of the initial mass.Panel (b) shows the evolution of a, a(1-e), and a(1+e) of all the planets.At the end of the simulation, the results are similar to that in Case 3. Panel (c) shows the evolution of the resonant angles.The resonant angles tend to be librate in a large range than in Case 1.Panel (d)  shows the evolution of eccentricities of planets in the system.The meanings of different colors are the same as in Figure9

Figure 11 .
Figure11.The evolution of the system in Case 6, A=0.1 and τm = 10 4 yrs.The system become unstable at about 0.35 Myrs.Panel (a) shows the evolution of the innermost planet.In this case, the innermost planet will lose about 0.3 Earth-mass which is about 10% of the initial mass.Panel (b) shows the evolution of a, a(1-e), and a(1+e) of all the planets.Panel (c) shows the evolution of the resonant angles.The resonant angles tend to be librate in a large range from 0 to 2π.Panel (d)  shows the evolution of eccentricities of planets in the system.The meanings of different colors are the same as in Figure9

Figure 12 .
Figure12.The evolution of the system with five planets, A=0.15 and τm = 10 4 yrs.Panel (a) shows the evolution of the innermost planet.In this case, the innermost planet will lose about 0.5 Earth-mass which is about 15% of the initial mass.Panel (b) shows the evolution of a, a(1-e), and a(1+e) of all the planets.At the last stage of the simulation, the planet pairs in the system tend to escape from the MMRs.Panel (c) shows the evolution of the resonant angles.The resonant angles tend to be librate in a large range at the end of simulation.Panel (d) shows the evolution of eccentricities of planets in the system.The meanings of different colors are the same as in Figure9.

Figure 13 .
Figure13.The evolution of the system with five planets, A=0.25 and τm = 10 4 yrs.The system become unstable at about 0.06 Myrs.Panel (a) shows the evolution of the innermost planet.After the system become unstable, the inner three planets collide into one larger planet which is about 9.33 M⊕.And then planet 4 and 5 collide into a planet with 6.67 M⊕.Panel (b) shows the evolution of a, a(1-e), and a(1+e) of all the planets.Panel (c) shows the evolution of the resonant angles.Panel (d) shows the evolution of eccentricities of planets in the system.The meanings of different colors are the same as in Figure9

Table 2 .
The initial parameters and results in the cases with high A (A=0.1) in section 2.2.The definitions of parameters listed in the table are the same as that in table 1.The cases shown in red color are stable for at least 10 7 yrs, while the cases in black represent the systems with orbital crossing happened before the end of simulations.

Table 4 .
The orbital parameters of the planets at the end of the simulation in the case in Group 2, which also represent the initial parameters of planets in the cases in section 3.2.The planets are named with their serial number from the inner to outer with a, e, i, Ω, ω, f , and Mass representing the semi-major axis, eccentricity, inclination, longitude of the ascending node, argument of periastron, true anomaly and mass of planets, respectively.