ON THE VERTICAL STRUCTURE OF THE PROTOLUNAR DISK

A useful basis of comparison is the vertical structure of a hypothetical isothermal, single-phase silicate gas disk with a surface density σ = 10⁷ g cm⁻². The equations that govern the structure are the ideal gas equation

$$\begin{equation} P = \rho RT/\mu, \end{equation} \tag{ 5 }$$

where P(z) is the pressure, T is the temperature, ρ(z) is the gas density, and R = 8.31 × 10⁷ erg mol⁻¹ K⁻¹ is the universal gas constant, and the equation of hydrostatic equilibrium,

$$\begin{equation} \rho ^{ - 1} {\it dP}/dz = - z\Omega ^2, \end{equation} \tag{ 6 }$$

where z is the vertical distance above the disk mid-plane. For a constant temperature T = T_* = 2000 K, the sound speed-like quantity c_* = (RT_*/μ)^1/2 = 7.44 × 10⁴ cm s⁻¹. Substituting Equation (5) into (6) and integrating leads to

$$\begin{equation} \rho = \rho _* e^{ - (z/H_*)^2 }, \end{equation} \tag{ 7 }$$

where ρ_* is the mid-plane density and

$$\begin{equation} H_* = \sqrt 2 c_* /\Omega = 2.40 \times 10^8 \,{\rm cm} \end{equation} \tag{ 8 }$$

is the scale height.

Integrating ρ vertically yields the surface density, i.e.,

$$\begin{equation} \sigma = \int_{ - \infty }^\infty {\rho dz = 2\int_0^\infty {\rho dz} }. \end{equation} \tag{ 9 }$$

Substituting Equation (7) for ρ yields $\sigma = \sqrt \pi \rho _{\ast} H_{\ast}$, implying ρ_* = 2.35 × 10⁻² g cm⁻³. Combining Equations (5) and (7) then gives

$$\begin{equation} P = \rho _ * c_ * ^2 e^{ - (z/H_ *)^2 } = P_ * e^{ - (z/H_ *)^2 } \end{equation} \tag{ 10 }$$

and a mid-plane pressure of P_* = 1.3 × 10⁸ dynes.

Applicability of this type of structure to the lunar disk would depend on two implicit assumptions: (1) an opacity so low as to preclude any significant temperature gradient and (2) a disk too hot for the stability of a condensed phase. Neither of these conditions is likely to have prevailed in the post-impact silicate disk.

When the disk has a finite opacity, the temperature is no longer constant and decreases with height above the mid-plane z = 0. If the opacity is high enough, the temperature profile needed to satisfy the luminosity equation would be super-adiabatic. In this case, it is likely that the vertical heat flux is carried by convective transport and that the adiabatic equation

$$\begin{equation} {\it dP}/dz = c^2 \;d\rho /dz \end{equation} \tag{ 11 }$$

applies as well, where again $c = \sqrt {\gamma RT{\rm /}\mu }$ is the sound speed. Combining with Equations (5) and (6) leads to the well-known adiabatic profiles where P∝ρ^γ. However, as stressed by Thompson & Stevenson (1988), if the both gas and condensed phase co-exist, there are additional constitutive equations to consider.

When two phases are present, there is a new state variable, x, representing the gas mass fraction; the condensed phase mass fraction then being 1 − x. One must now also distinguish between the gas density ρ_g and the total density ρ of the gas plus condensed phase mixture given by

$$\begin{equation} 1/\rho = x/\rho _g + (1 - x)/\rho _s,. \end{equation} \tag{ 12 }$$

where ρ_s denotes the density of the condensed phase. The total density must appear in the hydrostatic (Equation (6)) and adiabatic (Equation (11)) equations, but not in the ideal gas law (Equation (5)), 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. TS88 approximate this by the relationship

$$\begin{equation} P = P_o e^{ - (T_o /T)}, \end{equation} \tag{ 13 }$$

