## Abstract

The structure and viscous evolution of a post-impact, protolunar disc is examined. The equations for a silicate disc in two-phase (vapour–liquid) equilibrium are employed to derive an analytical solution to vertical structure. Both a vertically mixed phase disc and a stratified disc, where a magma layer exists in the mid-plane surrounded by a vapour reservoir, are considered. The former largely reproduces the low gas mass fraction, *x*≪1, profiles of the disc described in earlier literature that proposed that the disc would hover on the brink of gravitational instability. In the latter, the vapour layer has *x*∼1 and is generally gravitationally stable, while the magma layer is vigorously unstable. The viscous evolution of the stratified model is then explored. Initially, the disc quickly settles to a quasi-steady state with a vapour reservoir containing the majority of the disc mass. The magma layer viscously spreads on a time scale of approximately 3–4 years, during which vapour continuously condenses into droplets that settle to the mid-plane, maintaining the magma surface density in spite of disc spreading. Material flowing inwards is accreted by the Earth; material flowing outwards past the Roche boundary can become incorporated into accreting moonlets. This evolution persists until the vapour reservoir is depleted in approximately 50–100 years, depending on its initial mass.

## 1. Introduction

The leading theory of lunar origin is the giant impact hypothesis (GI), wherein a Mars-sized impactor collides with the Earth during its final stage of accretion [1,2]. The impact launches material into orbit around the Earth, and the Moon subsequently forms from this post-impact disc of debris (see also [3,4]). The collision event itself has been modelled extensively by smooth particle hydrodynamics (SPH) techniques [5–12]. These simulations generally find that torques exerted by the deformed target were instrumental in modifying the debris trajectories following their launch to prevent their re-impact with the Earth. The disc material is hot enough that a significant fraction of the material is vaporized.

