Self similar Shocks in Atmospheric Mass Loss due to Planetary Collisions

We present a mathematical model for the propagation of the shock waves that occur during planetary collisions. Such collisions are thought to occur during the formation of terrestrial planets, and they have the potential to erode the planet's atmosphere. We show that under certain assumptions, this evolution of the shock wave can be determined using the method of self similar solutions. This self similar solution is of type II, which means that it only applies to a finite region behind the shock front. This region is bounded by the shock front and the sonic point. Energy and matter continuously flow through the sonic point, so that energy in the self similar region is not conserved, as is the case for type I solutions. Instead, the evolution of the shock wave is determined by boundary conditions at the shock front and at the sonic point. We show how the evolution can be determined for different equations of state, allowing these results to be readily used to calculate the atmospheric mass loss from planetary cores made of different materials.


Introduction
The formation of terrestrial planets is thought to proceed in two phases. In the first phase, rocky cores form by accretion in the protoplanetary disc [3], and in the second stage, multiple planetary cores merge to form planets after the disc evaporates [2,6].
Theoretical models predict that when the disc evaporates, it will leave behind a large number of tightly packed planetary cores. Dynamical instabilities will excite eccentricities which will cause the orbits to intersect and the core to collide [7,33,18,21]. Among other things, these collisions can erode a planetary atmosphere.
It is difficult to observe planetary collisions directly, as they are expected, in the most optimistic scenario, to produce a faint, short (few hours) and hard (X-ray and EUV) transient [32]. Study of this phase is therefore based on theoretical or computational models, and indirect observational evidence, such as isotope ratios on Earth [23] or debris rings around exoplanetary systems [27]. Planetary collisions have been suggested as possible mechanism to explain systems like Kepler 107 [5], where a pair of neighbouring planets have very different average density. Kepler 107b has a low density, indicating the presence of an atmosphere, while Kepler 107c, while farther away from the host star, has a density almost three times larger than 107b, indicating a bare rocky core devoid of an atmosphere. If both planets started out with a thick atmosphere, other mechanisms, such as photoevaporation, could not have stripped Kepler 107c of its atmosphere while sparing 107b.
Simulations of giant collisions have been carried out by many authors in the past, e.g. [25]. However, these authors mostly considered the change in solid mass. For this purpose, authors ran many simulations, spanning a wide range of parameters. In contrast, for simulations of atmospheric mass loss due to giant collisions, most works focused on a single specific case, e.g. [14] which considered Kepler 36.
A systematic study of atmospheric mass loss for head on collisions has been carried out in [22,30]. The model they developed consists of two phases. In the first phase a shock wave emerges from the impact site and propagates through the planet. The second stage focuses on the interface between the ground and the atmospheric, and notice that when the ground moves due to the shock propagating through the target, another shock wave emerges from the gasground interface and moves through the atmosphere. Since the atmosphere has a declining density profile, the shock waves that travels through it accelerates as it moves up. At some point the shock wave exceeds the escape velocity, and all the gas above this point is lost. For this reason, some of the gas column above the ground is expelled even if the ground move slower than the escape velocity.
Propagation of a shock due to an impact event has been studied using simulations [20,17] and experiments [11,15,10]. However, under certain assumptions, the shock can be described analytically, using the theory of self similar solutions [4]. The most famous solution of this kind is the Sedov Taylor solution for the strong adiabatic explosion [24,26]. In this case, the motion of the shock front is determined using conservation of energy. Such cases where the self similar motion can be determined using conservation laws is called type I solutions. In contrast, the shock waves in planetary collisions are type II solution. In such cases, conservation laws cannot be used and the motion is determined by the behaviour near a singularity [28].
The intention of this paper is to explain how the evolution of the shock in the interior can be described using the theory of self similarity. The structure of this paper is as follows: in section 2 we describe the self similar solution of the shock propagating in the core. In section 3 we demonstrate how the atmospheric mass loss can be estimated. Finally. we discuss the results in section 4.
2 Core Shock

Impulsive Piston Problem
Let us consider an impact on a very large planet, so that the radius of curvature is unimportant and the ground can be considered flat. Such an impact will launch a shock wave that spreads as a hemisphere below the ground. When the radius of the shock wave is much larger than the size of the impactor, there are no other length scale in this problem and the shock radius grows as some power law in time R ∝ t α . The purpose of this section is to estimate α.
Even if we assume the shock to have azimuthal symmetry, the equations of motion still depend on two spatial coordinates and time. We therefore turn to an even simpler problem -the impulsive piston problem [1,31]. In this case, we consider a slab symmetric, one dimensional analogue to the self similar impact problem. In this scenario, a thin wafer hits a much thicker slab of material and both are perfectly cold prior to the collision. As a result of this collision a shock wave emerges from the contact surface and travels into the target. Once the shock wave has travelled a distance much larger than the width of the wafer, there is no other relevant length scale in this process other than the position of the shock front, and hence we expect the shock wave to evolve in a self similar way.
What's incredible about the impulsive piston problem is that the relation between the shock velocity and the swept up mass v ∝ m δ holds also in the three dimensional impact problem [30], i.e. a power law with the same exponent δ. In the remainder of this section we present the mathematical formulation of the impulsive piston problem.

