# Formation of Jupiter using opacities based on detailed grain physics

###### Abstract

Numerical simulations, based on the core-nucleated accretion model, are presented for the formation of Jupiter at 5.2 AU in 3 primordial disks with three different assumed values of the surface density of solid particles. The grain opacities in the envelope of the protoplanet are computed using a detailed model that includes settling and coagulation of grains and that incorporates a recalculation of the grain size distribution at each point in time and space. We generally find lower opacities than the 2% of interstellar values used in previous calculations [Hubickyj, O., Bodenheimer, P., Lissauer, J. J., 2005. Icarus 179, 415–431; Lissauer, J. J., Hubickyj, O., D’Angelo, G., Bodenheimer, P., 2009. Icarus 199, 338-350]. These lower opacities result in more rapid heat loss from and more rapid contraction of the protoplanetary envelope. For a given surface density of solids, the new calculations result in a substantial speedup in formation time as compared with those previous calculations. Formation times are calculated to be , , and Myr, and solid core masses are found to be , , and , for solid surface densities, , of , , and , respectively. For and , respectively, these formation times are reduced by more than 50% and more than 80% compared with those in a previously published calculation with the old approximation to the opacity.

###### keywords:

Planetary formation, Jovian planets, Jupiter: interior, Accretion## 1 Introduction

The core-nucleated accretion theory of giant planet formation was effectively originated by Safronov (1972), who realized that a major part of the process involved the accumulation of solid bodies. Consider the accretion of relatively small objects (planetesimals) onto a significantly larger embryo. If is the mass of the embryo, then the fundamental equation for its growth, in the absence of gas, is

(1) |

where is the effective geometrical capture cross-section, is the surface density of solid material (planetesimals), is the escape velocity from the embryo, is the orbital frequency, and is the relative velocity of embryo and planetesimal. The quantity in brackets is the gravitational enhancement factor , here given in the (2 + 2)-body approximation. Safronov mentioned the later capture of gas to form giant planets, and Cameron (1973) stated that Jupiter could form by gas accretion once a solid core of about 10 M had been accreted. Perri and Cameron (1974) then constructed an adiabatic, hydrostatic equilibrium model of a protoplanet with a core and a gaseous envelope, extending outward to the Hill radius. They were able to show that there was a critical mass for the core, above which the envelope was unstable, implying that rapid gas accretion could occur. However, the value of this critical core mass was at least 70 M, much higher than the deduced core masses of Jupiter and Saturn. The assumption of adiabaticity turns out to be incorrect (Stevenson, 1982), as more detailed calculations show that radiative transport with subadiabatic gradients is able to carry the energy outwards in significant portions of a growing giant planet.

Mizuno (1980) constructed models of a protoplanet including both radiative and convective energy transport in the envelope, along with an energy source provided by the planetesimals accreting onto the core at a constant rate, typically M/yr. The model assumed strict hydrostatic and thermal equilibrium, and the outer radius was assumed to be the protoplanet’s Hill radius. Again a critical core mass was found, above which the envelope was unstable, and the value was about 12 M, close to the estimated core masses of the giant planets. Furthermore, the critical value was practically independent of the distance to the Sun, partly as a result of the assumption of the same constant value of at all distances.

The assumption of strict hydrostatic equilibrium was relaxed by Bodenheimer and Pollack (1986), who allowed quasi-static contraction and calculated the first evolutionary sequences applied to the core-nucleated accretion model. The planetesimal accretion rate was assumed to be constant, also at M/yr. They found that gravitational contraction (but not dynamical collapse) began to become important at about the point where core and envelope masses were equal, at a value which later came to be known as the crossover mass. Rapid gas accretion also began near the crossover mass, which they found to be about 16.5 M.

The next major step was a series of calculations by Pollack et al. (1996) who allowed a variable core accretion rate given by eq. (1). The simulations contained three major elements: (1) three-body accretion cross-sections of solids onto the embryo (Greenzweig and Lissauer, 1992) to determine , (2) a modified stellar evolution code (Henyey et al., 1964) to follow the structure and evolution of the gaseous envelope, and (3) a calculation of the trajectories of planetesimals through the gaseous envelope (Podolak et al., 1988), to deduce the deposition of solid mass and energy in the envelope and to determine the protoplanet’s effective capture radius, taking into account gas drag. The often-quoted result of these calculations is that it took Jupiter 8 Myr to form in a solar nebula with g cm, three times the density in the minimum mass solar nebula, under the assumption that the opacity due to small grains in the cool outer envelope was provided by grains with an interstellar size distribution and with interstellar abundances with respect to the gas. It was also shown that if the grain opacities were arbitrarily reduced to 2% of interstellar values, the formation time was reduced to about 3 Myr.

Subsequent observations of protoplanetary disks and modelling of the interior of Jupiter revealed two major problems with the Pollack et al. (1996) baseline-case results for the formation of Jupiter. First, that case gave a formation time longer than the mean lifetime of protostellar disks, and, second, it produced a core mass of about 16 M, above the plausible range at the time of (0–11) M, derived from matching interior models of the present Jupiter with observations (Saumon and Guillot, 2004). Further calculations (Hubickyj et al., 2005) considered, with generally improved physics, the effects of the following parameters: (1) the opacity arising from small grains in the envelope, (2) the solid surface density in the protoplanetary disk, and (3) a possible cutoff in core accretion rate as a result of competition for solid particles by neighboring embryos. The reduction in opacity was justified by calculations (Podolak, 2003) showing that grains entering into the protoplanetary envelope would rapidly coagulate, settle, and eventually evaporate, resulting in an actual opacity well below interstellar values. In the actual protoplanetary formation calculation the grain opacity was simply reduced to 2% of interstellar values, with the result that the formation time was reduced to under 3 Myr, half that found for interstellar values of grain opacity. The final core mass, however, still remained near 16 M. A reduction in from 10 to 6 g cm, again with the opacity at 2% of interstellar, increased the formation time to 13 Myr and reduced the core mass to 8 M. It turned out that the constraints of short formation time and small core mass could be satisfied simultaneously only if the accretion rate of solids into the protoplanet were arbitrarily cut off at some time, simulating roughly the effect of planetesimal accretion by neighboring protoplanets. A cutoff at M led to a formation time of 4.5 Myr, while a cutoff at 10 M gave a formation time of 1 Myr.

