# Impact-induced Vaporization during Accretion of Planetary Bodies

Adrien Saurety, Razvan Caracas, Sean N. Raymond

► To cite this version:

Adrien Saurety, Razvan Caracas, Sean N. Raymond. Impact-induced Vaporization during Accretion of Planetary Bodies. The Astrophysical Journal Letters, 2025, 981, 10.3847/2041-8213/adb30e. insu-05080929

HAL Id: insu-05080929

https://insu.hal.science/insu-05080929v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Impact-induced Vaporization during Accretion of Planetary Bodies

Adrien Saurety¹, Razvan Caracas¹,², and Sean N. Raymond³

¹Institut de Physique du Globe de Paris, 1 rue Jussieu, Paris, 75005, France

²The Center for Planetary Habitability (PHAB), University of Oslo, Oslo, 0371, Norway

³Laboratoire d'Astrophysique de Bordeaux, Université de Bordeaux, CNRS B18N, Pessac, 33615, France

Received 2024 November 8; revised 2025 February 5; accepted 2025 February 5; published 2025 February 26

## Abstract

Giant impacts dominate the late stages of accretion of rocky planets. They contribute to the heating, melting, and sometimes vaporizing of the bodies involved in the impacts. Due to fractionation during melting and vaporization, planet-building impacts can significantly change the composition and geochemical signatures of rocky objects. Using first-principles molecular dynamics simulations, we analyze the shock behavior of complex realistic silicate systems, representative of both rocky bodies. We introduce a novel criterion for vapor formation that uses entropy calculations to determine the minimum impact velocity required to pass the threshold for vapor production. We derive impact velocity criteria for vapor formation—7.1 km s⁻¹ for chondritic bodies—and show that this threshold is reached in 61% and 89% of impacts in dynamical simulations of the late stages of accretion with classical and annulus starting configuration (respectively) for analogs of Earth. These outcomes should be nuanced by factors such as the impact angle and the mass of the impacting bodies, which further influence the vaporization dynamics and the resultant material distribution. Our findings indicate that vaporization was common during accretion and likely played a crucial role in shaping the early environments and material properties of terrestrial planets.

*Unified Astronomy Thesaurus concepts: Impact phenomena (779); Planet formation (1241); Mineral physics (2230); Carbonaceous chondrites (200); Solar system formation (1530); Planetesimals (1259)*

## 1. Introduction

The impacts are a central process in planet formation. During the growth of terrestrial planets, some models suggest that the impacts became increasingly energetic as 10–100 km scale planetesimals accreted into planetary embryos that were ~Mars-sized or larger (E. Kokubo & S. Ida 2000). Other models favor a first phase of pebble accretion followed by impacts between large planetary embryos at the end of accretion (A. Johansen et al. 2021). In any case, the final phase of terrestrial accretion was dominated by giant impacts between planetary embryos (S. N. Raymond et al. 2009; A. Morbidelli et al. 2012), culminating with the Moon-forming impact (R. M. Canup & E. Asphaug 2001; M. Ćuk & S. T. Stewart 2012). While this phase has been well studied with N-body simulations (C. B. Agnor et al. 1999; J. E. Chambers 2001; S. N. Raymond et al. 2009; A. Morbidelli et al. 2012; R. A. Fischer & F. J. Ciesla 2014; A. Izidoro et al. 2015; E. V. Quintana & T. Barclay 2016; M. S. Clement et al. 2018), including some studies that include the effect of collisional fragmentation (J. E. Chambers 2013; M. S. Clement et al. 2019; A. Emsenhuber et al. 2020), such simulations cannot currently capture the full, complex outcomes of giant impacts (H. Genda et al. 2012; Z. M. Leinhardt & S. T. Stewart 2012).

The energies reached during giant impacts could sometimes be sufficiently high to melt or even partially vaporize the impactor and target (E. J. Davies et al. 2020). Vaporization would enhance the loss of volatiles in space due to ionization and may significantly affect the geochemical composition of rocky objects in terms of volatile and moderately volatile elements and light isotopes (R. C. Hin et al. 2017; C. A. Norris & B. J. Wood 2017).