Self Similar Equations
The slab symmetric hydrodynamic equations in one dimension are given by [13] ∂ρ ∂t We assume an ideal gas equation of state Naturally, one would expect the material in a planetary core to have a more complicated equation of state. However, we argue that at very high pressures and temperatures all materials behave like ideal gases. We can replace the pressure by the sound speed using the relation of p = ρc 2 /γ. After this substitution we have a system of partial differential equations in v, c, ρ. We make the assumption that one dimensional shock has a self similar behaviour, so the position of the shock front X(t) evolves as a power law in time where α is some number which will determined later. We can use the conservation of energy or momentum to obtain limits on this parameter 2 3 > α > 1 2 [31]. The partial differential equations can be reduced to ordinary differential equations in the dimensionless position χ = x/X. We define dimensionless hydrodynamical variables V, C, D in the following way To get rid of d 2 X(t) dt 2 terms, we make the substitution: . Substituting equation 5 into this gives a relation between δ and α: δ = 1 − 1 α . The energy and momentum conservation limits on this parameter are − 1 2 > δ > −1. After making these adjustments and simplifying, the hydrodynamic equations equations 1, 2 and 3 become where a prime denotes derivation with respect to χ. These equations are now linear in V , C , D so these can be isolated. The relevant equations are V , C , and are given by where Dividing one equation by the other, we get a single ODE where: and

Boundary Conditions
In order integrate equation 14, boundary conditions and the value of α are needed. The boundary conditions at the shock front are given by the Rankine Hugoniot conditions for a strong shock [13,31] V f = 2 γ + 1 (17) and In some self similar problems, like the Sedov Taylor explosion, the parameter α can be be inferred directly from conservation laws. Such problems are known as type I solutions. In our case, however, α is determined by the condition that the hydrodynamic trajectory passes smoothly through some singularity. These are called type II solutions [28]. In such cases the self similar solution only applies to a finite portion of space, bounded by the shock front and the singularity. This singularity occurs at a point where information cannot propagate back to the shock front, and is therefore also referred to as the sonic point. Matter, momentum and energy continuously flow through the sonic point, so that they are not conserved in the self similar region. Equation 10 has a singularity when We call this curve the sonic line. On the sonic line, the denominator (equation 16) vanishes. To prevent a divergence, the numerator also has to vanish. This happens on a specific point on the sonic line, which we call the sonic point and Knowing the boundary conditions and dC dV , we can integrate from the sonic point to the shock front, as shown in Figure 1. The starting point is slightly shifted from the sonic point in the positive V direction by dV i , and in C direction by dV i dC dV sonic for some small dV i 1. The slope at the sonic point is given by the equation: The correct value of δ is such that the curve C (V ) satisfies both boundary conditions at the sonic point and the shock front. In practice, we use the shooting method to find the value of δ. We guess a value for δ, numerically integrate equation 14 w.r.t to V from the sonic point to the shock front, and note the value of C at the shock front. We then use the bisection method to refine the value of δ to minimise the distance between the value of C obtained from numeric integration and its theoretical value (equation 18). An example for some hydrodynamic C (V ) trajectories for the same value of γ but different values of δ can be seen in figure 1, for γ = 2.8.
In the case γ > 2, numerical integration is straightforward since both the sonic point and the shock front are on the same side of the impact site (i.e. both x > 0 when the impact site is x = 0). The case 2 > γ > 1 is more complicated because the integration goes through x = 0. Because of the way we defined the self similar variables (equation 6), they diverge at x = 0, even though the physical quantities v, c remain finite. We can circumvent this difficulty by noting that the Mach number, i.e. the ratio between the velocity and the speed of sound, remains finite and changes continuously across x = 0. The sonic point in the case 2 > γ > 1 lies in the fourth quadrant of the of the C − V plane (i.e. positive V and negative C). Numerical integration proceeds to even higher values of V , moving away from the shock front. Far away from the sonic point, the curve attains some asymptotic slope, and this slope is the Mach number at x = 0, which we'll denote by M 0 . Once we have this piece of information, we can restart the numerical integration in the second quadrant (C is positive and V is negative), beginning from some arbitrarily highly negative V = V r (subscript r for restart) and C r = V r /M 0 . From this point we can continue the numerical integration all the way to the shock front. An example for several hydrodynamic trajectories for γ = 7/5 can be seen in figure 2.
Thus, for every value of γ we can obtain the corresponding value of δ. The relation between them is shown in figure 3. These results are consistent with previous works [1]. We note that in the limit γ = 1, we get δ = −1, which corresponds to momentum conservation. However, in the limit γ → ∞, the solution does not converge to δ = −1/2, which corresponds to energy conservation, but to a different value. More on this in the next section.
We note that the results obtained in this section are in accord with the numerical simulations presented in [30], which found that for the three dimensional impact problem, the shock decelerates roughly asẊ ∝ m −2/3 (where m is the swept up mass) for γ = 5/3, where the model presented here predicts δ = −0.64. The model here also explains impact experiments, although the translation between the theoretical results and the experimental results is not straightforward. In experiments, like those described in [9], a projectile is fired at a slab of some material and the radius of the resulting crater is measured. In materials without shear strength, like fluids, the crater basin stops growing when the pressure inside the crater is comparable to the hydrostatic pressure. The analysis above gives us a relation between the shock velocity and the swept up massṘ ∝ m δ , where m is the swept up mass, and in a three dimensional case it scales as a cube of the radius m ≈ ρR 3 , where, for simplicity, ρ is the initial mass density of both projectile and target. If we denote the radius of the projectile by a, and we assume that the shock sweeps a mass comparable to the projectile before decelerating, then the relation between the shock velocity and the crater radius is roughly given bẏ where U is the initial projectile velocity. The pressure inside the crater is roughly given by ρṘ 2 , while the hydrostatic pressure is given by ρgR, where g is the gravitational acceleration. Equating the two pressures yields a relation between the volume of the crater R 3 c (where R c is the terminal size of the crater) and the properties of the projectile where we define the same dimensionless variables π v = R 3 c /a, π 2 = ga/U 2 as [9]. To produce a prediction, we need to estimate the effective adiabatic index. This can be done in different kinds of high pressure experiments, where a material is shocked and the material U m and shock velocity U s are measured. At high velocities, the ratio between the two tends to a constant, U s /U m = β, and for an ideal gas β = (γ + 1) /2, so γ = 2β −1. For water, equation of state experiments yields β = 1.78 [8], so the effective adiabatic index γ = 2.56 our model yields δ = −0.61 and d ln π v / ln π 2 = −0.643. The value obtained from experiments in water is d ln π v /d ln π 2 = −0.648 [9]. Overall, despite our crude approximations, we obtain a value relatively close to the experimental value.