In the planet formation simulations of Pollack et al. (1996) and Hubickyj et al. (2005), simplified surface boundary conditions were applied during the rapid gas accretion phase. First, the gas density in the protoplanetary disk was assumed to be constant in time. Second, the gas accretion rate was capped at M/yr, a rough estimate of the rate at which the disk could supply gas to the protoplanet as a result of its viscous evolution. Third, the outer boundary of the protoplanet was set at the modified accretion radius , defined by Bodenheimer et al. (2000) to be

(2) |

where is the Hill sphere radius, is the sound speed in the disk, is the planet’s mass, and are constants, both set to unity. Equation (2) is an interpolation between the Bondi accretion radius and the Hill sphere radius, such that when is large, reduces to the Bondi radius, and when is small, reduces to the Hill radius.

These assumptions were modified by Lissauer et al. (2009) to take into account the results of new three-dimensional numerical simulations of protoplanetary disks accreting onto protoplanets. These simulations employed the code developed by D’Angelo et al. (2003). The results from the 3-D simulations showed that , that is, the planet can collect material from the disk only within a quarter of its Hill radius, at best. The gas accretion rates from the 3-D simulations, which included the effects of gap opening, were used whenever they were less than the accretion rate needed to maintain the boundary of the planet at ; thus the arbitrary cap on this rate was removed. The parameters investigated in this set of simulations included (1) the size of the region from which the protoplanet can accrete, (2) the viscosity in the protoplanetary disk, which has a strong influence on the disk accretion rate onto the protoplanet once the envelope has shrunk and the planet is massive enough to begin clearing gas from nearby its orbit, and (3) the time-dependence of the gas density in the protoplanetary disk in the vicinity of the protoplanet. The grain opacity was set to 2% interstellar and g cm. The results of these improved simulations showed that Jupiter’s formation time was still less than 3 Myr, and that a relatively low-viscosity disk () is consistent with the formation of planets of about 1 Jovian mass, while higher viscosity produced planets that were more massive than most observed extrasolar giant planets.

In all of the simulations just mentioned, as well as in detailed formation simulations by other groups (Alibert et al., 2005; Fortier et al., 2007; Ikoma et al., 2000), the opacities arising from grains in the protostellar envelope have been approximated either by interstellar grain values (Pollack et al., 1985) or by a reduction of those values by some arbitrary amount, typically a factor 50. However, the simulations of Podolak (2003) and Movshovitz and Podolak (2008), in which grain coagulation and sedimentation were calculated for a few particular protoplanetary envelopes, showed that the ratio of actual grain opacity to interstellar grain opacity was highly variable, depending on the depth in the envelope. Near the surface of the planet the opacities were close to interstellar, while deep in the envelope they were reduced by a factor of up to 1000. Those opacity calculations, however, were not coupled to the evolution of the protoplanet, and they included only the outer radiative layer in the atmosphere, which typically has a temperature of 600 K at its inner boundary.

In the present paper we extend the basic grain settling and coagulation model of Movshovitz and Podolak (2008) by including the entire depth of the envelope down to a temperature of 1800 K, where practically all grains have evaporated. Since some layers of the outer part of the envelope are convective, we also modify the grain simulation to take into account convective effects. We then couple the grain calculation to the evolution of the protoplanet, recalculating the grain size distribution and opacity in every layer of the protoplanet at every time step of the evolutionary calculation. These significant improvements in the core accretion model allow us to investigate the following questions: (1) How is the formation time for a Jupiter-mass planet at 5.2 AU modified compared with the earlier approximate calculations regarding the opacity? (2) Most of our previous simulations for the formation of Jupiter required an initial solid surface density of 10 g cm in the protoplanetary disk. How far can this assumed density be reduced such that Jupiter still forms within a dissipation time of the disk – roughly 3 Myr – and what is the resulting solid core mass? The following sections describe the nature of the grain simulations, the basic idea of the core-nucleated accretion model, how the two calculations are integrated, and our new results for four formation runs, up to the onset of rapid gas accretion, for a Jupiter-mass planet at 5.2 AU.

## 2 Grain Physics

This section describes, for a protoplanetary envelope with a given distribution of temperature, density, and convective velocity as a function of distance to the center, how grain settling and coagulation modify the grain size distribution and opacity at each level. The size distribution of the grains as a function of height is determined as follows: At each atmospheric level the number of grains per unit volume, , evolves according to

(3) |

where the derivatives are Eulerian. Here is the mass of the grain, is the distance from the center of the planet, and is the time.

The first term on the RHS of Eq. (3) is a source term, and refers to those grains that are carried into the envelope either by the accreted gas or by planetesimal ablation. We do not consider the effects of evaporation and condensation on the grain size distribution. A grain is instantaneously evaporated when its temperature reaches 1800 K. (We assume that the grains are composed of silicates. The error introduced by this assumption is probably small compared with the overall uncertainty in the physics of grain structure and agglomeration. This point is furthered discussed below.)

The second term on the RHS of Eq. (3) refers to the evolution of the grain size distribution by coagulation. We treat this by solving the Smoluchowski equation

(4) |

The coagulation kernel, is essentially the probability of a collision between grains of mass and that results in the two grains sticking together to form a new grain. The details of how is computed are given in Podolak (2003) and Movshovitz and Podolak (2008).

In practice, the grain size distribution is approximated by a finite number of size bins, where the th bin represents a grain of radius

and is a scaling constant, set in this case to 1 m. It has been shown (Movshovitz and Podolak, 2008) that if the minimum grain size in a protoplanetary envelope is reduced to 0.1 m, the results are practically identical to those for a minimum size of 1 m, because of rapid coagulation of the small particles. The number of bins is chosen so that the size range of interest can be adequately represented. As a result, the integrals in Eq. (4) are replaced by sums.

The third term on the RHS of Eq. (3) is a transport term. In the absence of convection, the transport is given by

