The Astrophysical Journal, 856:175 (14pp), 2018 April 1

https://doi.org/10.3847/1538-4357/aab369© 2018. The American Astronomical Society. All rights reserved.

Particle Number Dependence of the N-body Simulations of Moon Formation
1

Takanori Sasaki1

and Natsuki Hosono2,3

Department of Astronomy, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan; takanori@kusastro.kyoto-u.ac.jp
2
Japan Agency for Marine-Earth Science and Technology, 2-15, Natsushima-cho, Yokosuka-city, Kanagawa 237-0061, Japan
3
RIKEN Advanced Institute for Computational Science, 7-1-26 Minatojima-minami-machi, Chuo-ku, Kobe, Hyogo, Japan
Received 2017 December 27; revised 2018 February 18; accepted 2018 February 27; published 2018 April 5

Abstract
The formation of the Moon from the circumterrestrial disk has been investigated by using N-body simulations with
the number N of particles limited from 104 to 105. We develop an N-body simulation code on multiple Pezy-SC
processors and deploy Framework for Developing Particle Simulators to deal with large number of particles. We
execute several high- and extra-high-resolution N-body simulations of lunar accretion from a circumterrestrial diskof debris generated by a giant impact on Earth. The number of particles is up to 107, in which 1 particle
corresponds to a 10 km sized satellitesimal. We ﬁnd that the spiral structures inside the Roche limit radius differ
between low-resolution simulations (N105) and high-resolution simulations (N 106). According to thisdifference, angular momentum ﬂuxes, which determine the accretion timescale of the Moon also depend on the
numerical resolution.
Key words: methods: numerical – Moon – planets and satellites: formation
the gravitational torque due to the spiral structure formed
within the Roche limit radius around the Earth and by thecollective motion of the particles associated with the structure.
They concluded that the formation and evolution of the spiralstructure were regulated only by the disk surface density suchthat the angular momentum transfer depends on only thesurface density but not on the particle number N whenN 103.
However, note that no simulations exist that focus on how
the satellite mass and growth speed depend on the number ofparticles. At the very ﬁrst stage of the accretion of the Moon,
the radius of debris particles is thought to be approximately
10 km, which corresponds to numerical simulations withapproximately 107 particles.
To address this shortcoming, we present herein the results ofN-body simulations based on the rubble pile model withN=104, 105, 106, and 107 particles, which allow us to follow
the lunar accretion process until Moon formation. We
investigate how the numerical resolution affects the formationand evolution of the disk structure, including the spiralstructure. We also compare the results to those of Kokubo
et al. (2000) and Takeda & Ida (2001). We ﬁnd that the spiralstructures inside the Roche limit radius are the same as thosefound in previous studies based on low-resolution simulations(N105) but differ from those found in high-resolution
simulations (N 106). On the basis of these differences, we
ﬁnd that the angular momentum ﬂuxes, which determine theaccretion timescales of the Moon, also depend on the numericalresolution of the simulation.
We assume herein that the circumterrestrial disk is anentirely particulate disk. SPH simulations suggest that mostparticles are in vapor or liquid phase just after the giant impact,with an average temperature exceeding 4000 K (Canup 2004).
In such cases, lunar accretion proceeds on the cooling timescaleof the entire disk (see, e.g., Thompson & Stevenson 1988;
Ward 2012) or by the slow spreading of the inner melt disk
(Salmon & Canup 2012), both of which are much longer(∼100 years) than the accretion timescale predicted by N-bodysimulations (∼1 month). While the disk’s thermal evolution

