Dust coagulation in the vicinity of a gap-opening Jupiter-mass planet
We analyze the coagulation of dust in and around a gap opened by a Jupiter-mass planet. To this end, we carry out a high-resolution magnetohydrodynamic (MHD) simulation of the gap environment, which is turbulent due to the magnetorotational instability. From the MHD simulation, we obtain values of the gas velocities, densities and turbulent stresses a) close to the gap edge, b) in one of the two gas streams that accrete onto the planet, c) inside the low-density gap, and d) outside the gap. The MHD values are then supplied to a Monte Carlo dust coagulation algorithm, which models grain sticking and compaction. We consider two dust populations for each region: one whose initial size distribution is monodisperse, with monomer radius equal to 1 m, and another one whose initial size distribution follows the Mathis-Rumpl-Nordsieck distribution for interstellar dust grains, with an initial range of monomer radii between 0.5 and 10 m. Our Monte Carlo calculations show initial growth of dust aggregates followed by compaction in all cases but one, that of aggregates belonging to the initially monodisperse population subject to gas conditions outside the gap. In this latter case, the mass-weighted (MW) average porosity of the population reaches extremely high final values of 98%. The final MW porosities in all other cases range between 30% and 82%. The efficiency of compaction is due to high turbulent relative speeds between dust particles. Future studies will need to explore the effect of different planet masses and electric charge on grains.
Subject headings:magnetohydrodynamics – planets and satellites: formation – protoplanetary disks – planet-disk interactions – turbulence
Numerical investigations of the dynamics of dust particles in the neighborhood of a gap-opening planet reveal that they tend to concentrate at gap edges, for planet masses between and , with the mass of Jupiter (Fouchet et al., 2007; Zhu et al., 2014; Owen, 2014; Gonzalez et al., 2015; Picogna & Kley, 2015). These accumulations are of great interest not only because they could facilitate planetesimal formation, but also because they could supply solid material to the planetary atmosphere through deposition of small grains, and hence contribute to the atmospheric opacity, although perhaps modestly (Ormel, 2014). In particular, a Jupiter-mass planet accretes predominantly particles smaller than 10 m (Paardekooper, 2007), which places a tight constraint on its solid enrichment.
One aspect of the evolution of dust in the vicinity of a gap-opening planet that has not been considered is the internal structure of solid aggregates. The porosity of dust particles affects their aerodynamic coupling to the disk gas. Their structure can also affect the ionization state of the protoplanetary disk: the increased surface area associated with internal voids within the dust aggregates allows for fast electron recombination, which can lower the ionization levels of the gas phase, and thus render the disk gas stable to the magnetorotational instability, or MRI (Balbus & Hawley, 1991), a likely candidate to produce turbulent viscosity in disks [although other mechanisms have been proposed to transfer a protoplanetary disk’s angular momentum, such as centrifugally-driven winds and jets (Blandford & Payne, 1982) and magnetic breaking (Matsumoto & Tomisaka, 2004)].
Dust in the vicinity of a gap-opening planet is subject to a dynamic and energetic environment. For example, inside the gap, shocks are generated when streams of gas flow through the L1 and L2 points towards the interior of the planet’s Hill sphere (Lubow et al., 1999). These shocks can modify the crystalline structure of silicates through thermal annealing. Knowledge of the porous structure of dust aggregates is essential to determine their heat conductivity during these events.
Previous studies of dust coagulation and dust charging in protoplanetary disks show that electric charges lead to larger, more massive and more porous aggregates than in the case of neutral coagulation (Matthews et al., 2012). In turn, the evolution of the aggregate porosity may define a region in weakly ionized disks where the growth of sub-micron-sized particles becomes stalled due to electrostatic repulsion between negatively-charged aggregates (Okuzumi et al., 2011).
The dusty ingredient in the MRI recipe could also be crucial for the dynamics of giant-planet circumplanetary disks. Turner et al. (2014) showed that, under some circumstances, the circumjovian disk has sufficient electrical conductivity for magnetic forces to drive accretion stresses. But if the disk contains enough dust, MRI turbulence is unlikely to occur, even in the presence of ionizing X-ray radiation. The circumjovian disk can develop a substantial magnetically inactive “dead zone”, where regular satellites could form.
It is thus evident that the size and internal structure of dust grains are key to understanding the birth environment of planets that are capable of opening gaps. To better characterize the spatial distribution of these dust properties around one such planet, in this study we combine a numerical model of MRI turbulence with a dust-coagulation algorithm. Our results will provide a clearer picture of the dust aggregate structure around a young, Jupiter-mass planet. In the interest of specificity, we use our models to invoke a scenario in which Jupiter underwent migration within the primitive solar nebula, the so-called Grand Tack hypothesis (Walsh et al., 2011). According to this hypothesis, Jupiter migrated from inside its current position (but beyond the Main Asteroid Belt) to a distance of about 1.5 AU from the Sun, where it encountered an orbital resonance with the trailing Saturn. Both planets then reversed their orbital motion outwards. The Grand Tack helps reproduce the formation of Mars analogs with the correct mass from a disk of planetary embryos.
2.1. Disk model
We use a magnetohydrodynamic (MHD) protoplanetary disk model in the local shearing box approximation (Hawley et al., 1995), in which the MHD equations are solved in a rectangular coordinate system that corotates with the disk at a fiducial orbital radius , with angular frequency . In this system, the axis is oriented along the radial direction, the axis along the azimuthal direction, and the axis is parallel to the disk’s angular momentum vector. Our solver is the Athena code (Stone et al., 2008), a grid-based algorithm that has been extensively tested and employed for various studies of protoplanetary disks.
Our numerical setup is similar to that of Zhu et al. (2013): the local, 3D disk model is isothermal, and does not include vertical stratification of the gas density. The box dimensions are , where is the disk scale height. We use a numerical resolution twice as large as that in Zhu et al. (2013), namely 64 grid cells per . In the code’s system of units, the gas sound speed , the initial gas density , and the angular frequency are all equal to 1. The initial magnetic field strength is given by the plasma beta, , and the initial field configuration corresponds to a non-zero net vertical magnetic flux.
To model the gravitational effect of a planet on the surrounding disk gas, we place a cylindrical potential at the center of the box, with axis coincident with the box’s vertical axis. As in Zhu et al. (2013), the planet potential is given by
where is the gravitational constant, is the planet mass, is the distance to the axis, and is a smoothing length to avoid small time steps close to the source. The potential is also smoothed out to take into account discontinuities at the box boundaries:
where is the cutoff distance to the axis beyond which the potential flattens out, and is a smoothing length at the cutoff radius. Following Zhu et al. (2013), we set and .
In Athena, we express the planet mass in terms of the so-called thermal mass , the mass at which the Hill radius and the Bondi radius of the planet are comparable to the disk scale height (Rafikov, 2006). For the minimum-mass solar nebula model [MMSN; Hayashi (1981)], the thermal mass is (Dong et al., 2011)
where is the planet’s semi-major axis, and is the mass of Earth. Within the Grand Tack scenario, Jupiter begins its inward migration at AU, reaches 1.5 AU, and reverses its motion towards its current orbital position at 5.2 AU (Walsh et al., 2011). Since we are not modeling the actual migration, but effectively only a “snapshot” of that process, we choose AU, because it is an interesting region of the solar nebula in terms of the dynamical evolution of chondritic parent bodies (Walsh et al., 2011). In that case, from relation (3) a value of gives our planet a mass .
We run this high-resolution setup for 51 orbits (265 yr at 3 AU in the MMSN), and use the resulting data from the last 12 orbits, after the gas flow reaches an approximate steady state, as input to the coagulation code described next.
2.2. Dust coagulation model
Our treatment of dust coagulation follows the implementation of Ormel et al. (2007), which in turn is based on the Monte Carlo model by Gillespie (1975). A population of particles (we refer to monomers and aggregates jointly as particles) is assumed to be uniformly distributed inside an abstract volume . The collision rate between particles and is calculated as
where is the collision cross-section, and are the respective particle radii, and is the relative speed between the two particles (described below). In this coagulation algorithm, we track the radius of any particle , , through the enlargement factor , where is the extended volume of the particle (i.e., the volume corresponding to the geometric cross-section of the aggregate), and is the volume of a monomer of radius . We then have . We also track the evolution of the particle mass .
The particle pair that will be involved in each collision is chosen using random numbers and partial sums of the collision rate . We treat two possible outcomes of a binary collision: particles either merely stick without further restructuring, or are compacted at the expense of internal voids. These possibilities occur when the kinetic energy of collision () is less than the critical energy to initiate compaction (), or exceeds it, respectively (Ormel et al., 2007). To conserve the total number of particles in the volume after a collision, one of the remaining particles is randomly chosen and duplicated. The duplication procedure requires that the volume be increased to preserve the spatial density of dust, .
In order to include the effect of the evolving porosity of aggregates formed in collisions, the enlargement factor of an aggregate is calculated according to Eqs. (15) or (17) of Ormel et al. (2007), depending on whether or , respectively.
The relative speed in Eq. (4) has contributions from Brownian motion and from turbulence. We use an expression for the turbulent relative speed based on Eq. (10) of Ormel et al. (2008). We insert into that expression values of the gas turbulent velocity, turbulent viscosity (given by magnetic and hydrodynamic stresses) and gas density, all obtained from the MHD simulation described above. These values are taken from four different regions in the vicinity of the planet, as shown in Fig. 1, which portrays the system at the end of our MHD simulation. The contours represent the gas surface density, obtained by averaging the gas density over the vertical extent of the box. The gap opened by the planet is clearly defined. The regions R1 through R4 are each divided into 128 subregions of grid cells each. We call these subregions macrocells. The gas variables are averaged inside each macrocell, and are then fed to the Monte Carlo coagulation code. Note that although the volume where each particle population resides is effectively associated to each macrocell, it is not the volume of the macrocell. The particles do not have position coordinates associated with them.
We set the number of particles in the volume to , and we use two different initial conditions for the particle sizes: =1 m, and distributed according to the Mathis-Rumpl-Nordsieck (MRN) distribution of dust grains in the interstellar medium, (Mathis et al., 1977). In the latter case, the range of monomer radii is 0.5 m. In both cases, the monomer bulk density is =3 g cm.
Figure 2 shows radial profiles of various quantities associated with the MHD flow. These quantities have been averaged over time (the last 12 orbits of the simulation) and over the and directions (regions for which are not included). The vertical gray bars mark the radial position of the reference regions. In Fig. 2a, the gas density, normalized by the initial density, exhibits a drop of a factor of inside the gap opened by the planet with respect to the surrounding gas. Panel c shows a corresponding factor-of-8 drop in the Maxwell stress, which is the main contributor to the turbulent viscosity in accretion disks (Hawley et al., 1995). Figure 2b shows the turbulent gas velocity in units of the gas sound speed. Overall, turbulent gas speeds are relatively high, reflecting the turbulent strength for our case of a non-vanishing initial magnetic flux. Finally, Fig. 2d reveals that the mean vertical magnetic field has two maxima at either side of the planet’s radial position, close to the gap edges. This is in contrast to the centrally-peaked profile of measured by Zhu et al. (2013), for a lower planet mass of ( at 3 AU in the MMSN), modeled with half the numerical resolution.
The coagulation evolution of the particle populations with an initial monomer radius of 1 m is shown in Fig. 3. Each curve color refers to one of the four regions of Fig. 1. Figure 3a shows the mean radius of the aggregates, normalized by the initial monomer radius m. Sticking without restructuring occurs more effectively than compaction for grains subject to the flow conditions in region R4 (green curve), where gas velocities are relatively low, as can be seen in Fig. 2b. On the other hand, growth is stalled for particles in regions R1, R2 and R3 (black, red and blue curves, respectively), which have higher turbulent gas velocities.
The evolution of porosity, defined by the mass-weighted average enlargement factor as , is shown in Fig. 3b. The compaction of aggregates in regions R2 and R3 is evident, and best understood by looking at the relative impact speeds of the aggregate pairs involved in each collision, as shown in Fig. 3c. The highest relative speeds are generated in region R3, as expected from the magnitude of the gas velocity there (Fig. 2b).
The case of an initial MRN-type size distribution is qualitatively similar to the monodisperse case, but with less overall growth. Figure 4a, which also displays the mean particle radius, shows that aggregates in all four regions have their growth reversed during the span of the calculation, with the largest growth achieved only by particles in region R4; they grow by % from the initial mean radius of 0.83 m. The resulting dusty structures are also more compact than in the monodisperse case, with final porosities in region R3 of 30%, compared to 48% in the same region when m.
Our results indicate that growth by sticking occurs most effectively outside the gap (region R4), where the gas density is higher and the gas turbulent velocity is lower (Figs. 2a,b). In this region, we also found that the mean vertical component of the magnetic field, , is relatively low (Fig. 2d).
Region R2, located in one of the two planetary wakes, is characterized by lower densities than R4, but higher gas velocities. It is interesting that R2 and R3 (the latter one located in the middle of the low-density gap) are two regions where has a local maximum and a local minimum, respectively, yet the level of aggregate compaction is similar.
The evolution of the aggregate radius exhibits different trends in the high-gas-density region R4, for the two dust populations with different initial size distributions. In the case of the initially monodisperse size distribution, the mean aggregate radius in R4 increases with time (although somewhat more slowly in the last third of the coagulation calculation), whereas in the case of the MRN-type initial distribution the mean radius decreases in the last four orbits. It is possible, then, that at least in some regions in the vicinity of the planet, different initial dust size distributions could lead to significant differences in aggregate growth rates. This point needs to be addressed in future investigations, using, for example, semi-analytical distributions obtained for coagulation/fragmentation equilibrium (Birnsteil et al., 2011).
Perhaps the most salient omissions in our MHD calculation are non-ideal effects caused by the varying ionization fraction of the disk. Magnetic field lines can be easily drawn into the planet-induced gap, and the Hall effect dominates the onset of the MRI depending on the relative orientation between the vertical component of the field and the rotation axis of the disk (Keith & Wardle, 2015). Nevertheless, certain combinations of column density, magnetic field strength, and dust content can render the gap susceptible to the MRI in the ideal-MHD regime, such as we assume here.
Our use of an orbital radius of 3 AU for a Jupiter-mass planet is consistent with the Grand Tack hypothesis. Our MHD gap model can be assumed to represent a snapshot of the Grand Tack scenario. In this setting, it is instructive to look at our calculations of aggregate porosities within the context of the formation of meteorite components. In the CV chondrite Allende, chondrule rims, composed of sub-micron grains, may have accreted with high porosities of 70–80% (Bland et al., 2011). If Jupiter’s migration did indeed occur and was contemporaneous with the formation of meteorite parent bodies, such high porosities may have occurred outside the gap, where turbulent gas velocities were lower. In any case, further exploration of our simulation parameters is needed to determine the effect of varying magnetic field strengths and geometries, a different equation of state, and electrical charging of dust grains. The latter effect could delay the growth of rim-forming grains due to electrostatic repulsion (Okuzumi et al., 2011).
Our Monte Carlo calculations of dust growth in and around a gap opened by a Jupiter-mass planet, in the presence of turbulence generated by the magnetorotational instability, indicate that aggregate compaction is effective inside the gap, near the gap edge, and in the planetary wake. This is due to the high turbulent relative speeds between aggregates (in the range 20 – 600 cm/s). In these regions, the relative kinetic energies frequently exceed the minimum energy required for restructuring. The lowest porosities occur inside the gap, with a value of if the initial size distribution of the population is monodisperse (with the radius of all initial monomers equal to 1 m), and if the initial size distribution follows a power law with exponent -3.5, as in the MRN distribution for interstellar dust.
The most porous aggregates occur in a high-gas-density region outside the gap, where the turbulent relative speeds are lower. In that region, porosities reach extremely high values of for the initially monodisperse population. In the case of the MRN-type population, the porosity reaches .
The outcome of our MHD model of a planet-induced gap, with the planet’s mass equal to the mass of Jupiter, differs noticeably from some of the results obtained by Zhu et al. (2013). In particular, the radial profile of the vertical component of the magnetic field is doubly peaked in our simulation, with one peak at either side of the planet’s radial position , whereas in the Zhu et al. (2013) study the profile exhibits only one peak at , for a planet 6.7 times less massive. Future investigations need to determine the effect of the planetary mass on this radial profile, as the vertical component of the magnetic field plays a key role in the susceptibility of a gap to the MRI (Keith & Wardle, 2015).
We will address the effect of grain charging in subsequent work, using a more sophisticated numerical scheme to treat the growth of dust aggregates.
- Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214
- Birnsteil, T., Ormel, C. W., & Dullemond, C. P. 2011, A & A, 525, A11
- Bland, P. A., Howard, L. E., Prior, D. J., Wheeler, J., Hough, R. M., & Dyl, K. A. 2011, Nat. Geosci., 4, 244
- Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
- Cuzzi, J. N., Dobrovolskis, A. R., & Champney, J. M. 1993, Icarus, 106, 102
- Dong, R., Rafikov, R. R., & Stone, J. M. 2011, ApJ, 741, 57.
- Fouchet, L., Maddison, S. T., Gonzalez, J.-F., & Murray, J. R. 2007, A&A, 474, 1037.
- Gillespie, D. T. 1975, J. Atmos. Sci., 32, 1977.
- Gonzalez, J.-F, Laibe, G., Maddison, S. T., Pinte, C., & Ménard, F. 2015b, Planetary & Space Science, preprint (arXiv:1506.00430).
- Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
- Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35.
- Keith, S. L., & Wardle, M. 2015, MNRAS, 451, 1104
- Lubow, S. H., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001
- Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
- Matsumoto, T., & Tomisaka, K. 2004, ApJ, 616, 266
- Matthews, L. S., Land, V., & Hyde, T. W. 2012, ApJ, 744, 8
- Okuzumi, S., Tanaka, H., Takeuchi, T., & Sakagami, M.-A. 2011, ApJ, 731, 96.
- Ormel, C. W., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 461, 215.
- Ormel, C. W., Cuzzi, J. N., & Tielens, A. G. G. M. 2008, ApJ, 679, 1588.
- Ormel, C. W. 2014, ApJ, 789, L18.
- Owen, J. 2014, ApJ, 789, 59.
- Paardekooper, S.-J. 2007, A&A, 462, 355.
- Picogna, G., & Kley, W. 2015, A&A, preprint (arXiv:1510.01498).
- Rafikov, R. R. 2006, ApJ, 648, 666.
- Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137.
- Turner, N. J., Lee, M. H., & Sano, T. 2014, ApJ, 783, 14.
- Walsh, K. J., Morbidelli, A., Raymond, S. N., O’Brien, D. P., & Mandell, A. 2011, Nature, 475, 206.
- Zhu, Z., Stone, J. M., & Rafikov, R. R. 2013, ApJ, 768, 143.
- Zhu, Z., Stone, J. M., Rafikov, R. R., & Bai, X.-N. 2014, ApJ, 785, 122.