EDP Sciences
Free Access
Issue
A&A
Volume 591, July 2016
Article Number A86
Number of page(s) 8
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201527732
Published online 20 June 2016

© ESO, 2016

1. Introduction

The conventional belief is that planetesimals are the building blocks of planets. However, the process of dust grains that grow into planetesimals in a protoplanetary disk is poorly understood. There are many theoretical difficulties in understanding the formation process of planetesimals.

The radial drift barrier is one of the most serious problems to understand in planetesimal formation. The disk is at hydrostatic equilibrium, in which the pressure and the centrifugal force balances the gravity of the central star. The pressure gradient causes the disk to rotate at a sub-Keplerian velocity. Since the dust particles tend to rotate at a Keplerian velocity, they feel a headwind from the gas flow. As a result, they lose their angular momentum through gas drag and drift inward (Whipple 1972; Adachi et al. 1976; Weidenschilling 1977). Because the dust growth timescale is generally significantly longer than the radial migration timescale, it is likely that the dust particles are lost from the disk.

Many ideas have been proposed to bypass the radial drift barrier. A straightforward solution is that planetesimals form rapidly before dust grains drift significantly. It has recently been proposed that the fluffy aggregates with a low bulk density grow more quickly and form icy planetesimals when they enter the Stokes regime (Okuzumi et al. 2012; Kataoka et al. 2013a,b). Another mechanism for rapid planetesimal formation is self-gravitational collapse to form massive clumps in dust-rich regions (the so-called streaming instability (SI); Youdin & Goodman 2005; Youdin & Johansen 2007; Johansen & Youdin 2007) and the following break-up of the massive clumps into smaller pieces. The SI is induced by inward radial drift of solid material that is due to gas pressure support. The instability modulates the spatial density pattern of solid material and forms dense particle clumps. When concentrated enough, these clumps may contract and become planetesimals as a result of self-gravity. We note that the SI does not work for small dust grains. The growth rate of SI is at its greatest when the particles are frictionally coupled with the gas on an orbital timescale. For instance, particles with a radius of 30 cm are most favorable for SI growth at the orbital distance r = 5 AU in the minimum mass solar nebula (MMSN; Hayashi et al. 1985) model (Johansen & Klahr 2011). For SI to be triggered, the particles must be decimeter-sized bodies in the inner part of protoplanetary disks.

The migration trap is another idea put forward to bypass the radial drift barrier. When a protoplanetary disk has a local pressure maximum, called a radial pressure bump, solid materials are trapped at the point where the pressure is radially maximized (Whipple 1972; Haghighipour & Boss 2003a,b). Radial pressure bumps are believed to form in various parts of protoplanetary disks, such as at the disk inner edge, at the dead zone inner edge, and at the snow line (Johansen et al. 2014). From recent submillimeter observations, it has been discovered that some protoplanetary disks have an asymmetrical dust profile (Casassus et al. 2013; Pérez et al. 2014), and some of these might possess a radial pressure bump (van der Marel et al. 2013).

A typical dust-to-gas mass ratio in a protoplanetary disk is ~ 0.01. Since the back-reaction is proportional to the dust-to-gas mass ratio, neglecting the back-reaction is a good approximation for dust and gas motion in a general situation. Haghighipour & Boss (2003a,b) considered the disk gas, and thus a pressure bump, to be a steady structure. When a pressure bump accumulates dust particles, however, the back-reaction from particles to gas becomes stronger (called the dust-dominated phase; Nakagawa et al. 1981, 1986), and eventually, the back-reaction alters the gas velocity and density distribution of the pressure bump (Kato et al. 2012), which affects its capability of trapping dust particles. It is important to take the back-reaction in the simulations of dust accumulation at the pressure bump into account.

