CIRCUMPLANETARY DISK FORMATION

The regular satellites of Jupiter and Saturn have orbits that are prograde, nearly circular, and approximately co-planar, suggesting that they formed within circumplanetary disks. Early studies of giant planet satellite formation employed a bottom-up or deconstruction approach in which a characterization of the satellite precursor disk was attempted by augmenting the satellite system mass to roughly solar composition (e.g., review by Pollack et al. 1991). This so-called minimum mass sub-nebula (MMSN) model is similar to that historically used in planetary formation models. However, over the intervening period our understanding of giant planet formation and disk–planet interactions has expanded considerably, and recently a new class of models has emerged utilizing a top-down approach that attempts to describe disk generation in the context of giant planet formation and contraction (Canup & Ward 2002, 2006, 2009; Magni & Coradini 2004; Alibert et al. 2005; Sasaki et al. 2010). This results in a dynamically evolving disk throughout most of the planet's contraction phase and a disk structure that can differ markedly from the earlier, more or less static, minimum mass models.

When a gas giant planet undergoes runaway gas accretion, it nearly fills the entire Hill sphere within which its gravity dominates over that of the Sun. The ratio of the Hill sphere radius, R_H, to the planet's radius, R, is ∼(4πρa³/9 M_⋆)^1/3, where M_⋆ is Sun's mass. For example, at Jupiter's distance, a = 5.2 AU, and current mean density, ρ = 1.3 g cm⁻³, the Hill sphere of a planet is more than 700 times the planet's radius; at Saturn's distance, this ratio rises to over a thousand. A forming gas planet will begin to contract within its Hill sphere once the rate of its gas accretion can no longer compensate for the increasing rate of its gravitational contraction due to the planet's growing mass and luminosity. As the planet contracts, two processes are thought to be effective means for producing circumplanetary disks of gas (H + He) and solids (rock + ice). (1) Accreted gas delivers angular momentum to the planet, and as the spinning planet contracts, conservation of angular momentum dictates that its rate of rotation increases. When this rate reaches the critical rate for rotational instability, the planet's outer equatorial layers begin to shed, forming a so-called "spin-out disk." (2) If the solar nebula is still present during the contraction phase, so that gas accretion continues, at some point the gas inflowing from solar orbit contains too much angular momentum to fall directly onto the planet as it becomes smaller and smaller. Instead, gas flows into orbit around the planet and directly creates an "accretion disk".

