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

TL;DR

This study investigates how atmospheric mass-loss and giant impacts influence the stability and evolution of closely packed multiple planetary systems, highlighting the importance of initial spacing and mass-loss limits for system stability.

Contribution

It provides new constraints on the extent of mass-loss in multiple planetary systems, especially in compact and resonant configurations, through idealized dynamical simulations.

Findings

Modest mass-loss can destabilize closely packed systems.

Wide initial separation or limited mass-loss is necessary for stability.

Giant impacts and >a few percent mass-loss induce instabilities.

Abstract

A gap in exoplanets' radius distribution has been widely attributed to the photo-evaporation threshold of their progenitors' gaseous envelope. Giant impacts can also lead to substantial mass-loss. The outflowing 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 which are closely spaced and/or contain planets locked in mean motion resonances (MMRs), 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.

Tables

Table 1. 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 𝑘 k is the initial planetary separation in units of mutual Hill radii, τ m subscript 𝜏 𝑚 \tau_{m} is the mass-loss timescale. The cases shown in red color are stable for at least 10 7 superscript 10 7 10^{7} yrs, while the cases in black represent the systems with orbital crossing happened before the end of simulations.

Case	A	Mass	$k$	${\tau }_m$	Mass	$k$	${\tau }_m$	Mass	$k$	${\tau }_m$
		$M_{\oplus }$	$R_{hill}$	yr	$M_{\oplus }$	$R_{hill}$	yr	$M_{\oplus }$	$R_{hill}$	yr
1	0.05	1	10	${10}^4$	3	10	${10}^4$	10	10	${10}^4$
2	0.05	1	10	${10}^5$	3	10	${10}^5$	10	10	${10}^5$
3	0.05	1	10	${10}^6$	3	10	${10}^6$	10	10	${10}^6$
4	0.05	1	12	${10}^4$	3	12	${10}^4$	10	12	${10}^4$
5	0.05	1	12	${10}^5$	3	12	${10}^5$	10	12	${10}^5$
6	0.05	1	12	${10}^6$	3	12	${10}^6$	10	12	${10}^6$
7	0.05	1	14	${10}^4$	3	14	${10}^4$	10	14	${10}^4$
8	0.05	1	14	${10}^5$	3	14	${10}^5$	10	14	${10}^5$
9	0.05	1	14	${10}^6$	3	14	${10}^6$	10	14	${10}^6$
10	0.05	1	16	${10}^4$	3	16	${10}^4$	10	16	${10}^4$
11	0.05	1	16	${10}^5$	3	16	${10}^5$	10	16	${10}^5$
12	0.05	1	16	${10}^6$	3	16	${10}^6$	10	16	${10}^6$

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

Case	Mass	A	$k$	${\tau }_m$	A	$k$	${\tau }_m$
	$M_{\oplus }$		$R_{hill}$	yr		$R_{hill}$	yr
1	10	0.1	10	${10}^4$	0.2	10	${10}^4$
2	10	0.1	10	${10}^5$	0.2	10	${10}^5$
3	10	0.1	10	${10}^6$	0.2	10	${10}^6$
4	10	0.1	12	${10}^4$	0.2	12	${10}^4$
5	10	0.1	12	${10}^5$	0.2	12	${10}^5$
6	10	0.1	12	${10}^6$	0.2	12	${10}^6$
7	10	0.1	14	${10}^4$	0.2	14	${10}^4$
8	10	0.1	14	${10}^5$	0.2	14	${10}^5$
9	10	0.1	14	${10}^6$	0.2	14	${10}^6$
10	10	0.1	16	${10}^4$	0.2	16	${10}^4$
11	10	0.1	16	${10}^5$	0.2	16	${10}^5$
12	10	0.1	16	${10}^6$	0.2	16	${10}^6$

