Dynamical Instability in Multi-Orbiter Systems with Gas Friction

Closely-packed multi-planet systems are known to experience dynamical instability if the spacings between the planets are too small. Such instability can be tempered by the frictional forces acting on the planets from gaseous discs. A similar situation applies to stellar-mass black holes embedded in AGN discs around supermassive black holes. In this paper, we use $N$-body integrations to evaluate how the frictional damping of orbital eccentricity affects the growth of dynamical instability for a wide range of planetary spacing and planet-to-star mass ratios. We find that the stability of a system depends on the damping timescale $\tau$ relative to the zero-friction instability growth timescale $t_{\rm inst}$. In a two-planet system, the frictional damping can stabilise the dynamical evolution if $t_{\rm inst}\gtrsim\tau$. With three planets, $t_{\rm inst} \gtrsim 10\tau - 100\tau$ is needed for stabilisation. When the separations between the planetary orbits are sufficiently small, $t_{\rm inst}$ can be less than the synodic period between the planets, which makes frictional stabilisation unlikely to occur. As the orbital spacing increases, the instability timescale tends to grow exponentially on average, but it can vary by a few orders of magnitude depending on the initial orbital phases of the planets. In general, the stable region (at large orbital spacings) and unstable region (at small orbital spacings) are separated by a transition zone, in which the (in)stability of the system is not guaranteed. We also devise a linear map to analyse the dynamical instability of the"planet + test-mass"system, and we find qualitatively similar results to the $N$-body simulations.

A widely-used criterion to evaluate the stability of a compact planetary architecture is to compare the distance between neighbouring planets to the mutual Hill radius, which scales with  1/3 , where  is the mass ratio between the planets and the central star.For two-planet systems with initially circular orbits, a "Hill instability" criterion for the critical planetary semi-major axis separation can be derived using the conservation of the Jacobi constant and a shearing box approach of the trajectories near the conjunctions (Henon & Petit 1986;Gladman 1993).On the other hand, an alternative approach focusing on the growth of orbital chaos due to resonance overlap leads to a critical separation proportional to  2/7 (Wisdom 1980;Duncan et al. 1989).The exact value of the exponent is still under debate even for the restricted three-body problem (i.e., one of the planets has a negligible mass), due to the chaotic nature of the system and the blurriness of the stability-instability divide (Petrovich 2015;Petit et al. 2020).
For a system comprising more than two planets, stability is never guaranteed, but the timescale for instability to arise depends on the mutual distance between the planets (e.g., Chambers et al. 1996;Smith & Lissauer 2009;Lissauer & Gavino 2021).Meanmotion resonances further prevents a clear determination of the stability boundary of multi-planet systems, leading to the recent development of alternative approaches, such as the use of machinelearning algorithms as a predictive tool (Tamayo et al. 2020).
Closely-packed, unstable multi-planet systems can be a natural product of planet formation and migration in protoplanetary discs.Planet-disc interaction typically damps the planet's eccentricity, thus prevents orbital crossings between planets and suppresses the dynamical instability.Several previous works have explored the interactions between planets embedded in gaseous discs in various scenarios of planet growth and migration, showing that a variety of planetary architecture can be produced (e.g. Lee et al. 2009;Marzari et al. 2010;Matsumura et al. 2010;Lega et al. 2013;Liu et al. 2022).However, a quantitative assessment of the effects of gaseous discs on the dynamical instability of multi-planet systems is lacking.The main goal of this paper is to systematically evaluate how the strength of planet eccentricity damping due to the disc affects the onset of the dynamical instability.To this end, we apply parameterized frictional forces acting on the planets, and use body simulations to quantify the onset and growth of instability as a function of the spacing between planets.
Although the problem studied in this paper pertains to planetary systems, it is also important for understanding the evolution of stellar-mass black holes (sBHs) embedded in AGN discs (e.g.McKernan et al. 2012;Bartos et al. 2017;Stone et al. 2017;Secunda et al. 2019;Tagawa et al. 2020;Li et al. 2022a).In particular, AGN discs can help bringing sBHs circulating around a supermassive BH into close orbits due to the differential migrations of the BHs and migration traps (Bellovary et al. 2016).Close encounters between such tightly-packed sBHs may lead to the formation and merger of binary BHs (Li et al. 2022a; see also Secunda et al. 2019;Tagawa et al. 2020).The effects of the AGN disc on the formation and evolution of the binary BHs are uncertain, and generally require hydrodynamical simulations for proper understanding (see Baruteau et al. 2011;Li et al. 2021bLi et al. , 2022b;;Li & Lai 2022;Dempsey et al. 2022).In Li et al. (2022a), we incorporated parameterised weak frictional forces (with the eccentricity damping time 10 5 times the orbital period) in the long-term -body integrations of multi-BH systems around a supermassive BH, and we found that these forces did not lead to enhanced formation of BH binaries.However, we did not explore how the initial instability growth is affected for a wide range of damping timescales.
In this paper, to systematically study the onset and growth of instability in multi-orbiter (planets or BHs) systems around a central massive body (star or supermassive BH) with and without frictional forces, we consider a wide range of "planet" to "star" mass ratios1 , from  = 10 −7 to 10 −3 .As noted above, because of the effect of mean-motion resonances, the stability property of multi-orbit systems does not just depend on the initial orbital spacings in units of the Hill radius (∝  1/3 ).
The rest of this paper is organised as follows.In Section 2, we study the restricted three-body problem (star, planet and test particle), and characterise the defining features of the instability as a function of , with and without damping force.In Section 3, we present an analytical model (algebraic map) based on the shearingbox approximation, which provides insights to the numerical results.In Section 4, we examine the multi-planet case and its relation to the restricted problem.We summarise our findings in Section 5.