Asymptotic Case
In this section we consider the case γ → ∞. In this case we cannot use the equations developed above. Part of the reason is that the integration domain V ∈ − 2 γ−2 , 2 γ+1 shrinks to zero width in this limit. To overcome this difficulty, we define a new variable W = γV . In this new variable, the integration domain   becomes W ∈ [−2, 2], and so remains finite in the limit γ → ∞. The expression for the derivative in the new variable is given by At the sonic point, the slope is given by .
We can determine δ using the shooting method described in the previous section. This case is actually simpler, because it does not depend on γ, so C = 1 at the sonic point and C = √ 2 at the shock front. A few hydrodynamic trajectories for different values of δ are shown in figure 4. Using numerical root finding methods we find that lim γ→∞ δ ≈ −0.557, in accordance with [1]. Interestingly, this result is different from the energy conservation limit δ = −1/2.

Analytic Case
For most values of γ, numerical integration has to be performed to obtain the power law index δ. However, for γ = 7 5 , it is possible to obtain a completely analytic solution [1]. Here we will briefly describe the analytic solution and the hydrodynamic profiles in the case γ = 7 5 . The dimensionless hydrodynamic quantities are We can see that at the shock front, where χ = 1, the dimensionless velocity and speed of sound are V = 5 6 , C = √ 7 6 , as verified by the coordinates of the black X in second graph of Figure 2. We also note that V, C diverge as χ → 0. From the shooting method, we obtain δ = − 2 3 , which is consistent with the results in Figure 2.

