Including Dust Coagulation in Hydrodynamic Models of Protoplanetary Disks: Dust Evolution in the Vicinity of a Jupiter-mass Planet

All the 2D (r + ϕ) hydrodynamic models are performed using a code that is a part of LA-COMPASS (Los Alamos CoMPutational AStrophysics Suite) and was described by Li et al. (2005, 2009). The protoplanetary disk is assumed to be geometrically thin so that the hydrodynamical equations can be reduced to two-dimensional Navier–Stokes equations by considering vertically integrated quantities. We adopt a locally isothermal equation of state:

$$\begin{eqnarray}&&{P}_{{\rm{g}}}(r)={{\rm{\Sigma }}}_{{\rm{g}}}{c}_{{\rm{s}}}^{2},\end{eqnarray} \tag{ 3 }$$

where P_g is the vertically integrated pressure, Σ_g is the gas surface density, and c_s is the sound speed, which only depends on distance r here (see Equation (2)).

Dust is treated as a pressureless fluid in a bi-fluid model that is governed by conservation laws (see, e.g., Paardekooper & Mellema 2006; Meheut et al. 2012). The gas and dust equations are coupled together through source terms that model the drag between the two fluids, i.e., including the backreaction from dust to gas. We include both Epstein and Stokes drag regimes. We have implemented a Godunov Riemann solver for dust equations. We also implement dust diffusion due to the turbulence and consistently combine it with the dynamic model of bi-fluid (see Fu et al. 2014b; Li & Li 2014). To deal with multiple timescales of coupling different dust species with gas dynamics, we develop an efficient and robust L-stable method to solve the coupled gas and dust equations (Li 2017).

The planet motion is governed by Newton's laws, whose equations are solved with an adaptive high-order Runge–Kutta method. The planet's gravitational potential is smoothed over 0.7 disk scale height H_g at the planet location. The disk self-gravity is not included in the models presented in this paper.

We used a linear polar grid with uniform spacing between the cells. The grid resolution is ${N}_{{\rm{r}}}\times {N}_{\phi }=1024\times 1024$. With this resolution, the disk scale height H_g is resolved with 17 cells and the Hill radius of the planet is resolved with 23 cells at the planet location.

The computational domain is set from 4 to 34 au. We keep the gas density constant at the inner and outer boundaries. This is justified because no significant viscous evolution is expected for the duration of the simulations (∼10⁵ yr). For dust, we have an open inner boundary to allow outflow. We keep the dust density at the outer boundary constant for the first 1000 planet orbits, which is equivalent to a steady inflow, and then we close the boundary, mimicking a decreasing flux of dust from the outer disk to the planet region. In the initial condition, velocities of gas and dust are set according to their equilibrium values derived by Nakagawa et al. (1986). At the inner and outer boundaries, the velocities (both radial and azimuthal) are kept constant during the simulation.

2.1.1. Fixed-size Dust

2.1.2. Full Dust Coagulation Treatment

The default code has been modified to include multiple dust fluids representing different dust sizes. Collisional evolution of dust is solved using explicit integration of the Smoluchowski equation. We include the Brownian motion, turbulence, differential radial and azimuthal drift, and vertical settling as sources of the collision velocities. The values of radial and azimuthal velocities for each dust species are obtained directly from the hydrodynamic solver and the other three sources are calculated using the same method as Birnstiel et al. (2010).

When calculating the collision probabilities, we take into account the midplane density of dust, which we calculate for each dust species i from the surface density ${{\rm{\Sigma }}}_{{\rm{d}}}^{i}$ used in the 2D version of LA-COMPASS:

$$\begin{eqnarray}&&{\rho }_{{\rm{d}}}^{i}=\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}^{i}}{\sqrt{2\pi }{H}_{{\rm{d}}}^{i}},\end{eqnarray} \tag{ 4 }$$

where we assume a Gaussian distribution of the grains around the midplane with scale height following Dubrulle et al. (1995):

