In pioneering work by Nelson et al. (2000), the midplane temperatures and surface densities from 2D thin-disk hydrodynamics simulations are used to fit simple vertical structures in hydrostatic equilibrium. The vertical fits then permit an estimate of the cooling rate for each column based on realistic opacities. In several simulations of Solar Nebula-sized disks, disk fragmentation, if it occurs at all, is much reduced over isothermal evolutions of the same disks. As an extension of Gammie's thin shearing box simulations, Johnson & Gammie (2003) use a similar approach by fitting what is effectively a "one-zone" Eddington-like vertical radiative equilibrium solution to every midplane cell of their 2D disks. Because their treatment is local, they are able to consider a wider range of parameter space than Nelson et al. (2000). They find that the Gammie fragmentation criterion applied to the initial cooling time is not a good measure for the likelihood of fragmentation, especially in parameter regimes where the opacity varies significantly with temperature due to evaporation of grain types. An average cooling time evaluated after the GIs become nonlinear is more reliable.
So far, there are only two treatments of radiative cooling for protoplanetary disks in full 3D for which I have sufficient information to review—those by Alan P. Boss (2001, 2002a,b, 2004a, 2005) and those by members of my own hydrodynamics group, especially
Annie C. Mejia (2004, 2005b) and Kai Cai (Cai et al. 2006). Although the work by Mejia and Cai is not yet published, radiative cooling is so important for understanding the behavior of GIs in real disks that I will unveil some interesting results that are currently in preparation. Several other research teams are working on this problem (e.g., Mayer, private communication), but I am not privy to any detailed preliminary information. Another student in my group, Aaron C. Boley, is collaborating with Ake Norlund to adapt a ray scheme for 3D radiative transfer (Heinemann et al. 2006) to disk geometry, but no disk simulation results are yet available. Because there are serious disagreements about the effects of radiative cooling between Boss and Mejia/Cai, it is worthwhile to lay out methodological similarities and differences in some detail.
4.2. Boss versus Mejia/Cai methodologies and results
Both Boss and Mejia/Cai use the same algorithm to do flux-limited diffusion in 3D within the optically thick part of the disk (Bodenheimer et al. 1990). The photosphere of the disk is determined by computing an optical depth t either inward along the spherically radial direction (Boss) or downward along the cylindrically vertical direction (Mejia/Cai). Boss and Mejia/Cai differ primarily in how they treat the optically thin regions outside the photosphere and how they match the optically thick and thin regions at the photosphere. Boss basically resets his optically thin cells either to their initial temperature or to a background temperature. The background temperature is intended to represent a thermal bath due to IR irradiation from an infalling envelope (Chiang & Goldreich 1997; Chick & Cassen 1997). He does not attempt to match fluxes at the thin/thick boundary. Mejia/Cai solve the energy equation in optically thin regions, allowing the gas to radiate freely. Without coupling to the upward flowing radiation from the optically thick disk, however, regions above the photosphere cool precipitously and do not approach radiative equilibrium. To prevent this from happening, some heating by the upward-moving radiation is included. As a boundary condition for the optically thick diffusion, the thin and thick solutions are fit together, albeit crudely, over one or two cells near the photosphere by an Eddington gray atmosphere. The Mejia/Cai algorithm takes so much care with the optically thin regions for two main reasons. First, the fractional volume of the optically thin part of the disk often exceeds that of the optically thick part, and it tends to grow as the disk cools and opacities decrease. Second, disproportionately large heating rates at higher altitudes in the disk due to shocks and irradiation can damp GIs (§4.3), and so deposition and removal of energy from these regions must be modeled.
Let us summarize significant differences between calculations by Boss and by Mejia/Cai:
1. Radiative Boundary Conditions (BCs) and the Optically Thin Regions, as detailed in the preceeding paragraph.
2. Equation of State. Mejia/Cai use a y = 5/3 monatomic gas throughout their disk, while Boss includes effects of the rotational levels and dissociation of molecular hydrogen.
3. Opacities. Boss uses Pollack et al. (1994) opacities with fixed dust grain sizes; Mejia/Cai adopt the same opacities used by Calvet et al. (1991) and D'Alessio et al. (1998), which allow variation of the maximum grain size.
4. Artificial Bulk Viscosity. Mejia/Cai always include ABV in their runs to simulate entropy production by shock waves, whereas Boss does not include ABV in the vast majority of his runs.
5. 3D Grid. The Boss code has a spherical grid; Mejia/Cai use a modified version of the Pickett et al. (2003) code with a cylindrical grid.
6. Initial Model. There are differences between the initial models used by the groups in terms of temperature T(r) and surface density distributions E(r).
7. Initial Perturbations. Mejia/Cai use a low-amplitude random density perturbation to the axisymmetric initial model; Boss gives his disk a strong hit with two, three, and four-armed modes plus a smaller random component.
Unfortunately, when we compare results from the two groups for disks with similar Md, Ms, Q, and radii, it is hard to imagine how the outcomes could be more disparate. In a series of papers over the last few years, Boss (2001, 2002a,b, 2004a, 2005) has presented simulations in which marginally unstable, radiative disks of Solar System size cool rapidly enough to fragment due to efficient convection. He finds this to be independent of metallicity and of IR irradiation by a circumstellar envelope. As Boss (2005) increases his azimuthal resolution, he tends to see the maximum density of clumps increase, which leads him to conclude that they are bound objects.
Mejia/Cai, on the other hand, generally find that, after a short adjustment period, cooling times become long, disks do not fragment, GIs' amplitudes increase with decreasing metallicity, and GIs tend to be damped by irradiation. Convection, if it occurs at all, seems to be localized and does not lead to rapid cooling. Cai, Boss, and myself have begun a collaborative effort to understand these severe differences by computing identical disks with the different methods using the same code. So far, we tentatively rule out items #2 to 7 above as the culprits. It seems probable, but not yet certain, that the primary cause of the differences is item #1, the treatment of the radiative BCs and the optically thin region. Of course, I am inclined to believe that my group's approach is the better one, but this is far from a firm conclusion at present. I will summarize some of the Mejia/Cai results (Mejia et al. 2005b; Cai et al. 2006) in the next section under the presumption that our simulations are more representative of real disk behavior, but the reader should remain respectfully skeptical until more groups present radiative hydrodynamics disk modeling in 3D.
4.3. Irradiation, metallicity, and grain size
Our first attempt to treat radiative cooling (Mejia 2004; Mejia et al. 2005b) is a simulation where the dust grains have a maximum radius of amax = 1 micron, where there is no external source of radiation shining down on the disk, and where the initial disk model is chosen to be the same as for the tcool = 2 orp simulation discussed in Section 3.1. We call this simulation "Shade," because the disk surface is not irradiated. Shade exhibits the same four phases of evolution—cooling, burst, transition, and asymptotic—as our global tcool = constant simulations, and mass transport also appears to be dominated by low-order GI modes. These characteristics are shared by all Mejia/Cai simulations, and the GIs appear to behave in a global manner.
An examination of tcool(r, t) throws considerable light on these results. Here, we define tcool in a column-wise sense, namely, tcool = the internal energy per unit projected disk area divided by the radiative flux out the top of the column. Initially, tcool is relatively low, an orp or less, and decreases with r. This tcool(r) is an artifact of our initial vertically isentropic disk structure, which is very far from radiative equilibrium. The temperature structure adjusts rapidly as the disk cools in such a way that tcool increases everywhere and becomes large, up to a few to 10 orps. It is relatively constant with r, although with very large local variations. Thus, after an initial radiative transient, Shade and all other Mejia/Cai simulations evolve toward a state best characterized as having a global tcool. This situation is drastically different from the tcool ~ r3/2 required for tcoolQ = constant. In this sense at least, GIs in real disks are a global, not local phenomenon. Typical mass inflow rates for Mejia/Cai disks in the asymptotic phase range between a few x10-7 to somewhat more than 10-6 M0/yr.
For reasons of computational cost, all Mejia/Cai calculations are run at an only moderate azimuthal resolution (128 cells from 0 to 2n), probably insufficient to permit fragmentation if it wants to occur (Pickett et al. 2003). From experience, we think we know how a disk that would fragment at high resolution behaves at low resolution, but, to be sure, we do rerun some stretches of some calculations with much higher resolution (512 cells from 0 to 2n). Although dense clumps do form during the burst of Shade, while the tcoois are still short by the Gammie fragmentation criterion, the clumps do not persist for more than a fraction of an orbit. So far, although we have not tested all cases, we have seen no evidence for significant fragmentation in any Mejia/Cai simulation with radiative cooling. The tcools in the asymptotic phase, typically of the order of five orps, are above Gammie's critical value of Prot/2 and are much higher than the initial tcools. We strongly agree with Johnson & Gammie (2003) that one cannot judge whether a disk will fragment by using initial values of tcoolQ.
Disks can be subject to irradiation by the central star (e.g., D'Alessio et al. 1999), by hot stars in a clustered environment (e.g., Johnstone et al. 1998), by a binary companion, and by a circumstellar envelope (e.g., Chiang & Goldreich 1997). The latter case is the easiest to handle, because an envelope will radiate in the far IR or at millimeter wavelengths where the opacities and optical depths are similar to those within our standard disk for r > few AU. To test the effects of envelope irradiation, Cai et al. (2006) rerun the Shade calculation with an IR flux of blackbody temperature T;rr shining down onto the disk from the grid boundary along the z-direction. Figure 7 compares the results at the same simulation time in the asymptotic phase for Tirr = 0 (Shade), 15 K, and 25 K. Irradiation suppresses the GI amplitudes. The difference between the 15 K and 25 K cases is a bit difficult to discern clearly in the figure. However, in extreme cases, for Tirr > about 40 K, the disk is not able to cool, and Q remains too large for instability to occur at all. Direct stellar irradiation is more difficult to treat, because the stellar photons are absorbed initially in a surface layer too thin for our code to resolve, but preliminary indications in Mejia (2004) are that it too tends to suppress GIs. Although other forms of irradiation have not yet been treated in GI simulations, I suggest that we will find that any external environmental factor that pumps heat into the disk will weaken if not prevent GIs. If this is correct, then comprehensive modeling of the radiative environment of real disks is required in order to determine whether GIs occur and how they behave.
Contrary to Boss (2002a, 2004a), Mejia/Cai find that metallicity matters. Cai et al. (2006) has repeated the Tirr = 15 K run with metallicities ranging from one quarter to twice the solar value. Equatorial plane density images from the two extreme cases are shown in the left and center panels of Figure 8. What the eye suggests when comparing Figure 7 (center panel) and Figure 8 (left and center panel) is confirmed by objective measures of GI amplitudes based on Fourier analyses, namely, GIs become stronger as the metallicity is lowered. Cai et al. (2006) do agree with Boss in one unexpected way. The optically thick parts of the disks converge to roughly similar behaviors and lose energy at similar rates regardless of metallicity. However, there is no efficient convection, and the cooling is slow, with tcool typically on the order of five to seven orps, except for localized regions which are now being tested for possible fragmentation at higher resolution. Major differences are found in the optically thin regions, which dominate the volume of the disk at low metallicity. If gas giant planets form directly, all at once, by f
Figure 8. Equatorial plane density grayscales in the asymptotic state for three calculations with Tirr = 15 K and varied metallicity and grain size. Left. One-quarter solar metallicity for amax = 1 micron. Center. Twice solar metallicity for amax = 1 micron. Right. Solar metallicity and amax = 1 mm. Images and simulations are similar in all other respects to those in Figure 7. Figures provided courtesy of K. Cai.
disk instability, then having stronger GIs at lower metallicity seems to run counter to the observed tendency for more planets to be found around high metallicity stars (Fischer & Valenti 2005). Anyhow, the tcools generally seem too long for fragmentation to occur.
Figure 8 shows that, at least in this disk model, use of a larger maximum grain size, amax = 1 mm, increases the opacity over the bulk of this rather cool disk, resulting in an effect similar to an increase in metallicity. Modeling of real disks does seem to require larger grain sizes (D'Alessio et al. 2001), even at early times (Osorio et al. 2003). The opacity variations are complex, however. In the hotter inner disk, increases in grain size and settling of grains should reduce rather than increase the opacity. Also, as shown by Nelson et al. (2000), as grains are transported by GIs through shocks and across temperature gradients, they may be vaporized. The coupling through opacity between GIs and grain growth, destruction, and settling will prove to be a serious challenge for future modeling. We will return to aspects of this in Section 5.
A great deal remains to be done before we claim to understand GIs in radiative disks with realistic opacities. Results available to date are almost as divergent as they possibly could be. Better treatments of both the optically thick and thin regions and their interface are required. If you believe the results of my own hydro group's calculations, then GIs are sensitive to their radiation environment and to variations in opacity. For the mass (Md/Ms = 0.14) and size (few to 40 AU) of disk considered, cooling is, by and large, too slow under realistic conditions to produce fragmentation, and GIs can be weakened or even suppressed by external radiation fields impinging on disks. Similar conclusions have been reached recently by others through analytic arguments (Rafikov 2005; Metzner & Levin 2005). The reader must remember, however, that the most extensive body of published work on this subject, by Boss, argues for fast cooling by convection regardless of metallicity or envelope irradiation.
Was this article helpful?