with P_o = 3 × 10¹⁴ dynes and T_o = μℓ/R = 6 × 10⁴ K, where ℓ = 1.7 × 10¹¹ erg g⁻¹ is the latent heat of vaporization of silicates. These are the data primarily employed by TS88 and are due to Krieger (1967). We shall use them here as well to facilitate comparisons. We should caution that Equation (13) is an approximation because "silicates" do not volatilize congruently, e.g., SiO is a major gas phase, but there is no SiO solid.

Equation (13) gives the pressure as a function of temperature alone. Combining Equations (5) and (13) also allows the determination of the gas density as a function of temperature alone,

$$\begin{equation} \rho _g = \left({\frac{{\mu P_o }}{{RT_o }}} \right)\frac{{T_o }}{T}e^{ - T_o /T} = \frac{{P_o }}{\ell }\frac{{T_o }}{T}e^{ - T_o /T}. \end{equation} \tag{ 14 }$$

Finally, the sound speed appearing in Equation (11) must now be replaced by the two-phase sound speed (e.g., Landau & Lifshitz 1959; TS88) when both species can be present

$$\begin{equation} c = \frac{{x + (1 - x)(\rho _g {\rm /}\rho _s)}}{{[(\mu /RT - 2/\ell + C_g T/\ell ^2)x + (C_s T/\ell ^2)(1 - x)]^{1/2} }}, \end{equation} \tag{ 15 }$$

where C_g = (3/2)RT = 4.16 × 10⁶ erg g⁻¹ K⁻¹, C_s ∼ 10⁷ erg g⁻¹ K⁻¹ are the specific heats of the vapor and condensed phases. Using Equation (14) to eliminate ρ_g results in an expression for the sound speed 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 2 shows the behavior of the sound speed as a function of the gas mass fraction x for some representative temperatures.

Zoom In Zoom Out Reset image size

Equations (5), (6), (11), (12), (13), and (15) close the system, and the hydrostatic and adiabatic equations can be reduced to two equations for dT/dz and dx/dz, namely,

$$\begin{equation} \frac{1}{T}\frac{{dT}}{{dz}} = {-} \left({\frac{\rho }{{\rho _g }}} \right)\frac{{z\Omega ^2 }}{\ell } \end{equation} \tag{ 16 }$$

$$\begin{equation} \frac{{dx}}{{dz}} = \!\left({\frac{{\rho _s }}{{\rho _s - \rho _g }}} \right)\!\left({\frac{\rho }{{\rho _g }}} \right)\!\left(\!{ - x\left[\! {1 + \frac{{(C_s - C_g)T}}{\ell }} \right] + \frac{{C_s T}}{\ell }} \!\right)\frac{{z\Omega ^2 }}{\ell }. \end{equation} \tag{ 17 }$$

The details of this derivation are given in Appendix A. Generally, ρ_g ≪ ρ_s and the first factor in Equation (17) can be set to unity. For the cases of interest here, it is also likely that x/ρ_g ≫ (1 − x)/ρ_s in Equation (12), 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 equations simplify to

$$\begin{equation} \frac{1}{T}\frac{{dT}}{{dz}} = - \frac{{z\Omega ^2 }}{{x\ell }} \end{equation} \tag{ 16′ }$$

$$\begin{equation} \frac{{dx}}{{dz}} = - \left({1 + \frac{{(C_s - C_g)T}}{\ell } - \frac{{C_s T}}{{x\ell }}} \right)\frac{{z\Omega ^2 }}{\ell }. \end{equation} \tag{ 17′ }$$

To obtain a system solution we begin by dividing Equation (17') by (16') to find

$$\begin{equation} \frac{{dx}}{{dT}} - x\left({\frac{1}{T} + \frac{{\Delta C}}{\ell }} \right) = \frac{{C_s }}{\ell },\quad \end{equation} \tag{ 18 }$$

where ΔC ≡ C_s − C_g. This can be integrated from the mid-plane values to give