Kato et al. (2009, 2010, 2012) partially confirmed the effects of the back-reaction. They investigated the inhomogeneous growth of a magneto-rotational instability (MRI) in a local area embedded in a protoplanetary disk. They found that the inhomogeneous growth of an MRI formed a quasi-steady radial pressure bump and meter-sized boulders were trapped in it. They found that the highest dust density was significantly lower with the back-reaction than without it. The back-reaction of dust onto gas, when included, brought the gas azimuthal velocity closer to Keplerian rotation in the dust-dense region, which was not favorable for a subsequent concentration of dust particles. Their simulations, however, were computationally expensive and included many complex settings. They did not follow the long-term evolution of the dust and gas distributions. In particular, because they also included MRI, the amount of gas surface density evolution that was contributed to by the back-reaction was not clear.

In this paper, we investigate the simultaneous evolution of the dust and gas density profile at a simplified radial pressure bump. For simplicity, we treat the disk gas as if it were in the (non-magnetized) hydrodynamic fluid. In addition, we set the pressure bump as an initial condition. The initial pressure balance is maintained by the Coriolis force from azimuthal velocity. Then the pressure bump is steady in the absence of a back-reaction from dust onto gas, even without the source that produced the bump (e.g., MRI zonal flows, edges of the MRI dead zone, a snow line, etc.). Kato et al. (2009, 2010, 2012) assumed that the MRI turbulence was due to a radially inhomogeneous external magnetic field and that the MRI does not occur anymore after a bump structure with local rigid rotation is established. The dust and gas motions are mutually coupled by the drag forces. The dust particles tend to accumulate at the point of highest pressure, while the dust accumulation alters the gas pressure profile. Since the altered gas profile also affects the dust motion, the simultaneous evolution of dust and gas dynamics must be investigated. In addition, the self-consistent approach enables us to consider the effects of SI. Since the SI becomes more effective at the dust-dense region (e.g., Youdin & Johansen 2007; Johansen & Youdin 2007), a pressure bump is expected to be a favorable location for the growth of SI.

This paper is organized as follows. In Sect. 2 we describe the equations and initial setup of our simulations. In Sect. 3 we show the simulation results. In Sect. 4 we discuss details of the mechanisms by which the pressure bump evolves. We also discuss the relevance of our results to the process of planetesimal formation. Our conclusions are presented in Sect. 5.

2. Equations and model

2.1. Equations

We considered the local Cartesian frame to rotate at Keplerian frequency Ω at a reference distance r0 from a central star to study dust and gas motions. The local Cartesian coordinates are (x,y,z), where x = rr0, y is the tangential distance from the origin, and z is the vertical distance from the disk midplane. We assumed an axisymmetric disk and performed the 1D (radial) and 2D (radial-vertical) simulations. Differences between the solid and the gas azimuthal velocities cause them to migrate radially as a result of angular momentum exchange through gas drag. We calculated the azimuthal velocities at each grid cell in 1D or 2D simulations. For the disk gas, we used isothermal hydrodynamic equations, where ρg is the density of the disk gas, v is the gas velocity, P is the gas pressure, and cs is the constant sound speed, ϵ = ρd/ρg is the dust-to-gas mass ratio, and τf is the friction (or stopping) time of solid bodies, respectively. The gas pressure gradient term is separated from the global one βcsΩ. We assumed that ρgrq, then the global pressure gradient term is (4)where H = cs/ Ω is the disk scale height. In our local model, βqH/r = −0.04 or 0 was adopted. The value β = −0.04 is slightly lower than the value at r = 1 AU in MMSN. Since we set up a pressure bump with a size similar to that found by Kato et al. (2012), we adopted the same value of β as they. On the other hand, β = 0 corresponds to the case without a global pressure gradient. Comparing our findings with the β = 0 result, we can define the effects of the global pressure gradient. The last term on the right-hand side of Eq. (2) is the back-reaction from dust onto gas. The description w is the allocated dust velocity at the grid point. We used the cloud-in-cell (CIC) model to allocate the average dust velocity. We assumed that the scaled stopping timescale Stokes number τsτfΩ is 1.0. Dust particles are only characterized by the stopping timescale. With τs = 1, the dust radial drift is fastest (Nakagawa et al. 1986).

