The Giant Impact Simulations with Density Independent Smoothed Particle Hydrodynamics

At present, the giant impact (GI) is the most widely accepted model for the origin of the Moon. Most of the numerical simulations of GI have been carried out with the smoothed particle hydrodynamics (SPH) method. Recently, however, it has been pointed out that standard formulation of SPH (SSPH) has difficulties in the treatment of a contact discontinuity such as a core-mantle boundary and a free surface such as a planetary surface. This difficulty comes from the assumption of differentiability of density in SSPH. We have developed an alternative formulation of SPH, density independent SPH (DISPH), which is based on differentiability of pressure instead of density to solve the problem of a contact discontinuity. In this paper, we report the results of the GI simulations with DISPH and compare them with those obtained with SSPH. We found that the disk properties, such as mass and angular momentum produced by DISPH is different from that of SSPH. In general, the disks formed by DISPH are more compact: while formation of a smaller mass moon for low-oblique impacts is expected with DISPH, inhibition of ejection would promote formation of a larger mass moon for high-oblique impacts. Since only the improvement of core-mantle boundary significantly affects the properties of circumplanetary disks generated by GI and DISPH has not been significantly improved from SSPH for a free surface, we should be very careful when some conclusions are drawn from the numerical simulations for GI. And it is necessary to develop the numerical hydrodynamical scheme for GI that can properly treat the free surface as well as the contact discontinuity.


Introduction
The origin of the Moon is one of the most important problems in the planetary science.
The giant impact (GI) hypothesis (Hartmann and Davis 1975;Cameron and Ward 1976) is currently the most popular, since it can solve difficulties that other models face, such as the current angular momentum of Earth-Moon system and Moon's small core fraction compared to the other rocky planets. According to the GI hypothesis, at the late stage of the terrestrial planet formation, a Mars-sized protoplanet collided with the proto-Earth and produced a circumplanetary debris disk, from which the Earth's Moon is formed. To examine whether this scenario really works or not, a number of numerical simulations of collisions between two planetary embryos have been carried out (e.g., Benz et al. 1986Benz et al. , 1987Benz et al. , 1989Cameron and Benz 1991;Cameron 1997;Canup and Asphaug 2001;Canup 2004;Nakajima and Stevenson 2014). Most of them were done by using the smoothed particle hydrodynamics (SPH) method which was a widely used particle-based fluid simulation method developed by Lucy (1977) and Gingold and Monaghan (1977). Recently, however, it is pointed out that the results of numerical simulation of GI by SPH method should be re-examined from the geochemical point of view.
Recent high precision measurement of isotope ratio revealed that it is not easy for the GI hypothesis to reproduce observed properties of the Moon. The Moon and the Earth have almost identical isotopic composition for oxygen (Wiechert et al. 2001), and isotopic ratios chromium (Lugmair and Shukolyukov 1998), titanium (Zhang et al. 2012), tungsten (Touboul et al. 2007) and silicon (Georg et al. 2007). This means that the bulk of the Moon should come from the proto-Earth, unless very efficient mixing occurred for all the isotopic elements (Pahlevan and Stevenson 2007). On the other hand, in previous numerical simulations of GI, the disk material comes primarily from the impactor, which is likely to have had the different isotopic compositions from that of the Earth. To solve this problem, several models have been proposed and studied numerically. These models have the total angular momentum significantly larger than that of the present Earth-Moon system. Models with a fast rotating proto-Earth (Ćuk and Stewart 2012), a hit-and-run collision (Reufer et al. 2012) and a massive impactor (Canup 2012) have been proposed. Although the excess angular momentum is assumed to be removed by the evection resonance with the Sun (e.g.,Ćuk and Stewart 2012), it may work only a narrow range of tidal parameters (Wisdom and Tian 2015). This means that the Moon was formed by a fortuitous event.
Recently, however, it is pointed out that the results of numerical simulations with the standard formulation of SPH (SSPH) is problematic. It turned out that SSPH has problems in dealing with a contact discontinuity and a free surface. It is pointed out that these difficulties result in serious problems, such as the treatment of hydrodynamical instabilities (e.g., Okamoto et al. 2003;Agertz et al. 2007;Valcke et al. 2010;McNally et al. 2012). This problem arises from the assumption in SSPH that the local density distribution is differentiable, though in real fluid, the density is not differentiable around the contact discontinuity. As a result, around the contact discontinuity, the density of the low-density side is overestimated and that of the high-density side is underestimated. Thus, pressure is also misestimated around the contact discontinuity and an "unphysical" repulsive force appears. This unphysical repulsive force causes a strong surface tension which suppresses the growth of hydrodynamical instabilities. In the GI simulation, since the core-mantle boundary is a contact discontinuity and the planetary surface also has a density jump, the accurate treatment of the contact discontinuities is very important.
We have developed a novel SPH formulation, density independent SPH (DISPH), to solve the problem for the contact discontinuity Hosono et al. 2013;Hopkins 2013). Instead of the differentiability of the density, DISPH requires the differentiability of the pressure. As a result, DISPH significantly improves the treatment of the contact discontinuity. Thereby, GI must also be re-investigated by DISPH. In this paper, we present results of GI simulations performed with DISPH and compare them with those obtained with SSPH by focusing on the treatment of contact discontinuity and its impacts on the results. We found that DISPH produced significantly different debris disk, which should lead to different moon forming process. We concluded that the results of GI is sensitive to the numerical scheme and previous numerical simulations of GI should be re-considered.
It is worth noting that Wada et al. (2006) reported the results of GI by a threedimensional grid-base method. They found that the post impact evolution of the disk is different from that of SSPH. They pointed out that it is due to the poor description of debris disk. They suggested that the difference may be due to low-resolution for the debris disk in SPH calculations. However, it could be rather due to their oversimplified polytropic-like EOS. Thus, it is not straightforward to compare their results with SSPH. Canup et al. (2013) reported the comparison of the results of GI between adaptive mesh refinement (AMR) and SPH and concluded that the predicted moon mass of two methods are quantitatively quite similar. Although we notice qualitative differences in disk spatial structures in some of these results (for example, the different clump structure between AMR and SPH in Fig. 4 of Canup et al. (2013)), our DISPH also predicts similar moon masses for the collision parameters that they tested, as we will mention in section 5.
Comprehensive code-code comparison is needed with grid codes, as well as between DISPH and SSPH. In this paper, we focus on the latter comparison.
Here we do not insist that the results of GI simulations by DISPH are much closer to realistic phenomenon than by SSPH. While DISPH has been improved for treatment of a contact boundary, both DISPH and SSPH have a problem to treat free surface, i.e., planetary surface. We here stress that only the improvement for treatment of a contact boundary significantly affects properties of circumplanetary disks generated by GI.
Therefore, we need to be very careful when some definitive conclusions are drawn from the current numerical simulations for GI. To clarify details of Moon formation, it is necessary to develop the numerical hydrodynamical scheme for GI that properly treats the planetary surface as well as the core-mantle boundary. This paper is organized as follows. In section 2, we briefly describe the numerical technique. We focus on the implementation of DISPH for non-ideal EOS. In section 3, we describe models of the GI simulations. In section 4, we show the results and comparisons of the GI simulations with the two methods and clarify the reason for the difference in the properties of the generated disks. We also show the results of single component objects, in addition to those of differentiated objects with core-mantle structure. The former and latter simulations discriminate the differences due to a free surface from a core-mantle boundary between the two methods. In section 5, we summarize this paper.