1. Introduction
The giant impact hypothesis is a widely accepted scenario
for the origin of the Moon (Hartmann & Davis 1975; Cameron& Ward 1976) and has been favored because it can explain the
Moon’s mass and iron depletion, as well as the angularmomentum of the Earth–Moon system (Canup & Asphaug2001). The giant impact process has been investigated bysmoothed particle hydrodynamics (SPH) simulations (see, e.g.,
Benz et al. 1986; Canup 2004), which showed that such giant
impacts usually result in the formation of a circumterrestrial
disk whose mass is usually smaller than 2.5 ML, where ML isthe present lunar mass (∼0.0123 M⊕, where M⊕ is the present
Earth mass). Most of the disk material is distributed near orinside the Roche limit radius RRoche (∼2.9 R⊕, where R⊕ is thepresent Earth radius), if the orbital angular momentum of the
impactor is about the same as that of the present Earth–Moon
system.
Lunar accretion from the circumterrestrial disk was investigated by using N-body simulations (Ida et al. 1997; Kokubo
et al. 2000). Ida et al. (1997) modeled the circumterrestrial disk
by 1000–1500 particles. However, this value is insufﬁcient toresolve the spatial structure of the disk, which is important forangular momentum transfer and mass transfer in the disk.
Kokubo et al. (2000) performed higher-resolution N-bodysimulations by using 10,000 particles for the disk and foundthat the results coincided qualitatively and quantitatively withthat of Ida et al. (1997) for the same surface density of the disk.
To more clearly see the evolution of the spatial structure of thedisk, Kokubo et al. (2000) also performed an N-bodysimulation based on the “rubble pile model,” in which no
mergers were allowed. To date, this 10,000-particle N-bodysimulation is the only rubble pile model simulation to follow
the process of lunar accretion until Moon formation.
Conversely, Takeda & Ida (2001) evaluated the angularmomentum transfer processes in the circumterrestrial disks, as
they followed the lunar accretion process halfway. They
performed global N-body simulations based on the rubble pilemodel with particle numbers up to 105 to represent the disk.
They found that angular momentum transfer was dominated by
1

The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

would depend on how the disk made of a vapor–liquid–solid
mixture cools and spreads, the resulting mass of the Moon islargely determined by the initial mass and angular momentumof the disk produced by the giant impact rather than thesubsequent disk evolution (Salmon & Canup 2012). The
purpose of the present study is to evaluate how the lunar
accretion process depends on the number N of N-bodysimulations. Therefore, we do not pursue such long-term disk
evolution and consider only the same situations adapted inprevious studies (Ida et al. 1997; Kokubo et al. 2000; Takeda &Ida 2001).
Section 2 describes the model and calculation method.
Section 3 presents the results of the N-body simulations.
Section 4 is devoted to a summary and some discussions.

whereas aggregates are tidally disrupted when they stray too
close to the proto-Earth. For simplicity, we do not consider
fragmentation, heat generation, or vaporization. We assume
that the central body is spherical and do not include the tidal
force due to tidal bulges raised on Earth, because the timescaleof the tidal orbital evolution is much longer than the accretion
timescale. When particles fall inside the central body’s radiusR⊕, we assume they are accreted onto the Earth and remove
them from the calculation.
In addition to self-gravity, we must resolve collisions
between particles. Toward this end, Ida et al. (1997) assumed
that when two particles overlap, the particles always andimmediately merge into a single particle. Kokubo et al. (2000)
allowed the merger of two particles only if the Jacobi energy of
the two particles was negative after the collision; otherwise,
they rebounded. Note that merging particles would lead to an
overestimation of the ﬁnal satellite mass because the mass loss
due to tidal stripping and/or collisions between rabble piles isignored. Thus, we do not allow particle mergers in this work;
instead, all collisions between two particles are resolved byrebound.
When the distance between two particles is less than the sumof the particle radii, we change the particle velocities according
to the restitution coefﬁcient ò (0ò1),

2. Method of Calculation
We performed several sets of high- and extra-high-resolution
N-body simulations that included physical collisions and selfgravity to simulate lunar accretion from a circumterrestrial
debris disk generated by a giant impact on Earth. To investigate
how the lunar accretion process depends on the number N of
N-body simulations, we adopt physical and numerical models
similar to those used in previous studies (Ida et al. 1997;
Kokubo et al. 2000), except for the number of particles.
We start the simulations of lunar accretion by assuming a solid
particle disk, which forms through the condensation of thesilicate vapor–liquid generated by the giant impact (Canup 2004).
We set the density ρ⊕ of the proto-Earth to 5.5 g cm−3 (equal tothe present bulk density of Earth) and the density ρp of the disk
particles to 3.3 g cm−3 (equal to the present bulk density of theMoon). We scale the orbital radii by the Roche limit radius asfollows:⎛ rp ⎞-1 3
RRoche = 2.456 ⎜ ⎟RÅ  2.9RÅ ,⎝ rÅ ⎠