We included solid materials as super-particles. Each super-particle consists of a large number of meter-sized boulders. The total number of super-particles is O(105−106) and each super-particle represents O(107) boulders. The equation of motion of the i-th super-particle is given by (5)where vi is the gas velocity at the location of the i-th particle, which is interpolated using v at the neighboring grid points.

In simulations, the gas equations are solved by the CIP scheme (Yabe & Aoki 1991), which is one of the methods of grid-hydrodynamics. The dust density and velocities are allocated to the closest four (2D) or two (1D) grid points in the computational region using cloud-in-cell interpolation. This algorithm strictly conserves angular momentum in the combined gas-dust system. We used a local simulation box, which means that the boundary conditions were periodic in all directions. For the radial boundary, however, we used the shearing-box approximation, which considers Keplerian differential rotation (Wisdom & Tremaine 1988; Hawley et al. 1995). These methods are similar to those used by Kato et al. (2012).

Table 1

Setup of individual runs.

2.2. Initial conditions

We conducted numerical simulations with three setups. The first was a 1D (radial) simulation with a bump. In this run, we tested the 1D evolution of the bump in the radial direction. The second was a 2D (radial-vertical) simulation with a bump. These are the main results of this paper. We investigated the pressure bump evolution and the growth of the SI at the same time. The third was a 1D simulation without the global pressure gradient. This setup was the same as the first, except that the global pressure gradient was included to check its effects on the time evolution of the radial pressure bump. Other parameters and settings for these cases are presented in Table 1.

In all types of calculations, we set the pressure bump as the initial condition. Our pressure bump was sustained by the azimuthal velocity profile given by Eq. (7). Therefore, it was steady if there was no back-reaction from solid materials onto gas. We note that we assumed that the pressure bump had already formed through some mechanism (e.g., MRI or baroclinic instability), and we followed the evolution of the bump without including its source mechanism in our simulation. This approximation assumes that the timescale of the bump structure production is longer than the deformation timescale that is due to back-reaction. The radial structure with equilibrium between dust and gas is found by solving a time-dependent initial-value problem.

The initial pressure bump is given by (6)where a and h are the arbitrary parameters for the height and width of a pressure bump. We set in our simulations. The initial azimuthal velocity profile is (7)which is derived from the radial force balances. In the 2D simulation, we assumed a uniform distribution vertically because of the narrow computational boxes here.

In all cases, the dust particles initially rotated with Keplerian velocity in the whole computational region. The initial particle distribution was uniform in our computational domain for simplicity, and the initial dust-to-gas mass ratio was ϵ = 0.1. The initial dust-to-gas mass ratio ϵ = 0.1 is higher than 0.01, which is the column dust-to-gas mass ratio in MMSN. Since we considered the solid materials to be meter-sized boulders, the scale height of the solid materials is more than 10 times smaller than the gas scale height (Youdin & Lithwick 2007). Therefore the dust-to-gas ratio is ~ 0.1 at the disk midplane. After a few frictional timescales, the dust particles attain terminal drift velocity given by the Nakagawa-Sekiya-Hayashi equilibrium solution (Nakagawa et al. 1986), except near the region of the bump, which is given by (8)To ensure that radial migration caused by gas drag provides a continuous supply of dust particles to the bump until run time Tend = 500, we set a sufficient radial width for the computational box. Since the radial drift velocity was ur,NSH = −0.018cs in our 1D simulation, the drift length is 9H. In our 2D simulation, the dust radial velocity decreases by ~ 40% because of the SI (Johansen & Youdin 2007), and the drift length is ~ 5.4H. Then the box size in the radial direction is Lx = 10H (1D simulations), 7.5H (2D simulation), which exceeds the drift length. The size of the vertical box is Lz = 0.25H, which is small enough compared to the disk scale height of the gas.

3. Results

3.1. 1D dust-density evolution