Overview of DISPH
In the SPH method, the evolution of fluid is expressed by the motions of fluid elements that are called SPH particles. The governing equations are written in the Lagrangian form of hydrodynamic conservation laws. The equations of motion and energy of the i-th particles are written as follows: where r, a, u and t are the position vector, the acceleration vector, the specific internal energy and the time, respectively. The subscript i denotes the value of i-th particle. The superscripts hydro, visc and grav mean the contributions of the hydrodynamical force evaluated by SPH, viscosity and self-gravity, respectively. The formal difference between SSPH and DISPH is in the form of a hydro i and (du i /dt) hydro . The other terms in the right hand side of Eqs. (1) and (2) have the same forms for both methods. Since the equations of SSPH can be found in the previous literature (e.g., Benz et al. 1986;Monaghan 1992;Canup 2004), here we will only describe those of DISPH.
DISPH is originally developed by Saitoh and Makino (2013) for the ideal gas EOS and then extended to an arbitrary EOS by Hosono et al. (2013). The main advantage of DISPH is the elimination of unphysical surface tension which rises at the contact discontinuity.
The unphysical surface tension in SSPH comes from the requirement of the differentiability of the density. Saitoh and Makino (2013) developed a new SPH formulation which does not require the differentiability of the density, but requires that of (a function of) pressure. As a result, their new SPH can correctly handle hydrodynamical instability.
Note that around the shock, neither the pressure nor the density is continuous. Thus the assumption of the differentiability of pressure and density is broken across the shock. Saitoh and Makino (2013) used the quantity p α where p is as pressure and α is a constant exponent for their formulation to improve the treatment of the shock 1 . Saitoh and Makino (2013) applied α = 0.1 and show that the treatment the strong shock is improved. However, they considered only ideal gas EOS. Here, we apply this formulation to fluids with non-ideal EOS. The choice of α is related with how to treat a free surface such as a planetary surface. DISPH overestimates the pressure gradient around the free surface, while SSPH underestimates it.
In Appendix A, we show the results of the 3D shock problem simulated with SSPH and DISPH. In the ideal gas EOS case (Fig. 17), the numerical pressure blip at the contact boundary is the smallest for α = 0.1, while the numerical density blip is relatively large for the same α = 0.1. Saitoh and Makino (2013) determined α = 0.1 as a best choice for ideal gas through more detailed discussions. On the other hand, in the Tillotson EOS case, it is not clear that which α is the best choice. In the case of 3D shock problem (Fig. 17), α = 1 is the best choice. However, in the case of strong shock test, α = 0.1 may be the best choice. Tillotson EOS (Fig. 18) shows the different dependence of the numerical pressure blip at the contact boundary on α from that with ideal EOS. This indicates that the smaller α works better to treat large pressure jump, namely, the strong shock and the free surface.
We adopt α = 0.1 also in this paper, although more comprehensive tests on the choice of α are needed in future study. We will show the effects of improved treatment of the free surface and the core-mantle boundary in section 4.1 and 4.2, respectively.
In Appendix B, we also show the results of the Keplerien disk test and Gresho vortex problem. Results of the Keplerien disk tests with SSPH and DISPH tell us that both schemes can maintain the disk structure during the first several orbital period, whereas catastrophic breakup takes place before 10 orbital period. Previous studies also reported the same results (Cullen and Dehnen 2010;Hopkins 2015 The essential difference between DISPH and SSPH is in the way to estimate the volume element of a particle, ∆V i . As the starting point, following Saitoh and Makino (2013) and Hosono et al. (2013), we introduce the physical quantity: where brackets mean "smoothed" values. Hereafter, we denotes p α as y. The value of y i is given as follows: where W and h i are the kernel function (see below) and the so-called smoothing length, respectively. By using y i and Y i , we derived the equations of motion and energy for our scheme as follows: where m i and v i are the mass and the velocity vector of particle i, respectively. Here Ω is the so-called "grad-h" term (e.g., Springel and Hernquist 2002;Hopkins 2013;Hosono et al. 2013); Here,ŷ is the value to determine the smoothing length; Note that the choice ofŶ and ŷ is arbitrary, as far asŶ /ŷ has the dimension of volume.
In this paper, following Saitoh and Makino (2013), we choseŶ = m and ŷ = ρ , where ρ is the density. Note that since the interactions between two particles are antisymmetric, our SPH conserves the total momentum and energy. The grad-h term improves the treatment of the strong shock (e.g., Saitoh and Makino 2013;Hopkins 2013). In Appendix C, we show the results of strong shock with DISPH both with and without the grad-h term. We show that our DISPH with the grad-h term works well for the strong shock. Our DISPH with the grad-h term has enough capability for the problems which include strong shock.
Note that in order to actually perform numerical integration, we need to determine new values of y i , by solving a set of implicit equations, Eq. (4) combined with the equation of state p = p(ρ, u). Thus, as in Hosono et al. (2013), we iteratively solve Eq. (4). The number of iterations is set to 3, following the previous works (e.g., Sec. 5.6 in Saitoh and Makino 2013). The iteration procedure is the same as that described in Hosono et al. (2013), except for the initial guess of Y i . The initial guess of Y i is obtained by the numerical integration of Y i using its time derivative: Note that these equations reduce to those of Hosono et al. (2013) in the case of α = 1.
For the kernel function W , we employ the cubic spline function (Monaghan and Lattanzio 1985). Note that the use of the cubic spline kernel for the derivative sometimes causes the paring instability (Dehnen and Aly 2012;Price 2012). In order to avoid this instability, we adopt a gradient of the kernel which has a triangular shape (Thomas and Couchman 1992).
For the artificial viscosity, with both methods, we use the artificial viscosity described in Monaghan (1997). Note that for both methods we use the smoothed density for the evaluation of artificial viscosity. The parameter for the strength of the artificial viscosity is set to be 1.0. In order to suppress the shear viscosity, we apply the Balsara switch (Balsara 1995) to the evaluation of the artificial viscosity.
The self-gravity is calculated using the standard BH-tree method (Barnes and Hut 1986;Hernquist and Katz 1989). The multipole expansion is calculated up to the quadrupole order and the multipole acceptance criterion is the same as Barnes and Hut (1986). The opening angle is set to be 0.75.

Timestep
In the SPH method, the timestep is usually determined by the Courant condition as follows (Monaghan 1997): where and c i is the sound speed of particle i. Here C CFL is a CFL coefficient which is set to 0.3 in this paper. In this paper, we consider three additional criteria: fractional changes in the specific internal energy, Y (DISPH only), and the accelerations. These are where C u , C Y and C a are dimensionless timestep multipliers. Throughout this paper, we set C u = C Y = C a = C CFL . Note that Eqs. (14) and (15) are applied when du i /dt < 0.
Throughout this paper, we use the Tillotson EOS (Tillotson 1962), which is widely used in the GI simulations. The Tillotson EOS contains 10 parameters, which we should choose to describe the given material. The material parameters of the Tillotson EOS for each material are listed in Melosh (1989), page 234, table AII.3. Note that in the very low density regime, the Tillotson EOS gives negative pressure which is unphysical on the scale of GI. To avoid numerical instabilities due to negative pressure, we introduce a minimum pressure p min for the Tillotson EOS. In the scale of GI, the typical value of pressure is order of ∼ 100 GPa. Throughout this paper, thus, we set p min = 0.1 GPa. Also, we impose the minimum timestep to prevent the timestep from becoming too small due to unphysical values of partial derivatives of EOS by density. We carefully determine the minimum timestep as one second not to cause poor description of the physical evolution of a system in this paper. In addition, in this case, we do not evaluate hydrodynamical terms in Eqs.
(1) and (2), since this small timestep is actually applied for particles with very low density.

Initial condition
We performed numerical simulations of GI from eight initial models. In this section, we briefly describe how we set up the initial conditions.
We first constructed two initial objects, the proto-Earth (target) and the impactor, which satisfy the given impactor-to-target mass ratio and total mass. We use ∼ 3 × 10 5 SPH particles in total. Following Benz et al. (1987), both objects consist of pure iron cores and granite mantles. First, we place equal-mass SPH particles in a Cartesian 3D-lattice.
Then, inner 30% of the object is set up as iron and the remaining outer part is set up as granite. The initial specific internal energy of particles is set to 0.1GM E /R E J/kg and the initial velocity of particles is set to zero. Here, G and M E are the gravitational constant and the current Earth's mass. We let the SPH particles relax to the hydrostatic equilibrium by introducing the damping term (Monaghan 1994) to the equation of motion. The end time of this relaxation process is set to t = 10000 seconds, which is about ten times of the dynamical time for the target. After this relaxation process, the particle velocities for each particle are ∼ 1% of the typical impact velocity (an order of 10 km/sec in the case of the Moon forming impact).
We constructed eight models. One of them, model 1.10, corresponds to run #14 of Canup and Asphaug (2001). They concluded that the Moon would form in this run. In this model, the impactor approaches the proto-Earth in a parabolic orbit. Another model, model 1.17, was close to run #7 of Benz et al. (1987). In this model, the initial relative orbit is hyperbolic. The remaining models have the same parameters as those of model 1.10 except for the initial angular momentum. High and low angular momentum models correspond to high-oblique and low-oblique collisions. In all models, initial objects are non-rotating. We integrated the evolution of these models for about 1 day. This duration time of the simulation is smaller than the time scale of numerical angular momentum transfer due to the artificial viscosity (for detail, see Canup 2004). When we present the result, we set the time of the first contact of two objects as the time zero. Table 1 shows the summary of the initial conditions. The columns indicate the initial separation between two objects (R init /R E ), initial angular momentum of the impactor (L init ), velocity at the infinity (v ∞ ), the total number of particles (N tot ), the number of particles of target (N tar ), the number of particles of impactor (N imp ), the mass of target (M tar /M E ), the mass of impactor (M imp /M E ) and the mass of core (M core /M tot ). Here, L EM is the angular momentum of the current Earth-Moon system, respectively. We set

Results
We first show the results of collisions of objects with a single component without core-mantle boundary structure in section 4.1. This set of runs is to discriminate the effects of a free surface from those of core-mantle structure and the core-mantle boundary. Then we show the results of collisions of two differentiated objects with core-mantle boundary. We overview the time evolution of eight models obtained with two different methods, DISPH and SSPH in section 4.2. In all runs, we found the differences between the results with two methods are rather significant. In section 4.3, we compare the predicted mass of the moon obtained with two methods. In section 4.4, we investigate the cause of this difference.

Collisions of single-component objects
We consider collisions between single-component planets consisting of only granite mantle. Here we performed two types of impacts; one is the collision between equal mass objects, the other is the same target-to-impactor mass ratio as described in section 3, but with single-component objects. Since the objects have no core-mantle boundary, the difference between the results of two methods should come from the treatment of the free surface.
The initial objects are constructed in a similar way to the prescription in section 3, although there is no iron core in this case. In the run with equal-mass objects, both objects have mass of 1M E and radius of 1R E , and the initial specific internal energy is set to be 0.1GM E /R E . The initial angular momentum is the same as the current angular momentum of the Earth-Moon system. The velocity of the impactor at infinity is zero. In this simulation we employ 300,754 particles in total. Figure 1 shows the radial profiles of mass and density of the final outcome of the collision of two equal mass objects for both methods. With SSPH, a gap in the particle distribution is formed around r ≃ 2.5R E , while with DISPH the radial distribution is continuous. This means that SSPH produces gap structure between the body and disk and more spreading disk than DISPH. The gap at r ≃ 2.5R E is also found in the snapshot on the x-y plane with SSPH ( Fig. 2). There is, however, no physical reason for the formation of this gap. It seems to be natural that the angular momentum distribution is continuous.
Why the gap is formed in the SSPH simulation is most likely the same as the gap formation at the contact discontinuity (for detail, see section 4.4). Since there is a density jump around the free surface, the free surface is a kind of contact discontinuity. Though there is no discontinuity in the density distribution, the slope is steep at 2-3 Earth radius and the density itself is low. Thus, the density difference between two particles radially separated can be very large, resulting in the problem similar to that in the contact discontinuity.
DISPH does not suffer from such a problem.
The result in Appendix D also suggests that DISPH is better than SSPH for the treatment of the free surface. However, since around the free surface the pressure is not continuous as well as density, it cannot be readily concluded that DISPH is sufficiently improved from SSPH for treatment of the free surface. Figure 3 shows the snapshots for a set of runs of single-component objects with the mass ratio of 10:1 given by Table 1. They show a similar trend on difference between the two methods to that found in Fig. 1. DISPH generally tends to produce more compact disks than SSPH does. Note that the sizes of the planet after an impact are roughly the same between SSPH and DISPH. Because DISPH produces more compact disks, the results by DISPH look as if the planet itself is inflated. The upper two panels show binned mass in log scale. The binned width is set to 0.01R E .
The mass is normalized in the current Earth mass and the distance is normalized by current Earth radius. The middle two panels show mean mass within the binned volume. The lower two panels show radial distance from central planet vs. density for each particle. Each axis is shown in log scale. The left column shows the results for DISPH and the right column shows those of SSPH. The left column shows the results of DISPH and right column shows those of SSPH. The upper row is shown in x − y plane and the lower row is shown in y − z plane. The color of each particle indicates the original objects of particle.
In order to compare the results between two methods quantitively, we employ the so-called "predicted moon mass" as just a reference value. First, we extract "disk particles" from the simulation results. Following Canup and Asphaug (2001), we classified SPH particles to three categories, namely, escaping particles, disk particles and planet particles.
Particles whose total (potential + kinetic) energies are positive are regarded as escaping. If the total energy of a particle is negative and its angular momentum is greater than that of the circular orbit at the surface of the planet, it is categorized as a disk particle. Then, other particles are categorized as planet particles, since these particles should fall back to the target. After the particles are classified, we predict the moon mass using information of disk particles. According to the N-body simulations by Ida et al. (1997) and Kokubo et al. (2000), the predicted moon mass, M M , is given by: where L disk , R Roche and M disk,escape are the angular momentum of the disk, the Roche radius of the planet, the mass of the disk, respectively. Note that M disk,escape is the total mass of disk particles that escape from the disk through scattering by accreting bodies. Following previous works (e.g., Kokubo et al. 2000;Canup 2004), we set M disk,escape to 0.05M disk .
Assuming that materials from the proto-Earth and the impactor are well mixed in the disk particles, we estimate the fraction of the moon materials originating from the proto-Earth.  However, DISPH produces more compact disks and accordingly smaller M M than SSPH does. Since the two objects have no core-mantle boundary, this difference should come from the treatment of the free surface. Figure 5 shows the distributions of specific internal energy for model 1.15 with both methods. The difference between two methods are clear. The first two snapshots for each method look fairly similar; shock heating and the arm-like structure can be clearly shown.
In the panels t > 2.3 hrs, however, clear difference between two methods can be seen. With DISPH, the arm re-collides to the proto-Earth and undergoes shock heating again, which results in hot and compact debris disk (panel t = 24.0 hrs). On the other hand, with SSPH, cold particles are ejected around the arm-like structure (panels t = 2.3 − 3.3 hrs). These ejected particles finally become the cold and expanded disk (panels t ≥ 4.7 hrs). Note that similar cold and expanded disk can be seen in previous studies with SSPH. This difference might come from the treatment of free surface or shock. Since in this paper we focused on the treatment of core-mantle boundary, further investigation of the origin of this difference is left for future works.
We will show a similar plot for collisions between differentiated objects. As we will show in Fig. 13, the dependance of M M on L init is different from that in Fig. 4. The difference is due to the contribution from core-mantle structure and its boundary. The core-mantle boundary has to be treated properly as well as the free surface in GI simulations.

Collision of objects with core-mantle structure
Next we discuss details of collisions between differentiated objects with core-mantle structure. We will show that the contact density discontinuity at the core-mantle boundary may be a problem in calculations with SSPH and it is improved with DISPH. Figure 6 shows the time evolution of model 1.10 obtained with two methods. What is shown here is the face-on view (in the x-y plane) of particles with negative values of the z coordinate, seen from z = ∞.
The two runs in the first frame (t = 0.5 hr) look fairly similar. In both runs, the core of the impactor (black) is significantly deformed. It pushes up the mantle of the proto-Earth (orange) and the mantle of the impactor (red) is left behind. In the first four frames, the results of the two models are qualitatively similar. In both runs, the core of the impactor forms an arc-like structure at t = 1 hr, which becomes more extended at t = 1.5 hr. However, this arc is more extended for the run with SSPH. In the frames of t = 2 hr through t = 7 hr, most of the mass of the impactor, which has not escaped from the planet gravitational potential, fall back to the proto-Earth in the run with DISPH, while some of the mass forms extended envelope and disk in the run with SSPH. As we will see in the next section, this difference in the structure causes a large difference in the predicted moon mass. Figure 7 shows the edge-on views of the two runs of model 1.10. The vertical distribution of ejected mantle material is also quite different. In the frames of t = 1 hr and 2 hr, the results are qualitatively similar. In the frames of t = 3 hr through 7 hr, however, SSPH produces vertically stretched disk. On the other hand, with DISPH, the number of the disk particles is much smaller than with SSPH. DISPH produces a thinner disk than SSPH. Figure 8 shows the distribution of the angular momentum. For t 2 hr, the ejecta with high angular momentum (> √ GM E R E ) is much abundant in the SSPH result (lower panels) than in the DISPH result. Figures 9 and 10 show the same plots as Fig. 6 for the models with low-oblique collision and high-oblique collision, respectively. From these figures, it is obvious that more extended debris disks are greatly formed in the runs with SSPH than in those with DISPH. Figure 11 shows the snapshots at the end time of simulations for all runs. The azimuthal distributions of the disks are different between SSPH and DISPH. With SSPH, particles are widely distributed. In particular, in model 1.17, the distribution is nearly azimuthally symmetric. On the other hand, with DISPH, the distribution of particles is clearly asymmetric and no large disk is formed.
The results of our SSPH runs are similar to those in the previous works. In both cases, disks extended to outside of the Roche radius are formed. On the other hand, the disk is very thin in our DISPH runs for these models. We quantify this difference more clearly in the section 4.4 in order to understand the cause of the difference.

Dependence on disk properties on impact angular momentum
In this section, we summarize dependence of disk properties on initial impact angular momentum. The conditions of the successful moon forming impact is defined by 1) the predicted moon mass is comparable to or larger than the present Moon mass, and 2) most of materials of the formed moon comes from the target planet (the proto-Earth).    with DISPH is similar to or rather smaller than that by SSPH, l disk is higher for DISPH, resulting in larger M M with DISPH than that with SSPH. In addition, SSPH produces larger amounts of escaping mass than DISPH does. These result in the larger M M with DISPH than with SSPH. Figure 13a shows the dependance of M M on L init for runs with SSPH and DISPH.
Generally, M M increases with L init for both methods as in Fig. 4 (runs with single-component objects). Notice that the dependence is much more sensitive in the differentiated objects impacts (Fig. 13a) than in single component objects impacts (Fig. 4). For a high-oblique collision, the impact momentum is transferred to ejecta from the outer parts of the impactor and the target. The volume of ejecta may be primarily regulated by a geometrical effect if the collision velocity is fixed. If the volume of ejecta and momentum transferred to the ejecta are the same for a fixed L init between an impact of differentiated objects and that of single component objects, the total ejecta mass from differentiated objects is smaller and its post-impact velocity is higher than that from single component objects. It results in formation of a more spread disk or a hit-and-run collision for the differentiated objects impact at high-oblique impacts.
With DISPH, M M is an order of magnitude smaller than those with SSPH for L init 1.1L EM , the trend of which is also found in Fig. 4. On the other hand, in the case of L init 1.15L EM , DISPH produces larger M M than SSPH. This is because in the runs with SSPH, more materials are ejected during the first contact event than in the DISPH runs ( Fig. 12d), probably by the artificial tension at core-mantle boundaries of the impactor and the target, as we discuss in details in section 4.4. In the panels of relatively high L init (models 1.17, 1.21 and 1.32) in Fig. 11 clearly show that much more materials are scattered away in the runs with SSPH, while more compact clumps remain in the runs with DISPH. As a result, DISPH produces abrupt transition of the predicted moon mass Since M M is sensitive to the distribution of the disk particles, the difference in M M between SSPH and DISPH is more pronounced than that in the distributions of formed disks.
In model 1.21 and model 1.32, the predicted moon mass with DISPH are greater than disk mass (see, Fig. 12a and c). Since Eq. (17)   Because L init in that case may be smaller than L EM , we can adopt higher impact velocity, which may further increase the disk fraction from the target. We will explore the parameters of the impact with M M ∼ M L and L init ∼ L EM that produces a disk mostly from the target (the proto-Earth) in a separate paper.