Table 3. Table 3: Parameters used in the simulations, Group 1, Group 2, and Group 3. k 0 = ( a i + 1 − a i ) / R H subscript 𝑘 0 subscript 𝑎 𝑖 1 subscript 𝑎 𝑖 subscript 𝑅 𝐻 k_{0}=(a_{i+1}-a_{i})/R_{H} represents the relative separation between two adjacent planets initially. r t subscript 𝑟 t r_{\rm t} is the transition radius, τ dep subscript 𝜏 dep \tau_{\rm dep} is the depletion timescale of the gas disk, c ​ o ​ e ​ f 𝑐 𝑜 𝑒 𝑓 coef is the coefficient in the corotation torque, h=H/r represents the disk’s aspect ratio, and f g subscript 𝑓 𝑔 f_{g} is the enhancement factor of the gas density.

	$k_0$	${\tau }_{\mathrm{dep}}$	$r_t$	$coef$	$h=H/r$	$f_g$
	$R_{\mathrm{hill}}$	Myr	AU
G1	8	1	0.5	10	$0.05r^{0.25}$	0.1
G2	16	1	3.0	10	0.05	0.1
G3-1	16	1/3	3.0	10	0.05	1
G3-2	16	1/3	3.0	10	0.05	0.5
G3-3	16	1/3	3.0	10	0.05	0.2
G3-4	16	0.1	3.0	10	0.05	0.2
G3-5	16	1	3.0	10	0.05	0.2
G3-6	16	3	3.0	10	0.05	0.2

Table 4. 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 𝑎 a , e 𝑒 e , i 𝑖 i , Ω Ω \Omega , ω 𝜔 \omega , f 𝑓 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.

Name	$a$	$e$	$i$	$\Omega$	$\omega$	$f$	Mass
	AU		rad	rad	rad	rad	$M_{\odot }$
1	0.753710	0.058854	0	0	0.387958	-1.367171	1e-5
2	0.987895	0.102676	0	0	-2.749728	-0.024010	1e-5
3	1.293961	0.077092	0	0	0.362406	-0.952595	1e-5
4	1.696034	0.064311	0	0	-2.754476	-1.830721	1e-5
5	2.222124	0.044478	0	0	0.591902	-2.323546	1e-5
6	2.912175	0.039683	0	0	-2.227325	-2.776581	1e-5
7	3.820460	0.024706	0	0	-0.013042	-0.058686	1e-5
8	4.623438	0.024523	0	0	2.482278	0.533084	1e-5
9	5.608244	0.011006	0	0	-0.760252	2.023750	1e-5

Table 5. Table 5: The parameters of mass-loss in the cases in section 3.2. A is the mass-loss fraction, τ m subscript 𝜏 𝑚 \tau_{m} is the mass-loss timescale, θ ¯ / 2 ​ π ¯ 𝜃 2 𝜋 \bar{\theta}/2\pi is the average libration amplitude of all the planet pairs, and e ¯ ¯ 𝑒 \bar{e} is the average eccentricity of all planets in the system.

Case	A	${\tau }_m$ (yr)	$\overline{\theta }/2\pi$	$\overline{e}$
1	0.05	${10}^6$	0.029	0.147
2	0.05	${10}^5$	0.029	0.150
3	0.05	${10}^4$	0.27	0.179
4	0.1	${10}^6$	0.17	0.195
5	0.1	${10}^5$	0.18	0.230
6	0.1	${10}^4$	0.36	-

Equations

$${\tau }_m=-\frac{m_p}{{\dot{m}}_p},$$

$$L=m_p\sqrt{G{M}_{\ast }{a}_p(1-e^2)},$$

$$\frac{{\dot{a}}_p}{a_p}=2\left(\frac{\dot{L}}{L}-\frac{{\dot{m}}_p}{m_p}+\frac{e\dot{e}}{(1-e^2)}\right).$$

$$m_p(t)=m_0[1+\mathrm{A}exp(-t/{\tau }_m)]/(1+\mathrm{A}).$$