$$\begin{eqnarray}&&{H}_{{\rm{d}}}^{i}={H}_{{\rm{g}}}\cdot \sqrt{\displaystyle \frac{\alpha }{\alpha +{\mathrm{St}}^{i}}},\end{eqnarray} \tag{ 5 }$$

where H_g is the gas scale height. For small grains, this equation is consistent with the work of Youdin & Lithwick (2007). In this approach, we assume that turbulent mixing is fast enough to always keep the vertical structure in the settling-mixing equilibrium. This assumption might break in a low turbulence case, when the interplay between settling and dust growth leads to the so-called sedimentation-driven coagulation (see, e.g., Zsom et al. 2011).

We assume that grains are compact spheres with internal densities of ρ_s = 1.2 g cm⁻³. Collisional outcomes include sticking for collisions with impact speeds below v_f = 10 m s⁻¹, fragmentation for collisions speeds above v_f, and erosion for collisions speeds above v_f when the mass ratio of colliding particles is greater than 10. Numerical implementation of the collisional evolution is the same as described in Birnstiel et al. (2010, their Section 2.6).

The dust-size distribution is resolved with 151 dust fluids covering sizes between 1 μm and 100 cm, which corresponds to 8.4 grid points per mass decade, a typical resolution used in dust coagulation models. In the initial condition, all the dust has a radius of 1 μm.

Due to the computational expense of solving dust coagulation, we call the coagulation solver every 50 time steps of the hydrodynamic solver. We tested that this sub-stepping routine does not impact the results significantly by running an analogical simulation where coagulation was solved at every time step (but with a shorter duration).

A more detailed description of the code will be given in a corresponding paper by S. Li et al. (2019, in preparation).

2.1.3. Simple Dust Coagulation Approach

We implemented a simple, sub-grid method for dust growth in the LA-COMPASS code, which is an updated version of the method adopted by Tamfal et al. (2018). In this method, dust is treated as a single fluid but its size is not fixed. The dust size is calculated at every time step and in every cell based on local conditions. In the initial condition, the size in all cells is set to a₀ = 1 μm. The initial Stokes number of grains is calculated as