Effect of the "Unphysical Surface Tension" in SSPH
In this section we discuss the origin of the difference between two methods by focusing on the treatment of the core-mantle boundary. Figure 14 shows the close-up view of SPH particles in model 1.17 at t = 6 min. Clear In Fig. 15, we show the acceleration per particle along the x-direction, y-direction and torque around the z-axis of the impactor's core and mantle particles. The hydrodynamical forces of SSPH and DISPH runs during the impact phase are different. In the first 10 minutes, the accelerations in both directions of SSPH are larger than those of DISPH. This difference results in the gap shown in Fig. 14. From t = 12 -17 minutes, the SSPH result shows larger torque than that of DISPH. Figure 16 illustrates the effect of this surface tension. In the lower left panel (SSPH, t = 3 min), none of the impactor's core particle suffers negative y-directional force, while in the corresponding snapshot with DISPH at t = 3 min (the upper left panel), some particles suffer negative y-directional force. The amplitude of the acceleration of impactor's core particles is much larger for the SSPH run. Thus, the impactor particles gain upward velocity (in the direction of y-axis), while losing the forward velocity (negative direction of x-axis), compared to the DISPH run. It is most likely that this difference is due to the numerical error of SSPH at the contact density discontinuity (the core-mantle boundary) and it results in the difference in the formed disks between the DISPH and SSPH runs.