1 3
1 ⎛ m ⎞
⎟ RRoche ,⎜2.456 ⎝ MÅ ⎠

(1 )

(2 )

following Kokubo et al. (2000).
The orbits of disk particles are calculated by numericallyintegrating the equation of motion,Nxi - xjdvi
x,= - GMÅ i 3 - å Gmi∣ xi ∣∣xi - xj ∣
dt
j¹i

(4 )

v¢t =  t vt ,(5 )where v and v¢ are the relative impact velocity and the relativerebound velocity, respectively; subscripts n and t indicate thenormal and tangential components, respectively. We assume
that the tangential restitution coefﬁcient is òt=1, neglecting
any exchange between orbital angular momentum and spin
angular momentum; furthermore, we assume a normal restitution coefﬁcient òn=0.1, independent of collision velocity.
Takeda & Ida (2001) showed that angular momentum transfer
rates are almost independent of òn except for highly elastic casessuch as òn 0.6. The velocity of each particle after a collision is
calculated on the basis of the conservation of momentum. To
avoid numerical difﬁculty with a detailed adjustment of collisionand rebounding, we follow Richardson (1994).
The initial condition of the disk is the same as that for run 25
in Kokubo et al. (2000): the initial disk mass is Mdisk=4.0 ML,
and the surface density distribution of the disk is modeled by apower-law distribution Σ ∝ a−3, where a is the semimajor axisof disk particles with inner and outer cutoffs. The inner cutoff
is the Earth radius ain=R⊕, and the outer cutoff is the Roche
limit radius aout=RRoche. All particles have an equal-mass m
determined by the number of particles N (i.e., m=Mdisk/N).
The semimajor axes of particles are randomly chosen between
the cutoffs to produce the disk’s power-law distribution.
Orbital eccentricities and inclinations are given by the Rayleighdistribution, root mean square (rms) eccentricity áe 2ñ1 2 =0.3, and rms inclination ái 2ñ1 2 = 0.15. In general, the initial
distributions of eccentricities and inclinations do not affect disk
evolution because they relax owing to collisional damping on atimescale of approximately one Kepler period. Eccentricanomaly E, longitude of ascending node Ω, and argument of
pericenter ω are randomly chosen between 0 to 2π. Initial
positions and velocities of the particles are uniquely determinedby the above parameters (a, e, i, E, Ω, and ω).

where R⊕ is the Earth radius. We use the proto-Earth mass as
the present Earth mass M⊕=5.97×1027 g, and we assume adisk of equal-mass particles. The physical radius of disk
particles with mass m isrp =v¢n = -  n vn ,(3 )where x and v are the position and velocity of disk particles,
respectively, and G is the gravitational constant. For numericalintegration, we follow Kokubo et al. (2000) and use a modiﬁed
fourth-ordered Hermite scheme for planetary N-body simulation with a shared timestep of 2−9. A full description of thescheme is given in Makino & Aarseth (1992).
We adopt the rubble pile model for physical collisions inwhich no mergers are allowed and instead particles simply
bounce inelastically off each other. In this model, gravitationally bound aggregates form outside the Roche limit radius,
2The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono
4

In this work, we vary the number of particles from 10 to107, in which one particle corresponds to 10 km sizedsatellitesimals. To execute a run with an extra-large numberof particles, we deploy Framework for Developing Particle
Simulators (FDPS; Iwasawa et al. 2016), which automaticallyadopts the tree method (Barnes & Hut 1986) and parallelizes anarbitrary particle-based numerical simulation method. Furthermore, we used a Pezy-SC processor, which contains a large
number of coprocessors to perform massively parallelized
numerical simulations (see the Appendix).

