EVOLUTION OF COATED GRAINS IN SPIRAL SHOCKS OF SELF-GRAVITATING PROTOPLANETARY DISKS

Protostellar disks likely undergo one or more phases of gravitational instability during the first 10⁵–10⁶ years of their life (Vorobyov & Basu 2010a, 2010b; Hayfield et al. 2010). Two factors determining their susceptibility to gravitational instability, as well as the ultimate outcome, are the mass loading from the molecular cloud envelope (Boley 2009; Rafikov 2009; Boley et al. 2010) and how fast the gas is able to cool and dissipate the energy generated by compression and shocks in the spiral arms (Rice et al. 2005; Boley et al. 2007; Mayer et al. 2007; Cossins et al. 2010; Meru & Bate 2010) either by radiative cooling alone or by a combination of radiative cooling and convective cooling (Boss 2003; Durisen et al. 2007; Mayer 2010). Whether gravitational instability results in a self-regulated, marginally unstable state with a long-lived spiral structure, or whether self-regulation is broken and fragmentation takes place within spiral arms, producing sub-stellar objects on brown dwarf or giant planet scales (Stamatellos & Whitworth 2009; Boley et al. 2010; Mayer 2010), depends on the ability to attain a gas cooling timescale comparable to, or shorter than, the local orbital timescale. Both radiative and convective cooling ultimately depend on the opacity of the gas. Therefore it is of pivotal importance to investigate in detail how the opacity might change within spiral shocks.