First, we describe the 1D view of the pressure bump evolution. The pressure bump halts inward, accreting dust particles at the point of highest pressure. If there is no back-reaction from dust to gas, the pressure bump is steady and the dust density increases until the gravitational collapse of solid components occurs, which is the expected result in a test-particle approach. When the dust drag is taken into account, the pressure bump is deformed by the angular momentum exchange and the dust particles are not concentrated at the point, but are distributed over a radial range. The timescale of the entire deformation of the pressure bump (i.e., there is no longer a super-Keplerian region) is ~ 150Ω-1. Since the radial dust accumulation causes the deformation of the pressure bump, this timescale depends on the global pressure gradient, the initial dust-to-gas ratio, and the stopping time of the particles. The radial width of the dust-dense region is about half that of the pressure bump. We assumed that the radial width of the pressure bump is ~ 2H. The width of the dust-dense region is ~ H.

Figure 1 shows the time evolution of the highest dust density, which is found in the cell with the highest density of the whole computational region. In Fig. 1 the green line represents the results of the run β04-1D. The highest dust density increases until tΩ ~ 10 because the dust is trapped in the pressure bump. At tΩ ≳ 10, the density stops increasing and the saturated dust density is ~ ρ0. The mechanism of this saturation is the same as with the 2D simulation (run β04-2D) and the simpler 1D simulation without global pressure (run β00-1D), which we discuss in subsequent sections.

The dust-dense region has a dust-to-gas mass ratio on the order of unity. Although a much higher dust-to-gas ratio is required for a simple Toomre-type self-gravitational instability to occur, Johansen & Youdin (2007) suggested that when two-dimensionality is introduced, the SI causes strong clumping of dust particles in such a dust-dense region. That is, a pressure bump is one of the most probable regions where the SI is excited and dust clumps are formed. In our simulations to be described in the next subsection, however, the behavior of SI is found to be different from previously studied cases because of the pressure bump deformation.

thumbnail Fig. 1

Time evolution of the highest dust density in run β04-2D (2D w/bump; red line) and β04-1D (1D w/bump; green dashed line). All the lines represent the dust density in the cell with the highest density of the whole computational region, normalized by the background gas density ρ0. The initial dust accumulation due to the bump structure is almost the same in the two runs (see inset).

Open with DEXTER

thumbnail Fig. 2

Dust-density pattern obtained from run β04-2D (2D w/bump). Each panel represents a snapshot at a) tΩ = 50, b)100; c)250; and d)500. Dust densities increase from black (zero density) to bright white (ρd = 2ρ0). Initially, the pressure bump concentrates the dust particles, with the dust-to-gas mass ratio being nearly unity (panel a)). Subsequent panels show the growth of SI in this dust-dense region. We note that the region near the pressure bump is zoomed-in and the whole computational domain is radially wider.

Open with DEXTER

3.2. 2D dust-density evolution

Based on 2D simulation results, we confirm that the formation of the dust-dense region is due to dust trapping and the growth of SI in the dust-dense region. In Fig. 1 the red line shows the time evolution of the highest dust density in the whole computational region from the 2D simulation. Until tΩ ~ 50, the highiest dust density evolves as in the 1D case. Then, the red line surpasses the green one, indicating that the linear growth of SI has begun. The dust concentration is saturated by tΩ ~ 200, with the highest dust density converging at ~ 5ρ0.

Figure 2 shows snapshots of the dust-density profile at the bump. In panel (a, b) of Fig. 2, which are the snapshots at tΩ = 50,100, the dust-dense region is formed by the radial pressure bump. The region has a nearly unity dust-to-gas mass ratio as in the 1D simulation. In panel (c), the density spatial pattern shows the linear growth of the SI. In panel (d), the saturated state of the instability is shown. There are some large dust clumps with densities higher than 2ρ0.

Based in Figs. 1 and 2, we confirm that SI forms dust clumps and that the highest dust density of these clumps is ~ 5ρ0 in the dust-dense region. On the other hand, Johansen & Youdin (2007) reported the highest dust density from clumping created by SI to be ~ 100ρ0 for τs = 1.0 and ϵ = 1.0. Johansen & Youdin (2007) employed an initially uniform dust and gas density of ~ ρ0 that was distributed throughout the whole computational domain. Our highest density is quite low compared with the value obtained from the other case, even though the dust-to-gas ratio that excites the SI is equally high ~ 1. This means that the SI with the spatial scale observed in previous studies is inefficient at the pressure bump. The reason for this weaker clumping in our case is discussed in subsequent sections.