(5) |

Here is the sedimentation velocity of the grain, and it is a function of the grain size and the local values of gravity, gas density, and temperature (Movshovitz and Podolak, 2008).

If convection is present, it affects the results in two ways.
First, small grains that are entrapped by a convective eddy reach
higher speeds than normal, while large grains are less affected by
these eddies. This effect causes the relative velocities between small and large
grains, and the resultant coagulation kernels, to change.
Weidenschilling (1984) gives expressions for grain relative velocities in
the limits where both grains are large, one is large and one is small,
and both small (his Eqs. 15, 16, and 18, respectively). More recently,
Ormel and Cuzzi (2007) give an expression for the relative velocity between
particles of arbitrary size (their Eq. 16), that matches the
Weidenschilling expressions in the relevant limits. We use
the Ormel and Cuzzi expression, modified to include the effects of the
systematic (sedimentation) velocity of the grains. The sedimentation velocity
shifts the boundary between small (coupled to the gas) and large (uncoupled)
grains, and the *difference* in sedimentation velocity between two
grains is also added in quadrature to the relative velocity due to turbulence.

In practice, however, the enhancement of collision probability due to turbulence, although fully included in the calculations, is not terribly important. For the most part, relative velocities between grains of different sizes are determined by the difference in sedimentation speed, and the relative velocities between small particles of similar size are determined only by their random thermal motion, because they are coupled to the gas. The turbulence does add collisions between large grains of equal size. These grains are too large for significant thermal motions, and with no difference in sedimentation speed they have no appreciable relative velocities, except those induced by turbulence. Because this enhancement is only apparent for large grains ( cm), the effect on the size distribution and opacity is negligible.

A second effect of convection is that the transport term becomes more complicated. Podolak (2003) treated convective transport with an eddy diffusion approximation. This leads to numerical difficulties and requires very short time steps. As a result, it is unsuitable for use in the current evolutionary models. We have therefore treated convective transport with a simplified algorithm that retains the major characteristics of convective mixing, explicitly conserves mass, and is very computationally efficient. For each layer in the envelope and for each grain size, we decide if convection is important or not based on whether two conditions are satisfied. The first is that the time it takes for a grain to become entrapped in the largest eddy be shorter than the lifetime of that eddy. If the grain is small enough so that the time to entrap it is short, the grain responds to the convective motion. For this calculation, we have taken the length of the largest eddy to be 1.5 times the pressure scale height. The eddy lifetime was taken to be this length divided by the speed of the eddy. The second condition that is relevant to convective transport is that the convective speed of the eddy be greater than the sedimentation speed of the grain. If this condition is not met, the grain sediments faster than the convection can mix it. Going back to the first condition, note that there are actually two requirements for a grain to be carried by an eddy. First, as stated above, the stopping time must be shorter than the eddy turnover time. Second, the stopping time must be shorter than the grain crossing time (eddy size/grain velocity). In our case, the crossing time is , and the turnover time is , where is the eddy size. Since the second condition requires , the crossing time is longer than the turnover time, which we have already required to be longer than the stopping time. Thus the requirement (stopping time crossing time) is automatically met.

For each layer in the envelope and for each bin size we determine whether or not these two criteria are met. If they are, convection will be the dominant mode of transport. We then determine, for each grain size, which regions (i.e., sets of contiguous layers) of the envelope will mix that grain size. In those regions we determine the total number of grains, and redistribute them so that the ratio of mass density of grains to mass density of gas is constant. Grain sedimentation velocities range from millimeters to meters per second. We find that almost all grains are mixed for convective velocities larger than cm s.

This algorithm has the disadvantage that it allows for either complete convective mixing or no convective mixing. Partial convective mixing is not considered. However, it has the advantage that it explicitly conserves the grain mass in the region, and it allows much larger time steps than the eddy diffusion algorithm. Numerical tests show that the resultant grain size distributions behave in accordance with what is expected for convective mixing, so this seems to be a useful approximation to convective transport.

Once the grain size distribution is established, the opacity is computed using an approximation to Mie scattering suggested by Cuzzi (personal communication; Movshovitz and Podolak, 2008). The effective cross-sections for absorption and scattering are calculated for each grain size and frequency, and combined to produce the opacity . This function is then integrated over the grain size distribution and over the blackbody spectrum to produce the Rosseland mean opacity for each level in the atmosphere.

We have assumed that the flux of grains into the uppermost layer of the envelope is 1% of the gas accretion mass flux, given a near-solar composition in the protoplanetary disk. This will cause us to overestimate the opacity for the following reason. It assumes that all of the condensed material in this region is in the form of small grains. This is certainly not true, since many of the grains will have accreted into larger planetesimals. These planetesimals are eventually accumulated by the planet, but their grains are released into the envelope at much deeper levels and thus have much less of an influence on the opacity of the outer layers.

Our model also assumes that all of the grains are composed of silicates. In fact, they are likely to be composed of a combination of silicates, organic material and water, in roughly equal fractions by mass. This has two major consequences. First, although the silicates and organics will remain solid at the temperatures relevant to the outer envelope, the water will evaporate. This means that we have overestimated the mass of grain material at all levels hotter than K and the true grain opacity will be correspondingly lower.

A second effect is that the bulk density of the grains in these upper layers will be lower than we have assumed, and this means that the grains will sediment more slowly and their number density will be higher than we have computed. As a result there should be some increase in the opacity. Below we estimate the magnitude of this effect.

In the uppermost layers of the envelope, there is very little grain growth and most of the grains are the same size (Movshovitz and Podolak, 2008). To first approximation, the time for such a grain to settle out of a layer of thickness is

(6) |

where is the sedimentation speed. The time between grain collisions is roughly

(7) |

Here is the mean free path between grains, and is the thermal speed of a grain.

If is short compared to , the grains will collide and grow before they sediment out of the region. The grain size can be estimated by setting . The mean free path between grain collisions will be

where is the number density of grains, and is the grain radius. The Knudsen number in this region is large so that the sedimentation speed can be written as

where is the acceleration of gravity, is the bulk density of a single grain, is the mass density of the gas, and is the thermal speed of a gas molecule. The mass flux through the layer is denoted by , where