INSTABILITY OF RESTRICTED THREE-BODY PROBLEM WITH FRICTION
In this section, we consider a co-planar system of a central star with mass , an inner planet with  1 =  on a circular orbit ( 1 = 0), and an outer test particle with  2 = 0.The test particle may experience a frictional force that tends to damp its eccentricity.
Using -body integration, we numerically determine the stability boundary of such a system for various values of  and the damping time.

Setup of the simulations
In our numerical simulations, the initial orbital separation between the planet and the test mass is set as where is the mutual Hill radius, and  is a dimensionless constant.We henceforth use the initial orbital period of the planet,  1 , as the unit of time.We consider different combinations of  and  and carry out 500 runs for each combination.The initial eccentricity of the test particle is  2 = 10 −5 , and the initial values of the argument of the periapsis, the longitude of the ascending node, and the mean anomaly are sampled randomly from the range [0, 2], assuming they all have uniform distributions.The systems are simulated using -body software REBOUND (Rein & Liu 2012) and the IAS15 integrator (Rein & Spiegel 2015).For each run, when the orbits of  1 and  2 overlap within their mutual Hill radius, i.e. when is satisfied by their real-time orbital elements, this system is considered unstable.If such instability is found, we stop the simulation immediately and register the time as  inst .Otherwise, the simulation ends when it reaches  = 10 5  1 and we consider the system "stable".

Systems with no friction
We first study the no-friction situation.With  = 10 −7 , 10 −6 , 10 −5 , 10 −4 , 10 −3 and initial  = 2.0, 2.1, 2.2, ..., 4.0, the whole experiment contains 105 suites of simulations (and 500 runs with randomised initial angles in each suite).Fig. 1 shows the fraction of stable runs in each suite of simulations.In all five panels (which present five different mass ratio ), there exist a "grey zone" of  between unstable systems for small 's and stable systems for large 's.The size of the "grey zone" depends on .respectively.Systems with  <  gz have more than 99% of chance to be unstable.Stability is guaranteed in the systems with  ≥  crit .
In the "grey zone" for  = 10 −5 , 10 −4 and 10 −3 , stable islands exist at  3.3, 3.1 and 3.3, respectively.We find that these islands roughly correspond to the  values of some mean-motion resonances (as marked in the figure).No stable islands are found for  = 10 −7 and 10 −6 because their  H is too small.
The "degree of instability" of a system can be characterised by  inst , the time for the instability to develop from two initially nearly circular orbits.Fig. 2 shows the  inst from our simulations.For each initial , the distribution of  inst results from the random selections of the initial argument of pericenter and mean anomaly.The horizontal spread is manually added for better display.The systems that remain stable at 10 5  1 are included in the plot and grouped in the  inst > 10 5  1 bins.
When  is small,  inst is almost always less than the synodic period of the test mass from the planet, where  1 and  2 are the initial orbital periods of  1 and  2 .This means that, except for a very small number of outliers, the mutual gravity between the planet and the test mass at their first orbital conjunction is strong enough to induce instability in a single shot.At least 80% of the instabilities are single-shot when  ≤  syn ∼ 2.5 (See Table 1).
Between  syn and  crit is a transition region, where the instability takes  inst  syn to develop.Note that  inst can vary by one to two orders of magnitude for the same  and .However, the "typical"  inst value is an exponential function of .We fit the ln ( inst )-vs- relation with ln ( inst ) = ln ( syn ) +   −  syn  syn .
(5) using the least squares method.The data points with  inst > 10 5  1 are excluded from the fitting.Table 1 gives the fitting parameters for different values of .
Table 1.Parameters for the instability timescale in the "planet + test mass" systems (of different mass ratio  =  1 / ) with no frictional force.(i)  syn , the largest  value that gives  inst <  syn for 80% of the runs; (ii)  gz , the minimum  value to have 1% systems to be stable (the left edge of the "grey zone" in Fig. 1); (iii)  crit , the minimum  value that ensures stability (the right edge of the "grey zone" in Fig. 1); (iv) , the best fit parameter  in equation ( 5) with one standard deviation error.