Summary
The giant impact (GI) is the most accepted model for the origin of the Moon. However, it is now being challenged. The identical isotope ratios between the Earth and the Moon found by recent measurement require a survey of new ranges of impact parameters, because the impact previously referred to as a "successful Moon forming impacts" produces a moon mostly consisting of materials from the impactor rather than those from the proto-Earth.
We have re-investigated GI by newly developed "Density Independent SPH" (DISPH) scheme with Tillotson EOS. Recently it is recognized that the standard SPH scheme (SSPH) has a serious problem in the treatment of contact discontinuities because SSPH assumes differentiability of density. The core-mantle boundary is a contact discontinuity and the planetary surface (free surface) also has a density jump. The errors result in the unphysical surface tension around the contact discontinuity and the free surface. Since In the case of collisions between single-component objects, compared with SSPH, DISPH always produces more compact disks, for which smaller moon masses are predicted.
This is because numerical repulsive force appears around the free surface, in SSPH runs (section 4.1). Note that since the predicted moon mass is sensitively dependent on the distribution of the disk particles, slight difference in the distribution between SSPH and DISPH can result in significant difference in the predicted moon mass.
On the other hand, in the case of collisions between differentiated objects with core-mantle boundary, DISPH predicts more massive moon masses than SSPH does for high-oblique impacts, while it still predicts lower mass moons for low-oblique impacts (section 4.2). The different dependence on the initial impact angular momentum from the single component objects would come from the transfer of impact momentum to the mantle layer with low density and numerical repulsive force at the core-mantle boundary (section 4.2 and 4.3). The overall trend that the predicted moon mass increases with the initial impact angular momentum is common between SSPH and DISPH.
Note that our result is consistent with the conclusion by Canup et al. (2013): SSPH and a grid code, AMR, produce disks that predict similar moon mass. They did a comparison for impacts with parameters similar to model 1.21. As we showed in Fig. 13, for model 1.21, SSPH and DISPH show similar results within 50% of the predicted moon mass. Thereby, DISPH produces consistent results with AMR for this parameter. The comparison with grid codes is necessary for other initial impact angular momentum for which DISPH and SSPH significantly differ from each other. However, clump structure looks different among AMR, SSPH and DISPH, suggesting that the angular momentum distributions are different among them.
What we want to stress in this paper is that properties of circumplanetary disks generated by GI are sensitive to the choice of the numerical scheme. Only the difference in the treatment of a contact discontinuity (core-mantle boundary) between DISPH and SSPH significantly affects the results. Other effects such as the treatment of free surface, shock propagation, heating so on are also likely to change the results. The results of GI also depend on the EOS, the initial thermal structure, density profiles, material strength and numerical resolution and so on.
We need be very careful when some conclusions are drawn from the numerical simulations for GI, because planets consist of solid layers with different compositions but not uniform gas and current numerical schemes have not been developed enough to treat planets. Thus, we need to develop numerical codes suitable for GI between planets, step by step. The next step of DISPH would be handling of free surfaces and shock propagation that currently has a free parameter α. For GI, code-code comparisons are now needed.
Comparison to experiments or other numerical schemes in the case of simple impact problem is also needed to calibrate the code. These are left for future works. A. 3D shock tube problem by DISPH with α To study the effect of choice of a parameter α in the DISPH scheme, we performed two calculation 3D shock tube problems with SSPH and DISPH with varying the parameter α. One calculation is with the ideal gas EOS and another is with Tillotson EOS. The parameter α is taken to be 1.0, 0.5 and 0.1.
The initial condition of the 3D shock tube problem for ideal gas EOS is set as follows: The velocity of each side is set to be 0. We employ a 3D computational domain, 0 ≤ x < 1, 0 ≤ y < 1/8 and 0 ≤ z < 1/8, and periodic boundary conditions are imposed in all directions. We employ ideal gas EOS with specific heat ratio γ = 1.4. The particle separation of high density side is set to 1/512; thus, the total number of particles is 1572864.
The initial condition of the 3D shock tube problem for Tillotson EOS is set as follows: (1.0, 3.5, 4.85) (otherwise). (A2) The parameters of granite are adopted for Tillotson EOS. The density and specific internal energy are normalized in reference density ρ 0 and reference energy E 0 (Tillotson 1962;Melosh 1989). The computational domain and number of particles are the same as those in the calculation with the ideal gas EOS. Figure 17 shows snapshots at t = 0.1 of the 3D shock tube problem for the ideal gas EOS with both SSPH and DISPH. Both methods show similar results, except for the contact discontinuity (at x ∼ 0.55). As expected, DISPH shows fairly smaller pressure blips than SSPH at the contact discontinuity for all values of α, while larger jumps are found in density with DISPH. These results are consistent with those shown in Saitoh and Makino (2013) and Hosono et al. (2013). Figure 18 shows snapshots at t = 0.05 of the 3D shock tube problem for Tillotson EOS with both SSPH and DISPH. Also in this case, DISPH shows smaller pressure blips around the contact discontinuity. For all values of α, DISPH shows better treatment of the contact discontinuity. The dependence of the magnitude of the pressure blip on α is different between the ideal gas EOS and Tillotson EOS. With Tillotson EOS, the pressure blip is higher for smaller value of α while the blip is still smaller than that with SSPH. However, as shown below, the treatment of a free surface is better for smaller α even with DISPH.
Note that this pressure blip can be a serious problem when we treat the contact discontinuity. Saitoh and Makino (2013) and Hosono et al. (2013) showed the consequence of this pressure blip (see Fig. 4 in Hosono et al. 2013). In their hydrostatic equilibrium tests, they put the high-density square in the low-density ambient in the pressure equilibrium. With SSPH, the high density square, which should remain its initial shape, quickly transforms into circular shape. With DISPH, on the other hand, the high density square remains its initial shape. This means that with SSPH, simulations suffer from the unphysical momentum transfer. Saitoh and Makino (2013) and Hosono et al. (2013) also showed the results of Kelvin-Helmholtz instability test, in which the contact discontinuity plays very important role (see Fig. 5 in Hosono et al. 2013). As expected, SSPH shows unphysical surface tension effect, while DISPH clearly eliminates it. Previous simulations should have suffer from this unphysical effect.
The treatment of the region with abrupt change in the pressure that corresponds to a free surface is improved by DISPH, in particular with small value of α. To show this, we show the pressure field around the strong shock region with Tilloston EOS. The initial pressure distribution is set as follows: The density is uniformly set to be ρ 0 . Figure 20 shows the pressure field and the error at the very first step, where the error is defined as follows: -47 -This figure clearly shows that taking the parameter α small improves the treatment of the large pressure jump, such as free surface.
In Fig. 21, we show the results of three runs of GI with α = 1. The results are somewhat different from those with α = 0.1. However, just like the case with α = 0.1, DISPH produces a smaller moon mass in the model 1.10 and a larger moon mass in the model 1.32.
Our modification for DISPH has good capability for both the shock and contact discontinuity, although DISPH includes the free parameter α. The tests of 3D shock tube and a free surface with ideal gas EOS show that α = 0.1 may be the best choice, similar to Saitoh and Makino (2013). However, with Tillotson EOS, different dependence on α can be seen. The results for GI with Tillotson EOS shows slightly different dependence on , though the number of runs is few. Thus, more careful calibration for the dependence on α should be done, which is left for future work. Unless otherwise specified, in the following, we adopt α = 0.1.

