Superparticle Method for Simulating Collisions
David Nesvorny, Andrew N. Youdin, Raphael Marschall, Derek C. Richardson
aa r X i v : . [ a s t r o - ph . E P ] A p r Superparticle Method for Simulating Collisions
David Nesvorn´y , Andrew N. Youdin , Raphael Marschall , Derek C. Richardson (1) Department of Space Studies, Southwest Research Institute, Boulder, CO, USA(2) Steward Observatory, University of Arizona, Tucson, AZ, USA(3) Department of Astronomy, University of Maryland, College Park, MD, USA ABSTRACT
For problems in astrophysics, planetary science and beyond, numerical simu-lations are often limited to simulating fewer particles than in the real system. Tomodel collisions, the simulated particles (aka superparticles) need to be inflatedto represent a collectively large collisional cross section of real particles. Here wedevelop a superparticle-based method that replicates the kinetic energy loss dur-ing real-world collisions, implement it in an N -body code and test it. The testsprovide interesting insights into dynamics of self-gravitating collisional systems.They show how particle systems evolve over several free fall timescales to formcentral concentrations and equilibrated outer shells. The superparticle methodcan be extended to account for the accretional growth of objects during inelasticmergers.
1. Introduction
Numerical simulations are often limited in their ability to simulate statistically largeparticle systems. For example, a protoplanetary disk simulation can account for ∼ -10 particles (and grid cells; e.g., Li et al. 2019), but the actual number of real-world elements(dust grains, boulders, etc.) is vastly larger. The simulated particles have more mass andrepresent a large number of real particles. We call this the superparticle method. Thesuperparticle method has general applicability in many areas of science. Here we discuss itin the context of planet formation (see Youdin & Kenyon (2013) for a review) mainly becauseour previous studies of the initial stages of planet formation directly motivated this work.During the earliest stages of planet formation, small grains condense in a protoplanetarynebula and grow to larger ice/dust aggregates by sticking to each other. The growth stallsnear cm sizes because the electrostatic forces are not strong enough to hold large particlestogether. Also, large grains feel strong aerodynamic drag from the surrounding gas and 2 –drift toward the central star on timescales too short for significant growth to happen (e.g.,Birnstiel et al. 2016, but see Michikoshi & Kokubo 2016). It is therefore quite mysterious howplanet formation bridges the gap between cm-size particle aggregates (hereafter pebbles ) and1–1000 km bodies ( planetesimals ; Chiang & Youdin 2010, Johansen et al. 2014). A growingbody of evidence now suggests that an aerodynamic interaction between pebbles and nebulargas (e.g., the streaming instability, Youdin & Goodman 2005) can collect particles in largeself-gravitating clouds that then directly collapse into planetesimals.The initial stages of particle concentrations in a gas nebula can be studied with spe-cialized hydrocodes (e.g., ATHENA with the particle module of Bai & Stone 2010). Simu-lations show that the streaming instability (hereafter SI) should be particularly important(Youdin & Johansen 2007; Johansen & Youdin 2007; Johansen et al. 2009; Nesvorn´y et al.2019; Li et al. 2019). The instability occurs because an initially small over-density of peb-bles accelerates the gas. This perturbation launches a wave that amplifies pebble densityas it oscillates. Strong pebble clumping eventually triggers gravitational collapse into plan-etesimals. Modern hydrocode simulations of the SI, however, do not have adequate spatialresolution to follow the gravitational collapse to completion. Moreover, once the pebbledensity within a collapsing cloud exceeds the gas density, aerodynamic effects of gas ceaseto be important, and detailed hydrodynamic calculations are no longer required (Nesvorn´y,Youdin & Richardson 2010; hereafter NYR10). Instead, one has to realistically model pebblecollisions that damp random speeds and stimulate growth (Wahlberg Jansson & Johansen2014).In our previous work (NYR10), we studied gravitational collapse with a cosmological N -body code known as PKDGRAV (Stadel 2001).
PKDGRAV is a scalable, parallel tree code thatis the fastest code available to us for these type of calculations. A unique feature of
PKDGRAV is the ability to rapidly detect and realistically treat collisions between particles (Richardsonet al. 2000). In NYR10, individual
PKDGRAV particles were artificially inflated to mimic avery large collisional cross section of pebbles. We simply scaled up the radius of
PKDGRAV particles by a multiplication factor that was the same for all particles and unchanging withtime. This is not ideal for several different reasons. Crucially, the method with fixed inflationfactor does not correctly account for the kinetic energy loss during inelastic collisions.Here we extend the superparticle method to be able to more realistically model plan-etesimal formation (Section 3). The new method is designed to reproduce the kinetic energyloss. It will be useful in many areas of science where particle collisions and energy dissipationare important. The method is tested for self-gravitating collisional systems in Section 4. Toprovide some background for these tests, we discuss several relevant timescales in Section 2.Section 5 summarizes our findings. 3 –
2. Collision and collapse timescales2.1. Collisions
Here we define the collisional timescale, t coll , as the timescale on which any particle in acloud would, on average, have one collision with another particle. The collisional timescaleis then t coll = 1 /ησv , where η is the number density of particles, σ = πr is the particlecross-section, r is the particle radius, and v is the relative speed. Adopting the virial speed v = ( GM/R ) / , where G is the gravitational constant and M and R are the mass and radiusof the cloud, we obtain t coll = 43 nr √ GM R / , (1)where n is the total number of cloud particles. The collisional timescale is therefore a steepfunction of cloud’s size. We assume that the cloud radius is some fraction of the Hill radius, R = f H R H , where R H = a (cid:18) M M ⊙ (cid:19) / , (2) a is the orbital radius, and M ⊙ ≃ × g is the solar mass. After substituting for R H and n = ( R eq /r ) , where R eq is the equivalent radius of a sphere with mass M and density 1 gcm − , we obtain t coll ≃
23 yr (cid:16) r (cid:17) (cid:18)
50 km R eq (cid:19) (cid:18) f H . a
45 au (cid:19) / . (3)For comparison, the orbital period at 45 au is P ≃
300 yr. This illustrates the importanceof particle collisions. Collisions will act to damp particle speeds and stimulate growth.
In virial equilibrium, T = − U/
2, where T = P i m i v i / U = − P i,j Gm i m j /r ij ,with r ij = | r i − r j | , are the total kinetic and potential energies of particles with masses m i ,positions r i and speeds v i . If a cloud is in virial equilibrium, it cannot collapse unless somedissipative process, such as inelastic collisions between particles, reduces the kinetic energy.In a single collision of two pebbles with masses m ′ i and m ′ j , the change in energy is δE = − µ ∆ v (1 − C ) , (4)where µ = m ′ i m ′ j / ( m ′ i + m ′ j ), ∆ v is the collision speed and C R is the coefficient of restitution.The above equation takes into account that collisions are not necessarily head-on, whichreduces the average dissipated energy by a factor of 1/2. 4 –Assuming that all particles have the same mass and ∆ v = √ v , where v = ( GM/R ) / is the virial speed, the total energy lost in the time interval ∆ t is∆ E = δE × n ∆ tt coll , (5)where n is the number of particles and t coll is given in Eq. (1). This leads to a differentialequation for the total energy of the cloud, E , or equivalently, for the cloud radius. WahlbergJansson & Johansen (2014) showed that E = E (cid:18) − tt vir (cid:19) − / , (6)where t is time and E is the initial energy. This assumes that frequent collisions instanta-neously virialize the cloud. The virial collapse timescale t vir is given by t vir = 0 . t coll (1 − C ) − (7)with t coll in Eq. (1). The virial collapse timescale is thus roughly four times shorter for C R = 0than the collisional timescale for t = 0. This is a consequence of the steep dependence of t coll on R . In the virial collapse, the cloud radius R is related to Eq. (6) via R = 3 GM E . (8)
The above analysis is valid when the collapse is “hot”, i.e., when particle speeds arevirial and remain virial during the whole collapse. An alternative to this is the “cold” collapsewith sub-virial particle speeds. In this case, the collapse timescale is related to the free falltimescale t free = (cid:18) π Gρ cloud (cid:19) / , (9)where ρ cloud = 3 M/ πR is the cloud’s mass density. This assumes that particles havenegligible initial velocities and collisions are ignored. If the mass is uniformly distributedwithin a spherical volume of radius R , the collapse is self-similar and all radial shells shrinkon the timescale given by Eq. (9). If the mass is initially concentrated toward the cloud’scenter, the inner shells will infall faster than the outer shells.Substituting ρ cloud and R = f H R H into Eq. (9), we obtain t free = 11 yr (cid:18) f H . a
45 au (cid:19) / , (10) 5 –which is independent of M (or R eq ), because more massive clouds have larger Hill radii andthese dependencies cancel out in ρ cloud . Thus, in our fiducial case with f H = 0 . a = 45 au, r = 1 cm, R eq = 50 km and C R = 0, we obtain t vir < t free , which shows that collisions dampthe random velocities faster than the cloud can collapse. The collapse will thus happen onthe free fall timescale. Given the different scaling of Eqs. (3) and (10) with R eq , only smallclouds with R eq <
25 km can, at least initially, contract on the t vir timescale. As t coll dropsfaster during the collapse than t free , however, all clouds will eventually free fall.
3. Superparticle method
Two particle systems are considered here: one that consists of real particles (RPs)and another one that consists of simulated superparticles (SPs). Inelastic collisions betweenparticles happen in both systems. They can result in inelastic bounces or mergers of particles,and the loss of kinetic energy. We assume that the number of SPs is much smaller thanthe number of RPs, and ask how to best deal with collisions in the SP system such as itstatistically reproduces the behavior (i.e., energy loss, particle growth) of the RP system.
We denote the number distributions of SPs as dN and RPs as dN ′ . In terms of theRP radius r (no prime here) we may have a size distribution dN ′ /dr ∝ r − q . The numberdistribution of SPs can (i.e., might need to) be different, e.g. dN/dr ∝ r − Q . For instance, Q = 1 would be a log-uniform distribution that might be a good choice to numericallysample.We define n ( r ) as the number of RPs replaced by a single SP, or (equivalently) the massratio of a SP and its RP. The mass in each bin is the same in both systems, dM = m ( r ) dN = m ′ ( r ) dN ′ = n ( r ) m ′ ( r ) dN , where m ′ ( r ) is the mass of a RP. Thus n ( r ) = dN ′ /dN = m/m ′ (11)and for the example of power-law distributions n ( r ) ∝ r Q − q and n ( r ) ∝ r − q for log-uniformSPs. With typical q ≈ n ( r ) should be much larger for The considerations ignore the accretional growth of particles during collisions, which acts to increase r , t coll and t vir . Particle fragmentation would have an opposite effect (Wahlberg Jansson et al. 2017). Note that dN/dr is not the same as the distribution dN/dR , where R denotes the SP size; the relationshipbetween RP’s r and its SP’s R is not specified at this point. m ′ ∝ r , our log-uniform example meansthat SPs for smaller r will be more massive for q > q > Q for other SP power-laws).We will refer to this pathological case as SP mass inversion.More generally, for non-constant n ( r ), the mass ratios of SPs will not be the same as theRPs they model. This may have some undesirable properties for gravitational and collisionaldynamics. In particular, both gravitational scattering and (inelastic) collisions tend towardsequipartition, i.e., equal kinetic energies for all species. If SPs do not have the same n , thenwe risk forcing equipartition of the artificial SP masses. Since varying n ( r ) may be difficultto avoid, the question is whether this tendency to equipartition is significant, or whetherit is over-ridden by other physics (e.g., the fact that collisions are inelastic). Moreoverthe timescale for equipartition could be long compared to timescales of interest (e.g., forgravitational collapse). Intuitively the system of SPs must have larger collisional radii than RPs to accountfor the reduced number. Specifically, we want the energy loss rate from collisions per unitvolume, ˙ E , to be the same in the RP and SP systems. For collisions between species i and j the differential loss rate is d ˙ E ij = 14 dN ( r i ) dN ( r j ) σ ij v ij µ ij f ( C R ) (12)for cross section σ ij , relative speed v ij , reduced mass µ ij = m i m j / ( m i + m j ) and coefficientof restitution C R . Averaged over impact angles, f ( C R ) = (1 − C ) / E is obtained by a double integral over the size distributions for species i and j . Inmore detail, the above equation is already averaged over a velocity distribution and impactangles, and the dN ’s are now per unit volume. The 1/4 in Eq. (12) represents 1/2 fromkinetic energy and 1/2 to avoid double counting of species pairs.We now equate our (unprimed) SP system with the (primed) RP system, d ˙ E i,j = d ˙ E ′ i,j .We assume v ij = v ′ ij and f ( C R ) = f ′ ( C R ), i.e., that the SPs have a similar velocity distribu-tion to the RPs and there is a similar distribution of impact angles. From Equations (11)and (12) we get the desired result: σ ij = σ ′ ij m i + m j m ′ i + m ′ j = σ ′ ij n i m ′ i + n j m ′ j m ′ i + m ′ j , (13)i.e., the cross section to mass ratio is the same for SP and RP collisions. 7 –If SPs i and j represent the same number of RPs, n i = n j , then σ = σ ′ n , i.e., the effectiveSP radii are R i = √ nr i and R j = √ nr j . For strongly unequal true masses m ′ i ≫ m ′ j , σ ij m ′ i ≫ m ′ j −−−−→ πr i m ′ i ( m i + m j ) . (14)If the SP masses are not inverted, such that m i > m j , we again get R i ≈ p m i /m ′ i r i = √ n i r i .Moreover, if we define (now relaxing the assumption of strongly unequal masses) σ ij ≡ π ( R i + R j ) ≡ π ( √ n i r i + R j ) , (15)then R j > R i ≈ √ n i r i is “safe” in that it won’t require negative radii for the collisional partner, i.e.the other SP. If SPs merge, this represents a merger of some large number of RPs. If SPs i and j merge to a new SP k , the masses simply add m k = m i + m j . However there is some choice asto how the underlying true particle masses combine. We would like to exploit this choice sothat SPs gradually approach RPs as coagulation proceeds. The final collapse outcome willthen have more realistic densities and a more meaningful final multiplicity (single, multiple,etc.). The obvious choice of simply merging the true masses, as m ′ k = m ′ i + m ′ j , doesn’t allowthis transition, as we now explain in more detail.Our desired transition requires a gradual decrease in the sampling rates n k . The obviousmerging strategy clearly doesn’t have this property. If all SPs merge they have a total mass m tot = P i m i = P i n i m ′ i , while the simple sum of RP masses is m ′ tot = P i m ′ i (where i is now a sum over each SP and not each size bin). Clearly, the sampling rate of this finalobject n tot = P i n i m ′ i / P i m ′ i will be large if the initial n i are large, as this final n tot is justa mass-weighted average of the initial n i . Again, a decrease of the sampling rate is neededto approach a true system.During a merger event, we thus choose to reduce the merged sampling rate n k, merged andthus increase the RP masses (since the SP masses must be conserved). From the merger of m k = m i + m j , the post-merger sampling rate must obey n k, merged < m k m ′ i + m ′ j ≡ n k, , (16)where n k, would arise from the obvious summation of RP masses. For instance, in Section4.3, we consider n k, merged = max( n k, (1 − m j /m k ) ,
1) for m i > m j . (Note that this correction 8 –uses the SP masses, but we also have m ′ i > m ′ j for the RP masses, provided we stick to thenon-inverted samplings.)With different merger prescriptions, we should consider how the mass growth timescalesof SP and RP systems are related in theory. Recall that we are considering an SP cross sectionthat correctly models energy dissipation but not necessarily the mass growth. To explorethis issue, we consider the mass growth rate of a given SP in bin i due to the SPs in bin j :˙ m i,j m i = m j m i dN j σ ij v ij (17)and the corresponding RP growth rate (for simple addition of RP masses, i.e. before anyrescaling) : ˙ m ′ i,j m ′ i = m ′ j m ′ i dN ′ j σ ′ ij v ′ ij = m j m ′ i dN j σ ij v ij m ′ i + m ′ j m i + m j = ˙ m i,j m i (cid:20) m ′ j /m ′ i m j /m i (cid:21) . (18)Above we use our previous result for conservation of mass per bin ( m j dN j = m ′ j dN ′ j ), andagain use the cross section relation of Equation (13) assuming accurate modeling of velocities v ij = v ′ ij . Equation (18) shows that the mass growth timescales for the RP and SP systemswill agree in the special case where the term in square brackets is unity. This special conditionoccurs (1) for uniform sampling, i.e. n i = n j , which again is often overly restrictive and (2)in the limit of very small accreted masses j , both for the SP and RP masses. Aside fromthese special conditions, mass growth rates will not agree exactly. We thus confirm that acost of non-uniform SP sampling is that energy dissipation and mass growth can’t both bereproduced exactly.In what direction should the inexact growth rates go? Consider m ′ i > m ′ j and thestandard case n j > n i for more efficient sampling. Then Equation (18) gives ˙ m i,j /m i > ˙ m ′ i,j /m ′ i . In words, the large SPs will grow at a statistically larger rate than they “should”if they were accurately modeling the RP system. In practice, differences in the SP vs. RPrelative velocities will also affect growth rates. Thus we must compare RP and SP growthrates numerically to asses these issues quantitatively. Our post-merger reduction of the sampling rates n k doesn’t affect the SP masses, but rather the RPmasses. With decreasing n k the SP growth rates approach those of the RPs, but since the RPs have beenmodified, this is not the same result as the (too difficult) simulation of RPs from the outset. PKDGRAV implementation
The method described above was implemented in the hard-sphere flavor of
PKDGRAV (Richardson et al. 2000). To detect collisions,
PKDGRAV sorts neighbors of each particle bydistance. The particle and each neighbor are then tested for a collision during the nexttimestep d T . The condition is simply d ij ( t ) ≤ R i + R j , where d ij ( t ) is the distance of i and j particles at time 0 < t < dT , and R i + R j is the sum of particle radii. We changed thecollision condition in PKDGRAV such that R i + R j = p σ ij /π , where σ ij is given by Eq. (13).When a collision is detected, PKDGRAV decides, based on user-defined parameters, collisionspeed, etc., whether it results in a merger or bounce.If a merger is applied, we use the method described in Section 3.3 to generate a new SP.
PKDGRAV then computes the velocity of the new SP from the linear momentum conservation.Particle bounces in
PKDGRAV are controlled by the normal coefficient of restitution, 0 ≤ C R ≤
1, where C R = 0 corresponds to a fully inelastic collision and C R = 1 to an ideally elasticbounce (i.e., no energy dissipation). The code computes post-collision velocities of SP i and j and moves to considering a new SP pair. PKDGRAV is also capable of producing particleaggregates, where interacting particles become rigidly locked as a group, but we do not usethis option here.Specifically, we implemented the following changes in
PKDGRAV :1. We modified the structure of the input/output files such that n i and r i of each SP isgiven in the last two columns. New variables, fNpebble and fRpebble were includedin the PKD (particle data) and
COLLIDER (collision data) structures. The new variablesare passed between these structures in pkdGetColliderInfo and pkdPutColliderInfo (functions to retrieve and store collision data).
SSDATA SIZE (solar system data size)in ssio.h was increased to make space for new variables. The function pkdReadSS now reads these variables from input and passes it to the
PKD structure.2. We changed the function
CheckForCollision (function that tests particle’s neighborsfor collision in the current timestep) in smoothfcn.c such that the sum of particle radiiis sr = p σ ij /π , where σ ij is given in Eq. (13). We also modified option OverlapAdjPos (a method to deal with overlapping particles by separating them along their line ofcenters) in pkdDoCollision , included in collisions.c , such that the displacement isaligned with the definition of sr in CheckForCollision.3. We implemented the merger algorithm from Sect. 3.3 in function pkdMerge , included in collisions.c . This was done in two steps to account for the merger of i and j SPs, andto reduce the number of RPs. In pkdBounce , the individual particle radii were replaced 10 –by R i = r i q ( n i m ′ i + n j m ′ j ) / ( m ′ i + m ′ j ) and R j = r j q ( n i m ′ i + n j m ′ j ) / ( m ′ i + m ′ j ). Thisassures that the linear momentum conservation part in pkdBounce operates with theright quantities. The SP bounces in PKDGRAV are treated as bounces of particles withinflated radii.4. The expression for the escape velocity of colliding SPs in pkdDoCollision (functionthat executes collisions) was changed such that v = 2( m i + m j ) / p σ ij /π with σ ij from Eq. (13). PKDGRAV uses the escape velocity and the parameter dMergeLimit todecide whether the SPs should be merged.The modified PKDGRAV code works with the pthread parallelization and we typically used5, 10 or 28 Broadwell cores in the tests described in Sect. 4.
4. Tests
The highest-resolution SI simulations published to date produce self-gravitating cloudswith up to ∼ SPs and cloud masses corresponding to solid planetesimals from tens tohundreds kilometers accross. To focus on a specific scenario, consider a pebble cloud with N = 10 SPs and total mass M ≃ × g (corresponding to a solid planetesimal with R eq = 50 km and density ρ = 1 g cm − ), and assume that all pebbles have the same radius r ≃
10 cm (this is merely a convenient choice for the tests described below; pebbles in theouter solar system are expected to be at least ∼
10 times smaller, Birnstiel et al. 2016). Thetotal number of pebbles is then 3 M/ πr ρ ∼ . We thus see that the number of pebblesexceeds the number of SPs by ∼
11 orders of magnitude. From Section 3.2, the collisionalradius of each SP would need to be R SP ≈ √ cm ≃
32 km to have a correct rate ofenergy dissipation.The tests described below were designed to mimic the situation in a self-gravitating,collisional cloud of particles. Since we are unable to simulate the actual number of pebbles,we consider a reduced system with 10 RPs distributed in a spherical volume with f H = 0 . a = 45 au and R eq = 50 km. If RPs were given a solid density 1 g cm − , then r ≃ . t coll ≃ yr according to Eq.(1)). We therefore scale up the radii to have r = 32 km. This is done to reproduce thecollisional timescale in the 10 pebble cloud case discussed above. With 10 particles and r = 32 km, we have t coll ≃ . t free ≃ The effects of aerodynamic gas drag are not considered here. See Nesvorn´y et al. (2010) for a justification
11 –To test the SP method, we keep f H = 0 . a = 45 au and R eq = 50 km, and representthe particle cloud by 10 –10 SPs. This means that each SP initially represents 10–10 RPs.The results of these tests were compared to the fiducial case with 10 RPs. This is a goodway to test the validity of the SP method because we know what the ground truth is (asgiven by the 10 RP simulation). Specifically, we considered the energy dissipation, radialand velocity distributions, angular momentum conservation, etc. If mergers were applied,we also compared the size distributions.Initially, particles were uniformly distributed within a sphere of radius f H R H ≃ ,
000 kmand were given slightly sub-virial speeds (see below). We also considered other initial con-ditions such the isothermal radial profile from Binney and Tremaine (1987) or locally virialconditions, where T = − U/ We first consider a case with C R = 1 (fully elastic collisions) to verify the energyconservation in PKDGRAV . Particles were given random speeds v rand = 0 . − , which isslightly lower than the virial speed ( v ≃ . − ). This corresponded to T / | U | = 0 . PKDGRAV timestep was set to d T = 0 .
001 yr and the whole simulation covered10 yr, which is about ten free fall timescales.We tested cases with different
PKDGRAV tree opening angles ( dTheta =0.2, 0.5 and 1)and found that the fractional change of the total energy was 2 × − for dTheta = 0 . − for dTheta = 1. This result is expected since for smaller opening angles,cells must be further away to be treated in the multipole approximation, and the forcecalculation is more accurate. The simulations with 10 RPs were completed in 7 (for dTheta =0 . dTheta = 0 .
5) and 2 ( dTheta = 1) hours of processing time on 28 Broadwell cores.Decreasing the timestep to d T = 0 . T = 0 .
001 yr and dTheta = 0 . of this assumption in more realistic pebble collapse simulations.
12 –The system remained near the virial equilibrium with
T / | U | monotonically increasing fromthe initial 0.47 to the final 0.51.Figure 1a shows the radial mass distribution of the cloud. We divided the integrationdomain, represented by a sphere of radius 40,000 km, into 25 concentric shells of equal radialwidth, d R , and plot the mass in each shell as a fraction of the total mass ( M ≃ × g).Mathematically, the plotted quantity is d M ( R ) /M = 4 πρ R d R /M , where ρ is the massdensity at radius R . For the singular isothermal profile with ρ ∝ R − (Binney & Tremaine1987), d M ( R ) = const. Thus our mass profiles highlight differences relative to ρ ∝ R − .By design most mass is initially contained in the radial shells near the outer edge ofthe cloud. Over one free fall timescale the cloud contracts and the radial mass distributionhas a maximum near 15,000 km. This trend continues during the subsequent evolution untileventually, by the end of the simulation, the maximum mass fraction is near 7,000 km. Thisis a by-product of particle inflation. The total volume of the inflated particles is equivalentto a sphere of radius ≃ v rand = 0 . − . The velocities increase in the inner region during the free fall stage (Fig.1b). After that, collisions act to virialize the system and the radial velocity profiles becomeflatter (dashed line in Fig. 1b). Eventually, the whole cloud, except for its crowded innerregion, evolves toward the virial equilibrium at each radius (Fig. 2). The slightly sub-virialconditions in the outer region ( T / | U | ≃ . t = 10 yr; Fig. 2) can be explained if particlesare preferentially located near their orbital apocenters. The escape speed from individual RPs in our tests is only v esc = 4 . − , i.e., ∼
20 times smaller thanthe random speed. The effect of gravitational scattering during particle encounters is therefore negligible.Note that even with radius inflation more massive SPs exaggerate the effect of gravitational scattering. Thestandard radius inflation by √ n gives an escape speed that increases as n / for SPs that are n times asmassive as their RPs. This is not a problem for the tests presented here because v esc < v rand and thegravitational collapse timescale is relatively short. Gravitational stirring of SPs may become a problem inapplications that will require large values of n . In such cases, it would be desirable to soften the gravitationalforce/potential for neighbor SPs, but this will not completely solve the problem, because of the substantialcontribution to stirring from distant SP encounters. This is an important limitation of the superparticlemethod.
13 –Additional simulations were performed with the SP method described in Section 3. Tostart with, we used 10 SPs each representing n = 10 RPs. The simulations were done inthe elastic regime and were thus not meant to test Eq. (13). The total energy change anddeparture from the virial equilibrium were both found to be slightly larger than in the fiducialcase, probably due to the increased graininess of particle interactions. We also monitored thebehavior of the SP system to see whether it reproduced Figs. 1 and 2. We found that, indeed,the behavior was very similar and the radial profiles, random velocities and energies closelyfollowed the trends described above. More information on these comparisons is providedbelow. The results with 10 and 10 SPs were found to be nearly identical to those obtainedwith 10 RPs. The 10 SP case, however, already significantly deviated from the referencecase in that it showed a far weaker particle concentration near the center. This happenedbecause SPs were strongly inflated and reached maximum packing in the center well beforethe collapse was completed.
The main goal of the tests reported here was to verify whether the SP method describedin Section 3 correctly reproduces the energy loss. To this end we adopted C R < q = 3 and distributed RPs within a factor 4 in size around r ≃
35 km (i.e., the value used inthe equal-size case). The most massive RPs thus had 64 times the mass of the least massiveRPs. The SP method was tested for power-law distributions with Q = 1, 2 and 3 (Sect 3.1).We experimented with C R = 0 .
5, 0.8 and 0.9. The total energy loss from Eqs. (4) and(5) is expected to be 1.5, 0.7 and 0.4 ( × in cgs units per year; Eq. (6) does not apply,as we discussed in Sect. 2.3, because the collapse happens on the free fall timescale). Forcomparison, the energy loss in a yr-long interval of the simulation with 10 RPs was 1.6, 1and 0.6 ( × in cgs units), respectively (Fig. 3). This is reasonable. An exact match isnot expected because the equations describe an average behavior that doesn’t account forspatial variations (or the time evolution of these variations) in the cloud.The superparticle method with 10 SPs replicated the energy loss reasonably well(Fig. 3). The small difference between the solid and dashed lines in Fig. 3 is a conse-quence of a relatively small number of SPs. The results with 10 and 10 SPs more closelyreplicated the reference case. The energy loss was more random and a factor of ∼ SP case in the following tests, butcomment on the results with other sampling rates as well. The radial mass profiles obtained 14 –in the case with C R =0.5 (maximum dissipation tested here; Fig. 4a) are similar to thoseobtained in the elastic case (Fig. 1a) except that the particles are even more crowded towardthe center, and this trend becomes stronger with time. This is a consequence of the energyloss in inelastic collisions that damp random velocities and facilitate collapse. Eventually, asthe number of collisions increased beyond reasonable limits in the central region, the codewas unable to make any progress. We stopped all simulations at t = 1 yr.The virial parameter T / | U | is shown in Fig. 4b. The T / | U | ratio is generally smallerthan that in Fig. 2 at least partly because here the kinetic energy decays by dissipativecollisions. This becomes especially clear near the center as time approaches the free falltimescale. All these trends are well reproduced with the SP method (Fig. 5). As the SPsoccupy a larger volume, however, they are not allowed to concentrate as strongly toward thecenter as in the reference case. This leads to the mass distribution that peaks at somewhatlarger radial distance (compare Fig. 4a and Fig. 5a). This “packing” problem becomes less(more) of an issue with > SPs ( < SPs). In addition, the radial profiles become verynoisy with only 10 SPs, suggesting that this sampling rate is not adequate to capture thestatistical behavior of the real system.We verified that the total angular momentum L is conserved by adding solid initialrotation of the cloud around the z axis. At the cloud’s outer radius of ≃ v circ = p GM/R ≃ . − . We thus added angular speed Ω = f Ω v circ /R ,where f Ω ≥ f Ω = 0 . ≃ . × − s − . This choice was motivated by rotation inferred forself-gravitating pebble clouds from our SI simulations (e.g., Nesvorn´y et al. 2019, Li et al.2019), but note that the initial structure of SI clouds is generally complex and cannot beidealized by solid rotation of homogeneous spheres. Here we used the idealized case to testthe SP method. More realistic initial conditions will be investigated elsewhere. We foundthat the angular momentum was very well conserved in both the reference case (relativechange δL/L ∼ × − ) and when the SP method was used (e.g., δL/L ∼ × − for thecase with 10 SPs).
Here we tested the merger algorithm described in Section 3.3. The general setup ofthe simulations was the same as the one used above. We tested cases with ( f Ω = 0 .
5) andwithout ( f Ω = 0) cloud rotation. The initial rotation provides centrifugal support to thecloud and reduces problems with particle crowding. We found that the case with f Ω = 0 . f Ω = 0 because the results are less sensitive to various 15 –integration details (e.g., the seed used to generate the initial conditions). The results for f Ω = 0 are reported below for completeness. In PKDGRAV , particle mergers are applied if thecollision speed between particles v col < f esc v esc , where v esc is the escape speed and f esc < v col > f esc v esc result is particle bounces.Figure 6 shows the size distributions obtained in the simulations with 10 RPs and 10 SPs. In both cases, we used v rand = 0 . − , C R = 0 . f Ω = 0 and f esc = 0 .
1. In thereference case, we simply plot the size distribution of RPs. In the SP simulation, the sizedistribution is constructed from the number and sizes of particles represented by each SP.We found that the SP algorithm with 10 SPs can reproduce the particle growth quite wellbut the results had a rather large statistical variability with minor changes of the input(e.g., when different seeds were used to generate initial conditions). We interpret this as aconsequence of stochastic interaction of large bodies that grow in the center of the collapsingcloud. The stochasticity of the results diminishes if > SPs are used. For this reason, inthe rest of this section, we discuss the results obtained with f Ω = 0 .
5, as rotation inhibitscollapse to the center of mass. For reference, the radii of the largest bodies shown in Fig. 6are ≃ ≃ f Ω = 0 . ≃ RPs) and ≃ SPs).Both simulations ended with a rounded cumulative distribution function (CDF) profile anddozens of bodies larger than 1000 km. For f Ω = 0 .
5, the SP method somewhat overshoots thenumber of the smallest objects but the difference is not large (the case with f Ω = 0 producedan opposite result; Fig. 6). We repeated the SP simulations with slightly changed initialconditions and found that the results were quite similar to those shown in Fig. 7b. Thismeans that the stochastic variations are not as large for f Ω = 0 . f Ω = 0.The radial mass distribution profiles are compared in Fig. 8. The SP method with 10 SPsreasonably well replicates the profiles obtained in the fiducial case. The results with 10 and10 SPs reproduce the fiducial radial profiles and size distributions even more closely.Figure 9 shows the evolving structure of the particle cloud with f Ω = 0 .
5. First, anagglomerate of particles formed in the center. The agglomerate had significant rotation andstretched along one direction to a ≃ t = 4 yr, the extremesdetached to release the angular momentum whereas the middle part of the agglomerateremained connected. Several small binaries formed at this point. Eventually, when thesimulation was extended to t = 100 yr, a massive, equal-size binary formed in the center. Thetwo components of the binary had radii R = 24 . R = 20 . a B = 4530 km. Thus a B / ( R + R ) / ≃
5. Conclusions
We developed a new method that will be useful for modeling statistically large colli-sional systems of particles. The method represents a significant computational speed-upby employing far fewer “superparticles” than the number of particles in the real system.It reproduces the energy loss in dissipative collisions of real particles and can account forparticle growth by mergers. We implemented the method in the
PKDGRAV code and tested it.We found that at least ∼ superparticles need to be used to reproduce the results of ourfiducial simulations with full resolution. The number of superparticles needed in real-worldapplications will have to be determined by convergence studies tailored to those applications.As the superparticles need to be inflated to represent a much larger collisional cross-section ofthe real system, particle crowding may become a problem in some applications. In the grav-itational collapse simulation performed here, this becomes apparent near the cloud’s center,where the code has difficulties to account for a vast number of superparticle collisions. Thisproblem can be mitigated by using a larger number of (less-inflated) superparticles, mergingsuperparticles and/or letting PKDGRAV form rigid particle aggregates.The work of DN was funded by the NASA EW and XRP programs. ANY acknowl-edges support from NASA Astrophysics Theory Grant NNX17AK59G and from NSF grantAST-1616929. RM acknowledges the support of the Swiss National Science Foundation(SNSF) through the grant P2BEP2 184482. We thank the anonymous reviewer for helpfulsuggestions.
Author’s contributions:
DN motivated this study, performed and analyzed all testswith
PKDGRAV , and prepared the article for publication. ANY derived Eq. (13) and wroteSection 3. RM helped to develop tests that are reported in Section 4. DCR provided supportrelated to the
PKDGRAV code. All authors contributed to the interpretation of the results andwriting this paper.
REFERENCES
Bai, X.-N., & Stone, J. M. 2010, ApJS, 190, 297 17 –Binney, J., & Tremaine, S. 1987, PrincetonBirnstiel, T., Fang, M., & Johansen, A. 2016, Space Sci. Rev., 205, 41Grundy, W. M., Noll, K. S., Roe, H. G. et al. 2019, Icarus, 334, 62Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493Johansen, A., & Youdin, A. 2007, ApJ, 662, 627Johansen, A., Youdin, A., & Mac Low, M. 2009, ApJ, 704, L75Johansen, A., Blum, J., Tanaka, H., et al. 2014, Protostars and Planets VI, 547Li, R., Youdin, A. N., & Simon, J. B. 2019, ApJ, 885, 69Michikoshi, S., & Kokubo, E. 2016, ApJ, 825, L28Nesvorn´y, D., Youdin, A. N., & Richardson, D. C. 2010, AJ, 140, 78 5 (NYR10)Nesvorn´y, D., Li, R., Youdin, A. N., Simon, J. B., & Grundy, W. M. 2019, Nature AstronomyNesvorn´y, D., & Vokrouhlick´y, D. 2019, Icarus, 331, 49Richardson, D. C., Quinn, T., Stadel, J., et al. 2000, Icarus, 143, 45Simon, J. B., Armitage, P. J., Youdin, A. N., et al. 2017, ApJ, 847, L12Stadel, J. G. 2001, Ph.D. ThesisWahlberg Jansson, K., & Johansen, A. 2014, A&A, 570, A47Wahlberg Jansson, K., Johansen, A., Bukhari Syed, M., et al. 2017, ApJ, 835, 109Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459Youdin, A., & Johansen, A. 2007, ApJ, 662, 613Youdin, A. N., & Kenyon, S. J. 2013, Planets, Stars and Stellar Systems. Volume 3: Solarand Stellar Planetary Systems, 1 18 –Fig. 1.— The mass (a) and mean velocity (b) profiles of a self-gravitating collisional cloud.Initially, 10 particles with r ≃
35 km were homogeneously distributed in a spherical volumewith f H = 0 . a = 45 au and v rand = 0 . − . Here we used C R = 1 (fully elastic case).The solid and dashed lines show the initial ( t = 0) and final ( t = 10 yr) profiles. The coloredlines follow profiles at 0.1 yr increments from t = 0 to 1 yr. 19 – This preprint was prepared with the AAS L A TEX macros v5.0.
20 –Fig. 2.— The ratio of kinetic ( T ) and potential ( U ) energies as a function of radial distance.The line styles and parameters used here are the same as in Fig. 1. The potential energywas computed as U = − GM ( R ) / R , where M ( R ) is the mass of particles inside a sphere ofradius R . 21 –Fig. 3.— The change of total energy for a case with equal-size particles (panel a) andunequal-size particles (panel b). Different lines correspond to C R = 0 . particles. The dashed lines show the results obtained with 10 SPs. For the unequal-sizecase, we used q = 2 and Q = 1 (Section 3.1). The results with Q = 2 and Q = 3 are nearlyidentical to the ones shown here for Q = 1. We were unable to determine the energy lossbeyond roughly one free fall timescale, because particles became crowded in the center andthe simulations stalled just past t = 1 yr due to an excessive number of collisions. 22 –Fig. 4.— The mass (a) and energy (b) profiles of a self-gravitating collisional cloud ofparticles. The line styles and simulation parameters used here are the same as in Fig. 1,except in this case C R = 0 .
5. The random velocities are damped by inelastic collisionsresulting in a stronger particle concentration in the inner region (compare panel (a) withFig. 1a). The simulation stalls just after t = 1 yr due to an excessive number of collisionsin the inner region. 23 –Fig. 5.— The mass (a) and energy (b) profiles of a self-gravitating collisional cloud. Thesimulation setup is the same as in Fig. 4, except that here we used 10 SPs each representing n = 10 RPs. The radial bin size was increased to compensate for poor resolution (10 binshere vs. 25 bins in previous figures). Panel (a) can be compared to Fig. 4a and panel (b) toFig. 4b. 24 –Fig. 6.— The cumulative size distributions of particles for the case with C R = 0 . f Ω = 0. Mergers between particles occur in PKDGRAV when the collision speed is less then0.1 of the particle escape speed. The size distribution is shown for 0 < t <
10 yr in 1 yrincrements (from thin black to thick blue lines). The reference case with 10 particles isshown in panel (a). The case with 10 SPs is shown in panel (b). 25 –Fig. 7.— The cumulative size distributions of particles for the case with C R = 0 . f Ω = 0 .
5. See caption of Fig. 6 for additional information. 26 –Fig. 8.— The radial profiles for the case in Fig. 7. The profiles are shown for 0 < t < t > particles and 10 SPs,respectively. 27 –Fig. 9.— Four snapshots showing particle distributions at t = 0 (top left), 1.2 (top right), 3(bottom left) and 7 years (bottom right). This is for the case with C R = 0 . f Ω = 0 .5.