Elemental and isotopic fractionation during partial vaporization might explain specific geochemical differences between chondrites and differentiated planets (H. S. O'Neill & H. Palme 2008; E. A. Pringle & F. Moynier 2017). It is, therefore, of great importance to estimate whether impacts leading to partial vaporization are frequent events during terrestrial accretion or whether they are rare and have a limited effect on the global geochemical evolution of rocky objects.

Shocks in geological materials have been studied extensively for decades. Using shock experiments, ab initio molecular dynamics (AIMD) and thermodynamic models such as ANEOS (H. J. Melosh 2007) or SESAME (C. E. Ragan 1982), the community proposed the shock equation of the state of several minerals and compounds such as forsterite (T. J. Ahrens & J. D. O'Keefe 1972; S. Root et al. 2018), MgO (C. A. McCoy et al. 2019), feldspars (T. J. Ahrens et al. 1969; P. D. Asimow & T. J. Ahrens 2010), silica (H. J. Melosh 2007; R. G. Kraus et al. 2012), or pure iron (Z. Li et al. 2020). Shocks on complex realistic materials, such as dunite (R. G. McQueen et al. 1967), pyrolite (S. Stewart et al. 2020), granite (J. L. Shang et al. 2000), etc., are less common, and they were often focused on exploring melting relations and the thermodynamic properties of silicate melts (C. W. Thomas & P. D. Asimow 2013). The reason is that shock experiments are expensive and time-consuming, providing limited points, while ANEOS and SESAME models rely on strong assumptions. Here, we use AIMD simulations to study complex silicate's Hugoniot equation of state (HEOS). This method has already been applied for feldspar by A. Kobsch & R. Caracas (2020), forsterite by S. Root et al. (2018), and pure iron by Z. Li et al. (2020). We computed the HEOS and the entropy production during shock for realistic silicate systems representative of planetary bodies: the average composition of a CM2 chondrule (CMc; P. Friend et al. 2018), the average composition of a CM2 matrix (CMm; P. Friend et al. 2018), and the bulk silicate Earth (pyrolite; W. McDonough & S.-s. Sun 1995). We

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
obtain an entropy production due to impact that is consistent with the entropy of shocked forsterite obtained experimentally by E. J. Davies et al. (2020).

## 2. Methods

### 2.1. Materials

We use the Jbilet Winselwan CM2 carbonaceous chondrite (P. Friend et al. 2018) as a model for the chondritic compositions. Carbonaceous chondrites are primitive meteorites that carry material representative of the protoplanetary disk (R. H. Hewins 1997; J. Villeneuve et al. 2015). P. Friend et al. (2018) analyzed the composition of the matrix (CMm) and of several chondrules (CMc) using a microprobe. The chondrules are mostly type I (FeO-poor), which are dominant in carbonaceous chondrites (R. Hewins et al. 2005).

The terrestrial composition, or pyrolite, closely follows the bulk silicate Earth model (BSE; W. McDonough & S.-s. Sun 1995) that we used in several previous studies (R. Caracas & S. T. Stewart 2023; R. Caracas 2024). This composition reflects the bulk terrestrial magma ocean and is a good approximation for the composition of the homogenized protolunar disk.

All three model compositions used here are realistic representations of solar system bodies. The main relative differences are in the (Mg+Fe)/(Si+Al) ratio, which can influence the degree of polymerization in the melt. Table A1 lists the three compositions adopted here.

### 2.2. Ab Initio Simulations

We perform ab initio (Al) molecular dynamics (MD) simulations (R. Car & M. Parrinello 1985) using the Vienna Ab Initio Simulation Package software (G. Kresse & J. Hafner 1993; G. Kresse & J. Furthmüller 1996a, 1996b). The interatomic forces are computed using the projector-augmented wave function (PAW) formalism (P. E. Blöchl 1994) of the density functional theory (P. Hohenberg & W. Kohn 1964; W. Kohn & L. J. Sham 1965). We employ the Perdew-Burke-Emrzerhof functional for the generalized gradient approximation (J. P. Perdew et al. 1996). Atomic cores and core electrons are treated as pseudopotentials. The kinetic energy cut-off of the plane waves outside the PAW spheres is set to 550 eV; the cut-off for the Hänkel functions, which describe the augmented charge inside the PAW spheres, is set to 800 eV. The Brillouin zone is sampled at the Gamma point. The electronic population of the energy states is treated using a Fermi distribution with smearing proportional to the temperature ($k_B T$; N. D. Mermin 1965).

Simulations are performed in the canonical ensemble with a constant number of particles (N), volume (V), and temperature (T). The temperature is controlled with the Nosé-Hoover thermostat (S. Nosé 1984; W. G. Hoover 1985). The time steps are 1 femtosecond long, the simulations are thermalized for about 1 picosecond, and the production runs between 2 and 10 picoseconds. Those values allow us to have both a good convergence of the pressure, temperature, and energy and enough simulation steps to make good statistics and get representative measures for the thermodynamic properties we are interested in. The timestep of 1 femtosecond allows the atoms to move slowly enough to have a realistic behavior. Too large a timestep will lead the simulation to crash because atomic cores overlap, but fast enough to explore several different configurations in the phase space; too small a timestep will increase the simulation time needed to let the atoms move

enough to encounter different neighbors. Those values of timestep and simulation time are classic in AIMD simulation and have been used in several studies (A. Kobsch & R. Caracas 2020; Z. Li et al. 2020; R. Caracas & S. T. Stewart 2023; R. Caracas 2024).

We use the Universal Molecular Dynamic (UMD)/magnetomix software (R. Caracas et al. 2021) to perform the postprocessing of the AIMD simulations. The entropy is computed from the simulation using thermodynamic integration methods (see Appendix for more details).

### 2.3. Hugoniot–Rankine Equation Equations of State

Shocks are rapid processes in which materials are subjected to extreme compression and heating during short-impact events. Due to the rapid and intense conditions generated by shocks, these processes are inherently out of equilibrium, presenting significant temperature and pressure gradients that cannot be adequately captured by equilibrium-based equations of state, such as the Vinet or Birch–Murnaghan models, which assume a uniformly equilibrated system. Instead, the Hugoniot–Rankine equations, specifically designed to address these rapid, nonequilibrium changes, provide a more accurate description of the material behavior under such extreme conditions.

The Hugoniot–Rankine equations express the conservation of the energy in the shocked material before and after the shock wave passes (B. Y. Zel’dovich & Y. P. Raizer 2002; W. M. Isbell 2005; J. W. Forbes 2012; A. Kobsch & R. Caracas 2020):

$$E - E_0 + \frac{1}{2}(P - P_0)\left(\frac{1}{\rho} - \frac{1}{\rho_0}\right) = 0, \quad (1)$$

where, E, P, and ρ stand for the internal energy, pressure, and density of the material attained during shock, respectively, and the 0 index represents the initial state.

The HEOS is sensitive to the initial density and temperature of the material. Consequently, we select three representative densities spanning the range measured by D. Ostrowski & K. Bryson (2020) for the chondrule and matrix phase (see Table A1). For each material, the three densities correspond to the average density reported by D. Ostrowski & K. Bryson (2020) ($ρ_{0a}$), the average with the standard deviation added ($ρ_{0h}$), and the average with the standard deviation subtracted ($ρ_{0l}$). This approach allows us to test representative preshock densities for our material. We consider two preimpact temperatures, 1500 and 10 K (undercooled melt). The higher temperature is realistic of what can be found in a molten body of the protoplanetary disk (A. P. Boss 2004; A. Johansen & S. Okuzumi 2018; B. A. L. Gaches et al. 2020), and the lower temperature is a proxy for a static simulation of a glassy state. We consider only 3.1 g cm⁻³ initial density for pyrolite at 300 K as representative of the cold state of the magma ocean.

We perform AIMD simulations along several isotherms to find the density where P, E, and ρ satisfy the Hugoniot–Rankine equation. Throughout this work, we call the HEOS of a material the curve that verifies the Hugoniot equation at all temperatures and pressures.

The behavior of a material during shock can also be expressed as a function of particle velocity in the target, $u_p$, and the impactor velocity, $U_p$. In this formulation, the pressure in the target, P, is related to $u_p$ (B. Y. Zel’dovich &
Y. P. Raizer 2002; J. W. Forbes 2012):

$$P = a \times u_p^2 + b \times u_p + c. \quad (2)$$

We can obtain a series of $u_p$ values at different states along the HEOS from AIMD simulations, using the energy conversion (W. M. Isbell 2005):

$$u_p = \sqrt{2(E - E_0)}, \quad (3)$$

which, eventually, allows us to fit the a, b, and c parameters from Equation (2).

Assuming that we are in the planar shock approximation, we can determine the peak pressure (J. W. Forbes 2012) during the shock by looking for the point where

$$P(u_p) = P(u_{\text{impactor}}) \quad (4)$$

with

$$u_{p, \text{impactor}} = U_{\text{impactor}} - u_p. \quad (5)$$

This means that the impactor particles will share their momenta with the impacted particles. Figure A1 shows a real case of determination of the maximum pressure for a CMm target and a CMc impactor, the latter with a speed of 10 km s⁻¹.

## 2.4. N-body Simulation Models and Classification of Planets

We analyzed the results of 100 N-body simulations of the late phases of terrestrial accretion (for a discussion of different models of terrestrial planet formation, see S. N. Raymond & A. Morbidelli 2022). Of these, 40 (taken from S. N. Raymond et al. 2009) were in the context of the so-called “classical model” of terrestrial planet formation (J. E. Chambers 2001), and contained a smooth disk of planetesimals and planetary embryos extending from 0.5 au out to close to Jupiter’s orbit. The other 60 simulations (taken from S. N. Raymond & A. Izidoro 2017) were in the context of the “empty asteroid belt” or “annulus” model (B. M. S. Hansen 2009), in which planetesimals and planetary embryos were initially distributed in a narrow ring between 0.7 and 1–1.5 au.

Each impact is recorded during those simulations, and the two impacted bodies fuse. At the end of the simulations, only a few bodies are left. Those bodies are the planets formed by the accretion simulation. All bodies remaining at the end of the classical and annulus are shown in Figure A2.

We extracted the bodies with similar masses and orbits to the actual Venus, Earth, and Mars, i.e., bodies in the yellow, blue, or red rectangles in Figure A2. For each of them, we extracted the impacts that this planet and all the bodies stuck to the main body have experienced during the accretion. This gives us the complete shock history of the constitutive material of each final planet.

# 3. Results

## 3.1. Behavior of the Constituent Materials during Impacts

During large shocks (i.e., supersonic), the impacted materials are shocked from their initial state on their Hugoniots. The kinetic energy associated with the shock is largely transformed into entropy. If the shock is sufficiently energetic, the target and/or impactor can be projected up to molten or even supercritical states (R. G. Kraus et al. 2015; R. Caracas & S. T. Stewart 2023). After the passage of the shock wave, the thermodynamic trajectories of the bodies involved in the

impact are largely quasi-isentropic, as there is no heat exchange with the environment. Depending on the amount of entropy generated during impact, the production of vapor may occur during release, when the materials reach the liquid-vapor dome (LVD). This makes entropy the most critical parameter in describing the thermodynamics of impacts.

The topology of the HEOS is such that its trajectory is highly dependent on the initial state, i.e., $ρ_0$, $T_0$, $E_0$, and $P_0$. To quantify this influence, we independently vary the initial densities and temperatures (see Figure 1) and compute the HEOS for the chondrule composition. We find that for a difference in initial density of 0.3 g cm⁻³, the equation of state diverges by 0.2 g cm⁻³ by the time the pressure reaches 200 GPa. In terms of temperature, this difference of 0.3 g cm⁻³ in the initial density leads to a difference of 2000 K at 200 GPa, i.e., the peak temperature will be 6000 K for a chondrule with $ρ_0 = 3.65$ g cm⁻³ and 8000 K for a chondrule with $ρ_0 = 3.35$ g cm⁻³. This divergence shows the importance of considering a material with a relevant initial density: the same impact may form a melt and vapor mixture in the first case and a supercritical fluid in the latter. Increasing the temperature from $T_0 = 10$ K to $T_0 = 1500$ K at the same initial density results in an upward shift in density along the entire Hugoniot curve, which is roughly constant: 0.11 g cm⁻³ at 80 GPa and 0.09 g cm⁻³ at 180 GPa. In the temperature-pressure phase space, the $T_0$ effect seems less important when the pressure increases and the two Hugoniot curves tend to converge. But it remains nonnegligible at moderate pressure with a difference of 1000 K at 150 GPa between the $T_0 = 10$ K and $T_0 = 1500$ K curves.

As $ρ_0$ and $T_0$ clearly have a strong influence on the resulting HEOS, we also try to grasp the effect of the chemical composition. We compute the equations of state for CMc, CMm, and pyrolite for various initial conditions. For the sake of completeness, we compare the theoretical equations of state with the HEOS of several representative materials involving planetary impacts obtained from previous experimental or computational studies. Figure 2 shows the differences between all these Hugoniot curves. Although the topologies of the different curves are similar, there is a large scatter in both temperature and pressure. For example, at 150 GPa, the equation of state of the Na feldspar, a continental crust mineral, lies at the lowest temperature, only 3000 K, and the equation of state of CMm reaches beyond 4500 K. Similarly, at 4000 K, the pressure reached by Na feldspar is 200 GPa, whereas the pressure on CMm is only 130 GPa. Compared to the HEOS mentioned above, the ANEOS curve of dunite (S. Root et al. 2018) intersects the equation of state of feldspar, but is about 50 GPa above the calculated one of pure forsterite (S. Root et al. 2018). The two are above CMc, CMm, and pyrolite in terms of pressure. The scatter between the different curves can have multiple causes, some intrinsic, such as initial densities, temperatures, degree of polymerization of the materials, the role of cations, changes in magnetism, etc., and some extrinsic, such as errors in measurements and approximations in calculations. Figure 2 also shows that for the same initial conditions of temperature (10 K) and pressure (3.5 g cm⁻³), the HEOS of CMc and CMm are almost overlapping. In contrast, those of pyrolite and CMm for slightly different initial densities are very similar. This suggests that the silicate component of chondritic-like planetesimals can be treated as a homogeneous phase within smooth particle
<figure><img src="image_3.png" /><figcaption>Figure 1. Pressure variation along the Hugoniot shock equations of state as a function of density and temperature, for CMc ((a) and (b)) and for CMm ((c) and (d)) compositions with various initial states. Higher initial densities yield higher peak pressures at a given temperature. Higher initial temperatures yield higher peak temperatures at a given pressure. The effect of the temperature is more important for CMm than for CMc, probably because of the larger relative thermal-induced depolymerization of the CMm melt.</figcaption></figure>

hydrodynamics (SPH) simulations (R. M. Canup & E. Asphaug 2001), representing a considerable simplification in the path to acquire realistic large-scale impact simulations.

Apart from the initial state, one potential explanation for the scatter and similarities in the pressure and temperature of the different HEOS is the difference or similarity in chemistry. The ratio between the number of interstitial cations, such as Na, Ca, Mg, or Fe, and the number of network-forming cations, such as Si or Al, induces differences in the degree of silica polymerization between these materials (B. O. Mysen et al. 1982). The substantial presence of these structural polyhedra could profoundly influence the dynamic response of the material by strongly structuring the melt.

## 3.2. Entropy Production Due to Impact

Passing a shock wave through a material is associated with an increase in entropy. We can link the entropy and the temperature along the Hugoniot using (B. Y. Zel'dovich & Y. P. Raizer 2002; J. W. Forbes 2012) the equation

$$S = A \times \ln(T) + B, \qquad (6)$$

where A and B are constants. Although we cannot directly compute the entropy and free energy from the AIMD simulations, we can use thermodynamic integration to calculate these numbers. This method provides a direct accurate value for the entropy and does not require indirect measurements. It is, therefore, more reliable than interpolating indirect experimental data on heat capacity, which relies on many

approximations (E. J. Davies et al. 2020)—or by integrating an ANEOS model, which assumes a set of constant parameters and often gives results significantly different from the experiment (R. G. Kraus et al. 2015). Thus, this is the first direct estimate of the entropy produced during a shock and the entropy of a material on the LVD using AIMD simulation and direct entropy measurements.

We calculate the entropy along the HEOS of pyrolite and CMc, starting with $ρ_0 = 3.5 \text{ g cm}^{-3}$ and $T_0 = 10 \text{ K}$. The results are plotted in Figure 3. The entropy of these materials is relatively similar, up to about 6000 K. Below 5000 K, shocks in pyrolite yield higher entropy production than in CMc, and above 5000 K, the tendency is reversed. Under huge shocks, at temperatures far exceeding 6000 K, these two materials show significantly different slopes and tend to diverge. This difference might be explained by the difference between the HEOS of CMc and pyrolite (see Figure 2). Indeed, for the same peak temperature, the peak pressure is higher on the CMc Hugoniot than on the pyrolite Hugoniot. As entropy is a measure of chaos, a significant pressure that slows down the movement of the atoms will also affect the entropy and make it more difficult to increase.

Next, we compute the entropy of the pyrolite at three different temperatures (2000, 4000, and 6000 K) along the liquid spinodal curve from R. Caracas & S. T. Stewart (2023). Using these values and the critical temperature estimated in R. Caracas & S. T. Stewart (2023), we can interpolate the shape of the spinodal line in the T-S diagram on the liquid side of the LVD. The spinodal line marks the limit of metastability: at
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Various HEOS computed for representative planetary-building silicate compositions (solid lines). All our computed equations of state start at low temperatures, either 10 K or 300 K. In the high-temperature and high-pressure regime, the equations of state are quasi-linear; this relation is less obvious at low temperatures and low densities. The HEOS exhibit similar topologies, but their position is highly dependent on the initial state and, in particular, on the initial density. The chondritic and pyrolitic compositions (this study), with a (Mg+Fe)/Si ratio close to mantle values, produce temperatures much higher than the other phases considered here, partly due to their high initial densities.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Entropy along the HEOS for CMc and pyrolite compositions. The liquid-vapor dome (LVD) is an approximation following the temperature along the liquid spinodal line (R. Caracas & S. T. Stewart 2023). The LVD line is a second-order polynomial fit that ends close to the critical point (R. Caracas 2024). The release follows a quasi-isentropic trajectory in the aftermath of an impact, represented by the vertical dashed arrows. Vapor forms if the entropy produced during the shock is larger than the threshold for vaporization (equivalent to the entropy of boiling). Note that the temperature scale is logarithmic.</figcaption></figure>

entropies higher than the spinodal, the pure liquid becomes unstable, and vapor forms.

## 4. Shock-induced Vaporization by Impacts during Accretion

The absolute value of the entropy produced during a shock is useful as a thermodynamic parameter per se. But its importance becomes evident when considering the postshock state. In the

planar shock approximation, we can link the impact velocity of the impactor to the maximum peak pressure and temperature reached in the impacted body. Then we can estimate the entropy production during shock in the *T–S* diagram (Figure 3). The shock wave passage causes a rapid increase in pressure and temperature. After the shock wave passes, the material undergoes a decompression. This decompression is a slow process compared to compression and, as such, can be assumed to be an isentropic decompression (B. Y. Zel'dovich & Y. P. Raizer 2002; J. W. Forbes 2012; K. Kurosawa et al. 2013). Although this assumption has some limitations (P. G. Heighway et al. 2019), it gives a realistic approximation of the conditions encountered during shock release.

Within the isentropic release assumption, the postshock evolution of the material in the *T–S* diagram is a vertical line along which the material cools down at constant entropy. If this vertical line crosses the LVD, partial vaporization occurs. Thus, if the entropy produced during a shock exceeds the entropy of boiling (i.e., entropy of incipient vaporization), vapor forms when the system temperature intersects the temperature on the LVD dome. This is our criterion for identifying vapor as a result of an impact. As such, it is a more accurate measure for the incipience of vapor than comparisons to the entropy required to vaporize 50% of forsterite (as in, e.g., R. G. Kraus et al. 2012; E. J. Davies et al. 2020). We then use the LVD calculated by R. Caracas & S. T. Stewart (2023) to obtain the density of molten pyrolite close to the equivalent of a triple point, at 2000 K, conditions that are similar to the prediction of the ANEOS model for forsterite (E. J. Davies et al. 2020). We calculate the entropy of pyrolite at 2000 K close to the spinodal line and find a threshold value of $S_m = 3590 \text{ J kg}^{-1} \text{ K}^{-1}$ (see Figure 3) to trigger vaporization. This value is close to the triple point entropy for forsterite ($3500 \text{ J kg}^{-1} \text{ kg}^{-1}$; E. J. Davies et al. 2020). Our work using a direct measure of the entropy during AIMD simulations gives credence to the results estimated by E. J. Davies et al. (2020) and R. G. Kraus et al. (2012).

The analysis of the HEOS curves (Figure 3) shows that the entropy of boiling at 2000 K is achieved for shock temperatures above $T_{mc} = 3631 \text{ K}$ for CMc and above $T_{mp} = 3882 \text{ K}$ for pyrolite. These minimal peak temperatures can be related to a minimal peak pressure of $P_{mc} = 122 \text{ GPa}$ for CMc and $P_{mp} = 98 \text{ GPa}$ for pyrolite (Figure 2). Finally, we relate the position reached on the HEOS curve with the impact velocity ($U_{mc}$), and find that vapor forms for impacts of at least $U_{mc} = 7.1 \text{ km s}^{-1}$ for two chondritic bodies and $U_{mp} = 6.1 \text{ km s}^{-1}$ for two pyrolitic bodies (Figure A3), which is consistent with the $6 \text{ km s}^{-1}$ critical impact velocity proposed by E. J. Davies et al. (2020) for cold silica. These values are considerably lower than the requirement to vaporize iron cores, at $U_{mc} = 11.5 \text{ km s}^{-1}$ (Z. Li et al. 2020).

We compare this criterion with the impact velocities recorded during dynamical simulations of terrestrial planet accretion, following two different scenarios (S. N. Raymond et al. 2009; S. N. Raymond & A. Izidoro 2017). Figure A4 shows that the annulus model leads to higher impact velocities than the classical model. Higher collision speeds dramatically increase the proportion of impacts that lead to partial vaporization.

In the approximation of undifferentiated cold silicate bodies, with starting temperatures at 10 K, we find that about 61.5% of the impacts are energetic enough to produce vapor when impacting chondritic bodies in the classical accretion simulations. This
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Histogram of the distribution of impact velocities in *N*-body simulations with a classical starting configuration (panels (a), (c), and (e)) and annulus starting configuration (panels (b), (d), and (f)) for Venus, Earth, and Mars analogs. For the sake of clarity, the vertical scale is not the same on the left and right panels. The dashed lines indicate the critical impactor velocities required to have impact-induced vaporization for chondrule material.</figcaption></figure>

proportion goes up to 90% of the impacts in the annulus accretion simulations. A similar analysis for BSE material leads to a proportion of vapor-forming impacts of 68% and 93% for classical and annulus cases, respectively. These estimations are valid for head-on planar collisions. More realistic impact geometries would reduce these numbers. Taking into account a maximum frequency of impact angles of 45° (E. Pierazzo & H. J. Melosh 2000), the threshold velocities are reduced to $U_{mc} = 10 \text{ km s}^{-1}$ and $U_{mp} = 8.6 \text{ km s}^{-1}$ for chondritic and pyro-litic bodies, respectively. This, in turn, reduces the number of impacts that can produce vapor to 43% and 51% for pyrolite and chondrule—respectively—in the classical model.

Furthermore, retracing the sequence of impacts that lead to the formation of Earth analogs (simulated planets with orbits and masses similar to Earth; see Section 2.4) allows us to

evaluate the significance of vapor formation during these events. We find that 61% and 89% of the chondritic impacts in the classical and annulus models, respectively, produced vapor (see Figure 4). These numbers are comparable to those obtained for Venus analog planets (65% and 92%, respectively). The proportion of vapor-forming impacts during the history of Mars analogs in the classical model is also close to that of Earth and Venus (56%). However, Mars experienced less vaporizing impacts in the annulus model than Earth and Venus (77%). Thus, an analysis of volatile and moderately volatile elements could give credit to the annulus accretion model if it shows that Mars is enriched in volatiles compared to Earth.

A complementary analysis showed that, within the results of *N*-body simulations, planetesimal impacts occurred in classical
and annulus models at higher velocities than embryo impacts. They then present more vaporizing impacts (see Figure A5). This is due to the higher eccentricity of planetesimals relative to embryos (E. Kokubo & S. Ida 2000). Moreover, impact velocities do not increase dramatically with time (see Figure A6). It is then unlikely that the dichotomy we see between impact velocities of embryo-embryo and embryo- planetesimal collisions is due to a timing effect. This higher vaporization rate for embryo-planetesimal impacts is consistent with E. D. Young et al. (2019) and R. C. Hin et al. (2017), suggesting that the enrichment of heavy Si and Mg isotopes of differentiated rocky bodies is due to vaporizing events during planetesimal growth.

## 5. Conclusion

Our findings suggest that vapor production during impacts between planetesimals could be a more common phenomenon than previously thought (Figure A4). Given the substantial number of vapor-producing impacts, they should be further considered in N-body simulations of impacts. Whether or not vapor is lost postimpact depends significantly on factors such as the mass, velocity, and impact angle of the colliding bodies. Impacts between small planetesimals might produce vapor, but probably only in small amounts, though this could be easily lost by ejection. As discussed above, large impacts would produce more vapor, but this might be trapped in the gravitational potential well that forms after the impact. In this case, the vapor would not escape but fall back into the condensing central body. SPH simulations that specifically use the entropy of the materials involved in the impact, like the ones we computed here, offer advanced predictive capabilities by accurately estimating the volume and fate of vapor produced during impacts. Our findings also suggest that impacts affecting the volatile budget of a planet could occur not only after giant impacts such as the Moon-forming impacts but also may be frequent in smaller-scale impacts. These events, often under-estimated, could significantly influence the mass-balance evolution and volatile distribution among accreting planetesimals. This significant proportion of vaporizing impacts we predict could affect our view of the formation of the terrestrial planets. E. A. Pringle & F. Moynier (2017) proposed that the rubidium isotopic signature of the Moon could be due to the loss of rubidium during vaporization events. They suggested that this evaporation could have occurred during the Moon-forming impact or earlier during the accretion of Theia. Potentially different ratios of volatiles and moderately volatile elements in the Martian moons could be related to an impact origin of Phobos. These insights into the dynamics of vaporization during planetary accretion offer pivotal corrections to existing planetary formation and evolution models, potentially altering our understanding of early solar system dynamics.

## Acknowledgments

A.S. acknowledges Charles Lelosq's help and constructive comments. R.C. acknowledges financial support from the UPC Labex UnivEarthS (project VADIS), as well as support from the Research Council of Norway, project numbers 332523 (PHAB), and 325567 (HIDDEN), and access to supercomputing facilities via eDARI stl2816, PRACE RA4947 and RA240046, and Uninet2 NN9697K grant. S.N.R. acknowledges funding from

the French Programme National de Planétologie (PNP) and in the framework of the Investments for the Future programme IdEx, Université de Bordeaux/RRI ORIGINS.

# Appendix
## Thermodynamic Integration and the Calculation of the Entropy

MD simulations cannot directly estimate thermal quantities, such as free energy and entropy, of a nonideal system (A. R. Leach 2001; D. Frenkel & B. Smit 2002). However, these quantities can be derived starting from the values of an ideal system, like an ideal gas.

We start with the free energy $F_0$ of an ideal gas system expressed as a function of its partition function as

$$F_0 = -k_B T \times \ln \frac{(N \times V)^N}{\Lambda^{3N} N!} = -k_B T \left[ N \ln \frac{V}{\Lambda^3} - \ln N! \right],$$

where $\Lambda$ is the Broglie wavelength, defined as

$$\Lambda = \frac{h}{\sqrt{2\pi mk_B T}}, \qquad (A1)$$

where $k_B$ is the Boltzmann constant, $h$ the Planck constant, $T$ the temperature, $V$ the volume, and $N$ the number of particles.

Then we express the free energy $F_1$ of the ab initio system starting from the ideal gas and adding a coupling term, $\Delta F_{0 \to 1}$, as in F. Dorner et al. (2018):

$$F_1 = F_0 + \Delta F_{0 \to 1}. \qquad (A2)$$

In practice, the coupling between the ab initio system and the ideal gas is realized using a linear mixing between the two descriptions of the dynamical states of the same system of particles. In this formalism, the forces, $f$, between the atoms can be evaluated as a linear combination stretching from the pure ab initio system to the pure ideal gas system:

$$f = \lambda f_1 + (1 - \lambda) f_0. \qquad (A3)$$

$\lambda$ takes into account the part of the dynamics that comes from the ab initio system and $(1 - \lambda)$ from the ideal gas. The coupling between the ab initio system and the ideal gas represents the integral along a thermodynamic path that couples these two systems:

$$\Delta F_{0 \to 1} = \int_{\lambda=0}^{\lambda=1} \frac{\partial F}{\partial \lambda}, \qquad (A4)$$

which is equivalent to (A. R. Leach 2001; D. Frenkel & B. Smit 2002)

$$\Delta F_{0 \to 1} = \int_{\lambda=0}^{\lambda=1} d\lambda \left\langle \frac{\partial U(\lambda)}{\partial \lambda} \right\rangle_{\lambda}, \qquad (A5)$$

and to (F. Dorner et al. 2018)

$$\Delta F_{0 \to 1} = \int_{\lambda=0}^{\lambda=1} d\lambda \langle V(\lambda) \rangle_{\lambda}, \qquad (A6)$$

where $U$ is the internal energy, $V$ is the potential energy, and the bracket denotes an average over a (well-thermalized and long-enough) simulation at constant $\lambda$.
<figure><img src="image_7.png" /><figcaption>**Figure A1.** Plot of $P(u_p)$ and $P(u_{p,imp})$ for a CMm target and a CMc impactor traveling at 10 km s$^{-1}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure A2.** Remaining bodies at the end of the classical simulations (a) and annulus simulations (b). The yellow, blue, and red rectangles indicate, respectively, Venus, Earth, and Mars analog planets.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure A3.** Determination of the peak impact pressure in the target in the planar shock collision approximation for chondrule (a) and pyrolite (b) material. The Hugoniot equations of the state of the two materials are plotted with opposite signs and cross at a given impact velocity $U_{mc}$ (respectively, $U_{mp}$). The corresponding peak pressure, $P_{mc}$ (respectively, $P_{mp}$), is read on the vertical axis (black horizontal arrow). Values are for impacts between two chondrule bodies and two pyrolitic bodies, respectively.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure A4.** Proportions of impacts producing vaporization that occurred in *N*-body simulations with a classical starting configuration (a) and annulus starting configuration (b). The color represents the type of material that is vaporized.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure A5.** Proportions of impacts producing vaporization that occurred in *N*-body simulations in between two embryos ((a) and (b)) and in between an embryo and a planetesimal ((c) and (d)) with a classical starting configuration ((a) and (c)) and annulus starting configuration ((b) and (d)). Colors represent the type of material that is vaporized.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure A6.** Relative velocities of impacts as a function of time for Venus ((a) and (b)), Earth ((c) and (d)), and Mars analogs ((e) and (f)) in *N*-body simulations with a classical starting configuration ((a), (c), and (e)) and annulus starting configuration ((b), (d), and (f)).</figcaption></figure>

<table><caption>**Table A1**

Chemical Composition of the Three Materials Described in This Study</caption><thead><tr><th rowspan="2"></th><th colspan="6">Oxide (wt%)</th><th rowspan="2">Molar Ratio<br>(Mg+Fe)/(Si+Al)</th><th colspan="3">Density (g cm<sup>-1</sup>)</th></tr><tr><th>FeO</th><th>SiO<sub>2</sub></th><th>MgO</th><th>Al<sub>2</sub>O<sub>3</sub></th><th>CaO</th><th>Na<sub>2</sub>O</th><th>ρ<sub>0l</sub></th><th>ρ<sub>0a</sub></th><th>ρ<sub>0h</sub></th></tr></thead><tbody><tr><td>Chondrule (CMc)</td><td>6.98</td><td>46.57</td><td>41.13</td><td>1.98</td><td>2.17</td><td>1.16</td><td>1.38</td><td>3.35</td><td>3.5</td><td>3.65</td></tr><tr><td>Matrix (CMm)</td><td>7.12</td><td>51.45</td><td>29.02</td><td>6.72</td><td>3.69</td><td>1.98</td><td>0.829</td><td>3.0</td><td>3.25</td><td>3.5</td></tr><tr><td>Pyrolite (BSE)</td><td>8.94</td><td>44.68</td><td>37.23</td><td>4.75</td><td>3.48</td><td>0.93</td><td>1.25</td><td>...</td><td>3.1</td><td>...</td></tr></tbody></table>

**Note.** CMc = chondrule composition, CMm = matrix composition, and BSE = pyrolite/bulk silicate Earth composition.
$\langle V(\lambda) \rangle_{\lambda}$ is estimated from a series of MD simulations. Once the free energy of the system is evaluated, the entropy becomes

$$S = \frac{U - F}{T}. \qquad (A7)$$

The Appendix presents one additional table and several additional figures that support the main findings of this study. Table A1 presents the composition of the materials used in this study and the densities of those materials before impacts. Figure A1 illustrates the determination of the maximum pressure for a CMm target and a CMc impactor. It highlights the range of impact pressures encountered in our simulations. Figure A2 shows the final distribution of planetary bodies in both the classical and the annulus N-body simulation setups, with Venus, Earth, and Mars analogs highlighted. Figure A3 details the determination of the peak impact pressures in the planar shock collision approximation for bodies with chondrule and pyrolite compositions, based on their respective Hugoniot equations of state. Figures A4 and A5 present the fraction of impacts that produce vaporization in different simulation configurations, distinguishing between classical and annulus scenarios, as well as between embryo-embryo and embryo-planetesimal collisions. Finally, Figure A6 tracks the evolution of the impact velocities over time for Venus, Earth, and Mars analogs, providing insights into how impact dynamics evolve during the late stages of planetary accretion.

## ORCID iDs

Adrien Saurety https://orcid.org/0000-0001-6967-9718

Razvan Caracas https://orcid.org/0000-0003-4586-776X

Sean N. Raymond https://orcid.org/0000-0001-8974-0758

## References

Agnor, C. B., Canup, R. M., & Levison, H. F. 1999, Icar, 142, 219

Ahrens, T. J., & O'Keefe, J. D. 1972, Moon, 4, 214

Ahrens, T. J., Petersen, C. F., & Rosenberg, J. T. 1969, JGR, 74, 2727

Asimow, P. D., & Ahrens, T. J. 2010, JGRB, 115, B10209

Blöchl, P. E. 1994, PhRvB, 50, 17953

Boss, A. P. 2004, ApJ, 610, 456

Canup, R. M., & Asphaug, E. 2001, Natur, 412, 708

Car, R., & Parrinello, M. 1985, PhRvL, 55, 2471

Caracas, R. 2024, E&PSL, 637, 118724

Caracas, R., Kobsch, A., Solomatova, N. V., et al. 2021, J. Vis. Exp., 175, e61534

Caracas, R., & Stewart, S. T. 2023, E&PSL, 608, 118014

Chambers, J. E. 2001, Icar, 152, 205

Chambers, J. E. 2013, Icar, 224, 43

Clement, M. S., Kaib, N. A., Raymond, S. N., Chambers, J. E., & Walsh, K. J. 2019, Icar, 321, 778

Clement, M. S., Kaib, N. A., Raymond, S. N., & Walsh, K. J. 2018, Icar, 311, 340

Čuk, M., & Stewart, S. T. 2012, Sci, 338, 1047

Davies, E. J., Carter, P. J., Root, S., et al. 2020, JGRE, 125, e2019JE006227

Dorner, F., Sukurma, Z., Dellago, C., & Kresse, G. 2018, PhRvL, 121, 195701

Emsenhuber, A., Cambioni, S., Asphaug, E., et al. 2020, ApJ, 891, 6

Fischer, R. A., & Ciesla, F. J. 2014, E&PSL, 392, 28

Forbes, J. W. 2012, Shock Wave Compression of Condensed Matter: A Primer, Shock Wave and High Pressure Phenomena (Berlin: Springer)

Frenkel, D., & Smit, B. 2002, Understanding Molecular Simulation: from Algorithms to Applications (2nd ed.; San Diego: Academic)

Friend, P., Hezel, D. C., Barrat, J.-A., et al. 2018, M&PS, 53, 2470

Gaches, B. A. L., Walch, S., Offner, S. R., & Münker, C. 2020, ApJ, 898, 79

Genda, H., Kokubo, E., & Ida, S. 2012, ApJ, 744, 137

Hansen, B. M. S. 2009, ApJ, 703, 1131

Heighway, P. G., Sliwa, M., McGonegley, D., et al. 2019, PhRvL, 123, 245501

Hewins, R., Connolly, H., & Lofgren, G. 2005, in ASP Conf. Ser. 341, Chondrites and the Protoplanetary Disk, ed. A. N. Krot, E. R. D. Scott, & B. Reipurth (San Francisco, CA: ASP), 286

Hewins, R. H. 1997, AREPS, 25, 61

Hin, R. C., Coath, C. D., Carter, P. J., et al. 2017, Natur, 549, 511

Hohenberg, P., & Kohn, W. 1964, PhRv, 136, B864

Hoover, W. G. 1985, PhRvA, 31, 1695

Isbell, W. M. 2005, Shock Waves: Measuring the Dynamic Response of Materials (London: Imperial College Press)

Izidoro, A., Raymond, S. N., Morbidelli, A., & Winter, O. C. 2015, MNRAS, 453, 3619

Johansen, A., & Okuzumi, S. 2018, A&A, 609, A31

Johansen, A., Ronnet, T., Bizzarro, M., et al. 2021, SciA, 7, eabc0444

Kobsch, A., & Caracas, R. 2020, JGRE, 125, e2020JE006412

Kohn, W., & Sham, L. J. 1965, PhRv, 140, A1133

Kokubo, E., & Ida, S. 2000, Icar, 143, 15

Kraus, R. G., Root, S., Lemke, R. W., et al. 2015, NatGe, 8, 269

Kraus, R. G., Stewart, S. T., Swift, D. C., et al. 2012, JGRE, 117, E09009

Kresse, G., & Furthmüller, J. 1996a, ComMS, 6, 15

Kresse, G., & Furthmüller, J. 1996b, PhRvB, 54, 11169

Kresse, G., & Hafner, J. 1993, PhRvB, 47, 558

Kurosawa, K., Kadono, T., Hironaka, Y., et al. 2013, LPSC, 44, 2537

Leach, A. R. 2001, Molecular Modeling: Principles and Applications (2nd ed.; New York: Prentice Hall)

Leinhardt, Z. M., & Stewart, S. T. 2012, ApJ, 745, 79

Li, Z., Caracas, R., & Soubiran, F. 2020, E&PSL, 547, 116463

McCoy, C. A., Marshall, M. C., Polsin, D. N., et al. 2019, PhRvB, 100, 014106

McDonough, W., & Sun, S.-s. 1995, ChGeo, 120, 223

McQueen, R. G., Marsh, S. P., & Fritz, J. N. 1967, JGR, 72, 4999

Melosh, H. J. 2007, M&PS, 42, 2079

Mermin, N. D. 1965, PhRv, 137, A1441

Morbidelli, A., Lunine, J., O'Brien, D., Raymond, S., & Walsh, K. 2012, AREPS, 40, 251

Mysen, B. O., Virgo, D., & Seifert, F. A. 1982, RvGeo, 20, 353

Norris, C. A., & Wood, B. J. 2017, Natur, 549, 507

Nosé, S. 1984, JChPh, 81, 511

O'Neill, H. S., & Palme, H. 2008, RSPTA, 366, 4205

Ostrowski, D., & Bryson, K. 2020, M&PS, 55, 2007

Perdew, J. P., Ernzerhof, M., & Burke, K. 1996, JChPh, 105, 9982

Pierazzo, E., & Melosh, H. J. 2000, AREPS, 28, 141

Pringle, E. A., & Moynier, F. 2017, E&PSL, 473, 62

Quintana, E. V., & Barclay, T. 2016, American Astronomical Society, AAS Meeting, 228, 404.05

Ragan, C. E. 1982, PhRvA, 25, 3360

Raymond, S. N., & Izidoro, A. 2017, SciA, 3, e1701138

Raymond, S. N., & Morbidelli, A. 2022, in Demographics of Exoplanetary Systems: Lecture Notes of the 3rd Advanced School on Exoplanetary Science, ed. K. Biazzo et al. (Cham: Springer), 3

Raymond, S. N., O'Brien, D. P., Morbidelli, A., & Kaib, N. A. 2009, Icar, 203, 644

Root, S., Townsend, J. P., Davies, E., et al. 2018, GeoRL, 45, 3865

Shang, J. L., Shen, L. T., & Zhao, J. 2000, Int. J. Rock Mech. Min. Sci., 37, 705

Stewart, S., Davies, E., Duncan, M., et al. 2020, in AIP Conf. Proc. 2272 (Melville, NY: AIP), 080003

Thomas, C. W., & Asimow, P. D. 2013, JGRB, 118, 5738

Villeneuve, J., Libourel, G., & Soulié, C. 2015, GeCoA, 160, 277

Young, E. D., Shahar, A., Nimmo, F., et al. 2019, Icar, 323, 1

Zel'dovich, B. Y., & Raizer, Y. P. 2002, Physics of Shock Waves and High-temperature Hydrodynamic Phenomena, Reprod edn. (New York: Dover)