Atmospheric Mass Loss
In this section we demonstrate how the atmospheric mass loss from a giant collision is estimated. The dependence of the mass loss on impactor size and velocity in head on collisions was studied in [30], so in this section we consider the effects of obliquity. For this purpose, we run simulations of collisions between planetary cores and planets. These simulations were run using the moving mesh numerical simulation RICH [29]. We use an ideal gas equation of state with an adiabatic index γ = 5/3. The planet and planetary core are described as blobs with uniform density. Since the simulation cannot handle actual vacuum, we fill the computational domain with a very tenuous gas, whose density is lower than that of the planet by a factor of 10 9 . We do not include gravity, and the entire domain has an initially uniform pressure such that the speed of sound in the planet is 10 −6 of the impact velocity. Lastly, to simplify the calculation, we perform these simulations in a 2D Cartesian geometry. This means that the planet and planetary cores are not represented by spheres, but by cylinders. As opposed to the other simplifications, this last one cannot be justified, and might lead to a large discrepancy with other results. The only reason we employ this assumption is that it simplifies the analysis greatly, and since the purpose of this section is mostly pedagogical. Figure 5 shows four log density snapshots from an oblique collision, where the radius of the impactor is 0.1 times the radius of the planet. The figure shows a shock wave emanating from the impact site that travels through the planet, and obliterates it, since there is not gravity to hold it together.
As the shock moves through the planet, it moves the ground, which pushes against the atmospheric gas column above it. In this sense, the ground acts like a piston, which sends a shock that moves up and away from the ground, accelerating due to the declining atmospheric mass density. The mass loss from each gas column as a function of ground velocity has been worked out in [22]. By integrating the local atmospheric loss across the entire surface we can obtain the total atmospheric mass loss. By running the simulation with different impact Figure 5: Snapshots from a simulation at different times of a collision between a planet and a planetary core whose radius is a factor of 0.1 smaller. We do not include gravity, so the planet is obliterated in the end. The simulation is performed using a moving mesh hydrodynamic simulation, so the polygons are actual computational cells.
parameter (or offset) and different impactor to target radius ratios (but fixed velocity, which we set to be 1.5 times the escape velocity), we obtain a map of atmospheric mass loss shown in Figure 6. One of the interesting features of this map is that for small impactors, the offset has very little influence on the outcome. This is because small enough impactors deposit all their energy in the target regardless of obliquity. This trend is in agreement with the findings of [30], where it was shown that when the shock radius is considerably larger than the size of the impactor, the shock wave tends to the self similar solution found in the previous section, and in the process loses some of the information about the initial conditions (i.e. fast and small impactors create the same shock wave as slow and big impactors). The loss of information about the initial conditions is a common consequence of self similar solution [4]. Figure 6: A map of the atmospheric mass loss as a function of the radius ratio between the impactor and target, and the offset between the centres normal to the initial velocity difference, or impact parameter. The velocity is held fixed at 1.5 of the escape velocity.

Conclusion
Planet formation involves a phase of giant collisions, when planets and planetary cores smash into each other. These collisions erode the planet's atmosphere. A key ingredient in modelling this mass loss is the evolution of the resulting shock wave in the interior of the planet.
In order to justify the use of self similar solutions, we had to make a number of simplifying assumptions. First, we had to neglect all dimensional parameters in the problem, so such solutions only apply for shock radii much larger than the impactor on the one hand, but much smaller than the planet on the other hand. Second, we had to assume a strong shock, so that the impact velocity has to be much larger than the escape velocity. Third, we used an ideal gas equation of state, whereas real materials have a more complicated equation of state, which includes phase transitions. We argue that at high pressures and velocities, all materials behave like an ideal gases. Moreover, high pressure experiments can be used to determine the effective adiabatic index for different materials [30].
Even with all these assumptions, the problem is still too complicated to be solved analytically, so we simplify it even further by considering the one dimensional analogue -the impulsive piston problem. This problem is a type II self similar shocks, where the evolution of the shock wave is not determined by some conservation law, but according to a singularity behind the shock front. By requiring that the hydrodynamic profiles satisfy both conditions at the shock front and the singularity, we can obtain δ = d lnẊ/d ln m, whereẊ is the shock velocity and m is the swept up mass. It turns out that this relation holds even in different geometries, i.e. in the three dimensional impact.
The impulsive piston problem has been originally solved by [1]. That work, however, used a different set of dimensionless variables, and so they could not reduce all equations to a single ordinary differential equation. They therefore had to use a multidimensional shooting method to determine δ. The variables we chose allow us to reduce the problem to a single ordinary differential equation, but it comes at a cost. The cost is that we introduce a coordinate singularity, i.e. a place where the dimensional variables remain finite, but the dimensionless variables diverge. The price that we have to pay for this choice is that for some range of values for the adiabatic index 2 > γ > 1 the numerical integration is not continuous.
It is surprising that despite the many simplifying assumptions made here, the amounts of atmospheric mass loss predicted from this model is similar to more complicated models, which take into account all the effects considered here [30]. In the future, it would be interesting to see if this model can be refined to include some of the effects previously neglected.