This paper is organized as follows. First, we describe the structures that are generated by giant impacts (sections 2.1 and 2.2 ). Next, we discuss the processes that dominate the evolution of a synestia as it cools (section 2.3 ). We present calculations of the pressure field of a synestia and estimate the mass and location of the moon that is formed (sections 2.4 and 2.5 ). Then, we present a calculation of the phase diagram for BSE composition material over the pressure and temperature range of the outer portions of the structure (section 3 ). In section 4 , we combine the results of the previous sections to propose a coupled dynamic and thermodynamic model for the formation of a moon from a terrestrial synestia and identify the pressure‐temperature‐spatial paths of condensates and growing moonlets. The chemical composition of the moon is estimated using the physical chemistry of the BSE system at the pressure and temperature predicted by our model. In section 5 , we discuss formation of our Moon from a synestia as a new model for lunar origin. We discuss the consistency of such a model with observations and examine the possible range of giant impacts that may generate post‐impact structures with the potential for forming our Moon. We present a short unified synopsis of our model for the origin of the Moon in section 6 . Finally, we draw our major conclusions and describe future tests of our model (section 7 ). This work includes supporting information .

At present, no single calculation can fully capture the dynamics, thermodynamics, and chemistry of lunar accretion. Therefore, our approach is to link the physics and chemistry of satellite accretion from a terrestrial synestia by understanding the processes that control the pressure and temperature paths of the material that forms a moon. First, we determine the pressure‐temperature conditions of a moon that grows by accretion of condensing silicate vapor. We then argue that the composition of the growing moon is set by equilibrium with BSE vapor over a specific range of pressures and temperatures determined by the structure and the phase relationships for material of BSE composition. Finally, we demonstrate that a variety of high‐energy, high‐AM giant impacts can generate initial conditions that can potentially lead to the formation of a lunar‐mass moon with the observed geochemical characteristics of our Moon. Our model provides a promising pathway to explain all the key observables of the Moon discussed above: the isotopic similarity between Earth and the Moon; the magnitude and pattern of moderately volatile element depletion in the Moon; and the large mass of the Moon. If Earth had a large obliquity after the giant impact, then a single event may also explain the inclination of the lunar orbit and the present‐day AM of the Earth‐Moon system (Ćuk et al., 2016 ).

Here we present a new model for lunar origin within a terrestrial synestia, an impact‐generated structure with Earth mass and composition that exceeds the corotation limit (CoRoL). Synestias are formed by a range of high‐energy, high‐AM collisions during the giant impact stage of planet formation (Lock & Stewart, 2017 , hereafter LS17). A synestia is a distinct dynamical structure compared to a planet with a condensate‐dominated circumplanetary disk, and, as a result, different processes dominate the early evolution of a synestia. Note that preliminary versions of this work (e.g., Lock et al., 2016 ; Petaev et al., 2016 ) used different nomenclature than is used here. In particular, synestias were referred to as continuous mantle‐atmosphere‐disk (MAD) structures.

A substantial constraint on the canonical Moon‐forming impact is that the AM of the Earth‐Moon system has not changed significantly since the formation of the Moon. Ćuk and Stewart ( 2012 ) showed that an evection resonance could drive significant AM loss from the Earth‐Moon system after the Moon‐forming impact. Since Ćuk and Stewart ( 2012 ), additional mechanisms have been found that could remove AM during lunar tidal evolution (Ćuk et al., 2016 ; Tian et al., 2017 ; Wisdom & Tian, 2015 ). Allowing for a change of AM after the impact significantly expands the range of possible impact parameters for the Moon‐forming collision. Ćuk and Stewart ( 2012 ) and Canup ( 2012 ) showed that high‐energy, high‐AM impact events can inject much more material into orbit than canonical impacts. In addition, Ćuk et al. ( 2016 ) showed that if the Earth after the impact had both high AM and high obliquity, an instability during the Laplace plane transition could both remove AM from the Earth‐Moon system and explain the Moon's present‐day orbital inclination. With these promising dynamical results, high‐AM giant impact scenarios for lunar origin warrant continued investigation.

Despite the fact that it has not yet explained major characteristics of the Earth‐Moon system, the giant impact hypothesis has not been rejected, primarily due to the lack of another viable mechanism for the origin of the Moon. A range of alternative impact models have been proposed (Canup, 2012 ; Ćuk & Stewart, 2012 ; Reufer et al., 2012 ; Rufu et al., 2017 ), but each calls upon an additional process or a fortunate coincidence to better explain the Earth‐Moon system. Hence, none of these recent variations on an impact origin have gained broad support.

The origin of the Moon's present‐day orbital inclination, which is about 5° from the ecliptic plane, has been a long standing problem in lunar tidal evolution. If a Moon‐forming giant impact also determined Earth's present obliquity, then lunar origin in an equatorial disk and subsequent tidal evolution through the Laplace plane transition, from an orbit that precesses in Earth's equatorial plane to one that precesses in the ecliptic plane, should have led to a lunar orbit with near‐zero inclination to the ecliptic. Therefore, dynamical processes subsequent to the impact are required to explain the present lunar inclination. Proposed solutions in the framework of the canonical model include a complex sequence of luni‐solar resonances (Touma & Wisdom, 1994 ), resonant interactions between the Moon and the circumterrestrial disk (Ward & Canup, 2000 ), and encounters between large planetesimals and the newly formed Earth‐Moon system (Pahlevan & Morbidelli, 2015 ). Recently, Chen and Nimmo ( 2016 ) and Ćuk et al. ( 2016 ) investigated inclination damping by lunar obliquity tides, and Ćuk et al. ( 2016 ) found that lunar inclination must have been large (∼30°) prior to the point in tidal recession where the lunar orbit transitions between Cassini states (distinct dynamical solutions that govern the alignment of the lunar spin axis and orbital plane, see Peale, 1969 ). Such a large inclination prior to the Cassini state transition defies explanation by any of the previously proposed mechanisms to raise lunar inclination after a canonical giant impact. Connecting the canonical giant impact to the Moon's current orbit remains an unsolved problem.

Predicting the final mass of satellites formed by giant impacts is challenging. The methods that are currently used for simulating giant impacts do not include the physics necessary for modeling lunar accretion; therefore, separate calculations of disk evolution are required to infer the mass of the satellite produced by a specific impact. Typically, scaling laws fitted to N ‐body simulations of idealized circumterrestrial debris disks (Ida et al., 1997 ; Kokubo et al., 2000 ) have been used to estimate the satellite mass from the total mass and AM of orbiting material (Canup, 2004 , 2008a , 2012 ; Canup & Asphaug, 2001 ; Ćuk & Stewart, 2012 ). Studies of canonical impacts have found that, over a narrow range of impact angles, sufficient mass is injected into orbit to produce a lunar‐mass Moon based on N ‐body scaling laws. Because N ‐body simulations do not include the multiphase physics of the lunar disk, they overestimate the efficiency of satellite formation. Simulations that include a simplified one‐dimensional model of Roche‐interior multiphase disks have inferred much lower accretion efficiencies (Salmon & Canup, 2012 , 2014 ). The Roche limit is the closest distance a satellite can withstand the tidal forces from the planet (about 18,500 km for silicate satellites orbiting Earth). Using the scaling laws produced by these most recent models, very few of the disks produced in published canonical giant impact simulations inject the required mass and AM into orbit to produce a lunar‐mass satellite ( supporting information section S1). Furthermore, the simple Roche‐interior disk model used by Salmon and Canup ( 2012 , 2014 ) likely overestimated the efficiency of the spreading of material beyond the Roche limit. By incorporating more multiphase physics, Charnoz and Michaut ( 2015 ) showed that viscous spreading of material beyond the Roche limit is slower than that calculated by Salmon and Canup ( 2012 , 2014 ) and that more mass from the disk is lost to Earth. Canonical impacts typically inject a large amount of mass directly beyond the Roche limit, and Charnoz and Michaut ( 2015 ) suggested that the Moon could have largely formed from this material. The efficiency of accretion has not been quantified and such a model would still need to explain the isotopic similarity and moderately volatile element depletion. Given the current results from giant impact calculations and available satellite accretion scaling laws, it is uncertain whether canonical giant impacts can form a sufficiently large moon.

