Pushing planets into an inner cavity by a resonant chain
aa r X i v : . [ a s t r o - ph . E P ] F e b Astronomy & Astrophysicsmanuscript no. discedge © ESO 2021February 18, 2021
Can planets be pushed into a disc inner cavity by a resonantchain?
S. Ataiee , ⋆ and W. Kley Institut für Astronomie & Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany Department of Physics, Faculty of Sciences, Ferdowsi University of Mashhad, Mashhad, 91775-1436, IranFebruary 18, 2021
ABSTRACT
Context.
The orbital distribution of exoplanets indicates an accumulation of compact Super-Earth sized planetary systems close totheir host stars. Assuming an inward disc-driven migration scenario for their formation, these planets could have been stopped andeventually parked at an inner edge of the disc, or be pushed through the inner disc cavity by a resonant chain. This topic has not beenproperly and extensively studied.
Aims.
Using numerical simulations we investigate how much the inner planets in a resonant chain can be pushed into the disc innercavity by outer planets.
Methods.
We perform hydrodynamical and N-body simulations of planetary systems embedded in their nascent disc. The inner edgeof the disc is represented in two di ff erent ways, resembling either a dead zone (DZ) inner edge or a disc inner boundary (IB), wherethe main di ff erence lies in the steepness of the surface density profile. The innermost planet has always a mass of 10 M Earth , with equalor higher mass additional outer planets.
Results.
A steeper profile is able to stop a chain of planets more e ffi ciently than a shallower profile. The final configurations in our DZmodels are usually tighter than in their IB counterparts, and therefore more prone to instability. We derive analytical expressions forthe stopping conditions based on power equilibrium, and show that the final eccentricities result from torque equilibrium. For planetsin thinner discs, we found, for the first time, clear signs for overstable librations in the hydrodynamical simulations, leading to verycompact systems. We also found that the popular N-body simulations may overestimate the number of planets in the disc inner cavity. Key words.
Hydrodynamics - Methods: numerical-Planetary systems - Protoplanetary disks - Planet-disk interactions
1. Introduction
The notable number of short period Super-Earths is one ofthe puzzles introduced by the exoplanet discovery surveys (e.g.Winn & Fabrycky 2015; Dressing & Charbonneau 2015). To an-swer the question of how these planets could accumulate at thesub-au proximity to their host stars, several scenarios have beenproposed, such as: a combination of planet trapping at the mag-netospheric truncation radius of the disc and the star-planet tidalinteraction (Lee & Chiang 2017), planet-planet interaction aidedby the eccentricity damping of the gas disc (Chatterjee et al.2008; Ogihara & Ida 2009), planet trapping at the disc inner edge(Benítez-Llambay et al. 2011; Miranda & Lai 2018), and push-ing of the inner planets in a resonant chain by the outer planets inthe system (Coleman & Nelson 2016; Carrera et al. 2019). Thecommon core idea in all of these scenarios is the existence of adisc inner edge, where the planets can be trapped. This trappingcan either lead to the formation of a resonant chain that mightbecome unstable or push the inner planets into the disc innercavity, or simply stall the migration a single planet.Soon after the discovery of the first exoplanet, 51 Peg b, itwas suggested that this short-period planet has been stopped bythe inner edge of its natal disc (Lin et al. 1996). By the meansof two-dimensional (2D) hydrodynamical simulations, Massetet al. (2006) showed that the migration of a low-mass planetcan be halted at a steep surface density transition thanks to thechange in the Lindblad torque and enhancement of the positive ⋆ [email protected] co-rotation torque. In order to sustain a transition in the sur-face density, they increased the viscosity inside the inner discin their models. This viscosity transition can itself modifies theco-rotation torque. Without assuming a transition in the viscos-ity, Romanova et al. (2019) supported the trapping of low-massplanet at the disc inner edge by three-dimensional (3D) hydrody-namical simulations. In their study, the disc inner edge is mod-elled as a region with transition in density and temperature suchthat the pressure equilibrium is sustained. Therefore, the trap-ping of a single low-mass planet at the inner edge of the discseems to be a robust phenomenon. However, whether an inneredge could also stop the migration of multiple planets is stillquestionable.Trapping of multi-planet systems at the disc inner edge isoften observed in N-body simulations (e.g Izidoro et al. 2017;Raymond et al. 2018). In these simulations, the migration andeccentricity damping of the planets are modelled using the for-mulae that are obtained from the results of hydrodynamical sim-ulations for a single planet in a power-law disc. It has not beeninvestigated if those formulae is applicable close to the inneredge of a disc, where a sharp gradient in the surface density ex-ists. On the other hand, the trapping of multi planets in studiesthat performed hydrodynamical simulations has not been vigor-ously investigated. Papaloizou & Szuszkiewicz (2005) and Cuiet al. (2019) inspected the convergent migration of low-masstwo-planet systems in a disc with a surface density profile thatresembles the disc inner edge. The surface density in the innerpart of their disc rises linearly from the inner boundary up to a Article number, page 1 of 22 & Aproofs: manuscript no. discedge given distance and then becomes constant. In the outer part itbecomes a decreasing power-law. There is no report of pushingthe inner planet into the inner region of the disc except in onefigure of Papaloizou & Szuszkiewicz (2005), where the planetsbecome very close to the inner boundary. Therefore, the trap-ping of multi planets using hydrodynamical simulation needs tobe studied more vigorously.Trapping of multi-planet systems at the disc inner edge candi ff er from a single planet in two aspects. Firstly, trapped multiplanets are usually in resonance, that results into the excitation oftheir eccentricities (e.g Papaloizou & Szuszkiewicz 2005). Ec-centricity damping of a single planet and multi-planet systemhave been investigated in power-law discs (Cresswell & Nelson2006, 2008; Bitsch & Kley 2010), while it is not obvious if atrapped eccentric planet at the disc inner edge behaves similarly.Ogihara et al. (2010) used analytical calculations along with or-bital integration and found that an inner eccentric planet in aresonant chain, that is close to the inner edge of a disc, can haltthe migration of the whole chain. This type of trapping, namedeccentricity trap, arises from di ff erent eccentricity damping rateof the inner planet at the apocentre and pricentre. Although, theydid not consider the change of the torque on the inner planetcorresponding to the surface density jump, their study hints ondi ff erent behaviour of the eccentric planets at the disc edge. Sec-ondly, multi planets in a resonant chain migrate as one entity.Whether trapping the inner planet can stop the whole chain orwhether a massive enough system would overcome and breakthe trap need to be probed properly.The outcome of N-body simulations of multi planets at thedisc inner edge highly depends on what equation of motion isused or how the inner edge is modelled (Brasser et al. 2018; Idaet al. 2020). One advantage of hydrodynamical simulations isthat the migration of a planetary system is directly calculatedthrough the disc-planet and planet-planet gravitational interac-tions. Therefore, the main complexity would refer to the mod-elling of the inner edge.The detailed modelling of the inner edge of a disc, where thematerial is transferred to the central star, can be very complicated(Romanova et al. 2002; Long et al. 2005; Romanova & Lovelace2006; Zanni & Ferreira 2009). However, we can simply defineit as where the gas leaves the disc and settles to the star. Nearto the inner edge of the disc can be other traps as well such as adead-zone edge. Flock et al. (2017) studied the inner region ofprotoplanetary discs around solar-like stars using radiative mag-netohydrodynamical simulations and show that the silicate sub-limation at ∼ .
08 au leads to the formation of a density jumpat ∼ . ff er in their steepness and long term evolution. InSec. 2, we describe how these two types of traps are modelled,which code is used, and how the disc and planets are set up. We performed simulations with two, three, and more planets withdi ff erent mass combinations, and also for a lower disc viscosityand smaller aspect ratio. The results are listed in Sec. 3. Afterthat, in Sec. 4, we check if the results depend on the width of ourdead-zone edge model. A power and torque analysis is presentedin Sec. 5, where we show the halting of a resonant system hap-pens where the total power of the system vanishes. In some ofour models, we observe overstability, that is discussed in Sec. 7.Finally, our findings are summarised in Sec. 9.
2. Modelling the inner edge: setup and numerics
As described in the introduction it is expected that a protoplane-tary disc disrupts near the stellar surface and has an inner edge,where the disc surface density drops significantly. Another loca-tion that is also close to the inner edge and has transition in sur-face density is the dead-zone edge. Both of them can trap planetsif the disc provides a strong enough co-rotation torque. However,the planets have to overpass the trap to get closer to the star. Us-ing 2D hydrodynamical simulations, we study the migration ofplanets and their trapping at these two types of density transi-tion. We analyse the final configuration with respect to their res-onances and monitor if the inner planet can be pushed into theinner cavity.Bringing the planets into a resonance, allowing the innerplanet to reach the trap, and properly modelling the migra-tion of planets at the disc inner edge is cumbersome and time-consuming. Even if a planetary system is initialised close to aresonance commensurability, it might need thousands of orbitsto reach a steady final state. On the other hand, properly resolv-ing the co-rotation torque for a low or moderate mass planetneeds at least six cells in a half horseshoe width of the planet(Paardekooper et al. 2011). Having such high resolution at thedistances very close to the star would decrease the time-step sodramatically that the simulations become infeasible. To evadethis issue, we shift the density transition to ∼ We construct two types of traps: a dead-zone inner edge (DZ)and a disc inner boundary (IB). To ease understanding of the re-sults, we stick to this colour code for the rest of the manuscriptwhen describing these two cases. These planetary traps are con-structed in the following way:Dead-zone inner edge: We model a dead-zone inner edgeby decreasing the viscosity over a distance of the order of discscale height, using the the method of Masset et al. (2006). Inthis model the inner part has a higher viscosity while the outerpart has a lower value. As a consequence the inflow velocity islarger inside the cavity and smaller outside, that gives rise to adensity maximum just outside of the viscosity transition region.The enforced viscosity determines the shape of the density jump.Because the viscosity profile is fixed, it leads to a surface den-sity profile that does not evolve much by time. This fixed surfacedensity profile is the advantage of this model. However, as itturned out, the strong viscosity variation around the location ofthe planetary trap and consequent vortex formation complicatesthe saturation of the co-rotation torque.Disc inner boundary: We assume that the region between thedisc edge and the star is emptied by some mechanism. Hence,when the gas reaches to the inner edge, it is taken out from the
Article number, page 2 of 22. Ataiee & W. Kley: Trapping planets at a disc inner edge simulations such that an extremely low surface density is createdat the inner boundary. The inner disc evolves correspondinglyand adopts its profile to this inner boundary condition. A sim-ilar method has been used by Romanova et al. (2019) to studytrapping of a planet in 3D hydro-simulations. The drawback ofthis model is that this boundary condition causes the surface den-sity profile, which determines both the Lindblad and co-rotationtorques, slowly evolves in time. However, the smooth viscositycan be considered as its advantage.To illustrate the di ff erent behaviours of these two traps, wedisplay the evolution of the surface density and the specifictorque on an imaginary 10 Earth-mass planet in Fig. 1. Here andin subsequent plots, the left panels (green frames) refer to theDead-zone inner edge (DZ) and the right panels (purple frames)to the Disc inner boundary (IB) model. The top row shows thatthe DZ model yield a very steep surface density transition thatbecomes stationary with time. In the bottom row, 2D maps ofthe specific torque, Γ , acting on a fictitious 10 Earth-mass planetat every location in the disc at a specific time during the disc’sevolution are displayed. The torque is calculated using the for-mulae in Paardekooper et al. (2010a, 2011). The positive den-sity slope at the inner disc’s edge produces a region where thetorque acting on the planet becomes positive. The black con-tours indicate where the torque is zero. The actual planetary trapis at the outer edge of the positive torque region. At the trap,the planet feels a negative torque if it is shifted outwards and apositive torque if slightly displaced inwards. In the DZ model,there is a sharp transition from the negative to positive torque atthe trap, whereas in the IB model, the gradient of the torque isvery shallow around the zero-torque location. Said di ff erently, ifthe planet is pushed slightly inwards, it would feel a very strongpositive torque in the DZ model and is sent back to the zero-torque location, while in the IB model, the planet would feel alittle positive torque which can be easily overcome by the pushof an outer planet. Therefore, based on the torque map of our discmodels, we anticipate that a resonant chain of planets would bemore successful in push a planet out of the planetary trap in theIB than the DZ model. In the next section, we will explain inmore detail the numerical setups of these two models. We used a modified version of the 2D hydrodynamical code
FARGO . This code simulates a locally isothermal disc in thecylindrical coordinate. We used a fixed coordinate frame cen-tred on the star and took into account the indirect terms both onthe planets and the disc.The aspect ratio of the disc is h = h ( r / r ) . with r = h is 0.05 for the reference model,but all parameter sequences include models with the smallervalue of 0.03. Our time unit, named orbit, defined as one Ke-plerian orbital time at r . To be initially in viscous equilibrium,we chose Σ = Σ ( r / r ) − . [M ⋆ / r ] with Σ = × − except forfew simulations which have double of this value. Σ = × − corresponds to ∼ / cm for a Solar-mass star M ⋆ = ⊙ and r = au. We note that this is the basic profile in the discthat is modified the edge at the inner part of the disc for eachdisc model as we explain later. Given this setup, the disc inneredge in our models will be located near 1 au, while in real discsit is approximately 10 times smaller. Our choice was done fornumerical reasons, because a smaller radius would require manymore time steps due to the shorter dynamical timescale, not al- http: // fargo.in2p3.fr / -Legacy-archive- −2 −1 Σ / Σ r [r ] t [ o r b i t ] Γ [ × − ] t = 0t = 2500t = 5000t = 7500t = 10000 r [r ] −1.0−0.50.00.51.0 Γ [ × − ] Fig. 1.
The evolution of the disc without any embedded planets forthe two types of trap dead-zone inner edge (left) and disc inner bound-ary (right).
Top : The 1D surface density evolution. The colour of linesrepresents the time in units of the orbital period at r =
1. Black denotesthe initial profile and lines pale as time proceeds.
Bottom : The analyticaltype I torque acting on a 10 Earth-mass planet calculated using the discprofiles. Every horizontal line in these plots shows the specific torque atevery distance from the star at a given time. The black contours repre-sent the zero-torque locations. lowing an extensive parameter study. As we deal with locallyisothermal discs the results should be scalable to smaller radii,using the value h = .
03 (Flock et al. 2017).For the disc viscosity ν , we use the α -viscosity model ν = α c s H , where c s = H / Ω is the sound speed with Ω being theKeplerian angular velocity, and H = hr is the disc scale height.In most of the simulations α = × − except in few models, inwhich the e ff ect of a ten times smaller viscosity is examined.In the DZ model, the viscosity becomes a hundred timeslarger inside the active zone. The transition happens over a dis-tance of λ = . H ( r t ) at r t = α = × − × r < r r − rr − r r ≤ r ≤ r r > r (1)with r = r t − λ/ r = r t + λ/
2. For the IB model, α isconstant through the whole discIn the DZ models, the disc is meshed into N r =
929 loga-rithmic segments in the radial, and N φ = r in = .
3. We set our outer boundary far enough at r out = r =
18 to 20, using the method in de Val-Borro et al.(2006).In the IB models, the outer boundary is identical to theDZ model. The inner boundary is located at r in = . r =
1. We changed the ra-dial resolution to N r =
712 while keeping the azimuthal grid thesame as the DZ model to have similar resolution in both models.Following the explanation in Sec. 2.1, we successively decrease
Article number, page 3 of 22 & Aproofs: manuscript no. discedge the surface density within the interval r ∈ [0 . − .
9] until itreaches a floor value of Σ floor = − . The radial velocity in thisregion is outflow. This slow decreasing of surface density pre-vents producing an instability by guaranteeing that the surfacedensity has enough time to be adopted to the new condition. Wewould like to point out that any other choice of inner boundary,either location or damping width, could produce a di ff erent pro-file and consequently dissimilar incommensurability than whatwe present here.The inner planet is planted at r = . r sm = ǫ H ( r ) with ǫ = . M i =
10 M ⊕ where M ⊕ = × − M ⋆ is the mass of Earth for a Solar-like star. The rest of the planetsare more massive than or as massive as the inner one. The maxi-mum number of the planets used in our models is 6 (See Table 1for the details).
3. Results
In this section we summarize the main results of our study. Toanalyse the impact of an inner edge on the final configuration ofembedded planets, we ran an extensive suite of simulation usingvarious disc parameter and planet masses. The set of models issummarized in Table. 1, and we refer to the quoted model identi-fiers in describing their outcome. The simulations are continueduntil all planet pairs settle into resonances or the system becomesunstable. As we found that the result of each simulation is indi-vidual, we present and explain them one by one in the following.To make the whole text more readable we display additional re-sults on the individual models in appendix B.
In the fist step, we checked trapping of a single planet in our discmodels. The top panels of Fig. 2 show the migration of a 10 M ⊕ (blue) and a Jupiter-mass planet (orange) in a dead-zone inneredge (left) and a disc inner boundary (right) model (1p10 and1p1Jup in Tab. 1). As pointed out, the colour of axes specifies thetype of disc model. In the bottom rows, the surface density pro-files are presented at t =
50 orbits (in black), when the planet’spotential is fully applied on the disc, and t = orbits.In the DZ model (with steep density slope), the low-massplanet is trapped at the zero torque location ( r ∼ .
14) and itssemi-major axis does not evolve, as the torque map in Fig. 1also suggests. The small oscillations in the semi-major axis isdue to the presence of some vortices outside of the edge. Thesevortices exert some torques on the planet as they pass by andproduce small wiggles. In contrast, the massive planet continuesits inward migration as it opens a planetary gap and deformsthe density profile which completely suppressing the positive co-rotation torque. The surface density in each timestep is multiplied by a factor of 0.999 t [orbit] a [ r ] t [orbit] r [r ] −2 Σ / Σ r [r ] Fig. 2.
Evolution of a single planet in a disc with a dead-zone inneredge (left) and at the disc inner boundary (right).
Top : Migration path, a ( t ). In each panel, every line belongs to a separate simulations withonly one planet. Bottom : The surface density profiles at t =
50 (black),when the full potential of the planet is applied on the disc, and t = orbit. The colour code is the same as for the top row and the finallocation of the planets are marked by circles. The low mass planet istrapped the DZ model and slowly migrates outwards in the IB model.In contrast, the massive planet (Jupiter-mass) continues its migration tothe inner disc thanks to the deformation of the bump by its gap. a i n [ r ] t [orbit] a i n [ r ] Fig. 3.
Semi-major axis evolution of the inner 10 M ⊕ planet for all mod-els with 2 planets, top for DZ models and bottom for the IB ones. Thesingle planet model (1p10) is added for comparison. In the IB model, the low mass planet moves slowly outwardas the trap regions widens and the zero torque location movesslowly outwards. The massive planet in this model behaves sim-ilar to the DZ case and continues its inward migration by re-shaping the surface density profile.These results imply that massive gap-opening planets couldpossibly push and shepherd small planets that are trapped at theedge toward the star. In the next section, we add a second planetwith various masses to these models and also examine the e ff ectof lower viscosity and smaller aspect ratio. Article number, page 4 of 22. Ataiee & W. Kley: Trapping planets at a disc inner edge h Σ α r p M p [M ⊕ ] Result1 1p10 0.05 2 × − × − × − × − × − × − × − × − × − × − ∗ × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − × − > × − × − × − × − Table 1.
The list of our main simulations. The information about the DZ and IB models are given by their corresponding colours. The data inblack are the same in both trap models. ∗ The planets are very close to 29:20 commensurability with T o / T i ≃ . In this section, we add a second planet just o ff the 2:1 commen-surability and monitor the evolution of the system. Mass of theinner planet in all cases is 10 M ⊕ but the outer one has di ff er-ent masses. We denote the inner planet by an index ’i’ and theouter by ’o’. Taking the case with M o =
20 M ⊕ (model 2p20)as the reference model, we additionally examined the e ff ect oflower viscosity, smaller aspect ratio, and higher surface densityfor this model.Figure. 3 embodies the migration of inner planet for all mod-els with two planets. The top panel contains the results for theDZ model (green frame) and the bottom panel (purple frame)demonstrate those of the IB model. The figure clearly showsthat the second planet undoubtedly has an a ff ect on the migra-tion of the inner planet specially if it is more massive than theinner one. Moreover, comparing the top and the bottom panelsimplies that it is easier to push the inner planet to the inner discin the IB models than the DZ models. In the DZ models, theinner planet insists on staying at the trapping location while inthe IB models, the outer planet in most cases is able to reversethe migration of the inner planet and pushes it inward after theyget into resonance. In the following, we describe the results ofthe models individually in more details. : In these models, the masses of both planets are 10 M ⊕ .As Fig. 4 shows, in the DZ model, the planets cross severalresonances and finally settle in the 5:4. The blue line in the toppanel of Fig. 3 indicates that the inner planet remains almost inits initial trapping location. In this case, the push by the outerplanet is not strong enough to overcome the positive torque onthe inner planet. On the contrary, in the IB model, after theplanets get into the 3:2 resonance, the outer planet pushes theinner planet inwards and reverses its slow outward migration. : Figure 5 shows semi-major axis, eccentricity, and orbitalperiod ratio for 2p20 models. In the DZ simulation, the outerplanet migrates inwards and gets into the 4:3 resonance with thetrapped inner planet. Because the outer planet is twice as mas-sive as the 2p10 model, the inner planet is pushed slightly intothe transition zone (top panel of Fig. 3). After forming the finalresonance, the migration of the planets stops. On the contrary,in the IB simulation, after the outer planet catches the outwardsmigrating inner planet in the 3:2 resonance, they continue theirmigration inwards but with a slower rate than the one of outerplanet before the resonance. : Fig. 6 shows the results for an outer planet of 100 M ⊕ .In addition to the previous figures, the bottom panel displays thesurface density of the disc at the end of simulations. Locations ofplanets are indicated by markers. The outer planet in these mod- Article number, page 5 of 22 & Aproofs: manuscript no. discedge a [ r ] e t [orbit] T o / T i t [orbit] M i = 10M ⊕ , M o = 10M ⊕ Fig. 4.
Evolution of semi-major axis, eccentricity, and orbital periodratio for 2p10 models. The inner planet is given by the blue line and theouter by the orange one. In the bottom panels resonances are marked byhorizontal lines. a [ r ] e t [orbit] T o / T i t [orbit] M i = 10M ⊕ , M o = 20M ⊕ Fig. 5.
Similar to Fig. 4 but for the 2p20 model. els opens a partial gap, not as strong as the one of the Jupitermass planet in 1pJup to destroy the bump, but can still a ff ectit. The left panels in this figure, that contain the results for theDZ model, show that the outer planet cages the inner planet be-tween the edges of the dead-zone and its gap. After temporar-ily settling into the 3:2 MMR, they rearrange and settle neara 29:20 commensurability. We checked the resonant argumentsbut found no definite pattern. In the IB model (right panels), theplanets get into the 2:1 resonance and keep this configuration un-til the end of the simulation. Both planets continue their inwardmigration and the inner planet, which becomes very eccentricwith the eccentricity of ∼ .
4, is pushed to the inner part of thedisc. : In these models the outer planet has one Jupiter massand opens a gap similar to the 1pJup model. For both edge mod-els, the inner planet traps into the 2:1 resonance with the outerone (see Fig. B.1 in the appendix B). The massive planet pushes a [ r ] e t [orbit] T o / T i t [orbit] r [r ] −2 −1 Σ / Σ r [r ] M i = 10M ⊕ , M o = 100M ⊕ Fig. 6.
The three top rows are the same as Fig. 4 but for the 2p100model. The forth row shows the final surface density with the locationsof planets indicated by markers, whose colours correspond to the toppanels. the inner planet into the inner part of the disc regardless of whichinner disc model is used. The low-mass planet in both modelsbecomes very eccentric. We stopped the simulations because forsuch large eccentricities ∼ .
3, the inner planet passes the areathat is close to the inner boundary of our computational domain.The overall behaviour of the planets in both of these models arevery similar and it seems that the type of the trap has only littleinfluence on the evolution. Therefore, we will not examine moremodels with a Jupiter-mass planet. : These models have identical setups as the reference2p20 models except that α is lowered to 5 × − . The lower vis-cosity can ease the partial gap opening and consequently reduceboth of the Lindblad and co-rotation torques. In addition, it canchange the saturation of the co-rotation torque (Masset 2001;Kley & Nelson 2012). In the bottom panels of Fig. 7 the sur-face density profiles at the end of the simulations are displayed.We can distinguish a shallow gap around the outer planet withthe mass of 20 M ⊕ . Comparing the two upper rows of Fig.7 withthose of Fig.5 shows that the planets in the DZ model man-aged to get closer and settle into the 6:5 resonance. However,they finally closely interact and the inner planet is sent to an or-bit outside of the bump, and the order of the planet is reversed.A comparison with the top panel in Fig. 3 shows that the innerplanet here is initially trapped more inside than the planet in the2p20 and 1p10 models with higher viscosity. This is consistentwith the idea that the lower viscosity reduces the positive co-rotation torque due to the saturation e ff ects.In the IB model, similar to its 2p20 counterpart, the planetspass the 2:1 resonance (around ∼ Article number, page 6 of 22. Ataiee & W. Kley: Trapping planets at a disc inner edge a [ r ] e t [orbit] T o / T i t [orbit] r [r ] −4 −2 Σ / Σ r [r ] M i = 10M ⊕ , M o = 20M ⊕ , α = 5 × 10 −5 Fig. 7.
Similar to Fig. 6 but for 2p20LA models which have a reducedviscosity. major axis of the inner planet (lower panel of Fig. 3) shows thatin contrast to the 2p20 model, the inner planet initially is migrat-ing inwards which indicates that the co-rotation torque is moresaturated. After the outer planet passes the 2:1 resonance, the in-ner planet changes its migration direction, moves outwards, andhalts until the outer planets catch it in the 3:2 resonance. There-after, its migration is governed by the outer planet. : These models have a smaller aspect ratio than2p20, here h = .
03. A smaller aspect ratio can change the torquein two ways. Firstly the torque from the disc on the planet is pro-portional to h − (Baruteau et al. 2014, and references therein),and therefore the torque is larger for smaller h . Secondly, theplanets are liable to open a partial gap and the smaller surfacedensity around the planet’s orbit decreases the torque. The brownlines in Fig. 3 shows that the semi-major axes of the inner planetsfor the 2p20LH models are smaller than those of the referencemodel 2p20.The left panels of Fig. 8 show that the planets in theDZ model finally settle into the 4:3 resonance and remain there.In the IB model, whose results are displayed in the right pan-els of Fig. 8, the inner planet migrates inwards until outer planetcaptures it in the 3:2 resonance at around t ≈ : The surface density in these models is twice thatof the 2p20 models. Hence, the planets feel larger torques com- a [ r ] e t [orbit] T o / T i t [orbit] r [r ] −2 −1 Σ / Σ t [orbit] r [r ] t [orbit] M i = 10M ⊕ , M o = 20M ⊕ , h = 0.03 Fig. 8.
Similar to Fig. 7 but for 2p20LH except that the panels of thelast row show three di ff erent times to demonstrate the gap opening ofthe planets. The markers display the location of each planet with theircolours correspond to the upper panels. pared to 2p20. The migrations and eccentricities are almost iden-tical to the 2p20 model, except the evolution is faster. It canbe seen in the lower panel of Fig. 3 that in the IB modelthe behaviour of the inner planet is very similar in the 2p20and 2p20HS models but the resonance capture happens earlierin 2p20HS. However, the results are slightly di ff erent for theDZ model. The planets are initially get into 5:4 resonance, thenare pushed into and stay in the deeper 6:4 resonance for about10 000 orbits. After that, they return back to 5:4 because of theplanet-planet gravitational interaction (see Fig. B.2 in the ap-pendix B). Despite changing of the resonance configuration, thesemi-major axis of the inner planet does not vary and remainsthe same as in the 2p20 model after the first resonance (upperpanel of Fig. 3). In this section, we will present the results of the simulations withthree planets. Mass of the outer planet is either 20 or 100 M ⊕ ,the middle 10 or 20 M ⊕ , and the inner planet is always 10 M ⊕ ,see Tab. 1 for details. The outer planet in these models is addedjust outside of 2:1 commensurability with the middle planet. Fig-ure 9 shows the migration of the inner planet for all models withthree planets. In the top panel, which shows the results for theDZ models, the inner planet in many simulations undergoes aninstability and is expelled. However, none of the inner planetsis ejected in the IB model (the bottom panel of Fig. 9). In thefollowing, we describe the results of each model individually. : These simulations have the same setup as 2p10 ex-cept that a third planet with the mass of 20 M ⊕ has been added. Article number, page 7 of 22 & Aproofs: manuscript no. discedge a i n [ r ] t [orbit] a i n [ r ] Fig. 9.
Semimajor axis evolution of the innermost planet. Similar toFig. 3 but for models with three planets. The inner planet has always10 M ⊕ and the labels give the masses of the other two in units of M ⊕ ,see Tab. 1. a [ r ] e t [orbit] T o / T i t [orbit] M i = 10M ⊕ , M o = 10, 20M ⊕ Fig. 10.
Top and middle panels display the semi-major axis and eccen-tricity for 3p1020 models. The bottom panels show the orbital periodratio for each consecutive pair of planets, where the colour correspondsto the outer planet. For example, the green curve demonstrates the or-bital period ratio of the third planet to the second planet when countingfrom inside out.
With the aid of the black lines in Fig. 3 and 9, that show the semi-major axis of a 10M ⊕ planet in the single-planet models, we cancompare the e ff ect of adding the third planet to the evolution ofthe inner planet. Comparing the blue lines in the upper panels ofthe mentioned figures indicates that, in the DZ model, the pres-ence of the third planet hardly a ff ects the migration of the innerplanet. When the third planet catches the other two in a 6:5:3resonant chain (with period ratios 6:5 and 5:3 from inside out)at ∼ a [ r ] e t [orbit] T o / T i t [orbit] r [r ] −2 −1 Σ / Σ t [orbit] r [r ] t [orbit] M i = 10M ⊕ , M o = 10, 100M ⊕ Fig. 11.
Similar to Fig. 10 but for 3p10100 model. the DZ model, the outer planet pushes the system into a tighterconfiguration, which is more prone to instability.In the IB model, the outer planet gets into 2:1 resonance withthe middle one but it has no e ff ect on the migration of the innerplanet (right panels of Fig. 10). As the bottom panel of Fig.9shows, the semi-major axis of the inner planet follows that ofthe single planet until it gets into 3:2 resonance with the middleplanet. Afterwards, it is pushed into the inner disc with a fastermigration rate than its 2p10 counterpart. : The trapping positions of the inner planets in thesemodels are very similar to their two-planet counterparts and theouter planet does not play a notable role. In the DZ model, al-though the outer planet does not push the most inner one, it candestabilise the system after it enters into 2:1 resonance with themiddle planet (Fig. B.3). After the two outer planets spend sometime in this resonance, the inner and middle planets switch theirorbits. In the IB model the middle planet enters 3:2 resonancewith the inner one at around 8000 orbits, and it causes the ec-centricity of the inner planet increases, exactly as in the 2p20model, see Fig. 5. Later, the outer planet catches up and a 3:2:1resonant chain forms with period ratios 3:2 and 2:1 (from insideout). Upon the capture, the eccentricities of the two inner plan-ets increase further. Otherwise no change in the migration rateof the inner planet happens. : These models are complementary to the 2p10 and3p1020 models in understanding the e ff ect of the third planet.Comparing the DZ simulations in the three mentioned models,shown in the left panels of Fig. 4, 10, and 11, indicates that thetwo inner planets form tight configurations in all three modelsand the role of the outer planet is mostly destabilising the sys-tem. The more massive the outer planet is, the earlier the systembecomes unstable. Di ff erently, in the IB simulations, the thirdplanet has a notable impact. When the third planet is more mas- Article number, page 8 of 22. Ataiee & W. Kley: Trapping planets at a disc inner edge a [ r ] e t [orbit] T o / T i t [orbit] M i = 10M ⊕ , M o = 20, 20M ⊕ , h = 0.03 Fig. 12.
Evolution of eccentricity and orbital period ratio for 3p2020LHwhich has a smaller disc thickness. sive, the reversal of the inner planet’s migration happens ear-lier, and the inward migration is faster after the reversal. In allof these IB models, the two inner planets are in 3:2 resonanceand when a third planet exists, the planets form a 3:2:1 resonantchain, which enhances the eccentricity of the lower mass planets.The last row of Fig. 11 shows the azimuthally averaged sur-face density of the 3p10100 models at given times. In both cases,the outer planet opens a gap. In the DZ model, the two innerplanets are sandwiched between the dead-zone edge and the par-tial gap of the outer planet. In contrast to the DZ model thatbecome unstable, the planets in the IB model migrate inwardswhile keeping their resonant chain. : Evolution of these models are very similar to3p2020 but faster due to the faster migration of the more massiveouter planet. Similar to 3p2020, the DZ model becomes unstableafter about 10 000 orbits. The planets in the IB model get into a3:2:1 chain. The only di ff erence is that the eccentricities of theplanets in this configuration are larger than those of 3p2020, dueto the larger mass of the outer planet. Looking at the migration ofthe inner planet for this model and 3p10100 in the lower panelof Fig. 9 indicates that the behaviour of the inner planet afterreversal is determined completely by the third planet which ismore massive than the second one. (The eccentricity and orbitalperiod ratio for these models are brought in Fig B.4.) : These models with a smaller α behave very sim-ilar to their high viscosity counterparts (as shown in Fig. B.5).The DZ model becomes unstable when the outer planet getsclose to the inner two, which are in 6:5 resonance. The planetsin the IB model, similar to the previous IB models with threeplanets, continue their inward migration as they stay in the 3:2:1resonant chain. : These models, with a smaller aspect ratio,show di ff erent behaviours than other three-planet models. TheDZ model here is the only one among other DZ models withthree planets that does not go unstable after being ina 4:3:2chain for more than 5000 orbits. Because the inner planet barelymoved from its trapping location (upper panel of Fig.9), the mi-gration of whole system halts as they form the resonant chain.One can see this in the left panels of Fig. 12. In the IB model, asthe lower panel of Fig.9 shows, the inner planet passes various a [ r ] e t [orbit] T o / T i t [orbit] a i n [ r ] M i = 10M ⊕ , M o = 20, 20, 20M ⊕ , DZ Fig. 13.
Results of onebyone simulation for the DZ model. The toprow displays the semi-major axes (left) and eccentricities (right) of allplanets. The bottom left panel contains the orbital period ratio for eachneighbouring pair. In the bottom right, the migration of the innermostplanet is shown for more clarity. phases of inward and outward migration depending on in whichresonance it is with the outer planet. Even after the final con-figuration forms, it does not migrate inwards. The right panelsof Fig.12 demonstrate the evolution of this system. The two in-ner planets get into various resonances, namely 3:2, 4:3, and 9:7,but these resonances are broken by the overstability. They finallydwell in 5:4 resonance. The outer planet, however, remains in the3:2 resonance with the middle planet after they enter into it. : The outcome of these models with higher sur-face density and 3p2020 models are alike, except that the oneshere have faster evolution. Therefore, we omit showing the fullresults. We can see in Fig. 9 that the inner planets in these mod-els behave very similar to their low surface density counterpartsbut their inward movements have started earlier. As before, theIB model ends up in the 3:2:1 chain and the DZ model becomesfinally unstable.
In this section, we examine the e ff ect of additional planets in thesystems and its e ff ect on the innermost planet. We take the 2p20as the reference model, add more planets with the mass of 20 M ⊕ ,and observe the behaviour of the innermost planet. We start withthe inner 10 M ⊕ planet and allow it to find its trapping location.Then the second planet with the mass of 20 M ⊕ is placed justout of 2:1 commensurability. We continue the simulation untilthe planets get trapped into a resonance. The next planets areadded one by one in a similar way when the inner pair form aresonance. We continue until the system becomes unstable orwe find that adding more planets has no e ff ect on the location ofthe inner planet. In the following the results of these simulationsfor discs with h = .
05 and h = .
03 for both inner edge modelsare presented.
Onebyone : The evolution of the system for the DZ modelis shown in Fig. 13. On the top left, the semi-major axes of theplanets are plotted. This panel, besides the migration of the plan-ets, demonstrate where and when the planets are added into thedisc. As the top right panel shows, the eccentricities are about0 .
07 unless when the planets closely interact. By the help ofthe bottom panels we can check how the location of the innerplanet is a ff ected by the outer planets. The second planet (orange Article number, page 9 of 22 & Aproofs: manuscript no. discedge a [ r ] e t [orbit] T o / T i t [orbit] a i n [ r ] M i = 10M ⊕ , M o = 20, 20, 20M ⊕ , IB Fig. 14.
Similar to Fig. 13 but for onebyone IB model. a [ r ] e t [orbit] T o / T i t [orbit] a i n [ r ] M i = 10M ⊕ , M o = 20, 20, 20, 20, 20M ⊕ , h = 0.03, DZ Fig. 15.
Similar to Fig. 13 but for onebyoneLH model which has a loweraspect ratio. curve) migrates inwards until it gets into 3:2 resonance with thefirst one. The bottom right panel reveals that the inner planet ishardly pushed inwards at this time ( t ∼
14 000 orbit). Shortly af-ter the third planet (green) forms 3:2 resonance with the secondplanet, the two inner planets go into the tighter 4:3 resonance.The forth planet (red) remains in 2:1 resonance with the thirdplanet and the whole system forms a 4:3:2:1 chain. As long asthe system abide in this configuration, the inner planet is mov-ing slowly inwards. But then the two inner planets go into thetighter 5:4 resonance and the system becomes unstable thence-forth. This instability leads to sending the inner planet out of thetransition zone into the inner disc.The IB model, unlike the DZ simulation, evolves verysmoothly. The two inner pairs get in 2:1 resonance and remain init. Thereafter, the migration of first planet reverses and becomesinwards. Later, its migration rate is not influenced by the pres-ence of the other planets. The third planet, however, ignites theeccentricities of the first and second planet after it get into 2:1resonance with the second planet. The last planet seems to haveno major e ff ect on the system. We did not continue this simula-tion because of the high eccentricity of the inner planet. OnebyoneLH : In these simulations the smaller disc aspectratio changes the migration rate compared to the two previouscases. Unlike the onebyone models, the DZ model here remainsstable even at the presence of six planets, while the IB model be-comes unstable after fourth planet. As Fig. 15 shows, we added a [ r ] e t [orbit] T o / T i t [orbit] a i n [ r ] M i = 10M ⊕ , M o = 20, 20, 20M ⊕ , h = 0.03, IB Fig. 16.
Same as Fig. 14 but for IB edge model. five 20 M ⊕ planets to the system in the DZ model. All of themstay in either 3:2 or 4:3 resonances but the inner planet is notpushed inwards. The only e ff ect of the outer planets on the firstplanet is producing some jitters in its semi-major axis by enhanc-ing its eccentricity.The IB model, on the other hand, as Fig. 16 demonstrates,becomes unstable after the forth planet comes into resonance bythe third one and forms a 5:4:3:2 resonant chain. From the topright panel, we can infer that the two inner pairs undergo over-stability. The inner planet shows absolutely no monotonic be-haviour. Depending on which resonance it forms with the otherplanets, it can move inwards or outwards.
4. Transition width and vortices
For our choice of the viscosity transition in the DZ models,some vortices are created and persist during the simulations.They form because of the narrow transition width and a smalldiscontinuity in α at the edges of the transition zone. In this sec-tion, we examine di ff erent widths and profiles in order to inspectwhether the vortices and transition width can a ff ect the migra-tion and final settling of the planets. For the models 2p20 and3p2020, we ran additional simulations using di ff erent viscosityprofiles, either by changing λ in Eq. (1), or by smoothing thetransition.The top panel of Fig. 17 shows four azimuthally averagedsurface density profiles for models with viscosity profile givenby relation (1) using λ = . , . , .
0, and one which has thesame transition width as λ = . λ = .
06 (our reference model),for which the outer planet could get closer to the inner planet.The corresponding two-dimensional surface densities aredisplayed in Fig. 18 at the time just before the planets getinto resonance. In the three simulations with our standard pro-file, there is at least one vortex visible, while the one with thesmoother transition shows none. In all of these simulations, re-gardless of the existence of a vortex / vortices, the inner planet istrapped at the edge of the transition zone. In the model with thesteepest transition, λ = .
06, there are many small vortices atthe transition edge. Although they can hardly be distinguishedin Fig. 18, their spirals are clearly visible (Paardekooper et al.2010b). Both the vortices and their spirals can interact with the
Article number, page 10 of 22. Ataiee & W. Kley: Trapping planets at a disc inner edge r [r ] −2 −1 Σ / Σ λ = 0.06λ = 0.3λ = 1.0λ = 0.3, smooth t [orbit] T o / T i Fig. 17.
Results for 2p20 model using di ff erent transition widths for theDZ case. Top : Azimuthally averaged surface density of the inner partof the discs with di ff erent transition widths, λ , see eq. (1). The profilelabelled with λ = . smooth has a transition width as λ = . Bottom : Orbital period ratio for the models in the top panelplotted with the corresponding colours. planets and produce small perturbation in their orbits. However,because these interactions are random, their torques / powers onthe planets do not influence the overall migration. For the mod-els with the wider transitions λ = . λ =
1, only a singlevortex is established and the inner planet is trapped azimuthallybetween the two ends of the vortex (Ataiee et al. 2014). In allof these models, after the planets get into resonance and theireccentricities are excited, the vortex / vortices fade.In the model with the smooth transition (labelled with λ = . λ = . λ = . . λ = .
06, wherethe system goes into 5:4 resonance and becomes unstable, theinner planet in the simulation with λ = . λ =
1) the outer planetpair managed to push the inner planet out of the trapping point(Fig. B.6). Therefore, it appears that, similar to the IB models,the inner planet can be pushed into the inner disc if the transitionzone in the DZ models are extremely wide.
5. Power and torque analysis of 2p20 models
The results of two-planet models show that, after the resonanceformation, the inner planet in the DZ simulations is slightlypushed inwards, where the torque is expected to be positive. Af-ter that the migration of the whole system halts. In the IB mod-els, in spite of the initial outward migration of the inner planet,the resonant system pushes the inner planet back to the innerdisc and the whole system migrates inwards with a di ff erent mi-gration rate than those of the individual planets. In this section r [ r ] λ = 0.06 λ = 0.3 θ r [ r ] λ = 1.0 θ λ = 0.3, smooth Σ/Σ Fig. 18.
Two-dimensional surface density for the discs in Fig. 17. Plan-ets are marked by white circles and the arrows show the centre of the(largest) vortex in the models with a vortex / vortices. −202 Γ d i s c + Γ i n d inner outer sum t [orbit] P d i s c + P i n d t [orbit] Fig. 19.
The left column shows the torque (top) and power (bottom)on a single planet of 1p10 DZ model. On the right, similar quanti-ties are plotted for 2p20 model. The dark green line represents the totaltorque / power on the system. The black horizontal lines in each panelmarks zero. We note that the contribution of the indirect term is includedin the plotted quantities. Torques and powers are averaged over 1000 or-bits to decrease the oscillations created by the vortices and eccentricitiesof the planets. The dashed lines in the right panels correspond to the val-ues of the eccentricity and semi-major axis damping timescales used inthe supporting N-body simulations, as displayed in Fig. 20. we analyse the torque and power for 2p20 model (displayed inFig. 5), and show that the halting position of the whole system isat that location where the total power vanishes.The energy E , and angular momentum L , of a planet is givenin terms of its semi-major axis and eccentricity as E = − GM ⋆ M p a and L = M p p GM ⋆ a (1 − e ) . (2)As the energy depends only on the semi-major axis of the planet,its migration rate is determined by the disc power, P disc , acting Article number, page 11 of 22 & Aproofs: manuscript no. discedge on it. The power gives the total energy change of a planet˙ E = P disc . (3)On the other hand, the disc torque, Γ disc determines the rate ofangular momentum change of a planet˙ L = Γ disc , (4)which is a combination of eccentricity and semi-major axischange (see e.g. Cresswell et al. 2007). For a single planet ona circular orbit, torque and power are equivalent. In this case, itis expected that the migration of planet halts where the torquevanishes. However, for a planet on an eccentric orbit, the haltinglocation will be given by vanishing power. In order to analysethe halting process of a planet pair, we monitored the disc torqueand power continuously during the evolution of the system . Fig-ure 19 shows the torque (top) and power (bottom) acting on theplanets as a function of time for a single planet model (1p10) onthe left and the reference two-planet model (2p20) on the right.In these plots, the contribution of the indirect term has been takeninto account in the calculation of the torque and power, as our co-ordinate system is centred on the star. Hence, the torque / poweris the sum of the disc contribution, denoted by subscript ’disc’,and that of the indirect term, labelled by ’ind’. The disc torqueand the disc power are calculated as Γ disc = N r , N φ X i , j = r p × F p i j , P disc = N r , N φ X i , j = r p · F p i j , (5)where r p is the location of the planet and F p i j is the gravitationalforce between the planet and the mass in the i j th cell of the disc.Since the calculations have been performed in the star-centrednon-inertial frame, there is an additional torque / power compo-nent due to the frame acceleration a ⋆ , which is the accelerationon the star exerted by the disc. The indirect torque / power van-ishes when the disc is completely symmetric and there is onlyone planet in the system. Otherwise, this term must be taken intoaccount. The indirect torque Γ ind and indirect power P ind are cal-culated using the following equations Γ ind = M p r p × a ⋆ , P ind = M p r p · a ⋆ . (6)The left column of Fig. 19 shows that the torque and powerare zero for the single planet of 1p10 DZ model. On the right,similar quantities are presented for the 2p20 DZ model with twoplanets (left column in Fig. 5). When the planets are trapped intheir final halting positions, at t ∼ orbits, the torque andpower acting on the inner planet are positive and those of theouter planet are negative, while their sum cancels out. If we con-sider each planet individually, we would expect the inner planetto move outwards because of the positive power, and oppositelyfor the outer planet. To maintain an equilibrium, the two powersneed to be equal but have opposite sign. This results in zero netmigration for the system as a whole. While, the migration of thesystem is governed by the total power on the system , the van-ishing total torque indicates that the eccentricities in the systemare also not changing anymore. The planet-planet interaction can be neglected since these quantities,which are very oscillatory, average out when the system has reachedequilibrium. If the calculation is carried out in the centre of mass frame, the sum ofthe power on all objects including the star should vanishes to maintainthe equilibrium. a [ r ] e t [orbit] T o / T i Fig. 20.
Results of the N-body simulation with τ a i = . × , τ a o = − . × , τ e i = τ e o = − × , and M i = M o =
20 M ⊕ are over-plotted those of 2p20 DZ model. The dark / pale colours refer to the N-body / hydrodynamical simulation. To ease the comparison, we shiftedthe time axis for the hydrodynamical simulation such that t = To examine our claim in more detail, we performed addi-tional N-body simulations using analogous parameters to the2p20 DZ model. To calculate an equilibrium situation, we con-sider a system of two planets, where the inner planet is migrat-ing outwards and the outer planet inwards. Upon capture in 4:3resonance (as observed in the hydrodynamical simulation), theyshould come to a halt and reach an equilibrium. When the mi-gration of the planetary system stops, its energy and angular mo-mentum stays constant, on average. This implies that the totalpower should vanish, P i + P o =
0. Here, we neglect the interac-tion energy between the two planets as it averages out after theyreach their parking position at the inner edge of the disc. Usingthis equilibrium condition and the relation P = | E | ˙ aa ≡ | E | τ a , (7)which follows from eq. (3), we obtain the ratio of the migrationrate of the inner to the outer planet as τ a o τ a i = − a i a o M o M i = − T o T i ! − / M o M i . (8)Because our system is in a 4:3 resonance, T o / T i = /
3, and themass of outer planet is twice of the inner, we get | τ a o /τ a i | ≃ . τ a i ≃ . × .The migration of the outer planet is simply 1 .
65 times slowerthan the inner planet and in opposite direction, meaning τ a o ≃− . × . We input these timescales in an N-body code (Kleyet al. 2005) that solves the equation of motion of the planets bythe method in Lee & Peale (2002). To assure that the system set-tles into the desired resonance, we place the planets just outside Article number, page 12 of 22. Ataiee & W. Kley: Trapping planets at a disc inner edge t [orbit] P d i s c + P i n d inneroutersum Fig. 21.
Time evolution of the power on each planet and total powerfor the 2p20 IB model. The pale green line is the total power withoutincluding the indirect term. The inward migration of the whole systemcan only be explained with considering the indirect term in calculatingof the torques and powers. of the 4:3 resonance. Besides the migration timescales, eccen-tricity damping timescales τ e i and τ e o are also needed. For sim-plicity, we assume they are equal and find τ e i = τ e o = − × to give a good match between the N-body and hydrodynamicalresults. From eqs. (4) and (2) we have Γ L =
12 ˙ aa − e − e ˙ ee ≡
12 1 τ a − e − e τ e . (9)Using the obtained values for the damping timescales, the corre-sponding torques are Γ i ≃ . × − and Γ o ≃ − . × − .These values are indicated in the upper right panel of Fig. 19with dashed lines, and agree well with the averaged torque onthe planets obtained from the hydrodynamical simulation.Figure. 20 demonstrates the results of the N-body comparedto the final part of the hydrodynamical simulation. The factthat in the N-body simulation, the final resonance and locationsmatch with the hydrodynamical results implies that the properratio of migration timescales has been used, in agreement witheq. (8). For a di ff erent ratio, the planets would not come to anequilibrium but would migrate jointly, either outwards or in-wards. One needs to keep in mind that the power and torque needto be calculated by including the indirect term when working inan accelerated reference frame.From the equilibrium of the torque, Γ i + Γ o =
0, we canderive another important relation for the eccentricity dampingtimescales at the final parking position of the planet pair. Usingeq. (2) for angular momentum and energy, eq. (8) from powerequilibrium, and assuming that e ≪ (cid:18) − e τ ai τ ei (cid:19)(cid:16) − e τ ao τ eo (cid:17) = T o T i . (10)This implies that the final eccentricities in equilibrium are de-termined by the ratio of migration over eccentricity timescale τ a /τ e . One should remember that τ a i is positive and τ a o nega-tive. The eccentricity damping timescales are always negativeand much shorter that the migration timescales. Using the spec-ified timescales and the equilibrium eccentricities of the N-bodysimulations, e i = .
05 and e o = .
6. Comparing N-body to Hydrodynamicalsimulations
In this section we compare the outcome of selected hydrodynam-ical models to customized N-body simulations, using the methoddescribed in appendix A. Adopting the density and temperaturedistributions from the hydrodynamical models the N-body simu-lations calculate the migration of the planets. As in our hydrody-namical simulations, we used ǫ = . β =
0. For calculating α Σ , wetook the surface density profile from the hydrodynamical simu-lation at two specific times. The viscosity transition is also takeninto account as in eq. (1).In Fig. 23, we display the results of two sets of N-body sim-ulations using the initial conditions of 2p20 and 3p2020 DZmodels, where the mass of inner planet is 10 M ⊕ and the outerplanet(s) 20 M ⊕ . Each set contains two identical setups and di ff eronly in their surface density profile, which is chosen at t =
50 or-bit (unperturbed) and t =
13 000 orbit (evolved). We chose twosnapshots to take a possible change of the surface density inthe hydrodynamical simulations into account. The left columnin this figure contains the migration of the planets, and on theright the used surface density profiles are plotted. The thickersections in the right panels mark the region where the torque ispositive. The upper / lower rows show the results for the modelwith two / three planets. In all of these models, the inner planetsin each model is pushed by the outer planets into the positivetorque region regardless of which surface density profile is used.The only di ff erence between the outcome of these simulationsis the slower migration of the planets in the model with evolvedsurface density profile. This slower migration is because of theflattened bump that is created by the planets in the hydrodynam-ical simulations.In the two-planet simulations, the inner planet is first trappedat the zero torque location and then pushed further in by the in-coming (more massive) outer planet. The migration of the outerplanet is halted when it reaches the trap at outer edge of the vis-cosity transition. By then the inner planet has been pushed intothe negative torque region again and migrates slowly inwards.However, in the three-planet cases, the two inner planets are firsttrapped in a resonant configuration while the outer planet is ap-proaching. Later, at t ≈
50 000, the outer planet catches the innerpair in resonance and pushes them out of the positive torque re-gion into the inner disc. The migration of the inner planets canonly be halted if they are still in the positive torque region at themoment the outer planet reaches its final trapping location. In allof our N-body cases, we saw that the outermost planet was ableto push the inner planets inward beyond the trap region into tothe negative torque region, such that they continue a very slowinward migration. In strong contrast, in the corresponding hy-drodynamical simulations, the inner planet was able to stop theinward migration of the planetary system, even though an unsta-ble evolution occurred in the three planet case.Our results hint that the strong positive torques experiencedby the planets at the trap in the hydrodynamical cases are not
Article number, page 13 of 22 & Aproofs: manuscript no. discedge -0.05 0 0.05 e i cos(ϕ i ) -0.0500.05 e i s i n ( ϕ i ) (a) 3:2 -0.05 0 0.05 e i cos(ϕ i ) (b) 4:3 -0.05 0 0.05 e i cos(ϕ i ) (c) 5:4 t [orbit] T o / T i Fig. 22.
Panel (a)–(c) : Evolution of the eccentricity and resonant angle for the 2p20LH IB model when the planets are engaged in 3:2, 4:3, and5:4 resonances, respectively. Panel (d) shows the evolution of orbital period ratio. The colour-coding corresponds to the right panel, and indicatesthe time evolution of the noted resonances in panels (a)–(c). Explicitly said, the blue colour in each panel (a)–(c) indicates the time before thesystem gets into the resonance, which is noted on top of each panel and also in panel (d), and the red colour represents the time after the system isin resonance. The exact time can be extracted from comparing the colour in panel (d) with the ones in panel (a)–(c). t [orbit] a [ r ] Σ/Σ a [ r ] t = 50t = 13000 Fig. 23.
The right panels show the surface density in the horizontal axisand distance from the star on the vertical. The surface density profilesat t =
50 and t =
13 000 orbit are denoted by blue and orange. Thethick sections represent where the torque on a 10 M ⊕ circular planet ispositive. The left panels show the results of N-body simulations for theinitial conditions as 2p20 (top) and 3p2020 (bottom) DZ models. Thecolours correspond to the surface density profiles shown on the rightpanels. One can see that the push by the outer planet in each pair canshove the inner planet out of the trapping point. fully reflected by the N-body simulations. This is depicted inFig. 24, which shows the torque from the hydrodynamical andN-body simulations (using the unperturbed surface density pro-file) against the semi-major axes of the planets. The grey solidand dashed lines are torques calculated using the formulae in theN-body code for a circular and an eccentric planet. For the outerplanet (lower panel), the torque agrees excellently between thehydrodynamical and N-body models down to a ∼ .
5. The dis-crepancy increases when approaching the transition edge. Spe-cially, in the hydrodynamical model, the torque on the innerplanet (upper panel of Fig. 24) is positive and relatively largecompared to the value given by formulae. It is why it can com-pensate the push by the outer planet and consequently halts at r ∼ .
4. In the N-body simulation, the torque on the inner planetis mostly negative. Even when the torque is positive inside the −1012 Γ i nn e r formulae (e = 0)formulae (e = 0.05)N − bodyhydro a [r o ] −505 Γ o u t e r formulae (e = 0)formulae (e = 0.025)N − bodyhydro Fig. 24.
Torque from the hydrodynamical and N-body simulations ver-sus the semi-major axes of the planet for the inner (top) and outer planet(bottom) in 2p20 DZ model. The grey lines denote the torque calculatedfrom the formulae used in the N-body code as a function of distancefrom the star for a circular (solid) and an eccentric (dashed) planet. Thevalues of the used eccentricities equal those of the hydrodynamical sim-ulation after equilibrium. Comparing the grey lines and the outputs ofthe simulations show that these formulae are capable of correctly calcu-lating the torque only for the outer planet before it approaches transitionzone. transition zone, it is much smaller than that of the hydrodynam-ical simulation. This indicates that the actual hydrodynamicalvalues for torque and eccentricity damping ( τ L , τ e ) inside andclose to a transition zone di ff er from those calculated from theformulae which are estimated from relatively smooth disc sur-face density profiles. We would like to mention that the N-bodytorque would be even smaller if the evolved surface density pro-file was used in Fig. 24 because both the smaller surface densityand also flatter profile would results in smaller torques.For illustrative purposes and to connect to the next exam-ple, we show the inverse of the migration timescales for the hy- Because the migration timescale is extremely large for the trappedplanets, showing its inverse is more proper for comparison.Article number, page 14 of 22. Ataiee & W. Kley: Trapping planets at a disc inner edge a [r o ] −2−1012 / τ a [ o r b i t − ] inner hydroouter hydro inner N − bodyouter N − body Fig. 25.
Inverse of migration timescale for the hydrodynamical and N-body simulation of 2p20 DZ model. drodynamical and N-body simulation in Fig. 25 . The migrationtimescale of the inner planet in the hydrodynamical simulation ismuch smaller than its N-body counterpart. For the outer planet,the migration timescale is very similar in both simulations until a ∼ .
5. As the planet approaches the transition zone, the migra-tion timescale in the hydrodynamical model deviates from thatof the N-body.To verify our suggestion that the di ff erence between hydro-dynamical and N-body simulations is hinged to a very small pos-itive torque near and inside of the trap region in the N-body runs,we performed a set of additional N-body simulations in whichthe migration timescale is given by an analytical model. For eachplanet, we choose a negative constant reference migration time, τ a , except inside a transition zone where it becomes positive inorder to mimic the planet trap. Mathematically speaking, the mi-gration timescale is given by τ a = τ a / f trans with f trans definedas f trans = + a (cid:20) tanh (cid:18) r − r ref ∆ (cid:19) − (cid:21) + (11) + b r − (1 + c ∆ ) r ref ∆ ! exp − (cid:18) r − r ref ∆ (cid:19) ! , where r ref is the location of the transition, ∆ determines thewidth, and a , b , c control the shape of the transition zone: a thedi ff erence between inner and outer inward migration, b strengthof the trap, and c the bump strength. The analytical form of thisfunction is motivated by the results from our hydrodynamicalmodelling. Di ff erent sets of these parameters corresponds to dif-ferent shape for the disc inner edge. For example, the third termadjusts the torque if a surface density bump exists before thezero-torque location. Attention should be paid that finding a re-verse process which gives the exact shape of the disc inner edgeis not trivial because of the torque calculation complexities. Infact, it is very hard, if not impossible, to find the shape and thedisc properties only by having the torque on the planet. This function is defined such that it gives a faster migra-tion in the outer than the inner disc plus a region where the mi-gration becomes very slow and outwards. Far from r ref , where f trans =
1, the migration of the planet is simply inwards with thetimescale of τ a . As the planet approaches r ref , f trans becomes ini-tially larger than unity. This is where the migration of the planetbecomes faster due to the existence of a surface density bump be-fore the disc edge. Afterwards, f trans decreases until it becomes Equation 11 is a general profile that can be used for any type of plan-etary trap only by proper adjustment of the paramters, either the trap isa dead-zone or disc inner edge. However, we note that in the IB mode,the paramters would be time-dependent due to the change of the surfacedensity profile by time.
Parameter DB D SB S ∆ r ref a b
20 8 3 1.4 c Table 2.
The list of parameters for f trans profiles shown in the top row ofFig. 26. The letters refer to, D: deep trap, S: shallow trap, B: bump. −10−505 f t r a n s DBD SBS a [r ] t [ o r b i t ] a [r ] Fig. 26.
The results of four N-body simulations in which a migrationtimescale transition zone, as depicted with the same colour on the toppanels, is used. The vertical dashed line corresponds to f trans =
0, andmarks the transition from positive to negative migration. zero and brings the migration of the planet to a halt. If by somemeans, for example the push of another planet, the planet movesto where r < r ref , a negative f trans causes an outward migrationsimilar to a disc inner edge. If the planet is pushed even moreinwards, such than the planet is out of the disc inner edge, weexpect the migration of the planet becomes again inwards butvery slow due to the low surface density. This is what f trans doesfor r ≪ r ref . We would like to insist that this factor does not pro-duce the torque from the asymmetric features such as vortices,which produce usually a stochastic torque.We examined two sets of models, one with a deep and onewith a shallow drop in f trans , resembling the hydrodynamical andN-body simulations, respectively. In each set, one has a bumpoutside of the drop and one does not. The latter models wereperformed in order to check the role of the bump, which speedsup the inward migration just outside of the transition zone, asseen in the hydrodynamical models. The upper panels of Fig. 26show f trans for the deep (left) and shallow (right) drop models.The used parameters in these profiles are listed in table 2.We used the same masses and initial locations as 2p20 andconsidered τ a o = − and τ a i = − × . The eccentricitytimescales are fixed to τ e i = τ e o = − × . These numbers aresimilar to those of Sec. 5. The result of the simulations are shownin the bottom panels of Fig. 26 with the corresponding coloursto each of the transition profiles.In all of the models, the planets get into 3:2 resonance afterabout 2800 orbits. In the deep models (left panel), the migrationof the planets are halted where the total power acting on bothplanets vanishes. This location is near f trans = ff erence be-tween the models with and without bump. In the shallow model Article number, page 15 of 22 & Aproofs: manuscript no. discedge t [orbit] ( τ e / τ a ) h Fig. 27.
Ratio of migration and eccentricity damping timescales scaledby the aspect ratio squared versus time. with a bump, the inner planet continues its migration to the tran-sition zone and is not able to stop the migration of the outerplanet. Finally, the inner planet leaves the transition region whilethe outer planet stays at the trapping location and therefore, theresonance breaks. However, in the model without bump, the mi-gration of the inner planet halts when the outer planet reachesthe trapping point and they remain in 3:2 resonance. This dif-ferent behaviour is due to the fact that, in spite of the S model,the inner planet has passed the transition zone before the outerplanet approaches the trapping point in the SB model.To summarise, based on the results in this section, we foundthat the N-body simulations favour resonant pushing, unlike hy-drodynamical simulations. This can be due to the fact that theformulae, which are used in the N-body calculations, do not giveproper values where the surface density deviates from the simplepower-law profiles. This point requires more detailed analyses ina separate study.
7. Overstability
Among the results of our simulations in Sec. 3, the IB modelswith smaller aspect ratio, h = .
03, show overstable behaviour(see runs 2p20LH, 3p2020LH, and OnebyoneLH, as displayedin Figs. 8, 12 and 16). Overstability for planets at the inner edgeof a protoplanetary disc was observed in N-body simulations byXiang-Gruess & Papaloizou (2015), but has not been reportedyet convincingly in hydrodynamical simulations .In an overstable system, the planets enter successively var-ious first-order resonances that are broken after a short time.In the transient resonances, the eccentricities oscillate aroundan equilibrium value with increasing amplitudes until the res-onance breaks and the planets enter an even tighter resonance.Overstability occurs when the eccentricity damping is not e ffi -cient enough. As a result, the eccentricities of the planets, thatare excited due to their mutual gravitational interaction, gradu-ally grow while the libration amplitude increases (Goldreich &Schlichting 2014). This growth continues until the resonance isbroken eventually.Figure 22 shows the overstable behaviour in more detailsfor the 2p20LH IB model, in which the system becomes over-stable during the 3:2 and 4:3 resonances (see also Fig. 8). Pan-els (a), (b), and (c) demonstrate how the eccentricity e and theresonant angle φ of the inner planet are evolving during the 3:2, Hands & Alexander (2018) attributed the mean-motion resonancebreaking that happened in their hydrodynamical simulations to oversta-bility. However, the eccentricity evolution of the planets in their modelsdo not resemble such a behaviour. ,
0) represents the eccentricity, and itsangle with the horizontal axis is a measure of resonant angle.The colour scheme indicates the time evolution as it is shown inpanel (d), where the time evolution of the orbital period ratio isplotted. The resonant angle for a first-order j + j resonanceis defined as φ i = ( j + λ o − j λ i − ̟ i , where λ and ̟ denotethe mean longitude and longitude of the pericentre, respectively.If φ i , which measures the angle between the conjunction and thepericentre of the inner planet, librates around a constant value,the system is in a resonance. Panel (a) implies that as the planetssettle into the 3:2 resonance, the eccentricity of the inner planetgrows to ∼ .
04 and the resonant angle librates around 3°. Theamplitude of the libration, as well as the eccentricity, increasesuntil it circulates and the system leaves the resonance. The sit-uation for the 4:3 resonance in panel (b) is very similar, exceptthat the equilibrium eccentricity and resonant angle are smaller.In contrast, the final resonance 5:4 seems to be stable around e ∼ .
025 and φ i = τ e /τ a issmall enough. They also present a condition between the equilib-rium eccentricity e eq and mass of the planets. They suggest thatif e eq . (M p / M ∗ ) / , the system remains stable. For the innerplanet in our 2p20LH IB model, (M p / M ∗ ) / ∼ .
03. There-fore, their stability criterion is only satisfied in our case for the5:4 resonance. They also found that e eq ∼ ( τ e /τ a ) / ∼ h thatpredicts e eq ∼ .
03 which roughly agrees with our results. How-ever, the study by Xu et al. (2018), that uses a migration modelbased on planet-disc interaction studies, suggests that the equi-librium eccentricity of the resonant planets are larger than thosegiven by the studies that used constant migration and eccentricitydamping rates. Our results seem to match better with the workof Goldreich & Schlichting (2014).It has been shown that τ e /τ a is of the order of h (e.g. Tanaka& Ward 2004). Therefore, we expect that overstability is moreprobable to occur in thicker discs while we observe it for ourthinner models. This may be explained by noting that the planetsin the models with h = .
03 open a partial gap around their or-bits. Although this partial emptiness around the orbit of the plan-ets slow down their migration, it also decreases the e ffi ciency ofthe eccentricity damping. Therefore, it seems that in our simu-lations the net e ff ect of the gap opening is more in the favour oflarger τ e /τ a . Figure 27 shows τ e /τ a scaled by h for the model2p20LH and 2p20 for when the planets are in a resonance. Wecan see that, in 2020LH model, τ e /τ a is larger at the overstableresonances than that of the stable resonance. Interestingly, thescaled ratio has similar values for the stable resonance in both2p20 and 2p20LH models. However, this speculation needs tobe justified by more rigorous studies.
8. Stability of the systems and comparing withobservation
Formation of planetary systems is a long journey from when theyform in a protoplanetary disc until the disc dispersal. Thus, it isuncertain if those of our systems that survive from instabilityduring the simulation time could remain stable afterwards. Be-cause continuing all of our simulations would be computation-ally very costly, we only examined whether DZ and IB 2p20 and2p20LH models retain their configuration during the disc dis-
Article number, page 16 of 22. Ataiee & W. Kley: Trapping planets at a disc inner edge a [ r ] −3 −1 Σ / Σ t [10 orbit] a [ r ] r [r ] −3 −1 Σ / Σ Fig. 28.
Time evolution of the models 2p20 after the photo-evaporationis activated. The left panels show the semi-major axes of the planets,and the right ones demonstrating the surface density evolution of thedisc. The vertical lines in the left panels correspond the time when thesurface densities are plotted. persal. We include the X-ray photo-evaporation model of Owenet al. (2012) with the X-ray luminosity of L x = erg / s andcontinue the mentioned simulations for over thousands of orbitsuntil the disc is dispersed around the planets. The outcome ofthese tests show that none of them become unstable (Fig. 28 andB.7). Similar outocmes can be seen in the results of the N-bodysimulations which include photo-evaporation (e.g. Migaszewski2015, 2016).To inspect the stability of other models, we used the criteriapresented by Chambers et al. (1996) and Gladman (1993), whichsuggest that a low-mass planetary system would remain stable if ∆ / R H > √ ∆ / R H <
10 for more than two planets. Here, ∆ is the initialdi ff erence of the semi-major axes of the planets, and R H is theirmutual Hill radius. For all of our two-planet models, this condi-tion is well satisfied. Among the multi-planet models, it is onlyIB 3p1020 that has ∆ / R H above 10. Therefore, it is the only two-planet system that would survive after a long-term evolution.We also compared these stable systems with observed exo-planets and found two systems, Kepler-804 and K2-189, thathave configurations close to our results. Kepler-804 has twoplanets revolving around a Solar-type star (Morton et al. 2016)with the size of the outer planet about twice of the inner one (thatis almost an Earth-size planet) and the outer to inner period ratioof ∼ . / . ≃ .
5. Using the mass-radius relation in Wolfganget al. (2016) for sub-Neptune-sized planets, the mass and periodratios of this system resemble our DZ 2p20 model (although themasses are smaller than ours). The other system, K2-189, hasalso two planets with the radii of about 2.48 and 1.51 times of theEarth that are revolving around a star with the mass of 0 . M ⊙ Mayo et al. (2018). The period ratio of the planets is about 1.3that is near 4:3 commensurability. Using the relation in Wolf-gang et al. (2016) gives a mass ratio of the outer to inner ∼ . This high value of L x is chosen to speed up the simulations.
9. Summary and conclusion
Using hydrodynamical simulations, we investigated whether atrapped planet at a dead-zone inner edge (DZ) or a disc innerboundary (IB) can be pushed into the inner cavity towards thecentral star by a resonant chain. We examined systems with twoand three planets where the planets were placed in the disc simul-taneously, and systems with more than three planets in which theplanets were added one by one after the previous planets formeda resonant chain. The DZ model has a very narrow and steeptransition in the surface density, and the IB model has a muchshallower profile. The mass of the inner planet was always 10 M ⊕ while the outer planets had either the same mass or were moremassive. We summarize and discuss our findings in the follow-ing.The steepness of the inner disc edge, i.e. the gradient of thesurface density plays a key role in the dynamical evolution ofour systems. The DZ model, where the sign of the power (andtorque) abruptly changes at the edge of the transition zone, inwhich the planet experiences positive power, acts as a dam thatcan stop an incoming resonant chain of planets. The inner planetcannot pass this dam unless an instability in the system hurls itto the inner disc. On the other hand, in the IB model, which hasa smoother transition and a less steep surface density profile, aresonant chain can push the inner planet to the inner cavity. How-ever, even in this case it requires multiple, and / or more massiveplanets to push the inner planet over the trap. Test calculationsusing a broader transition zone in the DZ setups produce widerresonant systems which are less prone to instability These testsindicate that if the viscosity transition zone is very wide, the in-ner planet might be pushed out of its trapping point, similar tothe IB models.In the two-planet cases, the final configurations in the DZmodels are usually tighter than in their IB counterparts. A typicaloutcome for a 10 and 20 M ⊕ combination is the 4:3 resonance.Most of the three-planet systems in the DZ models became un-stable and a scattering event between the two inner planets oc-curred, while in the IB models, the planets were trapped in a3:2:1 resonant chain. In the models where we added the planetsone by one, increasing the number of planets has no e ff ect on theinnermost planet and only rises the chance of instability.Planets in the models with h = .
03 behave di ff erentlythan those with h = .
05. In the cooler discs, the planets openpartial gaps around their orbits. The gap opening changes boththe migration rate and eccentricity damping such that overstableoscillations of the planets’ eccentricities occur (Goldreich &Schlichting 2014) and the planets reach successively tighter res-onances, jumping for example from 3:2 to 4:3 to 5:4. Our resultsare the first hydrodynamical simulations that clearly show theoverstable behaviour in resonant planetary systems, see Fig. 22.For computational reasons we focussed our study on transitionslocated at 1 au with typical aspect ratio of h = .
05. However, theinner edge of discs lies about 10 times closer to the star wherethe disc is much thinner h = .
02 to 0 .
03, and we expect thatunder these conditions overstability may occur frequently. Thisraises and highlights the need for more detailed studies in thistopic.In addition to the aspect ratio, we varied the disc mass (sur-face density, Σ ) and viscosity ( α ). Using a higher surface densityfor the disc, only speeds up the evolution while having no spe-cific e ff ect on the final results. This is in agreement with earlystudies of resonant capture by Lee & Peale (2002) who showedthat it is the ratio of migration over eccentricity damping that de-termines the outcome of capture, and not the absolute timescale. Article number, page 17 of 22 & Aproofs: manuscript no. discedge
Although the planets in the less viscous disc were able to slightlymodify the surface density profile of the disc through partial gapopening, the outcome of the models are not significantly di ff er-ent than those with higher viscosity.In the DZ simulations the steep positive surface density pro-file gives rise to the creation of vortices, that can potentially in-terfere in the migration of the planets. Our results show that thepresence or number of vortices at the edge of a dead-zone doesnot influence the evolution of the system greatly, besides a smallenhancement of the susceptibility to unstable evolutions.The eccentricity of the planets remain small as long as no in-stability happens or they are not pushed into the inner disc, wherethe surface density is very low. This indicates that the eccentric-ity damping timescales in our models are still much shorter thanthe migration timescales. The comparison of our hydrodynami-cal and N-body simulations show indeed a ratio of | τ a /τ e | ≈ ⊕ planet.We analysed in detail the conditions for an equilibrium park-ing of a system of two planets at the inner edge. In agree-ment with previous results for migrating planets (Cresswell et al.2007), we show that it is determined by vanishing total power,where all objects have to be taken into account. This could beformulated as a condition for the ratio of migration timescales,see eq. (8). The additional requirement of vanishing total torqueyields a condition for the equilibrium eccentricities, see eq. (10).This last equation implies particularly that resonant planetarysystems created by trapping must have a non-zero eccentricity.The validity of these relations was confirmed in our hydrody-namical as well as N-body simulations.As hydrodynamical simulations are very time-consuming,we performed additional customised N-body simulations in or-der to test the often-used approach of taking analytical formulaefor the planet’s torque and eccentricity damping. In the appendix,we demonstrated that the direct usage of the disc parameters,namely density and temperature profiles, from the hydrodynam-ical simulations as inputs for the forces in the N-body formulaedoes result in di ff erent evolutions of the planetary systems. Inparticular, the inner planet in the N-body simulations was notable to stop the resonant chain. This indicates that the positivetorque (power) contribution from the inner edge of the disc isunderestimated by the analytical formulae.The main conclusion of our study is that it is di ffi cult forthe planets to be pushed into the cavity created by a dead-zoneinner edge through resonant migration due to its steep densityslope. However, for a smoother planetary trap, such as a discinner boundary a resonant chain can help pushing the planetsinwards. Therefore, if we consider a dead-zone inner edge thatdoes not coincide the the inner boundary of the disc but lies far-ther out, the formed planets inside the dead-zone would havedi ffi culties to get into the inner disc, unless a the resonant chainbecomes unstable. Nonetheless, if they manage to pass the damof the viscosity transition, they can easily move towards the discinner edge by the help of resonant migration.In order to compare our results with the observed exo-planetary system, we also examined the stability of our systemsand found that the two-planetary systems are more probable toremain stable. Among our stable models, we only found two ob-served systems, Kepler-804 and K2-189, that are similar to theoutcome of our 2p20 models. Although, we can not make a con-clusion based on these two observed systems, we suggest thatthe more packed systems are more likely to form in a disc witha steeper inner edge.We would like to note that our results are obtained using lo-cally isothermal disc. The study by Faure & Nelson (2016), who investigated the planet trapping at a dead-zone edge in a non-isothermal disc using 2D and 3D magnetohydrodynamical sim-ulations, hints that the inner edge of a dead-zone can be a mass-dependent barrier for the migration of vortices, which formedinside the dead-zone, and also planets.As mentioned before, for computational reasons, we placedthe inner boundary at a distance of 1 au using typical values forthe disc surface density, aspect ratio, and viscosity at that lo-cation. The models with smaller aspect ratio or higher surfacedensity may resemble more the conditions closer to the star. Inparticular, the thinner disc models ( h = .
03) give an insight to-wards how the resonant pushing can be di ff erent in the inner discin the later stages of the disc evolution. Hence, we might expectthat overstable librations with subsequent tightening of the res-onant chains are a frequent outcome and may help in producingvery compact resonant systems such as Kepler-223 (Mills et al.2016).This study has been carried out using a background disc pro-file which dictates divergent migration in the absence of a plane-tary trap. If one used a di ff erent disc setup that governs the con-verging migration, instability would be more probable.One of the questions that remain open in this study is thatwhat happens to the resonant planets once they are pushed intothe inner cavity. In some IB models, where the inner planet ispushed all the way to the region with very low surface density, itseccentricity increased to large values. Whether the system willremain stable afterwards is a question we could not answer usingour hydrodynamical simulations.The tidal interaction between the star and the planets can bean important phenomenon for the planets located very close tothe star. In this work, as the first hydrodynamical study whichtries to simulate multi-planets at the disc inner edge, we ig-nored this e ff ect for simplicity. How much the evolution can bechanged by this e ff ect is a question for future works. Acknowledgements.
We acknowledge the support of the DFG priority programSPP 1992 "Exploring the Diversity of Extrasolar Planets under grant KL 650 / / his constructivecomments. References
Ataiee, S., Dullemond, C. P., Kley, W., Regály, Z., & Meheut, H. 2014, A&A,572, A61Ataiee, S. & Kley, W. 2020, arXiv e-prints, arXiv:2002.08077Baruteau, C., Crida, A., Paardekooper, S. J., et al. 2014, in Protostars and PlanetsVI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 667Benítez-Llambay, P., Masset, F., & Beaugé, C. 2011, A&A, 528, A2Bitsch, B. & Kley, W. 2010, A&A, 523, A30Brasser, R., Matsumura, S., Muto, T., & Ida, S. 2018, ApJ, 864, L8Carrera, D., Ford, E. B., & Izidoro, A. 2019, MNRAS, 486, 3874Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996, Icarus, 119, 261Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580Coleman, G. A. L. & Nelson, R. P. 2014, MNRAS, 445, 479Coleman, G. A. L. & Nelson, R. P. 2016, MNRAS, 457, 2480Cresswell, P., Dirksen, G., Kley, W., & Nelson, R. P. 2007, A&A, 473, 329Cresswell, P. & Nelson, R. P. 2006, A&A, 450, 833Cresswell, P. & Nelson, R. P. 2008, A&A, 482, 677Cui, Z., Papaloizou, J. C. B., & Szuszkiewicz, E. 2019, ApJ, 872, 72de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529Dressing, C. D. & Charbonneau, D. 2015, ApJ, 807, 45Faure, J. & Nelson, R. P. 2016, A&A, 586, A105Fendyke, S. M. & Nelson, R. P. 2014, MNRAS, 437, 96Flock, M., Fromang, S., Turner, N. J., & Benisty, M. 2017, ApJ, 835, 230Gladman, B. 1993, Icarus, 106, 247Goldreich, P. & Schlichting, H. E. 2014, AJ, 147, 32Hands, T. O. & Alexander, R. D. 2018, MNRAS, 474, 3998Ida, S., Muto, T., Matsumura, S., & Brasser, R. 2020, MNRAS, 494, 5666
Article number, page 18 of 22. Ataiee & W. Kley: Trapping planets at a disc inner edge
Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750Kley, W., Lee, M. H., Murray, N., & Peale, S. J. 2005, A&A, 437, 727Kley, W. & Nelson, R. P. 2012, ARA&A, 50, 211Lee, E. J. & Chiang, E. 2017, ApJ, 842, 40Lee, M. H. & Peale, S. J. 2002, ApJ, 567, 596Lin, D. N. C., Bodenheimer, P., & Richardson, D. C. 1996, Nature, 380, 606Long, M., Romanova, M. M., & Lovelace, R. V. E. 2005, ApJ, 634, 1214Masset, F. S. 2001, ApJ, 558, 453Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 642, 478Mayo, A. W., Vanderburg, A., Latham, D. W., et al. 2018, AJ, 155, 136Migaszewski, C. 2015, MNRAS, 453, 1632Migaszewski, C. 2016, MNRAS, 458, 2051Mills, S. M., Fabrycky, D. C., Migaszewski, C., et al. 2016, Nature, 533, 509Miranda, R. & Lai, D. 2018, MNRAS, 473, 5267Morton, T. D., Bryson, S. T., Coughlin, J. L., et al. 2016, ApJ, 822, 86Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123Ogihara, M., Duncan, M. J., & Ida, S. 2010, ApJ, 721, 1184Ogihara, M. & Ida, S. 2009, ApJ, 699, 824Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880Paardekooper, S. J., Baruteau, C., Crida, A., & Kley, W. 2010a, MNRAS, 401,1950Paardekooper, S. J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293Paardekooper, S.-J., Lesur, G., & Papaloizou, J. C. B. 2010b, ApJ, 725, 146Papaloizou, J. C. B. & Larwood, J. D. 2000, MNRAS, 315, 823Papaloizou, J. C. B. & Szuszkiewicz, E. 2005, MNRAS, 363, 153Raymond, S. N., Boulet, T., Izidoro, A., Esteves, L., & Bitsch, B. 2018, MNRAS,479, L81Rein, H. & Liu, S. F. 2012, A&A, 537, A128Rein, H. & Tamayo, D. 2015, MNRAS, 452, 376Romanova, M. M., Lii, P. S., Koldoba, A. V., et al. 2019, MNRAS, 485, 2666Romanova, M. M. & Lovelace, R. V. E. 2006, ApJ, 645, L73Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2002,ApJ, 578, 420Tanaka, H. & Ward, W. R. 2004, ApJ, 602, 388Winn, J. N. & Fabrycky, D. C. 2015, ARA&A, 53, 409Wolfgang, A., Rogers, L. A., & Ford, E. B. 2016, ApJ, 825, 19Xiang-Gruess, M. & Papaloizou, J. C. B. 2015, MNRAS, 449, 3043Xu, W., Lai, D., & Morbidelli, A. 2018, MNRAS, 481, 1538Zanni, C. & Ferreira, J. 2009, A&A, 508, 1117
Article number, page 19 of 22 & Aproofs: manuscript no. discedge
Appendix A: N-body simulations - methods
N-body simulations of migrating planets are a common alterna-tive for the lengthy and computationally expensive hydrodynam-ical simulations. In studies where the long term evolutions ofmany planetary systems are desired, using hydrodynamical sim-ulation is not only expensive but can even be unfeasible. There-fore, using N-body simulations augmented with a suitable mi-gration treatment is a good option (e.g. Papaloizou & Larwood2000; Ogihara & Ida 2009; Coleman & Nelson 2014; Izidoroet al. 2017). The basic ingredients for performing such a simula-tion is an N-body code, an expression for the disc torque actingon the planet, a relation for calculating the eccentricity dampingtimescale, and an equation of motion that modifies the accelera-tion of the planets accordingly.In order to examine how well the results of such simulationsagree with their hydrodynamical counterparts, we employed the
REBOUND code (Rein & Liu 2012) with the WHFAST integra-tor (Rein & Tamayo 2015) with a timestep of dt = − orbits,and implemented the formulae of Paardekooper et al. (2010a,2011) into it for calculating the torques on the planets. Thetorque has two components: Lindblad, Γ L , and corotation, Γ C .Both of these are altered if the planet is eccentric. The cor-responding correction factors for the Lindblad and corotationtorques, ∆ L and ∆ C , respectively, are taken from Cresswell &Nelson (2008) and Fendyke & Nelson (2014). The formulae weuse are identical to Coleman & Nelson (2014) and Izidoro et al.(2017) except that: (a) our simulations are 2D and therefore theinclination is set to zero, (b) we applied the additional correctionfactors for the smoothing length of the planet as in Paardekooperet al. (2010a) for consistency purpose, and (c) because the discis locally isothermal, we consider the thermal di ff usivity χ as in-finite and set the adiabatic index to γ =
1. For completeness andreference, we summarize the relevant formulae briefly. The totaltorque is given by
Γ = ∆ L Γ L + ∆ C Γ C , (A.1)where Γ C is the sum of the two components, barotropic andentropy-related co-rotation torque, and each of them has a lin-ear torque and a horseshoe drag (see Paardekooper et al. 2011).Therefore, the co-rotation torque reads Γ C = Γ c , baro + Γ c , ent , (A.2) Γ c , baro = Γ c , hs , baro F ( p ν ) G ( p ν ) + Γ c , lin , baro (1 − K ( p ν )) , (A.3) Γ c , ent = Γ c , lin , ent p (1 − K ( p ν )) . (A.4)The functions F ( p ν ), G ( p ν ), and K ( p ν ) regulate the saturationof the co-rotation torque due to the disc viscosity, where p ν isdefined as 2 / p a Ω x / πν , with x s being the half-horseshoewidth of the planet. For brevity, we avoid bringing the defini-tion of the functions and refer the readers to Paardekooper et al.(2011) or Izidoro et al. (2017). Please note that the entropy-related horseshoe drag vanishes due to the locally isothermalassumption. The Lindblad torque and the components of the co- The N-body code
REBOUND is freely available athttp: // github.com / hannorein / rebound. rotation torque are as following Γ L Γ = ( − . − . β + . α Σ ) . ǫ ! . , (A.5) Γ c , hs , baro Γ = . . − α Σ ) . ǫ ! , (A.6) Γ c , lin , baro Γ = . . − α Σ ) . ǫ ! . , (A.7) Γ c , lin , ent Γ = . β . ǫ ! . − . β . ǫ ! . , (A.8)where Γ is the torque normalization at the planet’s semi-majoraxis, a , which is given by Γ = (cid:18) qh (cid:19) Σ ( a ) a Ω . (A.9)The radial slopes of the density and temperature stratificationsare denoted by α Σ = − ∂ log Σ ∂ log r and β = − ∂ log T ∂ log r , (A.10)with T being the temperature. The eccentricity correction factorsfor the torque are calculated using ∆ L = − ( e . h ) + ( e . h ) . + ( e . h ) , (A.11) ∆ C = exp (cid:18) e . h + . (cid:19) . (A.12)The calculated torque along with the eccentricity and semi-majoraxis of the planets give the timescale of angular momentumchange τ L = L / Γ . (A.13)The second ingredient, the eccentricity damping timescale τ e , isobtained from eq. (11) in Cresswell & Nelson (2008) and reads τ e = . " − . (cid:18) eh (cid:19) + . (cid:18) eh (cid:19) h a ΩΓ . (A.14)And finally, we the use accelerations given by Papaloizou & Lar-wood (2000) as a = − u τ L − r · u ) r τ e r . (A.15)The connection of the actual migration timescale, τ a , and theangular momentum change, τ L , is given by1 τ a = τ L + e − e τ e , (A.16)see eq. (9). Only for a circular orbit, the migration timescaleequals half of the angular momentum loss timescale. Appendix B: Supplementary figures
In order to keep the main text more concise we moved supple-mentary figures to this appendix. Refer to the main text and Ta-ble 1 for details on the specific model setup.
Article number, page 20 of 22. Ataiee & W. Kley: Trapping planets at a disc inner edge a [ r ] e t [orbit] T o / T i t [orbit] M i = 10M ⊕ , M o = 333M ⊕ Fig. B.1.
Similar to Fig. 4 but for 2pJup models. t [orbit] T o / T i M i = 10M ⊕ , M o = 20M ⊕ , Σ = 4 × 10 −4 Fig. B.2.
Orbital period ratio for the 2p20HS DZ model, whose surfacedensity if twice that of 2p20. a [ r ] e t [orbit] T o / T i t [orbit] M i = 10M ⊕ , M o = 20, 20M ⊕ Fig. B.3.
Similar to Fig. 10 but for 3p2020 model. e t [orbit] T o / T i t [orbit] M i = 10M ⊕ , M o = 20, 100M ⊕ Fig. B.4.
Eccentricity and orbital period ratio of for 3p20100 models. e t [orbit] T o / T i t [orbit] M i = 10M ⊕ , M o = 20, 20M ⊕ , α = 5 × 10 −5 Fig. B.5.
Evolution of eccentricity and orbital period ratio for3p2020LA which has a smaller viscosity. r [r ] −2 −1 Σ / Σ λ = 0.06λ = 0.3λ = 1.0 T o / T i t [orbit] a i n [ r ] Fig. B.6.