Figure 3 shows a snapshot of the vertically averaged and highest values of dust and gas densities from the 2D result at tΩ = 500. The orange and magenta lines show the highest and vertically averaged dust density at each radial location, respectively. The deviation of the orange line from the magenta line shows that the SI is excited throughout the dust-dense region. The average dust density at the pressure bump is ~ ρ0 and is almost the same as the 1D results. This means that the basic dust-density radial structure is due to the 1D dynamics and that SI adds modulation along the vertical direction.

The blue and the cyan lines in Fig. 3 show the vertically averaged and highest gas densities at each radial location x, respectively. There is a small difference between the two, implying that SI does not result in a strong vertical variation of gas density.

3.3. Evolution of the gas pressure profile

In our simulations, we followed the time evolution of the local gas pressure, assuming that the global pressure gradient was constant. Here we study the evolution of the total pressure gradient force in the dust-dense region. The total gas pressure is a sum of the assumed constant global pressure Eq. (4) and the local gas pressure P calculated in our simulations.

Expressions for the total gas pressure and its gradient are Figure 4a shows snapshots of the total pressure at tΩ = 0,50,500 from the 2D case. At the pressure bump, the total pressure profile tends to become flat over time, and finally at tΩ = 500, the whole dust-dense region has an almost zero total pressure gradient. In Fig. 4b, to confirm pressure flattening, we check the evolution of the pressure gradient force in the dust dense region (−0.2H < x < 0.1H). The spatially averaged pressure gradient converges to ~ 0.0013, while the initial value of the global pressure gradient is | β | = 0.04. This means that the pressure gradient force at the dust-dense region is ~ 10−100 times lower than the value is initially set for the region outside the bump. As we discuss below, the significantly lower pressure gradient plays an important role in the size of the dense particle clumps formed by SI.

To investigate a simpler case without the global pressure gradient, we conducted a 1D simulation with β = 0. In this case, there was no accumulation of particles coming from regions outside the pressure bump, but dust was concentrated in the pressure bump toward the point of the pressure maximum. Figure 5 shows the snapshots of the gas profile evolution (tΩ = 50,250,500). Panels (a, b) shows the gas pressure profile and the gas azimuthal velocity profile (deviation from Keplerian rotation). Panel (c) shows the dust density profile. The system achieves steady-state by tΩ ~ 250. In the dust-dense region, the dust density is ~ ρ0, the gas azimuthal velocity is ~ 0, and the gas pressure profile is flattened. That is, the velocity distribution of gas is also forced to follow a Kepler rotation that is due to strong drag from dust particles in the area with a dust-to-gas ratio ~ 1. This clearly depicts the essential nature of what is seen in the 2D or 1D simulations with the global pressure gradient. The only difference is that in the simulations with a global pressure gradient, the completely flattened pressure profile (the dust-dense region) expands outward as a result of the continuous inward radial drift of dust.

thumbnail Fig. 3

Dust and gas density radial profile from run β04-2D (2D w/bump). Each line is a snapshot at tΩ = 500. Gray solid lines show the initial profile of gas (upper line) and dust (lower line) densities. The orange dotted line and the magenta solid lines show the highest and vertically averaged dust densities at each x. The vertically averaged dust density stays at ~ 1 as in the 1D run, while the upward deviation of the orange line from the pink one shows the growth of SI in the 2D case. The cyan solid and blue dotted lines showing the gas densities (highest and vertically averaged, respectively) show small differences that imply that the gas density has a nearly uniform distribution in the vertical direction.

Open with DEXTER

thumbnail Fig. 4

a) Vertically averaged gas pressure profile and b) time evolution of the radial pressure gradient force averaged spatially in the dust-dense region from run β04-2D (2D w/bump). As the dust particles accumulate at the bump, the total gas pressure profile is flattened. Then the pressure gradient force at the bump converges to ~ 0.0013, which is lower than the global pressure gradient β = −0.04 set initially (green dashed line).