(8) |

The thermal speed of a molecule (or grain) of mass is given by:

where is the temperature and is Boltzmann’s constant. Setting and using the above relations, we find

(9) |

The cross-section for scattering or absorption is the geometric cross-section, , multiplied by some efficiency factor, , that must be computed from Mie theory. Typically, for grains much smaller than the wavelength of the impinging radiation, and roughly equal to one for larger grains. Since we are interested in Rosseland mean opacities at temperatures of around 100 K, the wavelengths will be in the micron range and the grains are in the 1 micron range in the upper levels of the envelope. The corresponding size parameter is

where is the wavelength of the radiation.

For small the extinction is mostly due to absorption, and a useful approximation for is (Movshovitz and Podolak, 2008):

where and are the real and imaginary indexes of refraction, respectively. Typically and so . In our case, , so .

The opacity per gram of gas is

(10) |

All of these quantities are interconnected, since the grain size is also a function of grain density to the power . As a result, the opacity in these upper layers is expected to depend only weakly on the assumed grain density. Thus, even if the grains are made of pure silicates with a density of g cm or of pure ice with a density 1/3 as large, the opacity should change by roughly a factor of 1.6. In practice the change should be even smaller.

## 3 Core Accretion Model

The physics of the core–nucleated accretion model are for the most part described by Pollack et al. (1996), Bodenheimer et al. (2000), Hubickyj et al. (2005), and Lissauer et al. (2009). The calculation employs the three basic elements mentioned in Section 1. The first element is essentially Eq. (1), which is used to calculate the accretion rate of solids. The value of is obtained from the numerical fits of Greenzweig and Lissauer (1992), assuming a planetesimal radius of 100 km. Other more detailed calculations of the core accretion rate (Kokubo and Ida, 2000) show in general that the rates calculated by Greenzweig and Lissauer (1992) should be reduced because they do not completely take into account the stirring and increase in relative velocity of the planetesimals as a result of gravitational interactions with the embryo. However, the accretion rate for smaller planetesimals is faster than that for 100 km. Thus our accretion rate should be a reasonable approximation to the actual rate assuming that the planetesimals determining the velocity distribution are somewhat smaller than 100 km.

As in previous calculations, we take the feeding zone from which the embryo can accrete planetesimals to extend 4 Hill sphere radii on either side of the orbit, and assume that the surface density is spatially uniform within that zone. As the planet accretes material, in the feeding zone is recalculated, again assuming uniform density, taking into account depletion of planetesimals by accretion onto the embryo and expansion of the feeding zone into undepleted regions as the embryo’s mass increases.

The second element of the code calculates the interactions between planetesimals and the gaseous envelope of the protoplanet as they fall through it (Podolak et al., 1988). The equations of motion of the planetesimals are integrated, for a variety of assumed impact parameters, taking into account the gravitational fields of the core and the envelope as well as gas drag. The critical impact parameter inside of which the planetesimal is captured is determined, and the effective capture radius (Eq. 1) is assumed to be the periapsis altitude of the critical orbit. This radius can be several times larger than the actual solid core radius, once a significant gaseous envelope has accreted, so this effect is important in speeding up the accretion rate of solids.

For those planetesimals with impact parameters less than the critical value, trajectories are calculated to determine the deposition of mass and energy into the envelope. Effects included are ablation and evaporation of planetesimal material, heat associated with phase changes, and fragmentation when the dynamical pressure of the gas exceeds the compressional strength of the planetesimal. During the earliest phases, when the envelope has low mass, planetesimals plunge all the way to the core. At later stages they can be fully ablated or fragmented in the envelope. It is assumed that the remains of the planetesimal later sink to the core, releasing gravitational energy in the process (Pollack et al., 1996). This assumption is not entirely accurate: Iaroslavitz and Podolak (2007) have made a calculation of how much of the heavy-element material should actually dissolve in the envelope. Ice should be almost entirely dissolved, but most of the organics and rock would sink to the core or to the layers just outside it. They conclude that the overall evolution would not be substantially changed if this effect were included; the accretion of the protoplanet would be slightly sped up and the core mass would be reduced to about 2/3 the calculated value, with little change in the overall metal content.

The third element of the simulation is the solution of the four differential equations of stellar structure for the gaseous envelope, with the nuclear energy term replaced by the energy delivered by accreting planetesimals and with the gravitational energy terms included. The adiabatic temperature gradient is assumed in convection zones. At the core-envelope interface the radius is set to that of the outer edge of the core, calculated from and the assumed core mean density of 3.2 g cm. Also, the luminosity is set to the energy deposition rate for the planetesimals that hit the core. Outside the core, the energy supplied by ablated and fragmented planetesimals is included as a source term in the energy equation, along with the energy derived from gas accretion and compression.

At the surface, mass is added at a sufficient rate to maintain an outer radius , and Eq. (2) is used with , as reported in Lissauer et al. (2009) from the results of three-dimensional calculations of disk-planet interaction. Otherwise the density and temperature at the surface are set to assumed nebular values , , respectively. These values are constant in time. The simulations are stopped when the gas accretion rate reaches the disk flow limit d/d, as determined from three-dimensional calculations with disk viscosity parameter (Lissauer et al., 2009). The lower curve in figure 3 of that paper gives the limiting gas accretion rate as a function of planet mass. In the various cases of the present calculation, d/d M yr when the limit is first reached.

The equation of state of the gas is taken from Saumon et al. (1995), interpolated to our assumed composition of hydrogen mass fraction X = 0.74, helium mass fraction Y = 0.243, and metal mass fraction . Although the equation of state in the outer, low-density layers is essentially that of an ideal gas, the inner regions near the core can be significantly non-ideal once the envelope has acquired sufficient mass.

The Rosseland mean opacity calculation has three components. At temperatures above 3500 K the molecular opacities of Alexander and Ferguson (1994) are used. In practice the details of the opacities in this region are unimportant because the energy transport is almost always by convection. In the temperature range 200–3500 K the molecular opacities, without grains, of Freedman et al. (2008) are used. The grain opacity calculation described in section 2 is then included in the temperature range 100–1800 K. In the higher end of this range the grain opacity can be so low at times that the molecular opacity dominates. Grains are assumed to be completely evaporated above 1800 K.