low-resolution simulations show that lunar accretion is divided
into two stages: rapid growth (t40TK) and saturated growth(t>40TK). The proto-Moon is almost formed in the ﬁrst stage.
In the second stage, the proto-Moon not only terminates its
growth but also decreases its mass because of the collisional
and/or tidal strip caused by gravitational interactions with disk
material and/or other moonlets. Conversely, high-resolution
simulations show that the lunar seed grows continuously andgradually until the ﬁnal mass is reached at t∼100TK.
On the other hand, the accumulative masses accreted onto
the Earth are plotted as a function of time in Figure 8. The massof over 2 ML, which corresponds to a half of the initial disk
mass, has fallen to the Earth in 40 TK. The accreted masses
for low-resolution simulations (N105) are slightly larger
than those for high-resolution simulations (N 106). However, in contrast to the mass evolution of the proto-Moon, theconvergent behavior of accreted mass does not depend on the
numerical resolution.
The simulated mass of the largest aggregate, that of the sumof the three largest aggregates, and that fallen to the Earth att=110TK are summarized in Table 1 and Figure 9. Figure 9
indicates that both the Moon’s mass and the accumulative massfallen to the Earth would not depend on the numericalresolution. The Moon’s masses are smaller than those reported
by Kokubo et al. (2000) with the same disk parameters but
different accretion models (i.e., total and partial accretion
models). Considering that we can use the rubble pile model tosimulate processes such as collisional breakup and/or tidal
strip for a proto-Moon, the ﬁnal mass of the proto-Moon(Table 1) would be smaller than that obtained by Kokubo et al.
(2000) (Table 2). Although stripping the proto-Moon’s material
would not qualitatively change the mode of lunar accretion, it
can quantitatively change the ﬁnal lunar mass. Forming a Moon
body requires a longer time because of the re-accretion of the
stripped materials to the proto-Moon and/or a more massive
circumterrestrial disk made by the giant impact. However, notethat aggregates are too fragile to withstand collisional breakup
and tidal strip in the rubble pile model because the model
neglects the strength of the materials. Therefore, the actual
masses of the proto-Moon produced in the disk would beslightly greater than those listed in Table 1.

3. Results
3.1. Overview
We simulate the evolution of a circumterrestrial disk andlunar accretion and investigate how numerical resolution
affects these processes. First, we show the overall evolutionof the disk: Figures 1–4 show snapshots of the disk projected
onto the x−y plane at t=0, 1, 5, 10, 20, 40 TK, where we
scale the time by TK, which is the Kepler period at the Roche
limit distance. Figure 1 in this paper corresponds to Figure4 in
Kokubo et al. (2000) (i.e., N=104). For high-resolution
simulations (N 106), snapshots at t=60, 80, 100, 110 TK
are also shown in Figures 5 and 6. The point sizes in
Figures 1–6 are proportional to the physical size of disk
particles.
The evolution of the circumterrestrial disk for the simulation
with N=104 (Figure 1) is consistent with that found by
Kokubo et al. (2000). The disk contracts by the collisional
damping of particles, and particle clumps grow inside theRoche limit radius (t=TK). The clumps are elongated by
Keplerian shear and spiral arms are formed (t=5TK). Particles
are transferred beyond the Roche limit radius via thegravitational torque exerted by the spiral arms (t=10TK).
Small particles beyond the Roche limit radius collapse to form
small moonlets, the rapid accretion of which forms a lunar seed(t=20TK). The seed grows exclusively by sweeping up
particles transferred beyond the Roche limit radius, and theproto-Moon forms just outside the Roche limit radius(t=40TK). Although relatively large moonlets survive on
the horseshoe orbit with the proto-Moon, they collide with eachother over the long-term tidal interaction with Earth to accrete
to a single moon.

3.3. Disk Evolution
We now discuss in detail the evolution of the circumterrestrial disks. Figures 10–12 show snapshots of the surface density
Σ, the radial component vR of the velocity dispersion, andToomre’s Q value, respectively. The velocity dispersion is therms deviation velocity from the local Keplerian circularvelocity. Toomre’s Q value for a particle disk is approximately

3.2. Mass Evolution of the Proto-Moon
Figures 1–4 show the evolution of the disks, which are
qualitatively the same as that of Kokubo et al. (2000) regardlessof the numerical resolution. However, the spiral arm structure
differs between low-resolution simulations (N105) andhigh-resolution simulations (N 106), and this difference
results in a longer timescale for lunar accretion in the lattercase. For low-resolution simulations, individual spiral arms
appear clearly, and the proto-Moon is almost formed att=40TK. Conversely, for high-resolution simulations, thespiral arms connect with each other to form more complicatedstructures between the arms, and the proto-Moon forms at
approximately t=100 TK.
The mass evolution of the largest aggregate and that of thesum of the three largest aggregates are plotted as a function oftime in Figure 7. These graphs clearly show that the mass
evolution processes for low-resolution simulations (N105)
differ from those of high-resolution simulations (N 106). The