Open with DEXTER

thumbnail Fig. 5

Results from the model without the global pressure support (run β00-1D). a) Gas pressure; b) gas azimuthal velocity (the deviation from the Keplerian rotation); and c) dust density. The essential physics of the dust-enriched region can be learned from these simple simulation results.

Open with DEXTER

4. Discussion

4.1. Schematic of the pressure bump evolution

In this subsection we discuss the mechanism for dust and gas density evolution at the radial pressure bump. First, we describe a 1D picture of the evolution in the radial direction. The 1D evolution can be separated into two parts. One is the saturation of the dust density that is due to the pressure gradient evolution. The second is the deformation of the entire pressure profile as a result of the continuous trapping of particles at the outer edge of the pressure bump from the outer part of the disk.

The saturation of dust density occurs when the dust and gas achieve equilibrium near the center of the radial pressure bump. Initially, the dust density grows faster than the gas density evolution. This is because the inward drift velocity is faster than the outward flow velocity of gas when the dust-to-gas mass ratio is lower than unity. As the dust density increases near the center of the pressure bump, the dust-to-gas mass ratio becomes ~ 1. Then the gas radial velocity becomes similar to the dust radial drift velocity, meaning that the timescale of the gas density evolution is equalized with the timescale of the dust density increase. The outward gas flow flattens the pressure bump and continues until the total pressure gradient force becomes ~ 0. Upon reaching the state where the pressure gradient force is ~ 0, the dust and gas have almost zero radial velocity and achieve steady-state. The equilibrium is such that both dust and gas are mostly in the Kepler rotation state. Since additional density concentration cannot be expected, the highest dust concentration saturates at a dust-to-gas mass ratio of ~ 1.

Regarding the second point of the density evolution mechanism at the radial pressure bump, the entire pressure bump deformation is brought about by the continuous dust trapping at the outer edge of the pressure bump. When there is continuous dust trapping, the gas component continuously leaks outward. This outward-flowing mechanism continues until the dust trapping ceases or the pressure bump is destroyed entirely.

The simpler 1D calculation without global pressure gradient helps us understand the bump evolution very clearly since the final stage does not change in time and thus the interpretation becomes easier. Since there is no global pressure gradient β = 0, there is no dust trapping from outside the pressure bump, but the dust particles that initially existed inside the pressure bump accumulate at the super- or sub-Keplerian transition point. The results shown in Fig. 5 indicate that the highest dust density of ρd,max ~ 1.4ρ0, which is the same as for β = −0.04, is obtained in the region where the gas pressure gradient is nearly 0. The flattened pressure profile in the dust-dense region is in a steady state of equilibrium. That is, the dust concentration stops when both gas and dust achieve Kepler rotation. A dust-to-gas mass ratio of ~ 1 corresponds to a scenario where the back-reaction becomes strong enough to force the gas into Kepler rotation. The size of the dust-dense region is ~ 0.3H. This is determined by the initial total mass of the particles in the pressure bump.

Results of the 2D simulation are explained by the 1D picture with the addition of vertical modulation due to SI. After the formation of a dust-dense region with the dust-to-gas ratio of ~ 1 that is due to the 1D dynamics, SI is excited because of the high dust density. We note that the pressure gradient in the dust-dense region is significantly reduced compared to initial conditions, which affects the properties of SI that could be excited in this environment, as discussed below.

4.2. Properties of SI in the dust-dense region

In this section, we discuss the detailed properties of SI excited in the dust-dense region formed by the radial pressure bump. The dust-dense region can be regarded as the initial conditions of SI. We found that in the dust dense region, (i) the dust-to-gas mass ration is ~ 1 and (ii) η is reduced by 10−100 times from the nominal value of MMSN.