Study of lunar samples has revealed that the Moon is significantly depleted in moderately volatile elements (MVEs; e.g., K, Na, Cu, and Zn) relative to the bulk silicate Earth (BSE). For example, potassium and sodium are inferred to be depleted by factors of 5 to 10 compared to terrestrial abundances (e.g., Ringwood & Kesson, 1977 , see section 5.1 ). In a series of studies (e.g., Ringwood, 1986 ; Ringwood & Kesson, 1977 ), Ringwood and colleagues argued that the lunar composition could be explained if the Moon was a partial condensate of vapor derived from Earth's mantle. The MVE depletion of the Moon is a key constraint on lunar origin models. In addition, the volatile depletion of the Moon has been used to argue for a process‐based link between giant impacts and MVE loss. Indeed, the lunar depletion has been used to propose that Mercury would be depleted, if it was formed by a giant impact (Peplowski et al., 2011 ). We must understand the physical processes that led to volatile depletion on the Moon in order to place its data in the context of other bodies in the solar system. Few studies have attempted to combine the dynamics, thermodynamics, and chemistry of lunar origin, which is necessary to be able to test the proposed models. Recently, Canup et al. ( 2015 ) used the lunar disk models of Salmon and Canup ( 2012 ) and physical chemistry calculations to link the dynamics and thermodynamics of accretion from a canonical disk. They also suggested that the lunar volatile element depletion could be explained if the material that formed the observable Moon was a partial condensate of disk material. Wang and Jacobsen ( 2016 ) recently reported that the potassium isotopes of the Moon are heavier than BSE, which supports the idea of partial condensation. However, the model presented by Canup et al. ( 2015 ) does not quantitatively explain the magnitude, nor pattern, of moderately volatile element depletion observed for the Moon. Small isotopic fractions could be produced if the canonical disk had a period of hydrodynamic loss (Pritchard & Stevenson, 2000 ), but, given the mean molecular weight of vapor in the disk, hydrodynamic escape is unlikely (Nakajima & Stevenson, 2018 ). Further work is needed to fully integrate chemical and physical models of lunar origin.

Two classes of solutions to the problem of isotopic similarity have been proposed. First, post‐impact mixing between the planet and lunar disk could have erased initial isotopic heterogeneities (Pahlevan & Stevenson, 2007 ), but the extent of mixing required to explain the observations is a problem in the canonical model (Melosh, 2014 ). Second, the impactor and proto‐Earth could have formed from the same source material and thus shared nearly identical isotopic signatures (Dauphas, 2017 ; Dauphas et al., 2014 ; Jacobsen et al., 2013 ). There is evidence for a reservoir of terrestrial precursor materials with fractionation‐corrected isotopic ratios that are distinct from the meteorites and planetary samples in our collections (e.g., Drake & Righter, 2002 ). If the impactor and target accreted the majority of their mass from the same reservoir, the Earth and Moon would share similar stable isotopic ratios. Stable isotopic ratios are controlled by the source material only and not affected by processes within the body (e.g., O, Cr, Ti). Yet even if bodies in the inner solar system were formed from material with similar isotopic signatures, this explanation for the isotopic similarity between Earth and the Moon relies upon a coincidence to explain tungsten, which is sensitive to the conditions and timing of core formation on each of the colliding bodies (Dauphas, 2017 ; Dauphas et al., 2014 ; Kruijer & Kleine, 2017 ; Kruijer et al., 2015 ; Touboul et al., 2015 ). As discussed in Melosh ( 2014 ) and Kruijer and Kleine ( 2017 ), neither of these proposals on their own provides a satisfactory explanation for the isotopic similarity between Earth and the Moon. Nevertheless, the possibility of a more homogeneous inner solar system relaxes the isotopic constraint on Moon‐formation, as models need only produce similar enough tungsten isotopes in Earth and the Moon. The tungsten isotope constraint is weaker than that for stable isotopes due to the uncertainties in inferring the post‐impact composition before the addition of the late veneer (Kruijer et al., 2015 ; Touboul et al., 2015 ).

Numerical simulations of giant impacts predict that the canonical lunar disk is derived primarily from the impactor (Canup, 2004 , 2008a ; Canup & Asphaug, 2001 ). However, increasingly precise isotopic measurements of terrestrial and lunar samples have shown that Earth and the Moon share very similar initial isotopic ratios for a wide range of elements (Lugmair & Shukolyukov, 1998 ; Wiechert et al., 2001 ; Zhang et al., 2012 ). Because the isotope ratios of such elements are observed to vary significantly among planetary bodies (Clayton & Mayeda, 1996 ; Trinquier et al., 2007 ; Yin et al., 2002 ; Zhang et al., 2012 ), the impactor is generally expected to have had a distinct isotopic composition, resulting in a measurable isotopic difference between Earth and the Moon (Pahlevan & Stevenson, 2007 ; Melosh, 2014 ; Young et al., 2016 ).

Most studies of the origin of the Moon have focused on a narrow range of impact scenarios. Cameron and Ward ( 1976 ) proposed that the Moon‐forming giant impact could have prescribed the present‐day angular momentum (AM) of the Earth‐Moon system. Numerical simulations have shown that a grazing collision with a Mars‐mass impactor near the mutual escape velocity can impart the present‐day AM and generate a silicate‐rich disk composed of more than a lunar mass of material (Canup, 2004 , 2008a ; Canup & Asphaug, 2001 ). This scenario, which we refer to as the canonical giant impact, has become the de facto working model for lunar origin. However, studies of the canonical impact and its aftermath have difficulty explaining some key observables of the Earth‐Moon system, including: the isotopic similarity between Earth and the Moon; the lunar depletion in moderately volatile elements; the large mass of the Moon; and the present‐day lunar inclination.

In the giant impact hypothesis for lunar origin (Cameron & Ward, 1976 ; Hartmann & Davis, 1975 ), the proto‐Earth suffered a collision with another protoplanet near the end of accretion that ejected material into a circumterrestrial disk, out of which the Moon formed (see reviews by Asphaug, 2014 ; Barr, 2016 ; Stevenson, 1987 ). Giant impacts are highly energetic events that vaporize a portion of the impacting bodies. Hence, the disk is a multiphase mixture of liquid and vapor (Canup, 2004 , 2008a , 2012 ; Canup & Asphaug, 2001 ; Ćuk & Stewart, 2012 ). Modeling the formation of the Moon from such a disk is challenging and the details of lunar accretion are still uncertain (e.g., Carballido et al., 2016 ; Charnoz & Michaut, 2015 ; Gammie et al., 2016 ; Machida & Abe, 2004 ; Salmon & Canup, 2012 ; Thompson & Stevenson, 1988 ; Ward, 2012 , 2014 , 2017 ). To date, lunar origin studies have not demonstrated that a single giant impact can explain both the physical and chemical properties of our Moon (Asphaug, 2014 ; Barr, 2016 ).

2 Structure and Dynamics of a Synestia

In this section, we examine the physical processes that occur during the formation and evolution of synestias. First, we describe the physics controlling the structure of post‐impact states and demonstrate the magnitude of pressure support in synestias (section 2.1). Next, we look at the transition from the impact to a post‐impact state (section 2.2). In section 2.3 we discuss the dominant processes that drive evolution of the synestia and argue that condensates formed at the photosphere can transport mass radially in the structure. Based on these arguments, we construct a simple cooling model to calculate the temporal evolution of the vapor structure of a synestia (section 2.4). In section 2.5 we present an example calculation of the cooling of a potential Moon‐forming synestia. Additional examples of cooling of potential Moon‐forming synestias, formed by very different impact events, are provided in the appendix.