$$\frac{\dot{a}}{a}=\frac{2f_p{\dot{m}}_p}{m_p}=\frac{-2f_p}{{\tau }_m}\left(\frac{\mathrm{Aexp}(-t/{\tau }_m)}{1+\mathrm{Aexp}(-t/{\tau }_m)}\right)$$

$$\frac{d}{dt}{\text{V}}_i=-\frac{G(M_{\ast }+m_i)}{{r_i}^2}\left(\frac{{\text{r}}_i}{r_i}\right)$$

$$+\sum _{j\ne i}^{N}G{m}_j\left[\frac{({\text{r}}_j-{\text{r}}_i)}{|{\text{r}}_j-{\text{r}}_i{|}^3}-\frac{{\text{r}}_j}{r_j^3}\right]$$

$$+{\text{F}}_{\theta }+{\text{F}}_{\mathrm{damp}}+{\text{F}}_{\mathrm{migI}},$$

$$F_{\theta }=\frac{\dot{J}}{r_i}=\frac{J}{r_i}\left(\frac{{\dot{a}}_p}{2a_p}\right).$$

$$a_{i+1}=a_{i}+kR_{H},\ \ \ \ \ \ \ \ \$$

$${\Gamma }_{\mathrm{tot}}(r_p)=f_{tot}{\Gamma }_0,\mathrm{where}$$

$${\Gamma }_0=(\mu /h{)}^2{\Sigma }_p{r}_p^4{\Omega }_p^2/\gamma ,$$

$${\Sigma }_{\mathrm{g}}(r)=\frac{1700f_{\mathrm{g}}}{(r/1AU)}\exp \left(\frac{-\mathrm{t}}{{\tau }_{\mathrm{dep}}}\right)\mathrm{g}\mathrm{c}{\mathrm{m}}^{-2},$$

$$f_{CR}=\mathrm{coef}\times \left(1-\frac{2(r/r_t{)}^2}{1+(r/r_t{)}^2}\right),$$

$$f_{LB}=1-\frac{2(r/r_m{)}^4}{1+(r/r_m{)}^4}.$$

$${\tau }_a=\frac{a}{\dot{a}}=\frac{M_p\sqrt{G{M}_{\ast }r}}{2{\Gamma }_{\mathrm{tot}}}=\frac{{\tau }_0(r_p)}{f_{\mathrm{tot}}},$$

$${\text{F}}_{\mathrm{migI}}=\frac{{\Gamma }_{\mathrm{tot}}\hat{\theta }}{M_{p}r}\mathrm{and}$$

$${\text{F}}_{\mathrm{damp}}=-\frac{(\mathbf{v}\cdot \mathbf{r})\mathbf{r}}{r{\tau }_e}.$$

Full text

Exoplanet dynamics (490) — Exoplanet evolution (491) — Exoplanet formation (492)

1 Introduction

More than 5000 explanets around mature stars have been identified and confirmed (Howard et al., 2012; Batalha et al., 2013; Fressin et al., 2013; Petigura et al., 2013; Fabrycky et al., 2014; Silburt et al., 2015; Zink et al., 2019; Kunimoto & Matthews, 2020). In order to utilize their dynamical properties to extract clues and place constraints on scenarios of planet formation, we must take into consideration planets’ post orbital formation and structural evolution (Ida & Lin, 2010; Ida et al., 2013; Jin et al., 2014; Helled et al., 2014). In this quest, over 800 members of multiple planetary systems provide valuable clues and best tests for theoretical models (Lissauer et al., 2011; Lovis et al., 2011; Fabrycky et al., 2014; Weiss et al., 2018; Millholland & Winn, 2021).

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 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).