B. Tests for the conservation of angular momentum
For the calculation of GI problems, it is important to treat the angular momentum transfer correctly. In both SSPH and DISPH, since the interaction between two particles is pairwise, the global angular momentum is conserved. However, it is often said that SSPH does not treat the local angular momentum transfer correctly; unphysical angular momentum transfer due to the so-called zeroth order error and spurious viscosity appear.
In this appendix, to test whether SSPH and DISPH can treat the angular momentum transfer correctly or not, we performed two well-posed tests; one is the Keplerian disk test   Fig. 13, but shows the results with α = 1.0 (blue squares). Gresho and Chan 1990).
For the Keplerian disk test, we initialize two-dimensional disk whose surface density is set to 1.0. The inner and outer edges of the disk are set to 0.5 and 2.0. The initial pressure of the disk is set to 10 −4 and the heat capacity ratio of the ideal gas is set to 1.4. The self gravity between particles are ignored, while the gravity from the central star a = −r/r 3 acts on each particle. In this test we employ 48228 particles in total.
For the Gresho vortex test, we employ [−1, 1) 2 periodic boundary computational domain with the uniform density of unity. The initial pressure and azimuthal velocity distributions are as follows: where R = x 2 + y 2 . In this test we employ 16, 384 particles in total.  From these tests, we should conclude that both methods can handle the local angular momentum transfer to the same degree. Note that in GI simulations, we set the end time at t = 24 hrs, which corresponds to about 3.4 orbital time at the Roche limit. Figure 22 shows that both methods can treat the local angular momentum transfer until 3.4 orbital time. Overall, both SSPH and DISPH are capable of dealing with rotation disks with similar degree, as far as the simulation time is less than 3.4 orbital periods.