Qº

vR W
,
pG S

(6 )

where Ω=2π/TK is the Keplerian angular frequency(Toomre 1964). A linear stability analysis of a self-gravitating
disk shows that the disk is stable against axisymmetricperturbation when Q>1. The ﬁgures show that the evolutionof each parameter is consistent with the results of Kokubo et al.
(2000), regardless of the numerical resolution. Therefore, thedisk evolution process is essentially independent of numericalresolution.
3

The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

Figure 1. Snapshots of circumterrestrial disk projected onto the x–y plane at t=0, 1, 5, 10, 20, 40 TK. The solid blue circle centered on the coordinate origin
represents Earth, and the green circle is the Roche limit radius. The number of particles is 104.

4

The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

Figure 2. Same snapshots as in Figure 1, but the number of particles is 105.

5

The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

Figure 3. Same snapshots as in Figure 1, but the number of particles is 106.

6

The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

Figure 4. Same snapshots as in Figure 1, but the number of particles is 107.

7

The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

Figure 5. Same snapshots as in Figure 3, but at t=60, 80, 100, 110 TK.

When Q has its minimum value (Q;1), the circumterrestrial disk becomes marginally gravitationally unstable, and thecritical wavelength expected from the linear stability analysis(see, e.g., Toomre 1964) is
lc º

2p 2G SM
~ 2pR disk .
MÅW2of the disk. (Please see Michikoshi & Kokubo 2018 for details.)
Figure 13 shows that the number of spiral arms differs between
low-resolution simulations (N105) and high-resolution
simulations (N106). In the former case, the spiral arms
number are around 20, which is consistent with the above
analytical estimate and with the numerical results of Kokubo
et al. (2000). Conversely, in the latter case, the spiral arms
always number less than 10, and no spiral arms are recognized
for t>80TK and N=107.

(7 )Given the moderate pitch angle of the spiral arms, the numberof spiral arms may be estimated by
ns 

2pR
MÅ
~
.
lc
Mdisk

(8 )

4. Summary and Discussions
We executed N-body simulations of varying numericalresolution (N=104–107) of lunar accretion from a circumterrestrial disk. Irrespective of the numerical resolution, the proto-Moon
forms just outside the Roche limit radius. We ﬁnd that the spiralstructure inside the Roche limit radius differs between lowresolution simulations (N105) and high-resolution simulations(N 106). According to the difference, the angular momentum
ﬂuxes determined by the spiral arms and those determined by theaccretion timescales of the Moon also depend on the numericalresolution.

Our simulation (Mdisk=4ML) corresponds to ns∼20. Considering that the Moon is primarily formed by the materialtransferred beyond the Roche limit radius, the timescale ofMoon formation is almost equivalent to the timescale of themass or angular momentum transfer in the disk due to thegravitational torque exerted by the spiral arms.
We estimate the number of spiral arms for each simulation
using the mode m that has the maximum amplitude of theFourier coefﬁcient of the azimuthally averaged surface density
8

The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

Figure 6. Same snapshots as in Figure 5, but the number of particles is 107.

that angular momentum transfer becomes inefﬁcient when the
wake-like structures appear. Therefore, we cannot conclude
that the complicated structures that appear in the highresolution simulations are self-gravity wakes.
Another candidate for this complicated structure is acoexistence of nonaxisymmetric spiral arms and axisymmetricstructures. The axisymmetric structure is caused by overstableoscillations (Schmit & Tscharnuter 1995). The viscous overstability appears in very-high-resolution simulations and can be agood candidate to explain the small structures in the planetary
rings (Daisaka et al. 2001). Such overstable oscillations would
play an important role in the angular momentum transport in thecircumterrestrial disk and signiﬁcantly modify the angularmomentum ﬂux. However, although Salo (1992) advanced thestudy of the self-gravity and viscous overstability mechanisms
through numerical N-body simulations, no analytical theory cancurrently model the formation of overstability with self-gravity
spirals or wakes (Ballouz et al. 2017). Therefore, a comprehensive explanation of the interplay between the mechanisms ofself-gravity and viscous overstability remains elusive.
In this paper, we focus on evaluating how the lunar accretion
process depends on the number N of particles in an N-body