2.1 Structure of Post‐impact States (1) G is the gravitational constant, M bnd is the bound mass of the structure, r xy is the cylindrical radius, ρ is the density of the parcel, p is the gas pressure, and ω is the angular velocity. The vertical structure, in the direction parallel to the rotation axis, is also controlled by hydrostatic balance. Figures 1 a and 1 e show the approximate fluid pressure structure of two post‐impact states, calculated using a smoothed particle hydrodynamics (SPH) code (see section 2.4 for details of methods), which evolves the position and thermal properties of particles of fixed mass using the governing forces and material equations of state (EOS). The two examples illustrate a post‐impact structure that exceeds the CoRoL (LS17) and one that does not. Both impacts generate significant amounts of vapor. The fluid structures are controlled by a balance between the gravitational, pressure gradient, and centrifugal forces. For a parcel of material in the midplane, the force balance is approximatelywhereis the gravitational constant,is the bound mass of the structure,is the cylindrical radius,is the density of the parcel,is the gas pressure, andis the angular velocity. The vertical structure, in the direction parallel to the rotation axis, is also controlled by hydrostatic balance. Figure 1 Open in figure viewer PowerPoint S crit =5.40 kJ K−1 kg−1, p crit =25.5 kbar, T crit =8,810 K, ρ crit =1,680 kg m−3). Material to the left of the dome is liquid, material to the right of the dome and below the critical point is vapor, material above and to the right of the critical point is supercritical fluid (SCF), and material within the dome is a mixture of both liquid and vapor (blue points). The liquid‐solid phase boundary is neglected in this equations of state. The calculated pressure and specific entropy distributions (a–c and e–g) estimate the overall shape of gravitationally and thermally equilibrated post‐impact structures, but the pressure structure at the photic surface was not resolved. During thermal equilibration, the Roche‐exterior condensate fraction is removed and the post‐impact distribution of specific entropy (d, h) is averaged such that the middle of the structure has a mass‐weighted isentropic region (vertical set of points in (c) and (g)) which transitions to a saturated vapor region that follows the vapor side of the dome. In (g), the isentropic and saturated vapor regions of the thermally equilibrated structures are labeled. The inner edge of the saturated vapor region is identified by the black line in (a), (b), (e), and (f). The structure in (a)–(d) was generated by an impact of a 0.13 M Earth projectile onto a 0.9 M Earth target with an impact velocity of 9.2 km s−1 and an impact parameter of 0.75. For (e)–(h), a 0.47 M Earth projectile hit a 0.57 M Earth target with an impact velocity of 9.7 km s−1 and an impact parameter of 0.55. Midplane profiles are presented in Figure Examples of sub‐corotation limit (sub‐CoRoL) and super‐corotation‐limit (super‐CoRoL) post‐impact structures generated by canonical (a–d) and high‐energy, high‐AM (e–h) giant impacts, showing axisymmetric pressure (a, e) and silicate specific entropy (b, f) contours perpendicular to the spin axis, and the thermal state of the silicate material in the midplane in pressure‐specific entropy space (c, d, g, and h). The liquid‐vapor phase boundary is a dome‐shaped curve in pressure‐specific entropy space (black line). The black dot on the vapor dome (c, d, g, and h) is the critical point for the equation of state used in these simulations (=5.40 kJ Kkg=25.5 kbar,=8,810 K,=1,680 kg m). Material to the left of the dome is liquid, material to the right of the dome and below the critical point is vapor, material above and to the right of the critical point is supercritical fluid (SCF), and material within the dome is a mixture of both liquid and vapor (blue points). The liquid‐solid phase boundary is neglected in this equations of state. The calculated pressure and specific entropy distributions (a–c and e–g) estimate the overall shape of gravitationally and thermally equilibrated post‐impact structures, but the pressure structure at the photic surface was not resolved. During thermal equilibration, the Roche‐exterior condensate fraction is removed and the post‐impact distribution of specific entropy (d, h) is averaged such that the middle of the structure has a mass‐weighted isentropic region (vertical set of points in (c) and (g)) which transitions to a saturated vapor region that follows the vapor side of the dome. In (g), the isentropic and saturated vapor regions of the thermally equilibrated structures are labeled. The inner edge of the saturated vapor region is identified by the black line in (a), (b), (e), and (f). The structure in (a)–(d) was generated by an impact of a 0.13projectile onto a 0.9target with an impact velocity of 9.2 km sand an impact parameter of 0.75. For (e)–(h), a 0.47projectile hit a 0.57target with an impact velocity of 9.7 km sand an impact parameter of 0.55. Midplane profiles are presented in Figure 3 In debris disks, the effect of the pressure support term on the mass distribution is negligible as the gas fraction is small. The density of condensates is much larger than that of vapor, and consequently the effect of the vapor pressure term in equation 1 on condensed particles is small. In prior work on post‐impact states (e.g., Canup, 2004, 2008a; Canup & Asphaug, 2001; Ida et al., 1997; Kokubo et al., 2000; Nakajima & Stevenson, 2014, 2015; Salmon & Canup, 2012, 2014; Ward, 2012, 2014, 2017), the disk was assumed to be dominated by condensates with negligible pressure support. Figure 2 Open in figure viewer PowerPoint δ xy . Sections of profiles are not shown due to model artifacts associated with the averaging of the isentropic region of the structures. The profiles for the sub‐corotation‐limit (sub‐CoRoL) example are more variable due to fewer particles and hence lower resolution in the disk‐like region. The equivalent profiles for different example post‐impact structures are shown in Figures A3. The disk‐like regions of post‐impact states can be strongly pressure supported. We calculate the magnitude of the terms in equation 1 for the two example, post‐impact structures shown in Figures 1 and 3 . (a, d) Profiles of the specific gravitational force (solid green), specific pressure gradient force (yellow), and specific centrifugal force (navy) in the midplane. Quantities have been calculated by taking the average properties of particles in the midplane in 1 Mm bins and calculating the gradients of pressure and potential using these average properties. Dashed green lines show the gravity assuming that the body is spherically symmetric, neglecting higher‐order terms. (b, e) The ratio of the pressure gradient and gravitational forces. (c, f) The difference between the orbit of the pressure‐supported mass and a circular Keplerian orbit of the same angular momentum in the disk‐like regions,. Sections of profiles are not shown due to model artifacts associated with the averaging of the isentropic region of the structures. The profiles for the sub‐corotation‐limit (sub‐CoRoL) example are more variable due to fewer particles and hence lower resolution in the disk‐like region. The equivalent profiles for different example post‐impact structures are shown in Figures A1 The thermal structure of our two example post‐impact states are shown in Figures 1b and 1f. There is a strong gradient in entropy from the inner to the outer regions of the structure. Figures 1c and 1d present the thermal structure in the specific entropy‐pressure phase space with each SPH particle shown as a colored dot. The black curve is the liquid‐vapor phase boundary for the single‐component silicate EOS used in the simulation. Material below the black line is a mixture of liquid and vapor, and material above the curve is pure vapor, supercritical fluid (SCF), or liquid. Post‐impact structures are highly thermally stratified. Typically, the thermal profile is such that it does not intersect the liquid‐vapor phase boundary until low pressures and the silicate transitions smoothly from liquid to supercritical liquid to vapor (LS17). Thus, the post‐impact structures have no surface, and a liquid‐vapor mixture is initially restricted to the outer regions and near the photosphere (see section 2.3). In the examples shown in this work, the midplane is initially completely vapor to beyond the Roche limit (within the black lines in Figures 1a, 1b, 1e, and 1f). At low pressures, the structure intersects the liquid‐vapor phase boundary and follows a saturated adiabat. Beyond the black line in Figures 1a, 1b, 1e, and 1f, the saturated adiabat extends to the midplane and a fraction of the silicate is condensed. In most of the structure, the condensed mass fraction is small, and we have neglected the gravitational effects of the condensed mass in calculating the vapor pressure structure. Due to the dominance of vapor in post‐impact structures, there is substantial pressure support, and the pressure gradient term in equation 1 is not negligible. (2) Figure 2 shows the relative magnitudes of the gravitational, pressure gradient, and centrifugal terms in equation 1 for the example post‐impact structures in Figure 1 . In the corotating region, the pressure gradient term (equation 1 ) is comparable to gravity as expected, but the pressure gradient force can also be the same order of magnitude as gravity in the outer regions. Additional examples of the magnitude of pressure support in synestias are given in Figures A1 A3 . Equation 1 can be rearranged to give an expression of the angular velocity of the gas at a given cylindrical radius, A negative pressure gradient term reduces the angular velocity at a given radius. In other words, a parcel of material with a fixed specific AM may orbit at a larger radius with the addition of a pressure gradient force. If the pressure support was subsequently removed, this parcel of material would evolve to a Keplerian orbit closer to the rotation axis in order to conserve AM. Figures 2c and 2f show the difference between the cylindrical radius of vapor in post‐impact states and circular Keplerian orbits of the same AM, δr xy . The presence of the pressure gradient in post‐impact structures supports material many megameters farther away from the rotation axis than would be possible based solely on the specific AM of the structure and a balance between gravity and centrifugal forces. All post‐impact structures have a corotating inner region and a disk‐like outer region (Figure 3). The equatorial radius of the impact‐heated corotating region is larger than the radius of a non‐rotating, condensed planet of the same mass (LS17). Different post‐impact structures have corotating regions that rotate at different rates. In some cases, there is a difference between the corotating angular velocity and the angular velocity at the inner edge of the disk‐like region, which requires the presence of a transition region between the corotating and disk‐like regions (e.g., the portion of the gray line that lies above the dashed black line in Figure 3b). The transition region in the fluid has a monotonically increasing angular velocity and, in the midplane, is similar to the profile for shear between two cylinders (e.g., Chandrasekhar, 1969; Desch & Taylor, 2013). Structures with a significant transition region (such as the example in Figures 1a–1c and 3a–3c) tend to be below the CoRoL. The CoRoL is defined by where the angular velocity at the equator of a corotating planet intersects the Keplerian orbital velocity (LS17). The CoRoL is a surface that depends on thermal state, AM, total mass, and compositional layering. In contrast, the structure shown in Figures 1e–1g and 3d–3f is above the CoRoL. The corotating region rotates much more rapidly than in the sub‐CoRoL example. There is no transition region between the corotating and disk‐like regions, and the angular velocity profile decreases monotonically with radius. Bodies above the CoRoL are called synestias (LS17). Figure 3 Open in figure viewer PowerPoint The sub‐corotation‐limit (sub‐CoRoL) and super‐corotation‐limit (super‐CoRoL) post‐impact structures produced by canonical (a‐c) and high‐energy, high‐angular‐momentum (high‐AM) (d‐f) giant impacts have important differences. Midplane profiles of angular velocity (a, d), specific angular momentum (b, e), and surface density (c, f) correspond to the thermally and gravitationally equilibrated post‐impact structures shown in Figure 1 . Blue points indicate smoothed particle hydrodynamics particles that are within the vapor dome and are of mixed phase. Solid black lines denote the angular velocity or specific angular momentum (AM) of a circular Keplerian orbit around a point mass, and dashed black lines indicate the specific AM for corotation with the inner region of the structure. Gray lines show the mass‐weighted average of the angular velocity and specific AM of particles in the midplane. Both surface density and average properties are calculated in 2 Mm radial bins. Sub‐CoRoL structures have a substantial shear boundary between the corotating region and disk‐like region (points above the dashed line in b). Super‐CoRoL structures, synestias, have a monotonically decaying angular velocity profile between the corotating and disk‐like regions (d). Black arrows indicate the approximate outer edge of the corotating regions (8.5 Mm and 8.6 Mm, respectively). Red arrows indicate the inner edge of the isentropic regions (7.9 Mm and 9.9 Mm). Blue arrows indicate the outer edge of the isentropic regions (24.3 Mm and 23.6 Mm). The slight kinks in the profiles in (a), (b), (d), and (e) are caused by the entropy discontinuity at the inner edge of the isentropic region imposed during processing of post‐impact structures (section 2.4 ). The substantial pressure support in the disk‐like region of post‐impact structures leads to sub‐Keplerian angular velocities (points below the black line in Figures 3a, 3b, 3d, and 3e). The effect of pressure support extends farther away from the rotation axis in the example synestia. The surface density in the disk‐like region of the synestia is an order of magnitude higher than in the sub‐CoRoL case (Figure 3), and consequently there are much higher pressures in the midplane (Figure 1). The surface density in the disk‐like region of the sub‐CoRoL structure is approximately constant (Figure 1c) as reported in previous studies (Canup, 2004, 2008a, 2008b; Canup & Asphaug, 2001), but the constant surface density region begins at a larger radius than previously reported. The dynamic and thermodynamic structure of post‐impact states depends strongly on the parameters of the impact that produced them. The synestia example shown here is typical of the impact‐generated synestias found by LS17. All studies of giant impacts in the regime of the Moon‐forming event report substantial production of vapor. However, in previous work on post‐impact structures (e.g., Canup, 2004, 2008a; Canup & Asphaug, 2001; Ida et al., 1997; Kokubo et al., 2000; Nakajima & Stevenson, 2014, 2015; Salmon & Canup, 2012, 2014; Ward, 2012, 2014, 2017), it was assumed that the vapor in the structure would cool rapidly and the pressure gradient term in equation 1 could be neglected. Prior post‐impact analyses generally simplified the structure to a point source planet and near‐Keplerian disk. In general, the immediate postimpact structure cannot be analyzed in this manner (LS17). Our conclusion that pressure has a significant effect in post‐impact structure is not in conflict with the findings of prior work. For example, Nakajima and Stevenson (2014) note the importance of the pressure gradient force. We analyzed the post‐impact structures calculated by Nakajima and Stevenson (2014) (pressure contours received from M. Nakajima by personal communication) to derive the relative magnitude of the force terms in equation 1 and found that there are also strong radial pressure gradient forces. However, in calculating the surface density of their disk structures, Nakajima and Stevenson (2014) generally neglected the pressure support. Based on the results presented in this section, we find that post‐giant‐impact structures must be analyzed as continuous rotating fluids and cannot be separated into a planet and near‐Keplerian disk. The properties of the structure depend on the integrated effects of radial force balance (equation 1), vertical hydrostatic equilibrium, and the thermal, mass, and angular momentum distributions imparted by the impact event, over the entire structure.