The grain calculation and the stellar envelope calculation are combined as follows. After each evolutionary time step for the envelope, grains are added at the surface in accordance with d/d, assuming a dust-to-gas ratio of 0.01. Grains are also added where the ablation/fragmentation calculation for the incoming planetesimals deposits them, typically fairly deep in the envelope. For the grains, 34 size bins are used, with the smallest at m and the largest at 2.58 mm. Incoming grains, represented by the term in Eq. (3), are always placed in the smallest bin. The envelope calculation then provides to the grain calculation a table of temperature, density, and convective velocity as a function of radius, as well as the core and envelope masses. Velocities are estimated by the stellar structure code through the mixing-length theory. An arbitrary number of distinct convection zones is allowed.

The grain evolution, involving settling and coagulation, then proceeds in the region from the surface down to 1800 K. A large number of sub-timesteps is used, adding up to a total time of , as the number of grains in any size bin is not allowed to change by more than 4% per sub-timestep. The number of sub-timesteps can be 10,000 to 100,000, so this part of the calculation is quite time-consuming. Once the grain calculation has been completed, with a new opacity provided at each depth, the opacities are fed back into the envelope calculation in the form of a table of Rosseland mean opacity as a function of mass fraction, integrated inward from the surface. Mass fraction is used rather than temperature or density because, first, the temperature is practically isothermal in the outer layers, and, second, temperature and density values can fluctuate somewhat as a result of the numerical procedure for adding mass to the envelope. The current grain size distribution is stored at each depth for use in the next time step. The envelope calculation then proceeds to the next model. Once the model has been completed, the Lagrangian zones have moved in radius with respect to the previous model. The grain number densities, tabulated as a function of radius, are interpolated onto the new grid by a procedure that conserves the total number and mass of grains. The results for the grain opacity are qualitatively quite similar to those obtained for specific envelope models by Movshovitz and Podolak (2008). In the outer zones, where small grains are continuously added from the protoplanetary disk, the opacity generally falls in the range 0.1–1 cm g, not too different from opacities derived from grains with interstellar size distribution and abundance. In contrast, at the deeper levels the typical grain size increases sufficiently so that the mean opacity can drop by up to 4 orders of magnitude relative to interstellar values!

## 4 Calculations and Results

The planet initially consists of a core of 1 M and an envelope of about M. The protoplanet is located at 5.2 AU in a solar nebula disk, with the solid surface density treated as a parameter. Based on our past calculations for the time required to reach (Hubickyj et al., 2005; Lissauer et al., 2009), the starting time is set to yr, yr, and yr for , 6, and 4 g cm, respectively. No orbital migration is included. The quantity is set to 115 K, and , where . The scale height of the gas in the disk , where is the orbital distance from the star. Once started, the evolution usually consists of three main phases. The first involves primarily accretion of solids onto the core, with a relatively low-mass envelope and a low gas accretion rate. The solids accretion rate slows down significantly near the point where the isolation mass () for the core is reached; for at 5.2 AU this mass is about 11.5 M. During the second phase, the gas accretion rate is about 3 times as high as the core accretion rate, and both are nearly constant in time. The envelope mass builds up relative to the core mass, which grows slowly. The third phase, rapid gas accretion, starts just prior to crossover mass, . The calculations reported here stop at a point well beyond crossover, where d/dt = d/dt. Beyond that point, accretion up to M is very fast, and the core mass does not change appreciably. Examples of this final phase are given in Lissauer et al. (2009); note especially Cases 2l and 2lJ, calculated with g cm.

Three main cases were calculated, with 10, 6, and 4 g cm and with all other parameters identical. These will be referred to as runs 10, 6, and 4, respectively. A fourth case (10R) was calculated with 10 g cm but without convective effects included in the grain calculation, although they were included in the actual planetary model calculation. The most closely analogous comparison case from Lissauer et al. (2009) was their Case 2lJ with 10 g cm and the same outer boundary conditions as used here; however Lissauer et al. simply set the grain opacity equal to 2% of interstellar values. That calculation reached a crossover mass of 16.16 M at a time of 2.3 Myr, and a final total mass of 1 M at 2.6 Myr. The input parameters of the four runs are shown in Table 1, which includes , the gas surface density , the surface boundary temperature , a column “conv?”, which indicates whether convection was included in the grain calculation or not, and the isolation mass .

Run | (g/cm) | (g/cm) | (K) | conv? | (M) |
---|---|---|---|---|---|

10 | 10 | 700 | 115 | Y | |

10R | 10 | 700 | 115 | N | |

6 | 6 | 420 | 115 | Y | |

4 | 4 | 280 | 115 | Y |

The results of our calculations with convection are shown in Figs. 1 and 2. For the case, the formation of the planet is significantly speeded up in comparison to case 2lJ. The early core accretion phase is almost identical to that in case 2lJ as the envelope opacity has little effect on the rate of buildup of the core. The first luminosity peak is at a level of L at 0.4 Myr; the isolation mass is approached at about 0.45 Myr. However Phase 2, during which core and envelope accretion rates are relatively low, is significantly shorter in the present calculation. The main effect of the grain settling calculation is to reduce the average opacity in the envelope, allowing it to lose heat more readily, resulting in a higher rate of envelope contraction and therefore needing a higher gas accretion rate to keep the outer boundary at . The crossover mass is almost identical to that in case 2lJ, but it is reached at a time of only 0.97 Myr. Beyond that point the gas accretion rate increases rapidly, reaching the limiting rate of M yr at yr with M and M. By way of comparison, case 2lJ reached the same d/d with M and M, but the time was yr.

Thus the main results are that, as in previous work where we examine the effect of opacity on planetary formation (Pollack et al., 1996; Hubickyj et al., 2005), the final core mass is very insensitive even to large changes in the opacity. Furthermore, compared to calculations with full interstellar grain opacity, not including grain coagulation and settling (Hubickyj et al., 2005), the formation time for Jupiter is reduced by a factor 6. The main contributor to this result is the sharp reduction in the opacity that appears in the outer radiative zone at temperatures 500-600 K. This feature is visible in the structure plots shown in Figs. 3 and 4. Interior to this minimum, molecular opacities rather than grain opacities dominate, and the opacity increases with temperature.