At the relatively low temperatures of protoplanetary disks reached at ten or more AU from the star the major source of opacity is dust grains, and different prescriptions are used to calculate their contribution. Some of these prescriptions are normally used in the three-dimensional hydrodynamical simulations of gravitationally unstable protostellar disks (D'Alessio et al. 2001) but it is unclear whether they account for all the important effects on the grains as they undergo rapid temperature and density variations within spiral shocks.

In this work we explore in detail the effect of varying background temperature and gas density on the grain size distribution and its resultant opacity. In particular, we consider the time-dependent density and temperature evolution in strong spiral shocks appearing in gravitationally unstable protoplanetary disks that attain a minimum Toomre Q parameter just above the value needed for fragmentation (the necessary condition for fragmentation is that the minimum Q be 1.4 or lower—see, e.g., Durisen et al. 2007). We extract such information from one of the three-dimensional smoothed particle hydrodynamics (SPH) simulations presented in Mayer et al. (2007), which included radiative transfer with a flux-limited diffusion approximation.

The paper is organized as follows: first we describe the model adopted for the three-dimensional spiral shocks and for dust grains, then we describe our results on the evolution of dust grain opacity in our reference dust model, analyzing different locations throughout the spiral shock. We show results for the disk midplane first and then for regions well above the midplane including the effect of irradiation from the star. Furthermore, we investigate the effect of changing parameters in our model for dust, such as the dust-to-gas ratio in the spiral shocks, the water content of grains, and the distribution of grain sizes. Subsequently, we explore a model in which dust grains are composite. We conclude the paper with a discussion of the implications of our results for disk cooling and fragmentation into sub-stellar objects.

The thermodynamics of spiral shocks adopted throughout the paper is extracted from the simulation shown in the bottom left corner of Figure 1 of Mayer et al. (2007), which led to a long-lasting strong spiral structure (the minimum Toomre Q parameter was hovering around ∼1.5 for many orbital times) but not to fragmentation (see Mayer et al. 2007 for details on the parameters of this particular run). Fragmentation did not occur because the gas was unable to radiate on a sufficiently short timescale.

The profiles are computed in thick annuli centered on the density maximum of one of the spiral shocks (which also roughly corresponds to the temperature maximum). In the figure, the circle marks the radial distance at which all annuli are centered. However, in order to avoid diluting the effect of the spiral arm, the profiles are not computed by azimuthally averaging over a 2π annulus (i.e., over the whole disk); rather, they are azimuthally averaged across a three-dimensional rectangular box with an initial height of 0.75 AU and an initial azimuthal extent of 4 AU, as shown in gray. At any time, the box is chosen in such a way as to encompass the same particles within the chosen section of the spiral shock, exploiting the Lagrangian nature of SPH. Therefore, the box changes in size as time progresses because particles diffuse in space as the shock widens (see Figure 1). The different annuli have increasing radial widths up to the edges of the box, so that they probe conditions increasingly further away from the density maximum, and, in order to simplify the comparison, the width of each annulus is kept constant while the box changes in size (we verified that we are able to follow most of the particles initially identified within the box).

Zoom In Zoom Out Reset image size

We purposely consider a vigorously unstable but non-fragmenting disk case because this should represent a case in which a sudden variation in the cooling timescale, eventually triggered by opacity changes, might rapidly bring the disk from a self-regulated state to rapid fragmentation (see Meru & Bate 2010; Cossins et al. 2010). At the location of the spiral shock, the disk reaches the minimum Toomre Q of ∼1.5 (about 12 AU from the center). In the simulation, the disk was growing steadily in mass over hundreds of orbital times, a crude attempt to model accretion from the surrounding molecular cloud core. At the chosen stage of disk evolution, the disk has a total mass of about ∼0.12 M_☉ out to 30 AU, i.e., slightly above 10% of the mass of the central star (which has a mass of 1 M_☉).

The SPH simulations of Mayer et al. (2007) used the opacity model of D'Alessio et al. (2001), which was, in turn, based on the model of Pollack et al. (1994). Pollack et al. (1994) considered both spherical grains and spheroids where the axial ratio was 4:1. They also considered details of the chemical composition of the different components. They investigated the case where each component is segregated in a separate class of grains, as well as the case where the grains were composites of several species, produced by the coagulation of smaller grains. A possibility that was not considered, however, was the case where the grains are composed of a refractory core surrounded by an ice shell. Such a situation might occur if the grains were in a spiral shock where the heating would evaporate the water ice.

As the gas cools the water vapor would recondense on the silicate nuclei to form an ice mantle. In such a case, we might expect that the change of opacity with temperature would be different from the changes found by D'Alessio et al. (2001) and Pollack et al. (1994). This could have important consequences for the stability of the disk models. We investigate this case in the present paper.

To model this process, we consider a distribution of silicate cores. We assume that the number of cores per unit of volume with radii between a and a+da is given by a power law:

$$\begin{equation} n(a)da=a^{-\xi }da, \end{equation} \tag{ 1 }$$

where ξ is a constant. The standard case we examine will have a ratio of silicates to hydrogen/helium of 4.23 × 10⁻³ by mass. We also assume, following D'Alessio et al. (2001), that there is water present and that the mass ratio of water to silicates is 2.13.

We use a temperature–density relationship in the gas which is relevant for spiral shocks in self-gravitating protoplanetary disks, derived for the case shown in Figure 1, and follow the evolution of the grain opacity. Unless otherwise specified, the water is assumed to be in vapor form initially, even though, for most of the cases considered, this vapor will be supersaturated by many orders of magnitude. We consider only inhomogeneous nucleation, that is, that the vapor will only condense onto the silicate cores. Once the silicate cores are introduced, the condensation of water vapor is very fast. Except for a very short time near the beginning of the simulation, the difference between the case where all the water is initially in the form of vapor and the case where all the water is in the form of an ice shell around the silicate core is negligible.

Marginally unstable protoplanetary disks, which do not undergo fragmentation but rather self-regulate to a minimum Toomre Q of ∼1.5–1.6 exhibit a nearly stationary spiral pattern (Boley & Durisen 2008). In the Mayer et al. (2007) simulations, the dominant spiral mode was m = 3 over several tens of orbital times. Spiral structure is indeed long-lived in a marginally unstable, self-regulated state.

Zoom In Zoom Out Reset image size

Locally, the gas density and temperature undergo maxima and minima over time as trailing density waves propagate through the disk, with a timescale of order the local orbital time, and the gas shocks correspond to the density jumps. For our analysis, we choose to describe a time window also of order the orbital time at the radius where shocks normally reach their maximum strength, which is of order 12–15 AU (at such distances the orbital time is 30–45 yr). While we have chosen a particular phase of the disk evolution for computing the shock profiles, this choice is by no means special because, as we mentioned, the spiral pattern is nearly stationary in disks that are in a self-regulated state.

We utilize the structural data for one of the three spiral shocks in particular (see Figure 1) after having verified that the relevant conditions, namely temperature and density evolution, are very similar among the three. We note that disks with a more complicated, higher order spiral pattern might exhibit more differences between shocks as these might appear over a wider range of distances in the disk where both background density and temperature as well as shear can vary substantially (compare, for example, the case of lighter, colder disks considered in Mayer et al. 2004).

The temperature–density profiles of the shock as a function of time are shown in Figure 2 for three positions through the spiral arm at the disk midplane; distances are computed with respect to the density maximum along the spiral arm (we average between the conditions at either side of the spiral arm for simplicity since no remarkable temperature or density differences are seen).

For each of these cases, we have computed the evolution of the temperature of a grain as a function of time. At any instant the grain is heated by absorbing ambient radiation from the surroundings and by latent heat released from the condensation of water vapor onto the grain. It is cooled by re-radiation and by water ice sublimation. Contact with the gas is also included and acts to either cool or heat the grain depending on the sign of the temperature difference between them. Details of the energy balance calculation are given in Mekler & Podolak (1994), Podolak & Mekler (1997), and Podolak & Zucker (2004). In computing the cross section for absorption and radiation, we used the FORTRAN code DMILAY.F developed by W. Wiscombe of the Goddard Space Flight Center.⁴ It is based on the theory presented in Toon & Ackerman (1981). This code computes the scattering and absorption cross sections for a spherical grain composed of a core and a shell.

The grain size distribution in our baseline case has ξ = 3.5, a minimum radius of a_min = 10⁻⁶ cm, and a maximum radius of a_max = 0.1 cm as in D'Alessio et al. (2001). This distribution can be well represented with five size bins. The bins are chosen to divide the logarithm of the size range equally. The grain size associated with each bin is that of the (logarithmic) center of the bin, and the number of grains in the bin is taken to be the total mass of grains in the size range represented by the bin divided by the mass of the representative grain for that bin. In our case, the silicate core sizes are 2.21 × 10⁻⁶, 2.21 × 10⁻⁵, 2.21 × 10⁻⁴, 2.21 × 10⁻³, and 2.21 × 10⁻² cm, respectively.

When the silicate cores are initially inserted into the gas, there is a rapid condensation of the vapor onto the grains. The background gas is highly supersaturated, so the water condensation flux onto the grains is much higher than the water sublimation flux off the grains. The net growth rate of the grain is nearly independent of the grain temperature and the grain radius. As a result, all the grains acquire an ice shell of the same thickness. The behavior of the grain opacity (Rosseland mean) as a function of time at different parts of the density wave is shown in Figure 3.

Zoom In Zoom Out Reset image size

The ice shell forms within minutes and the increased radius of the grains causes the grain opacity to rise from around 0.5 to around 3 cm² g⁻¹. The saturation of the water vapor in the gas phase drops quickly to 1. At this point the grains are all nearly at the same temperature as the background gas, although a slight difference in temperature between the different grain sizes does exist. Typically the smallest grains are very near the temperature of the gas, while the larger grains are a fraction of a degree cooler than the gas. The difference is usually of the order of 10 mK although it can be as much as 0.3 K when rapid condensation of water vapor onto the grains heats them. A convenient visualization of this is given in Figure 4 where we show the energy flux transferred from the grain to the gas as a function of time for a particular choice of the shock profile (black dashed curve). Negative values correspond to colder grains being heated by the gas. The warmer grains sublimate their ice at a slightly higher rate than the background gas returns it, while the opposite occurs with the cooler grains. As a result, this temperature difference between different grain sizes causes the ice to slowly migrate from warmer grains to cooler ones while the background gas remains 100% saturated. The rate of migration is moderated by the number density of water molecules in the background gas, and will therefore be a sensitive function of the temperature of the gas.

Zoom In Zoom Out Reset image size

For a power-law size distribution with ξ = 3.5, most of the surface area is in the smallest grains. At the relevant temperatures (∼100 K), the Planck peak is at λ ∼ 30 μm. At these wavelengths, the size parameter of Mie theory for the smallest grains is

$$\begin{eqnarray*} &&x=\frac{2\pi a}{\lambda }=2\times 10^{-3}. \end{eqnarray*}$$

For such small values of x, the cross section of the grain for both scattering and absorption is much less than its geometric cross section. Thus despite the relatively large fraction of mass they contain, these small grains contribute very little to the grain opacity, even with their ice shells. As the ice migrates to the larger, more efficient scatterers, it contributes more efficiently to the cross section, and the opacity increases.

Eventually, the gas temperature in the density wave becomes high enough that the gas becomes undersaturated, and the ice shells quickly sublimate from all the grains. When the gas cools after reaching its maximum temperature, there is again a rapid condensation of the water vapor onto the silicate cores, and again, the initial thickness of the ice shell is independent of the grain radius. Since the gas temperature is still quite high, the ice quickly migrates to the coolest grains, and as the gas cools the ice shell thickness remains unchanged. The opacity still decreases with decreasing temperature, but this is again due to the temperature dependence of the Rosseland mean and not to changes in the grain size distribution.

The thickness of the ice shell for the grain distribution described above is shown in Figure 5 for the case of 0.2 AU in the density wave. The shell thicknesses all begin at the same value, but the ice migrates to the larger, cooler grains. Note that it is not the largest grain in the distribution that ends up with the thickest shell. This is fortunate because the largest grains contribute relatively little to the effective surface area of the distribution. Sequestering most of the ice on these grains would not allow it to contribute significantly to the opacity, and the variation of opacity with temperature would be much more moderate.

Zoom In Zoom Out Reset image size

This behavior leads to an interesting hysteresis effect. If a grain distribution develops ice shells at a low enough temperature, the ice will migrate very slowly or not at all. If the shell is formed at a higher temperature and then cooled, the ice will have the opportunity to migrate to cooler grains before being "frozen" in place. The point is that the opacity of a grain distribution is not only a function of temperature at a given time, but also of the thermal history of the gas and grains. The opacity behavior with time depends on the details of this thermal history. Below we present four examples from different regions of the density wave in more detail. The hysteresis effect is shown for these four cases in Figure 6 where we show the opacity as a function of temperature. Shown for comparison are the grain opacities for these particular conditions as taken from Pollack et al. (1985) (dotted curve in the upper left panel) and D'Alessio et al. (2001) (dashed curves). As can be seen, the agreement with published opacities is quite good before the heating event, but can differ by a factor of two or more after the heating event.

Zoom In Zoom Out Reset image size

3.1. 0.2 AU

3.2. 0.54 AU

3.3. 0.72 AU

3.4. 1.32 AU

At 1.32 AU in the density wave, heating is weak enough that the temperature never rises above 160 K. In this case, the entire region where the ice shell is completely evaporated is avoided. Instead, the opacity rises monotonically to its peak value as in the first two cases and then drops slowly as the temperature decreases again (see lower right panel in Figure 3). In the lower right panel of Figure 6, we see that the hysteresis effect in the opacity as a function of temperature is very similar to the first two cases except that now the drop at higher temperatures is absent since these temperatures are never reached.

Zoom In Zoom Out Reset image size

It is important to note that this hysteresis effect only occurs the first time that a shock wave encounters the grains. Once the grain distribution has experienced such a heating event, the ice is redistributed preferentially on to the colder grains. Future heating events may remove the ice shells, but as the gas cools they will reform on the colder grains as before. This is illustrated in Figure 7 where we show the opacity as a function of temperature for the 0.2 AU case, but where the grains have encountered a second shock immediately after the first one. The green triangles show the same case as Figure 6; the red squares show the opacity–temperature relation for the second shock. The profile is unchanged. The blue diamonds show the opacity–temperature relation for the case where the grains started with an ice-to-olivine ratio of 2:1 by mass. Since the ice is distributed differently among the grains in this case, the opacity originally differs as well. Once the shock passes, the ice rearranges itself into the same distribution independent of its original configuration.

In addition to the standard case, we have run several non-standard scenarios to test the dependence of the opacity on other parameters of the model. These additional cases were run for a point at 0.2 AU in the density wave. We summarize the results below.

4.1. High-Z Case: Dust Concentration in Spiral Arms

Haghighipour & Boss (2003b, 2003a) have studied the motion of solids in a disk. As Mayer et al. (2007) point out, shocks can cause solids to concentrate and thus raise the abundance of high-Z material with respect to hydrogen. We have run a case with the same size distribution as the standard case, but with the rock and water increased by a factor of 10 over the solar value. The result is shown in Figure 8. The panel on the upper left shows the evolution of the opacity as a function of time for this case. Although the number of silicate cores increases by an order of magnitude, the ratio of water to silicates remains the same. As a result, we expect that the opacity at any given time will simply be 10 times the opacity of the standard case. A comparison with Figure 3 shows that this is indeed the case. Small differences can probably be attributed to the fact that since the background water vapor can be 10 times higher than in the standard case, there will be some differences in the rates of sublimation and condensation of the ice shell, and this will lead to small differences in opacity.

Zoom In Zoom Out Reset image size

4.2. High Water Content in Grains

4.3. Different Power Law for Dust Size Distribution

4.4. Low Opacity Regions and the Effect of Stellar Irradiation

All the cases described above were run for the case of high opacity because we have considered shock profiles at the midplane of the disk. In this case, the radiative heating is due to the surrounding gas. Since this gas is at a temperature of the order of 100 K, the radiation it emits will be mostly in the infrared. Well above the disk plane, however, the opacity drops, and the visible radiation from the star can penetrate into the disk. Because the radiation from the star is concentrated at much shorter wavelengths, the size parameter of the grains with respect to this wavelength increases, and their efficiency for scattering and absorption changes. In addition, water ice is essentially transparent in the visible so that the albedo of the grains changes as well. The magnitude of any changes due to this effect is reduced by the fact that these density waves only appear at 10 AU or more from the star.

The total optical depth in the inter-arm region is of the order of 100, so the region where the optical depth is of order unity comprises only a few percent of the total mass in a column. In the arm itself, the total optical depth is a factor of 10 higher, so the low optical depth region is an even smaller fraction of the total mass. Nonetheless, it is of interest to see how the evolutionary behavior of the ice grains is affected by the addition of stellar radiation.

Figure 9 shows the opacity as a function of temperature at 0.2 AU in the arm (the same case as shown in the upper left panel of Figure 6). Here, however, we have added the stellar radiation due to a star at a distance of 10 AU. The parameters of the star (radius and surface temperature) are the same as for our Sun. The blue curve shows the case where the gas density is the same as the density for the optically thick case shown in the upper left panel of Figure 6 so that the only additional effect is the added stellar radiation. The grains in the first bin have a silicate core of radius 0.022 μm. At low temperatures with an ice shell they grow to a radius of ∼0.2 μm but most of this additional radius is composed of nearly transparent ice. As a result, their extinction cross section is small even in the visible, and the additional stellar heating does not appreciably affect their temperature. The grains in the larger bins absorb a larger percentage of the stellar heating, however, and reach higher temperatures. As a result, the ice migrates to the smallest grains rather than to the larger ones as in the high opacity case. At IR wavelengths that are relevant to the Rosseland mean, the smallest grains do not contribute much to the opacity, and the opacity is about half what it is in the high optical depth case.

Zoom In Zoom Out Reset image size

Since the optical depth depends on the gas density, we would expect that this density is much lower in the low optical depth region. The red curve in Figure 9 for a gas density of 10⁻⁴ of the midplane is therefore more appropriate. Here, because of low gas density, cooling by contact with the gas is inefficient. In this case, the grains reach much higher temperatures than the gas, and the coolest grains are those in bin 4, with core radii of 22 μm. Even these grains are too hot to retain ice shells until the gas temperature drops to below around 125 K and the opacity starts to rise. Even so, the opacity is still much lower than the corresponding value at midplane.

So far we have considered the grains to consist of a silicate core with an ice mantle, but this is likely a very simplified model of the actual structure of grains. In particular, suppose the grains consist of silicate micrograins embedded in an ice matrix. In this case, the effect on the opacity can be much more dramatic.

Consider first a case similar to that investigated above. If the grains are assumed to be composed of a mixture of rock and ice, and we assume that the ice-to-rock ratio, by mass, is 2:1, then M_i = 2M_r, where M_i is the ice mass and M_r is the rock mass. This means that the mass fraction of ice is X_i = 2/3 and the mass fraction of rock is X_r = 1/3. If we take the bulk density of ice to be ρ_i = 1 g cm⁻³ and the bulk density of the silicate to be ρ_r = 3 g cm⁻³, the density of the mixture is

$$\begin{eqnarray*} &&\frac{1}{\rho }=\frac{X_i}{\rho _i}+\frac{X_r}{\rho _r} \end{eqnarray*}$$

and ρ = 1.29 g cm⁻³. For a given mass of grain material, M, the number of grains, N is given by

$$\begin{equation} M=M_i+M_r=\frac{4\pi }{3}Na^3\rho \end{equation} \tag{ 2 }$$

and

$$\begin{equation} N=\frac{3M}{4\pi a^3\rho }. \end{equation} \tag{ 3 }$$

To a zeroth approximation, we can neglect the differences in the refractive index between water and rock, and assume that they all have a refractive index of m = 1.5 + 0.1i. Then

$$\begin{eqnarray*} &&\frac{m^2-1}{m^2+2}=0.296+0.05i \end{eqnarray*}$$

and, if the size parameter of the grain is given by

$$\begin{eqnarray*} &&x=\frac{2\pi }{\lambda }, \end{eqnarray*}$$

where λ is the wavelength of the incident radiation, the efficiency factor for extinction for x ⩽ 1 can be well approximated by (van de Hulst 1957)

$$\begin{equation} Q_e=Q_0x^4=x^4\Re \left[\frac{8}{3}\left(\frac{m^2-1}{m^2+2}\right)^2\right], \end{equation} \tag{ 4 }$$

where the symbol ℜ refers to the real part of the expression in the brackets. For our case, this gives Q_e = 0.0851x⁴. For a given mass of grains, therefore, the opacity will be proportional to

$$\begin{equation} N\pi a^2Q_e=\frac{3MQ_e}{4a\rho }=Q_0\frac{3M}{4a\rho }x^4. \end{equation} \tag{ 5 }$$

If the ice now evaporates, the new grain mass will be M* = X_rM. The new grain bulk density will be ρ* = ρ_r, and the new grain radius will be

$$\begin{eqnarray*} &&a^*=a\left(\frac{X_r\rho }{\rho _r}\right)^{1/3}. \end{eqnarray*}$$

If the original opacity is

$$\begin{eqnarray*} &&\kappa =\kappa _0NQ_e\pi a^2 \end{eqnarray*}$$

the new opacity will be

$$\begin{equation} \kappa ^*=\kappa _0Q_0\frac{3M^*x^{*4}}{4a^*\rho ^*}=\kappa \left(\frac{X_r\rho }{\rho _r}\right)^2. \end{equation} \tag{ 6 }$$

For the parameters we have chosen,

$$\begin{eqnarray*} &&\left(\frac{X_r\rho }{\rho _r}\right)\approx 0.15 \end{eqnarray*}$$

so the opacity should drop to around 0.02 of its original value. In fact, in the full model we see a smaller change, around 0.1, but this is partially due to the fact that we are using the same refractive index for water, rock, and the water + rock mixture. Since ice is more transparent than rock, and in this simple model is assumed to have the same extinction properties as rock, the removal of ice will overestimate the change in opacity. More important, however, is the fact that the higher opacity occurs at a lower temperature, where the peak of the Planck function is at longer wavelengths. In going from 130 K to 220 K the wavelength decreases by a factor of 0.6, so x increases by the same factor, and x⁴ increases by a factor of around eight, so the overall change is the product of these two factors ≈0.16, which, allowing for the previously mentioned refractive index effect, is about right. So, to zeroth order, the opacity change for grains with x ⩽ 1 is

$$\begin{equation} \kappa ^*=\kappa \left(\frac{X_r\rho }{\rho _r}\right)^2\left(\frac{T^*}{T}\right)^4, \end{equation} \tag{ 7 }$$

where T* is the high temperature at the place where the ice has evaporated and T is the low temperature at the place where the ice is in the grain.

Now suppose that the grains are actually composed of a conglomerate of micrograins held together in an ice matrix. As the ice evaporates, the micrograins are released. Suppose that each large grain of radius a contains micrograins of radius a'. Again the ratio of ice to rock is 2:1. Now when the ice disappears, the micrograins are released and the relevant parameters change as follows: M* = X_rM, and ρ* = ρ_r as before, but now a' = a and

$$\begin{eqnarray*} &&x^*=x\frac{a^{\prime }}{a}, \end{eqnarray*}$$

$$\begin{eqnarray*} &&N^*=\left(\frac{a}{a^{\prime }}\right)^3\frac{X_r\rho }{\rho _r}N. \end{eqnarray*}$$

As a result

$$\begin{equation} \kappa ^*=\kappa \frac{X_r\rho }{\rho _r}\left(\frac{a^{\prime }}{a}\right)^3\left(\frac{T^*}{T}\right)^4. \end{equation} \tag{ 8 }$$

Note that when a' = a* for the simple ice shell evaporation, Equation (8) reduces to Equation (7) as it should. If a ∼ 1 μm and a' ∼ 0.1 μm, for example, the reduction in opacity could be very dramatic.

Using the time-dependent spiral shock profiles from a three-dimensional hydrodynamical simulation of a gravitationally unstable disk, we have calculated the effects that such shocks have on the opacity of gas. Our results show that opacity variations of about an order of magnitude can take place on timescales of several months to years. Such timescales are much shorter than the orbital time, which is larger than 10 years at the distances of the shocks in the simulations (10–20 AU from the star). This means that the cooling time in the disk will change on timescales much shorter than the orbital time. In particular, when the opacity drops by nearly an order of magnitude on such short timescales it means that the cooling time will also decrease proportionally fast on such short timescales (see also Cossins et al. 2010). The cooling disk may not be able to self-regulate on such short timescales and therefore the likelihood of fragmentation should be substantially increased. Likewise, the temperature increase caused by the shock provides the right conditions to produce a rapid increase in molecular weight within the spiral shock, supporting the claim of Mayer et al. (2007), as long as rapid migration of solids toward the arms can occur simultaneously and increase the dust-to-gas ratio.

The water vapor that is released by the grains will increase the mean molecular weight of the gas. As Mayer et al. (2007) point out, this will make the disk more unstable. If the molecular weight of the background gas is μ_x and this gas is a mass fraction X of the total, then introducing a gas of molecular weight μ_z and mass fraction 1 − X will increase the mean molecular weight, $\overline{\mu }$, to

$$\begin{equation} \overline{\mu }=\mu _x\left[1+\frac{\mu _z-\mu _x}{\mu _z}\left(\frac{1-X}{X}\right)\right]. \end{equation} \tag{ 9 }$$

For solar composition gas, μ_x ≈ 2.4, while μ_z = 18 for water vapor. For the baseline model we investigated, the water vapor mass fraction is 9 × 10⁻³ so that $\overline{\mu }$ should increase by only 1%. But if grains are concentrated by the spiral arms and the mass fraction of water vapor increases by an order of magnitude, $\overline{\mu }$ can increase from 2.4 to 2.6. This is large enough to have noticeable effects on the stability of the disk (Mayer et al. 2007).

No published hydrodynamical calculations exist that can show the simultaneous effect of both variations. By combining the results of Mayer et al. (2007) and Meru & Bate (2010), one can expect that the formation of giant planets in the inner disk would be much more likely with a concurrent rapid drop in opacity and an increase in molecular weight, probably irrespective of whether convection occurs or not (see also Johnson & Gammie 2003).

Ruden & Pollack (1991) have pointed out that if the opacity, κ, has a temperature dependence of the form

$$\begin{eqnarray*} &&\kappa \sim T^{\beta } \end{eqnarray*}$$

then a gas with an adiabatic exponent of 1.45–1.5 (which is a good approximation in the temperature range of the shock at peak strength, between 150 and 240 K, for solar metallicity, see Boley et al. 2007) will satisfy the Schwarzschild criterion for convective instability if β > 0.5–1. The rising parts of the dashed curves in Figure 6 have slopes corresponding to β = 0.5. As can be seen from this figure, when the shock ramps up, in the region with T = 150 K where the ice quickly migrates to the largest cores, the value of β is temporarily >1, namely it is above the critical value expected at this temperature, for which γ = 1.5 according to Boley et al. (2007). The exact value of β depends on the location within the shock (i.e., it varies across different annuli because the opacity–temperature relation varies), but it is always high enough to allow, in principle, the development of convection. Later on, as the temperature increases and then decreases again, the adiabatic index also varies, thus changing the critical threshold, and β becomes sub-critical. The same is true for a second passage of the shock since the opacity has also changed in the meantime. Therefore, it appears that convection could be triggered initially but would be short-lived. However, since convection itself could change the density and temperature distribution through the shocking region, only a dynamical simulation that includes our coated grain model will be able to ascertain the thermodynamical evolution of the gas phase. Finally, in this discussion we have not considered how the conditions might change as a function of the altitude from the midplane since we have used azimuthally averaged quantities, but going beyond this will once again require dynamical (three-dimensional) simulations since convection, if it starts, is expected to change the physical conditions, particularly across the vertical extent of the disk.

Another possible effect could arise because molecular weight variations can, in principle, damp the convective instability. Indeed, since vertical convection requires ds/dz < 0, where s is the specific entropy, and, ignoring adiabatic index gradients, ds/dz = dln T/dz − (γ − 1)dlnρ/dz − dlnμ/dz for an ideal gas, one sees that a vertical gradient in the molecular weight can have a stabilizing influence on convection. Therefore, it will be important to understand whether or not significant μ gradients would develop once μ can be self-consistently calculated in simulations including both gas and dust. The important point to remember is that, for convection to be relevant for disk fragmentation, what matters is that it affects the cooling of the disk at the midplane, and not necessarily across the entire vertical extent of the disk. Therefore, as long as molecular weight gradients are negligible near the disk midplane, convection will not be affected. In practice, one will need to show that there is a large enough region around the midplane in which heat can be redistributed rapidly by convection but within which the molecular weight is increased nearly uniformly following dust grain collection and sublimation through the shock.

It is tempting to speculate that the gas becomes convectively unstable at this point, and this convection allows the gas to cool more quickly. If this is indeed the case, then the increase in opacity would be offset by convective cooling. Note that this high β occurs only during the heating of the gas, and β is negative during the cooling stage.

Certainly the influence of ice coatings on grains needs to be investigated directly with numerical simulations since these effects are not necessarily additive; for example, one can imagine that an increasing molecular weight would affect the temperature evolution. In addition, the opacity variation across the shock would also be affected if the timescales of the two processes are even slightly uneven.

In general, the present work shows the importance of modeling grain chemistry within the spiral shocks of self-gravitating protoplanetary disks. While at 50 AU < R < 200 AU fragmentation in disks loaded by high accretion rates from the molecular envelope is an almost unavoidable event that depends little on the details of the gas and dust chemistry because the cooling time is expected to be short (Boley 2009; Boley et al. 2010; Rafikov 2009; Vorobyov & Basu 2010a), fragmentation in the inner disk (<30 AU) depends very strongly on thermodynamics and therefore on the properties and response of dust grains that affect the cooling time via the opacity. Other thermodynamical factors, such as the effective adiabatic index (Boley & Durisen 2008) and the molecular weight (Mayer et al. 2007), are also bound to play an important role because they change the temperature evolution across the shock.

The models that we describe in this paper are not dynamical; the grains do not move in response to the gas, nor do they coagulate and change their size distribution. Recent calculations (Boley & Durisen 2010) show that grains significantly above interstellar size can migrate toward spiral arms, increasing the dust-to-gas ratio and even affecting fragmentation itself. In our calculations, we have assumed a size distribution appropriate for interstellar grains and no grain growth (other than the growth of an ice shell) so that the assumption of coupling between gas and dust is satisfied. In the spiral shock, the dust-to-gas ratio might go up by a factor of 10 if most grains have already grown to centimeter grains. This could open a new avenue for forming planetesimals rapidly via gravitational instability of the dust layer (Rice et al. 2006). Dust accumulation in spiral shocks might also be an additional effect promoting fragmentation in the gas phase (Boley & Durisen 2010), concurrent with those that we have just described.

The results of our paper show that the detailed modeling of the response of the dust grains to spiral shocks has an important effect on the evolution of the opacity, one of the crucial parameters controlling the development and outcome of gravitational instability in massive disks. In their present form our results should be regarded as a "proof of concept" since the temperature evolution of the shocks that we have assumed in the paper has not been recomputed using our dust model, but rather stems from the (D'Alessio et al. 2001) model adopted in the simulations. Directly incorporating the dust model that we have described in this paper in our three-dimensional hydrodynamical simulations of disks will be the first and most important next step. The stronger and faster opacity changes that we have found might produce fragmentation in a disk model that does not fragment with the standard opacity models used in Mayer et al. (2007) as well as in other published work. Likewise, since grain growth might have already taken place by the time the disk becomes spirally unstable, we will need to extend our model to larger grains with sizes comparable to those that can accumulate efficiently within spiral shocks.