The relatively high dust-to-gas ratio is favorable for SI (Youdin & Johansen 2007; Carrera et al. 2015). Especially when the stopping timescale τs is shorter than ~ 1, the pressure bump may be important for SI to occur. The bump deformation mechanism does not strongly depend on the frictional timescale. If we neglect gas turbulence, the bump deformation process becomes slower, but the saturation mechanism of dust density does not change even for millimeter-sized particles. While the dust-to-gas mass ratio also becomes ~ 1 for small particles in the dust-dense region, the growth rate of SI rises drastically at a dust-to-gas mass ratio 1 (Johansen & Youdin 2007). Therefore, the pressure bump may be a good location for planetesimal formation, even for small dust particles that have a short stopping time of τs ~ 0.003 corresponding to the mm-sized object at r = 1 AU in MMSN.

The evolution of pressure gradient, however, affects the properties of SI at the pressure bump. The spatial scales of SI are normalized by ηr, where η is another expression of the global pressure gradient η = −/ (2r) (Youdin & Goodman 2005). As seen in Sect. 3.3, the total pressure gradient at the pressure bump becomes ~ 10−100 times lower than the initial global value. Accordingly, the spatial scale of SI also becomes ~ 10−100 times smaller than the value of previous studies. This means that SI with the spatial scale observed in previous studies has a reduced growth rate. We did not observe dense particle clumps with a size similar to those found in Johansen & Youdin (2007). The most unstable scale of dense particle clumps at the pressure bump predicted by our results is also 10−100 times smaller than that given by Johansen & Youdin (2007). Johansen et al. (2012) conducted simulations including many effects (e.g., vertical gravity, particle self-gravity, and particle collision) to find that the SI forms gravitationally bound clumps corresponding to ~ 100−1000 km sized bodies with a normal MMSN value of β ~ 0.1. If the SI occurs in the same way as reported by Johansen et al. (2012) at the pressure bump, the predicted clump size is 10−100 times smaller than their results. The predicted size is ~ 1−100 km sized bodies, which is similar to the classically predicted planetesimal size.

We note that the spatial resolution of our calculations is insufficient to resolve this small-scale SI in the dust-dense region. We assumed that the width of the radial pressure bump is similar in size to the gas pressure scale height. To ensure a sufficient width for the dust drift and accumulation, our computational box has ~ 10 scale heights in the radial direction. However, the predicted short wavelength of the SI in the dust-dense region requires higher resolution. Since the wavelength of the SI becomes short but the growth rate does not change even if the global pressure gradient becomes small, a simple prediction is that the highest dust density is similar to the results of the case with τs = 1.0 = 1.0 in Johansen & Youdin (2007), which is almost the same setting as with the dust-dense region formed by the pressure bump. Thereby, it is not clear in our simulations whether the SI with a short wavelength does occur.

4.3. Future works on planetesimal formation at the radial pressure bump

Our calculations did not take the Kelvin-Helmholtz instability (KHI) by radial shear in azimuthal velocity and vertical shear in radial or azimuthal velocity into account. Johansen & Youdin (2007) and Bai & Stone (2010) found that SI dominantes the KHIs. In general, small-scale instability can be damped by other perturbations. At the radial boundary of the dust-dense region, there must be a shear stronger than the Keplerian one to cause strong KHI. The investigation of the effect of KHI is left for future work.

Dust self-gravity and vertical gravity of the host star are other important physics neglected in this study. Bai & Stone (2010) and Johansen et al. (2012) included them and confirmed that gravitationally bound clumps are formed by SI at the midplane dust layer. This may also be applied to the pressure bump.

We set the radial pressure bump as the initial condition of our calculations. Various formation mechanisms for pressure bumps have been proposed (e.g., MRI zonal flows, edges of the MRI dead zone, and a snow line). These mechanisms and the global disk evolution affect the bump deformation process. If the bump restoration or restructuring process is faster than the destruction process we discussed here, the radial pressure bump could accumulate many more dust particles.

5. Summary

We numerically studied the dust and gas density evolution at a radial pressure bump in a protoplanetary disk. Disk gas dynamics was treated with grid hydrodynamics and dust components with the particle-in-cell scheme. Computational regions were approximated by the local shearing-box model. The gas drag force to dust and its back-reaction to gas were consistently included in our simulations.

