Testing Disk Instability Models for Giant Planet Formation
Testing Disk Instability Models for Giant Planet Formation
Alan P. Boss
Department of Terrestrial Magnetism, Carnegie Institution of Washington, 5241 Broad
Branch Road, NW, Washington, DC 20015-1305
boss@dtm.ciw.edu
ABSTRACT
Disk instability is an attractive yet controversial means for the rapid formation
of giant planets in our solar system and elsewhere. Recent concerns regarding
the first adiabatic exponent of molecular hydrogen gas are addressed and shown
not to lead to spurious clump formation in the author’s disk instability models.
A number of disk instability models have been calculated in order to further test
the robustness of the mechanism, exploring the effects of changing the pressure
equation of state, the vertical temperature profile, and other parameters affecting
the temperature distribution. Possible reasons for differences in results obtained
by other workers are discussed. Disk instability remains as a plausible formation
mechanism for giant planets.
Subject headings: solar system: formation – planetary systems
1. Introduction
The disk instability mechanism for giant planet formation is based on the formation
and survival of self-gravitating clumps of gas and dust in a marginally gravitationally un-
stable protoplanetary disk (Boss 1997; reviewed by Durisen et al. 2007). In order for a
disk instability to succeed, the disk must be able to cool its midplane as the clumps form,
allowing them to continue to contract to higher densities, and the clumps must be able to
survive indefinitely in the face of Keplerian shear, tidal forces, and internal thermal pressure.
Considering the complicated physical processes involved during the time evolution of a three
dimensional disk instability, it is perhaps not surprising that the theoretical basis for the
disk instability hypothesis remains unclear even after a decade of work on the subject.
The Indiana University (IU) group has been active in studying disk instabilities, and has
generally found that disk instabilities are unable to lead to the formation of self-gravitating,
http://arxiv.org/abs/0704.1138v1
– 2 –
dense clumps that could go on to form gas giant protoplanets (e.g., Pickett et al. 2000;
Mej́ıa 2004; Cai et al. 2006a,b; Boley et al. 2006, 2007a,b). On the other hand, the
Washington-Zurich group (e.g., Mayer et al. 2002, 2004, 2007) has presented models that
support the hypothesis that disk instabilities can lead to the formation of long-lived giant
gaseous protoplanets. We present here several new calculations that attempt to understand
the reasons for some of these different outcomes regarding disk instability.
2. Energy Equation of State
Boley et al. (2007a) pointed out that uncertainties about the ortho/para ratio of molec-
ular hydrogen at low temperatures might lead to differences in the outcome of disk insta-
bility models. They suggested that a mixture intermediate between pure parahydrogen and
a 3:1 ortho/para ratio might be the most appropriate choice. Boley et al. (2007a) also
noted that discontinuities in the specific internal energy equation for molecular hydrogen
could lead to artificially low values of the first adiabatic exponent for a simple perfect gas
(Γ1 = γ = 1 + Rg/(µcV ), where Rg is the gas constant, µ is the mean molecular weight,
and cV is the specific heat at constant volume; Cox & Giuli 1968). Values of γ ≤ 4/3 can
lead to dynamical instabilities (either expansion or contraction) away from a configuration
of hydrostatic equilibrium (Cox & Giuli 1968). Boley et al. (2007a) suggested that γ ≤ 4/3
could artificially lead to clump formation.
Boley et al. (2007a) stated that the energy equation of state (EOS) used in disk instabil-
ity models by Boss (2001, 2002b, 2005) contained a discontinuity that might be responsible
for artificially lowering γ below 4/3, leading to spurious fragmentation. Boley et al. (2007a)
based this assertion on the equations of state described by Boss (1984). However, all of the
hydrodynamical models run by the present author since Boss (1989) have been based on a
different energy equation than that reported by Boss (1984), as a result of a direct com-
parison with the equation of state routines employed by Werner Tscharnuter. The models
of Boss (1989) revised the energy equation treatment to include an interpolation between
temperatures of 100 K and 200 K, but this revision was not explicitly stated in the Boss
(1989) paper as it seemed insignificant at the time.
Using the notation of Boss (1984), the specific internal energy for molecular hydrogen
for temperatures T < 100K is taken to be E∗H2 = 3/2 RgT/µ, while for temperatures between
100K and 600K it is E∗H2 = 5/2 RgT/µ. For intermediate temperatures, 100K < T < 200K,
the internal energy is interpolated according to the equation E∗H2 = 3/2 RgT/µ [1+2/3 (T −
100)/100]. The Erratum by Boley et al. (2007b) used this revised EOS.
– 3 –
Boley et al. (2007a,b) showed that for a pure parahydrogen or 3:1 ortho:para mix, γ
decreases significantly at ∼ 100K, but does not fall below 4/3. Figure 1 depicts the behavior
of γ calculated with the revised Boss equation of state and shows that γ drops below 4/3 for
135K < T < 200K. The densest clumps found in Boss’s disk instability models generally
have maximum temperatures below 135K, and even lower mean temperatures. Boss (2005)
presented a disk instability model with the highest spatial resolution computed to date, and
found that the densest clump that formed had a maximum temperature of 120K and a mean
temperature of 94K. Boss (2006c) found that clumps formed in disk instabilities around M
dwarf protostars had maximum temperatures less than 100K. Boss’s (2002a) Table 1 showed
maximum clump temperatures of 115K to 126K for a range of models with varied opacities.
Perhaps the most important point raised by Boley et al. (2007a, b) is that γ is likely to
decrease significantly around 100K, and this softening of the pressure EOS will enhance the
formation of dense clumps, as happens in locally isothermal disks with an effective γ = 1. In
fact, Figure 1 in Boley et al. (2007b) shows that the Boss γ is higher than that of either of
the preferred hydrogen mixtures for T < 100K, implying that clump formation is suppressed
somewhat in the Boss models as a result.
3. Pressure Equation of State
Several new models have explored using the same pressure EOS as is used in the IU
group models (e.g., Cai et al. 2006a,b; Boley et al. 2006): the gas pressure p is given by
(γ − 1)ρE, where γ = 5/3, ρ is the gas density, and E is the specific internal energy of
the gas. Three models were run with this pressure EOS, starting from the same initial disk
model as model HR in Boss (2001) – a 0.091M⊙ disk orbiting a 1M⊙ protostar with an
outer disk radius of 20 AU. The calculations were made with the same three dimensional,
gravitational, radiative hydrodynamics code as in Boss (2001) and in all subsequent Boss
disk instability calculations (see these earlier papers for more details about the calculational
techniques and initial conditions). The models had Nφ = 256 and NY lm = 32, though,
compared to Nφ = 512 and NY lm = 48 for model HR (Boss 2001).
The three new models varied the choice for the critical disk density (ρcr = 10
−13, 10−12,
or 10−11 g cm−3) below which the disk temperature was forced back to its initial value
(typically 40K in the outer disk). This artifice was employed in Boss (2001) and subsequent
models in order to maintain a reasonably large time step when low density regions develop
in the disk that are undergoing decompressional cooling.
In all three models, spiral arms and transient clumps form within 200 yrs of evolution,
– 4 –
similar to the behavior with the usual (Boss 1984) pressure EOS (e.g., Boss 2001). Relatively
high density clumps form, with maximum densities similar to those in a calculation with the
same spatial resolution but the usual pressure EOS. Evidently using the γ = 5/3 EOS
does not alter the results in a significant way because the Boss EOS also has γ = 5/3 for
T < 100K, the regime where clumps form. The choice of ρcr makes little difference as well, as
was found by Boss (2006b) for disk instability models in binary star systems. These results
suggest that the reason for differing outcomes must be sought elsewhere (see Discussion).
4. Varied Disk Temperature Parameters
We now present a set of four models varying several of the parameters that could affect
the temperature distribution in the disk models. Table 1 summarizes the four models, which
are all variations on model HR of Boss (2001). However, these models all had the same
spatial resolution (Nφ = 512) and number of terms in the spherical harmonic expansion for
the gravitational potential (NY lm = 48) as model HR in Boss (2001). The parameters that
were varied included: the temperature of the thermal bath (i.e., the envelope temperature
Te), the critical density in the disk below which the temperature was reset to the initial
disk temperature at that orbital radius (as in the previous section), the critical density
in the envelope below which the gas was assumed to be at a temperature equal to the
envelope temperature Te (for grid points at least 8 degrees above the disk midplane), and
finally, whether the temperature was forced to decline monotonically (Boss 2002a) with
vertical height inside the disk (mono) or not (free). The former constraint errs on the side of
artificially cooling the disk by removing local temperature maxima in the vertical direction.
All models started at a time of 322 yrs of evolution in model HR, and continued for at least
another 17 yrs of evolution (∼ 1 clump orbital period).
The results of these variations on the standard assumptions are shown in Figures 2 and
3. The variation that produced the largest deviation from the standard assumptions (model
H in Figure 2) was relaxing the constraint on the monotonic vertical (more precisely, in the
θ angle) decline of temperatures within the disk (model TZ in Figure 3), though even this
model led to the formation of well-defined clumps that were no more than a factor of two less
dense than in model H. Models T and TE led to evolutions that were very similar to that of
model H and so are not shown. The models show that these three variations in the details of
how the disk thermodynamics is treated in the Boss models are not particularly significant
for the outcome of a disk instability, presumably because of the thermal bath assumption.
– 5 –
5. Discussion
There are a number of possible reasons for different outcomes compared to other groups:
Spatial resolution – Clump formation is strongly enhanced as the spatial resolution
in the critical azimuthal direction is increased from Nφ = 64 to 512 (Boss 2000). Boss
(2005) presented a model with Nφ = 1024 and a locally refined radial grid (equivalent to
a calculation with over 8 × 106 grid points) that implied that in the continuum limit, the
outcome of a disk instability is the formation of dense, self-gravitating clumps. Cai et al.
(2006a,b) calculated models with Nφ = 128, increasing Nφ to 512 for two models only
after those models had entered a phase of evolution when nonaxisymmetry was no longer
growing. Boley et al. (2006) similarly calculated models with Nφ = 128, increasing Nφ for
some models to 512 for the earliest phase of evolution, leading to the formation of dense
clumps at the intersections of spiral structures. The clumps disappeared in a fraction of an
orbital period. Clumps typically last no more than an orbital period in even the highest
spatial resolution models of Boss (2005). Boss (2005) thus used virtual protoplanets to allow
the orbital evolution of these dense clumps to be followed further than is possible with even
a high spatial resolution calculation with a fixed Eulerian grid code.
Gravitational potential solver – The Boss models use a spherical harmonic (Ylm(θ, φ))
expansion to solve Poisson’s equation for the gravitational potential, with the accuracy of the
resulting gravitational potential being strongly dependent on the number of terms (NY lm)
carried along in the expansion. As we have seen, model HR in Boss (2001) used NY lm = 48.
Boss (2000, 2001) found that increasing NY lm led to the formation of significantly denser
clumps. Boss (2005) further explored the effects of using an enhanced gravitational potential
solver by replacing some of the mass in the densest regions of a clump with a point mass at
the center of the relevant grid cell, finding that this led to even better defined, higher density
clumps. In comparison, the IU group uses a direct solution of Poisson’s equation, with a
boundary potential employing terms up to l = m = 10 (Pickett et al. 2000), implying a
limited ability to depict small-scale gravitational forces in a strongly nonaxisymmetric disk.
Mej́ıa (2004) considered Fourier analysis of her disk models for m ≤ 6, while Boley et al.
(2006) considered m ≤ 63, finding increasingly little power in modes with m > 10, possibly
consistent with the cutoff at l = m = 10 in their boundary potential.
Artificial viscosity – Pickett & Durisen (2007) found that the inclusion of certain artificial
viscosity (AV) terms could enhance the survival of clumps formed in a disk instability,
and suggested that AV could thus explain the long-lived clumps found in SPH calculations
by Mayer et al. (2002, 2004, 2007). Pickett & Durisen (2007) further noted that even
calculations without any explicit AV (e.g., the Boss models) should be considered suspect,
given the intrinsic numerical viscosity of any numerical code. Boss (2006b) showed that with
– 6 –
a large amount of explicit AV, clump formation is suppressed, as was also found by Pickett &
Durisen (2007). However, the level of intrinsic numerical viscosity in Boss code models with
Nφ = 256 appears to be equivalent to an α-viscosity of ∼ 10
−4 or smaller (Boss 2004), a level
that appears to be negligible in comparison to typical explicit AV levels. Considering that
the virtual protoplanet (VP) models of Boss (2005) had Nφ = 256, the continued survival of
the VP for at least 30 orbital periods in these models is not likely to have been affected by
the intrinsic numerical viscosity of the Boss code.
Stellar irradiation – Mej́ıa (2004) considered the effects of stellar irradiation on the
surface of the disk as a means of heating the disk surface and thereby possibly suppressing
clump formation. The Boss models assume that the disk is immersed in a thermal bath
appropriate for backscattering from infalling dust grains in the protostellar envelope, with
an envelope temperature appropriate for a protostar that is not undergoing an FU Orionis
outburst (Chick & Cassen 1997). While the Boss models thus do not include the effects of
direct irradiation by the central protostar, the dynamical evolution of a three dimensional
disk leads to strongly variable vertical oscillations and structures (Mej́ıa 2004; Boley &
Durisen 2006; Jang-Condell & Boss 2007) that are not considered in simple theoretical
models of flared accretion disks. The highly corrugated inner disk surfaces (stretching at
least 29 degrees above the disk midplane; Jang-Condell & Boss 2007) will shield the outer
disk surfaces from the central protostar, eliminating this source of heating for much of the
outer disk.
Radiative transfer – Boley et al. (2006) have calculated disk instability models using the
flux-limited diffusion approximation (FLDA) along with a detailed treatment of the transi-
tion from the optically thick disk to the optically thin atmosphere of the disk. They suggest
using a plane-parallel (one dimensional) atmosphere as a test case. Boss (2001) investigated
the effects of using the FLDA instead of the standard diffusion approximation (DA) coupled
with a thermal bath for low optical depths, but did not find any significant differences. My-
hill & Boss (1993) showed the results of the fully three dimensional, standard nonisothermal
test case for protostellar collapse, calculated with their two independent Eddington approx-
imation (EA) codes, finding good agreement. Whitehouse & Bate (2006) found that their
FLDA models of the standard nonisothermal collapse problem led to temperature profiles
similar to those found by Myhill & Boss (1993), though with appreciably hotter gas tem-
peratures where the optical depth was ∼ 2/3. They attribute this difference to the FLDA
retarding the loss of radiation in these layers compared to the EA. Given that the choice of
the flux limiter can have an effect on the outcome (Bodenheimer et al. 1990), it is unclear
whether any particular implementation of the FLDA is superior to the EA. The standard
Boss models have used the DA coupled with a thermal bath to force the DA models to mimic
an EA calculation.
– 7 –
Numerical heating – In calculations by the IU group, nonaxisymmetric perturbations
tend to grow rapidly for a certain period of time and then begin to damp out (e.g., Cai et
al. 2006a,b). In the Boley et al. (2006) calculations, the disk starts out with a mass of 0.07
M⊙ and a radius of 40 AU, but then expands outward to a radius of ∼ 80 AU, leading to
the formation of rings inside 20 AU and a gravitationally stable region outside 20AU with
spiral arms that do not fragment. The latter behavior is roughly consistent with the models
by Boss (2003), who studied disks extending from 10 AU to 30 AU, and found fragments
to form at 20 AU but not at 30 AU. Similarly, Boss (2006a) studied disks extending from
100 AU to 200 AU, and found no evidence for fragmentation. Thus on scales larger than
∼ 20 AU, the results of Boley et al. (2006) and Boss (2003, 2006a) are in basic agreement.
The disagreement arises about what happens in the inner disks. Fragmentation typically
occurs at 8 to 10 AU in the Boss models, whereas the inner disk rings do not fragment in the
Boley et al. (2006) models. Boley et al. (2006) state that their inner disk is stable to ring
fragmentation because of “numerical heating” at distances out to ∼ 7 AU. This non-physical
heating appears to have affected the models of Cai et al. (2006a,b) as well as those of Boley
et al. (2006), making the IU results for inner disks difficult to accept.
6. Conclusions
While there are a number of potential reasons for the differences in disk instability
models calculated by the IU group and the present author, at this time the major sources
of these differences would appear to be some combination of several effects, namely spatial
resolution, gravitational potential solver accuracy, and numerical heating in the inner disk
of the IU models. Handling of the boundary conditions for the disk’s radiative energy losses
is another possibility that is still under investigation by the author, though the models
presented here suggest that this may not be a dominant effect. Given the current state of
knowledge, and the new results presented herein, it appears that reports of the death of the
disk instability model for giant planet formation have been greatly exaggerated.
I thank Aaron Boley, Kai Cai, and Megan Pickett for working with me to understand the
differences between their group’s results and my results, and the referees for their remarks.
This research was supported in part by NASA Planetary Geology and Geophysics grant
NNG05GH30G and is contributed in part to NASA Astrobiology Institute grant NCC2-1056.
The calculations were performed on the Carnegie Alpha Cluster, the purchase of which was
partially supported by NSF Major Research Instrumentation grant MRI-9976645.
– 8 –
REFERENCES
Bodenheimer, P., Yorke, H. W., Róz̀yczka, M., & Tohline, J. E. 1990, ApJ, 355, 651
Boley, A. C., & Durisen, R. H. 2006, ApJ, 641, 534
Boley, A. C., et al. 2006, ApJ, 651, 517
Boley, A. C., Hartquist, T. W., Durisen, R. H., & Michael, S. 2007a, ApJL, in press
——. 2007b, ApJL (Erratum), submitted
Boss, A. P. 1984, ApJ, 277, 768
——. 1989, ApJ, 346, 336
——. 1997, Science, 276, 1836
——. 2000, ApJL, 536, L101
——. 2001, ApJ, 563, 367
——. 2002a, ApJ, 567, L149
——. 2002b, ApJ, 576, 462
——. 2003, ApJ, 599, 577
——. 2004, ApJ, 616, 1265
——. 2005, ApJ, 629, 535
——. 2006a, ApJL, 637, L137
——. 2006b, ApJ, 641, 1148
——. 2006c, ApJL, 643, 501
Cai, K., et al. 2006a, ApJL, 636, L149
——. 2006b, ApJL (Erratum), 642, L173
Chick, K. M., & Cassen, P. 1997, ApJ, 477, 398
Cox, J. P., & Giuli, R. T. 1968, Principles of Stellar Structure (New York, Gordon and
Breach)
Durisen, R. H., et al. 2007, in Protostars and Planets V, B. Reipurth, D. Jewitt, & K. Keil,
eds. (Tucson, Univ. Arizona Press), p. 607
Jang-Condell, H., & Boss, A. P. 2007, ApJL, in press
Mayer, L, Quinn, T., Wadsley, J., & Stadel, J. 2002, Science, 298, 1756
——. 2004, ApJ, 609, 1045
– 9 –
——. 2007, astro-ph/0606361
Mej́ıa, A. C. 2004, PhD thesis, Indiana Univ.
Myhill, E. A., & Boss, A. P. 1993, ApJS, 89, 345
Pickett, B. K., Cassen, P., Durisen, R. H., & Link, R. 2000, ApJ, 529, 1034
Pickett, M. K., & Durisen, R. H. 2007, ApJL, 654, L155
Whitehouse, S. C., & Bate, M. R. 2006, MNRAS, 367, 32
This preprint was prepared with the AAS LATEX macros v5.2.
http://arxiv.org/abs/astro-ph/0606361
– 10 –
Table 1. Models with varied disk temperature parameters.
model Te (K) ρcr (disk) ρcr (envelope) T (θ)
H 50K 10−11 10−11 mono
T 100K 10−13 10−11 mono
TZ 50K 10−13 10−11 free
TE 50K 10−13 10−13 mono
– 11 –
1 1.5 2 2.5 3
Fig. 1.— Adiabatic exponent Γ1 = γ used in Boss (2001) and in all subsequent Boss disk
instability models (solid line). The short-dashed line shows γ = 4/3, while the long-dashed
line shows γ = 5/3, as used by Cai et al. (2006a,b) and Boley et al. (2006).
– 12 –
-1 -0.5 0 0.5 1
Fig. 2.— Equatorial density contours for model H after 339 yrs of evolution. The disk has
an outer radius of 20 AU and an inner radius of 4 AU. Hashed regions denote clumps and
spiral arms with densities higher than 10−10 g cm−3. Density contours represent factors of
two change in density.
– 13 –
-1 -0.5 0 0.5 1
Fig. 3.— Same as Figure 2, but for model TZ.
	Introduction
	Energy Equation of State
	Pressure Equation of State
	Varied Disk Temperature Parameters
	Discussion
	Conclusions