Systems with frictional force
Planets embedded in gaseous discs are subject to "frictional" forces that induce eccentricity damping and orbital migrations.For example, a low-mass (non-gap openning) planet (of mass  and semimajor axis ) experiences eccentricity damping on the timescale (e.g., Tanaka & Ward 2004) where Σ is the disc surface density and ℎ is the aspect ratio, and we have adopted some representative numbers for protoplanetary discs.
To assess the effect of the disc eccentricity damping to the dynamical instability, we re-run the simulations of Section 2.2 with an extra force (see Papaloizou & Larwood 2000), applied to the test mass, where  is the test mass velocity and r is the unit radial vector.For  1, this force leads to eccentricity damping  −/2.We have also experimented with a different frictional force model with where  K = √︁  /|| φ is the Keplerian velocity at the location of the test mass.This force also gives  −/2.We found that the effect of equation ( 8) on the stability of the "planet+test particle" system is similar to that of equation ( 7).In the following we describe our results based on equation (7).
Fig. 3 shows the fraction of stable runs when the frictional force (equation 7) is applied.The blue lines show the results from the weak friction case with  = 10 4  1 .The fraction of stable systems increases for all  values except around the stable islands.The friction force can protect a system from instability by reducing the orbital eccentricity of the test mass.However, the friction force can also affect the orbital energy of a test mass.Hence, inside the "greyzone" stable islands, the protective effect of the mean motion resonances is largely weakened.Overall, as the stable fraction changes from 0 to 1 with increasing , the weak friction force smooths the transition region by expanding the "grey-zone" to smaller values of  and reducing the non-monotonic features inside the "grey-zone".1, except showing the results that includes the frictional force (equation 7) on the test particle.The red lines show the cases with  = 10 2  1 and the blue lines with  = 10 4  1 .The black lines and the "grey-zone" regions are from the no-friction results (same as in Figure 1).Table 2. Critical initial  when the frictional force (equation 7) is applied in the "planet + test mass" systems for different the mass ratio  =  1 / .(i)  =14 , the largest  value that allows instability when  = 10 4  1 ; (ii)  =14 , the critical  when  = 10 4  1 as estimated by equation ( 10); (iii)  =12 , the largest  that allows instability when  = 10 2  1 ; (iv)  =12 , the critical  when  = 10 2  1 as estimated by equation ( 10).
From the fitting formula (equation 5), for each  and , there is a special initial orbital separation, obtained by setting  inst = .We define With initial  >   , a system is expected to have  inst > .That means the eccentricity damping is faster than the growth of instability and the system is expected to be stable.Table 2 gives the value of   from equation (10) and   , the true numerical value of the critical separation for 100% stability.Although equations ( 5) and (10) allow us to relate  and the "typical"  inst , we note that even systems with the same initial  can have a large spread of  inst 's (see Fig. 2).Whether a system can be stabilised by friction depends more fundamentally on its  inst than its initial .
Fig. 4 compares the  inst results from the with-friction runs (coloured dots) to those from the no-friction runs (black dots).We shade the area with  ∈ [2.0,   ] and  inst ∈ [0, ].Inside the shaded area, the instability growth is faster than the eccentricity damping ( inst < ), so the with-friction and no-frictions runs show the same distribution of  inst , and the coloured dots cover the black dots almost exactly.Outside of the shaded area, the instability is weaker ( inst > ), so the systems can be stabilised by the frictional force; thus, the black dots are not covered by the coloured dots as the coloured dots are pushed upward to the simulation time limit at  inst > 10 5  1 .
Due to the spread of  inst 's for runs with the same , a damping force with a constant  can only stabilise a fraction of the runs if  inst () ∼ .That creates a wider region of "grey zone" where the outcome of the evolution is undetermined and a more smooth transition from "all unstable" to "all stable" (see Fig. 3).The coloured outliers with  inst ∈ (, 10 5  1 ) in Fig. 4 are systems that would be protected by mean motion resonances but is now destabilised by the damping force.