$$\begin{eqnarray}&&{\mathrm{St}}_{0}=\displaystyle \frac{\pi }{2}\displaystyle \frac{{a}_{0}{\rho }_{{\rm{s}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}},\end{eqnarray} \tag{ 6 }$$

as all the micron-sized grains in our computational domain are in the Epstein drag regime. Dust growth is modeled as

$$\begin{eqnarray}&&{a}_{i}={a}_{i-1}+\dot{a}\cdot {\rm{\Delta }}t\end{eqnarray} \tag{ 7 }$$

where ${a}_{i-1}$ is dust size obtained in the given cell in the previous time step, Δt is the length of the time step, and the growth speed $\dot{a}$ is calculated as (see, e.g., Kornet et al. 2001)

$$\begin{eqnarray}&&\dot{a}=\displaystyle \frac{{\rho }_{{\rm{d}}}{\rm{\Delta }}v}{{\rho }_{{\rm{s}}}},\end{eqnarray} \tag{ 8 }$$

where ρ_d is the midplane density of dust (see Equation (4)), Δv is impact velocity between grains of Stokes number equivalent to St_i−1 and $0.5\cdot {\mathrm{St}}_{i-1}$ (where St_i−1 corresponds to the Stokes number of grains with size ${a}_{i-1}$), and ${\rho }_{{\rm{s}}}$ is the internal density of grains. When calculating the impact speed Δv, we take into account turbulence, radial drift, and azimuthal drift. The impact speeds are calculated from the radial and azimuthal velocities returned by the hydrodynamic solver assuming that the radial speed depends on the Stokes number as v_r ∝ St (correct for St < 1) and the azimuthal speed as ${v}_{\phi }\propto 1/(1+{\mathrm{St}}^{2})$.

Dust growth can be halted either by fragmentation or radial drift. To take this into account, we calculate the maximum dust size that could grow using the semi-analytic expressions derived by Birnstiel et al. (2012). The maximum Stokes number with turbulence-driven fragmentation is

$$\begin{eqnarray}&&{\mathrm{St}}_{\mathrm{frag}}={f}_{{\rm{f}}}\cdot \displaystyle \frac{{v}_{{\rm{f}}}^{2}}{3\alpha {c}_{{\rm{s}}}^{2}},\end{eqnarray} \tag{ 9 }$$

where the fudge factor ${f}_{{\rm{f}}}=0.37$, and v_f is the fragmentation threshold velocity, which we set to 10 m s⁻¹ in this paper. This equation was derived assuming that the turbulence-driven impact velocities scale as ${\rm{\Delta }}{v}_{{\rm{t}}}\propto \sqrt{\mathrm{St}}$, which applies for grains in the so-called fully intermediate regime of Ormel et al. (2007). In fact, grains that hit the fragmentation barrier are typically in this regime as shown by Birnstiel et al. (2011). Fragmentation can also be caused by the differential drift, and the maximum Stokes number for the drift-induced fragmentation is

$$\begin{eqnarray}&&{\mathrm{St}}_{\mathrm{df}}={f}_{{\rm{f}}}\cdot \displaystyle \frac{{v}_{{\rm{f}}}}{| \eta | {v}_{{\rm{K}}}},\end{eqnarray} \tag{ 10 }$$

where ηv_K is the maximum drift speed calculated using the midplane radial pressure gradient

$$\begin{eqnarray}&&\eta {v}_{{\rm{K}}}=\displaystyle \frac{1}{2{\rho }_{{\rm{g}}}{{\rm{\Omega }}}_{{\rm{K}}}}\cdot \displaystyle \frac{{{dP}}_{{\rm{g}},z=0}}{{dr}},\end{eqnarray} \tag{ 11 }$$

where ρ_g is the midplane gas density and P_g,z=0 is midplane gas pressure.

Fragmentation is the dominant factor in setting dust size when the coagulation timescale is shorter than the drift timescale. However, in a realistic disk, this is not always true. Particularly, in the outer part of the disk, radial drift may be faster than coagulation. This sets a limit on how far the growth can proceed before the grains are removed faster than they can grow. This effect is naturally recovered in the full-coagulation models, in which each dust fluid can be advected at its own speed. Although the radial drift is still accurately modeled by the hydrodynamic solver, in the simple coagulation approach the advection of dust does not have a direct effect on the representative size. Therefore, we must include the drift limit explicitly in the size calculation. The maximum Stokes number that can remain at a given location taking into account radial drift is

$$\begin{eqnarray}&&{\mathrm{St}}_{\mathrm{drift}}={f}_{{\rm{d}}}\cdot \displaystyle \frac{1}{2| \eta | }\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}},\end{eqnarray} \tag{ 12 }$$

where the fudge factor f_d = 0.55.

The values of the fudge factors ${f}_{{\rm{f}}}$ and ${f}_{{\rm{d}}}$ that we adopted were derived by Birnstiel et al. (2012) by comparing simple coagulation results to 1D coagulation in the framework of a global, smooth disk.

The new Stokes number is decided by choosing the minimum of the values calculated when taking into account growth (St_i, corresponding to the size a_i obtained in Equation (7)), and the possible barriers (Equations (9), (10), and (12)):

$$\begin{eqnarray}&&\mathrm{St}=\min \left({\mathrm{St}}_{i},{\mathrm{St}}_{\mathrm{frag}},{\mathrm{St}}_{\mathrm{df}},{\mathrm{St}}_{\mathrm{drift}}\right).\end{eqnarray} \tag{ 13 }$$

We found that, particularly in cases where the pressure gradient is briefly enhanced by the spiral wakes (see Figure 3), the Stokes number recovered from this treatment can be much lower than the one given by full-coagulation results. This is because the radial advection of size and the timescale needed to fragment all particles in a given cell are not taken into account. To minimize this effect, we limit the impact of fragmentation by comparing the Stokes number obtained in Equation (13) to the Stokes number from the previous time step, ${\mathrm{St}}_{i-1}$, and if $\mathrm{St}\lt {\mathrm{St}}_{i-1}$ we set