The properties of Runs 10, 10R, 6, and 4 at particular times in the evolution are summarized in Table 2. Here and . Case 10R is not included in the figures, because, as shown in Table 2, the results are very similar to those of 10. The small differences in the two cases are mainly a result of numerical fluctuations in the calculations. The properties are indicated at the following times: the maximum in luminosity that occurs during the main core accretion phase, the transition between Phase 1 and Phase 2 where the isolation mass of the core is approximately reached, a time in the middle of Phase 2, the time at which the crossover mass () is attained, and the time when the limiting gas accretion rate is reached. Following this last time, the planet quickly reaches its final mass through rapid gas accretion; the details of this final phase are treated in Lissauer et al. (2009).

10 | 10R | 6 | 4 | ||
---|---|---|---|---|---|

First Luminosity Peak | Time | 0.408 | 0.405 | 0.654 | 0.888 |

8.53 | 8.12 | 3.98 | 2.12 | ||

0.064 | 0.044 | 0.0061 | 0.0014 | ||

-4.85 | -4.89 | -5.84 | -6.53 | ||

28.8 | 27.7 | 17.3 | 9.86 | ||

End of Phase 1 | Time | 0.451 | 0.454 | 0.838 | 1.19 |

11.5 | 11.5 | 5.60 | 2.97 | ||

0.67 | 0.61 | 0.63 | 0.18 | ||

-5.46 | -5.52 | -6.64 | -7.49 | ||

34.7 | 35.2 | 23.1 | 15.7 | ||

Mid Phase 2 | Time | 0.775 | 0.781 | 1.295 | 2.738 |

14.0 | 14.0 | 6.50 | 3.54 | ||

7.0 | 7.0 | 3.25 | 1.77 | ||

-5.88 | -5.77 | -6.65 | -7.52 | ||

47.4 | 48.0 | 30.8 | 21.8 | ||

Crossover Point | Time | 0.968 | 0.947 | 1.625 | 3.618 |

16.09 | 16.11 | 7.50 | 4.09 | ||

-5.54 | -5.52 | -6.27 | -7.13 | ||

57.8 | 60.1 | 38.6 | 25.5 | ||

Onset of limited gas accretion | Time | 1.003 | 0.997 | 1.865 | 4.038 |

16.8 | 17.3 | 8.92 | 4.74 | ||

56.8 | 52.1 | 54.3 | 34.0 | ||

Time is in units of millions of years, Myr. Mass ( and ) is in units of Earth’s mass, M. The accretion rate ( and ) is in units of Earth masses per year, Myr. Luminosity () is in units of solar luminosity, L. Radius () is in units of Jupiter’s present equatorial radius, R.

Case 6 is calculated with , and as expected the formation time is longer, but still less than that for the calculation done by Hubickyj et al. (2005) for this level of . The final core mass is, however, significantly reduced, to 8.92 M. The first luminosity peak is reached in 0.65 Myr, somewhat later than in case 10. The core accretion rate is of course slower in case 6, but the available planetesimal supply is considerably less, a compensating effect. The maximum luminosity is typically reached when has been reduced (by accretion onto the planet) to 18% of its initial value. For case 6, the maximum is a factor 10 less than in case 10. Phase 2 is reached at a time of yr with a core mass of 5.6 M, half that of case 10 at the corresponding epoch. Mainly because of the reduced core mass (see Pollack et al., 1996), the luminosity is considerably lower in Phase 2 than it is for case 10, therefore energy is released less rapidly, contraction is slower, and the rate of gas accretion is reduced. Thus Phase 2 is lengthened, but its duration remains less than 1 Myr, giving a total formation time of 1.86 Myr.

The case 4 represents a disk with g cm at 5.2 AU. This value is only slightly greater than that of the minimum mass solar nebula, but note that our planetesimal accretion prescription neglects expulsion of solids from the planet’s accretion zone, so a somewhat higher is probably required in order to reproduce the core accretion rate calculated in this case. The first luminosity peak occurs at about the same time as in 6, but the value of the luminosity is a factor of 5 lower. Crossover is reached at a core mass of only 4 M and a luminosity of L. As shown in Figs. 5 and 6, the temperatures at a given density in the envelope structure are considerably lower than in case 10. The planet’s internal structure, both in Phase 2 and at crossover, is characterized by an inner convection zone and an outer radiative zone, with the boundary at a temperature of only 200 K. The deep minimum in opacity just outside this boundary allows for a high rate of envelope radiation and the very rapid addition of mass. The onset of limited gas accretion is reached after a time of 4.04 Myr, approximately the mean lifetime as observed for protostellar disks. The heavy element core has a mass of only 4.74 M.

## 5 Discussion and Conclusions

The formation phase of Jupiter, as determined by the concurrent accretion of solid particles and gas, has been calculated under the following major assumptions: (1) The protoplanet remains fixed in position at 5.2 AU. (2) The planet is assumed to exhibit runaway accretion of solids within a disk of 100 km planetesimals. The planetesimals are well mixed within the feeding zone of the planet; thus their surface density is constant in space but varies in time as the planet accretes. Migration of planetesimals into or out of the feeding zone is not considered. (3) The protoplanet is assumed to be the sole dominant mass in the region; other embryos which might compete for planetesimals are not considered. (4) The solids accretion rate is determined from Safronov’s formula (eq. 1) with gravitational enhancement factor obtained from the 3-body numerical fits of Greenzweig and Lissauer (1992). (5) All of the planetesimal material accreted by the protoplanet is assumed to fall into the core; thus what is referred to as ‘core mass’ in this paper is really the total heavy element excess over solar abundances. (6) The gas accretion rate is determined by the requirement that the outer boundary of the protoplanet remain at the modified accretion radius (eq. 2). Three-dimensional simulations of a planet embedded in a disk (Lissauer et al., 2009) determine in part the region of a planet’s gravitational influence. (7) The opacity caused by grains in the envelope of the protoplanet is determined by a detailed calculation of grain settling and coagulation, with the grain size distribution and the resulting grain opacity being recalculated at every level in the envelope at every time step. It is this last assumption that represents a major improvement over past calculations of Jupiter’s formation, and it results in a significant reduction in the formation time scale as compared to past calculations. There are approximations involved in the treatment of convection zones in the grain calculations, but the comparison of cases 10 and 10R shows that these have a negligible effect on the results. The grains in the incoming gas are composed of silicates only, with an abundance by mass of 1% that of the gas. Molecular opacities are included (Freedman et al., 2008).