Following emplacement, disc material interior to the Roche distance, *r*_{R}∼3*R* (where *R*=6.38×10^{8} cm is the Earth's radius), cannot undergo accretion owing to tidal shear. Instead, material will viscously spread, with some being re-accreted by the Earth while the rest diffuses outwards. Material flowing across the Roche boundary can be incorporated into accreting moonlets, which eventually coalesce into the Moon itself. The diffusion rate depends on the strength of the disc viscosity, *ν*, and early estimates of the spreading time scale for the disc were very short. The surface density of the disc is high; i.e. approximately two lunar masses distributed over an area ∼*π*(*r*^{2}_{R}−*R*^{2}) implies an average density of *σ*∼1.4×10^{7} g cm^{−2}. Ward & Cameron [13] pointed out that the disc could become susceptible to gravitational instability if the dispersion velocity of its constituent components damped below a critical value *c*_{crit}=*πGσ*/*Ω* [14], where *Ω*=(*GM*/*r*^{3})^{1/2} is the orbital frequency of the disc material, with *M*=5.97×10^{27} g being the mass of the Earth. Since tidal stresses shear out the unstable clumps as they try to form, the overall result is to generate an effective viscosity that from theoretical grounds [13] is expected to be of the order of
1.1
This functional form was largely confirmed by numerical simulation reported by Takeda & Ida [15] and Daisaka *et al.* [16]. This implies a short spreading time estimate of the order of approximately a year. As material crosses the Roche boundary, it can begin to quickly accrete into a moon. This rapid formation has been numerically simulated by, for example, Ida *et al.* [17] and Kokubo *et al.* [18].

This time scale estimate was brought into serious question by Thompson & Stevenson [19] (hereafter TS88), who pointed out that when so much dissipation occurs over such a short time it is likely that the entire disc would vaporize. The dissipation rate per unit area in a viscous Keplerian disc is of the order of
1.2
where equation (1.1) has been used in the final form. This would heat the disc until it was hot enough to radiate away the energy at the rate
1.3
where *σ*_{SB}=5.67×10^{−5} erg cm^{−2} s^{−1}K^{−4} is the Stefan–Boltzmann constant, *T*_{ph} is the photospheric temperature, and the lead factor of 2 counts both surfaces of the disc. Equating equations (1.2) and (1.3) implies a photospheric temperature of *T*_{ph}≈6700 K, far too high for condensed silicates, and implying a contradiction. TS88 argued that a much lower photospheric temperature (*T*_{ph}∼2000 K) is more consistent for a disc where vapour and condensed phases coexist, and that this should ultimately limit the allowed rate of spreading (see also [20]), resulting in a time scale of the order of several tens of years.

## 2. Vertical disc structure

To reconcile these two time scales, Ward [21] has recently re-examined the vertical structure of the protolunar disc. For a single-phase gas disc, the constitutive equations are: (i) the ideal gas law, , where *P* is the pressure, is the universal gas constant and *μ* is the molecular weight, for which I use 30 g mol^{−1}; (ii) hydrostatic equilibrium, *ρ*^{−1}d*P*/d*z*=*zΩ*^{2}, where *z* is the vertical height above the mid-plane; and, if it is convecting, (iii) the adiabatic equation, d*P*/d*z*=*c*^{2}d*ρ*/d*z*, where *c* is the gas sound speed.

On the other hand, when two phases are present, there is another state variable, *x*, representing the gas mass fraction, with the condensed-phase mass fraction then being 1−*x* (TS88). One must now also distinguish between the gas density *ρ*_{g} and the total density *ρ* of the gas plus condensed phase mixture given by: (iv) 1/*ρ*=*x*/*ρ*_{g}+(1−*x*)/*ρ*_{s}, where *ρ*_{s} denotes the *body* density of the condensed phase. It is the total density that must appear in the hydrostatic and adiabatic equations, but not in the ideal gas law, where only the gas density *ρ*_{g} is to be used. In addition, when both phases are present, the temperature and pressure must satisfy the Clausius–Clapeyron phase equilibrium equation, which TS88 approximate by the relationship: (v) *P*=*P*_{0}e^{−T0/T}, with *P*_{0}=3×10^{14} dyn, K, where ℓ=1.7×10^{11} erg g^{−1} is the latent heat of vaporization of silicates. These data are primarily due to Krieger [22]. Finally, the sound speed must now be replaced by: (vi) the two-phase sound speed (e.g. [23]; TS88) when both phases can be present,
where , are the specific heats of the vapour and condensed phases, assuming the vapour is predominantly diatomic.

Equation (v) gives the pressure as a function of temperature alone. Combining this with (i) also allows the determination of the gas density as a function of temperature alone,
2.1
Using equation (2.1) to eliminate *ρ*_{g} in the sound speed (vi) results in an expression that depends only on *x* and *T*. The two-phase sound speed is generally less than that of a single phase because the ability of gas to condense makes it more compressible. Figure 1 shows the behaviour of the sound speed as a function of the gas mass fraction *x* for some representative temperatures.

The hydrostatic and adiabatic equations can be reduced to two equations for d*T*/d*z* and d*x*/d*z*, *viz.*
2.2
and
2.3
To obtain a system solution, I first divide equation (2.3) by (2.2). Generally, *ρ*_{g}≪*ρ*_{s}, and the first factor in equation (2.3) can be set to unity. For the cases of interest here, it is also likely that *x*/*ρ*_{g}≫(1−*x*)/*ρ*_{s} in equation (iv), i.e. that gas occupies most of the volume even if the gas mass fraction is low (TS88). In that case, *ρ*/*ρ*_{g}≈1/*x* and the equation's quotient simplifies to
2.4
where Δ*C*≡*C*_{s}−*C*_{g}. This can be integrated from the mid-plane values to give [21]
2.5
where denotes the exponential integral. Equation (2.5) tells us, for a given set of mid-plane values {*x*_{c},*T*_{c}}, how the gas mass fraction varies as a function of temperature, *x*(*T*). We already know *P*(*T*) and *ρ*_{g}(*T*) from equations (v) and (2.1), respectively, but equation (2.5) can now be used to find *ρ*(*T*)=*ρ*_{g}/*x* and *c*(*T*,*x*(*T*)). All that remains is to determine the height *z*(*T*) at which a given temperature is reached. To do this, set *ρ*_{g}/*ρ*≈*x* in equation (2.2), eliminate *x* with equation (2.5), and integrate, leading to
2.6
where the integration limits are *y*≡Δ*CT*/ℓ and *y*_{c}≡Δ*CT*_{c}/ℓ. Integration by parts gives , and combining with (2.6) yields
2.7
which completes the solution.

Note that, for the values adopted here, Δ*CT*/ℓ=(*C*_{s}−*C*_{g})*T*/ℓ∼3.5×10^{−3}(*T*/2000 K) is significantly smaller than unity and suggests that we can expand equations (2.5) and (2.7) in that quantity. To lowest order, one finds [21]
2.8
In writing the r.h.s., a new length scale, , has been introduced. For the parameter values adopted here, *H*_{0}=1.33×10^{9} cm=2.08 *R*_{⊕}. To get comparable accuracy in expanding equation (2.5), one needs only to keep terms up to ln *y* because there is no factor of (Δ*C*)^{−1} in the coefficients, i.e.
2.9
An alternative form for *x* can be found by using (2.8) to eliminate the log term in (2.9), i.e.
2.10
Since the temperature decreases with height, the second term in equation (2.10) increases the gas mass fraction while the last term decreases it. Which dominates depends on how high one must go for the temperature to drop to the photospheric value, *T*_{ph}∼2000 K. Figures 2 and 3 show example disc profiles from Ward [21].

For the sound speed, the level of approximations used to derive (2.9) and (2.8) suggests that we set Δ*C*→0 and ignore the second term in the numerator of equation (vi) for consistency. The sound speed expression then reduces to . The adopted gas mass fraction at the mid-plane, *x*_{c}, is much smaller in figure 2 than in figure 3. Allowable choices of {*x*_{c},*T*_{c}} are constrained through the surface density, , so their mid-plane temperatures must differ as well.

Figure 2 represents a profile of a well-mixed two-phase, marginally unstable disc with most of the mass in the condensed phase, whereas figure 3 is appropriate for a gravitationally stable vapour layer surrounding an unstable magma layer in a stratified disc. In the above derivation, self-gravity of the disc has been ignored. This is only a fair approximation for a disc on the brink of instability, *Q*∼*O*(1), but becomes more justified as *Q* increases in the vapour atmosphere of the stratified disc. In fact, Thompson & Stevenson [19] cautioned that there was some uncertainty as to whether disc mixing could be maintained by convection, the strength of which was estimated from mixing length theory by requiring vertical energy transport sufficient to replenish radiation losses. However, this constraint could also be satisfied by the rain-out of a very small fraction of magma droplets carrying latent heat. Consequently, I feel that the stratified structure is probably more robust and I shall pursue this option henceforth. Such a model was also considered by Machida & Abe [24], who indeed concluded that an unstable layer of condensates would form quickly by settling, surrounded by a stable vapour layer. However, their treatment is actually similar to a two-phase *two*-component disc because the phase equilibrium physics, which is essential to the overall behaviour, was not included. Here, I describe a stratified model including this physics.

## 3. Disc viscous evolution

The governing equations for the surface density, *σ*(*r*,*t*), of the magma layer can be derived by modifying a standard approach for viscous discs, e.g. Lynden-Bell & Pringle [25], to include an active source of supply through the settling of condensate droplets. The resulting magma layer's evolution is characterized by a varying in-plane radial flux of material, *F*(*r*,*t*), that can be found from Euler's equation, *F*=−∂*g*/∂*h*, where *g*=3*πσνh* is the viscous couple, and *h*=(*GMr*)^{1/2} is the specific angular momentum at radius *r* for a Keplerian disc and is convenient to use as the independent variable. This is differentiated to find ∂*F*/∂*h*=−∂^{2}*g*/∂*h*^{2} and is combined with the continuity equation, which in the absence of active supply reads ∂*σ*/∂*t*=−(1/2*πr*)∂*F*/∂*r*=−(*G*^{2}*M*^{2}/4*πh*^{3})∂*F*/∂*h*. If, instead, the disc is surrounded by a condensing gas reservoir, material is continuously added to the magma layer so that −∂*σ*_{V}/∂*t* must be added to the r.h.s. This causes a modification in the divergence of the in-plane flux, *viz.*
3.1

The rate of change of the vapour surface density is controlled by the excess/deficit of energy radiated away at the gas disc surface compared to the energy input from viscous dissipation in the melt layer,
3.2
where is given by equation (1.2). This can be expressed in terms of *g* and *h* by the substitutions *σν*=*g*/3*πh*, *Ω*^{2}=(*GM*)^{4}/*h*^{6}. Finally, the magma surface density can then be expressed in terms of the couple as
3.3
and differentiating this with respect to time ultimately leads to the diffusion equation
3.4
with final source terms on the r.h.s. due to the material exchanged with the vapour reservoir.

### (a) Quasi-steady-state solutions

The active supply enables the magma layer to maintain a quasi-steady state where ∂*g*/∂*t*≈0 and the l.h.s. of equation (3.4) nearly vanishes. The layer tends to be driven to this state from its initial post-impact configuration. Of course, equation (3.2) tells us that this is not the case for the vapour layer. To solve for the magma steady-state profile, it is convenient to cast the defining equation in non-dimensional form by introducing the quantities *h*′=*h*/*h*_{p}, *g*′=*g*/*g*_{0}, *σ*′=*σ*/*σ*_{0}, *ν*′=*ν*/*ν*_{0} and *F*′=*F*/*F*_{0}, where their reference values are defined as
The normalized surface density, viscosity and in-plane flux can then be expressed in terms of *g*′ and *h*′ as *σ*′=*g*^{′1/3}/*h*^{′10/3}, *ν*′=*g*^{′2/3}*h*^{′7/3} and *F*′=−∂*g*′/∂*h*′, and the equation giving the steady state of *g*′ takes the compact form
3.5
Equation (3.5) describes a one-parameter set of curves where the parameter *η*^{2}≡3*GM*/ℓ*R* is proportional to the ratio of the specific potential energy at the planet's surface to the latent heat of vaporization. Normalizing time to years, one can also cast equation (3.2) in a non-dimensional form, ∂*σ*′_{V}/∂*t*′=−(1−*η*^{2}*σ*^{′3}*h*^{′3}).

I now make another change of independent variable to *z*≡*η*/*h*′ that converts equation (3.5) into a non-homogeneous, modified spherical Bessel equation of order zero,
3.6
which can be solved analytically by combining homogeneous solutions, *g*′_{1}=e^{z}/*z*, *g*′_{2}=e^{−z}/*z* with the particular solution found from , where *W*(*z*) is the Wronskian, *W*=*g*′_{1}d*g*′_{2}/d*z*−*g*′_{2}d*g*′_{1}/d*z*=−2/*z*^{2}. Combining homogeneous and particular solutions leads to
3.7
where I have defined . Free accretion of disc material by the planet requires *g*(*z*_{p})=0, implying *A* e^{zp}=−*B* e^{−zp}≡*C*.

### (b) Freely expanding discs

I will denote the outer edge of the disc by *r*_{d} for which *z*_{d}=*η*/*h*′_{d}. If the disc extends to the Roche boundary (*r*_{d}=*r*_{R}) and can freely expand across it, that material will begin to mutually accrete, converting to a low optical depth particle layer beyond *r*′_{R} that can in turn be accreted by the developing moon(s). For such a case, I set *g*′(*z*_{d})=0 at the outer disc edge as well, yielding
3.8
for which
3.9
At the inner boundary, *E*(*z*_{p})=*E*(−*z*_{p})=0, , and the flux *F*′=*η*^{−1}*z*^{2} d*g*′/d*z* onto the planet reads
3.10
while the flux across the Roche boundary is
3.11

For the protolunar disc parameters, *η*=3.33, and *z* ranges from *z*_{d}=*η*/*h*′_{d}=1.92 to *z*_{p}=*η*/*h*′_{p}=3.33. The solution using these values is shown in figure 4. The viscous couple is rather symmetrical around the mid-point of the disc. By contrast, the viscosity is strongly peaked at a larger distance, while the density is peaked at a smaller *r*. There is a stagnation point in the flux, *F*′=0, near the halfway mark, with a positive value exterior to that indicating material flowing towards the accretion zone of the Moon. Negative values of *F*′ in the inner portion of the disc indicate material flowing towards and destined to be accreted by the Earth. The fluxes at the planet boundary and at the Roche distance are computed by evaluating *F*′=−d*g*′/d*h*′ at *h*′=1 and , respectively, which gives *F*′_{p}=−0.610 and *F*′_{R}=1.129. Thus, the amount of magma leaving the disc through both boundaries is (*F*′_{R}−*F*′_{p})*F*_{0}=9.51×10^{16} g *s*^{−1}=0.0408 *m*_{moon} yr^{−1}. The mass of the magma layer is maintained at
3.12
It takes 3.4 years to completely replace a magma layer. However, at the beginning, most of the disc material is in the vapour state surrounding the magma layer and continues to condense,^{1} so it requires *m*_{disc}/*m*_{magma}∼7.5(*m*_{disc}/*m*_{moon}) ‘cycles’ to deplete the disc. Assuming a total vapour plus magma mass of approximately two lunar masses requires 15 cycles, and a time of approximately 50 years to deplete the disc. This is roughly the time scale for the disc to radiate all of the latent heat released by condensation at the two 2000 K photospheric surfaces.

### (c) Confined discs

Setting the viscous couples at both boundaries to zero implicitly assumes that the disc is subject to no external torques. However, many impact simulations launch significant material outside of the Roche zone at the outset [9,11,12]. Exterior condensed material can accrete into a lunar embryo on a very short time scale, i.e. days to months [26]. As the forming moon gains mass, it begins to exert a significant torque on the disc that opposes its outward spreading, and this can decrease material outflow. Although angular momentum is still transferred outwards, most of it is absorbed by the embryo, which in turn recoils from the disc. A satellite's torque on a disc is due to resonant interactions and can be characterized by a torque density, . If resonances are close enough that their wavelengths overlap, one can use [27]
3.13
as a good approximation, where *m*, *a*_{s}, *Ω*_{s} are respectively the moon's mass, semi-major axis and mean motion, and *f*∼2.5 is a constant. This modifies the in-plane flux to , where . Because the torque density is proportional to the surface density, *σ*∝*g*^{1/3}/*h*^{10/3}, the steady-state equation becomes nonlinear in *g* and no longer reduces to a modified Bessel equation. However, a solution is still easily found by numerically integrating the flux equation along with the steady-state version of equation (3.1), which in normalized form are
3.14a
and
3.14b
where . Figure 5 illustrates the disc structure for embryos with various fractional lunar masses, *m*/*m*_{moon}, located at *a*_{s}=4*R* as an example. One finds that the outward flux at the Roche boundary diminishes with increasing *m* and eventually shuts off near a tenth of a lunar mass. Further mass increase causes the effective edge of the disc to lie at smaller radii, but has only a minor effect on the in-plane flux profile.^{2}

The total reaction torque, , from the confined disc will cause an embryo to recoil from the disc at a rate , where *r*_{2:1}=*a*_{s}/2^{2/3} denotes the location of the innermost 2:1 resonance. diminishes with increasing embryo distance because fewer resonances fall in the disc. In a quasi-steady state, we can evaluate the torque indirectly by setting it equal to
3.15
where
3.16
are the rate of angular momentum being deposited into the magma layer from the vapour reservoir and the change in the angular momentum content of the quasi-steady-state layer configuration, respectively. One finds, for instance, that the time for a 0.5 lunar mass embryo to recoil to *r*=4.76*R*, where its 2:1 resonance falls just outside the Roche boundary, is of the order of a few years. During this interval, there is a hiatus in the lunar accretion process [26].

## 4. Summary

The evolution of an impact-generated protolunar disc has been solved under the assumption that it achieves quasi-steady state not long after its initial emplacement. It is found that the behaviour exhibits the following characteristics:

— The overall time scale for the protolunar disc evolution is set by the time to radiate its latent heat, as in Thompson & Stevenson [19].

— If it has a stratified structure [21], the magma layer spreads on a time scale of a few years, as in Ward & Cameron [13]. However, a quasi-steady-state magma layer contains only approximately 14% of a lunar mass, so most of the disc material is in the vapour reservoir once a steady state is achieved.

— A freely expanding magma layer can supply material to the accretion zone at the rate of approximately 0.027 lunar masses per year.

— However, as the lunar embryo accretes from this outflow (or from substantial material initially outside the Roche zone), it will shut off the flow before attaining lunar size.

— There is then a hiatus of several years during which the embryo recoils from the disc and disc material only flows onto the Earth.

— Once the resonances from the embryo are no longer in the disc, the disc can re-expand to the Roche boundary, material resumes flow into the accretion zone and the process continues.

— Eventually, the vapour reservoir is exhausted, and a quasi-steady state is no longer possible. The disc continues to expand, but the surface density diminishes and the disc cools below the phase equilibrium temperature. The viscosity also decreases and becomes progressively dominated by individual particle dynamics rather than by collective forces [28].

The accretion of the Moon seems to involve three rather distinct stages, as first identified by Salmon & Canup [26]. Following the giant impact, a mixed vapour–condensate disc is launched into orbit around the Earth. At first, the disc is neither in hydrostatic nor phase equilibrium, but begins to settle towards these states. If significant material initially resides outside the Roche limit in condensed form, it will rapidly accrete as it settles, forming a lunar embryo. If this embryo is massive enough, it will shut off the outward flow from the Roche interior disc until it recoils enough that its principal resonances lie outside of the disc. Alternatively, if most of the material external to the Roche zone is vapour, then formation of lunar embryo(s) will be on the longer time scale required for condensation to occur and some incorporation of interior disc material can directly occur. Either way, the ultimate time scale of complete lunar formation will then depend on the more extended interval required for the mutual accretion of any embryos spawned during the lifetime of the inner disc.

## Funding statement

The author acknowledges the support of NASA's LASER program and the NASA Lunar Science Institute (NLSI) during the performance of the supporting research.

## Acknowledgements

The author thanks the Royal Society and Prof. David Stevenson, Caltech, for his invitation to the conference on the *Origin of the Moon* and for publication support for this communication.

## Footnotes

↵† Dedicated to the memory of Prof. A. G. W. Cameron, Harvard University, framer of the giant impact hypothesis of lunar origin.

One contribution of 19 to a Discussion Meeting Issue ‘Origin of the Moon’.

↵1 Although it is highly unlikely that the initial post-impact disc will match this particular vapour/solid partitioning, it can be approached either by evaporation in an initially more condensed state (typical of canonical Mars-sized, low-velocity impactors, e.g. [9]) or by condensation from a more vaporized one (obtained by larger and/or higher-velocity impactors [11,12].

↵2 It is assumed here that only the magma layer is subjected to the density wave torque, and that the vapour continues to condense throughout the region. Any density waves excited in the vapour reservoir will be pressure waves that propagate towards the Earth and reflect at the Earth–disc interface, creating a standing wave pattern. Some net torque will occur if there is some wave dissipation, but will probably be relatively small compared to that felt by the magma.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.