$$\begin{eqnarray}&&\mathrm{St}=\min \left(1,{f}_{{\rm{n}}}\right)\cdot \mathrm{St}+\max \left(0,(1-{f}_{{\rm{n}}})\right)\cdot {\mathrm{St}}_{i-1},\end{eqnarray} \tag{ 14 }$$

where the fudge factor ${f}_{{\rm{n}}}={\rm{\Delta }}t/{t}_{\mathrm{coag}}$ is a ratio of the simulation time step ${\rm{\Delta }}t$ and the coagulation timescale calculated as ${t}_{\mathrm{coag}}=a/\dot{a}$ (see Equation (8)). This way we avoid any sudden, local drops of the Stokes number but let it decrease gradually. Finally, we limit the minimum value of the Stokes number to the Stokes number of the smallest, micron-sized grains:

$$\begin{eqnarray}&&\mathrm{St}=\max \left(\mathrm{St},{\mathrm{St}}_{0}\right).\end{eqnarray} \tag{ 15 }$$

The main difference between our implementation of the simple coagulation and the algorithm presented by Tamfal et al. (2018) is in the treatment of the initial dust growth phase. We changed the dust growth prescription from an exponential function implemented by Tamfal et al. (2018) to calculating the growth rate based on the local conditions (see Equation (8)). We have also introduced a limit on how much the size can decrease between two consecutive time steps (see Equation (14)).

To run the azimuthally averaged models, we used the DustPy code. The code, developed by S. M. Stammler and T. Birnstiel, is a Python-based version of the commonly used dust coagulation code described by Birnstiel et al. (2010). It solves dust coagulation and radial surface density evolution in an azimuthally and vertically averaged framework, performing implicit integration of the Smoluchowski equation and the advection-diffusion equation.

To test the impact of solving dust coagulation in an azimuthally averaged framework, and reproduce the approach previously used to study dust coagulation in the presence of a gap-opening planet, we set up a model where the DustPy code uses azimuthally averaged gas evolution obtained in the LA-COMPASS simulation as an input. This is done in the following way: the azimuthally averaged output of the 2D model is stored at every 10 planet orbits. An interpolation routine is used to generate input at time instances needed by the 1D model. Otherwise, we use the coagulation setup with the parameters given in Table 1.

