Here we examine the structure of terrestrial planets over a wide range of thermal and rotational states. We used smoothed particle hydrodynamics (SPH) simulations and developed a new code, the Highly Eccentric Rotating Concentric U (potential) Layers Equilibrium Structure (HERCULES) code, for calculating the equilibrium structure of bodies based on modeling the body as a series of overlapping spheroids (section 2 ). We find that there is a wide range of sizes and shapes for hot, rotating bodies (section 3 ). In particular, we show that there is a limit beyond which a body cannot have a single angular velocity. Beyond this limit, a body can exhibit a range of morphologies with disk‐like outer regions. The corotation limit (CoRoL) is a function that depends upon the composition, thermal state, AM, and mass of a body (section 4 ). We show that typical terrestrial planets experience multiple substantially vaporized states during accretion as a result of giant impacts. For the expected range of AM of growing planets, some post‐impact states exceed the CoRoL (section 5 ). The range of possible physical and thermal structures for terrestrial planets has significant implications for understanding physical processes during accretion, for example, differentiation, satellite formation, and interpreting the composition of exoplanets (section 6 ). The supporting information includes an extended description of the methods used in this work and data tables.

The need to understand the physical structures of rocky planets extends beyond accretionary processes. Internal structure models are the primary tool used to infer the possible compositions of exoplanets from mass and radius observations [e.g., Valencia et al. , 2006 ; Seager et al. , 2007 ; Swift et al. , 2012 ]. Most of the discovered exoplanets fall in a mass and radius range between the rocky and gaseous planets in our solar system [ Morton et al. , 2016 ]. Hence, most exoplanets are unlike any of the planets in our solar system. In addition, because of the biases in the astronomical techniques used to find them, most of the known exoplanets are close to their stars and have high equilibrium surface temperatures. However, current internal structure models do not consider a wide range of thermal or rotational states, especially for rocky planets. Most exoplanet structure calculations model non‐rotating, relatively cold bodies similar to the planets in the present‐day solar system [e.g., Zharkov and Trubitsyn , 1978 ; Valencia et al. , 2006 ; Swift et al. , 2012 ; Hubbard , 2013 ; Zeng and Sasselov , 2013 ; Zeng et al. , 2016 ; Unterborn et al. , 2016 ]. In order to understand the possible range of exoplanet compositions, rocky planet models must be extended to include hot and rotating structures.