For high-resolution simulations, the spiral arms are connected to each other, forming more complicated structures
between the spiral arms. The complicated structures resemble
the gravitational wake structures that appear in the numericalsimulations for a dense planetary ring (Daisaka et al. 2001;
Michikoshi & Kokubo 2011). The self-gravity wakes are small
nonaxisymmetric spiral structures formed by gravitational
instability. The wake structure is time dependent and transient,being created and destroyed on a timescale of the order of a
Keplerian period. Daisaka & Ida (1999) showed that as thewake structure grows, the radial velocity dispersion increases
with a large magnitude of ﬂuctuation and tends to maintain therelationship Q∼2. Daisaka et al. (2001) showed that when thewake structure is developed by gravitational instability, angularmomentum transport is signiﬁcantly enhanced not only bygravitational torque due to the wake-like structure but also by
the associated coherent motion of the wakes. Takeda & Ida(2001) also showed that the angular momentum transferred bynonaxisymmetric clumps is always as large as that due to thegravitational torque exerted by the spirals. Therefore, the selfgravity wakes cause efﬁcient radial diffusion in the circumplanetary disk. However, our high-resolution simulations indicate
9The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

Figure 7. Mass evolution for the largest aggregate (lines) and for the sum of the three largest aggregates (dotted lines) normalized by lunar mass ML for N=104, 105, 106,
and 107. Each panel, except for N=107, shows two sets of results with the same initial conditions but different random seeds to generate the initial particle positions andvelocities.

Figure 8. Accumulative masses accreted onto the Earth for N=104, 105, 106, and 107. Each panel, except for N=107, shows two sets of results with the same initialconditions but different random seeds to generate the initial particle positions and velocities.

10

The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

Figure 9. Plot of the data in Table 1 to see the convergent behavior with numerical resolution.

Figure 10. Surface density Σ plotted as a function of distance from Earth for N=104, 105, 106, and 107 at t=0, 1, 5, 10, 20, 40 TK.

11The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

Figure 11. Radial proﬁle of rms velocity dispersion in the radial direction of disk particles for N=104, 105, 106, and 107 at t=0, 1, 5, 10, 20, 40 TK(the legend is
the same as for Figure 10).
Table 1
Final Mass of the Proto-Moon and the Mass Fallen to the Earth
N
4

10
104
105
105
106
106
107

Table 2
Final Mass of the Proto-Moon Found by Kokubo et al. (2000)

M[ML]

M3[ML]

Mfall[ML]Run

Mf=0.7[ML]Mf=1.0[ML]

0.556
0.409
0.390
0.420
0.374
0.352
0.383

0.562
0.481
0.462
0.430
0.644
0.394
0.490

2.980
3.186
2.786
2.742
2.457
2.554
2.636

13
25
26
27
30

0.925
0.875
0.957
0.640
0.895

0.834
0.676
0.965
0.872
0.860

Note. The mass of the largest aggregate for the partial ( f=0.7) and total( f=1.0) accretion models in Kokubo et al. (2000). We only show limited
results from runs in which the initial conditions are the same as those in our
simulations. The masses are normalized by the present lunar mass ML.

Note. The mass of the largest aggregate M, the mass of the sum of the threelargest aggregates M3, and the mass fallen to the Earth Mfall at t=110TK forN=104, 105, 106, 107. The masses are normalized by the present lunar massML=0.0123 M⊕. Each row except for N=107 shows two sets of the resultswith the same initial conditions but different random seeds to generate theinitial particle positions and velocities.

valuable discussions. We also thank the people at Pezy
Computing for their help in coding on Pezy-SC. We used thesuper computer Shoubu at RIKEN. This work was supported
by JSPS KAKENHI grant number JP15K17750 and MEXT
grant number 26106006.

simulation and ﬁnd qualitative and quantitative differences in
the Moon formation process. In future works, we plan to
execute a large number of high-resolution N-body simulationsof Moon formation with a variety of disk parameters toinvestigate the physics that would cause such differences.

Appendix
Pezy-SC Processors
The Pezy-SC processor is a novel new architecture that
integrates 1024-core MIMD processors with a physicallyshared memory and hierarchical cache. It was developed by