2.2 Transition From the Impact to the Post‐impact Structure The methods that are currently used to simulate giant impacts do not include multiphase flow processes or thermal equilibration between parcels of material. In this section, we discuss the processes that govern the transition from the impact event to the post‐impact state, the starting point for subsequent evolution and satellite formation. During and after the impact, condensates separate from the vapor and fall radially inward and toward the midplane (section 2.3). The separation of the condensate from the vapor in the outer regions of the post‐impact structure changes the mass distribution and thermal structure. Condensates that fall into higher‐pressure, hotter regions of the structure, vaporize and transfer their mass to vapor (section 2.3). The fluid structure (Figure 1) adapts to the redistribution of mass and entropy. Condensates that have sufficient AM to remain beyond the Roche limit can accrete to form satellites. These processes occur on dynamical timescales of hours to days. Simultaneously, the vapor in the structure will convect, leading to rapid mixing vertically (parallel to the rotation axis). In all the example post‐impact structures considered in this work, the specific entropy of the structure within the Roche limit is such that the equilibrium phase in the midplane is pure vapor (inside black lines in Figures 1a, 1b, 1e, and 1f). The vertical thermal profile in these regions is adiabatic until it intersects the liquid‐vapor phase boundary at lower pressures, at which point it follows a saturated adiabat. Farther out in the structure (beyond the black lines in Figures 1a, 1b, 1e, and 1f), the midplane pressure is lower, and the whole column is on a saturated adiabat. As noted by Nakajima and Stevenson (2014), a portion of each post‐impact structure is approximately isentropic. Generally, the quasi‐isentropic region has a specific entropy that exceeds the critical point (LS17). This region corresponds to the approximately vertical subset of red SPH particles at pressures below the critical point in Figures 1d and 1h, which shows the thermal state of particles at the end of an SPH impact simulation. The methods that are currently used to simulate giant impacts do not model the condensate separation or fluid convection that occurs in the hours after the impact. It is necessary to process the output from impact simulations to mimic the effect of these processes on the post‐impact structure. We took the thermal profile from SPH impact simulations (Figures 1d and 1h), redistributed and removed condensate, and averaged the quasi‐isentropic region to a single isentrope (Figures 1c and 1g). The details of this processing step are described in section 2.4 and section S4. Following the rapid adjustment after the impact, there is a longer period of secular cooling of the post‐impact structure that we discuss in the following sections.