Prior disk formation models have typically assumed either that the accretion of gas ends before the planet contracts, so that as the planet contracts a pure spin-out disk is formed (e.g., Korycansky et al. 1991), or that the planet has fully contracted while it continues to accrete gas, resulting in a pure accretion disk (e.g., Canup & Ward 2002; D'Angelo et al. 2003). More likely, the disk produced around a growing gas giant may be due to a combination of both processes. Understanding how disks are created around gas giants is important not only for modeling a gas planet's growth, but also because the disk properties ultimately determine the nature of its final satellites, including their bulk composition, masses, orbital distributions, etc.

In the next section, the basic physical processes involved in disk formation and evolution are introduced, and the relatively simple case of a planet that has accreted all of its gas prior to its contraction is discussed. We estimate the properties of the resulting spin-out disk as a function of the planet's angular momentum, and find, in particular, that an inviscid spin-out disk is generally inconsistent with the properties of the regular satellite systems at Jupiter and Saturn. In Section 3, we proceed to the more complex situation when a planet is simultaneously contracting and accreting gas. We determine how inflowing gas is partitioned between the planet versus a circumplanetary disk as a function of the planet's contraction history, the mass inflow rate, and the angular momentum contained in the inflow. Section 4 then addresses the still open question of angular momentum accretion by the planet. In Section 5, the disk evolution model is merged with both inflow and planet evolution models to complete the solution. We consider several different descriptions for the planet's contraction: an existing numerical model (Papaloizou & Nelson 2005), a proxy model that describes the planet as a polytrope, allowing one to examine the effects of varying the adopted parameters of the system, and a model that includes a short duration collapse phase that could have occurred after the planet has contracted to a small fraction of its starting radius. We close the presentation in Section 6 with a discussion of the results and their ramifications. Section 7 then concludes the paper with a review of our overall findings and their implications. For the reader's convenience, a list of symbol definitions has been included in Table 1.

Table 1. Definition of Symbols

a	Planet semi-major axis
c	Gas sound speed
E, ER, ET	Polytrope total energy, rotational energy, thermal energy
$\mathscr{F}(t)$	Mass inflow rate
$\mathscr{F}_{{\rm shear}} {\rm }\mathscr{F}_{{\rm diff}}$	Shear- or torque-limited mass inflow rate, diffusion-limited mass inflow rate
finflow	Mass inflow rate per unit area
F(r)	In-plane mass flux, circumplanetary disk
Fd	Mass flux across outer disk edge
Fp	Mass flux from planet to disk
g(r), gd, gp	Viscous couple at distance r, at outer disk, at planet
h	Circumsolar disk half-thickness
Id	Moment of inertia, circumplanetary disk
J, JR	Cumulative angular momentum in inertial frame, in rotating frame
j(r), jd, jp	Specific angular momentum of circumplanetary disk at r, at disk outer edge, at planet
jR	Specific angular momentum of circumplanetary disk in rotating frame
jc	Average inflow specific angular momentum
Δj	Spread in inflow specific angular momentum
jrot	Specific angular momentum of inflow in frame rotating with planet
ℓ	Inflow angular momentum bias
L, Lc	Planet spin angular momentum, critical spin angular momentum for rotational stability
Ldisk	Circumplanetary disk angular momentum
LT	Total accreted angular momentum
$\mathscr{L}$	Luminosity
m	Molar molecular weight
M, Mp, MJ, MS	Planet mass, final planet mass, mass of Jupiter, mass of Saturn
M⋆ , M⊕	Solar mass, Earth mass
Md, MD	Circumstellar disk local and global masses (see Section 5.2)
Mdisk	Circumplanetary disk mass
Misol	Isolation mass
Mo	Planet mass at time of separation from Hill sphere
Mν	Viscous mass
Mtq	Transition mass between diffusion and torque limited regimes
MT	Total accreted mass
n	Polytropic index
P	Polytrope pressure
R(t), Rp	Planet radius, final planet radius
ℜ	Universal gas constant
RB	Planet Bondi radius
RH, RHP	Planet Hill radius, final planet Hill radius
Rrot	Planet radius at commencement of spin-out
r	Orbital radius, circumplanetary disk
rd	Outer radius of circumplanetary disk
rc, ri, ro	Centrifugal radius, inner radius of inflow region, outer radius of inflow region
rs, rz	Distance marker for stage (Equation (24)), distance marker for zone (Equation (25))
T	Polytrope temperature
Teff	Effective temperature of polytrope
Tc	Central temperature of polytrope
U	Polytrope potential energy
xo	Unperturbed streamline location at large distance from planet
xmax	Maximum streamline x-value that intercepts planet Hill sphere
xmin	Minimum streamline x-value that intercepts planet Hill sphere
α	Alpha viscosity parameter
γ	Outer radius of circumplanetary disk in Hill radii
γa	Adiabatic index
η	Gas enthalpy
κ	Opacity
λ	Planet moment of inertia coefficient
ν, νd	Viscosity, viscosity at circumplanetary disk outer edge
ρ	Planet density
Σ	Gas surface density, circumsolar disk
σ	Gas surface density, circumplanetary disk
σSB	Stephan–Bolztman constant
τneb	Nebular dissipation timescale
τν	Viscous timescale
τcol	Envelope collapse timescale
τKH, τKH,J	Planet Kelvin–Helmholtz timescale, Kelvin–Helmholtz timescale for MJ, RJ
φ	Planet gravitational potential
Ω	Planet orbital frequency
ω, ωc, ω'	Planet angular rotation rate, critical rotation rate, rotation rate scaled to critical rate

Download table as:  ASCIITypeset images: 1 2

2.1. Angular Momentum Bias

The planet of mass M acquires angular momentum L from the incoming material, and material with both positive and negative angular momentum will be accreted. The net magnitude is difficult to estimate without detailed hydrodynamic modeling, although some dimensional arguments have been employed to provide rough estimates. A natural characteristic length scale is the planet's Hill radius, R_H, while its mean motion Ω = (GM_⋆/a³)^1/2 provides a characteristic frequency. We define the angular momentum bias, ℓ, of the inflowing material by setting

$$\begin{equation} j_c \equiv {\it dL}/{\it dM} = \ell R_H^2 \Omega, \end{equation} \tag{ 1 }$$

where j_c is the average specific angular momentum of the inflow and prograde (retrograde) rotation corresponds to ℓ>0(<0). For example, if one considers a strictly two-dimensional (2D) Keplerian flow across the radius $\mathop R\nolimits_H$ surrounding the planet, then the angular momentum bias of the material crossing this boundary is ℓ = 1/4 (Ruskol 1982; Lissauer & Kary 1991; Mosqueira & Estrada 2003). In Section 4, we discuss other estimates for ℓ based on both analytic arguments and hydrodynamical simulation results. However, the bias is still a key uncertainty, and as such, we construct generalized models that treat ℓ as an unspecified parameter before returning to this issue in Section 4.

2.2. Rotational Instability

As the planet contracts, its spin frequency, ω ∼ L/C, increases, where C = λMR² is its moment of inertia, and we will assume for simplicity that convection maintains nearly uniform rotation within the planet. A moment of inertia parameter can be roughly estimated from a polytropic model (e.g., Chandrasekhar 1958). For example, for polytropic indices n = 1.5,  2,  2.5, the gyration constant is λ = 0.21, 0.15, 0.11. For comparison, Jupiter's current value is λ = 0.26.

When the planet's spin frequency reaches the critical value for rotational instability, ω_c, the planet starts to shed material from its equator into a disk. If it is assumed that the bias is independent of mass, the planet's cumulative angular momentum is found by integrating Equation (1),

$$\begin{equation} L = \int j_c {\it dM} = (3\ell /5) {\it MR}_H^2 \Omega. \end{equation} \tag{ 2 }$$

For a full mass Jupiter, this gives L = ℓ(5.41 × 10⁴⁷) g cm² s⁻¹ (i.e., likely much larger than the current angular momentum of Jupiter and its satellites augmented to solar composition, ∼8.6 × 10⁴⁵ g  cm² s⁻¹; e.g., Korycansky et al. 1991). An estimate of $\mathop \omega \nolimits_c \approx (GM/\mathop R\nolimits^3 \mathop)\nolimits^{1/2}$ is obtained by setting the spin rate equal to the Kepler frequency at the surface of a spherically symmetric planet, where both M(t) and R(t) are time dependent.¹ As the planet contracts, rotational instability sets in once L equals the critical value,

$$\begin{equation} L_c = \lambda {\it MR}^2 \omega _c = \lambda M\left({\it GMR} \right)^{1/2}, \end{equation} \tag{ 3 }$$

which determines the critical planet radius, R_rot, for commencement of spin-out from the planet,

$$\begin{equation} R_{{\rm rot}} = \frac{1}{3}\left({\frac{{3\ell }}{{5\lambda }}} \right)^2 R_H. \end{equation} \tag{ 4 }$$

Note, however, that if $\ell > 5\lambda /\sqrt 3,$ Equation (4) does not apply since R_rot cannot exceed R_H. In this situation, as long as the planet fills its Hill sphere and accretion is limited by its luminosity (see Section 5.1) we assume that the planet maintains itself at the brink of stability. Once the planet becomes luminous enough to separate from the Hill boundary and contract, it can immediately begin to shed a disk.

2.3. Contraction of an Accreted Planet

We first examine a case in which the bulk of the planet's mass has already been accreted before rotational instability sets in, so that there is no additional mass inflow to the planet as it contracts. Once R ⩽ R_rot, the planet begins to shed mass into a disk. The behavior of the disk material depends on its ability to transport angular momentum via stresses.

2.3.1. An Inviscid, Zero Stress Disk

If we consider the extreme of a zero stress disk, shed material remains in orbit where it first leaves the planet. Mass and angular momentum conservations require $M + \mathop M\nolimits_{{\rm disk}} = \mathop M\nolimits_T,\,\,L + \mathop L\nolimits_{{\rm disk}} = \mathop L\nolimits_T,$ where the subscripts disk and T refer to the disk and total values, respectively, and accordingly$\mathop {\dot M}\nolimits_{{\rm disk}} = - \dot M,\;\;\mathop {\dot L}\nolimits_{{\rm disk}} = - \dot L.$ Since material is spun out at the local mean motion, we have $\mathop {\dot L}\nolimits_{{\rm disk}} = \mathop {\dot M}\nolimits_{{\rm disk}} ({\it GMR}\mathop)\nolimits^{1/2}$ implying

$$\begin{equation} \dot L = \dot M({\it GMR}\mathop)\nolimits^{1/2} \end{equation} \tag{ 5 }$$

as well. However, requiring the planet to maintain L = L_c, and differentiating Equation (3) leads to

$$\begin{equation} \dot L = L_c \left({\frac{3}{2}\frac{{\dot M}}{M} + \frac{1}{2}\frac{{\dot R}}{R}} \right). \end{equation} \tag{ 6 }$$

Equating Equations (5) and (6) and rearranging yield

$$\begin{equation} \left({1 - \frac{{3\lambda }}{2}} \right) \frac{R}{M}\frac{{\it dM}}{{\it dR}} = \frac{\lambda }{2} \end{equation} \tag{ 7 }$$

which integrated gives M = M_T(R/R_rot)^q as the mass of the contracting planet, where q = λ/(2 − 3λ). Consequently, the mass left behind in a disk as a function of the planet's radius is

$$\begin{equation} \mathop M\nolimits_{{\rm disk}} = \mathop M\nolimits_T [1 - (R/\mathop R\nolimits_{{\rm rot}} \mathop)\nolimits^{ q} ] \end{equation} \tag{ 8 }$$

with the final disk mass obtained by setting R = R_p ≪ R_rot, with R_p being the final radius of the planet. For example, with λ = 0.15, q = 0.097 is small, which mitigates its influence somewhat. If we set $\mathop R\nolimits_{{\rm rot}} \approx \mathop {10}\nolimits^2 \mathop R\nolimits_p,$ then the final mass of the planet is, $\mathop M\nolimits_p \sim 0.64\, M_T,$ while $\mathop M\nolimits_{{\rm disk}} = 0.36M_T.$ In the case of a Jupiter mass planet, M_p = M_J implies $\mathop M\nolimits_T = 1.56\,M_J,$ with the mass of the disk being about ∼56% of the mass of Jupiter. Even for $\mathop R\nolimits_{{\rm rot}} \approx 25R_p$ (implying a very small bias ℓ ≈ 0.08), the disk mass decreases somewhat but is still ∼0.37M_J. These masses are vastly larger than that necessary to account for observed satellite systems, i.e., ∼fM_sat, where f is the gas to solid ratio of the infalling material and M_sat is the total mass of the satellite system. For the giant planets in our solar system, M_sat ∼ 10⁻⁴M_p, and for a solar composition infall, f ∼ 10², which implies an MMSN of order ∼1% of the planet's mass. In addition, massive disks with M_disk/M ≳ O(10⁻¹) would likely be gravitationally unstable, which would generate an effective viscosity and invalidate the starting assumption of a zero stress disk.

The massive spin-out disk calculated above is a result of the large angular momentum excess accreted by the planet compared to the giant planet plus satellite systems observed in our solar system. Since the fully contracted planet is limited in the amount of angular momentum it can contain, a significant amount of material must remain in orbit to store it. Indeed, the initially spun out material has an even larger specific angular momentum than the inflow average because $(\mathop {{\it GMR}}\nolimits_{{\rm rot}} \mathop)\nolimits^{1/2} = (3\ell /5\lambda)\mathop R\nolimits_H^2 \Omega = 3\mathop j\nolimits_c /5\lambda > \mathop j\nolimits_c,$ assuming a largely convective and uniformly rotating planet. Thus given the estimated angular momentum budget associated with accreting a gas giant, the angular momentum that would need to be contained in an inviscid spin-out disk is so large that the resulting disk is too massive and/or too radially extended compared to the MMSN construct.

2.3.2. Viscous Disks

A number of factors suggest that once gas is in circumplanetary orbit, it will not remain static but will instead radially spread due to angular momentum transport. A viscous protosatellite disk has been invoked by Coradini et al. (1989), Canup & Ward (2002, 2006), Alibert et al. (2005), and Sasaki et al. (2010), and perhaps the most compelling argument for this is that viscous diffusion allows most of the disk's mass to evolve inward and be accreted by the planet, while most of the disk's angular momentum (contained in a small fraction of its mass) is transported outward where it can eventually be removed from circumplanetary orbit. Thus, a much smaller disk mass can be consistent with the expected angular momentum budget than in the inviscid case. A circumplanetary disk also shares basic traits with the circumsolar disk, in which kinematic viscosity is commonly invoked to explain, e.g., observed mass accretion rates onto stars in extrasolar systems (e.g., Stone et al. 2000).

The source and magnitude of viscosity in cold disks are actively debated, and while a number of candidate mechanisms have been proposed, the mode and rate of angular momentum transport are quite uncertain. Given this uncertainty, a simplified parameterization is often used, i.e., the so-called alpha model (Shakura & Sunyaev 1973) that defines a viscosity ν = αch, where α is a constant, c is the gas sound speed in the disk, and h ≈ c/Ω is its vertical scale height (i.e., the disk's half thickness). The alpha model envisions angular momentum transport due to communication between turbulent disk eddies with velocity, $\mathop \alpha \nolimits^{1/2} c,$ and eddy length scale, α^1/2h (e.g., Dubrulle et al. 1995). Because the characteristic eddy velocity and length scale would be no larger than the disk sound speed and the disk thickness, respectively, the α parameter is a constant less than unity. Circumstellar disk models typically consider $\mathop {10}\nolimits^{ - 4} < \alpha < 0.1$ (e.g., Stone et al. 2000).

It is generally accepted that magnetorotational instability (MRI) in weakly magnetized disks drives turbulence and outward angular momentum transport, with effective α-values ∼ 10⁻²–0.1 (Balbus & Hawley 1991). To be operative, MRI requires a minimum ionization fraction, which can be produced thermally (for temperatures in excess of ∼10³K), by galactic cosmic-rays (for regions of low gas surface density, e.g., Dullemond et al. 2007), or by stellar X-rays (e.g., Glassgold et al. 1997). However, the presence of dust grains can deactivate MRI because they are effective charge absorbers and can cause the ionization fraction to plummet (Sano et al. 2000). It is thus plausible that both the protosatellite disk and the circumsolar disk at Jupiter's distance would have been "layered": comprised of MRI "active" regions at the disk surfaces, together with a neutral, dust-rich "dead zone" surrounding the disk mid-plane (e.g., Gammie 1996). This would result in a smaller effective α as a vertically averaged value than in a fully MRI active disk.

Non-magnetic hydrodynamic mechanisms have also been proposed as sources of turbulence and viscous transport in neutral disks, including notably those that rely on non-axisymmetric effects, including decaying vortices (e.g., Lithwick 2007), baroclinic instability (Klahr & Bodenheimer 2003), and shear instability (Dubrulle et al. 2005). Recent experimental results claim to rule out angular momentum transport via axisymmetric hydrodynamical turbulence at levels corresponding to α>10⁻⁵ (Hantao et al. 2006). Gravitational torques from objects accreting within the disk may provide an effective viscosity (Goodman & Rafikov 2001; see also Canup & Ward 2002). In the specific context of an actively supplied accretion disk, the shock front resulting from the difference in the free-fall velocity of mass infalling to the disk and that of the orbiting gas is another potential source of turbulence and therefore viscosity (Cassen & Moosman 1981). Such an inflow-driven viscosity would cease once the inflow had stopped completely, but it is possible that the last remaining orbiting gas could be then removed by other sources of viscosity, e.g., by MRI once small grains have been removed through accretion onto the satellites. Spiral shocks generated by the inflowing gas may also drive angular momentum transport (Lubow et al. 1999), although such shocks appear weaker in 3D simulations than in 2D simulations (Bate et al. 2003).

Viscous angular momentum transport is characterized by a torque due to viscous shear, g = 3πσνj,exerted by the inner disk on the outer disk across the circumference 2πr, where σ is the gas surface density, and j = (GMr)^1/2 is the specific angular momentum at orbital radius r; the quantity g is often referred to as the viscous couple (Lynden-Bell & Pringle 1974). The viscous couple causes the angular momentum of an annular ring of mass δm = 2πrσδr to change at a rate

$$\begin{equation} \delta m\frac{{dj}}{{dt}} = g\left(r \right) - g\left({r + \delta r} \right) \to - \frac{{\partial g}}{{\partial r}}\delta r. \end{equation} \tag{ 9 }$$

With dj/dt = (∂j/∂r)u, where u = dr/dt is the radial velocity of disk material with u > 0 defined as an outward velocity, Equation (9) gives

$$\begin{equation} 2\pi \sigma r\delta r\frac{{\partial j}}{{\partial r}}u = - \frac{{\partial g}}{{\partial r}}\delta r. \end{equation} \tag{ 10 }$$

Thus, the couple drives an in-plane radial mass flux in the disk, F = 2πrσu,that satisfies (Lynden-Bell & Pringle 1974)

$$\begin{equation} F\frac{{\partial j}}{{\partial r}} = - \frac{{\partial g}}{{\partial r}}. \end{equation} \tag{ 11 }$$

When the outer edge of the disk expands to distances comparable to the planet's Hill radius, at some point material escapes from the gravity of the planet and is returned to the circumstellar orbit. Accordingly a viscous circumplanetary disk will have an outer edge r_d = γR_H, where γ is a moderate fraction of unity, and a vanishing couple there,² with g(r_d = 0). Provided the viscous spreading time, τ_ν ∼ r²_d/ν, is less than the characteristic time scale over which processes supplying material to the disk change, a quasi-steady state is attained wherein the time variation in the flux F can be ignored when integrating Equation (11) so that F = F(r). When there is no inflow onto the disk, a steady-state F is constant throughout and integration of Equation (11) gives Fj + g = C. Evaluating C at r_d gives

$$\begin{equation} g\left(r \right) = F\left({j_d - j} \right) = - \dot M[ {\left({GMr_d } \right)^{1/2} - \left({GMr} \right)^{1/2} } ], \end{equation} \tag{ 12 }$$

where j_d = j(r_d).

To track the mass of the planet, we again employ conservation of angular momentum, but Equation (5) is now incomplete. The angular momentum of the planet not only decreases due to the shedding of material, but also due to the viscous torque, −g_p = −g(R) exerted on the planet by the disk, i.e.,

$$\begin{equation} \dot L = \dot M\left({{\it GMR}} \right)^{1/2} - g_p = \dot M\left({GMr_d } \right)^{1/2}. \end{equation} \tag{ 13 }$$

Equating this to Equation (6) yields

$$\begin{equation} \left({\left({\frac{{r_d }}{R}} \right)^{1/2} - \frac{{3\lambda }}{2}} \right)\frac{R}{M}\frac{{dM}}{{dR}} = \frac{\lambda }{2} \end{equation} \tag{ 14 }$$

where the quantity (r_d/R)^1/2 = (γR_H/R)^1/2 now replaces unity in the lead bracket. When $r_d > R_{{\rm rot}} \,\xyz\, R,$ the term 3λ/2 can safely be ignored and the equation integrated to give

$$\begin{equation} M/M_T \approx \left[ {1 - \frac{\lambda }{6}\frac{{R_{{\rm rot}}^{1/2} - R^{1/2} }}{{r_{dT}^{1/2} }}} \right]^6, \end{equation} \tag{ 15 }$$

where r_dT = r_d(M_T). Expanding to first order, the mass, ΔM, shed by the planet reads

$$\begin{equation} \frac{{\Delta M}}{{M_T }} \approx \lambda \left({\frac{{R_{{\rm rot}} }}{{r_{dT} }}} \right)^{1/2} \left[ {1 - \left({\frac{R}{{R_{{\rm rot}} }}} \right)^{1/2} } \right]. \end{equation} \tag{ 16 }$$

For R = R_p, R_HP = 744R_p, r_dT = γR_HP, λ = 0.15, the percentage shed is ∼(4.9, 2.2)/γ^1/2 for R_rot = (10², 25)R_p, respectively, which is much less than in the inviscid case.

Although Equation (16) gives us the mass processed through the disk, unlike the inviscid case, not all of the gas phase resides in the disk at the same time. Differentiating ΔM_g = fΔM/(1 + f) ≈ ΔM for f ≫ 1 with respect to time tells us the flux through the disk,

$$\begin{equation} F \approx - \frac{\lambda }{2}\left({\frac{R}{{r_d }}} \right)^{1/2} \frac{{\dot R}}{R}M_T, \end{equation} \tag{ 17 }$$

which is positive since $\dot R < 0.$ From Equation (12) and the definition of g(r), the surface gas density can be found

$$\begin{equation} \sigma = \frac{{\lambda M_T }}{{6\pi \nu}}\left({\frac{R}{{r_{dT} }}} \right)^{1/2} \left| {\frac{{\dot R}}{R}} \right|\left[ {\left({\frac{{r_d }}{r}} \right)^{1/2} - 1} \right] \end{equation} \tag{ 18 }$$

and integration over the disk yields its mass at any given time,

$$\begin{eqnarray} M_{{\rm disk}} &=& 2\pi \int_R^{r_d } {\sigma rdr} = \frac{\lambda }{3}M_T\left({\frac{R}{{r_{dT} }}} \right)^{1/2} \left| {\frac{{\dot R}}{R}} \right|\nonumber\\ &&\times \int_R^{r_d } {\left[ {\left({\frac{{r_d }}{r}} \right)^{1/2} - 1} \right]} \frac{r}{\nu}dr. \end{eqnarray} \tag{ 19 }$$

For instance, if ν = ν_d(r/r_d),

$$\begin{equation} M_{{\rm disk}} = \frac{\lambda }{3}M_T \left({\frac{{r_d^2 }}{{\nu_d }}} \right)\left| {\frac{{\dot R}}{R}} \right|\left({\frac{R}{{r_{dT} }}} \right)^{1/2} \left[ {1 - \left({\frac{R}{{r_d }}} \right)^{1/2} } \right]^2. \end{equation} \tag{ 20 }$$

Assuming R ≪ r_d, the ratio of this to ΔM is of order M_disk/ΔM ≈ (1/3)(R/R_rot)^1/2(τ_ν/τ_KH), where we have introduced the Kelvin–Helmholtz contraction timescale for the planet $\tau _{{\rm KH}} \equiv | {R/\dot R} |$. Since it was assumed at the outset that the disk's viscous time scale is much shorter than the planet's contraction time scale, only a small fraction of ΔM_g is actually in the disk at a given time. The gas that has left the system and gone into circumstellar orbit is of order ΔM_g − M_disk. For R_HP = 744R_p, λ = 0.15, M_disk/M_T ≈ 1.8 × 10⁻³γ^−1/2(τ_ν/τ_KH)(R/R_p)^1/2. On the other hand, solid particles do not necessarily follow this evolution if they coagulate and decouple from the gas phase.

So far we have considered a planet that had acquired its full mass prior to its contraction. However, estimates of maximum gas accretion rates, (∼10⁻² M_⊕ yr⁻¹, e.g., Hayashi et al. 1985) imply that, e.g., Jupiter would begin to contract to a size smaller than its Hill sphere once it exceeded ∼0.2–0.3 M_J (e.g., Tajima & Nakagawa 1997). Papaloizou & Nelson (2005) tracked a Jovian planet's evolution including ongoing gas accretion and predicted that a Jupiter-mass planet would contract to within the current orbit of Io in less than a million years. When accretion is ongoing during planet contraction, the system passes through stages in which the circumplanetary disk transitions from a spin-out disk to an accretion disk as we next describe.

3.1. Inflow Pattern

A quantity closely related to the angular momentum bias, ℓ, is the so-called centrifugal radius, r_c, of the inflowing material (e.g., Cassen & Summers 1983). This is the distance from the planet where a circular orbit would have the same specific angular momentum as the inflow average j_c, i.e.,

$$\begin{equation} r_c = \frac{1}{GM}\left(\frac{dL}{dM} \right)^2 = \frac{\ell ^2 }{3}R_H. \end{equation} \tag{ 21 }$$

If r_c exceeds the planet radius R, most of the inflowing material would have too much angular momentum to reach the planet and would instead accumulate in a disk. Setting ℓ ∼ 1/4 (Section 2.1) results in r_c ∼ R_H/48, which for Jupiter (R_HP ≈ 744R_J) and Saturn (R_HP = 1080R_S) implies r_c ≈ 15R_J and 23R_S, respectively, as the planets approach their final masses. This has led to the suggestion that the compactness of these planet's satellite systems compared to their Hill spheres is due in part to a small prograde angular momentum bias of incoming material (e.g., Stevenson et al. 1986; Canup & Ward 2002; Mosqueira & Estrada 2003).

Of course, it is unlikely that at any given time, all inflowing material has a specific angular momentum, j, equal to the average dL/dM = j_c. Both positive and negative angular momentum material would mix as it descends toward the planet. Numerical simulations have not yet determined this angular momentum distribution. In the absence of this information, we adopt a heuristic model in which an incoming mass δM is distributed evenly across a range of specific angular momenta j_c − Δj/2 < j < j_c + Δj/2 so that there is a mass dM = (δM/Δj)dj in any interval dj. When this material achieves centrifugal balance, dj = (GMr)^1/2dr/2r. For a mass inflow rate $\mathscr{F}$, this implies an inflow rate per unit area, f_inflow, to the planet's and/or disk's mid-plane of

$$\begin{equation} f_{{\rm inflow}} = \frac{\mathscr{F}}{{4\pi r_c^2 }}\left({\frac{{j_c}}{{\Delta j}}} \right)\left({\frac{{r_c }}{r}} \right)^{3/2} \end{equation} \tag{ 22 }$$

with inner and outer boundaries of r_i(o) = r_c(1 ∓ Δj/2j_c)². Assuming a broad distribution with Δj ≈ j_c, the inner and outer boundaries of the inflow upon reaching the planet or disk are r_i ≈ r_c/4, r_o ≈ 9r_c/4 (the latter of which is comparable to the outer limit assumed by Canup & Ward 2002).

3.2. Inflow Stages

As the planet contracts and the centrifugal radius increases due to its growing mass and Hill radius, the disk–planet configuration passes through three stages depending on whether the inflow falls on the planet or the disk or both. These are depicted schematically in Figure 1. The overall behavior hinges on which state is occupied as well as whether the planet rotates at a stable rate or at the critical rate, ω_c. In the former case, the viscous couple vanishes at the planet, g_p = 0, and the planet can accept material from the disk unimpeded while increasing its rotation rate. These stages will be prescripted with an "S'' to signify stable rotation. In stage S1, there is no disk; in stages S2 and S3, a disk forms by direct infall, and a stagnation point (F = 0) exists within the inflow boundaries (zone 2) separating inward and outward fluxes of disk material, with the flux at the planet always inward (F_p < 0). In contrast, if the planet rotates at the critical rate for stability, ω_c, the viscous couple at the planet–disk interface adjusts to the value needed to maintain this state. These stages will be prescripted with a "C'' to indicate critical rotation. This results in material shedding (F_p > 0) by the planet, and the development of a spin-out disk in stage C1. In stages C2 and C3, the disk is due to a combination of direct infall and the planet's tendency to shed, resulting in a value of F_p that can be either positive or negative, depending on which influence dominates. In this section, we solve for these cases in detail.

Zoom In Zoom Out Reset image size

3.3. In-plane Flux, F(r)

From Equation (22) with Δj = j_c we have an inflow per unit area of $f_{{\rm inflow}} = \mathscr{F}/({4\pi r_c^{1/2} r^{3/2} })$ between r_i = r_c/4 and r_o = 9r_c/4. This causes a divergence of the in-plane flux, F(r), such that (e.g., Canup & Ward 2002)

$$\begin{equation} \frac{{dF}}{{dr}} = 2\pi rf_{{\rm inflow}} = \frac{\mathscr{F}}{{2r_c^{1/2} r^{1/2} }}. \end{equation} \tag{ 23 }$$

Equation (23) can be integrated to obtain the flux variation with r for the various stages. In the following expressions, the flux F_n is subscripted by a zone number. Zone 1 is the region of the disk exterior to r_o, zone 2 is the disk region between r_i and r_o, while zone 3 is the disk region interior to the inflow pattern, r < r_i, as shown in Figure 1.

• 1.  
  Stage 1. Since R > r_o, there is no infall on the disk and only zone 1 exists for which the flux is constant throughout: F₁ = F_p = F_d, where F_d = F(r_d).
• 2.  
  Stage 2. A zone 2 appears from the planet's radius to r_o within which the flux increases (becomes more positive), while outside of that in zone 1, the flux is again constant:

  $$\begin{eqnarray*} && F_2 = F_p + \mathscr{F}\left[ {\left({\frac{r}{{r_c }}} \right)^{1/2} - \left({\frac{R}{{r_c }}} \right)^{1/2} } \right],\quad \nonumber\\ && F_1 = F_p + \mathscr{F} \left[ {\frac{3}{2} - \left({\frac{R}{{r_c }}} \right)^{1/2} } \right] = F_d \end{eqnarray*}$$
• 3.  
  Stage 3. There is now a zone 3 close to the planet where the flux is constant, then a flux that becomes more positive between r_i and r_o, followed by another region of constant flux:

  $$\begin{eqnarray*} F_3 &=& F_p, \nonumber\\ F_2 &=& F_p + \mathscr{F}\left[ {\left({\frac{r}{{r_c }}} \right)^{1/2} - \frac{1}{2}} \right],\quad \nonumber\\ F_1 &=& F_p + \mathscr{F} = F_d. \end{eqnarray*}$$

To economize on notation, we introduce two distance markers: the first changes with the stage,

$$\begin{eqnarray} r_s \equiv r_o \ &{\rm if }&\ r_o < R{\rm }\ ({\rm Stage\ 1}); \nonumber\\ {\rm }R\ &{\rm if }&\ r_i < R < r_o {\rm }\ ({\rm Stage\ 2});\quad \nonumber\\ {\rm }r_i\ &{\rm if }& \ R < r_i {\rm }\ ({\rm Stage\ 3}) \end{eqnarray} \tag{ 24 }$$

and the second changes with the zone,

$$\begin{eqnarray} r_z \equiv r_o\ &{\rm if }&\ r_o < r{\rm }\ ({\rm Zone\ 1}),\nonumber\\ {\rm }r\ &{\rm if }&\ r_i < r < r_o {\rm }\ ({{\rm Zone\ 2}}),\quad \nonumber\\ {\rm }r_i \ &{\rm if }& \ r < r_i {\rm }\left({{\rm Zone\ 3}} \right). \end{eqnarray} \tag{ 25 }$$

Using these we can write the flux in the generalized form

$$\begin{equation} F\left(r \right) = F_p + \mathscr{F}\left[ {\left({\frac{{r_z }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right] \end{equation} \tag{ 26 }$$

Equations (24)–(26) together with the obvious condition r ⩾ R can reproduce the flux for each case. In each expression, the flux is related to its value F_p at the planet–disk boundary. This flux can be positive with the planet shedding material to the disk, or negative with the planet accreting material from the disk. Evaluating Equation (26) at r_d tells us that the flux exiting the disk at the outer edge is

$$\begin{equation} F_d = F_p + \mathscr{F}\left[ {\left({\frac{{r_o }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right]. \end{equation} \tag{ 27 }$$

The value of F_d is always positive. With this, a useful alternative expression for F in terms of F_d can be found, namely,

$$\begin{equation} F\left(r \right) = F_d - \mathscr{F}\left[ {\left({\frac{{r_o }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_z }}{{r_c }}} \right)^{1/2} } \right]. \end{equation} \tag{ 28 }$$

3.4. Viscous Couple

To determine F_p (or F_d) we must also know the variation in the viscous couple, g(r). The couple is found by integrating Equation (11). Recall that we assume quasi-equilibrium states exist so that ∂σ/∂t = 0. In zones 1 and 3 where the flux is constant, the combination Fj + g will be constant as well. However, in zone 2, we use Fdj/dr = d(Fj)/dr − jdF/dr to rewrite Equation (11) as

$$\begin{equation} \frac{d}{{dr}}\left({Fj + g} \right) = j\frac{{dF}}{{dr}} \end{equation} \tag{ 29 }$$

and integrate with the help of Equation (23). The expressions below subscript the couple by the zone number and relate the flux at the planet to the value of the couple there.

• 1.  
  Stage 1: only zone 1 pertains for which F₁j + g = F_pj_p + g_p = F_dj_d, where in the final term the couple is assumed to vanish at the outer disk edge, g(r_d) = 0. We conclude that

  $$\begin{eqnarray*} && F_p ({j_d - j_p }) = g_p, \nonumber\\ && g_1 = g_p + F_p ({j_p - j}). \end{eqnarray*}$$
• 2.  
  Stage 2: in zone 2, Equation (29) is integrated from the planet's radius and Equation (26) used to eliminate F to yield,

  $$\begin{eqnarray*} g_2 &=& g_p + F_p \left({j_p - j} \right) - \mathscr{F}\left[ {\left({\frac{r}{{r_c }}} \right)^{1/2} - \left({\frac{R}{{r_c }}} \right)^{1/2} } \right]j \nonumber\\ &&+ \frac{1}{2}\mathscr{F}\left[ {\frac{r}{{r_c }} - \frac{R}{{r_c }}} \right]j_c, \end{eqnarray*}$$

  while in zone 1 we again write, F₁j + g = F_dj_d, but use F₁ from stage 2 to find,

  $$\begin{equation*} g_1 = \left[ {F_p + \mathscr{F}\left[ {\frac{3}{2} - \left({\frac{R}{{r_c }}} \right)^{1/2} } \right]} \right]\left({j_d - j} \right). \end{equation*}$$

  Requiring that g₁(r_o) = g₂(r_o) at the mutual boundary of the two zones then reveals that

  $$\begin{eqnarray*} F_p ({j_d - j_p }) &=& g_p - \mathscr{F}\left[ {\frac{3}{2} - \left({\frac{R}{{r_c }}} \right)^{1/2} } \right]j_d \nonumber\\ &&+ \frac{1}{2}\mathscr{F}\left[ {\frac{9}{4} - \frac{R}{{r_c }}} \right]j_c, \end{eqnarray*}$$

  $$\begin{eqnarray*} g_1 &=& g_p + F_p ({j_p - j}) - \mathscr{F}\left[ {\frac{3}{2} - \left({\frac{R}{{r_c }}} \right)^{1/2} } \right]j\nonumber\\ && + \frac{1}{2}\mathscr{F}\left[ {\frac{9}{4} - \frac{R}{{r_c }}} \right]j_c. \end{eqnarray*}$$
• 3.  
  Stage 3: following similar procedures one obtains

  $$\begin{eqnarray*} && F_p ({j_d - j_p }) = g_p - \mathscr{F}\left({j_d - j_c } \right), \nonumber \\ && g_3 = g_p + F_p \left({j_p - j} \right), \end{eqnarray*}$$

  $$\begin{eqnarray*} g_2 &=& g_p + F_p ({j_p - j}) - \mathscr{F}\left[ {\left({\frac{r}{r_c}} \right)^{1/2} - \frac{1}{2}} \right]j \nonumber\\ &&+ \frac{\mathscr{F}}{2}\left[ {\frac{r}{{r_c }} - \frac{1}{4}} \right]j_c \end{eqnarray*}$$

  $$\begin{equation*} g_1 = g_p + F_p ({j_p - j}) + \mathscr{F}\left({j_c - j} \right) \end{equation*}$$

  for stage 3 behavior.

As with the flux, these results can be written in generalized forms, namely,

$$\begin{eqnarray} F_p ({j_d - j_p }) &=& g_p - \mathscr{F}\left[ {\left({\frac{{r_o }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right]j_d \nonumber\\ &&+ \frac{\mathscr{F}}{2}\left[ {\frac{{r_o }}{{r_c }} - \frac{{r_s }}{{r_c }}} \right]j_c, \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} g &=& g_p - F_p ({j - j_p }) - \mathscr{F}\left[ {\left({\frac{{r_z }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right]j \nonumber\\ &&+ \frac{\mathscr{F}}{2}\left[ {\frac{{r_z }}{{r_c }} - \frac{{r_s }}{{r_c }}} \right]j_c. \end{eqnarray} \tag{ 31 }$$

3.5. Application to an Accreting/Contracting Planet

With inflow from the disk, the planet's growth rate is

$$\begin{equation} \dot M = \mathscr{F}\left[ {\left({\frac{{r_s }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_i }}{{r_c }}} \right)^{1/2} } \right] - F_p, \end{equation} \tag{ 32 }$$

where the first term is the infall directly onto the planet and the second is that accreted from (or lost to) the disk. In a quasi-steady state, conservation of mass dictates that the rate of inflow to the system must equal the rate of change of the planet's mass plus any mass loss at the outer edge, i.e., $\dot M = \mathscr{F} - F_d.$

3.5.1. Stable Rotation

If the planet rotates at a subcritical rate, it can accrete material from the disk even though it arrives with the critical specific angular momentum (GMR)^1/2, because of our assumption that nearly uniform planet rotation is maintained through convective mixing. Accordingly, we set the viscous couple at the planet to g_p = 0. Again for stage S1, no disk exists, but for stages S2 and S3 there is an in-plane flux, with part of the material falling on the disk being accreted by the planet (F_p < 0) as it spreads inward and the remainder spreading outward and eventually exiting the disk at r_d.

Setting g_p to zero, and using Equation (30) to eliminate F_p in Equation (31) leads to

$$\begin{eqnarray} g &=& \mathscr{F}\left({\frac{{j - j_p }}{{j_d - j_p }}} \right)\Bigg(\left[ {\left({\frac{{r_o }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right]j_d\nonumber\\ && - \frac{1}{2}\left[ {\frac{{r_o }}{{r_c }} - \frac{{r_s }}{{r_c }}} \right]j_c \Bigg)- \mathscr{F}\Bigg(\left[ {\left({\frac{{r_z }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right]j\nonumber\\ && - \frac{1}{2}\left[ {\frac{{r_z }}{{r_c }} - \frac{{r_s }}{{r_c }}} \right]j_c \Bigg). \end{eqnarray} \tag{ 33 }$$

As an example, let us consider a late stage S3 configuration for which r_s = r_i and Equation (30) gives $F_p = - \mathscr{F}({r_d^{1/2} - r_c^{1/2} })/({r_d^{1/2} - R^{1/2} }).$ In the limit that $\left({r_d /R} \right)^{1/2} \xyz 1,$ g/j = 3πσν reduces to

$$\begin{equation} 3\pi \sigma \nu \,{=}\, \mathscr{F}\! \left[ \!{1 \,{-}\, \left({\frac{{r_c }}{{r_d }}} \right)^{1/2} } \right]\left[\! {1 \,{-}\, \left({\frac{R}{r}} \right)^{1/2} } \right]\quad {\rm for}\ R \,{<}\, r \,{<}\, r_c /4, \end{equation} \tag{ 34 }$$

$$\begin{eqnarray*} &=& \mathscr{F}\Bigg(\left[ {1 - \left({\frac{{r_c }}{{r_d }}} \right)^{1/2} } \right]\left[ {1 - \left({\frac{R}{r}} \right)^{1/2} } \right]\nonumber\\ && + \frac{1}{2}\left[ {1 - \left({\frac{r}{{r_c }}} \right)^{1/2} - \frac{1}{4}\left({\frac{{r_c }}{r}} \right)^{1/2} } \right] \Bigg)\quad {\rm for}\ r_c /4 < r < 9r_c /4, \end{eqnarray*}$$

$$\begin{eqnarray*} &=& \mathscr{F}\Bigg(\left[{1 - \left({\frac{{r_c }}{{r_d }}} \right)^{1/2} } \right]\left[ {1 - \left({\frac{R}{r}} \right)^{1/2} } \right]\nonumber\\ && - 1 + \left({\frac{{r_c }}{r}} \right)^{1/2} \Bigg)\quad {\rm for} \ 9r_c /4 < r. \end{eqnarray*}$$

This is the analog of the quasi steady-state disk profile used by Canup & Ward (2002). The differences are due to their use of a constant inflow, $f_{{\rm inflow}} = \mathscr{F}/\pi r_o^2$ between r = 0 and an outer boundary of the inflow at 2r_c.

The strategy is then to track the quasi steady-state evolution of the disk viscous couple g together with M(t), R(t), which are needed to determine the appropriate values for r_c(M), r_d(M), and then to employ σ = g/3πνj to determine the disk profile for a given viscosity profile. For stable rotation, M(t) is found by setting g_p = 0 in Equation (30) and substituting into Equation (32),

$$\begin{eqnarray} \frac{{\dot M}}{M} &=& \frac{\mathscr{F}}{M}\Bigg(1 + \left({\frac{{j_p }}{{j_d - j_p }}} \right)\left[ {\left({\frac{{r_o }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right]\nonumber\\ && - \frac{1}{2}\left({\frac{{j_c }}{{j_d - j_p }}} \right)\left[ {\frac{{r_o }}{{r_c }} - \frac{{r_s }}{{r_c }}} \right]\Bigg). \end{eqnarray} \tag{ 35 }$$

The evolution of the planet's radius, R(t), on the other hand, must be obtained from a planet contraction model (see Section 5).

3.5.2. Critical Rotation

Once the planet is rotating at the critical rate, ω_c ≈ (GM/R³)^1/2, a combination of mass shedding and viscous torques (g_p ≠ 0) will maintain its spin angular momentum at L = L_c = λM(GMR)^1/2. Notwithstanding, the rate of change of L is

$$\begin{equation} \dot L = \frac{\mathscr{F}}{2}\left[ {\frac{{r_s }}{{r_c }} - \frac{{r_i }}{{r_c }}} \right]j_c - F_p j_p - g_p, \end{equation} \tag{ 36 }$$

where the first term is the angular momentum falling directly on the planet and the last two terms describe the angular momentum exchange with the disk. If Equation (30) is used to eliminate g_P, Equation (36) becomes

$$\begin{equation} \dot L = - F_p j_d + \mathscr{F}j_c - \mathscr{F}\left[ {\left({\frac{{r_o }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right]j_d. \end{equation} \tag{ 37 }$$

This must be equated to Equation (6) to find

$$\begin{eqnarray} \lambda M\left({\frac{R}{{r_d }}} \right)^{1/2} \left[ {\frac{3}{2}\frac{{\dot M}}{M} + \frac{1}{2}\frac{{\dot R}}{R}} \right] &=& \mathscr{F}\Bigg[ \left({\frac{{r_c }}{{r_d }}} \right)^{1/2} - \left({\frac{{r_o }}{{r_c }}} \right)^{1/2}\nonumber\\ && + \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} \Bigg ] - F_p. \end{eqnarray} \tag{ 38 }$$

Using (32) to remove F_P yields after some manipulation,

$$\begin{equation} \left[ {3 - \frac{2}{\lambda }\left({\frac{{r_d }}{R}} \right)^{1/2} } \right]\frac{{\dot M}}{M} + \frac{{\dot R}}{R} = - \frac{2}{\lambda }\left({\frac{{r_d }}{R}} \right)^{1/2} \frac{\mathscr{F}}{M}\left[ {1 - \left({\frac{{r_c }}{{r_d }}} \right)^{1/2} } \right]. \end{equation} \tag{ 39 }$$

In the limit $(r_d /R)^{1/2} \, \xyz\, 1,$ this expression reduces to

$$\begin{equation} \dot M \approx - \frac{\lambda }{2}\left({\frac{R}{{r_d }}} \right)^{1/2} \frac{M}{{\tau _{KH} }} + \mathscr{F}\left[ {1 - \left({\frac{{r_c }}{{r_d }}} \right)^{1/2} } \right]. \end{equation} \tag{ 40 }$$

Note that the disk outflow is

$$\begin{equation} F_d = \mathscr{F} - \dot M \approx \frac{\lambda }{2}\left({\frac{R}{{r_d }}} \right)^{1/2} \frac{M}{{\tau _{KH} }} + \mathscr{F}\left({\frac{{r_c }}{{r_d }}} \right)^{1/2}. \end{equation} \tag{ 41 }$$

The flux across the disk's outer edge is always positive and contains a contribution generated by the contraction and one by the inflow. However, the flux at the planet is

$$\begin{equation} F_p = F_d - \mathscr{F}\left[ {\frac{3}{2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right], \end{equation} \tag{ 42 }$$

which can be either positive (a spin out disk) or negative (an accretion disk). In state C1, F_p = F_d > 0, and the disk is due to spin-out. However, F_p progressively decreases from F_d in state C2, until in stage C3, $F_p = F_d - \mathscr{F},$ and the planet accretes material from the disk (F_p < 0) if

$$\begin{equation} \frac{\lambda }{2}\left({\frac{R}{{r_d }}} \right)^{1/2} \frac{M}{{\tau _{{\rm KH}} }} < \mathscr{F}\left[ {1 - \left({\frac{{r_c }}{{r_d }}} \right)^{1/2} } \right]. \end{equation} \tag{ 43 }$$

We now return to Equations (30) and (31), but this time g_P is eliminated to find

$$\begin{eqnarray} g &=& F_p \left({j_d - j} \right) + \mathscr{F}\left[ {\left({\frac{{r_o }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right]j_d\nonumber\\ && - \mathscr{F}\left[ {\left({\frac{{r_z }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right]j -\frac{\mathscr{F}}{2} \left[ {\frac{{r_o }}{{r_c }} - \frac{{r_z }}{{r_c }}} \right]j_c.\qquad \end{eqnarray} \tag{ 44 }$$

Using Equation (42), setting r_o = (9/4)r_c, g = 3πσνj, and dividing throughout by j gives

$$\begin{eqnarray} 3\pi \sigma \nu &=& ({\mathscr{F} - \dot M})\left[ {\left({\frac{{r_d }}{r}} \right)^{1/2} - 1} \right] + \mathscr{F}\left[ {\frac{3}{2} - \left({\frac{{r_z }}{{r_c }}} \right)^{1/2} } \right] \nonumber\\ &&- \frac{\mathscr{F}}{2}\left[ {\frac{9}{4} - \frac{{r_z }}{{r_c }}} \right]\left({\frac{{r_c }}{r}} \right)^{1/2}, \end{eqnarray} \tag{ 45 }$$

which completes the solution for the critical rotation case.

Equations (33–34) and (44–45) are among the main contributions of this paper. In the following sections, the critical rotation case will be applied to some specific examples. Before doing so, however, we will discuss the angular momentum of the inflow in more detail.

Disk solutions Equations (33) and (44) are given in terms of the centrifugal radius, r_c. This quantity is quite sensitive to the details of the inflow, depending on the square of the angular momentum bias, ℓ = j_c/R²_HΩ, which itself comes from averaging a specific angular momentum, j, that could vary widely across the inflow boundaries. Clearly, the compactness of the giant planet satellite systems could most economically be accounted for by an inflow that penetrates deeply into the planet's potential well and deposits presatellite material there directly. Support for this notion is provided by the calculations of Lissauer and co-workers that found ℓ = 1/4 for undeflected 2D Keplerian flow across a region of radius r = R_H (Lissauer & Kary 1991; Lissauer 1995), and was adopted as a fiducial value in some recent Galilean satellite formation models, e.g., Canup & Ward (2002) and Mosqueira & Estrada (2003) (although the latter applies this only to the inner three satellites). Undeflected flow can be viewed as a limiting case of streamline behavior found from integrating the planar Hill equations of motion

$$\begin{equation} d^2 x/dt^2 - 3\Omega ^2 x - 2\Omega dy/dt = - \partial \left({\varphi + \eta } \right)/\partial x, \end{equation} \tag{ 46 }$$

$$\begin{equation} d^2 y/dt^2 + 2\Omega dx/dt = - \partial \left({\varphi + \eta } \right)/\partial y, \end{equation} \tag{ 47 }$$

where x, y are local Cartesian coordinates in a frame attached to the planet and rotating at its circumstellar orbital frequency Ω, φ = −GM/(x² + y²)^1/2 is the planet's gravitational potential (ignoring the indirect part) and η is the enthalpy of the gas. The enthalpy tends to counteract the potential so that their sum, φ + η, varies less than φ alone (see, e.g., Paardekooper & Papaloizou 2009). In this context, Lissauer's problem can be considered a limiting case where φ + η = constant, and the solution reduces to pure Keplerian shear, with $\dot x = 0,$ x = x_o, $\dot y = - 3\Omega x_o /2,$ where x_o is the unperturbed streamline location at a large distance from the planet, namely, the dashed lines in Figure 2(a). In general, the specific angular momentum of the flow with respect to the planet at any point is given by

$$\begin{equation} j\left({x,y} \right) = x\dot y - y\dot x + r^2 \Omega = j_{{\rm rot}} + r^2 \Omega. \end{equation} \tag{ 48 }$$

The first term, j_rot, is the z-component of r × v while the second term compensates for the fact that the cross product is being taken in a rotating frame (e.g., Lissauer & Kary 1991). The evaluation of r²Ω = R²_HΩ is the same for a streamline crossing anywhere on the r = R_H circle, but j_rot is not, and finding its average for the entire flow, $\bar j_{{\rm rot}},$ requires weighting each streamline contribution by its flux fraction, $d\bar j_{{\rm rot}} = j_{{\rm rot}} d\mathscr{F}/\mathscr{F},$ where $d\mathscr{F} = \Sigma \left| v \right|dl$ is the flux along a "tube' of width dl, and we have here denoted the circumstellar disk surface density by Σ. If a steady-state flow is assumed, the flux $d\mathscr{F}$ is constant and can be evaluated far from the planet where the flow is Keplerian, so that Σ → Σ_o, v → −(3/2)x_oΩ, dl → dx_o, and $d\mathscr{F} = \left({3/2} \right)\Sigma _o \Omega x_o dx_o.$ Integrating $d\mathscr{F}$ across the range x_min < x_o < x_max of streamlines that eventually hit the Hill sphere, one finds $\mathscr{F}=\left({3/4} \right)\Sigma _o \Omega ({x_{\max }^2 - x_{\min }^2 }),$ and thus the quantity,

$$\begin{equation} \frac{{ {d{\bar j}} _{{\rm rot}} }}{{dx_o }} = \frac{{2j_{{\rm rot}} x_o }}{{x_{\max }^2 - x_{\min }^2 }} \end{equation} \tag{ 49 }$$

reveals from where contributions to $\bar j_{{\rm rot}}$ originate. For undeflected shear flow, i.e., neglecting the planet's gravity, j_rot = −(3/2)x²_oΩ, x_min = 0, x_max = R_H, $d\bar j_{{\rm rot}} /dx_o = - 3x_o^3 \Omega /R_H^2$, and $\bar j_{{\rm rot}} = - \left({3/4} \right)R_H^2 \Omega.$ The left curves in Figure 2(b) display j_rot, $d\bar j_{{\rm rot}} /dx_o$ as functions of x_o for the undeflected case in which all streamlines have j_rot ⩽ 0. Transferring to the non-rotating frame then gives the average specific angular momentum, $j_c = \bar j_{{\rm rot}} + R_H^2 \Omega = \left({1/4} \right)R_H^2 \Omega,$ implying a bias of ℓ = 1/4 (Lissauer & Kary 1991).

Zoom In Zoom Out Reset image size

It is instructive to compare this situation with the opposite extreme, φ + η ≈ φ, where parcels respond only to the planet's gravity; we will refer to this as the cold disk case corresponding to zero pressure and η ∼ O(c²) → 0. Figure 2(a) shows a numerical integration of trajectories in this limit (solid curves) and their corresponding contributions, j_rot, $d\bar j_{{\rm rot}} /dx_o,$ are also shown by the right curves in Figure 2(b). The resulting average is $\bar j_{{\rm rot}} = - 0.55R_H^2 \Omega,$ implying a bias ℓ ∼ 0.45, and would lead to a more extended disk, i.e., r_c ∼ 50R_J for Jupiter. Two effects contribute to a less negative value of $\bar j_{{\rm rot}}$ than in the pure shear case: (1) some of the streamlines now have positive values of j_rot because they acquire a velocity component toward the y-axis and (2) the range of x_o intersecting the r = R_H circle has shifted outward, increasing the flux fraction associated with less negative j_rot (Figure 2(b), right curves).

Under what conditions is φ + η ≈ constant a more appropriate approximation than φ + η ≈ φ? We note that φ ≈ η ≈ c² at r = GM/c² ≡ R_B, which is the Bondi radius. For r ≫ R_B, the enthalpy can counteract variations in the potential fairly well because as the density rises from gravitational focusing, the pressure rises as well, creating a pressure gradient that opposes the planet's gravitation (e.g., Paardekooper & Papaloizou 2009). However, when $r \ll R_B,\varphi \,\xyz\, \eta$ and pressure is no longer sufficient for this. This suggests that we might obtain a sense of the combined influence by a simple heuristic extension of Lissauer's approach that treats streamlines as undeflected until they come within the planet's Bondi radius and then follows the trajectories under the planet's gravity alone. Figure 3 displays as black diamonds the resulting ℓ's from such a procedure for selected values of R_B/R_H = 3(R_H/h)², where h ∼ c/Ω is the half-thickness of the circumstellar disk. For a constant Hill radius, there is a steady decrease of the bias with decreasing Bondi length or equivalently for a hotter, thicker nebula. What if R_B < R_H? In this circumstance, materials with impact parameters R_B < x_o < R_H also pass by the planet and R_B supercedes R_H as the accretion boundary. Consequently, j_c ∼ (1/4)R²_BΩ with an effective bias of ℓ ∼ (1/4)(R_B/R_H)²; this is indicated in Figure 3 by the black dotted curve.

Zoom In Zoom Out Reset image size

Although illustrative, its doubtful that the heuristic model above can provide accurate enough predictions to model the pre-satellite disk, and ultimately, numerical simulations will be required to ascertain the behavior of ℓ and r_c as functions of the planet's mass and surrounding circumstellar disk environment. Korycansky & Papaloizou (1996) numerically integrated the equations of motion along with the continuity equation, ∇ · (σv) = 0, for a 2D steady-state flow, assuming an isothermal disk and no accretion. They found that the flow pattern included a largely pressure supported region circulating the planet that has a radial extent almost equal to R_H. Since the material could not cool below the assumed isothermal temperature, this most closely corresponds to the early phase when the gas envelope fills the Hill sphere and cannot contract fast enough to detach from that boundary (i.e., designated Model A by Papaloizou & Nelson 2005, see Section 5). In that situation, the planet's accretion rate is regulated by its luminosity. Even so, the mass and angular momentum of the circulating region were determined, and we can use them to ascertain the specific angular momentum of the envelope destined to contract. These values are shown in Figure 3 as green circles for cases presented by Korycansky & Papaloizou (1996). The cases they chose are for rather small Bondi radii, namely, R_B/R_H ⩽ 3^1/3 = 1.44, which nevertheless could correspond to a Jovian size planet in a solar nebula of scale height h/a ∼ O(0.1). Interestingly, they are comparable to the simple estimate of our heuristic model when R_B < R_H.

D'Angelo et al. (2003; see also Estrada et al. 2009) used a nested grid hydrodynamics code to study gas flow in the vicinity of accreting and non-accreting planets. The surface density and azimuthally averaged specific angular momentum, j, were determined as a function of distance from the planet. Ayliffe & Bate (2009) have also computed the specific angular momentum inside the Hill sphere, but used a 3D smoothed particle hydrodynamics code. Both works adopt a thin circumstellar disk of scale height h/a = 0.05 so that for a Jupiter mass planet, the Bondi radius is 5.7 times the Hill radius. The peak values of j ∼ (4, 3) × 10¹⁷ cm² s⁻¹ for planet orbiting material found in these studies, respectively, correspond to peak biases of ℓ_o ≡ j/R²_HΩ ∼ 0.8, 0.6, which, if identified with the outer radius of the infall, imply that r_o = ℓ²_o/3R_H ∼ 0.21R_H, 0.12R_H. Although the disks actually extend farther than this in radius, they are not centrifugally supported structures in their outer portions.

To relate these findings to an average specific angular momentum, j_c for these disks, we have digitized the disk surface density and angular momentum curves found in Estrada et al. (2009) and Ayliffe & Bate (2009) and used them to compute the total disk mass and angular momentum; their ratio then gives the average specific angular momentum, j_c. The resulting biases are ℓ ∼ 0.69 and 0.56, giving centrifugal radii of 118R_J and 78R_J, respectively, and these are indicated in Figure 3 as gray and red points at R_B/R_H = 5.7. Ayliffe & Bate also consider two other planetary masses, 0.3M_J and 0.5M_J, for the same nebula scale height, although the smallest mass does not always form a circumplanetary disk (perhaps because the high accretion rate they adopt results in very hot structures). These correspond to R_B/R_H = 2.6 and 3.6, and the bias values obtained by digitizing their curves are ℓ ≈ 0.40 and 0.48, also included in Figure 3 as red dots. Since as Korycansky & Papaloizou point out, the problem is scalable by the single parameter R_H/h = (R_B/3R_H)^1/2, we can scale these results to a full Jupiter size, where the same R_B/R_H applies if h/a = 0.075, 0.063. It does not seem unreasonable that scale heights this large or even higher could have prevailed within a few Hill sphere distances of Jupiter given its energetic interaction with the solar nebula. These ℓ values yield centrifugal radii of ∼(39, 57)R_J.

Machida et al. (2008; see also Machida 2009) have also investigated angular momentum accretion onto a giant planet by means of a 3D nested grid code using a much higher spatial resolution inside the Hill sphere than other studies. They report the average specific angular momentum directly for several planetary masses and these are shown in Figure 3 by open blue circles. They find that the resulting values can be fit with the scaling j_c = 1.4 × 10¹⁷(M/M_J) cm² s⁻¹ for masses ≲M_J in a solar nebula with h/a = 0.05 at 5.2 AU. This gives a bias of ℓ = 0.29(M/M_J)^1/3, and since these authors also point out that the problem is scalable by R_H/h ∝ (R_B/R_H)^1/2 ∝ M^1/3, we can deduce that ℓ ∼ 0.12(R_B/R_H)^1/2 for R_B/R_H ⩽ 5.7. Although the slope of this curve is quite similar to that found for our heuristic model, the values fall below those estimates as well as the results of the other numerical simulations. However, Machida et al. (2008) and M. N. Machida (2009, private communication) did not correct for the rotating frame of the shearing sheet formalism employed. This can be inferred from their finding that the cumulative angular momentum, J, inside a sphere of their initial unperturbed distribution is negative (their Section 4.4). To correct for the fact that the coordinate frame is rotating at frequency Ω, one would then have to include an additional angular momentum of IΩ in the inventory, where I = (2/5)MR² is the moment of inertia of a constant density sphere, giving a total positive J = (1/10)MR²Ω.

To make an analogous rotational correction to their accreted disks, one must know the details of the disk's structure, since a constant density cannot be used. If we make the assumption that the disk is flattened and replace the spatial density with a surface density description, σ(r), the angular momentum is

$$\begin{eqnarray} J = {\int_0^{R_H } {2\pi r^2 } \sigma } v_\varphi dr &=& \int_0^{R_H } {2\pi r^2 \sigma \left({v_{\varphi,R} + r\Omega } \right)} dr\nonumber\\ & = & J_R + I_d \Omega, \end{eqnarray} \tag{ 50 }$$

where $v_{\varphi^\prime},v_{\varphi,R}$ are the azimuthal velocities of disk material in the inertial and rotating frames respectively, and $I_d = 2\pi \int_0^{R_H } {r^3 \sigma dr}$ is the moment of inertia of the disk. Dividing (50) by $M = 2\pi \int_0^{R_H } {r\sigma dr},$ the first term should recover the values, j_R = J_R/M, found by Machida et al. (2008), while the second term is the desired correction $j_{{\rm cor}} = \Omega \int_0^{R_H } {r^3 \sigma dr/\int_0^{R_H } {r\sigma dr.} }$ Assuming a power-law disk, σ ∝ r^−s, the correction factor reduces to (2 − s)/(4 − s)R²_HΩ (for s ≠ 2, 4). Unfortunately, an azimuthally averaged surface density profile is not given directly in Machida et al. (2008) although the cumulative disk mass as a function of radius is plotted logarithmically. There one finds a curve that is quite flat with d(log M)/d(log r) ≈ 0.3, implying that M ∝ r^0.3, dM ∝ r^−0.7dr. Since dM = 2πrσdr as well, we infer that σ ∝ r^−1.7, which is considerably steeper than that given by Ayliffe & Bate (2009), although the latter's curves do steepen quite a bit at smaller r. Setting s ∼ 1.7 yields a correction factor of ∼0.13R²_HΩ, and with this, the angular momentum bias in the inertial frame becomes ℓ ≈ 0.12(R_B/R_H)^1/2 + 0.13. This is also included in Figure 3 as the blue dashed curve and is remarkably similar in both values and slope to our heuristic values (indeed, a better concurrence than we have any right to expect). The rotation corrected individual values are also shown as filled blue circles.

We note in general that inferring the infall angular momentum from the accumulated circumplanetary disk in hydrodynamic calculations is complicated by the fact that some viscous evolution of the disks is likely in these simulations and this tends to increase the disks' average specific angular momentum because mass accreted by the planet has very low j = (GMR)^1/2, thus fractionally decreasing the mass of the disk more than its angular momentum. For example, D'Angelo et al. (2003) use a constant viscosity of 10¹⁵ cm² s⁻¹ for some of their modeling for which the time scale of viscous evolution, τ_ν ≈ r²_d/ν = 1.6 × 10⁻³(r_d/R_J)² yr, is only ∼O(10) yr for r_d ∼ 10²R_J. The effective viscosity for the SPH simulations of Ayliffe & Bate (2009) depends on the smoothing length of the particles which is not given, but some viscous evolution of their circumplanetary disks is likely to have occurred as well (M. R. Bate 2009, private communication). Machida et al. (2008) do not include a viscosity term in their equations of motion. There is, of course, a numerical viscosity associated with a grid code, although Machida (2009) argues that it is small enough to be ignored for their simulations.

Based on modeling efforts to date, we infer that the inflow could fall in a fairly compact pattern for a sufficiently small Bondi radius, helping explain the scale of the satellite systems. The calculation is challenging (e.g., accounting for the viscous evolution of the circumplanetary disk material may well be required and has not been done to date), and the ultimate values will be dependent on the local nebula conditions, which are likely to remain uncertain. However, the bias is probably larger than the pure shear estimate of 1/4, particularly for Jupiter. For a full mass Jupiter, 2 < R_B/R_H < 3 for 0.085>h/r > 0.069, suggesting a bias between 1/3 and 2/5; for Saturn, 1 < R_B/R_H < 2 for 0.079>h/a > 0.056 suggesting ℓ more in the range of 1/4–1/3. The lower of these Bondi values imply r_c ∼ 28R_J, 23R_S, while the higher values correspond to r_c ∼ 40R_J, 40R_S. The outer boundary of the inflow, r_o, would still be quite small compared to R_H for Δj/j ∼ 1 and could be similar in scale to the regular satellite systems if the range of incoming Δj/j was smaller than order unity used in Section 3. We look forward to further, high-resolution simulations of giant planet gas accretion for a variety of planet masses and disk conditions, which could help tighten the constraints on this important issue.

Returning to the circumplanetary disk model, Equations (39) and (45) can be used to determine the time evolution of the accompanying disk structure for a critically rotating planet as functions of the angular momentum bias, the total inflow rate, $\mathscr{F}\left(t \right),$ and a planetary contraction model. We begin in Section 5.1 by describing how to implement our disk formation machinery with a numerical planet contraction model such as that of Papaloizou & Nelson (2005). We then examine the case of a contracting polytrope in Section 5.2. This polytropic model retains much of the fundamental physics of the planet–disk system but has the advantage of furnishing analytical formulae for the accretion and contraction rates. Finally, we consider the possibility of a relatively brief envelope collapse phase and the subsequent post-collapse evolution.

5.1. A Numerical Contraction Model

5.2. An Inflow Model, $\mathscr{F}(t)$

To implement the polytropic and/or collapse approaches we need to furnish a specific model for the time-varying mass inflow rate, $\mathscr{F}(t)$. Early on, we assume that the planet fills its Hill sphere and accretion becomes regulated by its luminosity, i.e., a type Model A object. However, because our interest is in the establishment of a circumplanetary disk, we wish to follow models that track the system as the planet separates from this boundary, i.e., $\dot R < \dot R_H,$ until the inflow ends as the nebula dissipates. The following discussion describes three possible growth phases that we term shear-limited, diffusion-limited, and torque-limited and develops a simple parameterization for $\mathscr{F}(t)$.

We start by considering a planet orbiting within a locally uniform circumstellar disk with surface density Σ. At first, the accretion rate is only limited by the Keplerian shear of the gas that brings material azimuthally to the planet and, as given in Section 4, is of order $\mathscr{F}_{{\rm acc}} \approx 2 \times \left({3/4} \right)\Sigma ({x_{\max }^2 - x_{\min }^2 })\Omega = 3\Sigma \bar x\Delta x\Omega,$ where $\Delta x,\bar x$ are the width and mid-point of the accretion annulus, and the lead factor of 2 comes from assuming material is accreted from both exterior and interior to the orbit. Using Figure 2 as a guide, we estimate x_max ∼ 2.3R_H, $x_{\min } \sim 1.6 R_H \Rightarrow \bar x \sim 2R_H,$ Δx ∼ 0.7R_H, and set $\mathscr{F}_{{\rm acc}} \approx 4\Sigma R_H^2 \Omega.$ As the planet gains mass, R_H ∝ M^1/3increases and its accretional reach extends further into the disk. However the half-width, w, of the disk annulus that has been cannibalized increases as well, varying roughly linearly with the mass, dM = 4πΣadw, and eventually the rapid growth rate begins to stall when the planet reaches the so-called isolation mass for which these distances become comparable,³

$$\begin{equation} \frac{{M_{{\rm isol}} }}{{M_\star }} = \frac{8}{{\sqrt 3 }}\left({\frac{{x_{\max } }}{{R_H }}} \right)^{3/2} \left({\frac{{\pi \Sigma a^2 }}{{M_\star }}} \right)^{3/2}. \end{equation} \tag{ 51 }$$

To increase its mass beyond this point, the planet must rely on the ability of the disk material initially outside its accretion annulus to viscously diffuse into its vicinity. The planet's growth rate then becomes diffusion-limited, for which $\mathscr{F}_{{\rm diff}} \approx 3\pi \Sigma \nu_{{\rm neb}}$ is a more appropriate rate.⁴ Locally, the shear rate must still apply to the diffusively delivered material, implying that for $\mathscr{F}_{{\rm diff}} < \mathscr{F}_{{\rm acc}},$ the steady-state surface density right at the planet, Σ_p, falls below that of the ambient disk such that Σ_p/Σ ≈ 3πν_neb/4R²_HΩ. As the planet's growth continues, it may become large enough that tidal torques exerted directly on the disk by the planet begin to oppose the diffusion process. A balance between the planet's density wave torques and viscous stresses in the disk leads to a more local and steeper edged gap that further lowers the density being accreted by a factor of exp(−M/M_ν),where

$$\begin{equation} \frac{{M_\nu }}{{M_\star }} \approx \frac{{6\pi }}{5}\left({\frac{{\nu_{{\rm neb}} }}{{a^2 \Omega }}} \right)\left({\frac{{\bar x}}{{R_H }}} \right)^3 \approx 9.6\pi \alpha \left({\frac{h}{a}} \right)^2 \end{equation} \tag{ 52 }$$

is the so-called viscous mass (e.g., Ward & Hahn 2000) at which the planet's torques begin to create a significant density depression, O(1/e)Σ, and the final form in Equation (52) assumes an alpha-type viscosity law, ν_neb = αch (Shakura & Sunyaev 1973). As planetary torques become the bottle neck of the inflow, diffusion re-supplies the region outside the torque-induced gap and the density there rises again to more or less the disk ambient value. This is the torque-limited phase and continues until the disk is dissipated, which is assumed to occur on a characteristic time, τ_neb, i.e., $\Sigma \approx \Sigma _o e^{ - t/\tau _{{\rm neb}} }.$

Assembling these considerations into a composite, the following prescription for $\mathscr{F}(t)$ is adopted:

$$\begin{eqnarray} &&{\rm Shear - limited:}\nonumber\\ &&\mathscr{F}_{{\rm acc}} \approx 4\Sigma _o R_H^2 \Omega e^{ - M/M_{\nu} } e^{ - t/\tau _{{\rm neb}} } \quad {\rm for}\quad M \le M_{{\rm isol}} \nonumber\\ &&{\rm Diffusion - limited:}\nonumber\\ && {\mathscr{F}_{{\rm diff}} \approx 3\pi \Sigma _o \nu_{{\rm neb}} e} ^{ - t/\tau _{{\rm neb}} } \quad {\rm for}\quad M_{{\rm isol}} \le M \le M_{{\rm tq}} \\ &&{\rm Torque - limited:}\nonumber\\ && \mathscr{F}_{{\rm acc}} \approx 4\Sigma _o R_H^2 \Omega e^{ - M/M_v } e^{ - t/\tau _{{\rm neb}} } \quad {\rm for}\quad M_{{\rm tq}} \le M. \nonumber \end{eqnarray} \tag{ 53 }$$

Note that the same functional form is used in the first and third expressions, but this rate tends to be limited more by the R²_H factor at low mass (shear-limited), and by the $e^{ - M/M_{\nu} }$ factor at high mass (torque-limited). The ratio of the two expressions employed in Equation (53) is

$$\begin{equation} \frac{{\mathscr{F}_{{\rm acc}} }}{{\mathscr{F}_{{\rm diff}} }} = \frac{4}{{3\pi }}\!\left(\!{\frac{{R_H^2 \Omega }}{{\nu_{{\rm neb}} }}} \!\right)e^{ - M/M_\nu } \approx 6.15\left({\frac{{M_\star }}{{M_\nu }}} \right)^{1/3} \!\!\left({\frac{M}{{M_\nu }}} \right)^{2/3} \!\! e^{ - M/M_{\nu},} \end{equation} \tag{ 54 }$$

where Equation (52) was used to eliminate ν_neb in favor of M_ν. Figure 4(a) shows this for M_ν = 0.3M_J and indicates the domain of each growth phase. The transition mass, M_tq,between the diffusion and torque-limited phases is found by setting this ratio to unity to find

$$\begin{equation} \left({\frac{{M_{{\rm tq}} }}{{M_\nu }}} \right)^{2/3} e^{ - M_{tq} /M_{\nu} } = 0.162\left({\frac{{M_\nu }}{{M_\star }}} \right)^{1/3}. \end{equation} \tag{ 55 }$$

Note that $\mathscr{F}_{{\rm acc}} /\mathscr{F}_{{\rm diff}}$ and M_tq do not depend on the disk surface density, but the value of the isolation mass does. The more massive the disk, the further to the right the M_isol transition boundary to diffusion-limited growth falls, i.e., at larger planetary mass. In fact for heavy disks, the planet may achieve sufficient mass to initiate gap opening before it hits the isolation mass given by Equation (51). This motivates a generalization of M_isol, obtained by including the exponential term in the relationship $dM \approx 4\pi \Sigma e^{ - M/M_{\nu} } adw.$ Finally, if the predicted gap width w ∼ x_max for the isolation mass is smaller that the scale height of the nebula, Rayleigh instabilities in the steep walls of the narrow gap would vigorously fill it in again so it seems likely that the width should not be less than h. As a simple device to insure this, we replace w ∼ x_max with w ∼ x_max + h. Incorporating these changes leads to

$$\begin{equation} \frac{{e^{M_{{\rm isol}} /M_{\nu} } - 1}}{{\left({x_{\max } + h} \right)/a}} \approx \frac{{4\pi \Sigma _o a^2 }}{{M_\nu }}. \end{equation} \tag{ 56 }$$

In Figure 6, M_isol is plotted as a function of the disk surface density through the parameter πΣa² ≡ M_d. The diffusion-dominated regime is bounded by the M_isol and M_tq curves; elsewhere $\mathscr{F}_{{\rm acc}}$ applies.

Zoom In Zoom Out Reset image size

The inflow model depends on three parameters, τ_neb, Σ_o (or equivalently M_d), and ν_neb (or equivalently, M_ν). It will also prove convenient to introduce yet another parameter that involves their product, namely, $M_D = \int_0^\infty {\mathscr{F}_{{\rm diff}} dt = 3\pi \Sigma _o \nu_{{\rm neb}} \tau _{{\rm neb}} }$. This can be thought of as the total amount of mass that can be delivered to the planet's location via diffusion over the lifetime of the disk. In terms of M_d and M_ν,

$$\begin{equation} M_D = \frac{5}{{2\pi }}\left({\frac{{R_H }}{{\bar x}}} \right)^3 M_\nu \left({\Omega \tau _{{\rm neb}} } \right)\left({\frac{{M_d }}{{M_\star }}} \right) = 52.6 M_\nu ^\prime \tau _6 M_d, \end{equation} \tag{ 57 }$$

where τ_neb = τ₆ × 10⁶ yr and the prime denotes a normalization to a Jupiter mass, i.e., M'_ν = M_ν/M_J. Note that for modest values of M'_v and τ₆, M_D tends to be ≳O(10) larger than M_d because the former refers to a global supply of material while the latter to a more local inventory. Figure 4(b) displays the inflow rate as a functions of time for several different values of M'_d. For high enough surface densities, where Σ = M_d/πa² = 10²M'_d g cm⁻² at 5.2 AU, the diffusion-limited regime can be avoided. Recall, however, that $\mathscr{F}\left(t \right)$ is the rate supplied to the planet by the disk environment and is not always the same as the accreted mass because during the luminosity-limited phase, not all of the offered mass is accepted. On the other hand, a circumplanetary disk cannot exist during this phase.

5.3. A Polytropic Model

We now describe how to build a model for the contraction of a polytrope subject to a known time-dependent inflow $\mathscr{F}\left(t \right).$ This is a useful construct that can be solved in detail analytically and used to explore variations in the predicted planet–disk behavior that arise for different assumed values of the input parameters.

5.3.1. Energy

The potential and thermal energies of a polytrope of degree n with mass M and radius R are

$$\begin{equation} U = - \frac{3}{{5 - n}}\frac{{GM^2 }}{2},\quad \quad E_T = \frac{n}{{5 - n}}\frac{{GM^2 }}{R} \cdot \end{equation} \tag{ 58 }$$

With a rotating object, there is also a kinetic energy contribution, E_R = Iω²/2. For critical rotation, we will adopt E_R = λGM²/2R (see Appendix A for more discussion), and accordingly, set the total energy to be

$$\begin{equation} E \approx - \left({\frac{{3 - n}}{{5 - n}} - \frac{\lambda }{2}} \right)\frac{{GM^2 }}{2} \equiv - p\frac{{GM^2 }}{R} \cdot \end{equation} \tag{ 59 }$$

This can vary in time from energy delivered by accreting material and from energy lost through the planet's luminosity, $\mathscr{L}$. For simplicity, it is assumed that the accreting material falls from rest at infinity so that any falling directly on the planet has a kinetic energy equal and opposite to its potential energy.⁵ Most of the kinetic energy is converted to thermal energy at the shock front, which we assume is then virialized. We also ignore the thermal energy of the infalling material that is originating from the colder solar nebula. Consequently, at this level of approximation, direct infall to the planet does not add net energy. Mass falling onto the disk also arrives with essentially zero net energy, but much of the shock-induced conversion of kinetic to thermal energy is subsequently radiated away by the disk due to its large surface-to-volume ratio. By the time material viscously migrates to the planet, most of its kinetic energy is due to its orbital motion and its specific total energy is negative, i.e., −GM/2R. Thus, we set

$$\begin{equation} \dot E \approx - \mathscr{L} - F_p \left({ - GM/2R} \right), \end{equation} \tag{ 60 }$$

where F_p < 0 implies flow onto the planet. Finally, a non-zero couple between the disk and a critically rotating planet also does work −g_pω_c on the planet that decreases its rotational energy. However, this mechanical energy is converted to thermal energy, and if we again make the simplifying assumption that the preponderance of this is virialized and mixed throughout the planet by convection, the couple does not change the net energy, and Equation (60) is still approximately valid. Other partitioning assumptions could also be accommodated by the formalism developed here, although the computations would proceed somewhat differently.

For a polytrope, the pressure, P, density, ρ, and temperature, T, are related through

$$\begin{equation} P \,{=}\, K\rho ^{1 + 1/n},\quad T = mK\rho ^{1/n} /\Re,\quad T = mK^{n/(n + 1)} p^{1/(n + 1)} /\Re, \end{equation} \tag{ 61 }$$

where m is the molar molecular weight, ℜ = 8.31 × 10⁷ erg mol⁻¹ °K⁻¹ is the universal gas constant, and the proportionality constant is K = N_nGM^{(n − 1)/n}R^(3−n)/n, in which N_n is a numerical value obtained from the Lane–Emden function for polytropic index n (e.g., Chandrasekhar 1958). The pressure at the base of the photosphere is set to P ≈ GM/κR², where κ is the opacity (e.g., Schwarzchild 1958; Papaloizou & Nelson 2005). Substituting these into the third expression of Equation (61) yields the effective temperature,

$$\begin{equation} T_{{\rm eff}} = \frac{m}{\Re }\left({\frac{{GM}}{R}} \right)\left({\frac{{N_n^n R^2 }}{{\kappa M}}} \right)^{1/\left({n + 1} \right)}. \end{equation} \tag{ 62 }$$

From this the planet's luminosity can be estimated,

$$\begin{equation} \mathscr{L} = 4\pi R^2 \sigma _{{\rm SB}} T_{{\rm eff}}^4 = 4\pi R^2 \sigma _{{\rm SB}} \!\left({\frac{m}{\Re }} \right)^4 \!\left({\frac{{GM}}{R}} \right)^4 \!\! \left({\frac{{N_n^n R^2 }}{{\kappa {\rm }M}}} \right)^{4/\left({n + 1} \right)}\!\!, \end{equation} \tag{ 63 }$$

where σ_SB = 5.67 × 10⁻⁵ erg cm⁻² s⁻¹ °K⁻⁴ is the Stefan–Boltzman constant. Note that the luminosity is a power law of the form $\mathscr{L} \propto M^a R^b$ where a ≡ 4n/(1 + n), b = 2(3 − n)/(1 + n).

To determine the contraction rate, E is differentiated with respect to time,

$$\begin{equation} \dot E = - p\frac{{GM^2 }}{R}\left({2\frac{{\dot M}}{M} - \frac{{\dot R}}{R}} \right) \end{equation} \tag{ 64 }$$

and equated to the rate of energy change (60) leading to,

$$\begin{equation} {-} 2\frac{{\dot M}}{M} + \frac{{\dot R}}{R} = - \frac{{\mathscr{L}R}}{{pGM^2 }} + \frac{{1{\rm }F_p }}{{2p{\rm }M}}. \end{equation} \tag{ 65 }$$

The first term on the right-hand side would be −τ⁻¹_KH if $\dot M,$ F_p → 0. Evaluating this term for Jupiter's current mass and radius provides a convenient reference time scale, namely, $\tau _{{\rm KH}J} \equiv \mathscr{L}\left({M_J,R_J } \right)R_J /pGM_J^2.$

5.3.2. Stable Rotation

Recall that for ω < ω_c, the couple at the planet vanishes so that Equation (30) gives

$$\begin{eqnarray} F_p &=& {-} \mathscr{F}\left[ {\left({\frac{{r_o }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_s }}{{r_c }}} \right)^{1/2} } \right]\left({\frac{{j_d }}{{j_d - j_p }}} \right)\nonumber\\ && + \frac{\mathscr{F}}{2}\left[ {\frac{{r_o }}{{r_c }} - \frac{{r_s }}{{r_c }}} \right]\left({\frac{{j_c }}{{j_d - j_p }}} \right). \end{eqnarray} \tag{ 66 }$$

Equation (65) can be combined with Equation (32) to get

$$\begin{equation} \frac{{\dot R}}{R} = - \frac{{\mathscr{L}R}}{{pGM^2 }} + 2\frac{\mathscr{F}}{M}\left[ {\left({\frac{{r_s }}{{r_c }}}\right)^{1/2}}-\left(\frac{r_i}{r_c}\right)^{1/2} \right] + \left({\frac{1}{{2p}} - 2} \right)\frac{{F_p }}{M}, \end{equation} \tag{ 67 }$$

which together with Equation (35) can be integrated to give M(t) and R(t). (Here, we have ignored a small variation in p coming from the non-critical rotation.)

Since these equations are valid only so long as ω < ω_c, one must also monitor the behavior of ω' ≡ ω/ω_c = L/L_c. Differentiating with respect to time, $\dot \omega ^{\prime} = \dot L/L_c - \omega ^{\prime} \dot L{_c}/L_c.$ We find $\dot L$ from Equation (36) with g_p = 0, and $\dot L_C$ from Equation (6) yielding,

$$\begin{equation} \dot \omega ^{\prime} = \frac{1}{\lambda }\left({\frac{1}{2}\frac{\mathscr{F}}{M} \left[ {\frac{{r_s }}{{r_c }} - \frac{{r_i }}{{r_c }}} \right]\left({\frac{{r_c }}{R}} \right)^{1/2} - \frac{{F_p }}{M}} \right) - \omega ^{\prime}\left[ {\frac{3}{2}\frac{{\dot M}}{M} + \frac{1}{2}\frac{{\dot R}}{R}} \right]. \end{equation} \tag{ 68 }$$

Note that if it is assumed that the polytrope starts in an S1 stage, r_s = r_o, F_p = 0, and Equations (35) and (65) reduce to

$$\begin{equation} \dot M = \mathscr{F},\quad \frac{{\dot R}}{R} = - \frac{{\mathscr{L}R}}{{pGM^2 }} + 2\frac{\mathscr{F}}{M}. \end{equation} \tag{ 69 }$$

These are the same mass and radius rates one would get in the case of no angular momentum, but unlike a zero angular momentum case, the spin of the object is changing such that

$$\begin{equation} \dot \omega ^{\prime} = \frac{\mathscr{F}}{M}\left[ {\frac{1}{\lambda }\left({\frac{{r_c }}{R}} \right)^{1/2} - \frac{5}{3}\omega ^{\prime}} \right] + \frac{1}{2}\frac{{\mathscr{L}R}}{{pGM^2 }}\omega ^{\prime}. \end{equation} \tag{ 70 }$$

Exactly when the body becomes rotationally unstable will depend in part on its initial spin rate. If we make the additional assumption of near quasi-equilibrium, $\dot \omega ^{\prime} \approx 0,$ at separation, the starting value is

$$\begin{equation} \omega ^{\prime}_o \approx \frac{3}{{5\lambda }}\left({\frac{{r_c }}{{R_H }}} \right)^{1/2} = \frac{1}{{\sqrt 3 }}\left({\frac{{3\ell }}{{5\lambda }}} \right). \end{equation} \tag{ 71 }$$

This is just what one gets by equating Equation (2) to λMR²_HΩ. If such an object were to then contract without further mass gain, it would become rotationally unstable at R_rot given by Equation (4). However, because further accretion occurs, R_rot must be determined through integration of Equation (68).

5.3.3. Critical Rotation

Once the planet is rotating critically, we can no longer set g_p = 0. Using Equation (32) to eliminate F_p in Equation (65) gives

$$\begin{equation} \left({\frac{1}{{2p}} - 2} \right)\frac{{\dot M}}{M} + \frac{{\dot R}}{R} = -\frac{{\mathscr{L}R}}{{pGM^2 }} + \frac{1}{{2p}}\frac{\mathscr{F}}{M}\left[ {\left({\frac{{r_s }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_i }}{{r_c }}} \right)^{1/2} } \right]. \end{equation} \tag{ 72 }$$

Equations (39) and (72) can now be combined to give first-order equations for M and R separately, namely,

$$\begin{eqnarray} &&\left[ {\frac{2}{\lambda }\left({\frac{{r_d }}{R}} \right)^{1/2} + \frac{1}{{2p}} - 5} \right]\frac{{\dot M}}{M} = - \frac{{\mathscr{L}R}}{{pGM^2 }} + \frac{\mathscr{F}}{M}\Bigg(\frac{2}{\lambda }\left({\frac{{r_d }}{R}} \right)^{1/2} \nonumber\\ &&\quad \times \left[ {1 - \left({\frac{{r_c }}{{r_d }}} \right)^{1/2} } \right] + \frac{1}{{2p}}\left[ {\left({\frac{{r_s }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_i }}{{r_c }}} \right)^{1/2} } \right] \Bigg) \end{eqnarray} \tag{ 73 }$$

$$\begin{eqnarray} &&\left[ {\frac{2}{\lambda }\left({\frac{{r_d }}{R}} \right)^{1/2} + \frac{1}{{2p}} - 5} \right]\frac{{\dot R}}{R} = \left[ {\frac{2}{\lambda }\left({\frac{{r_d }}{R}} \right)^{1/2} - 3} \right]\nonumber\\ &&\quad \times \left({ - \frac{{\mathscr{L}R}}{{pGM^2 }} + \frac{1}{{2p}}\frac{\mathscr{F}}{M}\left[ {\left({\frac{{r_s }}{{r_c }}} \right)^{1/2} - \left({\frac{{r_i }}{{r_c }}} \right)^{1/2} } \right]} \right)\nonumber\\ &&\quad - \left({\frac{1}{{2p}} - 2} \right)\frac{2}{\lambda }\left({\frac{{r_d }}{R}} \right)^{1/2} \frac{\mathscr{F}}{M}\left[ {1 - \left({\frac{{r_c }}{{r_d }}} \right)^{1/2} } \right]. \end{eqnarray} \tag{ 74 }$$

If the inflow is known, these can be solved simultaneously with the luminosity law (63) to obtain the mass and radius of the polytrope with time. Equation (32) can then be used to back-out F_p. Finally, substituting M(t), $\dot M(t)$, and $\mathscr{F}(t)$ into Equation (45) gives the behavior of the couple from which the disk profile can be found for a given viscosity model.

5.3.4. Polytropic Index, n

We now need to consider what is a suitable polytropic index for our problem. In the case of convective-adiabatic equilibrium, n = 1/(γ_a − 1); γ_a = C_p/C_v, where γ_a is the adiabatic index and C_p, C_v are the specific heats at constant pressure and volume. For a diatomic gas, γ_a = 7/5,  n = 5/2; for a monoatomic species, γ_a = 5/3, n = 3/2. Obviously, since a contracting giant planet is not actually a polytrope, no index will be completely satisfactory. However, we do expect in the early stage of giant planet contraction that the predominant form of hydrogen will be molecular, with m ∼ 2 g mol⁻¹, and our calculations will employ an n = 2.5 polytrope throughout much of the evolution. For the luminosity Equation (63), a = 20/7, b = 2/7, and in terms of the opacity, κ, this implies τ_KHJ = 2.2 × 10⁷(κ/10⁻² cm² g⁻¹)^8/7 yr.

5.4. Collapse Phase

Our choice of an n = 2.5 polytrope is predicated on the assumption that molecular hydrogen is the dominant species. On the other hand, as the planet contracts, its temperature rises and hydrogen dissociation will occur at temperatures of ∼3 − 4 × 10³ °K depending on the pressure (e.g., Guillot 2005). The temperature of a polytrope, from Equation (61), can be written as

$$\begin{eqnarray} T &=& N_n \frac{m}{\Re }\left({\frac{{GM}}{R}} \right)\left({\frac{{3\rho }}{{4\pi \bar \rho }}} \right)^{1/n}\nonumber\\ & =& 2.3 \times 10^5 \;^{\circ}{\rm K }\,mN_n \left({\frac{{3\rho }}{{4\pi \bar \rho }}} \right)^{1/n} \left({\frac{{M^{\prime}}}{{R^{\prime}}}} \right), \end{eqnarray} \tag{ 75 }$$

where $\bar \rho \equiv 3M/4\pi R^3$ is the average density, and the mass and radius have been normalized to Jupiter values in the last expression. For n = 2.5, the structure is quite centrally condensed, with a central density $\rho _c = 23.4\bar \rho$ (Chandrasekhar 1958). This implies a central temperature T_c = 3.2 × 10⁵(M'/R') °K. We will see in the next section that for a planet that has attained more than a few tenths of a Jovian mass by the time it contracts within the satellite zone, dissociation can commence near its center, and an n = 2.5 polytrope begins to lose validity.

Dissociation of hydrogen could cause important structural changes. Absorbing molecular kinetic energy by breaking hydrogen bonds enables the body to contract without having to wait for luminosity losses. Early models, e.g., Bodenheimer (1974), predicted that envelope collapse occurred on a time scale, τ_col,of only a few years, which would lead to a massive disk similar to that discussed in Section 2.3.1. However, more recent work, e.g., Lissauer et al. (2009), finds a more gradual collapse over centuries to millenia. The longer of these timescales is marginally within the quasi-steady state regime, τ_col ≈ τ_ν, and can be roughly followed by replacing the polytrope contraction rate (Equation (74)) with $\dot R/R \approx 1/\tau _{{\rm col}}$ between the radius when T_c ≈ 3500 °K until the collapse halts near R ≈ 1.75R_J (Lissauer et al. 2009). Past this point, there is little change in R and we can set $\dot R \approx 0$ for the duration of the disk.

We are now in a position to examine circumplanetary disks for the above models at various points in their evolution and start again with the numerical model of Papaloizou & Nelson (2005) followed in turn by the polytropic and polytropic+collapse+post-collapse models. In each case, a viscosity law will be required to infer the surface density from the couple. For an alpha model, ν = αc²/Ω = αγℜT/μ_molΩ, which depends on temperature as well as distance, where Ω = (GM/r³)^1/2 now refers to the orbital frequency of disk material around the planet. If it is assumed that viscous dissipation per unit area, ∼(9/4)σνΩ², is the primary heat source in the disk and that this is balanced against the radiation losses from both of the disk surfaces, ∼2σ_SBT⁴, one can set T⁴ ≈ 9Ω²σν/8σ_SB = (3/8π)(Ω²/σ_SB)g/j to derive (e.g., Canup & Ward 2002)

$$\begin{equation} \sigma \approx \frac{m}{{3\pi \alpha \gamma \Re }}\left({\frac{{8\pi \sigma _{{\rm SB}} }}{3}} \right)^{1/4} \left({\frac{{GM}}{{r^3 }}} \right)^{1/4} \left({\frac{g}{j}} \right)^{3/4}. \end{equation} \tag{ 76 }$$

6.1. Numerical Model

Figures 5(a) and 5(b) show the planet mass and radius as functions of time from a Papaloizou & Nelson (2005) model. Since we know M, $\dot M,$R, $\dot R$, Equation (39) can be used to deduce an appropriate inflow $\mathscr{F}$. Equation (45) is then used to find 3πσν as a function of time and combined with Equation (77). Figure 5(c) shows the resulting behavior of σ as a function of r at various times during the evolution assuming ℓ = 1/3 and α = 3 × 10⁻³. Typical surface densities throughout the satellite zone are of order ∼10³ g cm⁻². With time, the inflow slows and the planet contracts, causing the inner disk edge to move inward and the peak value for 3πσν to decrease. For the case shown here, once the planet contracts to a radius smaller than about 13R_J, the flux at the planet becomes negative (F_p < 0) and the disk is increasingly well approximated as an accretion disk. Figure 5(d) shows an estimated temperature profile, T(r), at the same time intervals. Temperatures in the outer part of the satellite zone become low enough to permit ice condensation. At late times in a gas giant's growth, the circumplanetary disk model approaches that of a pure accretion disk (see also Canup & Ward 2009).

Zoom In Zoom Out Reset image size

6.2. The Polytropic Model

6.2.1. Final Planet Mass

For a polytropic model, the inflow (Equation (53)) is specified along with the luminosity, and M(t) and R(t) are to be derived. To a good approximation, the planet's mass as a function of time can be found by simply integrating the inflow model, although this ignores the mass in the disk and that lost from its outer boundary. However, Equation (53) is only applicable after the planet has detached from its Hill sphere boundary. To follow the planet as it separates, its mass, M_o (and spin rate, ω_o, for subcritical rotation) at that time are needed as initial conditions for the evolutionary equations. Separation of the polytrope radius from its Hill radius does not require that the body actually contract, but only that it expands less rapidly than R_H, i.e., $\dot R - \dot R_H < 0.$ For simplicity, we ignore F_p and set $\dot M = \mathscr{F}$ in Equation (65). Comparing the resulting $\dot R$ with $\dot R_H = \left({R_H /3} \right)({\dot M/M})$ leads to a criterion for separation

$$\begin{equation} \mathscr{L}R_H /pGM^2 \ge 5\mathscr{F}/3M. \end{equation} \tag{ 77 }$$

This can be cast in terms of a mass, M_o at separation using either $\mathscr{F}_{{\rm diff}}$ or $\mathscr{F}_{{\rm acc}},$ i.e.,

$$\begin{equation} M^{\prime}_d = 56M^{\prime}_o{} ^{16/7} /M^{\prime}_v \tau _{KH6} \quad {\rm for}\quad \mathscr{F}_{{\rm diff}}, \end{equation} \tag{ 78 }$$

$$\begin{equation} M^{\prime}_d = 0.91M^{\prime}_o{ }^{34/21} e^{M_o /M_v } /\tau _{KH6} \quad {\rm for}\quad \mathscr{F}_{{\rm acc}}, \end{equation} \tag{ 79 }$$

where again the primes denote normalization to M_J, and solutions to both are included in Figure 6 as dot-dashed and dashed curves. Self-consistency then dictates that only those portions of each curve that fall in the proper domain are valid; namely, if $M_o \left({\mathscr{F}_{{\rm diff}} } \right)\left[ {M_o \left({\mathscr{F}_{{\rm shear}} } \right)} \right]$ lies in (out of) the diffusion-limited domain. If neither of these conditions are satisfied, the separation mass is set equal to M_isol.

Zoom In Zoom Out Reset image size

Integrating the inflow model from M_o onward gives the mass versus time dependences shown in Figures 7(a) and 7(b), where the angular momentum has been omitted to obtain a base-line case. In Figure 7(a), the surface density of the disk is varied though the parameter M_d. Also shown is M_tq for the viscosity corresponding to M'_ν = 0.3, i.e., α = 2.4 × 10⁻³ for h/a = 0.07. For small disks, the final planet mass is limited by M_D, the amount of material that can be delivered before the disk is dissipated. A Jupiter mass is reached for M_d = 0.02M_J for which M_D = 0.95M_J. At large disk mass, it is gap opening that limits the planet's growth. Figure 7(b) illustrates the effect of lowering the viscosity, which makes gap formation easier. A Jupiter mass can become the limit for M_d = 1M_J, M_D = 15.7M_J if M'_ν = 0.1, α = 8.4 × 10⁻⁴, and it is clear that there are a continuum of possible disk masses and viscosities that could produce a one Jupiter mass end state. In what follows, we will follow two specific cases as representatives of diffusion-limited and torque-limited growth scenarios.

Zoom In Zoom Out Reset image size

6.2.2. Contraction Phase

Figures 8(a), 8(b) and 9(a), 9(b) show M(t) and R(t) curves for diffusion-limited and torque-limited examples respectively, assuming τ_neb = 3 × 10⁶ yr and τ_KHJ = 10⁷ yr. These now take into account the difference between the mass delivery rate obtained by integrating the inflow and the M(t) that occurs because of the presence of the disk physics. This requires a slight adjustment to the conditions described in Section 6.2.1 to reach a final state of M = M_J. In the diffusion-limited case, we require M'_d = 0.028, M'_ν = 0.3 to get a full Jupiter mass for the planet (Figure 8(a)); for the torque-limited case, we obtain a Jupiter mass by adopting M'_d = 1, M'_ν = 0.118 (Figure 9(a)). An angular momentum bias ℓ = 1/3 is employed for both cases. For this value, the n = 2.5 polytrope is rotating critically as it first separates from R_H. The polytrope radii are shown in Figures 8(b) and 9(b), along with the evolving boundaries of the inflow pattern, r_i, r_o and the Hill sphere R_H. The stages whose boundaries occur as the planet radius contracts first below r_o and then below r_i are labeled by C1, C2, and C3. Note that the polytrope contracts to the scale of a giant planet satellite system while the planet's mass is only a modest fraction of its final value. This is reasonably consistent with the numerical modeling results of Papaloizou & Nelson (2005) as well as a more recent work of Lissauer et al. (2009). Figures 8(c) and 9(c) show the time variations of total inflow, $\mathscr{F}$ (black solid), the portion, $\mathscr{F}_p$ (black dashed), falling directly on the planet, the planet's accretion rate, $\dot M$ (green dashed), and the in-plane disk fluxes at the planet, F_p (red dashed), and outer disk boundary, F_d (red solid). At first (stage C1, when R > r_o), F_p = F_d > 0, and material is shed by the planet and then migrates outward at a constant flux to escape from the outer disk edge, r_d. As the system enters stage C2 (r₁ ⩽ R ⩽ r_o), an increasingly larger fraction of the infall lands directly on the disk and the flux rate at the planet, F_p, starts to decrease, eventually switching sign to indicate that material is now being accreted from the disk. Early fluxes in the torque-limited case (Figure 9(c)) are much higher than for the diffusion-limited one (Figure 8(c)) because of the higher accretion rate. However, this period is short lived, and the precipitous drops in $\mathscr{F}$ and $\dot M$ in the torque-limited case signals that gap opening has commenced. Figures 8(d) and 9(d) plot the effective and central temperatures of the polytropes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.2.3. Circumplanetary Disk Profile

In Figures 10 and 11, the disk structure is examined for the diffusion-limited and torque-limited polytropic models, respectively. Figures 10(a) and 11(a) plot the ratio of the couple to the specific angular momentum, g/j = 3πσν. The curves are labeled by time and terminate at the surface of the planet, which becomes smaller as time elapses. In Figures 10(b) and 11(b), the behavior of the in-plane flux, F(r), throughout the disk is exhibited. At the earliest times, it is positive throughout, but as the planet contracts, the flux becomes negative at the planet while remaining positive at the outer disk edge. Once it appears, the stagnation point where F(r) = 0 remains relatively fixed. Figures 10(c), 10(d), 11(c), and 11(d) illustrate the resulting gas surface density and disk temperature profiles. These curves are determined assuming a viscosity strength of α = 3 × 10⁻³. For this α, the surface density of the disk in the satellite region spanning 5.9 < r/R_J < 26.4 is only a few hundred g cm⁻³ or less by the time the planet has contracted within this region in both the diffusion-limited and torque-limited examples. This is far less that the typical MMSN surface densities of O(10⁵–10⁶) g cm⁻² estimated by augmenting satellite masses to solar composition (e.g., Pollack & Consolmagno 1984), but is comparable to the values employed in the more recent starved disk models of Canup & Ward (2002, 2009). Along with low surface densities, one finds cool disk temperatures, because of the smaller viscous dissipation rate compared to an MMSN surface density, although the reader is cautioned that contributions from the planet's luminosity and nebular insolation have not been included in the graphs. Cool temperatures favor ice condensation late in the accretion process, perhaps helping to explain the progressively higher ice content and lower bulk densities with distance of the Galilean satellites' compositions. This issue is discussed in more detail in Canup & Ward (2002, 2009).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.3. Collapse Phase

The central temperatures of the polytropes are plotted in Figures 8(d) and 9(d) and exceed the dissociation temperature when R'/M' ≲ 10². Here we consider the effects of an envelope collapse phase once this occurs. If the collapse phase is fast compared to the viscous time scale, the resulting disk is similar to the inviscid case of Section 2.3.1 and will contain a fair fraction of the total planet's mass. This is dictated by the large angular momentum content of the critically rotating body. For solar composition material, this would be enough solids to form much larger and/or more numerous satellites than observed for the giant planets in our solar system. One could postulate a much lower solid-to-gas ratio in the planet's outer layers, perhaps due to settling, so that only the appropriate amount of satellite forming material was actually present. But in this case, the low solids-to-gas ratio in the disk leads to a satellite accretion time scale that is too slow to prevent the rapid loss of satellites from type I orbital decay before they reach their requisite size (Canup & Ward 2002, 2006).

If τ_col is comparable to the disk viscous time scale, a disk mass that is intermediate to the no-collapse and the rapid collapse cases results. A detailed description of the disk density and temperature would require the full solution of the time-dependent diffusion equation rather than a quasi steady-state approach, but we can venture an order of magnitude estimate of the amount of mass shed by the planet by the following argument: from Equations (6) and (13) we infer that $\left({1 - 3\lambda /2} \right)\dot M = \left({\lambda /2} \right)M\dot R/R \,+\, g_p /\left({{\it GMR}} \right)^{1/2}.$ Recall that in an inviscid disk, the couple is zero and the mass of the spin-out disk is given in Equation (8). The presence of a viscous torque, g_p ≠ 0, exerted by the disk on the planet drains angular momentum from the planet so that less mass needs to be shed to maintain the planet at the edge of rotational stability. If g_p is given by Equation (12), then $g_p /\left({{\it GMR}} \right)^{1/2} = - \dot M[ {\left({r_d /R} \right)^{1/2} - 1} ],$ which leads to Equations (14) and (16). However, this assumes that changes in the physical state of the planet are viscously communicated all the way to r_d on a time scale τ_ν that is short compared to changes in $\dot R/R,{\rm }\dot M/M.$ But for small τ_col < τ_ν, the collapse will be over before this can happen. This suggests that we replace r_d with an effective distance, r_eff,to which the collapse information has diffused in estimating the couple at the planet, i.e., $g_p /\left({{\it GMR}} \right)^{1/2} \approx - \dot M[ {\left({r_{{\rm eff}} /R} \right)^{1/2} - 1} ],$ which then leads to ΔM/M = O[λ(R_col/r_eff)^1/2] instead of Equation (16). (This effective "edge" can be found by replacing τ_ν with τ_col and r_d with r_eff in Equation (81) of the next section.) Assuming a roughly constant aspect h/r, we conclude that r_eff/r_d ≈ (τ_col/τ_ν)^2/3 and ΔM/M ≈ λ(R_col/r_d)^1/2(τ_ν/τ_col)^1/3. Setting R_col/r_d ∼ O(10⁻¹)  yields a total disk mass reminiscent of the MMSN approach if τ_ν/τ_col ∼ O(1), i.e., ΔM/M ∼ 3.5 × 10⁻². However, there is an important difference between this situation and the traditional MMSN disk. The regulation of the spin-out mechanism by the viscous couple applies to both gas and solids, but once shed, the solids may coagulate and decouple from the gas, remaining in the vicinity of where they were spun out. The gas, on the other hand, must spread out to nearly the disk edge to maintain the couple. Thus, the amount of gas remaining in the vicinity of the solids is more like ∼ΔM(R_col/r_eff)² ∼ 3.5 × 10⁻⁴M, resulting in comparable gas and solid surface densities of order ∼10³ g cm⁻².

We examine this situation in Figures 12 and 13, which repeat the planet/disk evolutions of Figures 8 and 9, but include a collapse phase, $\dot R/R = - 1/\tau _{{\rm col}}$ on a time scale comparable to the viscous time, τ_col ≈ τ_ν, that is initiated when the planet's central temperature reaches ∼3500 °K. The viscous time is ∼10⁴ yr for α = 0.003 (see Equation (81)), which would require a more protracted collapse than seen in recent models (Lissauer et al. 2009). The timing of this is displayed in Figures 12(b) and 13(b) where it is apparent that the planet passes through the C2 stage very quickly. Figures 12(d) and 13(d) display the various flux behaviors during the collapse episode on a finer resolution. (The artificially sharp discontinuities are artifacts of abruptly switching the $\dot R$ prescription from the polytrope description to the collapse description and then again to the post-collapse state.) At the onset of collapse, there is a sudden spike in the outward flux, while at the same time the planet's mass accretion rate, $\dot M,$ goes negative. The effect of this is just discernable on the M versus t plots shown in Figures 12(a) and 13(a). The structure during the collapse is primarily a renewed spin-out disk with a sudden increase in F_p > 0 that temporarily interrupts the disk's transition to that of an accretion disk. These rates then decay over the duration of the collapse, with F_p going negative again at about its mid-point. There is a second point of discontinuity in the $\dot R$ prescription that starts the post-collapse epoch, and over time scales much in excess of τ_col, the accretion disk behavior resumes except that the flux at the planet disk boundary F_p < 0, is more negative than was found in Figures 8 and 9. This is because the radius of the planet is already near its final value and has cleared the satellite zone having passed stage 2 quickly. The inflow now falls directly on the disk and must diffuse to the planet to be accreted. Figure 14 contrasts the resulting surface density and temperature profiles of the disk at three separate times selected to be prior to, during, and after the collapse for the diffusion-limited case (Figures 14(a) and 14(b)) and the torque-limited case (Figures 14(c) and 14(d)). The surface density reaches somewhat higher values, ∼10³ g cm⁻², during collapse, but the condition persists for an only relatively brief period. Were satellites to form during the collapse interval, their survival during the post-collapse phase would seem unlikely since several tenths of the planet mass has yet to be accreted.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have developed a self-consistent model for the growth of a gas giant planet and the development of an accompanying circumplanetary disk. The model is composed of three elements: (1) an inflow model describing the properties of the material inflowing to the planet, primarily the inflow rate and its specific angular momentum, (2) a viscous quasi steady-state disk model that describes the in-plane flux, surface density and temperature of the circumplanetary disk at any given time, and (3) a planet growth and contraction model that can follow the evolution of its mass and radius. These components perform as modules that are linked through appropriate boundary conditions and evolved in concert. Each could be swapped out for alternative versions if and when improvements in them become available.

The specific angular momentum of inflowing gas is estimated using both a heuristic analytical model and results from recent 3D hydrodynamical simulations. We describe the planet's contraction history both with published numerical works and with an analytical model that treats the planet as a uniformly rotating polytrope. To describe the time-dependent mass inflow rate in our analytical model, we consider three regimes of gas giant planet growth: shear-limited, diffusion-limited, and torque-limited. The accretion rates associated with these regimes and the conditions when each predominates are themselves functions of the planet's mass and the circumstellar disk mass, viscosity, and lifetime.

A suite of planet and nebula histories can yield a Jupiter-mass planet before the nebula dissipates and gas inflow ends. We consider two limiting cases in our polytropic model that both produce a final planet with M ∼ M_J. The first involves a circumstellar disk with a moderate viscosity and a relatively low initial mass. In this "diffusion-limited" case, the final stages of the planet's growth are regulated by the rate at which material can diffuse into its local region. The second case involves a circumstellar disk with a higher mass and somewhat smaller viscosity. This case results in a planet whose final growth is "torque limited", regulated by the strength of its gap-opening torques on the surrounding nebula.

The coupled evolution of the planet and its circumplanetary disk displays common features for a numerical planet contraction model from Papaloizou & Nelson (2005) and for both the diffusion-limited and torque-limited histories in our polytropic model. These include:

• 1.  
  As a critically rotating planet contracts, an initial circumplanetary disk is produced via spin-out from the planet's equatorial region. Once the planet contracts enough that some of the inflowing gas contains too much angular momentum to fall directly onto the planet, the disk begins to be supplied by a combination of spin-out from the planet and direct inflow from the nebula. Eventually the disk transitions to an accretion disk supplied solely by the inflow and material flows from the disk onto the planet at the disk's inner edge.
• 2.  
  For a circumplanetary disk viscosity parameter α ∼ 10⁻³ and a final planet contraction timescale ⩾10⁷ yr (corresponding approximately to a planet opacity ⩾10⁻²cm² g⁻¹), peak gas surface densities in the circumplanetary disk throughout most of its history are ∼10²–10³ g cm⁻². Higher viscosities and/or larger planet opacities typically yield even less massive circumplanetary disks by the time the planet contracts to a radius smaller than the satellite zone. These low densities are similar to those invoked in "gas-starved" disk models (Canup & Ward 2002, 2006, 2009), and are orders-of-magnitude lower than that of the canonical MMSN, in which σ ∼ 10⁵–10⁶g cm⁻². To reach MMSN disk densities, either a very small circumplanetary disk viscosity (e.g., α ≲ 10⁻⁶) compared to the circumstellar disk is required, or else a fairly massive solar nebula must be invoked, although in this latter situation, the high disk densities occur when the planet is only a small fraction of its final mass (see below).
• 3.  
  Temperatures in the circumplanetary disk as the planet nears its final mass are typically low enough for ices at orbital radii greater than about 10–20 planetary radii, consistent with the formation of ice-rich satellites in these regions.
• 4.  
  The history of the system may contain a relatively brief collapse phase during which the radius of the planet shrinks through the satellite zone to nearly its final radius. Such an event is punctuated by a somewhat higher density spin-out disk, but this decays to the more typical long-term densities on the viscous time scale of the disk. For the models considered here, the collapse phase occurred well before the planet had achieved its final mass, which is consistent with recent numerical modeling, for instance Lissauer et al. (2009).

The basis for the prevailing low circumplanetary disk surface densities can easily be identified using the polytropic model. From Equation (45), it is apparent that the magnitude of 3πσν scales with $\mathscr{F}$. It follows that if the inflow rate is diffusion-limited, $\mathscr{F}_{{\rm diff}} = 3\pi \Sigma \nu_{{\rm neb}},$and σ will be of order Σ(ν_neb/ν). For comparable alpha values, ν_neb/ν = T_neb/T_disk × Ω_p/Ω(r), where Ω_p is the planet's heliocentric orbit frequency and Ω(r) is the circumplanetary orbit frequency of the disk material. Evaluating the circumplanetary disk frequency at the centrifugal radius gives Ω_p/Ω(r_c) = 9/ℓ³. For ℓ = 1/3, this reads 2.4 × 10² while the temperature of the cicumplanetary disk is a factor of 2 or more higher than the nebula implying that σ ≈ 10²Σ ≈ 10⁴M'_d g cm⁻² At first, it might appear that setting M'_d ∼ 10–10² would give MMSN-type values, but there is a constraint on this choice. If the planet's mass is diffusion-limited, we must also satisfy $M_D = 3\pi \Sigma_o \nu_{{\rm neb}} \tau _{{\rm neb}} \approx M_J - M_o \left({\mathscr{F}_{{\rm diff}} } \right).$ Using this to eliminate $\mathscr{F}_{{\rm diff}},$ and assuming M_o ≪ M_J leads to

$$\begin{eqnarray} \sigma &\approx& \frac{{M_J e^{ - t/\tau _{neb} } }}{{3\pi v\tau _{neb} }} \approx \frac{1}{3}\left({\frac{{M_J }}{{\pi r_d^2 }}} \right)\left({\frac{{\tau _\nu }}{{\tau _{neb} }}} \right)e^{ - t/\tau _{neb} } \nonumber\\ &\approx& 10^2 \left({\frac{{3 \times 10^{ - 3} }}{\alpha }} \right)\left({\frac{r}{{10h}}} \right)^2 \left({\frac{{e^{ - t/\tau _{neb} } }}{{M^{\prime 2/3} }}} \right)\;{\rm g}\;{\rm cm}^{ - 2}, \end{eqnarray} \tag{ 80 }$$

where in evaluating (80) we have set

$$\begin{eqnarray} \tau _\nu &\equiv& \frac{{r_d^2 }}{\nu } \approx \frac{1}{\alpha }\left({\frac{{r_d }}{h}} \right)^2 \frac{1}{{\Omega \left({r_d } \right)}} = \frac{{\gamma ^{3/2} }}{{3^{1/3} \alpha }}\left({\frac{{r_d }}{h}} \right)^2 \frac{1}{{\Omega _p }} \nonumber\\ &\approx& 10^4 \left(\frac{\gamma }{0.5}\right)^{3/2}\left({\frac{{3 \times 10^{ - 3} }}{\alpha }} \right)\left({\frac{{r_d }}{{10h}}} \right)^2 \;{\rm yr}{\rm.} \end{eqnarray} \tag{ 81 }$$

Because σ varies inversely with the viscosity of the circumplanetary disk, lowering its alpha value compared to the circumstellar disk would increase the circumplanetary disk surface density, but values as low as α ≲ O(few × 10⁻⁶) would be necessary to achieve MMSN values. This would imply quite different turbulence properties for the circumplanetary environment versus those typically assumed for circumstellar disks. Even then, the disk is not really equivalent to the more or less static MMSN paradigm because it is continuously being accreted by the planet and supplied by the inflow.

If, on the other hand, the inflow rate is shear/torque-limited, we would instead expect

$$\begin{eqnarray} \sigma &\approx& \Sigma \left({\frac{{\nu_{{\rm neb}} }}{\nu}} \right)\frac{{\mathscr{F}_{{\rm acc}} }}{{\mathscr{F}_{{\rm diff}} }}\nonumber\\ &=& \left[ {\frac{4}{{3\pi }}\left({\frac{{M_J }}{{\pi a^2 }}} \right)({\Omega _p \tau _\nu })} \right]M^{\prime}_d \gamma ^{ - 2} e^{ - M/M^{\prime}_\nu } e^{ - t/\tau _{{\rm neb}} }. \end{eqnarray} \tag{ 82 }$$

The bracketed quantity has a value of 2.4 × 10⁵(τ_ν/10⁴ yr)g cm⁻² so again one could seemingly pick values of M'_d, M'_ν for which the resulting circumplanetary disk surface density is of order the MMSN values except that these parameters are still constrained by the requirement that the final planetary mass be M_J. Let us assume that the disk mass is great enough that the diffusion-limited regime is completely avoided during accretion. Integrating $\dot M = \mathscr{F}_{{\rm acc}}\ {\rm to }\ t \to \infty,$ and requiring the final mass be M_J leads to the constraint

$$\begin{eqnarray} \frac{{M^{\prime}_d }}{{M^{\prime}_{\nu} }} &=& \frac{\pi }{4}\left({\frac{{3M_\star }}{{M_J }}} \right)^{2/3} \frac{{M_{\nu}^{ - 1} }}{{\Omega _p \tau _{{\rm neb}} }}\int_0^{M_J - M_e } M^{ - 2/3} e^{ - M/M_\nu } dM \nonumber\\ &=& \frac{\pi }{4}\left({\frac{{3M_\star }}{{M_J }}} \right)^{2/3} \frac{{e^{x_f} {I\left({x_f} \right)} }}{{\Omega _p \tau _{{\rm neb}} }}, \end{eqnarray} \tag{ 83 }$$

where $I\left({x_f } \right) = x_f^{2/3} \int_0^{x_f } {e^{ - \left({x_f - x} \right)} x^{ - 2/3} } dx$ is a slowly varying function of x_f = (1 − M'_o)/M'_ν whose value ranges from 1.03 to 1.68 over 1 < x_f < 25 during which $e^{x_f }$ changes by several orders of magnitude. Substituting Equation (83) into (82) gives

$$\begin{eqnarray} \sigma &\approx& \left({\frac{{M_\star }}{{\pi a^2 }}} \right)\left({\frac{{M_J }}{{3M_\star }}} \right)^{1/3} \left({\frac{{\tau _\nu }}{{\tau _{{\rm neb}} }}} \right)M^{\prime}_{\nu} e^x I({x_f })e^{ - M/M_\nu } e^{ - t/\tau _{{\rm neb}} } \nonumber\\ &\approx& 300M^{\prime}_\nu e^{1/M^{\prime}_v } e^{-(M + M_{\rm o})/M_{\nu}} \,{\rm g}\,{\rm cm}^{ - 2} \end{eqnarray} \tag{ 84 }$$

where we have set α = 3 × 10⁻³, r/10h = 1, and γ = 1/2. This can still get quite large for small enough M'_ν, namely, $300M^\prime _\nu e^{1/M^\prime _\nu } = ({6.7 \times 10^5,6.4 \times 10^6 })$ for M'_ν = (0.1, 0.08) respectively. However, these peak surface density values occur early in the accretion, when the planet mass is still rather small. Indeed, it is likely that the planet has not yet contracted within the satellite zone. As the planet approaches it final mass, the inflow wanes and the circumplanetary disk density decreases with the factor $e^{ - \left({M + M_o } \right)/M_\nu }.$ For a planet at 50% of its final mass, these densities have dropped to (4.4 × 10³, 1.2 × 10⁴) g cm⁻²; at M + M_o = 0.9M_J, the densities are only of order 10² g cm⁻². Thus, low circumplanetary disk gas surface densities are an inevitable consequence of the progressive decay of the inflow rate near the planet's terminal mass in either the diffusion-limited or torque-limited regimes.

Is there a potential to create an MMSN-style disk through the onset of a rapid collapse phase? We remarked in Section 6.3 that hydrogen dissociation temperatures are reached at the center of a polytrope when M'/R' = 10². If this occurs when the planet's radius is comparable to the satellite system (R' ∼ 25), the planet would have achieved only a fraction of its final mass at the time of the collapse generated spin-out disk. The later arrival of several tenths of the planet's mass onto the disk would completely alter the outcome, with structure relaxing back to an accretion disk in its post-collapse stage. If, on the other hand, the opacity of the planet was so high that its contraction occurred subsequent to its accretion, i.e., M' ∼ O(1), the collapse would ensue when R' ≈ 10², which is rather large compared to the observed scale of the satellite systems. In addition, we find that the polytropic model requires a contraction time scale parameter, τ_KHJ, some ∼ 300 times τ_neb, for accretion to be ∼97% complete before collapse. This implies such high planet opacities, i.e., κ ∼ O(few) cm² g⁻¹, that it is problematic as to whether the planet could form within the nebula's lifetime (Hubickyj et al. 2005). Thus, while a spin-out disk of sufficient solids to form the satellites appears just attainable for rather extreme values of the model's parameters, the range of parameter space resulting in a gas starved environment at the end of planet formation is much larger. In addition, the only partially differentiated states of Callisto (Anderson et al. 2001) and Titan (Iess et al. 2010) argue against their formation in an environment where all the necessary solids are supplied in a collapse time scale, because this would allow their accretion to proceed very quickly, resulting in their complete differentiation.

A key parameter throughout the modeling of the evolution of the planet–disk system is the angular momentum of the inflowing material. A prior estimate (Lissauer 1995) for a purely Keplerian flow predicted an angular momentum bias of ℓ = 1/4, corresponding to centrifugal radii at Jupiter and Saturn of 15R_J and 23R_S, respectively. Based on results of recent hydrodynamic simulations (D'Angelo et al. 2003; Machida et al. 2008; Ayliffe & Bate 2009) and a heuristic model described here that roughly accounts for both pressure forces and the planet's gravity, we estimate an average bias in the range of ℓ = 1/3–2/5 for a fully formed Jupiter, and l = 1/4–1/3 for a fully formed Saturn. Such values would imply centrifugal radii at Jupiter between 28 and 40R_J, and between 23 and 40R_S at Saturn. These radii are vastly smaller than the Hill radius of either planet, and are roughly comparable to the compact scale of the regular satellite systems. Determination of the centrifugal radius is a challenging calculation: it is a sensitive function of the inflow specific angular momentum, which in turn depends on somewhat uncertain local nebular conditions (e.g., the nebular scale height and temperature). Future high-resolution hydrodynamical simulations that explore a wider range of conditions are needed to provide improved constraints on this quantity.

Finally, the coupled planet–disk evolution model described here assumes that the viscous torque between the disk and the surface of the planet will maintain the planet's rotation rate at the limit of rotational stability as the planet contracts. Yet ultimately, Jupiter and Saturn must be left with less angular momentum than this to account for their current ∼10 hr rotational days, which are about a factor of 3 longer than their critical rotation periods. Accounting for sub-critical rotation is a well-known issue for Jupiter and Saturn, as well as for protostars. Protostars, also believed to grow through mass delivered via a viscous accretion disk, have observed rotation rates that are often much slower than breakup (e.g., Herbst et al. 2002). Proposed solutions to account for protostar angular momentum loss may also apply to gas giant planets, including (1) "disk-locking", or magnetic coupling between the central object and its disk that results in angular momentum transferred to the disk beyond the co-rotation radius (where the Kepler velocity equals the primary rotational velocity; e.g., Koenigl 1991; Takata & Stevenson 1996) and (2) magnetocentrifugally driven "X-winds" that originate from the inner accretion disk, diverting mass and angular momentum that would otherwise be delivered to the primary (e.g., Shu et al. 2000). A full treatment of this issue in the context of the current model will involve, e.g., modifying the disk profile to include the magnetic torque and assessing the expected degree of ionization in the disk, which are planned topics of future work.

This work was supported by grants from NASA's Outer Planets Research Program and NSF's Astrophysics Program (WRW) and by a grant from NASA's Origins of Solar Systems Program (RMC).

A rapidly spinning object will flatten into the form of a spheroid. The general form for the potential of a spheroid in its midplane is (e.g., Murray & Dermott 1999)

$$\begin{equation} U = - \frac{{GM}}{r}\left[ {1 - \Sigma _{j = 1}^\infty J_{2j} \left({\frac{R}{r}} \right)^{2j} } \right], \end{equation} \tag{ A1 }$$

where R is the equatorial radius of the body and the coefficients, J_2j involve integrations over the mass distribution. The orbital frequency of the disk material is given by

$$\begin{equation} \Omega ^2 = \frac{1}{r}\frac{{dU}}{{dr}} = \frac{{GM}}{{r^3 }}\left[ {1 - \Sigma _{j = 1} \left({2j + 1} \right)J_{2j} \left({\frac{R}{r}} \right)^{2j} } \right] \end{equation} \tag{ A2 }$$

where the lead term gives the Keplerian value, while the other terms arise from the non-sphericity of the object. Evaluating at the radius r = R, gives an expression for the critical rotation frequency, ω_c = (GM/R³)^1/2[1 − Σ_{j = 1}(2j + 1)J_2j]^1/2. However, the values of J_2j depend on the rotation rate so the solution has to be done self-consistently. The expression for the viscous couple in the disk, g = −2πσνr³dΩ/dr would also be altered for non-Keplerian rotation, modifying the implied surface densities somewhat. These changes are strongest in the vicinity of the planet, but are minor corrections in the satellite zone. On the other hand, the mass exchange rates between planet and disk and the strength of the couple at the interface, g_p, should still be comparable to those used in the text because the critical angular momentum of the planet continues to scale with M(GMR)^1/2. Furthermore, if there is a magnetic coupling with the disk (see Section 6.3) the rotation may actually be sub-critical anyway. Since we are not interested in the planet per se, but rather in its role of providing an inner boundary condition for the disk, we chose to not include rotational flattening in our analysis. Nevertheless, it is interesting to estimate how much the critical frequency might differ from (GM/R³)^1/2.

To estimate this, we can employ a rotating polytrope. The equation of hydrodynamic equilibrium reads

$$\begin{equation} \nabla \left[ {\varphi + \eta - \frac{{\omega ^2 }}{2}r^2 \sin ^2 \theta } \right] = 0, \end{equation} \tag{ A3 }$$

where φ is the gravitational potential, η = (n + 1)Kρ^1/n is the enthalpy, and ω is the rotation rate. Taking the divergence of Equation (A3) yields

$$\begin{equation} \nabla ^2 \eta \,+\, \nabla ^2 \varphi \,-\, 2\omega ^2 = \left({n + 1} \right)K\nabla ^2 \rho ^{1/n} \,+\, 4\pi G\rho - 2\omega ^2 = 0, \end{equation} \tag{ A4 }$$

where Poisson's equation ∇²φ = 4πGρ has been used to eliminate φ. Setting ρ = ρ_cSⁿ, where ρ_c denotes the central density of the polytrope and defining a radial scaling factor ξ = r/a where a² = (n + 1)Kρ^{1/n − 1}_c/4πG, Equation (A4) can be rewritten as

$$\begin{equation} \nabla _\xi ^2 {\rm S} + {\rm S} ^n = \frac{{\omega ^2 }}{{2\pi G\rho _c }} = \alpha. \end{equation} \tag{ A5 }$$

For a non-rotating object, α → 0, there is no angular dependence and ∇²S − ξ⁻²∂/∂ξ(ξ²∂σ/∂ξ), reducing Equation (A5) to the Lane–Emden equation for the density of a polytrope of degree n. For a rotating polytrope, α remains small even near rotational instability, ∼O(GM/R³),because the central density is significantly larger than the average density $\bar \rho$, namely, $\rho _c /\bar \rho = 6.0,{\rm }11.4,{\rm }23.4$ for n = 1.5, 2, 2.5, using non-rotating models as a guide. Thus, perturbation techniques can be used to construct rotating models (e.g., Chandrasekhar 1933; Monaghan & Roxburgh 1965). It is generally found that $\alpha _c = 0.36\bar \rho /\rho _c$ is the maximum value corresponding to equatorial breakup, regardless of polytrope index, so that α_c ≈ 0.060, 0.032, 0.015 for n = 1.5,  2, 2.5. This implies a critical rotation rate of $\omega _c = \left({2\pi G\rho _c \alpha _c } \right)^{1/2} \approx 0.85\left({\pi G\bar \rho } \right)^{1/2}$ compared to $({GM/R^3 })^{1/2} = \left({4\pi G\bar \rho } \right)^{1/2} = 1.16\left({\pi G\bar \rho } \right)^{1/2}$ used in the text.

We have used the escape of material at the outer disk boundary as the primary mechanism of angular momentum loss. However, our model of the outer edge of the disk, r_d, is not well defined other than it be much further out than the centrifugal radius, r_c of the inflow. R. Nelson (2009, private communication) has suggested to us that solar tides may also contribute to the angular momentum loss. If the solar tides could remove all of the excess angular momentum, it would obviate the need for an outward flux at r_d. This would not change the disk profile too much because the edge flux is only a small fraction of the inflow anyway. However, we have envisioned the disk as not extending to R_H because as long as accretion is ongoing, incoming material that penetrates the Hill sphere in the near equatorial region interferes with the disk's expansion. Rather, we have assumed that disk material is entrained in a portion of this heliocentric flow that passes into and then exits the Hill sphere. Under these conditions, the solar tidal torque should be less important and we have not included it in this initial model. Nevertheless, models including the solar torque could be constructed and will be pursued in future work.