C. Sedov-Taylor blast wave test
Here we show the results of Sedov-Taylor blast wave test, which shows the capability of the strong shock. The initial condition of this test is the same as that used in Saitoh and Makino (2013). We employ 3D computational domain [0, 1) 3 . We first place equal-mass 2, 097, 152 particles in the glass distribution with the uniform density of unity. Then, the explosion energy is added to the central 32 particles. In this test we use the ideal gas EOS with γ = 5/3. Figure 25 shows the results of this test with DISPH with and without grad-h term.
Without grad-h term, clearly, the shock propagates slower than the semi analytic solution.
DISPH with the grad-h, on the other hand, shows better results than DISPH without the grad-h term. Our DISPH with the grad-h can treat the strong shock. Note that these results are consistent to Saitoh and Makino (2013) and Hopkins (2013).

D. Test for the treatment of the free surface
To check the capability for the free surface, we performed three simple 2D tests which include free surface with both DISPH and SSPH. The first test is hydrostatic equilibrium test, which is carried by Monaghan (1994). The second test is the vertical impact of aluminium-to-aluminium test and the third test is the glass-on-water test. The latter two tests are performed by Pierazzo et al. (2008).
For the first test, we employ 2D computational domain [0, 40) × [0, 40) and 4, 096 particles in total. The periodic boundary condition is imposed on the x-direction. We set up the fluid which is initially in pressure equilibrium under the constant gravitational acceleration g = −10 along the y-direction. We fix the positions and internal energies of all particles with y < 4. Around y = H, there is a free surface. In this test we use the following EOS for linear elastic material: The material parameters A and ρ 0 are set to the granite's values in the Tillotson EOS. Since this system is in a hydrostatic equilibrium, particles should maintain their initial positions.
For the second test, we first placed the target particles in [0, 10 4 ) 2 . Then, we placed projectile particles with the impact velocity 20 km/sec. We employ 67, 597 particles in total. We set the impact angle to 0 • (vertical impact). The radius of projectile is set to 10 3 m. In this test we use the Tillotson EOS and the material parameters are set to the aluminium's values, similar to Pierazzo et al. (2008). Note that in this test, we omit the material strength.
In the third test, we followed the early time evolution of the glass-on-water test, following Pierazzo et al. (2008). We first place target particles in [0, 0.05) 2 . Then, we placed projectile particles with the impact velocity 4.64 km/sec. We employed 65, 336 particles in total. We set the impact angle to 0 • (vertical impact) The radius of projectile is set to 1 mm. We used the Tillotson EOS and the material parameters are set to water for target and wet tuff for projectile (Melosh 1989). indistinguishable between SSPH and DISPH. DISPH has similar accuracy/errors for free surface as SSPH does. Note that there are several differences between two results, e.g., the height and expansion of impact jetting. Figure 29 shows the results of the glass-on-water test with both methods. Unlike the aluminium-to-aluminium test, this test contains the contact discontinuity between water and wet tuff. Similar to the aluminium-to-aluminium test, the height and expansion of the ejecta curtain is different between two methods, which could be due to the unphysical surface tension between two different materials arising in SSPH calculations. The target particles are pushed up by the projectile particles at the early step of the impact (t = 0.6µs -2.0µs). This results in the higher crater rim with SSPH than DISPH. At t = 13.9µs, SSPH produces oblate projectile, while with DISPH, the projectile and target are mixed.
This clearly due to the unphysical surface tension term which results in an underestimate of material mixing. This difference may be related to the difference in impact-generated disks in the GI simulations between DISPH and SSPH. Figure 30 shows the cumulative mass of ejecta M ejecta (v) with a vertical velocity greater than a given velocity v. According to the previous works (e.g., Holsapple 1993; Housen and Holsapple 2011), the results should have a power-law form with a slope of 3µ: where v imp , ρ imp and ρ tar are the velocity of the impactor, the density of the impactor and the density of the target. Here, ν and µ are material parameters which are set to be 0.4 and 0.55 (Housen and Holsapple 2011). Both methods reproduced roughly similar results to the experiments shown in Housen and Holsapple (2011). The power-law regime with a slope of −3µ is well reproduced with both methods. However, SSPH produces high speed jetting component (v 0.5), which can hardly be seen in the experimental results. This difference should come from the fact that the target particles feel unphysical surface tension from the penetrating projectile, as stated in the previous paragraph. The target particles are pushed up by the projectile particles to acquire high vertical velocity. This could result in the difference of the results of GI between two methods.
It is not clear which SPH scheme is more correct, especially for the free surface. The tests carried out in this Appendix are performed using a 2D cartesian geometry. The results may differ from 2D cylindrical or 3D geometries. Thus, it is not straightforward to compare these results with Pierazzo et al. (2008) and the experiments. To carry out appropriate comparison, we need to perform 3D impact tests or use 2D axisymmetric domain. However, note that DISPH does not show unphysical behavior compared to the results with grid code. We need further investigation to find an appropriate treatment of the free surface, which is left for future work.  The velocity is normalized by v imp (ρ imp /ρ tar ) (3ν−1)/3µ while the mass is normalized by the impactor mass. The red curve indicates the result with DISPH, while the blue curve indicates that with SSPH. The black line indicates reference line for a power-law with a slope of v −3µ .