2.3 Cooling a Synestia In sections 2.1 and 2.2, we considered the structures of both sub‐CoRoL structures and synestias. Henceforth, we focus solely on synestias and describe their evolution in the years following the impact. In particular, we consider the evolution of terrestrial synestias, those with Earth‐like mass and composition. Some of the processes we consider are universal to all post‐impact structures, but the evolution of sub‐CoRoL structures is left to future work. A terrestrial synestia, as shown in Figures 1e–1g, cools by radiation from the photosphere, where the structure is optically thin. We estimated the optical depth of the outer edge of the structure and found that the photic surface is at low pressures (10−6 to 10−2 bar, supporting information section S2) and radiates at a temperature of about T rad =2,300 K, determined by the liquid‐vapor phase boundary of BSE composition material (section 3). At the photosphere, the majority of the energy lost by radiation is compensated for by condensation of vapor, and the material state is a mixture of vapor and condensates (liquid droplets and/or solid dust). Since condensates are not supported by the pressure gradient of the vapor structure, they are not dynamically stable at the photic surface and will fall. Because the photospheric temperature is far above equilibrium with the incoming solar radiation, radiative cooling is efficient and drives a torrential rain of condensates into the higher‐pressure regions of the structure. Initially, radiative cooling produces on the order of a lunar mass of silicate condensate per year, and the production rate of silicate rain near the photosphere of a synestia (a few centimeters an hour) is about an order of magnitude greater than heavy rainfall during hurricanes on Earth today. (3) F D is the drag force from the gas, m cond is the mass of the condensate, and g is the gravitational acceleration. The droplet size is then given by the balance between the drag force and surface tension over the droplet, (4) ρ cond is the density of the condensate, R cond is the radius of the droplet, and γ is the surface tension. Boča et al. ( 2003 −1. From the balance of the forces, the radius of the droplet is (5) The size of droplets forming at the photosphere and falling is controlled by a balance between shear from the vapor and surface tension. We approximate the drag force on condensates by assuming that the droplets are in free fall such thatwhereis the drag force from the gas,is the mass of the condensate, andis the gravitational acceleration. The droplet size is then given by the balance between the drag force and surface tension over the droplet,whereis the density of the condensate,is the radius of the droplet, andis the surface tension. Boča et al. () found that the surface tension for silica melts at 2,000 K is 0.3 N m. From the balance of the forces, the radius of the droplet is The gravitational acceleration at the photosphere ranges from 10 to 0.1 m s−2 due to the large spatial scale of the synestia. The corresponding range in droplet sizes is a few millimeters to a few centimeters. As discussed in section 2.1, there is a substantial pressure gradient force perpendicular to the rotation axis acting on the vapor in synestias. When condensates form from the vapor at the photosphere, they do not have sufficient AM to remain in a circular orbit at the same location. In the absence of vapor, the condensates would fall rapidly on significantly elliptical orbits in a plane through the center of mass. For the synestia shown in Figures 1-3, the initial eccentricity varies from 0.6 at the inner edge of the disk‐like region to 0.27 at a cylindrical radius of 25 Mm. As they fall into the synestia, condensates interact with the vapor and experience drag which perturbs their motion. (6) ρ vap is the density of the vapor, C D is the gas drag coefficient, v is the velocity vector of the condensate, and v vap is the velocity vector of the vapor. A set of simple classical dynamics equations were then integrated to find the position and velocity of the particle as a function of time. Full details of the calculation are given in Appendix To calculate the motions of condensates in the synestia, we used a simple orbital evolution model, including an acceleration due to gas drag of the formwhereis the density of the vapor,is the gas drag coefficient,is the velocity vector of the condensate, andis the velocity vector of the vapor. A set of simple classical dynamics equations were then integrated to find the position and velocity of the particle as a function of time. Full details of the calculation are given in Appendix A Figure 4 presents the velocity and distance traveled in the direction perpendicular to the rotation axis for small, isolated condensates starting at different points in the synestia shown in Figures 1e–1g. We initialized the condensates at a velocity equal to the gas velocity, as if the particle had just condensed from the vapor. The gas velocity in the structure was calculated using the midplane angular velocity profile of the SPH structure, assuming that the angular velocity did not vary with height above the midplane in accordance with the Poincaré‐Wavre theorem. We assumed the vapor motion was purely azimuthal. In these calculations, the gas density was constant. In the outer regions of the structure the scale height is large (on the order of megameters), and over the timescales shown in Figure 4, the particles do not typically move more than a scale height. The assumption of constant density does not have a significant effect on our results. Figure 4 Open in figure viewer PowerPoint δr xy , g–l) as a function of time after formation, for condensates in the post‐impact synestia in Figures −4 kg m−3, corresponding to vapor on a saturated adiabat at 10−3 bar. The second row shows the dynamics at 5 Mm above the midplane at pressures interpolated from the post‐impact synestia (23 and 5 bar). The corresponding gas densities were 2.1 and 0.6 kg m−3. The bottom row is for condensates formed in the midplane at pressures of 46 and 7 bar, and gas densities of 3.7 and 0.8 kg m−3. The condensates were initialized at the vapor velocity, and the gas density was assumed to be constant in each calculation. We considered condensates with radii of 1 mm (blue), 1 cm (green), and 10 cm (yellow). The density of the condensates was assumed to be 3,000 kg m−3. Condensates formed at the photosphere of a synestia can rapidly spiral inward in cylindrical radius, toward the rotation axis. Shown are the velocities of individual condensates toward the rotation axis (a–f) and the corresponding change in cylindrical radius (, g–l) as a function of time after formation, for condensates in the post‐impact synestia in Figures 1 and 3 . The condensate behavior is demonstrated at two different cylindrical radii (18 and 25 Mm) and at different heights in the structure, with correspondingly different gas densities. The top row shows the case of condensation at the approximate height of the photosphere of the structure (10 and 13 Mm for the two radii) with a gas density of 1.6 × 10kg m, corresponding to vapor on a saturated adiabat at 10bar. The second row shows the dynamics at 5 Mm above the midplane at pressures interpolated from the post‐impact synestia (23 and 5 bar). The corresponding gas densities were 2.1 and 0.6 kg m. The bottom row is for condensates formed in the midplane at pressures of 46 and 7 bar, and gas densities of 3.7 and 0.8 kg m. The condensates were initialized at the vapor velocity, and the gas density was assumed to be constant in each calculation. We considered condensates with radii of 1 mm (blue), 1 cm (green), and 10 cm (yellow). The density of the condensates was assumed to be 3,000 kg m At the photosphere, the residual vapor after condensation is low density (∼10−4 kg m−3) and offers little resistance to the falling particles. Condensates of the size calculated above (a few millimeters to a few centimeters) rapidly accelerate to velocities of several tens to hundreds of meters per second and both fall vertically toward the midplane and spiral in toward the rotation axis (Figure 4). The velocity initially increases rapidly and then plateaus at later times. Condensates originating closer to the midplane behave similarly to those formed at the photosphere, but, due to the higher gas density and hence amplified drag, reach lower terminal velocities (on the order of a few meters per second, Figure 4). The rate of radial infall is greater for larger condensates that experience less gas drag. The AM of the condensates formed from the vapor is such that they begin falling on highly elliptical orbits. The orbits of smaller condensates are rapidly circularized and the particles spiral inward. Particles that experience less gas drag are less perturbed from their original elliptical orbits and fall more rapidly toward the rotation axis. (7) We pause to confirm the typical droplet size using our condensate orbit calculations. Using the form in equation 6 , we once again assume a balance between gas drag and surface tension (8) Rearranging, the droplet size is given by In our calculations, the differential velocity varies depending on droplet size, with larger droplets having greater shear. At the later times shown in Figure 4, the 10 cm bodies reach velocities at which they would be sheared apart by the gas. For the velocities at 104 s in Figure 4, the largest droplets that would travel slowly enough to not shear apart are on the order of a few millimeters to a few centimeters, in good agreement with our earlier size estimate. At earlier times, when the differential velocities are smaller, condensates could have been much larger. The equilibrium droplet size does not vary significantly with height in the structure as increased gas density is compensated for by lower differential velocities. Falling condensates are a radial mass transport mechanism in synestias. Condensates originating at the photosphere rapidly accelerate to velocities comparable to the convective velocities of the vapor (hundreds of meters per second at the photosphere, see section S3), and are likely to avoid entrainment by gas convection. As a result of rapid radial motion, condensates fall at relatively shallow angles to the photosphere. Due to the large‐scale height of the outer regions, condensates can move considerable distances at relatively low pressures and hence low vapor densities. In higher‐density regions of the structure, the relative velocity between condensate and gas is lower, and it is possible that condensates could be entrained in the turbulent fluid. As a result, the bulk of the radial mass transport by condensates likely occurs near the photosphere of the synestia. (9) σ is the Boltzmann constant, T vap and T cond are the temperatures of the surrounding vapor and condensate, respectively, and R cond is the radius of the condensate. Assuming a homogeneous condensate, the energy balance is given by (10) M is the mass, l is the latent heat, c p is the specific heat capacity, and t is time. Assuming that vaporization occurs linearly between the initial temperature of the condensate, , and the temperature of complete vaporization, , (11) , where M 0 is the initial mass of the condensate. This is a reasonable approximation for moderate pressures and initial condensate temperatures close to the liquid‐vapor phase boundary (see section (12) (13) The efficiency of condensates as a radial mass transport mechanism is dependent on how long condensates can survive in the synestia before vaporizing. At the photosphere, condensates are in thermal equilibrium with the vapor and could persist indefinitely. However, as they fall into higher‐pressure regions of the structure, condensates are heated and begin to vaporize. Here we approximate lower limits on the evaporation timescales for individual, isolated condensates. We assume that condensates are heated by blackbody radiation from the surrounding vapor. The net power gained by a spherical particle is given bywhereis the Boltzmann constant,andare the temperatures of the surrounding vapor and condensate, respectively, andis the radius of the condensate. Assuming a homogeneous condensate, the energy balance is given bywhereis the mass,is the latent heat,is the specific heat capacity, andis time. Assuming that vaporization occurs linearly between the initial temperature of the condensate,, and the temperature of complete vaporization,for, whereis the initial mass of the condensate. This is a reasonable approximation for moderate pressures and initial condensate temperatures close to the liquid‐vapor phase boundary (see section 3 ). With this assumption, the temperature increases linearly with mass loss,Using equations 10 and 12 , the condensate mass changes at a rate We solve equation 13, using equations 9 and 11, for the mass of the condensate as a function of time. We present two example calculations of the vaporization of condensates: at pressures close to the photosphere at 10−3 bar and at a midplane pressure of 20 bar. In section 3, we calculate the multicomponent phase diagram for silicates in a terrestrial synestia, and find that, at 10−3 bar, vaporization of initially fully condensed material occurs over a temperature range between 2,300 and 2,700 K (Figure 9a). Based on this calculation, we used ,300 K and ,700 K for our low‐pressure calculations. Similarly, we use ,700 K and ,200 K for the midplane (Figure 9b). We do not consider the time for the condensate to heat up to since it is negligible compared to the vaporization time as the specific heat capacity is 4 orders of magnitude less than the latent heat of vaporization. In the outer regions of a structure during cooling, the vapor is expected to be close to the liquid‐vapor phase boundary (Figure 1). We hence assume reference vapor temperatures similar to , 2750 K at 10−3 bar and 4,250 K in the midplane. We used a latent heat of l = 1.7 × 107 J K−1 kg−1, which has been widely used in lunar disk studies and is consistent with the EOS used in the SPH simulations (e.g., Thompson & Stevenson, 1988; Ward, 2012), a condensate density of 3,000 kg m−3, and a specific heat capacity of 1,000 J K−1 kg−1, a typical value for silicate melts (Stebbins et al., 1984). Figure 5 Open in figure viewer PowerPoint Small, isolated droplets can survive for minutes to tens of minutes in the pure‐vapor region of the synestia. Panels show (a) the loss of mass over time for small condensates of a given initial radius at low pressure and (b) midplane pressures. Mass loss is more rapid at the midplane due to the higher vapor temperature and larger temperature difference between the condensate and vapor. Isolated condensates, of the typical sizes we calculated, can survive on the order of minutes to tens of minutes in the vapor of the synestia (Figure 5). Our estimate for the vaporization timescales of condensates are lower limits, and a number of processes could increase the survival time of condensates. Energy exchange by black body radiation with pure vapor is the most efficient energy transfer process possible. The vapor produced from the condensate is colder than the gas of the surrounding structure. The colder gas will partly insulate the condensate from radiation and thermal exchange with the hotter vapor of the synestia. We did not account for the heating of condensates by thermal diffusion, due to the difficulty of accounting for this blanketing effect. Collisions and combination of small condensates would lower the effective surface area to mass ratio, again reducing the efficiency of vaporization. The short lifespan of small, isolated condensates in the synestia limits their ability to transport mass radially. An isolated condensate at the photosphere could only move on the order of 10 km toward the rotation axis before vaporizing. However, there is a substantial mass of condensate being continuously formed at the photosphere. When a fluid parcel rising in the synestia reaches the photosphere, a large mass fraction condenses almost instantaneously. Due to the high radiative temperature, a 1 km thick photosphere would condense in about a second. Condensates forming at the photosphere are not isolated but are falling as part of larger, condensate‐rich downwellings. Condensates dominate the optical depth of silicate vapor‐melt mixtures (section S2), and the large mass fraction of condensates in downwellings makes them initially optically thick. Thus, the radiation field within the downwelling is dominated by emission from the lower temperature condensates and not the hot vapor of the synestia. Heating from radiation within the mixture is much less efficient than for isolated condensates. Heating can also occur by thermal conduction between the condensates and the vapor that they are falling through. As condensates vaporize, they produce gas that is in thermal equilibrium with the condensate. The gas flow then advects the lost vapor away from the particle. In a downwelling, subsequent condensates are falling into vapor produced by partial vaporization of the leading condensates. The temperature difference between the condensates and the vapor is smaller and the condensates are not efficiently heated. The effects of self‐shielding and self‐buffering by condensates in the downwelling last until the downwelling is broken up by eddies in the fluid. The timescale for this process is uncertain but is likely to be on convective timescales. Thus, the lifetime of condensates in downwellings would be significantly longer than isolated condensates. Due to their longer lifetimes, condensates in condensate‐rich downwellings from the photosphere could transport mass over substantial radial distances. As the radial velocity of condensates increases rapidly with time (Figure 4), a relatively small increase in survival time can significantly increase the radial distance traversed. For example, assuming the infall rates are the same for isolated condensates and those in downwellings, a survival time of 104 s (a few hours) would allow for radial mass transport at the photosphere over 105–106 m. The survival time for isolated condensates that we calculated above cannot be applied to groups of particles. The survival time of condensates in condensate‐rich downwellings is much longer, and we expect that condensates can transport mass over megameter length scales before vaporizing. Condensates could also transport mass in regions where condensates are stable in the midplane (region outside the black lines in Figures 1e and 1f). These condensates would also experience drag and spiral inward toward the rotation axis. As the midplane pressure is higher than at the photosphere (0.1–10 kg m−3), the infall velocities are somewhat smaller (several to tens of meters per second for centimeter‐sized bodies). However, where condensates are thermodynamic stable, they can fall for substantial distances. Outside the Roche limit, mutual collisions between condensates would lead to the accretion of larger bodies. The radial infall of moderate‐sized bodies can be more efficient than for smaller condensates depending on the gas density. Condensates must accrete onto moonlet sized bodies in order to avoid spiraling into higher‐density regions of the structure and revaporizing (section 4.3). A more sophisticated model of condensation and gas drag is needed to fully examine the dynamics of condensates in the synestia. However, our simple calculations suggest that falling condensates are a significant radial mass transport mechanism. During the evolution of synestias, several lunar masses of material condenses at the photosphere (section 2.5) and could be advected substantial distances (e.g., hundreds of kilometers to megameters). Significantly, condensates can move mass perpendicular to the rotation axis, whereas mass transport radially in the vapor is difficult due to the substantial Coriolis force. As condensates are dragged by the gas, they also deposit AM into the vapor. Large‐scale mass and AM transport radially by condensates is a process that has not previously been considered in post‐impact evolution models. The condensation of mass in the synestia is also a major component of the energy budget. The energy transported by condensates and lost by radiation from the photosphere is redistributed vertically by convection. We estimated the vertical convective mixing timescale for the structure using mixing length theory by making an analogy to purely thermal convection in rotating systems (section S3). We found that convective velocities in the outer regions of the synestia are on the order of tens of meters per second in the midplane and hundreds of meters per second at the photosphere. The corresponding mixing timescales are on the order of days, much shorter than the cooling timescale for the disk‐like region, which is on the order of years (section 2.5). The mass in each column of vapor is cooled simultaneously by radiation and condensate transport. Radiation from the photosphere leads to rapid evolution of a synestia. As post‐impact structures cool, the entropy of the vapor decreases, a fraction of the vapor condenses, and the pressure support is reduced. The production of condensates leads to a reduction in the vapor surface density. Material is no longer supported at such large radii and the structure radially contracts. Next, we describe a model for calculating the cooling of synestias and present example cooling simulations.

2.4 Cooling Calculation: Methods Based on the results presented in previous sections, post‐giant‐impact structures must be analyzed as continuous rotating fluids. For the early time evolution after the impact event, a fluid code calculation that approximates the redistribution of condensing material can capture the basic physics of the system. Here we studied the early evolution of post‐impact structures using SPH methods. We constructed a simple model to calculate the physical structure (e.g., mass, pressure, and entropy distribution) of a synestia during radiative cooling. This calculation is intended to assess the timescales for cooling, the potential mass and orbit of a primary satellite, and the magnitude of the vapor pressure around the growing satellite. In this work, we neglect or simplify a number of physical processes and, therefore, take a conservative approach and attempt to estimate the fastest possible timescale for cooling and condensing the Roche‐exterior mass of a synestia. In this section, we describe the methods used in the calculation. Complete details of the code implementation of the model are provided in the supporting information (section S4). We adapted the GADGET‐2 SPH code to calculate the cooling of post‐impact synestias. In LS17, GADGET‐2 was used to calculate the equilibrium structure of corotating planets and the results compared well with a potential field method. LS17 also used GADGET‐2 to generate synthetic synestias, formed by heating isolated planets, as well as impact‐generated synestias. Based on this previous work, GADGET‐2 is able to solve for the pressure field of a synestia of a given distribution of mass, AM, and thermal energy. Because synestias evolve both by mass redistribution by condensates and viscous spreading, they are not static structures. The structure of a synestia with a given mass and AM is not unique and is dependent on the physical processes that generated the synestia and acted during its evolution. Therefore, our calculation includes a series of steps that model the major physical processes controlling the creation and evolution of a terrestrial synestia. First, we generate a synestia by a giant impact that is calculated until a near‐axisymmetric structure is achieved (24 to 48 h). Second, motivated by our discussion of the multiphase dynamics in the synestias during and immediately after the impact in section 2.2, the outer portion of the post‐impact structure is thermally equilibrated. In this step, the fraction of condensates with sufficient AM to be rotationally supported in circular orbits beyond the Roche limit is removed from the SPH calculation under the assumption that this material quickly accretes onto a body that we refer to as the seed of the moon. Third, we mimic radiative cooling by decreasing the thermal energy of SPH particles while accounting for condensation and redistribution of mass and AM. We estimate the mass and orbit of the moon that is formed and determine the range of vapor pressures surrounding the moon from the calculated pressure structure of the cooling synestia. In the first step, giant impacts were modeled in the same manner as in Ćuk and Stewart (2012) and LS17. For this study, we drew examples of impact‐generated synestias from the database presented in LS17. The impacting bodies were differentiated (2/3 rocky mantle, 1/3 iron core by mass) with forsterite mantles and iron cores modeled using M‐ANEOS equations of state (Canup, 2012; Melosh, 2007). The forsterite EOS is a single‐component model, which provides a simple treatment of the liquid‐vapor phase boundary. In section 4, we discuss the effect of a multicomponent phase boundary on the evolution of the structure and the formation of a moon. In the second step, the outer portion of the post‐impact structure was thermally equilibrated. The adjustment of the structure to the changes in entropy is calculated in a manner identical to the third step (cooling) but with no radiative heat loss. Based on the unique post‐impact entropy distribution in the structure, a value of specific entropy was chosen to demarcate between an inner stratified region and a well‐mixed outer region. SPH particles in the well‐mixed region are divided into two groups: an isentropic pure‐vapor group and a vapor‐dome group. The three thermal groups are shown schematically in Figure 6a. Particles in the isentropic vapor group were all assigned the same entropy given by the mass‐weighted mean value of the particles in that group (yellow particles). In the vapor‐dome group, the mass fraction of condensate was removed from each particle using the lever rule. The remaining mass was assigned the density and specific entropy of vapor on the phase boundary at the same pressure (blue particles on the saturated adiabat). Examples of the change in the entropy distribution in modeled synestias are shown in Figures 1d, 1c, 1h, 1g, and S3c–S3f. The position and velocity of the particles were not changed, and hence, they retained the same specific AM, j. Figure 6 Open in figure viewer PowerPoint a R . (a) At the beginning of each time step, particles are assigned to thermodynamic groups. (b) Every particle i in bin k is cooled in proportion to radiative heat loss over the surface area 2A k . Cooling is calculated for each radial bin. (c) The mass falling into bin l is added to all isentropic particles and the enthalpy of the particles is reduced by the latent heat of vaporization. Mass addition and revaporization is repeated for each radial bin. The variables are defined as follows: m i is the original mass of particle i, dm i is the change in mass of particle i upon condensation or upon addition of mass to a bin, is the new mass of particle i upon cooling or addition of mass, S liq and S vap are the specific entropies on the liquid and vapor side of the liquid‐vapor phase boundary, respectively, at the pressure of particle i, is the updated entropy of particle i, dQ k is the energy lost by bin k due to cooling, σ is the Stefan‐Boltzmann constant, dt is the time increment, T eff is the effective radiative temperature, T i is the temperature of particle i, is the radius of the circular Keplerian orbit corresponding to the specific angular momentum of particle k, is the number of bins into which mass is being redistributed from particle i, and is the number of isentropic particles in bin l. A schematic of the radiative cooling model for a synestia. The left column shows the position of a select number of smoothed particle hydrodynamics (SPH) particles in pressure‐specific entropy space (colored points). The liquid‐vapor phase boundary is shown in black. The right column presents a spatial schematic of the cooling synestia as a cross section parallel to the rotation axis. SPH particles in the inner group (gray), isentropic group (yellow), and vapor dome group (blue) are given as colored circles, and the division of radial bins is shown by the black lines. The dashed line indicates the Roche limit,. (a) At the beginning of each time step, particles are assigned to thermodynamic groups. (b) Every particlein binis cooled in proportion to radiative heat loss over the surface area 2. Cooling is calculated for each radial bin. (c) The mass falling into binis added to all isentropic particles and the enthalpy of the particles is reduced by the latent heat of vaporization. Mass addition and revaporization is repeated for each radial bin. The variables are defined as follows:is the original mass of particleis the change in mass of particleupon condensation or upon addition of mass to a bin,is the new mass of particleupon cooling or addition of mass,andare the specific entropies on the liquid and vapor side of the liquid‐vapor phase boundary, respectively, at the pressure of particleis the updated entropy of particle, dis the energy lost by bindue to cooling,is the Stefan‐Boltzmann constant, dis the time increment,is the effective radiative temperature,is the temperature of particleis the radius of the circular Keplerian orbit corresponding to the specific angular momentum of particleis the number of bins into which mass is being redistributed from particle, andis the number of isentropic particles in bin The extracted condensate with specific AM greater than that of a body in a circular, Keplerian orbit at the Roche limit, j Roche , was assumed to accrete into a single body. The mass of condensate with specific AM less than j Roche was distributed evenly in radial bins between the point of origin and the radius corresponding to a circular, Keplerian orbit for the specific AM of the condensate. The redistribution is shown schematically in Figure 6c. The mass of condensate that fell into a given bin was divided equally between each particle in that bin belonging to the isentropic group. The motivation for this simple mass‐distribution function is the expectation that condensates would be vaporized and reincorporated into the structure over a range of radii. Since the details of the dynamics of condensates in the synestia are not yet known, we chose a simple mass‐distribution function in order to investigate the role of condensates in the evolution of the structure. For the purposes of this section, we refer to condensates with j < j Roche as falling condensates. The thermal effect from falling condensates is discussed below. After the extraction of condensates and the redistribution of mass, the physical structure was evolved using the forces calculated by the SPH code under the constraint of the EOS. At each time step the identification of groups and condensate extraction and redistribution were repeated. Typically, after a few time steps the production of condensate became negligible and the structure attained the desired thermal profile. The SPH calculation was continued for a few dynamical times (typically several hours) to form a quasi‐steady structure. In the third step, the specific entropy of the well‐mixed region was reduced in a process that approximates radiative cooling (Figure 6b). The time step in SPH codes is limited by the sound speed and the Courant criterion. In addition, the artificial viscosity in the SPH code causes unrealistically rapid viscous spreading of the structure. Thus, an SPH code cannot be used to directly simulate the 10 to 100 year evolution of a synestia. Here our goal is to estimate the temporal distribution of condensates and the approximate pressure field of the cooling vapor structure. To overcome the timescale issue, radiative cooling was implemented using a large effective radiating temperature, T eff , to determine the energy loss per time step. Then the calculation time was scaled by a factor of (T eff /T rad )4, where T rad is the true radiating temperature, to obtain the corresponding cooling time. The evolution of the synestia was broadly similar for a range of T eff , and we typically used 15,000 or 20,000 K. In each cooling time step, we identified the particle groups and recalculated the entropy of the isentropic group (Figure 6a). Next, we decreased the specific entropy of each well‐mixed particle in each 1 Mm radial bin such that the total enthalpy removed equaled the radiative energy lost from the surface area of that bin. The process is illustrated by the cooling of a single bin, indexed by k, by dQ k in Figure 6b. For example, the ith particle initially on the saturated adiabat will cool, reducing its specific entropy by dS k . As a result, a portion of the particle condenses, dm i . This process is repeated for each particle in each bin. The condensate fraction was removed in the same manner as in the thermal equilibration step. The mass of condensates falling within the Roche limit, a R , was redistributed by adding mass to isentropic group particles and reducing the enthalpy of those particles by an amount determined by the latent heat of vaporization. In Figure 6c, an example particle j in bin l increases in mass by dm j and decreases in entropy by dS j . This procedure is repeated for every particle in bin l. The addition of mass from falling and vaporizing condensates also mimics the transfer of the original condensate's AM to the gas over a range of radii. If the structure cooled to the point where a given Roche‐interior bin no longer contained pure vapor particles, the condensates falling into that bin were removed. In this work, we did not attempt to model the evolution of any thermodynamically stable condensates within the Roche limit, and removal of this material follows our conservative approach of estimating the fastest possible cooling time for the vapor structure. Based on the results of the thermal equilibration and cooling steps, we estimated the mass and angular momentum of a growing moon. We estimate the range of mass and AM for a primary satellite using two different assumptions. The first estimate (A) only includes condensates with specific AM larger than j Roche . The second estimate (B) adds falling condensates that fell beyond the Roche limit according to our simple condensate redistribution scheme. Both estimates of the satellite mass assume perfect accretion of Roche‐exterior condensates into a single body. Perfect accretion is consistent with N‐body simulations that show efficient accretion of condensates beyond the Roche limit (e.g., Ida et al., 1997; Kokubo et al., 2000; Salmon & Canup, 2012, 2014). The growing moon has a large collisional cross section and would accrete some falling condensates, which motivates our inclusion of them for the moon B estimate. The addition of falling condensates increases the total mass and decreases the specific AM of moon B compared to moon A. We calculate the radii of the circular orbits of moons A and B and record the midplane pressure of the vapor structure at those radii. In our model, we make a number of simplifications. Condensates are assumed to perfectly separate from the vapor structure and the gravitational effects of condensates are neglected. Thus, neither the gravitational perturbations from the growing moon on the vapor structure nor the effect of gas drag from the structure on the satellite orbit are included. The gravitational field of the growing moon would increase the local vapor pressure, and the vapor pressure around the growing moon in our calculation is a hard lower limit on the pressure around the moon (section 4.4). Conversely, the pressure at the Roche limit provides a weak upper limit on the pressure around the moon, assuming the gravitational effect of moonlets are negligible. We neglect tidal forces and dynamical resonances, such as Lindblad resonances, as discussed in the supporting information. Internal heating by viscous dissipation is not included in the energy budget in order to determine the fastest cooling time. The contribution from viscous heating to the energy budget is uncertain as discussed in Charnoz and Michaut (2015) who showed that the viscous heating rate is much less than the radiative cooling rate in canonical disks.