As expected, the massive planet placed in a viscous disk quickly clears a gap; however, some accretion flow through the gap is retained (Lubow & D'Angelo 2006), which is visible in Figure 1. The initial power-law density profiles are modified not only by planetary gap opening, but also by the planet-induced spiral density waves and by dust drift. Dust evolution depends strongly on grain sizes, which influence their aerodynamic interaction with gas.

Zoom In Zoom Out Reset image size

We run the simulation for 4000 planet orbits (corresponding to t = 1.26 × 10⁵ yr). Figure 2 presents the time evolution of the gas and dust surface density, as well as characteristic dust size. As can be seen, the gap profile is practically saturated at the end of the simulation. As the dust growth timescale is of the order of 100 local orbits (or, to be more precise, the growth timescale is orbital timescale divided by the dust-to-gas ratio; see, e.g., Birnstiel et al. 2012), the dust sizes quickly reach a steady state, therefore they do not change significantly after the gap is fully open.

Zoom In Zoom Out Reset image size

Interaction between the planet and the disk causes formation of a pressure bump at the outer edge of the planet gap. Figure 3 shows the η parameter (see Equation (11)), which defines the maximum possible drift speed of dust and its direction. The red line corresponds to azimuthally averaged value of η, which determines the overall evolution of dust. Negative values of η translate into inward drift and positive values mean outward drift of dust. As η > 0 between 11 and 15 au, the inward drift is reversed, thus we can expect that the radial dust drift is halted around 15 au. Indeed, as visible in the middle panel of Figure 2, dust is radially concentrated around 15 au. However, trapping is not 100% efficient and the inner region of the domain is not completely cleared, but some population of grains is retained throughout the simulation. This is because the effect of trapping is compromised by viscosity (see, e.g., Pinilla et al. 2012; Ataiee et al. 2018). Because the radial drift speed in the pressure bump is directed toward η = 0 and it increases with size, there is a critical size for which the particle will always drift back to the pressure bump after being displaced by random turbulent movements. Thus, small grains are expected to pass through the gap, while large grains are expected to stay outside of the gap. Gas flows quickly through the gap and small grains, which are well-coupled, are carried along (Rice et al. 2006; Zhu et al. 2012). This is confirmed in the lower panel of Figure 2, where the density-averaged size at each location is plotted. The typical size of grains in the gap is much smaller than that outside of the gap. We will discuss this effect in more details in the subsequent section, where we focus on comparing the full-coagulation run to models employing the fixed-size approach.

Zoom In Zoom Out Reset image size

It is worth noting that the radial profiles of the η parameter plotted for different disk sectors (the gray lines in Figure 3) display significant variations, driven mostly by the spiral wakes. These wakes sweep the disk, causing temporary, small-scale pressure bumps. However, in the full-coagulation approach, these do not seem to modify dust evolution considerably.

Figure 4 compares dust distribution in the protoplanetary disk obtained in the series of models assuming fixed dust size and in the model with full coagulation. In agreement with previously published results (see, e.g., Paardekooper & Mellema 2004; Rice et al. 2006; Fouchet et al. 2010; Weber et al. 2018), we find that large grains are trapped outside of the gap opened by the planet and cannot pass it. The larger the grains, the larger and deeper the gap they open in dust density distribution. On the other hand, small grains that are coupled to the gas pass through the gap. The critical size of grains that can pass through the gap is about 1 mm in our setup. The millimeter-sized grains open a clear gap but at the same time do not form a distinct peak outside of the planetary gap, which is characteristic for simulations including larger dust sizes.

Zoom In Zoom Out Reset image size

Figure 5 compares the azimuthally averaged dust density profiles obtained in the series of fixed-size models (upper panel) to the results of the full-coagulation model (lower panel). We find that the results obtained when applying the full dust coagulation treatment cannot be adequately fitted using a single fixed-size model. Dust distribution resulting from the interplay between multi-size dust advection and coagulation shows both confined peak outside of the planetary gap, characteristic for centimeter-sized and larger grains, and the partially filled planetary gap, characteristic for submillimeter grains. Overall, the slope of dust density through the outer edge of the planet gap is much shallower in the full-coagulation model than in most of the fixed-size models.

Zoom In Zoom Out Reset image size

This can be understood when considering what are the contributions to the total dust density from grains of different sizes (see the lower panel of Figure 5). While the maximum dust sizes that can be obtained in the pressure bump outside of the planet orbit are on the order of few centimeters, the gap is filled exclusively with grains smaller than 300 μm.

In models that do not solve dust coagulation but need input on size distribution (mostly for comparison to observations), it is often assumed that the dust size distribution follows the so-called MRN distribution with $n(a)\propto {a}^{-3.5}$ (Mathis et al. 1977). In the upper panel of Figure 5, we summed up contributions from each single-size model assuming the MRN distribution. The resulting total density profile (red dashed line) is relatively similar to the density profile obtained in the full-coagulation model presented in the bottom panel. In Figure 6 we compare the global size distribution obtained in the full models to the assumed MRN profile. At 1000 orbits of the planet, it does indeed match the power-law distribution reasonably well, although the slope of the distribution is generally shallower and there is significantly less of the largest grains. The MRN profile may be a reasonable assumption for the overall size distribution, although a good estimate of the maximum dust size is necessary, as most of the mass is contained in the largest grains. The size distribution is significantly different at the beginning of the simulation, when the grains have not reached their maximum sizes yet. Also, the maximum size of grains decreases toward the end of the simulation, and falls from 10 cm to about 3 cm at 4000 planet orbits. This is caused by the decreasing dust influx at the outer boundary (see Section 2.1). With a lower dust-to-gas ratio, the drift barrier affects smaller grains (Birnstiel et al. 2012). Thus, the MRN size distribution with a fixed maximum dust size may not be very useful to model disk evolution over a long period of time.

Zoom In Zoom Out Reset image size

A significant difference between the full-coagulation and fixed-size models is that in the full models fragmentation constantly replenishes small grains that can pass through the planetary gap. The impact of fragmentation is noticeable when comparing the upper and lower panels of Figure 5. In the full-coagulation simulation (lower panel), the small grains follow the density profile of larger grains outside of the planet orbit, differing from what is expected from the fixed-size simulations (upper panel). This is because the small grains are constantly produced in collisions between larger aggregates. Our results suggest that, if fragmentation is efficient, the density of small aggregates should be enhanced in dust traps, despite them not being trapped themselves. Indeed, a look into the size distribution of grains presented in Figure 7 reveals that outside of the planet orbit, the size distribution is close to coagulation-fragmentation equilibrium (see, e.g., Birnstiel et al. 2011).

Zoom In Zoom Out Reset image size

Another effect that distinguishes the full-coagulation simulation from the fixed-size models is the growth of small particles after they passed the gap. Dust density in the gap is too low to allow for efficient coagulation. But in the inner region of the simulation domain, where dust density increases, grains larger than 3 cm are present, which would not be expected from the fixed-size models.

It is worth noting that the dust distribution obtained outside of the planet gap (the upper panel of Figure 7) is remarkably symmetric, with little deviation when considering different azimuthal angles. This is not true in the planet co-orbital region (lower panel of Figure 7). Some of the small grains that pass through the gap are trapped either in the direct vicinity of the planet (we do not consider accretion onto the planet) or in the Lagrange points (this is visible in Figure 1). Due to this asymmetric nature of the co-orbital region, the size distributions sampled around the planet orbit exhibit different profiles.

Because the density and size distribution profiles are generally symmetric (outside of the planetary gap region) and the size distribution generally follows the expected coagulation-fragmentation equilibrium profile, we can expect that it would be possible to recover results of the full-coagulation model using less computationally expensive methods. In the subsequent section, we compare results of three methods including dust coagulation.

We aimed to reproduce the results of the full-coagulation run using two less computationally intensive methods. The first method relies on a semi-analytical estimate of the coagulation outcome and is based on the work of Birnstiel et al. (2012). The second method relies on extracting the gas evolution from the full 2D simulation and running a 1D azimuthally averaged dust evolution model in a post-processing step. The estimated computational expense of the three methods, per 1000 planet orbits, is as follows:

• 1.  
  Full coagulation: 27 684 CPU hours (12 hr on 2304 CPUs).
• 2.  
  Simple coagulation: 192 CPU hours (20 minutes on 576 CPUs, very similar to an analogical fixed-size simulation).
• 3.  
  1D coagulation: 78 CPU hours.

Figure 8 compares azimuthally averaged dust densities and dust sizes obtained using the three methods. As can be inferred from this plot, the simple coagulation method is generally better at reproducing the full results than the azimuthally averaged approach.

Zoom In Zoom Out Reset image size

The main problem with 1D coagulation is that azimuthal averaging of gas density "kills" the effects of spiral wakes (see Figure 3): they induce additional impact speeds, limiting the maximum size possible to obtain, but they also induce extra mixing. The 1D coagulation predicts almost one order of magnitude higher peak density in the trap outside of the planet orbit. The peak is also narrower than that in full-coagulation results. This effect has multiple reasons: (1) due to averaging out of the spiral wakes, the 1D coagulation predicts larger particles, that are trapped more efficiently, and (2) in the full-coagulation results, the exact position of the trap may change with the azimuthal angle because of the asymmetric nature of planet–disk interaction. Thus, the trap is radially "smeared" in the full-coagulation results. Additionally, as we mentioned before, the spiral wakes induce additional mixing that is not taken into account in the 1D model. The full-coagulation model includes the effect of backreaction of dust to gas, which additionally increases the width of the dust ring (although this is not a significant contribution in our setup; see Section 4.3). It is worth noting that similar results for widening the peaks in 2D versus 1D models were found by Weber et al. (2018) for fixed-size grains.

Neglecting the 2D effects of planet–disk interaction has the most significant outcome in the planetary gap region. The 1D coagulation predicts significantly more dust inside of the gap than the full and simple coagulation models. In 2D simulations, dust can only flow through the gap if it enters the streamline around the plane (see Figure 1). The 1D model cannot take this subtlety into account. In our case, using the azimuthally averaged gas density information to calculate pressure gradient and drift speed leads to an increased dust flux through the gap. This is the opposite conclusion of presented by Weber et al. (2018), who also compared a 2D and 1D approach to dust dynamics in the vicinity of Jupiter-mass planet, using a fixed-size model. However, they did not extract the density information from a 2D simulation, but ran a self-consistent 1D model with the planet, which led to a different density profile in the 1D and 2D runs.

The simple coagulation results do not reproduce the full coagulation perfectly either. The main problem of the simple model is that the size is calculated locally, without input from neighboring cells, thus the effect of dust mixing is not taken into account in all aspects. As shown by Birnstiel et al. (2012), the simple method reproduces results of full coagulation very well in the case of a smooth, axisymmetric disk. We find that including a massive planet that induces spiral wakes, locally enhancing the pressure gradient and thus impact speeds, leads to violent fragmentation events. In the current simple coagulation approach, we only track one "representative" particle size per cell. If this size suddenly drops due to the spiral-wake-induced fragmentation, dust will take a relatively long time to regrow at this position. In the full-coagulation method, this effect is significantly reduced as the fluids representing different sizes mix, leading to a similar size distribution in neighboring cells. This is why we had to limit the effect of fragmentation in the simple coagulation model (see Section 2.1.3).

Despite these difficulties, the simple model reproduces the full-coagulation results reasonably well. Outside of the region where the spiral wakes have the strongest effect (∼10 to 25 au, the dust size calculated in the sub-grid method fits the density-averaged size obtained in the full-coagulation model almost perfectly (see the lower panel of Figure 8). It is worth noting that the implementation of the simple coagulation practically does not increase the computational cost of the 2D hydrodynamic model, so this calculation is as fast as a fixed-size run.

We presented results of computational models utilizing state-of-the-art methods for modeling dust evolution in protoplanetary disk. However, our models are not free from limitations, which we discuss in this section.

We performed 2D models, solving for radial and azimuthal structure of the protoplanetary disk, assuming that the vertical density distribution is Gaussian, and depends on dust size in a simple way (see Equations (4) and (5)). Thus, we neglect the potential effect that sedimentation-driven coagulation could have on dust growth (Zsom et al. 2011). It is known that in some cases, the 3D effects may change the conclusions of 2D hydrodynamic models (see, e.g., Lyra et al. 2018).

We adopted a simple, isothermal protoplanetary disk model with a fixed temperature structure. Thus we do not take into account the effects of a planet heating the protoplanetary disk (Pohl et al. 2015; Szulágyi 2017), which could potentially change the outcome of dust coagulation as the collisional speeds are dependent on the sound speed (see Equation (9)).

Our computational domain covers a patch of the protoplanetary disk ranging from 4 to 34 au. While this domain allows us to cover most of the physics connected to dust drift, it is still a relatively small fraction of the global disk, which could extend to several hundreds au. One of the problems associated with not including the outer parts of the protoplanetary disk directly is that, without a proper boundary condition, we would quickly "run out" of dust. In the models presented in this paper, we adopted an outer boundary condition that allows an inflow of gas and dust, thus preventing the density at the outer edge from dropping significantly. However, particularly for the long runtime of the simulation, this condition cannot adequately account for an evolving pebble flux that is expected from global disk simulations (see, e.g., Birnstiel et al. 2012).

In the full-coagulation and 1D coagulation models, we adopted a relatively simple collision model with only two possible collision outcomes: sticking and fragmentation. We have assumed a single fragmentation threshold value in the whole domain (v_f = 10 m s⁻¹). While more complex collision models can be developed based on the results of laboratory experiments (see, e.g., Blum 2018), these are much harder to implement in the Smoluchowski equation solver (Windmark et al. 2012). We have also neglected the evolution of the porosity of dust aggregates, which can potentially lead to a different coagulation pattern (Ormel et al. 2007; Okuzumi et al. 2012; Krijt et al. 2016).

Despite these limitations, our results may have implications for the theory of pebble isolation mass. This concept assumes that delivery of solids to a growing gap-opening planets is halted if grains are large enough to be trapped in the pressure maximum outside of the planet orbit (see, e.g., Lambrechts et al. 2014; Ataiee et al. 2018; Bitsch et al. 2018). However, our results suggest that those large grains will fragment and constantly replenish the population of small grains, which are able to pass through the gap and potentially regrow in the planet co-orbital region (although the resolution of our models does not allow us to resolve the potential circumplanetary disk, in which the growth would be most efficient, see Shibaike et al. 2017; Drażkowska & Szulágyi 2018). Dust coagulation and fragmentation could thus increase the pebble isolation mass.

Similarly, our results cast doubt on the efficiency of dust filtration by growing Jupiter which is postulated to explain some features of the solar system. Efficient isolation of different reservoirs by a gap-opening planet, as postulated by, e.g., Kruijer et al. (2017), Alibert et al. (2018), Haugbølle et al. (2019), would only be possible if particles do not fragment during collisions (but at the same time are large enough to undergo efficient trapping). Because of the efficient fragmentation outside of the planet, small grains passing the gap around a growing Jupiter could still transport water into the inner solar system, in contrast to the idea proposed by Morbidelli et al. (2016), where the proto-Jupiter blocks the delivery of water when it reaches a mass of about 20 M_⊕.

The 2D hydrodynamic models include the effect of backreaction of dust on gas. We tested the importance of this effect by running a setup analogous to the full-coagulation run but with backreaction switched off. The comparison of these two runs is presented in Figure 9. The gas density is not modified significantly, but the effect of including backreaction is visible in the dust distribution. The dust ring formed outside of the planetary gap is placed a little bit farther away and it is slightly wider in the run including backreaction. This is because the large grains in the overdense region push on the gas, leading to a slight modification of the pressure gradient. The outward drift of dust in the pressure bump is sped up, leading to the wider ring profile. This effect was also observed by Kanagawa et al. (2018). Gonzalez et al. (2015) suggested that backreaction may lead to formation of a second pressure maximum, and as a consequence, a second dust ring caused by a single planet. We do not observe such an effect in our results, but this may be due to a difference in our setup: our planet is significantly less massive than the 5 M_J implemented by Gonzalez et al. (2015).

Zoom In Zoom Out Reset image size

The limited effect of backreaction we observe is a consequence of assuming a viscous disk with α = 10⁻³. Viscosity prevents the dust-to-gas ratio from becoming high: in the full-coagulation model, the maximum vertically integrated dust-to-gas ratio stays below 10% (see Figure 2). In disks with lower viscosity, planet–disk interactions lead to development of a vortex outside of the planet orbit (Li et al. 2009; Fu et al. 2014a; Hammer et al. 2017). The vortices are able to significantly concentrate dust and the effects of backreaction are more pronounced, including destruction of the vortex (Fu et al. 2014b; Crnkovic-Rubsamen et al. 2015; Raettig et al. 2015; Surville et al. 2016), although this effect is mitigated in 3D models (Lyra et al. 2018). We plan to study the effects of dust coagulation inside of a vortex in a future paper.