Formalism of the linear map
In this section, we present an analytical map that can qualitatively reproduce the simulation results of Section 2. Our approach is modified and generalised from the previous work by Duncan et al. (1989, hereafter DQT89).
We map the evolution of the outer test particle's orbit to a succession of conjunctions with the inner planet.The leading-order change of the orbital elements of the test mass during a conjunction can be derived using the shearing box framework (Henon & Petit 1986).We define where  0 and  1 are the modified Bessel functions of the second kind, and the subscript  denotes the massive planet.We define   to be the complex eccentricity just after the  th conjunction, and  − the complex eccentricity just before the  th conjunction (see Fig. 5).We use the same notations for the relative semi-major axis  and the longitude  of the test mass at the conjunction.During the  th conjunction, the complex eccentricity is changed by the amount (see DQT89) and the following quantity (Jacobi integral) is conserved In DQT89,  is fixed to  − in equation ( 14), an approximation that becomes incorrect close to the instability boundary, where  can vary significantly.In this work, we use a leap-frog-like approach to evaluate  ,mid , the average of the pre-conjunction and the estimated post-conjunction -values: The actual post-conjunction orbital elements are In the no-damping case, the orbital elements are conserved between the conjunctions, i.e.
In the damping case, with the frictional force described as in equation (6) (corresponding to an exponential eccentricity damping with timescale ), we have at first order: where   is the synodic period after the  th conjunction.It can be derived from the orbital periods of the planet   and test particle : Orbital crossings occur when  >  and equation ( 3) is met when   −  H / p .However, both criteria lie outside the validity range of the map.In the following, we adopt the stability condition Empirical tests show that our results are insensitive to the precise numerical threshold (0.5 can be replaced by 0.4 or 0.6).
Starting from  = 0, equating |Δ| (equation 14) to  gives the Hill radius scaling ( ∝  1/3 ).Alternatively, equating the variable  part of |  −  − | to  gives the resonance overlap criterion scaling ( ∝  2/7 , Wisdom 1980, DQT89, see Appendix A).As we will see below however, these scalings do not account for the complexity of the instability threshold.