This theoretical prediction has been confirmed by an observed gap (based on Kepler’s data) in the planets’ size distribution between $\sim 1.3-2R_{\oplus}$ over a range of periods (Fulton et al., 2017; Petigura et al., 2017; Berger et al., 2018; Fulton & Petigura, 2018; Berger et al., 2020; Van Eylen et al., 2018; Petigura et al., 2022). The masses of planets on either side of the gap have been determined from the transit timing variation (TTV) (Lithwick et al., 2012; Hadden & Lithwick, 2014, 2017; Wu & Lithwick, 2013; Jontof-Hutter et al., 2016) and radial velocity measurements (Weiss & Marcy, 2014; Hadden & Lithwick, 2014; Dressing et al., 2015; Wolfgang et al., 2016; Chen & Kipping, 2017). Their nearly uniform masses match well with the theoretically predicted threshold transition across an “evaporation valley” (Lopez & Fortney, 2013; Owen & Wu, 2013, 2017; Weiss et al., 2013; Wu & Lithwick, 2013; Jin et al., 2014; Owen & Adams, 2019; Mordasini, 2020).

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}\sim 10M_{\oplus}$) and a modest fraction ($f\sim$ 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_{\oplus}$) 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 millimeter- to meter- sized pebbles; (Ormel & Klahr, 2010)), until they attain the pebble-isolation mass, typically a few to 10$M_{\oplus}$ 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 $\sim 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 endure intense stellar irradiation. In principle, photo-evaporation 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).

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.

2 Stability of Systems with Mass-Loss of the Innermost Planets

In this section, we illustrate the consequence of mass-loss 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.

2.1 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 $\tau_{m}$ such that

[TABLE]

where $\dot{m}_{p}$ is the mass-loss rate of planet. The planet’s orbital angular momentum $L$ is

[TABLE]

such that