Exploring the range of physical structures after giant impacts is particularly important for understanding the mechanisms for terrestrial satellite formation, including the origin of our Moon and its chemical relationship with Earth [e.g., Canup et al. , 2015 ; Lock et al. , 2016 ; Wang and Jacobsen , 2016 ]. Previous studies of lunar origin have focused on the conditions after the canonical Moon‐forming impact, which we summarize here. After a canonical impact event, the Earth's spin period would have been ∼5 h, which would have produced minor rotational flattening for a fully condensed body. There would have been a large angular velocity discontinuity between the corotating planet and a sub‐Keplerian disk. The dynamics of this shear boundary have not been well studied, except to consider possible mixing between the planet and disk [e.g., Pahlevan and Stevenson , 2007 ; Desch and Taylor , 2013 ; Melosh , 2014 ]. The material that is injected into orbit in canonical style impacts originates primarily from the antipode hemisphere of the impactor and is less shocked than the impacted hemisphere [ Canup and Asphaug , 2001 ; Canup , 2004 ; Nakajima and Stevenson , 2014 ]. As a result, the disk material would have been less vaporized than the planet's predominantly silicate vapor atmosphere, creating a thermodynamic discontinuity between the Earth and the disk [although the disk vapor may have been continuous with the planet's atmosphere; see Pahlevan and Stevenson , 2007 ; Desch and Taylor , 2013 ]. Thus, studies of the system after Earth's last giant impact assume that a dynamically and thermodynamically distinct disk orbits a corotating, nearly spherical central body [e.g., Thompson and Stevenson , 1988 ; Canup and Asphaug , 2001 ; Canup , 2004 ; Machida and Abe , 2004 ; Ćuk and Stewart , 2012 ; Canup , 2012 ; Nakajima and Stevenson , 2014 ; Charnoz and Michaut , 2015 ]. Given the possible range of thermal and rotational states produced by giant impacts, it is necessary to explore satellite formation from a range of different post‐impact structures.

The physical structures of bodies after giant impacts have not been well studied with the exception of the Earth after the proposed Moon‐forming event [ Hartmann and Davis , 1975 ; Cameron and Ward , 1976 ]. Most studies have assumed that this event set the present AM of the Earth‐Moon system, as proposed by Cameron and Ward [ 1976 ]. Detailed studies of giant impact outcomes find that this constraint strongly limits the mass ratio and impact parameters of a potential Moon‐forming event [ Canup , 2004 , 2008a ]. Hence, the canonical model for the Moon‐forming giant impact is that of an approximately Mars‐mass body obliquely colliding with the proto‐Earth near the mutual escape velocity [ Canup and Asphaug , 2001 ; Canup , 2004 , 2008a ]. This specific giant impact is often used to represent all giant impacts. However, impacts similar to the canonical Moon‐forming impact represent only a small fraction of the impacts that occur during accretion [ Stewart and Leinhardt , 2012 ]. Studies of terrestrial planet formation with more realistic collision outcomes have found that giant impacts with much higher specific energies [e.g., Quintana et al. , 2016 ] and AM [ Kokubo and Genda , 2010 ] than the canonical impact are common. The resulting post‐impact structures have a correspondingly wide range of thermal and rotational states that can be substantially different from that produced by the canonical Moon‐forming impact. Furthermore, the canonical impact model for Moon formation is currently under scrutiny because of the difficulty of explaining the isotopic similarity between Earth and the Moon [e.g., Melosh , 2014 ; Burkhardt , 2014 ] and the discovery that AM could have been transferred away from the Earth‐Moon system since the last giant impact [ Ćuk and Stewart , 2012 ; Wisdom and Tian , 2015 ; Ćuk et al. , 2016 ; Tian et al. , 2017 ]. If the AM of the Earth‐Moon system immediately after the event was different than the present day, then a wider range of giant impact scenarios are possible [ Ćuk and Stewart , 2012 ; Canup , 2012 ; Lock et al. , 2016 ]. Thus, a quantitative study of the possible physical structures of hot, rotating rocky bodies is needed.

Most planet formation studies have implicitly assumed that the lifetime of the transiently hot states generated by giant impacts [hundreds to thousands of years, Zahnle et al. , 2007 ] is short enough to neglect a detailed treatment of partially vaporized rocky bodies. For some aspects of terrestrial planet formation, this assumption may be justified, and most investigations of post‐impact processes begin after the silicate fraction of a body is fully condensed and the planet has a distinct magma ocean and volatile‐dominated atmosphere [e.g., Solomatov , 2000 ; Lebrun et al. , 2013 ; Hirschmann , 2012 ]. However, the physical structure of a body in the time period between an energetic impact and a magma ocean can influence a number of different processes, including core formation, the reaccretion of impact debris, and the formation of satellites. Furthermore, recent studies have suggested that the Earth's lower mantle records chemical signatures that predate the last giant impact [ Rizo et al. , 2016 ; Mukhopadhyay , 2012 ; Tucker et al. , 2012 ; Parai et al. , 2012 ; Pető et al. , 2013 ]. Ascertaining how such reservoirs survived at least the Moon‐forming impact event and persisted to the present day requires an understanding of the range of pressure‐temperature conditions generated by giant impacts.

The physical structure of a planet is essential information that is needed to investigate all major physical processes during planet formation and evolution. The planet's mass, size, and shape govern the efficiency of accretion onto the body, and the internal pressure and temperature profiles control the conditions for differentiation and core formation. The vigor of convection and the resulting extent of mixing depend upon the physical structure of the body as well as the driving forces. The energy deposited by giant impacts, a key stage of terrestrial planet formation in our solar system and exoplanetary systems [ Chambers , 2010 ], can radically change a body's physical structure. Significant portions of the silicate mantles of the impacting bodies melt or vaporize [ Melosh , 1990 ]. In addition, planets can acquire significant angular momentum (AM) via one or more giant impacts [ Agnor et al. , 1999 ; Kokubo and Ida , 2007 ; Kokubo and Genda , 2010 ].

In general, HERCULES is a more efficient and accurate method to model planetary structure compared to SPH. However, the current version of HERCULES can only calculate the structure of corotating bodies. A version of HERCULES that may solve for an imposed angular momentum structure is under development. As a result, we utilize SPH to study the structure of non‐corotating bodies.

The limited resolution of the SPH simulations leads to an error in resolving the outer boundary of bodies, particularly for silicate specific entropies above the critical point. The SPH resolution is not sufficient to capture the mass of material in the low‐density, partially vaporized regions of hot bodies. For the example planet with an isentropic profile below the critical point entropy (Figures 2 a– 2 c), the SPH calculation does not resolve pressures below about 2 × 10 4 bar. The scale height at the edge of a mostly condensed body is small, so the equatorial radius is only slightly smaller than that calculated with the HERCULES code with bar, about 6% in our example case. However, for substantially vaporized bodies with specific entropies above the critical point value (e.g., Figures 3 d– 3 f with S lower =6.5 kJ K −1 kg −1 ), the scale height of the low‐pressure regions is larger. Although the minimum pressure of the SPH planet is lower, at 1257 bar in our example case, the difference in equatorial radius is greater, about 15% smaller than that calculated with the HERCULES code with bar. To mimic the effect of the resolution limit in the SPH simulations, we also calculated a HERCULES body using a bounding pressure of bar (blue lines in Figures 3 d– 3 f), the lowest pressure in the midplane of the SPH planet. In this case, the equatorial radius agreed with the SPH to within about 0.2%. The mass of the HERCULES body that lies outside of the surface of the SPH structure (outer red line in Figure 3 d) and within the HERCULES 10 bar contour (black dashed line) is only about 10 −3 M Earth . The mass of this unresolved outer, partially vaporized layer is well below the resolution of the SPH simulation, where each particle has a mass of about 10 −5 M Earth . For all the isolated, corotating planets considered here, the differences between HERCULES and SPH calculations typically involved less than the mass of 500 SPH particles (<0.5 wt %), resulting in errors in the SPH equatorial radii less than about 15%. Since there is little mass in the region not resolved by the SPH planets, the corotating angular velocity of the bodies calculated using SPH and HERCULES is not substantially different (Table 1 ). SPH does not define the polar radius as well as the equatorial radius due to the difference in particle resolution along each axis of an oblate planet. This leads to a small error in the aspect ratio, the ratio of the polar to equatorial radii (Table 1 ).

A comparison of equilibrium planetary structures calculated using the GADGET‐2 SPH code and the HERCULES code for a largely condensed body ((a–c) silicate specific entropy=5.0 kJ Kkg) and a substantially vaporized body ((d–f) silicate specific entropy=6.5 kJ Kkg). All bodies are Earth mass, have an angular momentum equal to that of the present‐day Earth‐Moon system () and an isentropic thermal profile (class I). SPH bodies are shown in red and HERCULES bodies with a bounding pressure,, of 10 bar are in black. In Figures 3 d– 3 f, a HERCULES body with a bounding pressure of 1257 bar, the lowest pressure in the midplane of the corresponding SPH body, is also shown in blue. Figures 3 a and 3 d show the axisymmetric SPH and HERCULES pressure contours on a quarter plane through the rotation axis. The pressure contours correspond to specific layers in the HERCULES code and are not the same pressures in Figures 3 a and 3 d. The outer contour is always the bounding surface of the body. In Figure 3 d, only the surface of the outermost layer is shown for the HERCULES body withbar. The lower panels display midplane pressure (Figures 3 b and 3 e) and density (Figures 3 c and 3 f) radial profiles.

The GADGET‐2 and HERCULES codes implement fundamentally different approaches for modeling the equilibrium structure of isolated bodies. Nevertheless, the two techniques produce very similar structures for corotating bodies. Two examples of isentropic Earth‐mass planets with a total angular momentum equal to the present‐day Earth‐Moon system ( L EM =3.5 × 10 34 kg m −2 s −1 ) are shown in Figure 3 . The SPH calculation is shown in red and the HERCULES calculation, using a bounding pressure of bar, is plotted in black. The well‐known issues with resolving boundaries with high‐density contrasts in SPH can be seen at the core‐mantle boundary, with a layer of anomalous density particles on either side of the boundary. However, the error at the boundary does not propagate to the rest of the structure. Similarly, cooler SPH planets that have a sharp boundary with vacuum have an outer layer of particles with an increased radial separation compared to interior particles (Figures 3 a– 3 c with S lower =5 kJ K −1 kg −1 ), but this layer captures the correct density and pressure at that radius.

We expect that the accuracy of the EOS model for silicates will be improved as new data are acquired on forsterite and other silicate chemical compositions [e.g., Bolis et al. , 2016 ; Sekine et al. , 2016 ; Root et al. , 2016 ; Davies et al. , 2016 , 2017 ]. In previous studies, the M‐ANEOS model has been shown to underestimate the gain in entropy at high shock pressures [e.g., Kraus et al. , 2012 ]. The model shock temperatures for forsterite overlap with the error bars reported by Bolis et al. [ 2016 ] up to about 500 GPa, which corresponds to a specific entropy slightly greater than the model critical point when shocked from standard pressure and temperature. As a result, the EOS model is inferred to underestimate the entropy gain at shock pressures >∼500 GPa. Thus, the impact outcomes from the highest velocity events are likely to shift to a greater degree of vaporization with improving EOS models, but our general conclusions about changes in planetary structure relative to the degree of vaporization are robust.

We have used specific entropy, rather than temperature, as the natural intensive variable to describe the thermal state for each layer within a planet. For liquid‐vapor mixtures, the temperature is insufficient to determine the relative proportions of each phase at a given pressure. With specific entropy as the independent variable, the lever rule can be applied to determine the mass fraction of each phase using the material's liquid‐vapor phase boundary, which is a dome in specific entropy‐pressure space (Figure 2 ). In rocky planets, the thermal state of the mantle is often represented by the potential temperature, which is the temperature on the mantle adiabatic at a reference pressure (usually 1 bar for terrestrial applications). The potential temperature is degenerate for thermal structures that intersect the liquid‐vapor phase boundary (Figure S7 ), so it is not a useful parameter in this work. In addition, partially vaporized planets do not have a well‐defined surface, and we must report the radius of the body at a specific pressure contour.

A stratified profile (class III) divides the silicate portion of a body into two separate material layers containing a specified mass fraction of the silicate. The lower layer is isentropic with specific entropy S lower . The upper layer has the same thermal structure as a vapor atmosphere profile: isentropic at higher pressures with specific entropy S upper and following the saturated vapor curve at lower pressures (blue line in Figure 2 ). The pressure at the boundary between the upper and lower layers varies between different bodies since the masses of the two silicate layers are dictated, not the pressure of transition. Here we considered stratified profiles with an upper layer that contains 50 or 25% of the mass of the silicate. These stratified profiles are intended to emulate the thermal structures that are typically produced in hydrodynamic simulations of giant impacts. The impacting hemispheres of the colliding bodies are more highly shocked than the antipodal hemispheres. After gravitational equilibration, the post‐impact body is thermally stratified with lower entropy material from the antipodal hemispheres at higher pressures and higher entropy material from the impacted hemispheres at lower pressures. Post‐impact thermal profiles are much more variable than a simple two‐layer structure, but the model stratified thermal profile allows us to examine the general effect of thermal stratification on planetary structure.

Vapor atmosphere profiles (class II) are isentropic at high pressure with specific entropy, S lower , but if a profile intersects the liquid‐vapor phase boundary, it follows the vapor side of the mixed phase region (red line in Figure 2 ). In this case, regions with pressures below the intersection point are assumed to be pure vapor on the phase boundary. The outer layers of such structures represent a saturated silicate vapor atmosphere and have a greater specific entropy compared to an isentropic thermal profile defined by the same S lower (green line). The vapor atmosphere profile approximates the scenario where condensate efficiently rains out from the low‐pressure regions of a body and the vapor is stratified. In a hot rocky body with a turbulent atmosphere, the rain out of condensates is unlikely to be wholly efficient, but the vapor atmosphere profile provides a useful end member case.

Examples of the different silicate thermal profiles considered for the calculation of the structure of isolated bodies. The silicate portion of each body was modeled as a simple system with two condensed phases and vapor using an equation of state for forsterite. The colored dashed lines correspond to the different classes of thermal profiles described in section 2.3 : (I) isentropic (green), (II) vapor atmosphere (red), and (III) stratified (blue). The solid black curve is the model liquid‐vapor phase boundary for forsterite with the critical point (=5.40 kJ Kkg=25.5 kbar,=8810 K,=1680 kg m) shown by the black dot. Material in a thermal state plotting within the liquid‐vapor phase boundary is a mixture of liquid and vapor with thermal states on each side of the boundary. Isentropes in pressure‐temperature space are shown in Figure S7

For the silicate portion of bodies, we used three different classes of thermal profiles: (I) isentropic, (II) vapor atmosphere, and (III) stratified (Figure 2 ). For isentropic profiles (class I, green line in Figure 2 ), the silicate portion of the body has a constant specific entropy, S lower , which is the simplest possible thermal state with which to make comparisons. Except for the coldest bodies considered in this study, the outer layers of these structures intersect the liquid‐vapor phase boundary. In these cases, the low‐pressure portions of the isentropic profile describe an ideal mixture of liquid and vapor, without any phase separation.

We considered isolated bodies with a range of thermal profiles. The variations in the thermal state of the core during accretion is not well understood. Here we largely neglect these variations to focus on the thermal state of the silicate component. Each isolated body had an iron core of fixed specific entropy S core =1.5 kJ K −1 kg −1 . This core isentrope has a temperature of ∼3800 K at the pressure of the present‐day core‐mantle boundary, similar to the present thermal state of Earth's core [e.g., Anzellini et al. , 2013 ]. The effect of varying core entropy is discussed in section 4 .

As in the SPH calculations, we modeled Earth‐like bodies with 2/3 silicate and 1/3 iron by mass and used the same equations of state derived from the M‐ANEOS model. A range of different thermal profiles were considered for the silicate (section 2.3 ). Unless noted otherwise, 100 constant‐density layers were used: 20 for the iron core and 80 for the silicate. When the silicate was divided into two layers with different thermal states, 40 layers were used for each layer. The equipotential surfaces were defined using N μ =1000 points, and terms up to spherical harmonic degree 12 were included in the iterative equation. The HERCULES code requires a bounding pressure for the planet at the surface of the outermost layer. Unless otherwise stated, we assumed a bounding pressure of bar.

The HERCULES code cannot reach the same level of precision as the CMS model. The varying integration domains in the full solution to the Poisson's equation prevent the use of Gaussian quadrature, resulting in less precise numerical integration. However, the HERCULES code is capable of efficiently calculating a variety of planetary structures, including rapidly rotating planets, making it a versatile tool to use in a wide range of problems in planetary science.

Schematic of how an axisymmetric planetary structure is modeled in the HERCULES code. A body is described as a superposition of a number of constant‐density spheroids of density δ ρ i (shown in gray). The superposition of the spheroids gives a body with increasing density with depth. The volumes between successive spheroids are called layers. Each layer has a constant density, ρ i , given by the sum of the densities of all the spheroids larger than the spheroid that defines the inner edge of the layer. Each layer belongs to a material layer which determines the relationship between pressure and total density in that layer by use of the material's equation of state. Two material layers are shown separately in blue and orange.

In order to be able to model rapidly rotating bodies, we built on the work of Kong et al. [ 2013 ] and Hubbard [ 2012 , 2013 ] to develop a new code. Here we summarize the method, and full details are presented in Text S1 . The HERCULES code models a body as a superposition of a number of overlapping, concentric, constant‐density spheroids (Figure 1 ). The surface of each spheroid is defined by an equipotential surface. We refer to the region between two equipotential surfaces as a layer, and each layer is defined as being composed of a single material (e.g., silicate or iron). The density distribution in the planet is given by the superposition of the density of each of the spheroids. The total density of a layer in the body, ρ i , is given by the sum of the density of all the spheroids, δ ρ i , that have an equatorial radius larger than the inner edge of the layer. The code iteratively solves for the equilibrium structure and shape of the body while conserving the mass of each of the materials and the total AM. During each iteration, the density of the spheroids is altered so that the total density is consistent with a hydrostatic pressure profile, the EOS of the materials, and an imposed thermal state for each of the material layers. Each material's mass is conserved by altering the radius of each of the concentric layers. Although we refer to the volumes as spheroids, in our formulation the equipotential surfaces can be any rotationally symmetric surface whose radius decreases monotonically from equator to pole.

We developed a new code for calculating the equilibrium structure of rotating planetary bodies based on describing a body as a series of overlapping, concentric, constant‐density spheroids. A model based on the same principle, the concentric Maclaurin Spheroid (CMS) model, was originally developed by Hubbard [ 2012 , 2013 ] for planets with slow rotation rates. The CMS model is a self‐consistent field method that iteratively solves for the equilibrium shape of the surfaces of each of the constant‐density spheroids. However, the iterative equation used by Hubbard [ 2013 ] is based on an incomplete solution of the Poisson equation that diverges for high degrees of rotational flattening [ Kong et al. , 2013 ; Hubbard et al. , 2014 ]. An iterative approach for finding the shape of a single Maclaurin spheroid based on the full solution to the Poisson's equation was formulated by Kong et al. [ 2013 ] but was not extended to consider planetary structures with multiple layers and variable densities. Although the issue of non‐convergence is of little concern for slowly rotating planets such as Jupiter, for which the CMS model was designed, it is a significant issue for the rapidly rotating bodies that we consider in this work.

Post‐impact states were generated by simulations of giant impacts in the same manner as described in Ćuk and Stewart [ 2012 ]. Briefly, bodies were initialized in hydrostatic equilibrium by forming each body in isolation with isentropic silicate thermal profiles (section 2.3 ) with silicate specific entropy, S silicate =3.2 or 4 kJ K −1 kg −1 , corresponding to 1 bar potential temperatures of ∼1900 and ∼3300 K, respectively. Each impacting body was composed of 10 5 to 5 × 10 5 particles. Impact simulations were calculated for 24 to 48 h of simulation time, until the change in bound mass was negligible and the post‐impact structure reached a quasi‐equilibrium shape. The parameters for each simulation are given in Table S4 . Part of this suite of simulations (116 impacts) was calculated for Ćuk and Stewart [ 2012 ] and is primarily composed of small impactors (of mass M p ≤ M Mars ) onto targets with preimpact rotation. We extended the suite of giant impacts (by 46) to include other proposed Moon‐forming scenarios [ Canup and Asphaug , 2001 ; Canup , 2004 , 2012 ] and to sample the full range of giant impacts found in the later stages of planet formation [ Raymond et al. , 2009 ; Quintana et al. , 2016 ]. Our collection of impact outcomes span a range much broader than those proposed for the terminal impact event on Earth. The final bound masses range from 0.45 to 1.1 Earth masses ( M Earth ) but are typically between 0.8 and 1.1 M Earth .

Non‐rotating, isolated bodies were initialized with a given mass and an approximate isentropic thermal profile. The bodies were then equilibrated for 24 h with the entropy of the layers imposed at each time step and any particle velocity damped. To produce cold isolated bodies with a mantle entropy of 4 kJ K −1 kg −1 , a constant angular velocity was imposed to particles in a non‐rotating body of the same thermal structure. The structure was then equilibrated with the entropy of the layers imposed at each time step but without any damping of velocities. Rotating bodies with hotter thermal states were produced by starting from a colder planet of the same AM and increasing the specific entropy of the mantle at a rate of 0.25 kJ K −1 kg −1 d −1 until the desired specific entropy was reached. Isolated bodies had 10 5 particles, which is comparable to much of the recently published literature on giant impacts using the SPH technique.

The GADGET‐2 SPH code [ Springel , 2005 ], modified to include tabulated equations of state (EOS) [ Marcus et al. , 2009 ; Marcus , 2011 ], was used to calculate both the structure of isolated bodies and to simulate giant impacts. We restricted this work to Earth‐like planets, and all planets were modeled as differentiated bodies with 2/3 silicate and 1/3 iron by mass. As in prior work, the metal core and silicate were modeled as pure iron and pure forsterite, respectively. We used the M‐ANEOS model [ Melosh , 2007 ] to generate the tabulated EOS, which includes two condensed phases and vapor (refer to the supporting information of Canup [ 2012 ]).

We modeled the equilibrium structure of isolated bodies and the physical states obtained after giant impacts. Isolated bodies were modeled using both SPH simulations and the HERCULES code, a new method for modeling rapidly rotating bodies based on overlapping, concentric, constant‐density spheroids. Impact simulations were performed using SPH in a manner similar to prior studies of giant impacts [ Ćuk and Stewart , 2012 ]. In this section, we describe both methods and compare the results for isolated, corotating bodies.

The shapes of synestias can vary from bodies with pinched equators to bodies with flared disk‐like regions. Flaring typically only occurs at low pressures in the structure but can be substantial. The most extended structures in Figure 4 have a vertical scale height in the low‐pressure regions on the same order as the polar radius of the corotating region. The equatorial radius of synestias can extend beyond the Roche limit, which is the closest distance a satellite can withstand tidal forces from the planet (∼18 × 10 6 m for silicate satellites orbiting an Earth‐mass body).

We name structures beyond the CoRoL synestias , after the Greek syn for connected and Hestia for the goddess of architecture. The traditional definitions of mantle, atmosphere, and disk are not applicable in synestias (Figures 5 d– 5 i). In this work, we refer to the corotating portion of both super‐CoRoL and sub‐CoRoL structures as the corotating region and the Keplerian or sub‐Keplerian portion of structures as the disk‐like region or connected disk. The region between the corotating and disk‐like regions we call the transition region. In typical synestias, there is a monotonic angular velocity profile connecting the corotating and disk‐like regions (Figures 5 e and 5 h) and the transition region is very narrow. Synestias can be attained by a variety of celestial bodies and are not restricted to rocky bodies (section 6.8 ).

For the isolated SPH bodies shown in Figure 4 , the silicate particles were incrementally heated in an attempt to obtain as small a structure as possible (see section 2.1 ). Incremental heating limited the mass and AM transported into the sub‐Keplerian region of the structure during gravitational re‐equilibration. Other processes that drive a planet beyond the CoRoL, such as giant impacts, can introduce more mass and AM farther out in the structure.

We define the limit for planetary bodies with constant angular velocity as the corotation limit (CoRoL). The CoRoL is a surface that depends upon the mass, compositional layering, thermal profile, and AM of a body. Above the CoRoL, there is no solution for the planetary structure that is both hydrostatic and corotating. For bodies below the CoRoL, there is a unique solution for a perfectly corotating structure. However, above the CoRoL, there is no unique solution to the structure as it depends on the spatial distribution of mass and AM. The range of possible super‐CoRoL structures depends on the mechanism that drove the body beyond the limit.

The structures of rocky bodies with high angular momentum spans corotating, oblate bodies to a regime with a connected disk. Each body has an Earth mass with total angular momentum of 2.45and silicate specific entropies of (a–c) 4.0, (d–f) 5.5, and (g–i) 6.25 kJ Kkg. Structures were calculated using the GADGET‐2 SPH code with isentropic thermal profiles (class I) imposed. Figures 5 a, 5 d, and 5 g present pressure contours in a plane through the rotation axis. Details of the plotting technique are described in section 3 . The angular velocities (Figures 5 b, 5 e, and 5 h) and specific AM (Figures 5 c, 5 f, and 5 i) in the midplane of the silicate portion of the body are shown by black dots. Solid red lines denote the specific AM of a circular Keplerian orbit around a point mass, and dashed red lines indicate the specific AM for corotation with the inner region of the structure.

To illustrate the dynamics of extended planetary structures, we compare the rotational profiles of bodies with different silicate specific entropies at a total AM of 2.45 L EM . In Figure 5 , the pressure field and size of the structure are represented by the colored contours in the top row. Pressures were interpolated in the same manner as in Figure 4 . The radial angular velocity and specific angular momentum in the equatorial plane are shown by the black dots in the Figures 5 b, 5 e, and 5 h and Figures 5 c, 5 f, and 5 i, respectively. High‐AM bodies with low specific entropies have a constant angular velocity (i.e., are in solid body rotation) and consequently adopt very oblate equilibrium structures to conserve AM (Figures 5 a– 5 c). At higher specific entropies, the fraction of vapor in the outer regions is greater and the structure thermally expands. To maintain a fixed total AM, increasing the specific entropy is balanced by lower angular velocities. At a specific combination of AM and specific entropy, the angular velocity at the equator of the body intersects the Keplerian angular velocity. The bodies in Figures 5 a– 5 c and 5 d– 5 f are before and after this intersection. Beyond the intersection, it is not possible to attain an equilibrium corotating structure. The centripetal force required to remain bound and corotating is greater than the gravitational force, and a negative pressure gradient would be required to keep the equator at the same angular velocity as the interior of the body. A negative pressure gradient is non‐physical, and instead, the outer edge of the body must adopt a Keplerian or sub‐Keplerian angular velocity. With increasing specific entropy, the pressure gradient support of the structure increases, and the sub‐Keplerian region of the structure expands substantially (Figures 5 g– 5 i).

The equilibrium structures of rocky planetary bodies are highly variable. Each panel displays a one‐Earth‐mass body with an isentropic thermal profile (class I). Colors denote pressure contours in a plane through the rotation axis for axisymmetric bodies, with different angular momenta (rows) and silicate specific entropies (columns), calculated using the GADGET‐2 SPH code. Details of plotting technique are given in section 3 . Heating the body significantly inflates the radius, and increasing AM produces synestias, structures with connected disk‐like regions. Bodies to the left of the red line are corotating, and bodies to the right of the red line are above the CoRoL and are synestias.

We calculated the structure of rocky bodies with varying thermal profiles and AM. To illustrate the range of solutions, we determined the equilibrium structures for one Earth‐mass bodies for a range of AM and silicate isentropes (thermal profile class I) using GADGET‐2. We find that there are a wide variety of possible structures, ranging from compact to highly extended, as shown by the pressure contours in Figure 4 . Pressures were interpolated between SPH particles using a Delaunay triangulation. To avoid spurious contour lines in the extended regions with low particle density, the outermost few particles in each radial bin were not included when calculating the contours. The width of the radial bins and the number of particles excluded were varied depending on the structure being plotted. The gaps in some structures in Figure 4 are due to the lack of resolution in those areas. In hot, high‐AM structures the plotted pressure contours do not close as the scale height at the upper and lower surfaces of the structure is not resolved. Increasing the specific entropy of the silicate leads to an increase in the fraction of vapor in the low‐pressure regions of the body and a substantial increase in radius. At low AM, the structures are nearly spherical, oblate spheroids rotating with a single angular velocity, regardless of thermal state. However, bodies with higher AM do not always maintain a simple oblate spheroidal shape. The equilibrium structure for an isolated body of a fixed mass and composition can attain a large range of morphologies and sizes, depending on thermal state and total AM.

(a) The angular momentum required to exceed the CoRoL and (b) the angular velocity of bodies at the CoRoL for Earth‐like bodies of varying masses calculated using the HERCULES code. Symbols are for bodies with isentropic (class I, plus) and stratified (class III, triangle) silicate thermal profiles. Stratified profiles with 50% by mass cold lower silicate layer=3 kJ Kkgwere used. The lines are power law (a) or log linear (b) fits to the calculated points. Colors indicate the specific entropy of the silicates or upper silicate layer for the isentropic and stratified profiles, respectively. The black dashed line in Figure 8 a indicates the average AM expected for terrestrial bodies, extrapolated from the results of]. The solid black line indicates the range of planetary masses formed in the simulations of].