$$\begin{eqnarray} x &=& \frac{T}{{T_c }}\bigg({x_c e^{\Delta C(T - T_c)/\ell } + \frac{{C_s T_c }}{\ell }[E_1 (\Delta CT/\ell)} \nonumber\\ &&- {E_1 (\Delta CT_c /\ell)]e^{\Delta CT/\ell } } \bigg), \end{eqnarray} \tag{ 19 }$$

where $E_1 (y) \equiv \smallint _y^\infty (e^{ - y} /y)dy$ denotes the exponential integral (e.g., Abramowitz & Stegun 1966). Equation (19) 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 (13) and (14), respectively, but Equation (19) 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, we substitute Equation (19) into Equation (16'), and integrate leading to

$$\begin{eqnarray} &&\frac{{\ell x_c }}{{\Delta CT_c }}(e^{\Delta C(T - T_c)} - 1) + \frac{{C_s }}{{\Delta C}}\bigg({\int_{y_c }^y {e^y E_1 (y)dy}} \nonumber\\ &&\quad - {{E_1 (\Delta CT_c /\ell)(e^{\Delta CT/\ell } - e^{\Delta CT_c /\ell })} } \bigg) = - \frac{1}{2}\frac{{z^2 \Omega ^2 }}{\ell }, \end{eqnarray} \tag{ 20 }$$

where the integration limits are y ≡ ΔCT/ℓ, y_c ≡ ΔCT_c/ℓ. Integration by parts gives

$$\begin{equation} \int_{y_c }^y {e^y } E_1 (y)dy = e^y E_1 (y) - e^{y_c } E_1 (y_c) + \ln (y/y_c) \end{equation} \tag{ 21 }$$

and combining with Equation (20) yields

$$\begin{eqnarray} &&\frac{{\ell x_c }}{{\Delta CT_c }}(e^{\Delta C(T - T_c)/\ell } - 1) + \frac{{C_s }}{{\Delta C}}({e^{\Delta CT/\ell } [E_1 (\Delta CT/\ell)} \nonumber\\ &&\quad -{ E_1 (\Delta CT_c /\ell)] + \ln (T/T_c)} ) = - \frac{1}{2}\frac{{z^2 \Omega ^2 }}{\ell }, \end{eqnarray} \tag{ 22 }$$

which completes the solution.

We note that for the values adopted here, ΔCT/ℓ = (C_s − C_g)T/ℓ = 0.0687(T/2000 K) is significantly smaller than unity and suggests that we can expand Equations (19) and (22) in that quantity. Setting e^y ≈ 1 + y + ⋅⋅⋅, E₁(y) ≈ −γ − ln y + y + ⋅⋅⋅, where γ = 0.5772 is Euler's constant, and keeping terms up to order y, Equation (22) reduces to

$$\begin{equation} \left({x_c + \frac{{C_s T_c }}{\ell }} \right)\!\left({\frac{{T - T_c }}{{T_c }}} \right) - \frac{{C_s T_c }}{\ell }\left({\frac{T}{{T_c }}} \right)\ln \left({\frac{T}{{T_c }}} \right) = - \left({\frac{z}{{H_o }}} \right)^2. \end{equation} \tag{ 23 }$$

In writing the right-hand side, a new length scale, H_o ≡ (2ℓ)^1/2/Ω =(2RT_o/μ)^1/2/Ω has been introduced. For the parameter values adopted here, H_o = 1.33 × 10⁹ cm = 2.08 R_⊕. Note that to this level of accuracy, ΔC does not appear. To get comparable accuracy in expanding Equation (19) we need only keep terms up to ln y because there is no factor of (ΔC)⁻¹ in the coefficients. This results in

$$\begin{equation} x = \frac{T}{{T_c }}\left[ {x_c - \frac{{C_s T_c }}{\ell }\ln \left({\frac{T}{{T_c }}} \right)} \right]. \end{equation} \tag{ 24 }$$

An alternative form for x can be found by using Equation (23) to eliminate the log term in Equation (24), i.e.,

$$\begin{equation} x = x_c - \frac{{C_s T_c }}{\ell }\left({\frac{{T - T_c }}{{T_c }}} \right) - \left({\frac{z}{{H_o }}} \right)^2. \end{equation} \tag{ 25 }$$

Since the temperature decreases with height, the second term in Equation (25) 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 3–5 show example disk profiles constructed from Equations (13), (14), (23), and (24) and ρ ≈ ρ_g/x. For the sound speed, the level of approximations used to derive Equations (24) and (25) suggests that we set ΔC → 0 and ignore the second term in the numerator of Equation (15) for consistency. The sound speed expression then reduces to

$$\begin{equation} c \approx \frac{{x\sqrt \ell }}{{[(T_o /T - 2)x + (C_s T_{} /\ell)]^{1/2} }}. \end{equation} \tag{ 26 }$$

The adopted gas mass fraction at the mid-plane, x_c, is progressive larger from Figures 3–5. However, for a given surface density, the choice of {x_c, T_c} is constrained through Equation (9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The surface density in Equation (9) involves an integration over the altitude, z. However, in our solution for the profile, the temperature has been employed as the independent variable. One can switch to an integration over T as follows:

$$\begin{eqnarray} \sigma &\approx & 2\int_0^{z_{ph} }\! {\rho dz} = 2\int_{T_c }^{T_{ph} } \!{\rho \frac{{dz}}{{dT}}} dT = 2\int_{T_c }^{T_{ph} } \!{\frac{{\rho _g }}{x}} \left(\!{ - \frac{{x\ell }}{{zT\Omega ^2 }}} \right)dT \nonumber\\ &=& \frac{{2\ell }}{{\Omega ^2 }}\int_{T_{ph} }^{T_{} c} {\frac{{\mu P}}{{RT}}} \frac{{dT}}{{zT}} = \frac{{2P_o T_o }}{{\Omega ^2 }}\int_{T_{ph} }^{T_c } {\frac{{e^{ - T_o /T} }}{{z(T)}}} \frac{{dT}}{{T^2 }}. \end{eqnarray} \tag{ 27 }$$

In the first step, we stop the integration at the photospheric height z_ph = z(2000 K) with little error since the density has become very low and there is little mass beyond that point. The second step switches to T as the integration variable, while step three replaces ρ = ρ_g/x and uses Equation (16') to remove the derivative. Interestingly, the gas mass fraction cancels out at this point. In the final two steps, Equations (1) and (13) to used to eliminate ρ_g and P, respectively. In the integration, Equation (23) is to be used for z(T) which contains the mid-plane gas mass fraction, x_c, as part of the expression so that ultimately σ = σ(x_c, T_c). Requiring σ(x_c, T_c) = 10⁷ g cm⁻² then constrains the allowable combinations as shown in Figure 6. In particular, assumed values of x_c =(0.01, 0.1, 1) require mid-plane temperatures of T_c = (3644 K, 3898 K, 4218 K) respectively, and these were used in the construction of the disk profiles.

Zoom In Zoom Out Reset image size

Although the roughly factor of two change in temperature from the mid-plane to the photosphere occurs over an altitude z_ph comparable to or greater than H_*, the e-folding distance for the pressure (and density) is much less, especially at low values of x_c. This is due in part to phase equilibrium as mandated by the Clausius–Claperyon equation. The pressure at mid-plane is $P_c = P_o e^{ - T_o /T_c }$, while the temperature, T_e, at which the pressure drops by a factor e⁻¹ is given by −T_o/T_e = −T_o/T_c − 1, implying a temperature change of $\Delta \tilde{T} _e \equiv (T_e - T_c)/T_c = - T_c /(T_o + T_c),$ which is typically of order ∼0.06. This suggests that we can set $T/T_c = 1 + \Delta \tilde{T}$ where $\Delta T \equiv (\tilde{T} - T_c )/T_c$ and Taylor expand Equation (23) to find

$$\begin{equation} \Delta \tilde{T} \left({x_c - \frac{{C_s T_c }}{\ell }\Delta \tilde{T} } \right) \approx - \left({\frac{z}{{H_o }}} \right)^2, \end{equation} \tag{ 28 }$$

which, solving for T, gives

$$\begin{equation} T \approx T_c - \frac{{x_c \ell }}{{C_s }}\left({\left[ {1 + \frac{{2C_s T_c }}{{x_c \ell }}\frac{1}{{x_c }}\left({\frac{z}{{H_o }}} \right)^2 } \right]^{1/2} - 1} \right). \end{equation} \tag{ 29 }$$

Using this in Equation (25) then implies

$$\begin{equation} x \approx x_c \left[ {1 + \frac{{2C_s T_c }}{{\ell x_c^2 }}\left({\frac{z}{{H_o }}} \right)^2 } \right]^{1/2} - \left({\frac{z}{{H_o }}} \right)^2. \end{equation} \tag{ 30 }$$

In practice, however, the second term in the bracket of Equation (28) is of order O(10⁻²) for $\Delta T\sim O(\Delta \tilde{T} _e)$ so that keeping only the leading linear term is reasonably good unless x_cis very small. In this case, we can write Equations (28) and (29) to lowest order as

$$\begin{eqnarray} T &\approx& T_c \left[ {1 - \frac{1}{{x_c }}\left({\frac{z}{{H_o }}} \right)^2 } \right],\nonumber\\ x &\approx & x_c \left[ {1 + \left({\frac{{C_s T_c }}{{\ell x_c }} - 1} \right)\frac{1}{{x_c }}\left({\frac{z}{{H_o }}} \right)^2 } \right]. \end{eqnarray} \tag{ 31 }$$

From these it is easy to see that while the temperature decreases with height, whether the gas mass fraction decreases or increases from the mid-plane depends on whether x_c is greater or less than C_sT_c/ℓ.

We return now to the Clausius–Clapeyron equation and write

$$\begin{eqnarray} P &=& P_o \exp (- T_o /T) = P_o \exp [ - (T_o /T_c)(1 + \Delta \tilde{T} )^{ - 1} ] \nonumber\\ & \approx & P_o \exp [(T_o /T_c)(1 - \Delta \tilde{T} )] = P_c \exp [(T_o /T_c)\Delta \tilde{T} ] \nonumber\\ &\approx & P_c \exp \big [ - x_c^{ - 1} (z/H_c)^2 \big], \end{eqnarray} \tag{ 32 }$$

where another new length scale H_c ≡ (2RT_c/μ)^1/2/Ω is introduced. The pressure scale height is of order $H_p = \sqrt {x_c } H_c$ so that decreasing the gas mass fraction thins the disk vertically as is evident in the disk profiles of Figures 3–5. Since (H_c/H_o)² = T_c/T_o, x⁻¹_c(z/H_o)² can still be small when (z/H_p)² ∼ O(1). An estimate of the density is now obtained from ρ = μP/xRT, namely,

$$\begin{eqnarray} \rho & \approx & \frac{{\mu P_c }}{{x_c RT_c }}\nonumber\\ &\times &\frac{{\exp \big[ - x_c^{ - 1} (z/H_c)^2 \big]}}{{\left({1 + \left({\dsty\frac{{C_s T_c^2 }}{{\ell x_c^{} T_o }} - 1} \right)\dsty\frac{1}{{x_c }}\left({\dsty\frac{z}{{H_p }}} \right)^2 } \right)\left({1 - \dsty\frac{{T_c }}{{T_o }}\left({\dsty\frac{z}{{H_p }}} \right)^2 } \right)}}.\nonumber\\ \end{eqnarray} \tag{ 33 }$$

Equation (33) can be used to derive an approximate analytic expression for the surface density $\sigma = 2\smallint _0^\infty \rho dz.$ This integral can be estimated by the method of steepest descent (Appendix B),

$$\begin{equation} \sigma \approx \left({\frac{\pi }{{x_c }}} \right)^{1/2} \frac{{\mu P_c H_c }}{{RT_c }}\left[ {1 + \left({\frac{{C_s T_c }}{{x_c \ell }} - 2} \right)\frac{{T_c }}{{T_o }}} \right]^{ - 1/2}, \end{equation} \tag{ 34 }$$

which can be rearranged to give x_c = x_c(σ, T_c),

$$\begin{eqnarray} x_c & \approx & \left[ {\left({\frac{{2\pi }}{\ell }} \right)\left({\frac{{P_o }}{{\sigma \Omega }}} \right)^2 \frac{{T_o }}{{T_c }}e^{ - 2T_o /T_c }}\right. \nonumber\\ &&\left.- {\left({\frac{{C_s T_c }}{\ell }} \right)\frac{{T_c }}{{T_o }}} \right]\bigg/\left({1 - 2\frac{{T_c }}{{T_o }}} \right). \end{eqnarray} \tag{ 35 }$$

This is compared with the numerical integration in Figure 6. The agreement is reasonably good but becomes less accurate at low x-values where the linear expressions of Equation (31) are less suitable.

Figure 3 illustrates a disk with x = 0.01 and mid-plane temperature T_c = 3644 K. Figure 3(a) shows the temperature (abscissa) normalized to T_* = 2000 K versus height (ordinate) normalized to the scale height, H_* of the isothermal reference solution. The temperature decreases with z at a nearly constant rate. Figure 3(b) displays the gas mass fraction that increases steadily with height until reaching between 7% and 8% at the photosphere. In contrast, the pressure and gas density in Figure 3(c) decreases much more precipitously, as does the full density in Figure 3(d). However, whereas the gas density at the mid-plane is smaller than that of the isothermal model, the total density is larger. This latter reflects the fact that the effective scale height is much less than H_* as discussed in the last section. For comparison, the density for the isothermal disk is included as a dotted line in the graphs. Figure 3(e) plots the sound speed behavior (normalized to c_*) with height. The sound speed increases through most of the disk, but hits a maximum near z/H_* ∼ 0.8. Figure 3(f) displays the cumulative surface density with altitude. The bulk of the surface density resides within ∼0.1H_*, consistent with the total density profile.

Thompson & Stevenson (1988) concentrate on models with low gas fractions x ⩽ O(10⁻²), which they liken to a foam. Their strategy is to adopt a gas fraction low enough to trigger marginal gravitational instability. This is schematically represented by the TS-label in Figure 2. The mid-plane sound speed when 1 − x ∼ O(1) reads $c_c \approx x_c \ell /\sqrt {C_s T_c } \approx 9x_c \times 10^5 \,{\rm cm}\,{\rm s}^{ - {\rm 1}}$. This is below the critical sound speed equation (1) for x_c ⩽ 5.3 × 10⁻³. For this value, the temperature is T_c = 3588 K and the central pressure and gas density are $P_c = P_o e^{ - T_o /T_c } = 1.64 \times 10^7 \,{\rm dynes}$ and ρ_{g, c} = μP_c/RT_c = 1.61 × 10⁻³ g cm⁻³, respectively. The length scale H_c = 3.22 × 10⁸ cm, but the pressure scale height is only $H_p = \sqrt {x_c } H_c = 2.35 \times 10^7 \,{\rm cm}$. In this case, the corresponding central total density is ρ_c = ρ_{g, c}/x_c = 0.303 g cm⁻³. Comparing central quantities to the isothermal reference model, T_c/T_* = 1.79, P_c/P_* = 0.126, ρ_g/ρ_* = 0.0685, ρ_c/ρ_* = 12.9, and H_p/H_* = 0.098. As in Figure 3, the gas density is less, but the total density is over an order of magnitude greater than the reference profile. Again, the preponderance of the mass is in the liquid phase.

To understand how TS88 propose a metastable state can lower the effective viscosity, consider the dispersion relationship for horizontal oscillations of a self-gravitating, Keplerian disk (e.g., Goldreich & Ward 1973),

$$\begin{equation} \omega ^2 = \Omega ^2 - 2\pi G\sigma k + c^2 k^2, \end{equation} \tag{ 36 }$$

where ω is the oscillation frequency and k is the wave number of the mode. The minimum of this expression is where dω²/dk = 0 for which k_min = πGσ/c² and ω²(k_min) = Ω²[1 − (πGσ/cΩ)²]. Marginal stability occurs when ω²(k_min) = 0, i.e., at c_crit → πGσ/Ω and k_crit = Ω²/πGσ. For c < c_crit, ω² = −Ω²[(c_crit/c)² − 1] and ω = ±iΩ[(c_crit/c)² − 1]^1/2 becomes imaginary, implying a growing mode on a timescale τ_inst = Ω⁻¹[(c_crit/c)² − 1]^−1/2. This is longer than the timescale ∼O(Ω⁻¹), to shear out clumps, suggesting that the perturbation velocities, v_inst, generated by the instability achieve a value less than c_crit by a factor ∼1/Ωτ_inst. This implies a viscosity, ν_eff ∼ (v_inst)²/Ω, less than the nominal value of Equation (2) by a factor 1/(Ωτ_inst)² = (c_crit/c)² − 1. Setting this equal to 8.0 × 10⁻³ results in a viscous spreading timescale comparable to the radiation timescale and implies c = 0.996c_crit.

This is an elegant concept, but one caveat of the model is that the liquid and vapor phases must be well mixed throughout the structure. Thompson & Stevenson (1988) caution that this may not be satisfied and that some settling could occur. This could lead to another class of models that are stratified with a highly unstable magma layer surrounded by a stable vapor reservoir. We turn to such a possibility in the next section.

Figure 4 presents the disk profile for x_c = 0.1,  T_c = 3898 K. The temperature again deceases with height at a fairly constant rate, but the photosphere is higher. Interestingly, the gas mass fraction at first increases, but reaches a maximum of ∼13% near z/H_* ∼ 1.35. Pressure and densities still drop quickly, but with an e-folding distance about three times that of Figure 3 consistent with the ∼$\sqrt {x_c }$ dependence of the effective scale height. As a result of the increased disk thickness, the gas density and total density, although still less than and greater than ρ_* respectively, are closer to that value. In Figure 4(d), we see that the sound speed is no longer monolithic with height. Although it increases during most of the cumulative surface density (shown in Figure 4(f)), it decreases in the higher altitude, more rarefied portion of the structure.

In Figure 5, the mid-plane gas mass fraction is assumed to be unity and the temperature T_c = 4218 K. For high x the sound speed is approximately

$$\begin{eqnarray} c &\approx & \left({\frac{{xRT}}{\mu }} \right)^{1/2} \approx \left({\frac{{x_c RT_c }}{\mu }} \right)^{1/2} \left[ {1 - \frac{1}{{x_c }}\left({\frac{z}{{H_o }}} \right)^2 } \right]^{1/2} \nonumber\\ &\approx & (x_c \ell)^{1/2} \left({\frac{{T_c }}{{T_o }}} \right)^{1/2} \left[ {1 - \frac{1}{{2x_c }}\left({\frac{z}{{H_o }}} \right)^2 } \right]. \end{eqnarray} \tag{ 37 }$$

The mid-plane sound speed is c_c = 1.09 × 10⁵cms⁻¹, which is far too high for gravitational instability. However, now both the gas mass fraction (Figure 5(b)) and sound speed (Figure 5(e)) monotonically decrease with altitude. The scale height H_p ∼ H_c = 3.49 × 10⁸ cm, which is over an order of magnitude greater than the low x disk considered in Section 3.1. The corresponding central pressure is P_c = 1.99 × 10⁸ dynes and the central gas density is ρ_{g, c} = 1.67 × 10⁻²gcm⁻³. Since x_c = 1, this is also the total density ρ_c. Compared to the isothermal reference model, T_c/T_* = 2.11, P_c/P_* = 1.53, ρ_{g, c}/ρ_* = ρ_c/ρ_* = 0.427,  H_p/H_* = 1.45. These are not too different from the isothermal case, although the two-phase disk is a little thicker. Note that in this model, most of the mass resides in the vapor state.

If the disk is stable but continues to radiate, the energy losses are compensated by latent heat release through condensation. Assuming that the droplets can settle, a magma layer will begin to collect in the mid-plane. The magma layer can itself be susceptible to gravitational instability. This is not just because the sound speed is low if the gas mass fraction is very small in the layer, but also because a typical parcel of the magma can be large enough that it is not well coupled to the gas through drag forces. The layer is not metastable but undergoes vigorous dissipation with an effective viscosity given by the original Ward–Cameron expression, Equation (2). The dissipation obtained from Equation (3) is at first insufficient to satisfy radiation losses, but as the layer density builds up, the dissipation rate climbs. The critical surface density at which these rates are comparable is found be equating dE/dt = F to find

$$\begin{equation} \sigma _{{\rm crit}} \approx \big(\sigma _{{\rm SB}} T_{{\rm ph}}^4 \Omega \big/\pi ^2 G^2\big)^{1/3}. \end{equation} \tag{ 38 }$$

At a = 2 R_⊕,  T_ph = 2000 K this reads σ_crit ≈ 2 × 10⁶ g cm⁻², which is about ∼1/5 of our assumed starting surface density value. The remaining ∼4/5 of the material resides in the stable atmosphere above the mid-plane. This configuration is represented schematically in Figure 2 by the labels A and M indicating the atmosphere and magma layer. Substituting Equation (38) into (2) and using this in the estimated spreading timescale yields, τ_ν ∼ 50 yr. As material spreads across the Roche boundary and into the Earth at the inner boundary, the surface density of the magma layer goes down and so does its energy dissipation. As it drops below the radiation rate, more material condenses and settles to replenish the magma layer. This cycle continues until the vapor atmosphere is exhausted in ∼250 yr. This constitutes an alternative way in which viscous and radiation time scales can be reconciled.

An additional complication arises for a stratified disk where there is an interface between the magma layer and atmosphere. This can introduce boundary layer drag that can decay the orbits of the magma layer on a timescale as short as a few years. This is due to loss of angular momentum to the vapor atmosphere which can be sub-Keplerian if an outward pressure gradient partially supports it. This issue was also considered in the treatment of Machida & Abe (2004) assuming a constant ratio of the radial pressure force to Earth's gravity. However, since there is no external torque, the entire disk cannot destroy itself. Material lost to the Earth has a lower specific angular momentum than the rest of the disk so that the specific angular momentum of the latter must increase. To accommodate this, the disk restructures itself with an outward migration of vapor that tries to remove the pressure gradient and shut off the boundary layer drag. To illustrate, consider a power-law disk, σ∝r⁻ⁿ, of mass M_n that extends from r = R_⊕ to 3R_⊕. The disk angular momentum is $L_n = 2\pi \smallint _{{\rm disk}} \sigma (GM_ \oplus r)^{1/2} rdr =$ M_nj_⊕(2 − n)(3^{5/2 − n} − 1)/(5/2 − 1)(3^{2 − n} − 1), where j_⊕ ≡ (GM_⊕R_⊕)^1/2 denotes the specific angular momentum if orbiting at Earth's radius. Suppose the initial disk has a 1/r dependence for which L₁ = 1.40M₁j_⊕. If the disk restructures itself to a constant surface density, its angular momentum will be L₀ = 1.46M₀j_⊕. During the process, the Earth absorbs an angular momentum ΔL = (M₁ − M₀)j_⊕ carried by the disk material that decays into it. Conservation of angular momentum requires that ΔL = L₁ − L₀ from which we conclude that M₀ = 0.87M₁, indicting that ∼13% of the disk is sacrificed to flatten the density gradient.