Results
Using the algebraic map described in Section 3.1, we compute the "eccentricity destabilisation time"   , which is the time for a system to reach the limiting eccentricity to break equation (27).A system is considered stable if it satisfies equation ( 27) for 1000 conjunctions (equivalent to 10 3 -10 5  1 depending on the mass ratio and the initial -value).The initial conjunction longitude is randomised, and the initial eccentricity of the test mass is set to .
Fig. 6 shows the results from the map without friction (i.e., with equations 20-22).The behaviour of   with respect to  (for a given ) resembles that of  inst from the -body simulations.Fig. 6 shows that the   - "curve" exhibits three branches.The leftmost which has the same form as equation ( 5) for  inst .Table 3 gives the values of  syn, ,  crit, and   for different 's.
To see why the stability islands exist, we analyse the map by assuming / 2 .Using equation ( 17 That means that  will not change if the conjunction happens exactly when the test mass is at its pericenter or apocenter ( −  = 0 or ).When the planet and the test particle are at conjunction, the mean longitude of the two objects are the same.Hence,  ≡  −  corresponds to the resonance angle of the planet and the test mass if they are in a first-order two-body mean motion resonance.The cumulative change of  will be zero if  "librates" (discretely at each conjunction) around 0 or .
From equations( 17) and ( 19) we find that the changes in  and  in a conjunction are given by where we wrap large angles into [0, 2).If  is a small constant because  remains close to 0 or , equation ( 18) suggests that  is which means the map allows Δ and Δ to be zero at the same time for some specific value of .The equilibrium is stable ( 2 Δ/ 2 < 0) at  = , allowing a long-term "libration" of  around  = .This is why the systems with some special initial  can have long-term stability according to the map.This stabilising mechanism can also be applied to the -body simulations and explain the stable islands found in Section 2 (see Figs. 1 and 2).Fig. 7 shows the distribution of  and Δ at conjunctions in the "planet + test mass" -body simulations.We select the cases with  = 10 −5 and  = 3.2 (red), 3.3 (blue) and 3.4 (green), where  = 3.3 is in a stable island while  = 3.2 and 3.4 are outside of the island.All simulated conjunctions approximately follow Δ ∝ − sin  as equation ( 29) predicts.Near the island (blue),  always stays between /2 and 3/2 and distributes symmetrically around  = , which agrees with the stabilising mechanism derived using the map.

Map with friction
Fig. 8 shows the eccentricity destabilisation timescale   from the map with the frictional damping included (see equation 23) for two damping timescale  = 10 4  1 and  = 10 3  1 .We have also experimented with  = 10 2  1 and found that all runs are stable.Similar to equations ( 9) and ( 10) for the -body simulations, we define the special initial spacing  *  via   () : Following equation ( 28), we use  *  to derine  , , the map estimate of the largest -value that allows instability when an eccentricity damping of timescale  is applied.In Fig. 8, we shade the area with  ∈ [2.0,  , ] and   ∈ [0, ], where the instability is expected to be stronger than the frictional damping.The result shows that the eccentricity can only grow in systems inside the shaded area.
Outside of the shaded area, where   > , the damping is too strong for the eccentricity to grow.Fig. 9 compares the size of the transition zone [ syn, ,  crit, ] from the map to  crit from the -body simulation (see Tables 1 and  2) in the no-friction and friction case.Better agreement between the stability thresholds is found when the damping is applied.Although the map does not exactly reproduce the results of the -body simulations, it correctly predicts the orders of magnitude of the instability timescale and the effect of friction on stability.The grey zone in Fig. 9 is impacted by the damping term in the same way as the effect of frictions in the -body simulations (see Section 2.3).The advantage of the map is twofold: it is much quicker to compute than -body simulations, which allows a better sampling of the parameter space, and its constitutive equations shed some light on the instability mechanisms.

INSTABILITY OF MULTI-PLANET SYSTEMS WITH FRICTIONS
In this section, we consider systems with multiple (non-zero mass) planets.Our simulated systems consist of a central star with mass Transition zone between stability (all particles are stable) and instability (particles are unstable in less than two synodic periods) as a function of the mass ratio .For each , 100 test particles are evolved for 1000 conjunctions using the map.The green line is a fit to the end of the transition zone without friction ( ∝  2/7 , as obtained by DQT89).The blue line is the predicted effect of the friction using the  -body simulations (Table 1 and equation 10).The discrepancy is large between the  -body and map results for the no-friction case (left), mostly due to the different criteria for instability.On the other hand, the results are much more similar when we add friction, as both versions of instability behave alike.
and planets with mass   = 2 2−   for  = 1, 2, .., .The initial orbital separation between two neighbouring planets is set as where  the same dimensionless orbital separation for all neighbouring pairs and is the mutual Hill radius of   and  +1 .The innermost (also the most massive) planet is given initial eccentricity  1 = 0 and other planets (  ≥ 2) have initial eccentricities   = 10 −5 .All planets are co-planar and have zero mutual inclination.For some experiments, we apply equation ( 7) to all planets to model the frictional damping.Same as in the previous sections, we consider different combinations of ,  and .We carry out 50 runs for each combination with randomised initial angles, anomalies and longitudes.Each simulation runs up to 10 6  1 .
It has been well recognised that stability in two-planet systems and that in the systems with more than two planets are substantially different: the Hill stability criterion only applies to two-planet systems; systems with three or more planets can become unstable even when the initial planetary separations are large, although the instability growth time increase rapidly with the separation (see Pu & Wu 2015, and refs therein).Hence, we will discuss the results from the simulations with  = 2 and  = 3 separately.
Fig. 10 shows the fraction of stable runs in the two-planet simulations.We give the values of  gz and  crit for the "grey zone" (where a system can be either stable or unstable) in the "no-friction"  runs in Table 4. Unlike the test-mass case (see Fig. 1), the "grey zones" for two-planet systems do not contain any stability island, except at  = 2.7 to 3.2 for  = 10 −3 .The strong instability at  ∼ 3.3 for  = 10 −3 almost coincides with the location of the 5 : 3 mean-motion resonance in the system.When the friction is applied, the stable fraction increases almost everywhere, and decreases inside the stability islands for for  = 10 −3 .Fig. 11 shows the instability growth time,  inst , from the nofriction simulations (with the systems that remain stable in 10 6  1 being added as the  inst > 10 6  1 dots).Similar to the test-mass cases (see Fig. 2), when  is small (≤  syn given in Table 4),  inst is almost always less than the synodic period (equation 4) of the two planets: the mutual gravity between the planets is strong enough to trigger immediate instability at their first orbital conjunction.
Between  syn and  crit , the instability growth time  inst can vary by several orders of magnitude for the same  and .The range of the instability time is roughly an exponential function of , despite the two outliers in the  = 10 −4 case at  = 3.5 and 3.9.We fit the data points in the transition region that has  inst < 10 6 with equation (5).Table 4 lists the fitting results.
Fig. 12 compares  inst 's from the with-friction runs (coloured dots) to the results from the no-friction runs (black dots).Similar to Fig. 3, instability can only happen inside the shaded area with  ∈ [2.0,   ] and  inst ∈ [0, ] (See Table 5 for the values of   , as well as the numerical   ).Outside of the shaded area,  inst >  and instability is suppressed.However, because  inst spreads a couple orders of magnitude for runs with the same  and , there is no 10 0 10 1 10 2 10 3 10 4 10 5 > 10 6 = 1e 6 10 0 10 1 10 2 10 3 10 4 10 5 > 10 6 = 1e 5 10 0 10 1 10 2 10 3 10 4 10 5 > 10 6 Instability timescale t  Table 4. Parameters for the instability timescale in the two-planet systems (of different mass ratio  =  2 / ) with no frictional force.(i)  syn , the largest  value that gives  inst <  syn for 80% of the runs; (ii)  gz , the minimum  value to have 1% systems to be stable (the left edge of the "grey zone" in Fig. 1); (iii)  crit , the minimum  value that ensures stability (the right edge of the "grey zone" in Fig. 1); (iv) , the best fit parameter  in equation ( 5) with 1 error.
10 −6 3.4 3.4 2.9 2.8 10 −5 3.3 3.3 3.1 2.9 10 −4 3.2 3.7 3.1 2.9 10 −3 3.4 3.4 2.7 2.7 with more than two planets.Fig. 13 summarises the fraction of systems that remain stable after 10 6  1 in our simulations.Fig. 14 shows the instability growth time  inst in the no-friction runs: most of the runs become unstable within 10 6  1 .The exceptions are at  ≥ 4.6 for  = 10 −3 , where  inst is too long due to the large initial planet spacings, and similarly at  = 3.8 for  = 10 −3 and  = 4.7 and 4.9 for  = 10 −4 , where  inst are larger than at their neighbouring 's possibly due to mean motion resonances.We expect the stable fraction (Fig. 13) to converge to zero for all   considered  values if we evolve the systems for more orbits.We find the value of  syn for each  using the synodic period of the two inner planets and fit the data points with  >  syn and  inst < 10 6  1 using equation ( 5).Table 6 lists the results.
When the frictional force with  = 10 4  1 is applied, the survival rate increases at large .Notably, around half of the systems with  = 10 −5 are stable at  = 5.0.Systems with  = 10 −4 have a 75% chance to be stable at  = 5.0 and more than 50% chance to survive at  ≥ 4.8.In the  = 10 −3 case, 100% of our simulations with  ≥ 4.6 are stable.The stable fraction increases dramatically from ∼ 0 to 90% at  = 3.5.
Fig. 15 displays  inst when the frictional force is applied.Unlike in Figs. 4 and 12, we see a significant amount of coloured dots (from the with-friction runs) outside of the shaded area (i.e. [2.0,   ] × [0, ]).The stability in three-planet systems cannot be simply determined by the condition  inst >  like in the two-planet case.Nevertheless,  inst still provides an estimate on what value  is needed to stabilise a three-planet system.Our results suggest that, a planetary system may not be stabilised by the frictional force even if  inst 10 − 100.Therefore, the critical initial  for stabilisation is larger than   that one would estimate using equation 10 (see Table 7).

Summary of key results
In this paper, we have studied the dynamical instability of closelypacked multi-planet systems in the presence of frictional forces that arise from planet-disc interactions.The systems have initially coplanar, near circular orbits, with the orbital separation characterised by the dimensionless ratio  = ( +1 −   )/ H (where   is the semi-major axis of the -th planet, and  H is the mutual Hill radius).Instability occurs due to planetary orbital crossing (see equation 3) as the system evolves.The goal of this paper is to evaluate how frictional forces (of various strengths) affect the growth of instability as a function of  and the planet-to-star mass ratio .Note that although we use the terms "planet" and "star" throughout this paper, our results are also relevant for understanding the evolution of stellar-mass black holes embedded in AGN discs around supermassive black holes.We consider both "2 planets" and "3 planets" systems using numerical -body integrations, and carry out theoretical analysis of the restricted three-body problem (star, planet and test particle) to gain analytical understanding.For each planetary architecture, we adopt a range of planet-to-star mass ratios (), initial dimensionless orbital separations () and frictional damping timescales (), and determine the fraction of stable systems and instability growth time  inst (i.e., the time to reach the first orbital crossing from initially circular orbits).In general, we find that the stable (large  inst and ) and unstable (small  inst and ) regimes are separated by a grey/transition zone, which can have complicated stable fraction-vs- dependence because of mean-motion resonances; this transition zone becomes "smoother" for more extreme mass ratios ( 10 −6 ).Frictional forces tend stabilise a system, pushing the unstable regime towards smaller 's.
For systems with no frictional forces, our key results can be summarised as follows: • In the "planet + test-mass" and "2 planets" simulations, the fraction of stable systems is 100% for  >  crit 2 √ 3 and is less than 1% for  <  gz , where  gz ∈ [2.6, 3.1] depending on the mass ratio  (see Tables 1 and 4).The region  ∈ [ gz ,  crit ] is a grey zone where there is no guarantee whether dynamical instability will occur or not (see Figs. 1 and 10).Inside the grey zone, the stable fraction does not always increase monotonically with  because of mean-motion resonances.In the test-mass cases with  ≥ 10 −5 , we find stable islands at some particular 's where the stable fraction is enhanced (see Fig. 1).
• The instability timescales  inst for different planetary architectures are shown in Figs. 2, 11 and 14.For systems with the same  and initial ,  inst can vary by 2 to 3 orders of magnitude due its dependence on the initial orbital phases of the planets.In the transition zone, the averaged  inst generally has a non-monotonic dependence on .
We then re-run the simulations with a frictional force that damps the planetary eccentricity (see equation 7).There are two main effects associated with the frictional force: • The frictional force generally increases the stable fraction in all simulations.To be stabilised, a two-planet system needs to have  inst  (see Figs. 4 and 12) and a three-planet system needs to have  inst 10 − 100 (see Fig. 15).
• In a two-planet system, the frictional force may lower the planet survival fraction in the stable islands.A weak friction force can smooth the stable fraction vs  curve inside the transition zone.When there are three planets, the frictional force can make the systems with relatively large  to achieve long-term stability.Due to the spread of  inst at the same , a three-planet system may acquire "grey zone" when the frictional force is applied (see Fig. 13).
We have also devised a linear map to analyse the dynamical instability of the "planet + test-mass" system (see Section 3).This map serves as a useful model that captures the key features of the stability, especially the effect of the damping force.

Discussion
Our results are useful for understanding the evolution of multi-planet systems born in protoplanetary discs, and the evolution of multiple stellar-mass black holes (BHs) in ANG discs (see Section 1).Multiple planets can be brought into closely-packed orbits due to differential migration or migration traps in the disc.The frictional forces from the disc acting on the planets can help maintain the stability of the system even for small planetary spacings -"how small" depends on the frictional damping timescale, which in term depends on the disc density.
In this paper, we have focused on the onset of dynamical instability, signalled by close encounters between two planets (see equation.3).Once the instability is initiated, the system will experience chaotic orbital evolution.In the absence of gas friction, close encounters keep recurring until one planet is ejected, or two planets collide with each other, or one planet crashes into the star.The branching ratio of each outcome depends on the planetary mass, radius and orbital velocity around the host star (see Li et al. 2021a, and references therein).In the case of stellar-mass BHs around a supermassive BH, the large mass ratio implies that ejection is highly unlikely or impossible; and the chaotic orbital evolution is terminated when two BHs undergo an extremely close encounter and form a bound, merging binary due to gravitational radiation (see Li et al. 2022a, and references therein).
When dynamical instability develops in the presence of gas discs, there exists another possible outcome in which the system regains its stability after the chaotic evolution: the combined effects of close encounters and frictional forces may push the planets (or BHs) into well separated orbits to ensure dynamical stability with friction (Li et al. 2022a).In particular, we expect that some planetary ejections and star crashings may be prevented by the frictional forces because they typically require a large number of closer encounters than planet-planet collisions.

Figure 1 .
Figure 1.Fraction of stable runs for systems with different  and initial  in the "planet + test mass" simulations with no frictional force.The "greyzone" regions start at the minimum  with stable fraction > 1% and end where the stable fraction reaches 100%.The dashed red lines mark the  values that correspond to (  + ) :  mean-motion resonances for  ≤ 20 and  ≤ 3.

Figure 2 .
Figure 2. Instability timescale  inst (black dots) for the simulations shown in Fig. 1.The data points for each  are manually spread on the horizontal axis for better display.The vertical dashed lines indicate  syn and  crit , and the "grey-zone" regions are shaded.At  ≤  syn , the magenta triangles mark the synodic period between the test mass and the planet.The magenta lines show the fitting results using equation (5) in the transition regions between  syn and  crit .

Figure 3 .
Figure3.Same as Fig.1, except showing the results that includes the frictional force (equation 7) on the test particle.The red lines show the cases with  = 10 2  1 and the blue lines with  = 10 4  1 .The black lines and the "grey-zone" regions are from the no-friction results (same as in Figure1).

Figure 4 .
Figure 4. Instability timescale  inst for "planet + test mass" systems with frictional force applied.The red dots show the results for the  = 10 2  1 case and the blue dots show the results for the  = 10 4  1 case.The black dots are the no-friction results (same as in Figure 2).The area within [2.0,   ] × [0,  ] are grey shaded.
Figure 5. Setup of the linear map and the definition of the key parameters and variables.

Figure 6 .
Figure 6.Similar to Fig. 2, except showing the eccentricity destabilisation timescale   (black dots) from the map (see Section 3.2.1).No damping force is applied.The data points are manually spread on the horizontal axis for better display.The vertical dashed lines indicate  ≤  syn and  crit, .The magenta triangles mark the synodic period between the test mass and the planet for  syn .The magenta lines show the fitting results of   between  syn and  crit, using equation (28).

Figure 7 .
Figure 7. Distribution of  =  −  and Δ at conjunctions in the  -body simulations with  = 10 −5 and  = 3.2 (red), 3.3 (blue) and 3.4 (green).This suite of  -body simulations finds a stability island at  = 3.3.Only the conjunctions during the first 1000 orbits are shown.

Figure 8 .
Figure 8. Similar to Fig. 4, but showing the eccentricity destabilisation timescale   (which is the time when equation 27 breaks down) from the map.

Figure 9 .
Figure9.Transition zone between stability (all particles are stable) and instability (particles are unstable in less than two synodic periods) as a function of the mass ratio .For each , 100 test particles are evolved for 1000 conjunctions using the map.The green line is a fit to the end of the transition zone without friction ( ∝  2/7 , as obtained by DQT89).The blue line is the predicted effect of the friction using the  -body simulations (Table1 and equation 10).The discrepancy is large between the  -body and map results for the no-friction case (left), mostly due to the different criteria for instability.On the other hand, the results are much more similar when we add friction, as both versions of instability behave alike.

Figure 10 .
Figure 10.Similar to Fig. 3, except showing the fraction of stable runs for systems with two finite-mass planets after 10 6  1 .The black lines are from the "no-friction" runs.The red lines show the cases with  = 10 2  1 and the blue lines show the cases with  = 10 4  1 .

Figure 11 .
Figure 11.Similar to Fig. 2, except showing the instability timescale  inst (magenta dots) for systems with two finite-mass planets.No friction force is applied.The vertical dashed lines indicate  syn and  crit , and the "greyzone" regions are shaded.At  ≤  syn , the magenta triangles mark the synodic period between the the two planets.The magenta lines show the fitting formulae using equation (5) in the transition region (which span from  syn to  crit ).

Figure 12 .
Figure 12.Similar to Fig. 4, but showing the instability time  inst for twoplanet systems with frictional force applied.The red dots show the results from the  = 10 2  1 cases and the blue dots show those from the  = 10 4  1 cases.The black dots are from the no-friction runs (same as in Figure 11).The area within [2.0,   ] × [0,  ] are grey shaded.
Table 1 lists the values for the left and right boundary of the "grey zone", which we denote as  gz and  crit ,

Table 3 .
Parameters for the transitional branch of the eccentricity destabilisation time from the no-friction map calculations: (i)  syn,e , same as the  syn in Table2; (ii)  crit, , the largest  -value that allows instability; (iii)   , the best fit for the slope in equation (28) with one standard deviation errors.is very short and slightly decreases as  increases.The rightmost region (except for the  = 10 −7 case) covers  >  crit, , where  crit,e is the largest  that allows instability.The middle branch is the transition region, where the "averaged"   increases exponentially with .Note that the actual   has a wide spread for given  and ; e.g.,   can differ by one to two orders of magnitude for two systems with the same  and .Quantitatively, even if they are closely related, the  syn, and  crit, do not always equal to  syn and  crit (see Section 2) since the map and the -body simulation use different instability criteria.For the map,   is always greater than or equal to the synodic period  syn (see equation 4).We fit   in the transition branch with ln (  ) = ln ( syn ) +    −  syn,  syn,

Table 6 .
Parameters for the instability timescale in the three-planet systems (of different mass ratio  =  2 / ) with no frictional force.(i)syn, the largest  value that gives  inst <  syn for 80% of the runs; (ii) , the best fit parameter  in equation (5) with 1 error.