Our conclusions are summarized as follows:

  • The pressure bump accumulated dust particles, but the outwardflow of gas created by the back-reaction smoothed it out after thedust-to-gas mass ratio reached ~ 1. However, this only occurred if bump restoration or restructuring was slower than the destruction process.

  • The pressure bump created a long-lived dust-dense region where the vertically averaged dust-to-gas mass ratio was ~ 1. When the inward drift of dust continued to deliver dust particles to the outer edge of the pressure bump, the flattened pressure gradient region in the dust-dense region expanded outward until the pressure bump was completely destroyed.

  • As a result of bump destruction, direct self-gravitational instability at the site of the pressure bump was inhibited.

  • The gas flow in pressure bumps was modulated to be close to the Keplerian one and the deviation from Keplerian velocity (η) was ~ 10−100 times smaller than the global value. Because the SI wavelength is scaled by η, the clumps formed by SI, even if the SI occurs, should become 1−100 km sized bodies, which are much smaller than the previously proposed size (~ 100−1000 km). This means that the SI forms planetesimals with sizes similar to those formed by the self-gravitational instability of a thin dust layer (Safronov 1969; Goldreich & Ward 1973) at the radial pressure bump.

Acknowledgments

We thank an anonymous referee for detailed comments. We are grateful to Taishi Nakamoto, Taku Takeuchi, Satoshi Okuzumi and Akimasa Kataoka for comments. This research was supported by a grant for JSPS (23103005) Grant-in-aid for Scientific Research on Innovative Areas.

References

All Tables

Table 1

Setup of individual runs.

All Figures

thumbnail Fig. 1

Time evolution of the highest dust density in run β04-2D (2D w/bump; red line) and β04-1D (1D w/bump; green dashed line). All the lines represent the dust density in the cell with the highest density of the whole computational region, normalized by the background gas density ρ0. The initial dust accumulation due to the bump structure is almost the same in the two runs (see inset).

Open with DEXTER
In the text
thumbnail Fig. 2

Dust-density pattern obtained from run β04-2D (2D w/bump). Each panel represents a snapshot at a) tΩ = 50, b)100; c)250; and d)500. Dust densities increase from black (zero density) to bright white (ρd = 2ρ0). Initially, the pressure bump concentrates the dust particles, with the dust-to-gas mass ratio being nearly unity (panel a)). Subsequent panels show the growth of SI in this dust-dense region. We note that the region near the pressure bump is zoomed-in and the whole computational domain is radially wider.

Open with DEXTER
In the text
thumbnail Fig. 3

Dust and gas density radial profile from run β04-2D (2D w/bump). Each line is a snapshot at tΩ = 500. Gray solid lines show the initial profile of gas (upper line) and dust (lower line) densities. The orange dotted line and the magenta solid lines show the highest and vertically averaged dust densities at each x. The vertically averaged dust density stays at ~ 1 as in the 1D run, while the upward deviation of the orange line from the pink one shows the growth of SI in the 2D case. The cyan solid and blue dotted lines showing the gas densities (highest and vertically averaged, respectively) show small differences that imply that the gas density has a nearly uniform distribution in the vertical direction.

Open with DEXTER
In the text
thumbnail Fig. 4

a) Vertically averaged gas pressure profile and b) time evolution of the radial pressure gradient force averaged spatially in the dust-dense region from run β04-2D (2D w/bump). As the dust particles accumulate at the bump, the total gas pressure profile is flattened. Then the pressure gradient force at the bump converges to ~ 0.0013, which is lower than the global pressure gradient β = −0.04 set initially (green dashed line).

Open with DEXTER
In the text
thumbnail Fig. 5

Results from the model without the global pressure support (run β00-1D). a) Gas pressure; b) gas azimuthal velocity (the deviation from the Keplerian rotation); and c) dust density. The essential physics of the dust-enriched region can be learned from these simple simulation results.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.