Giant impact simulations

[maxpkatz]

In this thread we will discuss using Castro for giant impact (GI) simulations, i.e. collisions between protoplanets, planetesimals, and similar objects.

[maxpkatz]

Living literature review:

The giant impact (GI) hypothesis has been invoked to explain several aspects of the solar system, including the formation of Earth's moon (Hartmann and Davis (1975); Cameron and Ward (1976)), the retrograde rotation of Uranus (Safronov (1966)), and the high density of Mercury via stripping of its mantle in a collision (Benz, Slattery, and Cameron (1988); Benz et al. (2007)).

Simulations of GI scenarios have been done since the 1980s, starting with Benz, Slattery, and Cameron (1986) investigating the Moon formation question using the SPH method. Most simulations to date have been performed with SPH, including Canup and Asphaug (2001), Reinhart and Stadel (2017), Reinhard et al. (2020), and Meier, Reinhardt, and Stadel (2021). Typical SPH particle counts are ~10^6 in recent literature, but some recent work (Genda et al. (2015)) noted that simulations may be unconverged at low particle counts and Kegerreis et al. (2019) used ~10^8.

A few studies have used grid-based codes, particularly the Sandia code CTH: Melosh and Kipp (1989), Leinhardt and Stewart (2009), and Canup, Barr, and Crawford (2013). Additionally, Wada, Kokubo, and Makino (2006) used an unnamed code that uses the AUSM method, and Mullen and Gammie (2020) used Athena++. Typical spatial resolutions are 100-200 km in this literature.

The composition profiles used in these simulations are fairly simple. For collisions among rocky objects, the typical choice in recent literature is a two-component system where the interior of the object is iron (1/3 by mass) and the exterior is a silicate material like granite or dunite (2/3 by mass). Some papers just used a uniform composition of granite for the whole object. For the Uranus case, the typical assumption is a rocky core, an ice mantle, and a hydrogen-helium atmosphere (see e.g. Kegerreis et al. (2018)).

The physics relevant for simulating these events includes: hydrodynamics, self-gravity, and an equation of state. Even though the target and impactor are initially partially or fully solid, these simulations generally neglect viscosity, so we can treat the objects as fluids, at the cost of being unable to represent permanent effects on the target such as an impact crater (the force of self-gravity will eventually cause the object to become spherical again in the absence of any shear strength). Self-gravity is important for understanding the long-term evolution of the system, in particular for understanding questions like what happens to an accretion disk formed by the impactor, although simulations focusing only on the time period surrounding the initial collision can possibly get away with simpler approximations.

Two equations of state are in common use in the literature. While these are not the only ones used, our initial investigation into this would probably stick with well-studied physics as we validate our implementation, and then we could discuss more recent work (e.g. Wissing and Hobbs (2020)).

The first common EOS is the Tillotson EOS and was described in Tillotson (1962), with a useful summary in Appendix A of Reinhart and Stadel (2017). This EOS breaks down into two general regimes, one where material has been compressed relative to a reference density and one where material has expanded relative to that reference density. The compressed regime is of the form of a Mie-Gruneisen equation of state, which describes the thermodynamic properties of a shocked material by calculating the R-H jump conditions relative to to the reference state (density and internal energy). Tillotson fit the numerical parameters of the model in the low-pressure limit to experimental data and in the high-pressure limit (pressures above what could be achieved in lab experiments) to calculations from the Thomas-Fermi model. In the expanded regime there are three cases depending on the internal energy. For sufficiently low internal energy, none of the material has vaporized so the pressure formula is the same as the compressed state. For sufficiently high internal energy, all of the material has vaporized, and a different pressure formula is given, corresponding to an ideal gas plus a correction term. In between these internal energy thresholds, which is where there would be a mixed solid/liquid and vapor phase, Tillotson interpolates the pressure between the cold and hot limits, which is commonly cited as a drawback of the model because it is physically motivated. Some notes about this EOS are warranted. First, it is fully analytic, but Tillotson only provides pressure given (rho, e). Other derived quantities such as sound speed (dP / drho), gamma (d ln P / d ln rho), and dP/de need to be derived from the analytic form. Reinhardt and Stadel give the derivation of sound speed in Appendix B, it should be straightforward to determine the others (and they also have a repo providing their implementation). Second, temperature is not involved in this EOS, we can determine pressure entirely from internal energy. The lack of temperature is not a serious problem as long as no chemistry is involved, but if we do want an approximate temperature, we can use the construction in Reinhardt's repo where we get T from u(T) = u(0) + c_v * T. This specific heat is usually provided in the relevant papers. Third, the Tillotson EOS is parameterized for each material. For SPH simulations you can assign a material type to each particle and just evaluate the EOS accordingly for each particle. But for a grid code with a continuous composition distribution, a given cell is likely to have some mix of materials in it, and it is not obvious how to compose the EOS for multiple materials in the same grid element (at least, I haven't seen any papers explicitly address this, even the grid code papers). Since these are not ideal gases a linear weighting by mass fraction is not obviously correct. Fourth, the pressure can be negative in the analytic formula. This is intentional because it represents a surface tension, however it obviously causes problems for a hydro code, so most papers set the pressure to zero or a floor when that happens, assuming that gravity will hold the object together.

The second commonly used EOS is called ANEOS, which is also an analytic equation of state (Thompson and Lauson (1972)). ANEOS is much closer to the kind of EOS we're familiar with: it constructs the Helmholtz free energy as the sum of three terms: F from the (zero-temperature limit of the) atomic lattice, F from the thermal vibrations of the atoms, and F from the ionization of electrons. The various thermodynamic terms are constructed as derivatives of F as usual. If we come in with (rho, e) then we do a Newton-Raphson iteration to get T. Several improvements to ANEOS have been suggested over time, particularly the M-ANEOS variant proposed in Melosh (2007) and later work. M-ANEOS points out that the original ANEOS was constructed for metals, not minerals, and in particular does not get the vapor phase correct for materials like SiO2 because it neglects the presence of molecules in the vapor (the energy of the vapor is too high so there's a higher barrier to forming it through vaporization in a simulation). While ANEOS and the later variants are too complicated to implement directly based on the description in the paper, M-ANEOS is available on GitHub in a Fortran code that is permissively licensed.

There is a fair bit of discussion on the initial conditions for these simulations, which breaks down along two lines. First is the construction of equilibrium profiles of the individual objects (target and impactor). As far as I can tell, this is very similar to the problem we have in astro for constructing equilibrium profiles. We may have the same problem we have in wdmerger which is that it's hard to actually damp something out to a stable long term equilibrium, and we will definitely have the same problem that it's hard to maintain that equilibrium when the object is advecting across the grid. It seems that the reason so much discussion goes into this is for the general reason that it's non-trivial to construct equilibrium distributions of particles in SPH simulations, and most of these papers use SPH. But constructing the 1D profiles to sample from is still the HSE problem we are familiar with. The second concern is the velocity and relative orientation of the objects. The angular momentum of the orbit of the Moon about the Earth is actually higher than the rotational angular momentum of the Earth, which is not true for other planets with moons in the solar system. This fact is one of the interesting facts about the system and is one reason why a giant impact scenario is favored, so it means that for Moon-forming simulations getting the angular momentum of the final configuration correct is important. Typically this is done by defining an impact parameter and a velocity at infinity; Meier, Reinhardt and Stadel (2021) has useful discussion (Section 3 and Appendix A) specifying the exact formulas for constructing the initial locations and velocities of the two objects.