The CoRoL is also dependent on the mass of the body (Figure 8 ). We modeled the relationship between mass and the AM required to exceed the CoRoL for Earth‐like bodies by a simple power law of the form , where M is the mass of the body. The value of the exponent, γ , only varies slightly between different thermal structures. We fit the exponent for the three example thermal profiles shown in Figure 6 . For the 5.5 and 6.0 kJ K −1 kg −1 isentropic profiles (class I, solid lines), γ = 1.69 and 1.73, respectively. For the example stratified planet (class III) with 50 wt % at 6.0 kJ K −1 kg −1 (dash‐dotted line), γ = 1.77. These fits are good for larger mass planets but miscalculate the AM to exceed the CoRoL by 10 to 25% for 0.5 M Earth bodies and produce a misfit on the order of a few percent for M Earth bodies. The angular velocity of a body at the CoRoL also increases with mass and varies significantly for different thermal structures. The increase in angular velocity scales linearly with the logarithm of the mass of the body for the whole range of mass we considered.

For the expected range of values for Earth, the specific entropy of the iron core has little effect on the CoRoL for Earth‐like bodies. For the iron EOS model used here, the present thermal state of Earth's core corresponds to a specific entropy of ∼1.5 kJ K −1 kg −1 . The thermal history of the core is debated, so we considered a range of specific entropies, S core =1 to 2 kJ K −1 kg −1 , corresponding to temperatures of 700 to 11,000 K at the pressure of the present‐day core‐mantle boundary. This temperature range covers estimates for the early Earth's core that were calculated by requiring a core dynamo of its present strength for all of Earth history [e.g., O'Rourke and Stevenson , 2016 ]. A higher specific entropy core slightly increases the AM required to exceed the CoRoL because the increased moment of inertia of the extended core reduces the corotating angular velocity at a given AM. We find that for the range of core entropies considered, the equatorial radius of the structure and the AM required to exceed the CoRoL vary by only a few percent due to the low thermal expansivity of iron at high pressures.

The different thermal profiles have substantially different equatorial radii at the CoRoL (Figure 7 b) because of the strong sensitivity of the radius to the thermal state of the outermost layers. However, in terms of specific entropy, the CoRoL is similar for the different stratified thermal profiles (Figure 7 a). Consequently, the total angular momentum and specific entropy of the outer silicate layers are reliable metrics to determine if a particular stratified body is above the CoRoL.

We also considered bodies with stratified thermal profiles (class III) with a colder lower silicate layer ( S lower ) and hotter upper silicate layer ( S upper ). For S upper >∼5.5 kJ K −1 kg −1 , such structures cross the CoRoL at a lower AM compared to bodies with the vapor atmosphere profile (class II) with a silicate specific entropy equal to S upper . The cold lower silicate layer is more compact and closer to the rotation axis than the equivalent mass fraction of bodies with vapor atmosphere profiles or isentropic profiles with high specific entropies. Accordingly, the structure has a lower moment of inertia, and the cold lower silicate layer of a stratified body does not accommodate as much AM for a given angular velocity. Thus, stratified bodies must rotate faster, compared to an isentropic body with the same AM, and the equator of the planet intersects the Keplerian orbit at a lower total AM. This effect increases with a larger difference in specific entropy between the upper and lower silicate layers. The CoRoL does not turn over for thermally stratified structures in the range of specific entropies considered because the cold dense cores of such bodies have a comparatively low moment of inertia. Because the isentropic case requires a higher specific entropy to reach the CoRoL compared to stratified structures, the suite of isentropic structures in Figure 4 has fewer synestias over this range of AM and specific entropy than would be achieved by stratified bodies. The fact that thermal stratification causes bodies to cross the CoRoL at lower AM is significant because giant impacts create stratified structures (section 5 ).

The CoRoL depends on the thermal structure of a body. Using the HERCULES code, we calculated the CoRoL for bodies with isentropic (black line, class I), vapor atmosphere (red line, class II), and stratified silicate (class III) thermal profiles. We present the CoRoL in terms of (a) the silicate specific entropy used to define the thermal profile and (b) the equatorial radius as a function of angular momentum. For stratified profiles only the specific entropy of the upper silicate layer was varied, and the lower layer was held constant with=4 kJ Kkg. Stratified profiles (section 2.3 ) with 25% (orange lines) and 50% (blue lines) by mass upper silicate layers are shown.

The CoRoL is also sensitive to the thermal profile within the body (section 2.3 ). Figure 7 presents the CoRoL in terms of (a) the specific entropy of the silicate layer (upper layer for stratified profiles (class III)) and (b) the 10 bar equatorial radius, as a function of AM. The AM required to exceed the CoRoL with an isentropic profile (class I) is offset from the other thermal profiles because the density profile assumes an ideal mixture of liquid and vapor in the mixed phase region. Vapor atmosphere profiles (class II), with pure gas densities in the lowest pressure layers, produce bodies with larger equatorial radii and reach the CoRoL at lower AM for the same value for S lower . For low specific entropies, the CoRoL for bodies with isentropic (class I) and vapor atmosphere (class II) thermal profiles tend toward each other because the mixed phase region of the structure is absent or negligible. At high specific entropies, the CoRoL for bodies with isentropic (class I) and vapor atmosphere (class II) thermal profiles converge as the mixed phase regions of isentropic bodies become increasingly dominated by vapor. At these high specific entropies, the CoRoL turns over as the slower rotation rate of the more extended structures begins to overcome the effect of an increased equatorial radius.

The radius of Earth‐like bodies varies strongly with specific entropy and angular momentum. Each colored line is the range of equatorial radii of bodies with isentropic thermal profiles (class I) of fixed silicate specific entropy and varying angular momenta, calculated using the HERCULES code. The black line denotes the limit for equilibrium corotating structures (CoRoL). Above the CoRoL, structures are synestias and a portion of the body must form a disk‐like region.

The equatorial radius, and hence the CoRoL, is very sensitive to the specific entropy of the outer regions of the body and the total AM. Figure 6 presents the calculated equatorial radius for Earth‐mass bodies with isentropic silicate layers (thermal profile class I), where each line indicates the radius for a constant specific entropy and varying total AM. Each line terminates at the limit of corotating equilibrium structures, and the locus of such points defines the CoRoL (black line) for this specific planetary mass, compositional layering, and thermal profile. Hotter bodies have significantly expanded equatorial radii and cross the CoRoL at lower AM. Even below the CoRoL, the equatorial radius varies by a factor of 3 depending on thermal state and AM.

We have determined the CoRoL as a function of mass, AM, and thermal state for Earth‐like planets using the HERCULES code. As the current version of HERCULES cannot calculate the structure of non‐corotating bodies, including synestias, we find the CoRoL by extrapolation from corotating bodies. In AM increments of 0.05 or 0.01 L EM , we calculated the equilibrium structure of corotating bodies with a given thermal profile, compositional layering, and mass. Above a certain AM, close to the CoRoL, HERCULES is no longer able to find a physical solution to the structure as, during at least one iteration, the total potential gradient at the edge of the body changes sign. This causes the equipotential surfaces to cross at higher latitudes, breaking the assumptions of the model. In order to find the CoRoL, we linearly extrapolate to higher AM to find the point at which the corotating and Keplerian angular velocities at the edge of the structure intersect (Figure S8 ). This point is the CoRoL and is generally close to the highest AM structure found by HERCULES.

The prevalence of both hot and high‐AM states for terrestrial bodies during accretion makes the formation of synestias highly likely. In Figure 11 , we find that thermal states with dominantly vapor outer layers (e.g., S upper > S crit ) reach the CoRoL for AM greater than about 1.5 L EM . Kokubo and Genda [ 2010 ] found that a majority of Earth‐mass planets exceed this AM at the end of accretion. Our examination of impact energies during accretion find that nearly all Earth‐mass planets experience one or more collisions that produce dominantly vapor outer layers. Therefore, we conclude that it is common for Earth‐mass bodies to exceed the CoRoL and form synestias one or more times during accretion.

The rotational states of terrestrial bodies during accretion is difficult to track, especially because partitioning of AM between the growing bodies, collision ejecta, and a planetesimal population is poorly understood. N ‐body studies have found that each giant impact can dramatically alter the spin state of a growing planet [e.g., Agnor et al. , 1999 ]. Kokubo and Genda [ 2010 ] have conducted the best assessment to date of the spin state of rocky planets during accretion. They used an N ‐body simulation of planet formation with bimodal impact outcomes, either perfect merging or hit‐and‐run, and tracked the AM of each of the bodies in the simulation. Kokubo and Genda [ 2010 ] found that the mean AM of rocky planets is large, e.g., 2.69 L EM for Earth‐mass planets and that the distribution of AM is wide. The mean angular velocity of final planets was similar for the entire mass range in their study, covering planetary masses between about 0.25 and 1.75 M Earth . We infer that that the AM of terrestrial bodies roughly scales as L = 2.69 L EM ( M / M Earth ) 5/3 (black line in Figure 8 a). These results indicate that the AM of planets during accretion is expected to be large and that the present‐day AM of planets have been modified by satellite and solar tides.

N ‐body simulations of terrestrial planet formation have recently incorporated more realistic collision outcome models including fragmentation [e.g., Chambers , 2013 ; Carter et al. , 2015 ; Quintana et al. , 2016 ]. Quintana et al. [ 2016 ] used the same modified specific impact energy, Q S , to evaluate the collisions in their simulations and graciously provided us with the collision histories for the final planets. We calculated the number of high‐energy impacts that near Earth‐mass planets experience as they accrete. We find that most near Earth‐mass planets experienced several giant impacts with Q S >2 × 10 6 J kg −1 during accretion and over half also experienced at least one impact with Q S >5 × 10 6 J kg −1 (Figure 13 ). Assuming that the scaling of thermal state with impact energy holds for smaller bodies (Figure 12 ), we expect rocky bodies to form highly vaporized structures multiple times during accretion. Other N ‐body studies without fragmentation and with different initial conditions and alternative configurations of the giant planets [e.g., O'Brien et al. , 2006 ; Hansen , 2009 ; Raymond et al. , 2009 ; Walsh et al. , 2011 ; Levison et al. , 2015 ] all produce high‐energy impacts in the final stages of accretion although the frequency of different energy impacts may change between these different scenarios.

The specific entropy of the upper silicate layers of post‐impact structures scales well with the specific impact energy,. The scaling for both the specific entropy of the upper (a) 50 wt % and (b) 25 wt % is shown. Symbols indicate the dynamical group of the structure: sub‐CoRoL (square), super‐CoRoL (circle), and co‐CoRoL (triangle). The colors indicate the average specific entropy of the corresponding lower 50 wt % (Figure 12 a) or 75 wt % (Figure 12 b) of the silicates. Filled symbols indicate structures with little or no mass in the disk‐like region of the structure. Only post‐impact structures with bound masses between 0.9 and 1.1are shown. Red symbols correspond to the example cases plotted in Figure 9

We find that the specific entropy of the outer silicate layers scales linearly with the logarithm of specific impact energy for post‐impact bodies with masses between 0.9 and 1.1 M Earth and Q S >∼10 6 J kg −1 (Figure 12 ). Q S does not account for any entropy increase caused by secondary impacts in graze and merge collisions or by reimpacting debris. Such effects can be significant, particularly for low‐energy impacts, and are likely responsible for the apparent plateau in entropies at low Q S in our suite of impacts as these low energies are dominated by graze and merge, canonical‐style Moon‐forming impacts. Impacts with Q S >2 × 10 6 J kg −1 generally deposit sufficient energy such that the upper 25 wt % of the silicates have an average specific entropy that exceeds the critical point value for the equation of state, S crit =5.4 kJ K −1 kg −1 (Figure 12 b, solid line). For impacts with Q S >5 × 10 6 J kg −1 , the upper 50 wt % of silicates typically attain mean specific entropies above the critical point (Figure 12 a, dashed line). Such hot post‐impact structures are dominantly vapor at low pressures. They have much larger radii compared to entirely condensed planets of the same mass and AM. The correlation between post‐impact thermal state and specific impact energy provides a straightforward means to determine which collisions during accretion lead to highly vaporized structures.

Next, we examine the likelihood of generating hot, rapidly rotating post‐impact states during the giant impact stage of terrestrial planet formation. First, we examine the specific entropy of post‐impact structures. We find that the specific entropy of the outer silicate portions of post‐impact structures scales well with a modified specific impact energy, Q S . Q S is a variation of the parameter developed in Leinhardt and Stewart [ 2012 ] that adjusts for the geometry of collisions between similarly sized bodies to estimate the relative deposition of impact energy into the post‐impact body for different collision scenarios ( supporting information Text S4 ).

A significant implication of this analysis is that some post‐impact structures cannot be analyzed as a separate planet and disk, as is common in studies of post‐impact states. Instead, the post‐impact structures are synestias, which must be treated as a single extended structure. The impact‐generated synestias have a large variation in the mass of the disk‐like region. The proportion of mass and AM in the disk‐like region is determined by the impact conditions [with consistent results between different types of codes, e.g., Canup et al. , 2013 ]. Because these structures exceed the CoRoL, any process that transfers mass from the disk‐like region to the corotating region (and preserves the mean specific entropy and AM) would not be able to transfer all the mass to the corotating region. Or, more colloquially, the entire disk‐like region cannot fall down. Because the structure is expected to cool by radiation (perhaps mitigated by internal dynamics such as viscous heating and continuing accretion), an impact‐generated synestia will evolve to a state below the CoRoL, but the timescale will depend on the evolutionary path for that particular body. This investigation of post‐impact states is a snapshot of the (near) initial structure after the event.

Next, we consider the structure of post‐impact bodies in the context of proposed Moon‐formation scenarios. The example cases shown in Figures 9 are examples of proposed Moon‐forming giant impacts. Figures 9 a– 9 d show a typical structure formed by an impact in the style of the canonical Moon‐forming impact [e.g., Canup and Asphaug , 2001 ; Canup , 2004 , 2008b ]: sub‐CoRoL structures with relatively cold disk‐like regions. Figures 9 e– 9 h and 9 i– 9 l show example high‐AM, high‐energy Moon‐forming impacts from Ćuk and Stewart [ 2012 ] and Canup [ 2012 ], respectively. These high‐AM impacts tend to produce super‐CoRoL or co‐CoRoL structures with substantially vaporized disk‐like regions. Note that both the high‐AM studies proposed a range of impact energies and impact parameters as potential Moon‐forming events. Previous studies of post‐impact states after the canonical impact [e.g., Canup and Asphaug , 2001 ; Canup , 2004 , 2008b ] did not observe synestias as they only considered bodies with an AM similar to the present‐day Earth‐Moon system. Such bodies would require very high specific entropies in their outer layers to exceed the CoRoL (Figure 11 ). Such hot thermal states are not reached in canonical style impacts.

Giant impacts can produce structures that are above the CoRoL. The characteristics of the post‐impact structures can be assessed by their total angular momentum and average specific entropy of the outer (a) 50% or (b) 25% by mass of silicates. The symbols indicate the dynamical group of the structure: sub‐CoRoL (square), super‐CoRoL (circle), and co‐CoRoL (triangle). The colors indicate the average specific entropy of the corresponding lower 50% (Figure 11 a) or 75% (Figure 11 b) by mass of silicates. Filled symbols indicate structures with little or no mass in the disk‐like region of the structure (see text for discussion of such cases). Only impacts that have a total bound mass of the post‐impact structures between 0.9 and 1.1are plotted, and the lines denote the CoRoL for Earth‐like composition bodies of 0.9 (dashed), 1.0 (solid) and 1.1 (dotted). The CoRoL lines indicate the upper silicate layer specific entropy and AM limit for stratified thermal profiles (=4 kJ Kkg) with 50 wt % (Figure 11 a) or 25 wt % (Figure 11 b) upper silicate layers. Red symbols correspond to the example cases plotted in Figure 9

In Figure 11 , we compare the post‐impact structures to the CoRoL for stratified thermal profiles. As shown in Figure 7 , the CoRoL is defined by similar AM and outer specific entropy for a variety of stratified thermal profiles. The majority of our post‐impact states that we identified as super‐CoRoL or co‐CoRoL exceed the CoRoL calculated for stratified bodies when considering the specific entropy of the outer 50 or 25 wt % of the body. All our sub‐CoRoL structures fall below the CoRoL for stratified bodies. There are some structures that we identified as super‐CoRoL structures that fall below the CoRoL plotted in Figure 11 . These structures typically have thermal profiles where the top few percent of the mass has very high specific entropies. As a consequence, the structures intersect the CoRoL but with very little mass in the disk‐like region (filled circles in Figure 11 ). The existence of very hot outer regions is not captured in the average entropy of the structure and comparing such structures with the CoRoL for stratified planets with 50 or 25 wt % hot outer layers is not a good indicator of whether the structure is above the CoRoL.

The approximate vapor structure for the three example post‐impact structures shown in Figure 9 . The post‐impact states were processed by removing condensed material and re‐equilibrating the vapor structure as described in section 5.1 . (a, d, and g) Pressure and (b, e, and h) entropy contours were plotted as described in section 3 . (c, f, and i) Angular velocity for the silicate SPH particles in the midplane. Solid red lines denote the angular velocity of a circular Keplerian orbit around a point mass. Pressure and entropy scale bars in this figure differ from elsewhere in the paper.

We would like to compare the shape of post‐impact structures to those of the isolated synestias; however, the inability of impact codes to correctly model phase separation makes it difficult to calculate the pressure field that would be expected immediately after an impact. To approximate the vapor pressure field in each of our example cases, we have modified the post‐impact states produced by SPH to remove the condensed material. We removed unbound SPH particles and the condensed mass fraction of bound SPH particles in the disk‐like region from the simulation at 24 h after the impact. The structure was then re‐equilibrated for a further 24 h with the mass of any condensate removed at each time step. The resulting structure approximates the continuous vapor pressure field at several dynamical times after the impact if the condensate was fully decoupled from the vapor and the effects of cooling were neglected. The calculated pressure fields, and corresponding specific entropy fields, for the three examples cases are given in Figure 10 . The pressure fields and angular velocity profiles for the super‐CoRoL and co‐CoRoL cases are similar to those found for isolated synestias (Figures 4 and 5 ). Impacts emplace significant amounts of mass and AM far from the rotation axis and the vapor structure can extend beyond the Roche limit, with significant vapor pressure (tens of bars) at several Earth radii. The larger masses in the disk‐like regions and the increase of specific entropy with decreasing pressure in the post‐impact states (Figures 10 b, 10 e, and 10 h) leads to more exaggerated flaring than in the isentropic isolated bodies in Figures 4 and 5 . The concentration of cold material in the midplane is due to buoyancy driven settling and is partly a result of the lack of thermal equilibration in the SPH code. For the particular example of a co‐CoRoL case in Figure 9 , the structure after condensate extraction is above the CoRoL and the resulting angular velocity profile structure becomes monotonically decreasing with radius. The sub‐CoRoL example also shows a flared pressure field similar to the other cases but with lower pressures. The angular velocity profile in the sub‐CoRoL case still has a shear boundary after removing the condensate, and the profile is similar to the expected solution for a shearing fluid [ Chandrasekhar , 1969 ; Desch and Taylor , 2013 ].

We examined whether post‐impact structures are below or above the CoRoL. Sub‐CoRoL structures (e.g., Figures 9 a– 9 d) do not have sufficient AM for the hot, inner regions of the structure to extend out to the intersection between the corotating and Keplerian angular velocities (black dashed line in Figure 9 ). Consequently, the vapor at the edge of the hot inner regions is sheared as the angular velocity profile transitions from corotating to sub‐Keplerian. The disk‐like regions of sub‐CoRoL structures tend to have a low vapor fraction in the midplane. In SPH post‐impact, sub‐CoRoL structures there are colder particles in the transition region that contribute to the shear, but the angular velocity profile is still dominated by the shear in the vapor. Conversely, in super‐CoRoL structures (synestias), the hot inner regions extend out to the intersection between the corotating and Keplerian angular velocities. There is a smooth transition in the angular velocity of the vapor from the corotating inner region to the sub‐Keplerian disk‐like region. In some super‐CoRoL structures, there are colder particles in the transition region that have different angular velocities than the vapor, leading to an apparent shear. However, this shear does not significantly affect the overall smooth angular velocity profile. In a few cases, it is difficult to determine whether the structure is above or below the CoRoL. An example of such a structure is shown in Figures 9 i– 9 l. The vapor of the inner region only just extends to the intersection between the corotating and Keplerian angular velocities, and there is a substantial mass of colder particles in the transition region. As a result, the mass averaged angular velocity profile has a small shear boundary. It is not clear from the SPH post‐impact structure alone whether such structures are above the CoRoL nor whether they would remain so following thermal equilibration. We classify such structures as co‐CoRoL as they are at the border between sub‐CoRoL and super‐CoRoL states.

The thermal state of the silicates farther out in the structure is more variable. Some structures have disk‐like regions that have high specific entropy and are mostly vaporized even to beyond the Roche limit (Figure 9 g), while others have disk‐like regions that are much colder and have a large mass fraction that is condensed (Figures 9 c and 9 k). In the latter case, the transition between the high‐entropy inner regions and the colder disk‐like regions typically occurs in the dynamical transition region. The parcels of material represented by SPH particles cannot thermally equilibrate with each other. As a result, colder, more condensed particles can coexist with hotter, mostly vapor particles in the transition region. The colder particles are less pressure supported than the vaporized particles and rotate at a lower angular velocity at the same radius. These colder particles complicate the dynamics in the transition region for the SPH post‐impact structures (as shown in Figure 9 ). In reality, the material represented by the colder and hotter SPH particles would thermally equilibrate, and in most cases in our suite of post‐impact states, the transition region would evolve to be dominantly vapor. In all cases, there may be some (comparatively colder) material in orbit at larger distances (e.g., a distal debris disk) that is not continuous with the central structure.

Impact‐induced heating is very heterogeneous and post‐impact structures are far from isentropic. The inner regions of the structure are thermally stratified with a monotonically increasing specific entropy with radius (Figures 9 c, 9 g, and 9 k). This general profile is produced by buoyancy and rotational forces during the impact acting to impose a monotonically decreasing density profile. Such gravitationally equilibrated profiles neglect any effects from chemistry, thermal equilibration, or shear strength that may act on the dynamical timescale of hours. For impacts between initially condensed bodies, the post‐impact lowermost mantle is a high‐pressure liquid or a liquid‐solid mixture (see discussions in Stewart et al. [ 2015 ] and Nakajima and Stevenson [ 2015 ]) with specific entropies typically between 3 and 5 kJ K −1 kg −1 . The specific entropy of the silicates in the inner regions increases substantially with radius (Figures 9 c, 9 g, and 9 k), typically far exceeding the specific entropy of the critical point ( S crit =5.40 kJ K −1 kg −1 ). In most cases we considered, the corotating region does not intersect the liquid‐vapor phase boundary. At high pressure in the body, the silicate grades smoothly from a liquid into a supercritical fluid and then vapor. In the cases where the thermal profile does intersect the liquid‐vapor phase boundary for a range of pressures within the corotating region, the condensate fraction is usually small.

Giant impacts produce a range of post‐impact structures, with varying dynamic and thermodynamic properties (Table S4 ). Dots represent silicate SPH particles in the midplane for three example impact scenarios (one in each column) that produce substantial disk‐like regions at 48 h after first contact. The dynamical structures identified as: (a–d) sub‐CoRoL, (e–h) super‐CoRoL, and (i–l) co‐CoRoL structures. Solid red lines denote the angular velocity or specific AM of a circular Keplerian orbit around a point mass. Dashed red line indicates the specific AM of material corotating with the inner region. The black particles are either at pressures above the critical point (25.5 kbar) or are pure vapor, and the blue particles are mixtures of liquid and vapor. The super‐Keplerian particles in Figures 9 a and 9 b are the remnants of a disrupting moonlet.

We find that there is a broad range of post‐impact structures, in terms of both thermal and dynamical state. Three example post‐impact structures with substantial disk‐like regions are shown in Figure 9 . All structures have a corotating region close to the rotational axis and a Keplerian or sub‐Keplerian disk‐like region farther out. The transition region varies substantially between different structures. A large fraction of our suite of post‐impact states have very narrow transition regions, and the corotating region grades smoothly into the disk‐like region, as is the case for isolated synestias (Figure 5 ). But in some cases there is a shear boundary (an increasing angular velocity with radius) in the transition region which varies in magnitude and width (Figures 9 a and 9 i). The prevalence of structures with narrow transition regions is likely biased by the many high‐AM impacts with fast, small impactors in our suite of impacts. The nature of the transition region is determined both by the AM and thermal state of the structure, as we discuss below.

Rocky bodies are naturally forced into hot, rotating states by giant impacts during the later stages of accretion. We seek to understand the significance of hot, rotating planetary structures by assessing the range of possible structures that could be attained during planet formation and the frequency of occurrence. First, we examine and attempt to classify the dynamical and thermal structures of post‐giant impact states. In order to do this, we calculated the outcome of impacts between Earth‐like bodies with a range of collision parameters using SPH, as described in section 2.1 , and combined our results with the simulations from Ćuk and Stewart [ 2012 ]. We examined the instantaneous structures of post‐impact structures at 24 to 48 h after first contact, when the structure had reached a quasi‐steady mass and AM distribution. This suite of impacts covers a range of impact energies but focuses on those that have high enough specific impact energy to be described as giant impacts by Quintana et al. [ 2016 ] (section 5.2 ). The full set of impacts are summarized in Table S4 .

6 Discussion

6.1 Dynamic Stability of Hot, Rotating Planetary Structures Rapidly rotating, axisymmetric bodies can be susceptible to dynamic instabilities, such as the bar instability that causes the breakdown of axial symmetry and drives reorganization of the structure [e.g., Chandrasekhar, 1969]. Constant‐density, rotating spheroids are also known to transition from axisymmetric Maclaurin spheroids to non‐axisymmetric Jacobi spheroids above a critical degree of rotational flattening, because Jacobi spheroids are lower energy states [Chandrasekhar, 1969]. If rapidly rotating axisymmetric planetary structures were susceptible to such instabilities, the structure could rearrange or break up. Here we examine the stability of the rapidly rotating planetary structures considered in this work. The stability of rapidly rotating fluid bodies is typically evaluated by the ratio of the rotational kinetic energy, T, to the gravitational potential energy, W. Rotating bodies are unstable to non‐axisymmetric modes if the ratio T/|W| is above a critical value, which varies depending on the body in question. Constant‐density Maclaurin spheroids [Chandrasekhar, 1969], stellar systems [Ostriker and Peebles, 1973], rapidly rotating white dwarfs [Ostriker and Tassoul, 1969], and polytropic stars [Ostriker and Bodenheimer, 1973] have critical values in the range 0.14 to 0.27. All the pre‐CoRoL and super‐CoRoL structures considered in this work have an energy ratio T/|W| below this range of critical values. For example, the highest value for the structures shown in Figure 4 is 0.07. For the post‐impact structures we studied for this paper (Table S4) the maximum value is 0.085 and the average value is 0.042. We also did not observe any instabilities or triaxial structures in our isolated body SPH calculations, where the equilibrium shapes were calculated over many dynamical timescales. Therefore, based on the energy stability criteria and our empirical observations, the rapidly rotating planetary structures presented in this work are likely to be dynamically stable.

6.2 Lifetimes of Hot Planetary Structures Although hot, extended planetary structures are dynamically stable, they evolve by radiating energy. In some cases, AM may be reduced by tides or resonant interactions. For synestias, a sufficient decrease in AM or thermal energy will cause the body to fall below the CoRoL and adopt a more compact, corotating structure. Bodies that are far enough away from their host stars will eventually cool to a silicate magma ocean overlain by a volatile‐dominated atmosphere. Calculating the timescales for the collapse of synestias and the cooling of post‐impact states to magma oceans is key to understanding the influence of post‐impact states on the evolution of terrestrial planets. Synestias formed by impacts are substantially vapor in their disk‐like regions. Unless the body is very close to its parent star, the outer regions of the structure are expected to cool quickly by radiation. When the photosphere is controlled by silicate condensation, the radiative temperature is high [∼2300 K, Lock et al., 2016], but the energy that needs to be radiated to fall below the CoRoL is substantial. Condensing the vapor requires significant energy loss, both to cool to the phase boundary and to release the latent heat of vaporization of silicates. The contraction of the structure upon cooling also releases a comparable amount of potential energy. The amount of energy that needs to be radiated to fall below the CoRoL depends on how the internal energy is redistributed during cooling because the thermal structure affects the location of the CoRoL. We can estimate a lower limit on the timescale to cool below the CoRoL by considering just the potential energy difference between a post‐impact structure and a corotating body of the same AM, with a thermal state slightly below the CoRoL. This approach neglects the loss of energy required to reduce the specific entropy of the outer layer of silicates (for a fixed AM) to below the CoRoL, which could be substantial. We assume a radiative temperature of ∼2300 K for silicate vapor. However, the surface area of the structure will change with time, and the difference in radiative surface area between the post‐impact state and the sub‐CoRoL structure is an order of magnitude. As a result, there is a substantial uncertainty in cooling time. Considering the full range of possible radiative surfaces, typical Earth‐mass post‐impact synestias would require a minimum time of order 101 to 103 years to cool to a state below the CoRoL. Alternatively, a synestia could be brought below the CoRoL by removing AM from the structure by tides, interactions with the host star [Murray and Dermott, 1999], or three body interactions with the host star and a moon [Goldreich, 1966; Touma and Wisdom, 1994; Ćuk and Stewart, 2012; Wisdom and Tian, 2015; Ćuk et al., 2016; Tian et al., 2017]. A significant reduction of AM by tides takes 103–109 years, depending on the mechanism, which is generally greater than the radiative cooling timescale. AM may also be redistributed by viscous spreading, but the whole structure would remain above the CoRoL (for a fixed thermal state). Thus, an impact‐generated Earth‐mass synestia would be stable for greater than tens to thousands of years. The timescale for a post‐impact structure to cool to a fully condensed magma ocean is longer. To estimate the timescale, we calculate the difference in potential, kinetic, and internal energy between post‐impact states and corotating, isentropic condensed planets of the same AM (first column in Figure 4). We again assume a radiative temperature of ∼2300 K and consider radiative surface areas ranging from the post‐impact state to a fully condensed body. For highly extended post‐impact states that are above the CoRoL, the cooling time to a magma ocean state is on the order of 102 to 103 years. For post‐impact states produced by less energetic impacts that do not begin above the CoRoL, the timescale is shorter but still on the order of 102 to 103 years. For all partially vaporized states, the cooling time may be substantially extended if the chemistry of the system produces a lower temperature photosphere. For example, refractory components in the photosphere may form relatively low temperature hazes. The chemistry of the silicate atmosphere will also evolve as it cools, potentially changing the effective radiative temperature. In addition, the energy input from continued accretion will increase the cooling timescale of post‐giant impact structures, especially if the accreting bodies were small and primarily deposited their energy in the high‐entropy outer layers of the body. The lifetime of hot, extended structures affects the accretion and evolution of terrestrial planets (section 6.6).

6.3 CoRoL Versus Solid Body Rotational Breakup The rotational breakup of small rocky bodies has been extensively studied due to its importance in bounding the possible rotational states of minor planets and in understanding the origin of multiple asteroid systems [e.g., Harris, 1996; Pravec and Harris, 2007; Richardson and Walsh, 2006; Holsapple, 2004; Richardson et al., 2005; Ćuk, 2007; Walsh et al., 2008] [see review by Walsh and Jacobson [2015]]. Solid bodies with a range of material strength properties have been considered. Cohesionless asteroids share similarities with fluid bodies, and the spin stability limit for self‐gravitating solid bodies with deformation has been studied in detail [e.g., Holsapple, 2004]. The CoRoL for the larger rocky bodies studied in this work is different than the critical spin limit for small bodies. The typical assumption of incompressibility used in the small body literature does not hold for planet‐sized bodies. In addition, planetary bodies with hot thermal states intersect the liquid‐vapor phase boundary, leading to a significant decrease in average density. Accordingly, the CoRoL is typically reached at a lower AM compared to the critical spin limit for a rigid, condensed body. The outcome of exceeding spin stability is also different for asteroids and hot planets. Upon exceeding the critical spin limit, an asteroid will shear into two or more bodies with relative motions that conserve AM [Ćuk, 2007; Walsh et al., 2008; Walsh and Jacobson, 2015]. In contrast, a hot planetary body that exceeds the CoRoL can remain a single structure with the excess AM forming a disk‐like region. Therefore, the CoRoL is a new dynamical transition that is fundamentally different from the well‐studied critical spin limit for smaller solid bodies.

6.4 Analysis of Post‐impact Structures The structure and evolution of post‐impact states can have a substantial impact on planet formation and the final properties of terrestrial planets (section 6.5). In particular, the properties of the disk‐like regions of post‐impact structures control the mechanisms and efficiency of satellite formation (section 6.6). The numerical methods used to model giant impacts cannot be used to directly model satellite formation. Hence, dedicated disk structure and evolution models have been used to study the formation of moons from the disk‐like regions of post‐impact structures, with initial conditions based on the results of impact simulations. Informed by work on astrophysical disks, most studies treat post‐impact structures as consisting of a distinct planet and disk, both for modeling the disk evolution and for analyzing post‐impact states. However, as discussed in section 5 and shown in Figure 9, above the CoRoL the entire post‐impact structure can be continuous between the corotating and disk‐like regions. Even below the CoRoL, there may be significant vapor pressure in the transition region between the corotating and disk‐like regions (e.g., Figure 10a), and the disk‐like region cannot be treated in isolation from the rest of the structure. Here we summarize previous work on the structure of circumterrestrial disks formed by giant impacts and identify some issues with common approaches taken in these studies. Published studies of giant impacts have made different assumptions on how to divide post‐impact structures into a planet and a disk. Since most impact studies have utilized SPH, in this section we refer to a parcel of material as a particle. The majority of studies [e.g., Canup et al., 2001; Canup and Asphaug, 2001; Canup, 2004, 2008b, 2012; Nakajima and Stevenson, 2014, 2015] have used an iterative routine to divide the structure, e.g., as described in Canup et al. [2001]. First, an initial guess is made to estimate the mass and equatorial radius of the planet. Particles are defined as being part of the planet or disk based on their AM. A particle is classified as being in the disk if it has sufficient AM such that the semi‐major axis of a circular Keplerian orbit with the same AM is greater than the planet's equatorial radius. Bound particles that are not classified as being in the disk are considered to be part of the planet. The equatorial radius of the planet is recalculated based on the total mass and AM of all the particles classified as being in the planet, assuming that the planet has a similar bulk density to the present‐day Earth and limited rotational flattening. The particles are then reclassified based on the new mass and radius of the planet, and the procedure is repeated until convergence. Alternatively, Ćuk and Stewart [2012] used a density contour of 1000 kg m−3 to define the planet's equatorial radius due to the large rotational flattening of the planet. As in the iterative routine, material with sufficient AM to be in a circular Keplerian orbit above the equator of the planet was considered to be part of the disk. Low‐density material that did not meet the AM requirement to be in the disk was considered to be part of a vapor atmosphere around the planet. The combined mass of the planet and atmosphere was used for calculating the orbit of particles. In both methods, once the structure is divided into a planet and disk, the surface density of the disk is calculated by moving the mass of each disk particle to the semi‐major axis of the equivalent circular Keplerian orbit for the AM of that particle. Although the details differ slightly, previous work applied the same principle of dividing the structure based on AM, and we refer to both methods as the conventional analysis. The conventional analysis neglects the connection between the disk‐like region and the corotating region of the structure through the vapor phase. Also, the radius of the planet defined in the conventional analysis does not account for the large mass of high specific entropy, low‐density material in the corotating and transition regions of the structure because the planet is assumed to have a bulk density corresponding to condensed material (of order 103 kg m−3). The assumptions made in analyzing post‐impact states are important as the properties of the post‐impact structure (e.g., the disk mass, AM, and surface density) calculated using the conventional analysis and reported in impact studies are used to inform the initial conditions of disk evolution models, either directly or indirectly. Most studies of disks in the aftermath of giant impacts have focused on the Moon‐forming event [e.g., Thompson and Stevenson, 1988; Ida et al., 1997; Kokubo et al., 2000; Ward, 2012, 2014, 2017; Salmon and Canup, 2012, 2014; Nakajima and Stevenson, 2014; Canup et al., 2015; Charnoz and Michaut, 2015]. In the canonical giant impact, ∼20 wt % of the disk would have been initially vapor [Canup and Asphaug, 2001; Canup, 2004], but in the recently proposed high‐AM models, the disk would have been largely vaporized [Ćuk and Stewart, 2012; Canup, 2012; Nakajima and Stevenson, 2014]. Modeling the structure and evolution of partially vaporized disks is challenging since the material in the disk is composed of multiple phases and chemical components. It is necessary to make a number of simplifying assumptions to make the problem tractable. The simplest disk models are N‐body simulations that neglect the vapor in the disk and treat the disk material as fully condensed particles [e.g., Ida et al., 1997; Kokubo et al., 2000]. As the disk evolves under gravitational forces, material with an orbit that intersects the planet is assumed to be lost from the disk and is removed from the simulation. More advanced disk models have made approximations of the multiphase physics [e.g., Thompson and Stevenson, 1988; Ward, 2012, 2014, 2017; Salmon and Canup, 2012, 2014; Nakajima and Stevenson, 2014; Canup et al., 2015; Charnoz and Michaut, 2015]. Typically, particles of condensate are assumed to experience mutual collisions that damp the eccentricity and inclination of their orbits. Outside the Roche limit, the vapor is assumed to condense quickly and the condensate rapidly collides together to form moonlets. Inside the Roche limit, the condensate separates from the gas and forms a liquid layer [Ward, 2012] or froth [Thompson and Stevenson, 1988] in the midplane of the disk, surrounded by a vapor atmosphere. Viscosity in the disk causes the inner disk to spread and the surface density distribution to evolve. As in N‐body simulations, material that spreads inside the radius of the planet is assumed lost from the disk, although some preliminary work has considered a more complex inner boundary [Desch and Taylor, 2013; Charnoz and Michaut, 2015]. Although the above assumptions are the most important for our discussion here, a number of other assumptions are made in different models of the inner disk, including a constant surface density [e.g., Salmon and Canup, 2012, 2014], neglecting the radial pressure gradient [e.g., Charnoz and Michaut, 2015], and using a single‐component equation of state [e.g., Ward, 2012]. All the disk models are initialized with conditions that are based on the conventional analysis described above and evolve under the assumption of unfettered mass flow to the planet. Our work identifies several issues with the assumptions made in previous studies of post‐impact structures and circumterrestrial disks. First, post‐impact structures cannot always be divided into a planet and disk that can be modeled separately. For synestias, the corotating region and the disk‐like region are continuous, as demonstrated by the angular velocity and thermodynamic profiles of the structure shown in Figures 9e–9h. In such cases, calculating the structure or evolution of a region of the structure without considering the whole will lead to substantial errors. Even for sub‐CoRoL structures, there can be significant pressures (approximately thousands of bars in the transition region between the corotating and disk‐like regions, e.g., Figure 10a), due to the high specific entropy of the material. In such cases, the structure and evolution of the disk‐like region cannot be considered in isolation from the rest of the structure because the influence of the vapor in the corotating and transition regions is not negligible. Due to the hot thermal state of the corotating and transition regions, the definition of the planet in the conventional analysis does not correspond to the inner edge of the disk‐like region in the post‐impact structure. For example, for the structure shown in Figures 9a–9d, in the conventional analysis the mass of the planet is 0.98 M Earth using the iterative routine of Canup et al. [2001] and 0.95 M Earth using the approach of Ćuk and Stewart [2012]. The combined mass of the planet and atmosphere using the Ćuk and Stewart [2012] approach is 0.99 M Earth . In contrast, the mass of the corotating region in the SPH post‐impact structure is only 0.88 M Earth but the radius of the corotating region extends out to ∼6500 km, similar to the radius of the planet defined in the conventional analysis, 6532 km or 6691 km. The hot vapor in the transition region extends out to beyond 8000 km, much greater than the conventional radius of the planet. Thus, the inner edge of the disk‐like region is much farther out than assumed in the conventional analysis. The thickness of the transition region may be dependent on the viscosity of the vapor [Desch and Taylor, 2013], and the thickness may not be captured accurately in SPH. However, even if the transition region was much thinner, the large radius of the hot inner portion of the structure cannot be ignored. Second, mass and AM cannot simply be lost from the disk‐like regions of a structure to the planet. The pressure at the interface between the corotating and disk‐like regions in most structures is high (103–104 bar), and these regions tend to be majority vapor. Vapor cannot simply flow across the boundary because of the pressure gradient. The fate of condensed mass in the disk that falls or spreads inward is more complicated. Falling condensate would encounter the high specific entropy of the corotating and transition regions, thermally equilibrate, and vaporize. Vaporizing condensates could increase the pressure support in the corotating region and, in turn, the disk‐like region. Increasing the vapor pressure could add vapor mass to the disk‐like region, particularly in the case of synestias. The magnitude of the effect of the infalling condensate depends on the balance between the potential energy released by the falling condensates and the redistribution of thermal energy. Crucially, mass cannot flow irretrievably from the disk‐like region to the corotating region, and previous studies have artificially depleted the mass in the disk‐like region. Third, the redistribution of the disk mass onto Keplerian orbits in the conventional analysis produces unrealistically compact surface density distributions for multiphase disks. In post‐impact structures the disk‐like regions are typically a variable mixture of condensate and vapor (for example, see the range of specific entropies in Figures 9c, 9g, and 9k). The condensates are not supported by the pressure gradient in the vapor. If the condensates decouple from the vapor, then they would fall inward via buoyancy forces. The codes that are typically used in giant impact studies do not include the separation of vapor and condensed phases. This inability of the impact code to model multiphase behavior is the rationale behind the redistribution of mass in the conventional analysis. The conventional analysis assumes that the mass fraction of vapor is small and so the effect of vapor pressure is neglected. However, for disk‐like regions that are substantially vaporized, the vapor pressure cannot be neglected. Even assuming condensate is decoupled from the vapor, the pressure gradient can support the fraction of mass that is vapor on orbits at greater semi‐major axes than purely Keplerian orbits. If the condensates are coupled to the vapor, then the whole mass of the disk can be supported by pressure gradient forces. If there is a substantial fraction of vapor, the conventional analysis produces a surface density distribution that is more compact than a pressure‐supported disk. In impacts where the majority of the disk‐like region is vaporized (e.g., the example in Figures 9e–9h), the structure of the disk‐like region is correctly modeled by the impact code, given sufficient resolution and time for the structure to reach quasi‐equilibrium (e.g., Figures 9e–9h and 10d–10f). In these cases, no post‐simulation redistribution of mass is necessary. Even so, there is a technical difficulty in calculating the post‐impact structure for a particular impact even in completely vapor structures. The numerical viscosity in the SPH code leads to some redistribution of mass and AM on the dynamical timescales required for the structure to gravitationally equilibrate. For the example in Figure 10d–10f, the mass outside of Roche increased by a few tenths of a lunar mass over days of simulation time. Thus, the equilibrium structure evolves over the calculation time. Given that the real viscosity in the structure is poorly constrained, there is uncertainty in the predicted post‐impact structure for a specific impact scenario. In a mixed phase disk, more careful analysis is needed to correctly model the dynamics and thermodynamics of liquid and vapor both during the impact and in the immediate aftermath. It is possible to quantify the error that is introduced in analyzing the output of an impact simulation using the conventional analysis in the case that the disk‐like region is majority vapor. There is initially no phase separation, and the code used to simulate the impact can accurately model the hydrostatic post‐impact structure. For cases similar to the Moon‐forming impacts proposed by Ćuk and Stewart [2012], the conventional analysis artificially decreases the mass of the disk‐like region by several lunar masses. The surface density of the disk‐like region is correspondingly reduced as shown in Figure 14. The black line in Figure 14 shows the true surface density of the structure, and the surface density of the disk in the conventional analysis is given in red. The blue dashed line shows the original position of the disk material in the conventional analysis. The mass depletion of the disk‐like region is strongest close to the planet and far from the rotational axis. Thus, the mass in the disk‐like regions can be substantially reduced by the conventional analysis. For the example in Figure 14, the mass beyond Roche is reduced from 1.7 to 0.5 lunar masses. Note that the surface density of the disk‐like regions within a synestia (black line in Figure 14) is not constant with radius as assumed in some disk evolution models. In most impact cases, some mass is injected into orbits far from the central mass. This mass will cool rapidly and can be modeled as pure condensate. Figure 14 Open in figure viewer PowerPoint Ćuk and Stewart [ 2012 The conventional analysis for processing the output of giant impact simulations artificially depletes the mass in the disk‐like regions for substantially vaporized structures. For example, the surface density in the disk‐like region of the structure shown in Figures 9 e– 9 h (black) is significantly higher than the surface density calculated using the conventional analysis in the style of] (red). The disk mass is moved inward from its original position (blue dashed). Nakajima and Stevenson [2014] have previously noted the significance of pressure gradients in the structure of highly vaporized post‐impact structures. They analyzed the disk structures generated by specific examples of canonical, high‐energy and high‐AM, and intermediate giant impacts. They found that in the disk‐like regions of high‐energy, high‐AM structures, the specific AM was nearly constant with radius, due to strong pressure support. Such structures are likely to be unstable by the Rayleigh criteria. As a result, Nakajima and Stevenson [2014] redistributed the mass in the high‐energy, high‐AM cases assuming an arbitrary stable surface density profile that conserved the mass and AM of the disk. However, they calculated the disk mass and AM to be conserved using the conventional analysis. The disk mass was significantly depleted by the use of the conventional analysis, and the calculated disks have less mass than found in the disk‐like region of synestias formed under the same impact conditions. Given the significant errors in the structure calculated using the conventional approach, it is important that better techniques are developed to analyze the results of giant impact simulations. For vapor‐dominated synestias, we suggest that removing the escaping mass from an SPH simulation and calculating the hydrostatic structure of the bound mass provides a reasonable estimate of the post‐impact structure, as shown in Figure 10. In addition, the treatment of the interface between the corotating and disk‐like regions in disk evolution models requires significant technical developments (as discussed in Charnoz and Michaut [2015]).

6.5 Synestias and Satellite Formation In the canonical disk model (e.g., Figures 9a–9d), recent work suggests that satellite accretion occurs in a multistage process [Salmon and Canup, 2012]. The material outside the Roche limit condenses and accretes quickly (on a timescale of weeks) to form a proto‐Moon. The proto‐Moon temporarily confines the edge of the Roche‐interior liquid‐vapor disk via resonant interactions. As the Roche‐interior fluid disk cools and viscously spreads both inward and outward, moonlets whose orbits are raised beyond the Roche limit are accreted onto the proto‐Moon on a timescale of hundreds of years [Machida and Abe, 2004; Salmon and Canup, 2012; Charnoz and Michaut, 2015]. The multiphase dynamics of a canonical circumterrestrial disk are challenging to model, and our understanding of the physical processes in the disk is incomplete. For example, most studies of the canonical disk evolution make an assumption of energy balance between viscous spreading of the Roche‐interior disk and radiative cooling; however, Charnoz and Michaut [2015] have shown that radiative cooling dominates the system. Furthermore, this work demonstrates that canonical disk models must include the vapor pressure support from the corotating region. Satellite accretion from the various post‐impact structures described here can be fundamentally different than accretion from a canonical circumterrestrial disk. Post‐impact synestias can have a much higher surface density in the disk‐like region (e.g., Figure 9). Consequently, the cooling time of vapor beyond the Roche limit reaches timescales relevant for satellite formation. For example, Lock et al. [2016] estimate that substantial vapor pressure (greater than bars) can persist at the Roche limit for tens of years after the impact for the synestia example in Figures 9e–9h. Synestias do not have a shear boundary between the corotating region and disk‐like region, and therefore, the whole structure can mix more easily than in the sub‐CoRoL structures formed in canonical giant impacts. Condensates that decouple from the vapor in the lowest surface density regions beyond the Roche limit have angular momenta such that they orbit within the vapor structure. In some cases, the proto‐Moon may form beyond the Roche limit and grow to near its final mass while remaining within the Earth‐composition vapor of the synestia. Lock et al. [2016] proposed that accretion within an impact‐generated synestia can explain the isotopic and chemical properties of our Moon. A future paper will investigate the evolution of synestias and the formation of our Moon.

6.6 Planet Formation With Hot Rocky Bodies We have shown that partially vaporized post‐impact bodies are common during accretion (section 5); however, the implications of the highly variable thermal and physical structures of such bodies on the formation of terrestrial planets have not been considered. Here we comment on a few key areas of planet formation affected by the possible range of planetary structures: accretion efficiency, core formation, and chemical evolution. At present, models of the accretion of terrestrial planets do not account for changes in thermal state. The physical size of planetesimals and planets is usually calculated by assuming a fixed bulk density over the entire calculation [typically 3000 kg m−3 in N‐body codes, e.g., Quintana et al., 2016]. For rotating rocky bodies with the highest specific entropy cases considered here, the bulk density could be an order of magnitude lower. Figure 15 presents the bulk density and areal cross sections for the suite of Earth‐mass bodies shown in Figure 4. Partial vaporization is accompanied by a significant decrease in bulk density for all rotational states. Post‐impact states can be significantly larger than isolated bodies with even lower density and larger cross sections. All giant impacts results in some inflation of the radius of a body. The larger sizes of hot rocky planets (Figure 4) translate to larger collisional cross sections and an increased probability of impacts, and hence more efficient accretion (Figures 15b and 15c). In particular, the expanded cross section of a body after an impact will increase the fraction of impact debris reaccreted in the early time period before the debris can be perturbed from crossing orbits. In proposed Moon‐forming giant impacts, typically a few to several lunar masses of material is ejected [Canup, 2004; Ćuk and Stewart, 2012; Canup, 2012], which is accreted onto the Earth and nearby planets within ∼10 Myr under current assumptions [Jackson and Wyatt, 2012; Bottke et al., 2015]. It has been suggested that the warm debris disks by giant impacts would be observable around other stars [Jackson and Wyatt, 2012]. But recently, Kenyon et al. [2016] noted that the observed occurrence rate of warm debris disks around young stars is much lower than what would be expected from the prevalence of extrasolar rocky planets around mature stars and current planet formation models. Kenyon et al. [2016] inferred that terrestrial planet formation must be quick and neat, with efficient reaccretion of impact‐produced debris. Realistic capture cross sections of planets could be a factor in explaining such observations and should be included in planet formation studies. Figure 15 Open in figure viewer PowerPoint . Colored lines show bodies of constant angular momentum. The extended structures of hot, rotating rocky bodies lead to very low bulk den