Numerical aspects of Giant Impact simulations

1.5 MB·

by  Christian Reinhardt

MNRAS000,1–??(2016) PreprintJanuary31,2017 CompiledusingMNRASLATEXstylefilev3.0
Numerical aspects of Giant Impact simulations
Christian Reinhardt,1(cid:63) Joachim Stadel2
Institute for Computational Science, University of Zurich, 8057 Zurich, CH
AcceptedXXX.ReceivedYYY;inoriginalformZZZ
ABSTRACT
7
In this paper we present solutions to three short comings of Smoothed Particles Hy-
1
drodynamics(SPH)encounteredinpreviousworkwhenapplyingittoGiantImpacts.
0
2 FirstweintroduceanovelmethodtoobtainaccurateSPHrepresentationsofaplanet’s
equilibrium initial conditions based on equal area tessellations of the sphere. This al-
n
lows one to imprint an arbitrary density and internal energy profile with very low
a
noise which substantially reduces computation because these models require no relax-
J
ation prior to use. As a consequence one can significantly increase the resolution and
8
more flexibly change the initial bodies to explore larger parts of the impact param-
2
eter space in simulations. The second issue addressed is the proper treatment of the

] matter/vacuumboundaryataplanet’ssurfacewithamodifiedSPHdensityestimator
P that properly calculates the density stabilizing the models and avoiding an artificially
E low density atmosphere prior to impact. Further we present a novel SPH scheme that
. simultaneously conserves both energy and entropy for an arbitrary equation of state.
h
This prevents loss of entropy during the simulation and further assures that the ma-
p
- terial does not evolve into unphysical states. Application of these modifications to
o impact simulations for different resolutions up to 6.4·106 particles show a general
r agreement with prior result. However, we observe resolution dependent differences in
t
s the evolution and composition of post collision ejecta. This strongly suggests that the
a
use of more sophisticated equations of state also demands a large number of particles
[
in such simulations.

1
Key words: hydrodynamics – planets and satellites: formation – planets and satel-
v
lites: terrestrial planets
6
9
2
8
0
1 INTRODUCTION tool to study the outcomes of such violent collisions. While
.
1 most of the earlier work (e.g. Benz et al. (1987), Cameron
0 Collisions are one of the most fundamental processes in & Benz (1991), Canup & Asphaug (2001), Canup (2004b),
7 planet formation as they are a key mechanism to explain
Wada et al. (2006)) concentrated on the giant impact hy-
1 growthofinitiallymicronsizeddustgrainstomassiveplan-
pothesisoftheoriginofthemoon(Cameron&Ward1976)
: ets like the ones we observe today in our solar system. The
v morerecentsimulationsexploreavarietyofpossiblescenar-
last stage of terrestrial planet formation is dominated by
i ios. Some recent examples include satellite formation in a
X very energetic collisions between roughly Mars-sized plane-
moregeneralsense(e.g.Canup(2011),Citronetal.(2015)),
r taryembryoscalledGiantImpacts(GI)(Agnoretal.1999). mantlestrippingevents(Benzetal.2007)orderivingmerg-
a SuchGIplayakeyroleintheprocessofplanetformationas
ing criteria to improve N-Body simulations (Genda et al.
they determine and influence many important properties of
2012; Leinhardt & Stewart 2012).
the final planetary system like the final number, mass and
spin state (Agnor et al. 1999), chemical and physical com- MostpriorworkusestheSmoothedParticlesHydrody-
position (Benz et al. 2007; Wiechert et al. 2001; Genda & namics (SPH) method to solve the equations of motion of
Abe 2005; Tonks & Melosh 1992; Jutzi et al. 2013) and the thematerial.SinceSPHinterpolatesthephysicalquantities
presenceofsatellites(e.g.Benzetal.(1986),Canup(2004b), using a discrete set of particles that evolve with the flow, it
Citron et al. (2015)). cansimulateverydeformedgeometriesandtracktheejected
Starting with Benz and Cameron’s (Benz et al. 1986) materialoveralargedynamicrangemakingitanidealtool
pioneering work in the mid 80’s simulations of GI between forsuchsimulations.Fortheresearchpresentedinthispaper
planetary embryos have became an increasingly important we modified the cosmological hydrodynamics code GASO-
LINE (Wadsley et al. 2004) that has been extensively used
in many astrophysical applications and applied it to such
(cid:63) [email protected] GI simulations. The code uses a modern SPH implementa-
(cid:13)c 2016TheAuthors
2 C. Reinhardt and J. Stadel
tion, a tree to calculate gravity and has a modular design toGIsimulationsarepresentedinSection5.Finallyasum-
thatallowsonetoimplementnewphysicswithouttouching mary of the present work and conclusions are discussed in
the complicated parallel layer. We have implemented the Section 6.
Tillotson equation of state (EOS) (Tillotson 1962) within
GASOLINE to model the condensed materials that planets
are made of. To facilitate the use of the Tillotson EOS we
developedalibrary1 thatprovidesfunctionstocalculatethe 2 CREATING QUIET PARTICLE
REALIZATIONS OF PLANETS
pressure, sound speed and temperature allowing its use in
anysimulationsoftware.BesidestheSPHmodificationsde- InSmoothedParticlesHydrodynamics(SPH)thedensityof
scribedinthispaperwealsoaddedanewclassofsmoothing eachparticle(aswellaspressureforcesandchangesininter-
kernels (Dehnen & Aly 2012) that do not suffer from parti- nal energy) is derived from the positions of its neighbours.
cleclumpingandallowincreasingthenumberofneighbours Thusitiscrucialtohaveanaccurateparticlerepresentation
thereby reducing the noise of the method. oftheinitialconditions.Aninhomogeneousdistributionwill
In this paper we present solutions to three short com- causelargespuriousradialandangulardensityfluctuations
ings of SPH that previous work encountered. First of all it that cause oscillations and requires careful relaxation. Hav-
is very difficult to get an accurate SPH representation of ing Poison noise in the IC not only means that we waste
a planet’s equilibrium initial condition (IC) so one has to computationaltimeonrelaxingthemodelsbutisingeneral
evolve (”relax”) the system carefully to reach a true equi- problematic for SPH as it affects the method’s convergence
librium state before doing any impacts. This is very time (Monaghan 1985). For larger simulations with several mil-
consuming and can, in the worst case, use more computa- lion particles, relaxation can take a substantial part of the
tionalresourcesthantheactualsimulations.InSection2we total simulation time and doing such simulations becomes
presentanewmethodtogenerateextremelylownoiseplan- prohibitiveiftheICistoonoisy.Thismakesthetraditional
etarymodelsthatdonotneedrelaxationpriortoacollision approach laborious for parameter studies where the target
which drastically reduces the total simulation time. and impactor properties are varied.
Another issue is the proper treatment of the material– Previousworkplacedtheparticlesonauniform3Dlat-
vacuumboundaryataplanet’ssurface,wherestandardSPH ticeeitherwithconstantinternalenergy(Gendaetal.2012)
tendstoseverelyunderestimatethematerial’sdensitycaus- or e.g., an isentropic thermal profile (Canup 2012), inside a
ing it to jump from the condensed to an expanded state. sphere to have a very uniform and thus low noise particle
Thiscreatesan”atmosphere”oflowdensitymaterialaround distribution. To obtain a density gradient that is consistent
a planet prior to the impact. A solution to this problem is with the thermal profile and get rid of the Poison noise at
discussedinSection3wherewepresentamodificationofthe theboundary,thesemodelsarethenevolvedwiththehydro-
standardSPHdensityestimatorthatproperlycalculatesthe dynamicscode(usuallywithavelocitydamperthatreduces
densityformaterial–vacuumboundariesencounteredatthe the particles velocities by a given fraction after each step)
surfaceoftheplanets.Anotherpossibleremedywouldbeto untilanequilibriumrepresentationisreached.Thismethod
apply DISPH (Saitoh & Makino 2013) to GI simulations as has two short comings: 1) it does not allow one to imprint
donebyHosonoetal.(2016).BecausethisSPHschemedoes a (radial) density gradient and 2) also has severe noise at
not use the density but pressure as a smoothed (and thus the planet’s surface where the grid is not adapted to the
differentiable)quantity,itperformswellatcontactdisconti- spherical symmetry of the profile. Some authors have tried
nuitiesandfreesurfaces,e.g.,encounteredatthecore-mantle to model a density gradient by distorting the uniform grid
boundary and the planet’s surface. In the present work we radially (Woolfson 2007) but the resulting models were out
focusonundifferentiatedbodies.However,thepropertreat- of equilibrium and still needed relaxation.
mentofthecore–mantleboundaryisnotaproblemandthe In the present paper we suggest a different approach
improvements described here have already seen application that respects the spherical symmetry of the problem and
in realistic models with an iron core and basalt mantle. providesaveryuniformradialparticledistributionbyusing
Finally,SPHdoesnotnecessarilyconservebothenergy an equal area tessellation of the sphere. This method not
and entropy for an arbitrary equation of state. In the case only produces very low radial and angular density fluctu-
ofanidealgasonecandefineanentropyfunction(Springel ations but it also allows one to imprint a density gradient
& Hernquist 2002) which allows one to conserve entropy thatcloselyfollowstheequilibriumsolutionthusmakingre-
forpurelyadiabaticflows.Foramorecomplex,maybeeven laxation of the models obsolete.
tabularEOS,thisisingeneralnotpossible.InSection4we PriortogeneratingtheSPHrepresentationoftheplane-
present a new, entropy conserving SPH scheme that can be tarybodiesonehasto(numerically)solvetheusualinternal
applied to any equation of state. This prevents loss of en- structureequations(e.g.,Alibert(2014)).Forourboundary
tropy during the simulation and assures that the material conditions(e.g.,M(r=R)=M andρ(r=R)=ρ )itis
tot 0
does not evolve to unphysical states thus also eliminating convenient to solve these for ρ(r), M(r) and u(r). Besides
code crashes resulting from particles reaching negative in- the equation for hydrostatic equilibrium
ternal energies. Code tests and applications of the methods
∇P(ρ,u) M(r)
=−G (1)
ρ r2
1 The TillotsonEOS library is freely available under https:// oneneedstospecifyanequationofstateandaninternalen-
github.com/chreinhardt/tillotson.git. We also provide our ergyprofiletohaveaclosedsetofequations.Forthepresent
codetocreatelownoiseparticlerepresentationsofplanets,bal- workweusedtheTillotsonequationofstate(EOS)(Tillot-
lic,athttps://github.com/chreinhardt/ballic.git. son 1962) that was originally developed to model hyper-
MNRAS000,1–??(2016)
Numerical aspects of Giant Impact simulations 3
velocity impacts and has been used in many prior simula- procedure results in somewhat more or somewhat less par-
tions (e.g. Benz et al. (1987), Canup & Asphaug (2001), ticles than specified, meaning that the total mass is deviat-
Marinova et al. (2011), Genda et al. (2012) and Jutzi et al. ing from the desired value. We then recalculate the particle
(2013)). Despite the simple analytic form the results are massgiventheobtainednumberofparticlesandrepeatthe
in good agreement with measurements (Benz et al. 1986; above until we converge. Our final model has the exact to-
Brundage 2013) and its ability to properly reproduce the tal mass and equilibrium density profile (Figure 4) for the
materialsHugoniotcurve,thusaccuratelymodellingshocks, planet as well as having equal mass for all particles. There
is excellent (Brundage 2013). More details on the EOS and are somewhat larger density fluctuations in the inner–most
the material parameters used in this study can be found in shells because the number of particles, and thus the radial
Appendix A. size of the cells, can not be adapted as flexibly due to con-
Thesimplestinternalenergyprofileisuniformandeach straintsfromtheHEALPixgrid.HEALPixuses12“rectan-
particle has the same internal energy. This works well for gular”patches on the sphere each of which can be further
bodies smaller than a few hundred kilometres where the subdivided into n×n pixels. An Icosahedral grid uses 20
densityisroughlyconstantbecauseself–gravityistooweak “triangular”patches,meaningthattheinner–mostshellhas
to significantly compress the material (e.g. Genda et al. 20particleswhichissomewhatmorerestrictivethanthe12
(2015)). More massive planets have increasingly steep den- particle shell that is possible with HEALPix.3 We also no-
sity gradients, where the use of a uniform thermal profile ticethatthesmootheddensitydeviatessignificantlyforthe
wouldbeunphysicalandfurthermoresomeparticlestendto imprinted profile in the outer most shells. There the result-
end up below the minimal allowed energy for the material ing density as determined by the SPH code is lower than
(thesocalledcoldcurve)2 whenstartingfromasuchinitial that originally imposed, which causes the bodies to be out
conditions.Suchaprofileisalsonotisentropicmeaningthat of equilibrium. For a condensed material EOS this is even
themodelisexpectedtoexhibitconvectionuntilaninternal more problematic as the estimated surface density is lower
energy profile with constant entropy is reached. than the material’s reference density (ρ , the density of the
0
WeimprintedanisentropicthermalprofileinourIC,so materialatzerotemperatureandpressure)whichcausesthe
thatthesurfacetemperaturesarebelow2000K.Thisseems materialtofallintoanunphysicalregime.Foracoldermodel
sensible in light of the fact that we expect a ”cold”surface withalowsurfacetemperaturethisislessproblematicasthe
because the time between GI is typically large enough to TillotsonEOSinthisregionhasthesameanalyticformasin
allow the planet’s surface to cool. Of course, it is possible the condensed states. For hotter models, however, the par-
to solve for higher surface temperature models in the case ticlesendupintheintermediateexpandedstateswherethe
thatthetimebetweencollisionsisshortorthecoolingtime EOSisdifferentandthemodelsbecomeveryunstable.Itis
is long. thereforecrucialtocorrectthesurfacedensityascalculated
To create the particle realization of an equilibrium by the SPH code and is made possible by the method we
modelwedividedthesphereintoconcentricsphericalshells present in the following section.
anduseanequalareatessellationofthesurfaceofeachshell. Inordertoassestoperformanceof ballic,wecompare
We experimented with two tessellation methods, the Icosa- it to two particle distributions of the same model, that are
hedronpackage(Tegmark1999)andHEALPix(Go´rskietal. obtained using more conventional methods (Figure 2). Fol-
2005)thatwasoriginallydevelopedforanalysingdatafrom lowingpreviouswork,weuniformlydistributedtheparticles
cosmicmicrowavebackgroundmeasurements.Bothmethods on a Cartesian grid inside of a sphere, so that a particle’s
in principle suit our purpose. Icosahedral tessellation pro- densityisequaltotheuncompresseddensity,ρ ,ofthema-
0
duces a nearly perfect tangentially homogeneous grid but terial. Of course such a uniform density is only physical for
suffersfromlargerradialdeviations.HEALPixontheother smaller bodies, as we expect larger bodies to compress un-
handproducesabetterradialdistributionbutshowsasmall dertheirowngravityresultinginadensitygradient,andan
amountoftangentialdensityfluctuationsduetoartefactsre- Earth-mass body is thus poorly modelled. We then imprint
sulting from the grid. Since, in our case, radial fluctuations either a uniform (Genda et al. 2012) or isentropic (Canup
are more problematic and the HEALPix grid adapts more 2004b) internal energy profile and let them evolve with our
flexibly in the inner part of the sphere (see below), it has codeuntilaselfconsistentequilibriumisreached.After30h
been used exclusively in all of the work presented here. (in simulation time) we determine each particle representa-
We begin particle generation by first solving the above tion’s root mean square velocity, which measures the parti-
structure equations (Equation 1). Our method then itera- clenoise.Asinpreviouswork(Canup2004b),(Gendaetal.
tivelyminimizestheaxisratiosofthefinitevolumessothat 2012), (Marinova et al. 2011) we consider a model to be
the tangential to radial size become almost identical by ad- relaxed, if its root mean square velocity is below a small
justing the number of particles for each shell and the spac- fractionoftheimpactvelocity,inourcasev <50ms−1.
RMS
ing between shells keeping the particle mass constant. This We also compare the relaxed body with the desired pro-
file, to see how close the different methods are to the input
model. The results are striking: the initial distribution ob-
2 Itshouldbenotedthatunlikeinthecaseofanidealgas,isother- tained from ballic can be considered to be relaxed right
mal does not mean that the internal energy is constant as there from that start and even after 20h in simulation time still
isanadditionalcontributionfromthecoldcurvetothetotalen-
followstheimprinteddensityprofile.Theothertwomethods
ergy in case of a condensed material EOS (see Appendix A for
details).Issueswithparticlesfallingbelowthecoldcurvedueto
anunphysical,uniforminternalenergyprofilewouldalsoreadily
provide an explanation for the problems with ANEOS reported 3 We also include a single central particle, although this is op-
inCanup(2004b) tional.
MNRAS000,1–??(2016)
4 C. Reinhardt and J. Stadel
Figure 1. A slice (2.2 x 2.2 x 0.2 R⊕) through the SPH repre-
sentation of an equilibrium model (0.997M⊕ and 105 particles)
where the colours represent the internal energy of the material.
Thesmallplotshowsazoomatthesurfaceofthemodel,where
wemarkedthesmoothingkernel(inred)ofaparticleintheouter
mostshell.Obviouslysuchoutermostparticleshaveonlyoneside
of their kernel sampled by neighbouring particles while the rest
isvacuum.
performsignificantlyworsebothintermsofrelaxationtime,
which is about 10h, and match to the original model. This
is not surprising as these particle representations both suf-
fer from Poison noise and start from unphysical initial con-
ditions. We also experimented with re-imprinting an isen- Figure 2.Inthisfigurewecompareballic(red)toamorecon-
tropic internal energy profile after every few steps of the ventional method that places particles on a Cartesian grid with
a uniform (green) or isentropic internal energy profile (blue). In
simulation to enforce a closer match to the desired model
bothcasesweevolvetheparticlerepresentationofanEarth-mass
and randomly distributing particles in concentric shells as
granitemodelusingISPH(Section4)andthedensitycorrection
inMarinovaetal.(2011)butbothprovedtobedifficultand
proposedinSection3.Thetopframe(a)showshowtherootmean
time consuming as particles would enter unphysical states
squarevelocityofthemodelschangeswithtime.Thedashedgrey
(see Section 4 for details) which caused the hydro code to line marks vRMS =50ms−1, below which models have typically
crash. beenconsideredtoberelaxed(vRMS <5·10−3vimp,withanas-
sumedimpactvelocityofvimp=10ms−1).Theplotclearlyshow
thataninitialdistributionmadewithballicisrelaxedrightfrom
thestart,whiletheothertwomodelsrequireabout10hinsimu-
3 A MODIFIED DENSITY ESTIMATOR FOR lationtimetoreachthispoint.Evenafterthattheyremainnoisy
FREE SURFACES and show oscillations in vRMS. In the bottom frame (b) we plot
thebinneddensityoftherelaxedbodies(t=20h)andcompareit
InSmoothedParticleHydrodynamicsthedensityofaparti-
to the model (dashed grey line). The Cartesian initial distribu-
cleisderivedfromthemassesandpositionsofitsneighbour- tionsclearlydeviatefromthemodel,expandandhavealayerof
ing particles (Monaghan 1992; Springel 2010; Price 2012) low density material at the surface. Our initial conditions show
an excellent match to the imprinted density profile and have no
N
(cid:88) suchlowdensity“atmosphere”.
ρ = m W(|r −r |) (2)
i j i j
j=1
wherethekernelfunctionW(r)weightseachparticle’scon-
tributionaccordingtoitsdistancefromthecentralparticle.
thestandardSPHdensityestimatorassumesthatthewhole
ModernSPHimplementationsuseacompactkernelanden-
kernel is (more or less homogeneously) filled with particles.
force a fixed number of nearest neighbour particles (Mon-
One possible treatment of vacuum/material interfaces
aghan 1992) or enclosed mass (Springel & Hernquist 2002)
is thus to correct the density using a correction factor
toautomaticallyadapttheresolutionofthemethod.When
thedensityisestimatedinthiswayaparticleattheplanet’s V
V = (3)
surfacehasaboutonehalfofitskernelsampledwithneigh- fac V
eff
bours while the rest of the volume is empty, or rather filled
with vacuum. This means that the mass to volume ratio, that only accounts for the effectively sampled volume V .
eff
and thus the resulting density, is underestimated because OnewaytoestimateV istoconsidertheimbalanceofthe
eff
MNRAS000,1–??(2016)
Numerical aspects of Giant Impact simulations 5
particle distribution in the kernel
|f |
fimb,i = 2h (cid:80) miW (h ) (4)
i j j ij i
where we used
(cid:88)
f = (r −r )m W (h ) (5)
i i j j ij i
j
to measure the asymmetry of the particle distribution in-
side the kernel. We also tried to estimate the imbalance as
suggested by Woolfson (Woolfson 2007) but found that the
lower order of this estimator results in density fluctuations
larger than the required correction. By using a kernel aver-
age (we use the Wendtland C , a fourth order kernel, with
2
80nearestneighbours)asshownintheaboveequationsour
resultingdensityatthesurfaceismuchlessnoisy.Themod-
ified density is then obtained from
ρ
ρ = SPH (6)
eff V
fac
whichobviouslygivesthestandardSPHdensityiftheparti-
cle distribution inside of the kernel is symmetric (V =1).
fac
For a given filling factor of the kernel, where we assume a
planeboundarybetweenmassandvacuum,wecancalculate
the corresponding f by integrating Equation 4 numeri-
imb,i
cally over the filled region. This is shown by the blue dots
in figure 3.
Sincetheresultwastoosensitivetonoiseforverysmall
andlargeimbalancesweusedalinearapproximation(shown
in green in Figure 3) and kept V constant in these cases.
fac
Since the proposed density correction only depends on Figure3.AplotthatshowsV (f )forthecubicspline(top)
fac imb
a given particle’s kernel volume, it is not symmetric (e.g., and Wendland C2 kernel (bottom). Rather than using the ob-
Equation 4 and 5) and we have to calculate the standard tained numerical result (blue dots) we fit a linear function (red
SPH density from andgreenlines)andkeepthevolumecorrectionconstantforsmall
andlargeimbalancestoavoidartificiallytriggeredimbalancedue
(cid:88)N to noise or large density gradients. The coefficients we obtained
ρ = m W (h ) (7)
i j ij i for the linear fit are a = −2.39701, b = 1.02918 for the cubic
j=1 splineanda=−2.52444,b=1.0259fortheWendlandC2kernel.
which differs from that typically used in GASOLINE Wad-
sley et al. (2004). It is not clear how to relate imbalance
to a V in a consistent way when using GASOLINE’s
fac
symmetrized density, a formulation which improves the be-
haviourforverystrongdensitydiscontinuitiesbecauseapar-
ticleinthehighdensityregionstillgetsacontributionfrom morethantheirequilibriumdistance(Melosh2007).Sincea
lowdensityparticles.Correctlyhandlingmaterialinterfaces, fluid at these low densities (and temperatures) is expected
e.g., at core–mantle boundaries, is of greater relevance for tofragmentand,forexample,formdropletsratherthanbe-
GIs than this refinement to the density definition. Further- havinglikeacontinuumwefollowMelosh(Melosh1989)and
more,weuseamuchlargernumberofneighbours,forwhich set the pressure to zero if it would become negative.
itisevenlessofanissue.However,wedostilluseasymmet- The resulting particle representations using the modi-
ricformulationtoevaluatetheequationofmotionsuchthat fied density estimator are in excellent agreement with the
both momentum and energy are conserved (e.g. Monaghan model (Figure 4) and the particles in the outer most shell
(1992),Wadsleyetal.(2004),Springel(2010),Price(2012)). remaininthecondensedstate.Thesmootheddensityprofile
SincetheusuallyusedB-splinekernels(seee.g.Wadsley exhibits a general trend to over estimate the density com-
etal.(2004)orMonaghan(1992))cantriggeraninstability paredtoclassicSPH,whichismostpronouncedintheouter
that causes particles to clump when they are getting close most shell. This occurs because the proposed estimator ar-
weimplementedtheWendlandkernelsproposedbyDehnen tificially triggers an imbalance due to a density gradient on
&Aly(2012)thatarebothstableandallowalargereffective the kernel scale. For less steep density gradients, e.g., for
numberofparticlesthusdecreasingthenoiseforagivenres- lower mass models, the corrected density lies even closer to
olution.Anothercauseforclumping,independentofthecho- thedesiredvaluesatthesurface.Inordertoproperlyrepre-
senkernelfunction,isanegativepressurebetweenparticles. sentverymassiveplanets,likesuper–earths,onewouldhave
This isproblematicas theTillotsonEOScalculatesanega- to modify the procedure that relates V to f to ac-
fac imb,i
tivepressureinthelowdensity,coldregion(AppendixA)to count for a density gradient. At this point we leave further
mimicatensionalforcebetweenparticlesthatareseparated improvements of the estimator to later work.
MNRAS000,1–??(2016)
6 C. Reinhardt and J. Stadel
(cid:18) (cid:12) (cid:19)v (cid:18) (cid:12) (cid:19)
ui,j+1, ∂∂uρ(cid:12)(cid:12)i,j+1 ui+1,j+1, ∂∂uρ(cid:12)(cid:12)i+1,j+1
(cid:18) (cid:12) (cid:19) 3 u(ρ,v) (cid:18) (cid:12) (cid:19)
ui,v, ∂∂uρ(cid:12)(cid:12)i,v Hermitefunction ui+1,v, ∂∂uρ(cid:12)(cid:12)i+1,v
1 2
ρ
(cid:18) (cid:12) (cid:19) (cid:18) (cid:12) (cid:19)
ui,j, ∂∂uρ(cid:12)(cid:12)i,j ui+1,j, ∂∂uρ(cid:12)(cid:12)i+1,j
Figure 5. A schematic explaining the 2D interpolation method
used in ISPH. The values of u(ρ,v) and ∂u/∂ρ| are calculated
v
from Equation 13 while their derivatives with respect to v are
obtainedbyfittingcubicsplinesinv.Ateachgridpointu(ρ,v),
∂u/∂ρ|v, ∂2u/∂v2(cid:12)(cid:12)ρ and ∂3u/∂ρ∂v2(cid:12)(cid:12)ρ are stored in the lookup
Figure 4. A density profile of the particle representation of an table. The value of u(ρ,v) is determined by (1) doing a spline
M = 0.997M⊕ protoplanet using N = 105 particles. The solid interpolation in v at ρi to calculate ui,v and ∂u/∂ρ|i,v, then
red line shows ρ(r) for the calculated equilibrium. The dots are (2) doing the same at ρi+1 to obtain ui+1,v and ∂u/∂ρ|i+1,v.
the smoothed density for the classic SPH (light blue) and our Step(3)istointerpolateu(ρ,v)betweenthesevaluesusingthird
new density estimator (light green). While the results agree for order Hermite functions. Since ∂u/∂v|ρ and ∂2u/∂ρ∂v can be
both methods in the inner part of the sphere, the particles in determinedfromthesecondderivativesusedforthecubicspline
thelastshellmatchtheimposedsurfacedensitymuchbetterfor method ∂u/∂v|ρ,vcanbecalculatedinthesamelookupifneeded.
our improved estimator. At this material/vacuum interface only
afractionofthekernelissampledwithSPHparticlessowecor-
rect the standard SPH density accordingly. Besides avoiding an
where
artificial ”atmosphere”of low density material this also provides
innecerdeeadsedprsiotarbtiolityusteoitnheaesqimuiulilbartiiuomn.mWoedeallss,osoobnsoerrveelaxthaatitonouisr dudtad = Pρ2i (cid:88)n mjvij·∇iWij (9)
proposed modification tends to slightly overestimate the density i j=1
compared to the desired values because imbalance of the parti-
is the work from adiabatic compression or expansion (PdV
cle distribution associated with the density gradient triggers an
over–correction. work) and
n
duΠ = 1(cid:88)m Π v ∇ W (10)
dt 2 j ij ij i ij
j=1
4 ENTROPY CONSERVING SMOOTHED isthecontributionduetoirreversible(ds>0)shockheating.
PARTICLES HYDRODYNAMICS The artificial viscosity term Πij (AV) is non zero only for
convergingflowsandcapturesshocks(e.g.Monaghan(1992),
The results with the modified density were improved but Wadsleyetal.(2004),Springel(2010),Price(2012)).Follow-
some simulations would crash, e.g., due to unphysical en- ingpreviouswork(e.g.Canup(2004b))weusedα=1.5and
tropyloss,eveniftheoriginalICsweresetupinaphysically β =2α for the viscosity parameters.
consistentwayasdescribedabove.Forexample,ejectedma- In the case of an ideal gas it is possible to define an
terialthatre–impactswiththetargetwouldallofasudden entropyfunctionA(s)=(γ−1)uρ1−γ (Springel&Hernquist
see a high density region when its kernel starts to overlap 2002)sothatonecanalternativelyevolvetheentropyrather
with the target, causing a sudden and unphysical increase than the internal energy given by
intheparticle’sdensitywhiletheinternalenergystayscon-
dA γ−1 γ−1du
stant (Figure 6). This would not happen, if the method i =− L(ρ ,u )+ Π, (11)
would exactly conserve the particle’s entropy since its in- dt ργ i i ργ−1 dt
i
ternal energy would increase accordingly.
where L is a source term. This has the advantage that in
Clearly evolving the entropy rather than the internal
the absence of dissipative processes the entropy is exactly
energy as a fundamental thermodynamical variable would
conserved(L=0).Thisisnotthecaseiftheinternalenergy
bedesirable.Aparticlewillfollowanisentropeforapurely
is evolved (Springel & Hernquist 2002). The main difficulty
adiabatic flow thus making it impossible to ever fall below
of this entropy formulation is, however, that it requires the
the cold curve (if we start from physical initial conditions)
EOS to be both analytic and thermodynamically complete
becausetheonlypossibilityistomovetoahigherisentrope
inordertodefineanentropyfunction.Thisisoftennotthe
as a result of shock entropy generation.
case for a more complex EOS like the Tillotson EOS that
Initsstandardversion,GASOLINEevolvestheinternal
weuse.Itispossiblehowevertocalculatetheisentropesfor
energy as (Wadsley et al. 2004)
an arbitrary EOS from the first law of thermodynamics
du du du P
i = ad + Π (8) du= dρ+Tds (12)
dt dt dt ρ2
MNRAS000,1–??(2016)
Numerical aspects of Giant Impact simulations 7
where ds=0 for constant entropy. So one can solve wherea =a(x )anda =a(x )aretheaccelerations
i i i+1 i+1
due to the forces acting between the particles.
du P
= (13) Implementing such an integrator for a fluid is slightly
dρ ρ2
morecomplicatedbecausetheforcesnotonlydependonthe
numerically and store the results for different initial values positionbutalsothevelocityandinternalenergyofthepar-
(taken at u(ρ=ρ0)) in a lookup table. Since this does not ticles(Quinnetal.1997).Thesequantitiesareknownatthe
requireananalyticexpressionofthepressurebutjustanu- startbutnotattheendofeachtimestepsoinordertocal-
mericalvalueitcanbedoneforanyEOS,evenifitistabular. culatethefluidforcesonehastouseapproximatepredicted
Intheabsenceofentropychangingprocesses,likeshocks,the values to update the velocity from v to v .
i+1/2 i+1
evolution of the internal energy for a given change in den- WenowdescribeasingleISPHtime–step.1)Thegrav-
sitycanthenbesimplydonebyinterpolatingbetweenthese itational and SPH forces on all particles a(t ) and contri-
i
values. butions from shock heating (Equation 10) are calculated
To setup the lookup table we integrated Equation 13 at the beginning of the time interval. 2) Then the veloci-
fordifferentstartingvaluesvi =ui(ρ=ρ0)inthedirection ties and internal energies are updated for half a time step
of both increasing and decreasing densities. The variable vi (”kick”). 3) Now all particle positions are ”drifted” from
thus can be used to label an isentrope and v0 = 0 corre- x(ti) to x(ti+∆t). Before calculating the updated density
sponds to the cold curve. For a given isentrope a particle’s ρ(t +∆t) at the end of the interval its old value ρ(t ) is
i i
internalenergythusdependsonlyonit’sdensityu=u(ρ,v). stored. 4) Using this value and the new density, calculate
So in order to evolve a particle’s energy from (ρ1,u1) to the change in internal energy due to adiabatic expansion
(ρ2,u2)wefirsthavetodetermineonwhichisentropeitlies. or compression using the lookup table as described above.
For this we calculate the root of 5) In order to do the closing ”kick” the predicted velocity
and predicted internal energy at the end of the interval are
F(v)=u −u(ρ ,v) (14)
1 1 determined. 6) The forces a(t +∆t) and shocking heating
i
which is bracketed by v and v using the Brent algo- are calculated.7)Finallyboththevelocityandinternalen-
0 max
rithm. To find u(ρ,v) we do a two dimensional bi-cubic in- ergyareupdatedusingthenewlycalculatedforceandshock
terpolation in the look-up table as described in Figure 5. heating term (closing ”kick”). This way of evolving the en-
Since the derivatives of u with respect to ρ are known an- ergy equation matches very well with our time integration
alytically we only have to fit cubic splines along the direc- scheme and can be integrated rather easily into the code.
tionofv(atconstantdensity).Onceaparticle’sisentropeis The actual algorithm is more complex than this because
foundwedothefinalinterpolationofu =u (ρ ,v)forρ . the particles have their own unique time steps based on a
2 2 2 2
Todorealistic(impact)simulationsonehastoaccount block time stepping scheme. While particles are all drifted
forthechangeofentropyduetoshocksthuswestillneedto simultaneously on the smallest time step, such that the po-
accountforthedu /dtterm.InSPHitispossibletoobtain sitions, and hence densities, are always time synchronized,
π
the internal energy contribution due to shocks from Equa- thoseonlargertimestepswillnothavetheirthermalenergy
tion10.Thisallowsonetosplitthe(temporal)evolutionof advanced(seestep4above)untiltheyreachtheendoftheir
the internal energy into an adiabatic and a shock heating time step.
part, enforcing the conservation of entropy for an adiabatic
flow and including the contribution of shock heating to the
internal energy using the standard SPH formalism. This is
5 TESTS
very much in the spirit of a leapfrog–like evolution of the
thermal energy and meshes well with the existing adaptive Asafirsttestwesimulatedtheequilibriummodelswithour
time–stepping leapfrog scheme used within GASOLINE. code to check if they truly stayed in equilibrium or would
One time step from ti to ti+1 in the ISPH adaptation expandorcollapseafteracertaintime.ForclassicSPHthe
of the leapfrog algorithm used in GASOLINE4 is given by, modelswouldremainstablebutoscillate,astherewassome
intrinsic noise in our model and the density at the surface
1
kick: vi+1/2 =vi+ 2ai∆t of the body was significantly underestimated as discussed
in Section 3. After several oscillations the particles would
1du
kick(entropy): u(i1) =ui+ 2 dtΠ∆t settleintoanequilibriumconfigurationbecausetheartificial
viscosity in SPH would dampen their motion since it does
drift: x =x +v ∆t
i+1 i i+1/2
not distinguish between a converging flow and a shock. As
du
drift(adiabatic): u(2) =u(1)+ ad∆t in previous studies (e.g. Canup & Asphaug (2001), Genda
i i dt
et al. (2012)) we used a velocity damper that reduced the
1
kick: v =v + a ∆t particle’s velocity each step to speed things up but this has
i+1 i+1/2 2 i+1
to be done carefully to not over–damp the system as the
1du
kick(entropy): u =u(2)+ Π∆t models will become unstable again as soon as the damping
i+1 i 2 dt
term is removed. For simulations that used our proposed
density correction the models remain stable right from the
startanddonotchangeatall.UsingISPHdoes,asexpected,
4 HeretheKick-Drift-Kick(KDK)versionofthealgorithmisde-
not affect the result in any way.
scribed.Thismeansthatboththepositionsandforcesareevalu-
atedatthebeginningandtheendwhilethevelocitiesareupdated Wethencheckedthecode’sabilitytoproperlymodelan
in the middle of the time step. This is convenient as the forces adiabatic flow which is not only a key requirement for any
areevaluatedatpositionsthataresecondorderaccurate. hydrodynamics code but also most important to test how
MNRAS000,1–??(2016)
8 C. Reinhardt and J. Stadel
theparticlesintheISPHsimulationcloselyfollowtheirini-
tial isentrope over more than 17 hours in simulation time
showing excellent entropy conservation.
Next we verified that we can recover the Rankine-
Hugoniot jump conditions (e.g. Melosh (1989)) relating the
fluidquantitiesinfrontandbehindthepassingshocktotest
thecode’sabilitytocorrectlycaptureshocks.Forthiswelet
twouniformgraniteslabs(15×15×8R )thatareinitially
⊕
incontactcollidewithoppositevelocities.Theinitialdiscon-
tinuityinthevelocityfieldthensendsashockwavethrough
the material causing the fluid variables to jump to the post
shockstateattheshockfront.Thesevaluescanthenbeob-
tainedfromthesimulationsandcomparedtothetheoretical
resultsusingtheRankine-Hugoniotequations.Forbothclas-
sic SPH and ISPH (with and without the modified density
estimator) we obtain the correct values and thus conclude
that shocks are adequately resolved in our code.
Figure 6. An example of a simulation (b = 0, v = 10km/s, As a last test we simulated collisions between granite
N = 105) between an Earth-mass target and a Mars-sized im- spheres for different impact parameters and compared our
pactor where the code crashed because a particle starting from results to prior work done in the field. Since there are not
(S)enteredtheunphysicalregionbelowthecoldcurve(C)ofthe
many publications on GI between undifferentiated bodies
TillotsonEOS.Thishappensbecauseanejectedclumpofmate-
and the early work (Benz et al. 1986) suffers from (for to-
rialreimpactswiththetargetandtheparticle’skernelsuddenly
day’sstandards)lowresolutionandshortsimulationtime,a
overlaps with the target, which causes an unphysical increase in
direct comparison of the results was complicated. However
the density while the internal energy stays constant. If ISPH is
applied to such a simulation, a particle can never evolve in such weexpectthebasicfeaturesofsuchacollisiontobesimilar
an way and the simulation does not crash. We also marked the fordifferentiatedandundifferentiatedbodies.Firstweruna
fourregionsoftheTillotsonEOS(I,II,IIIandIV)definedbythe simulation of an oblique impact (b=0.71, v =10km/s,
imp
Tillotsonparametersρ0,us andu(cid:48)s asdiscussedinAppendixA. Mt =0.997M⊕ andMi =0.099M⊕)liketheonethatmight
have created our moon (Canup & Asphaug 2001; Canup
2004b) using a total number of 105 particles for both clas-
u
sic SPH without any of the modifications presented in this
paper and ISPH (including the density correction discussed
inSection3).Priortotheimpactbothbodiesareplacedin
ui+1
∆uπ the xy-plane at a given separation so that they collide with
the desired impact angle and velocity. As in previous sim-
ulations the impactor is deformed due to tidal forces while
approaching the target. The impact sends a strong shock
through both bodies destroying most of the impactor and
the ejected material forms an arm like structure for both
∆u
ad
methods which is consistent with prior results (Canup &
∆unum Asphaug2001).Sinceourplanetarymodelsrequirenorelax-
ationpriortotheimpactsimulationwerunthesamesimula-
tionforanincreasingnumberofparticlesforN =105,8·105
∆uπ
ui up to 6.4·106 to check the resolution dependence of the
results.Whilethegeneralmorphology(massivetarget,arm
like structure and impactor remnant) is the same the den-
ρ
ρi ρi+1 sity and internal energy of the material and thus the phys-
ical state of the material clearly depends on the resolution.
Thisismostpronouncedforthearmlikestructureasshown
Figure 7.Theevolutionoftheinternalenergyduringone(sub)
in Figure 9 and Figure 10. This makes the use of a highly
stepinISPH.Thedashedblueline∆unum representstheexact sophisticated EOS in low density, poorly resolved, regions
numerical solution of evolution of the the internal energy as a
questionable as the values of the fluid variables there are
weak shock occurs. ISPH divides the total change in u into an
highly dependant on the resolution and the treatment of
adiabaticpart∆uad andshockheating∆uπ.
the hydrodynamics itself.
Thesecondsimulationweperformedwasamorehead–
accuratelyISPHconservesentropy.Forthiswecutasphere on impact (b = 0.5) of an 0.2M impactor onto a 0.9M
⊕ ⊕
from a cubic uniform density and internal energy particle target at 10.62 km/s with 105 SPH particles using classic
grid and let it evolve only due to pressure forces (neglect- SPHandISPH.Weletthesimulationevolveforonedayand
ing gravity). We ran two simulations, one with classic SPH again find that both methods initially produce very similar
and one with ISPH both using our proposed modified den- results. They do however differ in how the ejecta evolves
sityestimator.Thedifferencesarestriking(Figure8).While on a longer time scale, as the impactor remnant leaves the
classic SPH suffers from huge scatter in the internal energy system in the classic SPH simulation but remains bound
MNRAS000,1–??(2016)
Numerical aspects of Giant Impact simulations 9
Figure 8. The adiabatic expansion of a uniform granite sphere. The left plot shows the results for classic SPH using our improved
treatmentofthematerial/vacuumboundaryfordifferenttimes(t=0sviolettot=0.7hred)whiletherightplotdoesthesameforISPH.
We clearly notice that for ISPH the particles follow their original isentrope, showing very little scatter in internal energy compared to
classic SPH. Furthermore, all scatter that is observed in the ISPH case originates from the IC (the violet points at internal energy of
2.5erg/g).
Figure 9. Simulation 1 is an example of an oblique (b = 0.71, vimp = 10km/s) collision between a M = 0.99M⊕ target and an
M =0.099M⊕ impactor (both bodies composed of granite only). The three top panels show the simulation 2.66h after the impact for
differentresolutions(N =105,8·105,and6.4·106 fromlefttoright).Thecoloursshowthedensityfrom0(deepred)to2.7g/cm3 (dark
blue)correspondingtoρ0 sodarkbluematerialiscondensedwhileeverythingelseisintheexpandedcoldorintermediatestatesofthe
EOS.AsinprevioussimulationsonMoonformingimpactsapartoftheimpactorsurvivesthedirectcollisionandformsadenseclump
attheendofan”arm-like”structure.AsshowninFigure10theimpactorremnant’sphysicalstateisverysimilarforallresolutionsand
iscomposedonlyofabout1%targetmaterial.Thethreebottompanelsshowazoomintotheabovesnapshotstoshowthe”arm”that
connectsthetargetandtheimpactorremnant.Forincreasingnumberofparticlesthearmismuchbetterresolvedresultinginobvious
differences in its structure and density profile (Figure 10). A movie showing the early phase of the impact for the 6.4·106 particle
simulationisavailableonline(https://youtu.be/QjpSg-gZl8U).
whenISPHisusedandre–impactswiththetargetafter4.4 we obtain a slightly more massive planet and less massive
days. In both cases the impact produces a disk of hot, par- disk than for classic SPH, the two methods produce very
tially vaporized material extending to about 3R from the similar disk profiles up to 24 hours in simulation time. To
⊕
target. Already after one day it assumes a flattened, rota- studythelongtermevolutionofthediskandthefateofthe
tionally supported profile which is more pronounced in the impactor remnant in the ISPH collision we continued both
case of ISPH. We also observe a steep density gradient at simulations for additional seven days beyond impact. As in
the interface between the planet and the disk. To analyse previous simulations the Balsara switch (Balsara 1995) was
disk properties we define it as the bound orbiting material used to reduce unwanted dissipation in the disk. After 2.5
thathasadensitylowerthatthereferencedensityρ atthe daysthedifferencesinthedisk’sstructureduetothemethod
0
surface of the planet. Despite visible differences, for ISPH aresignificant.ThediskobtainedfromISPHisrotatingsig-
MNRAS000,1–??(2016)
10 C. Reinhardt and J. Stadel
Figure 10. Results from a b = 0.71, v = 10km/s impact of a
Mars-sized (M = 0.099M⊕) impactor onto an Earth-mass tar-
get (M = 0.99M⊕) shown in Figure 9 for different resolutions
of 105 (green), 8·105 (blue) and 6.4·106 (red) particles. The
top plot shows the internal energy vs density of the particles in
the ”arm”like structure that forms between the target and the
surviving part of the impactor. For increasing resolution it be-
comes denser and more particles have an internal energy larger
than uIV meaning that the material is partially evaporated. So
thephysicalcomposition(andthusthepressure)inthislowden-
Figure 11.Thecircum–planetarydiskresultingfromab=0.5,
sity region clearly depends on the resolution. The bottom plot
shows the same for the impactor remnant that forms the end v = 10.61km/s impact between a M = 0.9M⊕ target and an
of the ”arm”and consist mostly of impactor material that was M = 0.2M⊕ impactor about three days after the collision. The
colour shows the materials temperature from 4500K (blue) to
not dispersed during the impact. The light orange region marks
40’000K(red).Thetopframe(a)showstheresultforclassicSPH
where we did the pressure cut–off because it would be negative
withoutanyofthemodificationspresentedinthispaperandthe
fortheTillotsonEOS.Inthisregionthematerialissupposedto
bottom frame (b) shows the outcome when ISPH is used. For
formsmalldropletsandwouldideallybetreatedbyamulti–phase
classicSPHthediskcooler,slightlymoreexpandedandlessflat.
fluidprescriptionwithinthehydrodynamics.
nificantly faster (Figure 12) and is less extended than the much better result and introduces less noise compared to
oneresultingintheclassicSPHsimulation.Wealsoobserve previouswork.Applyingthisdensitycorrectiontoourequi-
that for ISPH the disk remains flatter and is slightlyhotter librium models results in stable initial conditions that do
(Figure 11). These differences remain until the simulations not have to be relaxed prior to the simulation, significantly
ends eight days after the impact. At this point we leave a reducing the total simulation time allowing us to reach so
detailed investigation to later work. farunprecedentednumberofparticles.Finallywepresentan
entropyconservingformulationofSPHthatdoesnotrelyon
an analytic or thermodynamically complete EOS. The evo-
lutionoftheinternalenergyissplitintoacontributionfrom
6 SUMMARY AND CONCLUSION
shockheatingusingthestandardSPHartificialviscosityand
In this paper we present three solutions to long standing an exact adiabatic evolution that is calculated by interpo-
problems of GI simulations with the SPH method. First we lating between isentropes. For standard SPH the particles
introduce a method to generate low noise particle repre- significantly deviate from their original isentrope. Our new
sentation of planets. Comparing the particle’s (smoothed) method conserves entropy and the particles never fall be-
density to the numerical solution of the equilibrium model low their adiabatic evolution curves. At the same time, the
shows excellent agreement and a low spread around the de- Rankine-Hugoniotjumpconditionsforthefluidvariablesat
sired value. We then propose an improved treatment of a the discontinuity are successfully recovered, so our scheme
material/vacuuminterfaceasitoccursataplanet’ssurface. properly captures shocks as expected.
There the classic SPH density estimator fails, calculating a AcomparisonbetweenthenewmethodandclassicSPH
significantly lower value. As a result the models obtained when applied to GI simulations shows that the methods
withballicdeviatefromequilibriumandstarttooscillate. mostly agree on the general outcome within the first 24
Thisisevenmoreproblematicwhenthesurfacetemperature hours after the impact. In one simulation (Simulation 2 in
ishighasinthiscasethematerialispartiallyvaporized.As Section5)theimpactproducesaflattened,rotationallysup-
a solution we propose a modification of the standard SPH ported disk of hot, low density material around the target.
densityestimatorthataccountsforonlythepartoftheker- While the total bound mass is very similar, ISPH produces
nel volume that is sampled with particles, which gives a a slightly more massive central body and less massive disk
MNRAS000,1–??(2016)