[TABLE]

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 ${\dot{e}}$. 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. $\dot{L}=0$) and the planets’ orbits expand on a time scale $\tau_{a}=a_{p}/\dot{a}_{p}=\tau_{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 $\tau_{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

[TABLE]

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 $\tau_{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

[TABLE]

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

[TABLE]

where $\textbf{F}_{\theta}={{F}_{\theta}}{\hat{\bf\theta}}$ is being applied at the radius $r_{i}$ in the azimuthal (${\bf\theta}$) direction with

[TABLE]

In the above equation, $J=\sqrt{GM_{*}a_{p}}$ is the planets’ specific angular momentum at the guiding center of their epicycles. Note that $F_{\theta}$ vanishes for $t>>\tau_{m}$. Two additional forces $\textbf{F}_{\rm damp}+\textbf{F}_{\rm 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_{\theta}$.

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_{\ast}=1M_{\odot}$ 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.

2.2 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

[TABLE]

where $R_{H}=(2\mu/3)^{1/3}(a_{i}+a_{i+1})/2$ is the Hill radius, $\mu=M_{p}/M_{\ast}$, 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 $\sim-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_{\oplus}$, respectively. For each group, 12 models are computed with different values of $k$ (initial separation) and mass-loss time scale ($\tau_{m}$).

We list the initial parameters in Table LABEL:allcase. 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 LABEL:allcase are also shown in the upper panels of Figure 3. The low-mass ($1M_{\oplus}$) group have relatively small $R_{H}$ such that some compact systems (with modest $k$) may have overlapping MMRs. Only systems with large initial $k(\geq 16$) and long mass-loss time scale $\tau_{m}\geq 10^{6}$ remain stable for more than $10^{7}$ yrs.

In comparison with the $M_{p}=1M_{\oplus}$ planets, the magnitude of $R_{H}$ is larger for the $M_{p}=3M_{\oplus}$ planets. Adjacent planets with $k\geq 16$ are stable because they avoid most MMR’s, regardless the mass-loss time scale $\tau_{m}$. In Figure 2, we show the evolution of a typical stable (for $>10^{7}$ yrs) system (with $M_{p}=3M_{\oplus}$, $A=0.05$ and $\tau_{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_{\oplus}$), $R_{H}$ is sufficiently large, all systems with $k\geq 10$ are stable for more than $10^{7}$ yrs. In general, systems with larger $M_{p}$, longer $\tau_{m}$, larger initial $k$ values are more stable.

These results are generalized to cases with non-conservative 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 mass-losing 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, non-conservative 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 LABEL:allcase.

The gaseous envelope contributes $\sim 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_{\oplus}$, $f_{p}=-1$ (outward migration), a range of values for $k$ and $\tau_{m}$ (Table LABEL:high). 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.

3 Planetary systems with MMRs

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.

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.

3.1 Multiple planetary systems with MMRs

Resonant systems are common, $\sim 20\%$ and $\sim 10\%$ of Kepler’s multiple-planet systems are near their 2:1 and 3:2 MMRs, respectively (Wang et al., 2012; Lee et al., 2013; Liu et al., 2015; Wang & Ji, 2014, 2017; Wang et al., 2021). Chain resonances with period ratio in 4:2:1 MMRs or 6:3:2 MMRs have been found in some triple systems. Two sets of 8:6:4:3 and 6:4:2:1 resonant chains are found among four super-Earths around Kepler-223 (Delisle, 2017; Mills et al., 2016) and Kepler-79 (Wang et al., 2012), respectively. They have masses $\sim 3.5-7~{}M_{\oplus}$. Five planets, separated by their mutual 3:2 MMR’s, are found with radii ($1.6-3.2~{}R_{\oplus}$) in K2-138 system on either side of size gap (Christiansen et al., 2018). Six out of seven planets in the TRAPPIST-1 system are in near chain resonances (Gillon et al., 2017). 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. In-situ 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.

3.2 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, 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

[TABLE]

where $\gamma=1.4$ is the adiabatic index, $h$ is the aspect ratio, and $\Omega_{p}$ is the Keplerian frequency at the location of the planet. We adopt the empirical minimum mass solar nebular model (MMSN; Hayashi 1981) for $\Sigma_{p}=\Sigma_{g}(r_{p})$ with

[TABLE]

where $f_{\rm g}$ is the enhancement factor. Herein, we consider the gas depletion time scale $\tau_{\rm dep}$ in the range of [0.1, 3] Myrs (Haisch et al., 2001).

The total-torque coefficient $f_{\rm tot}=f_{CR}+f_{LB}$ is the sum of corotation ($f_{CR}$) and Lindblad ($f_{LB}$) torque on the planets (Paardekooper et al., 2010, 2011). The magnitude and sign of $f_{CR}$ and $f_{LB}$ are functions of the evolving distribution of $s=\partial{\rm ln}\Sigma_{g}/\partial{\rm ln}r$, $\beta=\partial{\rm ln}T/\partial{\rm ln}r$, viscosity $\nu$ and thermal diffusivity $\xi$. 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

[TABLE]

where we consider a range for the torque coefficient ${\rm coef}$ (2, 5, or 10). While $r_{t}\sim$ 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

[TABLE]

Tidal torque due to disk-planet interaction leads to type I migration on a time scale

[TABLE]

where $\tau_{0}=M_{\ast}(h^{2}\gamma/2\mu\Omega_{p})/(\Sigma_{p}r_{p}^{2}$) (Tanaka et al., 2002; Kley & Nelson, 2012). The associated time scale for eccentricity damping is $\tau_{e}=e/{\dot{e}}=h^{2}\tau_{a}$. These contributions are incorporated into the Equation of motion (Eq 6) through

[TABLE]

[TABLE]

Aside from $\textbf{F}_{\rm migI}$ and $\textbf{F}_{\rm 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).

3.3 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_{\odot}$) planets designated to be P1 to P9 in order to increasing $a$ around a $M_{\ast}=1M_{\odot}$ (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 semi-major 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 $\tau_{\rm dep}(=0.1,1/3,1,3$ Myr). For the planets, we adopt the same value for ${\rm coef}(=10)$ and $r_{m}(=0.05$ AU) for two initial sets of $k(=8$, 16) and $r_{t}(=0.5,3$ AU) (Table LABEL:run1). 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/1\rm AU)^{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.2 $R_{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\leq 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\times 10^{5}$ yrs with P8 firstly, and into a 4:3 MMR at $\sim 5\times 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 the exception of that between P7 and P8 ($\sim 83^{\circ}$). Their eccentricities are $\sim 0.1$ (see Table LABEL:run2). 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 $\sim 5$ Myrs. In a less-massive disk (Model G3-2 with $f=0.5$ in panel b), the planets are captured into similar MMR’s in $\sim 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 $\tau_{\rm 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.

3.4 Island of stability near the MMRs

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.

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 constant 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}_{\theta}$, $\textbf{F}_{\rm damp}$, and $\textbf{F}_{\rm 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\leq f_{k}\leq 1.01$) or a large ($f_{k}\geq 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 $\tau_{\rm 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 LABEL:run1). 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 $\tau_{\rm 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}\leq 0.9$) of its original stable configuration can lead to quick orbital crossing, the crossing time shorter than $10^{4}$ yrs.

3.5 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 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 LABEL:6case. The parameters are adopted for the fractional amount ($A$) and time scale ($\tau_{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 LABEL:run2).

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 $\pi$. 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 change due to the mass-loss from the innermost planet. Their MMR configuration is maintained. Although the magnitude of $\tau_{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 $\tau_{\rm dep}=10^{4}$ yrs) the planets’ eccentricities $e$ and libration amplitude $\theta$ of the resonant angles between adjacent MMR pairs grow substantially (Figure 10). The planets’ average eccentricities $\bar{e}$ and average libration amplitude $\bar{\theta}$ in Case 3 are much larger than those in Case 1 (Table LABEL:6case). 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 $\tau_{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 LABEL:6case, are similar to that in Case 3. In Case 6, $\tau_{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 ($\bar{\theta}$) grow to $\pm\pi$ 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 ($\sim 5\%-10\%$) fractional mass-loss, compact systems with multiple MMRs have a tendency 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 LABEL:run2. 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 $\tau_{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 configuration (with period ratio of $\sim$ 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 $\sim$ 1.46 and the period ratio is $\sim$ 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.

4 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 ($\sim 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 after 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 $\tau_{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 ($\geq 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 mass-loss process and the libration of the resonant angles. When the libration timescale is larger than the mass-loss timescale, the amplitude of the resonance angles become large enough to break the resonance making orbital crossing happened in the system. With the timescale of mass-loss less than $2\times 10^{4}$ years, planet pairs may escape from MMR. The system can be destroyed if the innermost planet lost $\sim 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 $\sim$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), and TOI-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), and TOI-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 mass-loss 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 six-planet 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 mass-loss process as concluded in the section §3. If the planets in multiple planetary system sustain mass-loss process for shorter than $\sim~{}2\times 10^{4}$ yrs, the system will be faced with the challenge of being unstable or being destroyed. Therefore, for planets undergoing a core-powered 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 $\sim~{}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 $\sim$ 5% of their total mass, the system may remain stable, even with a ”boil-off” phase.

We thank the referee for a thorough and constructive report that significantly improved the manuscript. This work is supported by the B-type Strategic Priority Program of the Chinese Academy of Sciences, Grant No. XDB41000000, National Natural Science Foundation of China (Grants No. 12033010, 11633009, 12111530175), the Natural Science Foundation of Jiangsu Province (Grant No. BK20221563), the Strategic Priority Research Program on Space Science of the Chinese Academy of Science (Grant No. XDA15020800), the China Manned Space Project with No. CMS-CSST-2021-B09, the Foundation of Minor Planets of Purple Mountain Observatory, and Youth Innovation Promotion Association.