Interactions Among Non-Interacting Particles in Planet Formation Simulations
DDraft version July 24, 2020
Typeset using L A TEX twocolumn style in AASTeX62
Interactions Among Non-Interacting Particles in Planet Formation Simulations
Shirui Peng and Konstantin Batygin Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA
ABSTRACTOver the course of the recent decades, N -body simulations have become a standard tool for quanti-fying the gravitational perturbations that ensue in planet-forming disks. Within the context of suchsimulations, massive non-central bodies are routinely classified into “big” and “small” particles, wherebig objects interact with all other objects self-consistently, while small bodies interact with big bodiesbut not with each other. Importantly, this grouping translates to an approximation scheme wherethe orbital evolution of small bodies is dictated entirely by the dynamics of the big bodies, yieldingconsiderable computational advantages with little added cost in terms of astrophysical accuracy. Herewe point out, however, that this scheme can also yield spurious dynamical behaviour, where even in ab-sence of big bodies within a simulation, indirect coupling among small bodies can lead to excitation ofthe constituent “non-interacting” orbits. We demonstrate this self-stirring by carrying out a sequenceof numerical experiments, and confirm that this effect is largely independent of the time-step or theemployed integration algorithm. Furthermore, adopting the growth of angular momentum deficit as aproxy for dynamical excitation, we explore its dependence on time, the cumulative mass of the system,as well as the total number of particles present in the simulation. Finally, we examine the degree ofsuch indirect excitation within the context of conventional terrestrial planet formation calculations,and conclude that although some level of caution may be warranted, this effect plays a negligible rolein driving the simulated dynamical evolution. Keywords:
Planet formation (1241), Solar system formation (1530), N -body simulations (1083) INTRODUCTIONThe past three decades have seen staggering advancesin computation, and few sub-fields of astrophysics havebenefited from these developments as much as the studyof planet formation. Having remained elusive for cen-turies, the chaotic evolution inherent to the coalescenceof planetary building blocks into bona fide planets isnow within reach of modern GHz-grade machines (Dun-can et al. 1998; Chambers 1999; Rein & Liu 2012). Itis with the detailed numerical modeling of this processthat we will concern ourselves in this letter. In partic-ular, here we point out that the conventional approachto modeling quasi-Keplerian, large- N -body systems issusceptible to spurious excitation of the constituent or-bits, although we find that this effect is negligibly smallwithin real astrophysical applications. Let us begin bybriefly outlining the context of our calculations.Crudely speaking, the process of planet formationcan be sub-divided into two temporal scales: the disk-bearing phase (during which the central star is encircledby an extensive disk of gas and dust) and the post-nebular epoch (Morbidelli et al. 2012; Lissauer 1993; Armitage 2020) which takes place after the large-scaledepletion of H and He from the system. In terms of gov-erning physics, the former is subject to a multitude ofcomplex gravito-hydrodynamic processes, while the lat-ter is governed primarily by purely gravitational dynam-ics (see for example, Raymond & Morbidelli 2020, for arecent review). For definiteness, here we will restrictourselves to consideration of the latter, more qualita-tively simple mode of planet formation (correspondingto the post-nebular epoch), where the buildup of plan-etary bodies proceeds primarily via pair-wise collisionsamong planetesimals.Modeling of post-nebular evolution of planetary sys-tems is typically carried out by breaking up the cal-culation into three types of constituents: big objects,small bodies, and test particles. Big objects interactwith all other bodies in a self-consistent N -body fash-ion. Small bodies (sometimes also called semi-activeparticles) interact with big bodies but not with eachother. Finally, test particles merely track the dynamicsfacilitated by the gravitational field of the big bodies,exerting no back-reaction. The reason for this catego-rization is two-fold. First, without this treatment, the a r X i v : . [ a s t r o - ph . E P ] J u l computational burden of a typical N -body simulationwould scale as O ( N ), with N being the total numberof bodies in the problem. The big-small categorization,however, alters the computational cost of the simulationto O ( N b ) + O ( N b N s ), with N b and N s being the num-ber of big and small bodies, respectively. Because in realplanet-formation systems N b (cid:28) N s ∼ O ( N ), this treat-ment translates to a drastic reduction of computationalcosts.Second, ignoring self-gravitational stirring amongsmall bodies circumvents unphysical excitation of theorbits. This is because in an effort to keep N s to areasonably low number, the population of solid debris isroutinely modeled as a swarm of super-particles objectsthat trace the dynamics of planetesimals but containmuch more mass than the individual bodies they rep-resent. In turn, suppression of self-interactions withinthe planetesimal swarm prevents the Safronov numberΘ = ( v esc / (cid:104) v (cid:105) ) which regulates the efficiency of ac-cretion (see e.g., Safronov 1972; Lissauer 1993) fromdecreasing artificially. In other words, the big-smallparticle characterization mimics the un-modeled effectof dynamical friction.Owing to these advantages, the classification of non-central bodies into big and small particles is widely uti-lized, with important examples set within the solar sys-tem itself. In particular, over the last decade, conglom-eration of Mercury, Venus, Earth, and Mars has beenmodeled by various groups as the gravitational evolu-tion of a ∼ M ⊕ annulus of debris where the initialmass-fraction of big planetary embryos to small plan-etesimals is taken to be approximately unity (Jacobson& Morbidelli 2014; Walsh et al. 2011; see also Hansen2009). Within the context of the outer solar system,a transient (Nice-model) instability is believed to havebeen sparked early in the solar systems life time by inter-actions among big planets and a ∼ M ⊕ disk of smallplanetesimals (Tsiganis et al. 2005; Nesvorn´y & Mor-bidelli 2012). Similarly, recent simulations of the for-mation of Galilean moons (Batygin & Morbidelli 2020)treat the satellite seeds as big objects, while modelingthe much more numerous aggregate of satellitesimals assmall bodies.In this paper, we show that even if direct interactionsare turned off, some degree of self-stirring within thedisk may be unavoidable. More specifically, we carry outtests with different disk to star mass ratios and varyingnumbers of non-interacting planetesimals. The resultsgenerally show a growing trend of angular momentumdeficit, indicating a gradual increase of average eccen-tricity in time. Numerical tests using finer time steps orintegrators with higher accuracy give essentially identi- cal results. Nevertheless, our simulations also show thatthe disk of debris responsible for the generation of ter-restrial planets is not sufficiently massive for this effectto meterialize in any appreciable manner.The remainder of this letter is organized as follows:Section 2 provides a description of the setup of our nu-merical experiments. Section 3 illustrates the key resultsof our simulations. In Section 4, the formation of the ter-restrial planets is examined as an illustrative example.Our findings are summarized in section 5. NUMERICAL EXPERIMENTSSwarms of planet-forming debris are routinely envi-sioned to emerge from their natal protoplanetary neb-ulae, possessing negligible eccentricities and inclina-tions. Evolving under self-gravitation over timescalesmuch longer than an orbit, massive objects perturbone another, causing the effective velocity dispersionof the system to increase. However, this process ismarkedly not uniform, as dynamical friction is exertedupon the more massive objects by less massive bodies,causing the largest members of the planet-forming ag-gregate to circularize at the expense of further excita-tion of their smaller, but more numerous counterparts(Safronov 1972; Lissauer 1993).In the language of standard N -body simulations, thispicture can be summarized in a straight-forward man-ner: in an initially quiescent disk, big bodies heat eachother as well as the small bodies, while small bodies onlycool the big bodies. Correspondingly, if no big bodiesare present in the simulation, the perfectly circular andcoplanar architecture of the system should be preserved.Let us check this assertion with the aid of a state-of-the-art N -body code REBOUND (Rein & Liu 2012).The basic setup of our numerical experiments drawsupon standard planet formation simulations. In ourfiducial model (M3 in Table 1), we represent the diskof planetesimals in orbit of a single central body of mass M (cid:63) as N s = 1000 super-particles, which cumulativelycomprise a disk with mass M disk /M (cid:63) = 10 − , confinedto a radius range between 0 . M (cid:63) = 1, G = 1, and the time andspace units are chosen such that an orbit with 1 semi-major axis unit has a period of T = 2 π time units, whichwe define as a single “year”. The radius of each body isset to zero to suppress any collisions. The semi-majoraxes are spread within the disk uniformly from 0 . .
0, and all orbits are assumed to be initially circularwith zero inclination. The default time step is taken tobe 1 / (21+ ε ) of the orbital period with a semi-major axisof 0 . ε = 10 − is an arbitrary smallquantity. The baseline N -body integrator is the hyr- Table 1.
Model ParametersModel M disk N s ∆ t − Integrator( M (cid:63) ) ( T − min )(1) (2) (3) (4) (5)M1 10 − MERCURIUS
M2 5 × − MERCURIUS
M3 10 − MERCURIUS
M4 5 × − MERCURIUS
M5 10 − MERCURIUS
M6 10 −
100 21
MERCURIUS
M7 5 × −
100 21
MERCURIUS
M8 10 −
100 21
MERCURIUS
M9 5 × −
100 21
MERCURIUS
M10 10 −
100 21
MERCURIUS
M11 10 − MERCURIUS
M12 5 × − MERCURIUS
M13 10 − MERCURIUS
M14 5 × − MERCURIUS
M15 10 − MERCURIUS
M16 10 − MERCURIUS
M17 10 − IAS-15
M18 10 − EOS
M19 10 − JANUS
M20 10 − LEAPFROG
A summary of simulations carried out in this work.The default model is M3. In the column of ∆ t − , asmall correction of εT − min = 10 − T − min is dropped.Our baseline integrator is MERCURIUS , which em-ploys the Wisdom & Holman (1991) mapping.
IAS-15 is a high accuracy non-symplectic integra-tor with adaptive timestepping (Rein & Spiegel2015);
EOS corresponds to the Embedded Opera-tor Splitting Methods;
JANUS is an integer basedintegrator (Rein & Tamayo 2018); and
LEAPFROG is a symplectic integrator that does not require aKepler solver. bid symplectic intergator
MERCURIUS (Rein et al. 2019),based upon the widely-used
Mercury6 software package(Chambers 1999). By default, the systems are run for6 × time units ( ∼ . M disk and N s . We also test therobustness of this self-string for different time steps andintegrators. A representative list of models is summa-rized in Table 1.In order to avoid examining the dynamics of eachsimulated particle individually, we quantify the results of our simulations in terms of the angular momentumdeficit (AMD) (Laskar 1997, 2000):AMD = N (cid:88) j =1 ( L j − G j ) = N (cid:88) j =1 L j (cid:16) − (cid:113) − e j (cid:17) , (1)where L j = µ j (cid:112) G ( M (cid:63) + m j ) a j and G j = L j (cid:113) − e j .In our situation, the reduced mass µ j ≈ M disk /N s (cid:28) M (cid:63) , and L j ≈ M disk (cid:112) G M (cid:63) a j /N s . Importantly, low val-ues of this quantity (1) correspond to near-circular or-bits, while high eccentricities give large AMD and smallΘ, implying unfavorable conditions for planet formation. RESULTSLet us begin with an illustrative example. That is,while all of our models display a certain degree of self-stirring among planetesimals in the disk, the tendencytowards self-excitation is more pronounced in more mas-sive disks. Correspondingly, as a demonstration of thedescribed collective behaviour, in Figure (1) we show aninitial as well as an evolved orbital states of a M disk =10 − M (cid:63) , N s = 1000 disk. More specifically, Figure (1)depicts the starting state of the disk with purely circularorbits on the left panel and a final state with overlappedeccentric orbits on the right panel, where the default in-tegration timescale was increased by a factor of eight toaccentuate the growth of eccentricities.Of course, disks of solid debris as massive as that con-sidered in Figure (1) are unlikely to be physical withinthe broader context of planet formation. Given thatgravitational stability limits the mass of quasi-Kepleriandisks from above to a value smaller than their aspectratio M disk /M (cid:63) (cid:46) h/r ∼ .
05 and that the typical dust-to-gas ratio of circumstellar nebulae is of order 1%, weadopt M disk /M (cid:63) = 10 − as a reasonable mass-scale forour fiducial experiment, M3. The scaled AMD evolution(see below) obtained within this simulation is shown inthe top left panel of Figure (2) as a black curve withdash-dotted line. Dependence on Time — An immediately notable feature ofthe depicted time-series is that the growth of the scaledAMD is approximately linear. In fact, such behavior canbe expected if the individual eccentricities themselvesincrease due to stochastic forcing (see e.g., Puranam &Batygin 2018). This can be understood as follows. First,we note that for small eccentricities, the simplification1 − √ − e ≈ e / e is “dif-fusive” (see e.g., Øksendal 2003) such that d e ∝∼ d W (Wiener process) yielding e ∝∼ √ t . This presumptionimmediately gives AMD ∝∼ e ∝∼ t for small eccentrici-ties. x y Initial Condition: M disk = 10 M , N s = 1000 x y Snapshot at t = 0.8 Myr Figure 1.
The initial (left) and final (right) states of model M1. The small bodies start with perfectly circular orbits, andend with overlapped eccentric orbits, despite the fact that direct gravitational interactions among the particles are suppressed.990 of the small bodies and their orbits are displayed in lightgray, and 10 in black with orbits whose eccentricities evolve above99 th percentile in the final stage. For the purposes of this example, the integration time was extended to ∼ . Dependence on Mass — The dependence of AMD growthupon disk mass can be reasoned out in a similar fashion.Particularly, if we postulate that d e/ d t ∝ M disk /M (cid:63) ,then e ∝∼ ( M disk /M (cid:63) ) t and a reasonable choice ofAMD. evolution scaling would be AMD ∝∼ M disk e ∝∼ M disk ( M disk /M (cid:63) ) t . If correct, then by defining thescaled AMD asAMD s = AMD G (0) (cid:18) M (cid:63) M disk (cid:19) , (2)it should be possible to collapse the time evolution ofmodels with identical N s but with different M disk ontoa common curve. In other words, to remove the en-visioned cubic dependence of AMD growth on M disk ,we divide the AMD by the initial angular momentum G (0) = (cid:80) Nj =0 G j ( t = 0) ∝ M disk , and further scale it bythe square of disk-to-star mass ratio, rendering AMD s dimensionless.To test this assertion, we show AMD s growth for aseries of N s = 1000 models (Figure 2 top right panel),spanning M disk = 10 − − − M (cid:63) (M1 – M5). The indi-vidual numerical experiments are marked with differentline styles. By and large, these N s = 1000 simulationsexhibit approximately linear growth in scaled AMD, andthe depicted curves have comparable slopes (to withina factor of ∼ s growth. In- terestingly, the time series of the experiments with twosmallest disk masses (M4, M5) also tend to overlap fairlywell, but have growth rates that are notably smaller.Although of some interest, chasing down the associatedcorrection to equation (2) is beyond the scope of ourexploratory paper. Instead, we now turn our attentionto the dependence of this effect upon the “resolution” ofour experiments, N s . Dependence on N s — Although N s ∼ N s = 100 (M6 – M10) and N s = 10 , N s = 100 exhibitmore rapid and more uneven AMD s growth than theirhigher- N s counterparts. To this end, the concavity ofthe highest-AMD s N s = 100 models can likely be at-tributed to the fact that the attained eccentricities areso high that the reasoning behind quasi-linear growthoutlined above no longer applies. Conversely, curvescorresponding to N s = 10 ,
000 are smooth, linear, and A M D s Fiducial simulation M3 M disk / M = 1 × 10 N s = 1000 e ( M disk M ) ( N s ) ( t kyr ) A M D s M disk / M = 10 M disk / M = 5 × 10 M disk / M = 10 M disk / M = 5 × 10 M disk / M = 10 A M D s N s = 1000 N s = 100 N s = 10000 0 20 40 60 80Time (kyr)0100200300400500600700 A M D s MERCURIUSMERCURIUS finerIAS-15EOSJANUSLEAPFROG0 20 40 60 80Time (kyr)10 | E / E ( ) | | G / G ( ) | Figure 2.
Results of our numerical experiments. The time series correspond to the models listed in Table 1. Different linestyles indicate different M disk /M (cid:63) , and different line colors indicate analytic approximation, different N s or numerical setups.The first two rows correspond to dependencies on time (upper left), mass (upper right), N s (middle left), and numerical setup(middle right). The third row shows energy (left) and angular momentum (right) errors. overlap one-another very well, implying that equation(2) constitutes a better approximation for simulationswith higher N s .In addition to the aforementioned experiments,we have also measured the characteristic Lyapunovtimescale of N s = 1000 experiments and found that itdecreases approximately as the inverse square root of thedisk mass. Specifically, for a M disk = 10 − M (cid:63) system, T l ∼
250 yr; for M disk = 10 − M (cid:63) , T l ∼
80 yr; and for M disk = 10 − M (cid:63) , T l ∼
20 yr. Our M disk = 10 − M (cid:63) sim-ulations further indicate that the Lyapunov timescaleexhibits ancillary dependence on the particle count,with N s = 100, N s = 10000 runs yielding T l ∼
250 yrand T l ∼
130 yr, respectively.Cumulatively, the results of these experiments areconsistent with an interpretation wherein the spuriousgrowth of the angular momentum deficit is driven byperturbations that small particles exert upon the cen-tral body, which are then transmitted to other membersof the system. In other words, the gravitational cou-pling we observe in our numerical experiments is likelyfacilitated in full via the indirect terms of the disturb-ing Hamiltonian (see Ch. 6 of Murray & Dermott 1999),since there are no other interaction terms in the code. Itfurther worth noting that all indirect terms of the dis-turbing function average out to zero in the secular limit(where perturbations are taken to be phase-averaged),and indeed, this is the limit we approach as N s → ∞ ,which explains why the rate of AMD s growth dimin-ishes with increasing N s . Specifically, we found thatAMD s ∝∼ N − s via regression for N s ∼ − Dependence on Timestep & Integration Method — While theaforementioned dependencies of the collective disk be-havior on t , M disk , and N s appear sensible, the abovediscussion falls short of addressing the possibility thatthe dynamical behavior observed in our simulations isnothing more than a numerical artifact. Thus, as a fi-nal check on our results, we have repeated our fiducial M disk = 10 − M (cid:63) , N s = 1000 simulation employing avariety of numerical setups. Specifically, we test the de-pendence of the observed behavior on timestep (M16,M17), as well as integration method (M17 – M20). Thegrowth of AMD s is depicted on the middle right panelof Figure (2). Importantly, all of these numerical ex-periments yield consistent results, insinuating that theobserved dynamical excitation is genuine, and is not afeature of any specific algorithm. We have further usedthe Mercury6 software package (Chambers 1999) to re-produce some of our results using the hybrid and Bu-lirschStoer algorithms (Press et al. 1992), as well as the
IAS-15 integrator to verify the dependence on M disk il- lustrated in the left panel and got good agreement in allcases. A HEURISTIC EXAMPLEHaving demonstrated that indirect gravitational cou-pling among “small” non-interacting particles is ageneric feature of N -body simulations, we are now in aposition to inquire if this effect is of appreciable prac-tical importance in real planet-formation calculations.To answer this question, we proceed by considering aspecific example of post-nebular dynamical evolutionalready mentioned in the introduction: the assembly ofterrestrial planets from a narrow annulus of rocky de-bris. In particular, we follow Hansen (2009) and buildour terrestrial planet formation experiment by initializ-ing a M disk = 6 × − M (cid:63) disk of planetesimals confinedbetween 0 . ∼ . π × time units.The relevant time-series of this simulation are sum-marized in Figure (3). Intriguingly, these results showno sustained self-stirring in the system. Instead, con-trary to the numerical experiments reported in the pre-vious section, AMD s exhibits only low-amplitude fluctu-ations, and stays below AMD s (cid:46) − (upper panel), asdoes the eccentricity distribution (middle panel), with (cid:104) e (cid:105) (cid:46) − . We solidify this conclusion by repeating theexperiment with different random configurations and al-ternate integrators, as well as with a smaller number ofbodies ( N s = 400).In our interpretation, the disparity between this ex-periment and those described above lies in that here, M disk is so low, that the indirect gravitational stirringfalls below machine precision. In other words, the effectwe describe herein operates only above a threshold massof the small particle swarm. To test this assertion, wehave repeated the Hansen (2009) once again, boostingthe disk mass to M disk = 1 × − M (cid:63) as in Figure 1,and observed growth of the angular momentum that isfully consistent with the results depicted in Figure 2. In-deed, it is likely that the precise value of the thresholdmass above which indirect self-excitation ensues is botha function of N s as well as other details of the physicalsetup of the simulation such as the radial extent of thedisk, surface density profile, etc. SUMMARYIn this work, we have considered the dynamical conse-quences of big-small particle categorization scheme em- A M D s e ( ) th percentile75 th percentile50 th percentile0 200 400 600 800 1000time (Myr)10 | E / E ( ) | Figure 3.
Evolution of terrestial planet forming annulus of debris, with direct interactions suppressed. The upper panel showsthe evolution of AMD s with no clear growth. The middle panel shows the evolution of three percentiles (100 th , 75 th and 50 th )in the eccentricity distribution. The lower panel gives the energy error of the integration. ployed in conventional N -body simulations of planet-forming disks. To this end, we have carried out a seriesof numerical experiments that demonstrate that evenin absence of any big non-central bodies, interactionsamong massive small particles can still yield self-stirringwithin the system. We argue that this mode of dynami-cal excitation arises from indirect gravitational coupling,wherein perturbations are transmitted among particlesvia the barycentric reflex motion of the central star, in-duced through a superposition of individual Keplerianorbits.Collectively, our simulation suite shows that the afore-mentioned effect yields a growth of the systems angularmomentum deficit that is approximately linear in time,and scales roughly as the cube of the cumulative diskmass. These results are consistent with a picture wherethe evolution of the individual eccentricities is drivenby stochastic fluctuations, whose amplitude scales lin-early with the disk mass (such that the diffusive progressof the eccentricity dispersion has the approximate form (cid:104) e (cid:105) ≈ ( M disk /M (cid:63) ) (1000 /N s ) − / (cid:112) . t/ N s disks displaying less rapid AMD s growth.Finally, we have examined the role played by spuriousexcitation of the orbital dispersion within the contextof the solar systems terrestrial planet formation simu-lations (Hansen 2009; Walsh et al. 2011). Remarkably,we found no sustained growth of the velocity disper-sion arising from indirect interactions among small par-ticles, further demonstrating that this effect only op-erates above a certain threshold mass-scale which theterrestrial planet-forming annulus does not reach. As aresult, we conclude that although some caution may bewarranted in simulations of massive planetesimal disks,it is unlikely that interactions among non-interactingparticles within N -body simulations constitute a signif-icant source of uncertainty in numerical models of plan-etary assembly.We are thankful to Matthew J. Holman, Juliette C.Becker, Marguerite Epstein-Martin, Max Goldberg, To-bias Koehne and Elizabeth Bailey for insightful dis-cussions. Additionally, we would like to thank HannoRein for providing a thorough and insightful referee re- port which led to a considerable improvement of themanuscript, as well as his assistance with implemen-tation of numerical experiments. K.B. is grateful tothe David and Lucile Packard Foundation and the Al-fred P. Sloan Foundation for their generous support.Simulations in this paper made use of the REBOUNDcode which is freely available at http://github.com/hannorein/rebound.REFERENCES-body simulations constitute a signif-icant source of uncertainty in numerical models of plan-etary assembly.We are thankful to Matthew J. Holman, Juliette C.Becker, Marguerite Epstein-Martin, Max Goldberg, To-bias Koehne and Elizabeth Bailey for insightful dis-cussions. Additionally, we would like to thank HannoRein for providing a thorough and insightful referee re- port which led to a considerable improvement of themanuscript, as well as his assistance with implemen-tation of numerical experiments. K.B. is grateful tothe David and Lucile Packard Foundation and the Al-fred P. Sloan Foundation for their generous support.Simulations in this paper made use of the REBOUNDcode which is freely available at http://github.com/hannorein/rebound.REFERENCES