We acknowledge the reviewer’s helpful comments. We
thank Shugo Michikoshi for providing an analysis code and
12

The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

Figure 12. Toomre’s Q value as a function of distance from Earth for N=104, 105, 106, and 107 at t=0, 1, 5, 10, 20, 40 TK(the legend is the same as for Figure 10).

Figure 13. Evolution of the number of spiral arms for N=104, 105, 106, and 107. The analysis code was provided by S. Michikoshi.

13

The Astrophysical Journal, 856:175 (14pp), 2018 April 1

Sasaki & Hosono

Figure 14. Measured wall-clock time of our code on a CPU and PEZY-SC for N=104, 105, 106, and 107 with a line of N log8 N.

the Japanese venture company Pezy Computing. The Pezy-SC
processor has achieved impressive performance per watt: the
ExaScaler system based on the ﬁrst-generation Pezy-SC ranked
second in the Green500 list in 2014 November, and three
ExaScaler systems occupied the top three ranks in theGreen500 list in 2015 June. Please see Yoshifuji et al. (2016)
for details on the Pezy-SC processor.
In order to show the speed-up efﬁciency of the PEZY-SC
device, we brieﬂy show the measured wall-clock times of our codeby changing the number of particles deployed in a run in
Figure 14. For a comparison, we also show the wall-clock timemeasured on a CPU (Intel Xeon E5-2618L v3). The wall-clock
time of our code becomes roughly N log8 N. The speed-up
efﬁciency of PEZY-SC compared to CPU reaches to about
20 times faster. In the case of N=107, the wall-clock time with
the number of processes of 32 is about 9 s. In our simulation, theend time is set to 3500 and the timestep is 2−9. In this case, if thenumber of particles is constant, the estimated wall-clock time toﬁnish the run is 3500×512×9 s which corresponds to about
half a year. Recall that the number of particles in one run will
decrease to a third of the initial number of particles. Thus, theactual wall-clock time is less than the estimated value. In this case,the actual value is about 3 months. Note that our code surely has
room for further speed-ups. Since the main purpose of this paper isthe Moon formation, detailed measurement and breakdowns (e.g.,
wall-clock time on the host vs. the device and strong/weak scaling
to the number of devices) should be given in another context.

ORCID iDs
Takanori Sasaki

https://orcid.org/0000-0003-1242-7290
References

Ballouz, R.-L., Richardson, D. C., & Morishima, R. 2017, AJ, 153, 146
Barnes, J., & Hut, P. 1986, Natur, 324, 446
Benz, W., Slattery, W. L., & Cameron, A. G. W. 1986, Icar, 66, 515
Cameron, A. G. W., & Ward, W. R. 1976, LPI, 7, 120
Canup, R. M. 2004, Icar, 168, 433
Canup, R. M., & Asphaug, E. 2001, Natur, 412, 708
Daisaka, H., & Ida, S. 1999, EP&S, 51, 1195
Daisaka, H., Tanaka, H., & Ida, S. 2001, Icar, 154, 296
Hartmann, W. K., & Davis, D. R. 1975, Icar, 24, 504
Ida, S., Canup, R. M., & Stewart, G. R. 1997, Natur, 389, 353
Iwasawa, M., Tanikawa, A., Hosono, N., et al. 2016, PASJ, 68, 54
Kokubo, E., Ida, S., & Makino, J. 2000, Icar, 148, 419
Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141
Michikoshi, S., & Kokubo, E. 2011, ApJL, 732, L23
Michikoshi, S., & Kokubo, E. 2018, ApJ, submitted
Richardson, D. C. 1994, MNRAS, 269, 493
Salmon, J., & Canup, R. M. 2012, ApJ, 760, 83
Salo, H. 1992, Natur, 359, 619
Schmit, U., & Tscharnuter, W. M. 1995, Icar, 115, 304
Takeda, T., & Ida, S. 2001, ApJ, 560, 514
Thompson, C., & Stevenson, D. J. 1988, ApJ, 333, 452
Toomre, A. 1964, ApJ, 139, 1217
Ward, W. R. 2012, ApJ, 744, 140
Yoshifuji, N., Sakamoto, R., Nitadori, K., & Makino, J. 2016, arXiv:1612.
00530

14