The main parameter varied in this set of calculations is the solid surface density at 5.2 AU in the primordial disk; values considered are 10, 6, and 4 g cm. As in previous simulations of Jupiter’s formation using the prescription of Lissauer (1987) for the runaway accretion of solids, the results show that there are three main phases. Phase 1 involves primarily the accretion of solid particles to build up the core; the mass of the gaseous envelope is fairly negligible. Phase 2 begins when the isolation mass is reached; accretion of gas takes place at about three times the rate of accretion of solid particles, and the envelope builds up in mass relative to the core. Phase 3, rapid gas accretion, begins close to the point where core and envelope masses are equal. The calculations in this paper are taken beyond the crossover mass up to the point where the gas accretion rate becomes limited by the rate at which the disk can supply the mass, as determined by three-dimensional simulations ( M yr). The accretion up to the final mass of Jupiter takes place quickly after that point (Lissauer et al., 2009).

The formation time to crossover is dominated by Phase 2 in all simulations. Case 10 may be compared with a very similar calculation (Lissauer et al., 2009) in which the grain opacity was simply assumed to be 2% that of interstellar grains. Phase 1 in the present case was completed in about the same time as in the comparison case, but Phase 2 (0.52 Myr) was a factor of 3.6 shorter, giving a total formation time of about 1 Myr, compared with 2.6 Myr in the comparison case. Correspondingly, the mean luminosity in phase 2 was a factor of 3.6 higher in the present run than in the comparison case. The main cause of the increase in luminosity is the very low opacity in the region around 600 K in the present model (Fig. 3), with a minimum of opacity of about cm g as compared with cm g at the same temperature in the comparison case. Some layers in the present case have opacities higher than 2% interstellar, but they are mainly in the outermost layers that are either convective or rather optically thin, so they do little to influence the flow of radiation. In spite of the difference in formation time for the two cases, the core mass at onset of limiting gas accretion was practically the same in both cases, about 17 M. This mass is higher than that usually derived for the present Jupiter (Saumon and Guillot, 2004) but it is close to that predicted from an independent calculation of the equation of state by Militzer et al. (2008).

The question then arises whether a reduction in , which is known to result in a smaller core mass, still leads to a reasonable formation time. Cases 6 and 4 show that this is indeed the case. The formation times are 1.9 and 4.0 Myr, respectively. The ratios of Phase 2 times for the three cases scale as /, as derived by Pollack et al. (1996) under the assumption that planetesimal accretion provides most of the luminosity during Phase 2. The value 4 Myr is within the range of expected lifetimes of observed circumstellar disks around young stars (Haisch et al., 2001). Note that the case has a formation time only 14% as long (and Phase 2 time only 8% as long) as the corresponding case in Hubickyj et al. (2005) with grain opacities set at 2% of interstellar values, while the case took 40% as long as the corresponding case () in that earlier work.

The main conclusions of this work are: (1) Within the basic framework of a core-nucleated accretion model for Jupiter’s formation, the effects of grain settling and coagulation in the gaseous envelope of a protoplanet are significant in determining the gas accretion rate and formation time. (2) The formation time for a Jupiter at 5.2 AU with a solid surface density of 10 g cm (about 3 times that in the minimum mass solar nebula) is about 1 Myr, significantly shorter than that in comparable comparison cases. If is taken to be 4 g cm, the formation time is still reasonably short, about 4 Myr. (3) In the case with g cm, Jupiter’s core mass is calculated to be 4.74 M, consistent with the Jupiter structure models computed by Saumon and Guillot (2004). No arbitrary cutoff in the core accretion rate was needed to obtain this result.

Further work is needed, however, to confirm this result. The main points which need to be improved include, first, the planetesimal accretion model, which should be replaced by a statistical calculation taking into account the evolution of the planetesimal size distribution and the effects of the protoplanet itself on the planetesimal velocity distribution (Weidenschilling et al., 1997). Numerous estimates of the Phase 1 accretion time have been published. One example is the (statistical + N-body) calculation of Inaba et al. (2003) which takes into account planetesimal fragmentation as well as the effect of the gaseous envelope in increasing the effective cross-section. A planetesimal size of 10 km is assumed. If the opacity in the planetary envelope is reduced by a factor 100 compared with interstellar values, a core approaching 10 M is reached in 5 Myr at 5.2 AU; however is required to be about 5 times that in the minimum mass solar nebula. This and other simulations, along with the present results, suggest that the formation time for Jupiter may be controlled by Phase 1, not Phase 2. Second, the effects of the ices should be included in the calculations of grain settling and coagulation. Third, the effects of the dissolved accreted solid material, primarily ices, should be included in the equation of state and opacity of the envelope and in the calculation of the solid core mass.

## Acknowledgments

This work was supported in part by NASA Origins grant NNX08AH82G. M. P. acknowledges the support of the Israel Academy of Sciences. N. M. acknowledges assistance from the Sackler Institute of Astronomy.

## References

- Alexander and Ferguson (1994) Alexander, D. R., Ferguson, J. W., 1994. Low-temperature rosseland opacities. Astrophys. J. 437, 879–891.
- Alibert et al. (2005) Alibert, Y., Mordasini, C., Benz, W., Winisdoerffer, C., 2005. Models of giant planet formation with migration and disk evolution. Astron. Astrophys. 434, 343–353.
- Bodenheimer et al. (2000) Bodenheimer, P., Hubickyj, O., Lissauer, J. J., 2000. Models of the in situ formation of detected extrasolar giant planets. Icarus 143, 2–14.
- Bodenheimer and Pollack (1986) Bodenheimer, P., Pollack, J. B., 1986. Calculations of the accretion and evolution of giant planets: the effects of solid cores. Icarus 67, 391–408.
- Cameron (1973) Cameron, A. G. W., 1973. Accumulation processes in the primitive solar nebula. Icarus 18, 407–450.
- D’Angelo et al. (2003) D’Angelo, G., Kley, W., Henning, T., 2003. Orbital migration and mass accretion of protoplanets in three-dimensional global computations with nested grids. Astrophys. J. 586, 540–561.
- Fortier et al. (2007) Fortier, A., Bienvenuto, O. G., Brunini, A., 2007. Oligarchic planetesimal accretion and giant planet formation. Astron. Astrophys. 473, 311–322.
- Freedman et al. (2008) Freedman, R. S., Marley, M. S., Lodders, K., 2008. Line and mean opacities for ultracool dwarfs and extrasolar planets. Astrophys. J. Suppl. 174, 504–513.
- Greenzweig and Lissauer (1992) Greenzweig, Y., Lissauer, J. J., 1992. Accretion rates of protoplanets. ii. gaussian distribution of planetesimal velocities. Icarus 100, 440–463.
- Haisch et al. (2001) Haisch, Jr., K. E., Lada, E. A., Lada, C. J., 2001. Disk frequencies and lifetimes in young clusters. Astrophys. J. 553, L153–L156.
- Henyey et al. (1964) Henyey, L., Forbes, J., Gould, N., 1964. A new method of automatic computation of stellar evolution. Astrophys. J. 139, 306–317.
- Hubickyj et al. (2005) Hubickyj, O., Bodenheimer, P., Lissauer, J. J., 2005. Accretion of the gaseous envelope of jupiter around a 5–10 earth-mass core. Icarus 179, 415–431.
- Iaroslavitz and Podolak (2007) Iaroslavitz, E., Podolak, M., 2007. Atmospheric mass deposition by captured planetesimals. Icarus 187, 600–610.
- Ikoma et al. (2000) Ikoma, M., Nakazawa, K., Emori, H., 2000. Formation of giant planets: dependences on core accretion rate and grain opacity. Astrophys. J. 537, 1013–1025.
- Inaba et al. (2003) Inaba, S., Wetherill, G. W., Ikoma, M., 2003. Formation of gas giant planets: core accretion models with fragmentation and planetary envelope. ICarus 166, 46–62.
- Kokubo and Ida (2000) Kokubo, E., Ida, S., 2000. Formation of protoplanets from planetesimals in the solar nebula. Icarus 143, 15–27.
- Lissauer (1987) Lissauer, J. J., 1987. Timescales for planetary accretion and the structure of the protoplanetary disk. ICarus 69, 249–265.
- Lissauer et al. (2009) Lissauer, J. J., Hubickyj, O., D’Angelo, G., Bodenheimer, P., 2009. Models of jupiter’s growth incorporating thermal and hydrodynamic constraints. Icarus 199, 338–350.
- Militzer et al. (2008) Militzer, B., Hubbard, W. B., Vorberger, J., Tamblyn, L., Bonev, S. A., 2008. A massive core in jupiter predicted from first-principles simulations. Astrophys. J. 688, L45–L48.
- Mizuno (1980) Mizuno, H., 1980. Formation of the giant planets. Prog. Theor. Phys. 64, 544–557.
- Movshovitz and Podolak (2008) Movshovitz, N., Podolak, M., 2008. The opacity of grains in protoplanetary atmospheres. Icarus 194, 368–378.
- Ormel and Cuzzi (2007) Ormel, C. W., Cuzzi, J. N., 2007. Closed-form expressions for particle relative velocities induced by turbulence. Astron. and Astrophys. 466, 413–420.
- Perri and Cameron (1974) Perri, F., Cameron, A. G. W., 1974. Hydrodynamic instability of the solar nebula in the presence of a planetary core. Icarus 22, 416–425.
- Podolak (2003) Podolak, M., 2003. The contribution of small grains to the opacity of protoplanetary atmospheres. Icarus 165, 428–437.
- Podolak et al. (1988) Podolak, M., Pollack, J. B., Reynolds, R. T., 1988. Interactions of planetesimals with protoplanetary atmospheres. Icarus 73, 163–179.
- Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., Lissauer, J. J., Podolak, M., Greenzweig, Y., 1996. Formation of the giant planets by concurrent accretion of solids and gas. Icarus 124, 62–85.
- Pollack et al. (1985) Pollack, J. B., McKay, C. P., Christofferson, B., 1985. A calculation of the rosseland mean opacity of dust grains in primordial solar system nebulae. Icarus 64, 471–492.
- Safronov (1972) Safronov, V. S., 1972. Evolution of the protoplanetary cloud and formation of the earth and planets. In: Safronov, V. S. (Ed.), Evolution of the protoplanetary cloud and formation of the earth and planets., by Safronov, V. S.. Translated from Russian. Jerusalem (Israel): Israel Program for Scientific Translations. Keter Publishing House.
- Saumon et al. (1995) Saumon, D., Chabrier, G., van Horn, H., 1995. An equation of state for low-mass stars and giant planets. Astrophys. J. Suppl. 99, 713–741.
- Saumon and Guillot (2004) Saumon, D., Guillot, T., 2004. Shock compression of deuterium and the interiors of jupiter and saturn. Astrophys. J. 609, 1170–1180.
- Stevenson (1982) Stevenson, D. J., 1982. Formation of the giant planets. Planet. Space Sci. 30, 755–764.
- Weidenschilling (1984) Weidenschilling, S. J., 1984. Evolution of grains in a turbulent solar nebula. Icarus 60, 553–567.
- Weidenschilling et al. (1997) Weidenschilling, S. J., Spaute, D., Davis, D. R., Marzari, F., Ohtsuki, K., 1997. Accretional evolution of a planetesimal swarm. 2. the terrestrial zone. Icarus 128, 429–455.