A Subgrid Model for the Growth of Dust Particles in Hydrodynamical Simulations of Protoplanetary Disks

In this work we employ the finite-volume code RoSSBi (Surville & Barge 2015; Surville & Mayer 2018), which uses the fluid approximation to treat both the gas and the dust component of the disk, solving the relevant equations in two dimensions and in cylindrical coordinates. The evolution of the dust and gas surface density (Σ_d, Σ_g) is described by the inviscid Euler equations in a cylindrical coordinate system. For the gas surface density one obtains

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{{\rm{\Sigma }}}_{g}+{v}_{g,r}\displaystyle \frac{\partial }{\partial r}{{\rm{\Sigma }}}_{g}+\displaystyle \frac{{v}_{g,\phi }}{r}\displaystyle \frac{\partial }{\partial \phi }{{\rm{\Sigma }}}_{g}+{{\rm{\Sigma }}}_{g}{\rm{\nabla }}\cdot {{\boldsymbol{v}}}_{g}=0,\end{eqnarray} \tag{ 1 }$$

with v_g,r, v_g,ϕ, and v_g denoting the radial, azimuthal, and total gas velocity, respectively. The evolution of the dust surface density reads

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{{\rm{\Sigma }}}_{d}+{v}_{d,r}\displaystyle \frac{\partial }{\partial r}{{\rm{\Sigma }}}_{d}+\displaystyle \frac{{v}_{d,\phi }}{r}\displaystyle \frac{\partial }{\partial \phi }{{\rm{\Sigma }}}_{d}+{{\rm{\Sigma }}}_{d}{\rm{\nabla }}\cdot {{\boldsymbol{v}}}_{d}=0,\end{eqnarray} \tag{ 2 }$$

with v_d,r, v_d,ϕ, and v_d denoting the radial, azimuthal, and total dust velocity, respectively.

In absence of additional force, the evolution of the gas velocity field follows

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}+({\boldsymbol{v}}\cdot {\rm{\nabla }}){\boldsymbol{v}}=-\displaystyle \frac{1}{{\rm{\Sigma }}}{\rm{\nabla }}P-{\rm{\nabla }}{\rm{\Phi }},\end{eqnarray} \tag{ 3 }$$

with Φ denoting the gravitational potential of the central star and P denoting the pressure. The interaction between gas and dust due to aerodynamical friction is implemented as a drag force

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{\mathrm{aero}}=-{{\rm{\Sigma }}}_{d}{{\rm{\Omega }}}_{K}(r){\mathrm{St}}^{-1}({{\boldsymbol{v}}}_{d}-{{\boldsymbol{v}}}_{g}),\end{eqnarray} \tag{ 4 }$$

with ${{\rm{\Omega }}}_{K}^{2}={{GM}}_{\star }/{r}^{3}$ describing the Keplerian orbital frequency and St being the Stokes number (see Weidenschilling 1977). The Stokes number is the ratio of the stopping time t_stop and the turnover time of the largest turbulent eddy t_turnover. The latter may be written as ${t}_{\mathrm{turnover}}={{\rm{\Omega }}}_{{\rm{K}}}^{-1}$ (see Cuzzi et al. 2001), and assuming the Epstein drag law, an isothermal volumetric gas density profile with the gas density in the midplane ${\rho }_{g}\,={{\rm{\Sigma }}}_{g}{{\rm{\Omega }}}_{k}/\sqrt{2\pi }{c}_{s}$ and spherical particles, we can write the Stokes number as

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

with a denoting the particle radius and ρ_s the particle density. In this work we always use 1 g cm⁻³ for the particle density ρ_s.

Combining Equations (3) and (4) in polar coordinates leads to the radial gas velocity evolution,

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{v}_{g,r}+{v}_{g,r}\displaystyle \frac{\partial }{\partial r}{v}_{g,r}-\displaystyle \frac{{v}_{g,\phi }^{2}}{g,r}+\displaystyle \frac{{v}_{g,\phi }}{r}\displaystyle \frac{\partial }{\partial \phi }{v}_{g,r}\\ &&\quad =\,-\displaystyle \frac{1}{{{\rm{\Sigma }}}_{g}}\displaystyle \frac{\partial }{\partial r}P-r{{\rm{\Omega }}}_{K}^{2}-\displaystyle \frac{{{\boldsymbol{f}}}_{\mathrm{aero}}\cdot {{\boldsymbol{e}}}_{r}}{{{\rm{\Sigma }}}_{g}},\end{eqnarray} \tag{ 6 }$$

and for the azimuthal gas velocity component we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{v}_{g,\phi }+\displaystyle \frac{1}{r}{v}_{g,\phi }\displaystyle \frac{\partial }{\partial \phi }{v}_{g,\phi }+\displaystyle \frac{{v}_{g,r}}{r}{v}_{g,\phi }+{v}_{g,r}\displaystyle \frac{\partial }{\partial r}{v}_{g,\phi }\\ &&\quad =\,-\displaystyle \frac{1}{{{\rm{\Sigma }}}_{g}}\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial \phi }P-\displaystyle \frac{{{\boldsymbol{f}}}_{\mathrm{aero}}\cdot {{\boldsymbol{e}}}_{\phi }}{{{\rm{\Sigma }}}_{g}},\end{eqnarray} \tag{ 7 }$$

with e_r and e_ϕ being the radial and azimuthal unit vectors, respectively. The corresponding equations for the dust can be written as

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{v}_{d,r}+{v}_{d,r}\displaystyle \frac{\partial }{\partial r}{v}_{d,r}-\displaystyle \frac{{v}_{d,\phi }^{2}}{r}+\displaystyle \frac{{v}_{d,\phi }}{r}\displaystyle \frac{\partial }{\partial \phi }{v}_{d,r}\\ &&\quad =\,-r{{\rm{\Omega }}}_{K}^{2}+\displaystyle \frac{{{\boldsymbol{f}}}_{\mathrm{aero}}\cdot {{\boldsymbol{e}}}_{r}}{{{\rm{\Sigma }}}_{d}}\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{v}_{d,\phi }+\displaystyle \frac{1}{r}{v}_{d,\phi }\displaystyle \frac{\partial }{\partial \phi }{v}_{d,\phi }+\displaystyle \frac{{v}_{d,r}}{r}{v}_{d,\phi }+{v}_{d,r}\displaystyle \frac{\partial }{\partial r}{v}_{\phi }\\ &&\quad =\,\displaystyle \frac{{{\boldsymbol{f}}}_{\mathrm{aero}}\cdot {{\boldsymbol{e}}}_{\phi }}{{{\rm{\Sigma }}}_{d}}.\end{eqnarray} \tag{ 9 }$$

The evolution of pressure is obtained by solving the adiabatic energy equation. If the energy conservation is solved, the corresponding pressure can be calculated after each time step from the total energy of the adiabatic gas:

$$\begin{eqnarray}&&E=\displaystyle \frac{P}{1-\gamma }+\displaystyle \frac{1}{2}{{\rm{\Sigma }}}_{g}{{\boldsymbol{v}}}_{g}^{2},\end{eqnarray} \tag{ 10 }$$

with γ = 1.4. Hence, the pressure is updated according to the energy conservation

$$\begin{eqnarray}&&P=(1-\gamma )\cdot \left(E-\displaystyle \frac{1}{2}{{\rm{\Sigma }}}_{g}{{\boldsymbol{v}}}_{g}^{2}\right).\end{eqnarray} \tag{ 11 }$$

We use the fluid approximation for the gas and a pressureless fluid model for the dust content. This assumption is justified because in our model the dust particles stay small, within the Epstein drag regime. Consequently, particles within one grid cell have nearly the same properties and can therefore be computationally modeled as a fluid. However, particles in the Epstein regime mainly interact with the gas molecules and not with other dust particles; hence, the pressure-less fluid assumption has to be used.

The conservative form of this system of coupled equations is solved by the RoSSBi code using a well-balanced finite-volume method; see Surville & Barge (2015) and Surville et al. (2016). The time integration is based on a second-order Runge–Kutta scheme, and the flux reconstruction is third order in space using parabolic interpolation. The numerical fluxes of the gas fluid are obtained by an exact Riemann solver, and the ones of the dust fluid are obtained by a pressureless Roe solver used in Paardekooper & Mellema (2006).

The boundary conditions implemented in the code RoSSBi are based on zero gradient conditions, where ghost cell variables are reconstructed to follow the steady-state profiles of the disk. These free conditions account for radial flux of gas and dust. However, in order to keep stability of the dust fluid (in particular at the outer disk boundary), only the dust density is damped toward the initial profile.

Additionally, the gravity of an embedded planet is implemented in the code RoSSBi. In some of the runs presented in this study, a Jupiter-mass planet is orbiting on a circular orbit r_p around the star. As a simplification, the center of mass of the star/planet system is kept at the star center, which is the origin of the reference frame.

The field of gravitational force exerted by the planet of mass M_p on the disk is given by

$$\begin{eqnarray}&&{{\boldsymbol{g}}}_{\mathrm{pla}}({\boldsymbol{r}})=-\displaystyle \frac{{{GM}}_{p}}{{[| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{p}| +l({r}_{p})]}^{3}}({\boldsymbol{r}}-{{\boldsymbol{r}}}_{p}),\end{eqnarray} \tag{ 12 }$$

where the gravity is modified using a well-known softening length $l({r}_{p})=0.6{H}_{0}({r}_{p})$ for a planet potential. The isothermal disk scale height at the planet orbit is ${H}_{0}({r}_{p})\,={[{P}_{0}({r}_{p})/{{\rm{\Sigma }}}_{0}({r}_{p})]}^{1/2}/{{\rm{\Omega }}}_{K}({r}_{p})={c}_{s}({r}_{0})/{{\rm{\Omega }}}_{r}({p}_{0})$, with c_s denoting the sound speed. The actual profiles of the gas background pressure and density will be given in Section 2.3. Finally, the mass of the planet is loaded directly from the beginning of the run, which has no critical influence on the evolution of the disk later on. The drag interaction between dust and gas fluids is solved using an implicit method described in Surville & Mayer (2018). Finally, simulations including a planet are done using a thermal relaxation term in the energy equation to avoid shock heating of the disk. This additional source term is computed implicitly. Thus, the time step is determined by the CFL condition, with a factor of 0.5 needed by the parabolic reconstruction of the RoSSBi scheme. For these simulations, except in the ZHU and DZHU simulations (see Section 3.4), we use an average grid size of 0.0244 au in radius and an average aspect ratio of 0.285. The usage of a cell aspect ratio (CAR) of ∼0.25–0.3 has been tested and robustly confirmed in several papers using the RoSSBi code (Surville & Barge 2015; Surville et al. 2016; Surville & Mayer 2018) for the evolution of vortices, and also the convergence at higher resolutions for these CAR has been tested during the preparation of the paper. Therefore, we can resolve the Hill radius, in the simulations containing a planet, with 78 grid cells in radius and 69 cells in azimuth.

We base our approach for modeling dust sizes on the so-called two-population algorithm proposed by Birnstiel et al. (2012). The idea behind that model was to reproduce the general pattern of dust surface density evolution obtained with the dust coagulation code of Birnstiel et al. (2010) at much reduced computational cost, without actually solving the dust coagulation equation. In each grid cell the dust size a_max is chosen to represent the full dust size distribution, based on a semianalytical approach. Making the right choice of a_max is possible thanks to the comprehensive understanding of the processes governing dust evolution presented by Birnstiel et al. (2012).

The original dust evolution model of Birnstiel et al. (2012) was 1D. Here we extend this method to two dimensions, by applying the subgrid dust growth and fragmentation prescription and keeping the original advection scheme of the RoSSBi code. In this model, the representative dust particle size, determined with the help of local gas and dust variables (see Equations (14) and (17)), is calculated before each Runge–Kutta time step and in each single grid cell. This representative grain size is then used to compute the aerodynamical friction source terms in the RoSSBi code. As the required calculation is local to each grid cell, this subgrid method does not change the parallelization model of the hydrodynamical method.

The representative size a_max is found in each cell by comparing the maximum aggregate size that could be obtained taking into account various physical processes: dust growth (Equation (17)), fragmentation (Equations (14) and (15)), and the loss of large aggregates due to radial drift (Equation (16)). We pick the smallest of these sizes as the representative size characterizing the local dust population. In what follows we describe how we model such individual processes.

The growth timescale can be written as

$$\begin{eqnarray}&&{\tau }_{\mathrm{grow}}=\displaystyle \frac{a}{\dot{a}}\approx \displaystyle \frac{1}{Z{{\rm{\Omega }}}_{K}},\end{eqnarray} \tag{ 13 }$$

with Z standing for the vertically integrated dust-to-gas ratio. In the inner parts of protoplanetary disks, this timescale is typically faster than the global dust redistribution timescale (Birnstiel et al. 2012; Dra̧żkowska et al. 2016). In such a case, the dust size distribution is governed by a coagulation–fragmentation equilibrium. Since the impact speeds of dust aggregates increase with size, there is a maximum size that can be obtained by dust growth before it is halted by fragmentation. If the dominant source of the impact speeds is turbulence, the representative size is

$$\begin{eqnarray}&&{a}_{\mathrm{frag}}=\displaystyle \frac{{f}_{{\rm{f}}}}{3}\displaystyle \frac{2{{\rm{\Sigma }}}_{g}}{\pi {\rho }_{s}}\displaystyle \frac{{u}_{f}^{2}}{\alpha {c}_{s}^{2}},\end{eqnarray} \tag{ 14 }$$

where f_f is an order of unity constant, α is the turbulence parameter (see Shakura & Sunyaev 1973), ρ_s is the internal density of dust particles, and u_f denotes the threshold fragmentation speed, which we set to u_f = 10 m s⁻¹. The turbulence that drives impact speeds has typical eddy overturn timescales comparable to a typical stopping time of dust grains; these correspond, in our models, to eddies of roughly the size of one grid cell. This implies that we could not resolve them and therefore the α parameter is describing a subgrid turbulence. In this paper we chose α = 10⁻³ because it is the standard value used in dust coagulation models. The maximum particle size possible to obtain with respect to the impact speeds triggered by the differential drift can be written as

$$\begin{eqnarray}&&{a}_{\mathrm{df}}=\displaystyle \frac{2{{\rm{\Sigma }}}_{g}}{\pi {\rho }_{s}}\displaystyle \frac{{u}_{f}{V}_{k}}{{c}_{s}^{2}(1-N)}{\left|\displaystyle \frac{d\mathrm{ln}P}{d\mathrm{ln}r}\right|}^{-1},\end{eqnarray} \tag{ 15 }$$

where N is defined as a typical ratio between the Stokes numbers of two colliding particles and V_k is the Keplerian velocity. Following Birnstiel et al. (2012), we use N = 0.5, since this gives the best fit to complete models.

In the outer parts of the protoplanetary disk, particles' growth timescale becomes longer than the radial drift timescale. Since in our approach the advection of dust does not have a direct effect on dust size, this process cannot be explicitly modeled. Therefore, we need to include a drift limit, which takes this effect explicitly into account. Leaving such a limit out of consideration would lead to an overprediction of the grain size because if the drift timescale is shorter than the growth timescale, dust grains should be removed by the radial drift before they grow to the size limited by fragmentation. The maximum size that can be kept at a given orbital distance before it would be removed by the drift can be written as

$$\begin{eqnarray}&&{a}_{\mathrm{drift}}={f}_{{\rm{d}}}\displaystyle \frac{2{{\rm{\Sigma }}}_{g}}{\pi {\rho }_{s}}\displaystyle \frac{{V}_{k}^{2}}{{c}_{s}^{2}}{\left|\displaystyle \frac{d\mathrm{ln}P}{d\mathrm{ln}r}\right|}^{-1},\end{eqnarray} \tag{ 16 }$$

with f_d being the numerical constant for drift. Moreover, we use an initial growth limit

$$\begin{eqnarray}&&{a}_{\mathrm{ini}}={a}_{0}\cdot \exp \left(\displaystyle \frac{t}{{\tau }_{\mathrm{grow}}}\right),\end{eqnarray} \tag{ 17 }$$

with a₀ denoting the initial particle size and τ_grow denoting the growth timescale (described by Equation (13)). This limit takes into account that the growth timescale significantly increases with the orbital distance, so the large particles occur in the outer part of the disk much later than in its inner part (Birnstiel et al. (2012) called this effect "the delayed drift effect," and Lambrechts & Johansen (2014) called it "the pebble formation edge").

The algorithm of Birnstiel et al. (2012) considered two characteristic sizes of dust particles in each cell: the minimum and the maximum size. In our code, we are restricted to a single size per grid cell. Hence, the smallest dust size in each cell is neglected, and the whole dust surface density is assumed to be generated by the maximum-sized particles. This is a good estimate if the particle size is limited by radial drift, because 97% of the dust surface density is determined by the mass of the a_max-sized particles. As for the turbulence barrier, 75% of the dust mass is in the maximum-sized particles.

In this work we chose the minimum-mass solar nebula (MMSN) model (see Hayashi 1981) as the initial dust surface density model, but the configuration of the RoSSBi code allows us to change it to another model in the future. This particular choice is justified with the MMSN being a well-studied model for which it is relatively simple to find simulations to be compared with this work. Within the MMSN model the initial gas surface density can be written as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{g}(r\gt 2.7\,\mathrm{au})=1700{\left(\displaystyle \frac{r}{\mathrm{au}}\right)}^{-3/2}\left[\displaystyle \frac{{\rm{g}}}{{\mathrm{cm}}^{2}}\right].\end{eqnarray} \tag{ 18 }$$

The assumption of an adiabatic ideal gas combined with a power-law surface density profile leads to the following simple form of the pressure profile of the disk (see also Surville et al. 2016):

$$\begin{eqnarray}&&P={P}_{0}\cdot {\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{-2},\end{eqnarray} \tag{ 19 }$$

with r₀ the reference radius, which corresponds to the planet location in the runs with a planet (r₀ = 10 au). In order to normalize the pressure and set the temperature profile, we normalize the disk scale height to H₀(r₀) = 0.05r₀. Supplementary information can be found in Surville et al. (2016). The domain size ranges always from 5 to 30 au, and simulations are stopped after 400 orbits at 10 au. For the simulations with an embedded planet we consider always a Jupiter-mass planet at 10 au. We note that in runs containing a planet we apply a fast cooling function to the energy equation, which makes the simulations effectively isothermal. The fast cooling function can be written as

$$\begin{eqnarray}&&{f}_{\mathrm{cool}}=-{{\rm{\Sigma }}}_{g}(T-{T}_{0})\displaystyle \frac{{{\rm{\Omega }}}_{k}}{{\tau }_{\mathrm{cool}}},\end{eqnarray} \tag{ 20 }$$

with τ_cool denoting the cooling constant, T the temperature, and T₀ the reference temperature at the planet position. In this study we choose τ_cool = 10⁻⁴.

Many hydrodynamic grid codes rely on the assumption that dust size (or, alternatively, the Stokes number) is fixed across the simulated domain. Often this size is set to represent large grains, which are expected to grow in some regions of protoplanetary disks. The top panel in Figure 1 presents the results of such a simulation, with the dust size fixed to 3 cm (RN32). The gas phase exhibits no significant evolution past the initial condition of this run. But the dust component is axisymmetrically piling up in the inner part of the disk. This is a result of the inward drift and the inflow at the outer boundary. Since the particle size is fixed, the Stokes number of particles increases with radial distance. Additionally, the radial drift velocity increases with the Stokes number (until St < 1); hence, particles in the outer part of the domain drift faster than those in the inner part, causing the increase of dust density at the inner edge of the domain.

Zoom In Zoom Out Reset image size

The bottom panels of Figure 1 present results of the corresponding run including our subgrid model for dust evolution (DN32). For consistency, we initialized this simulation with 3 cm grains, which means that the initial growth stage is bypassed, as dust size is immediately adjusted to the maximum possible so that at each radial distance fragmentation and radial drift become the governing mechanisms. The panels clearly illustrate that the pattern of surface density evolution is quite different from that in the fixed-size simulation. Particularly, no pileup of dust in the inner disk arises. The dark annulus marks depletion of dust that develops at the boundary between the regions dominated by fragmentation, which is triggered by turbulence and radial drift.

The first conclusion to be drawn is that there is a qualitative change in the dust evolution pattern depending on whether we assume fixed or variable dust size even if we start from large grains in both cases. We can now proceed to compare the dust evolution pattern when we start from small grains and consider the initial growth stage. Figure 2 presents results of two runs (RN1-42 and DN1-42) that begin with dust grains having a size of 1 μm, which is typically assumed as an initial size in works modeling dust growth in protoplanetary disks (Windmark et al. 2012). In the top panels of Figure 2, this size is fixed throughout the evolution, and thus the dust surface density practically does not evolve, as the small grains stay well coupled to the gas. The bottom panels of Figure 2 display the inside-out evolution pattern commonly obtained in 1D dust growth models (Birnstiel et al. 2012; Krijt et al. 2016). The particles grow faster at smaller orbital distances, where the Keplerian frequency is higher (see Equation (17)). Those particles then decouple form the gas disk and drift inward, before the particles in the outer part of the domain manage to reach their maximum size.

Zoom In Zoom Out Reset image size

Figure 3 illustrates the results of two runs including a Jupiter-mass planet at 10 au: one with the dust size fixed to 3 cm (RY32; top panels), and one including our dust evolution model starting from grains having a size of 3 cm (DY32; bottom panels). Adding a planet to the simulations affects not only the dust component but also the gas. This is explained by the additional gravitational potential of the planet, as well as by the back-reaction of the dust fluid to the gas component. Our results for the gas surface densities and gap opening are comparable to previous studies (e.g., Rosotti et al. 2016). As seen in the top panels of Figure 3, the planet-induced pressure bump stops the migration of particles, and therefore a ring structure behind the planet is formed (see, e.g., Paardekooper & Mellema 2004, 2006; de Val-Borro et al. 2007; Fouchet et al. 2010; Gonzalez et al. 2012; Zhu et al. 2012, 2014). Particles in the inner part of the planetary orbit can simply drift toward the central star and vanish from the simulated disk. Since in the fixed-size RY32 simulation the 3 cm particles have high drift velocities, particles from the outer part of the disk arrive at the pressure bump quickly.

Zoom In Zoom Out Reset image size

The bottom panels of Figure 3 show results of the run with a planet and the variable dust size algorithm (DY32). The evolution of the disk is less dramatic and hence slower than the evolution of the fixed-size simulation. This is because our dust size routine returns sizes that are, with the exception of the very inner part of the domain, significantly lower than the 3 cm used in the fixed-size case. The presence of the planet triggers strong pressure gradients, and thus particles in the disk are now also controlled by the radial drift. In the variable size simulation we do not obtain the pronounced dust annulus at the outer edge of planetary gap as observed in the fixed-size simulation. Also, the depletion of the inner disk is significantly weaker.

The first part of this section exemplified the results of each individual method. In order to compare them in a more qualitative approach, in this subsection we investigate histograms of azimuthally averaged and normalized values for the dust and gas surface density component.

In the absence of the planet, the azimuthal averaging is straightforward and does not reduce the amount of information. However, using one analysis tool for all simulations and applying it to axisymmetric solutions may be used as a test for its accuracy. Because the simulations start with the same conditions, a change in the results is a product of the different particle sizes.

The simulations with a 3 cm fixed dust size (RN32), a 10⁻⁴ cm fixed dust size (RN1-42), and a variable dust size, starting with 3 cm sized particles (DN32), without a planet show no significant change of the gas surface density profiles, as presented in Figure 4. This is expected since a change of the gas surface density occurs only via the back-reaction from dust to gas. Because of the low dust-to-gas ratio, the impact of the back-reaction during the simulated timescales in an axisymmetric simulation should be negligible, and therefore the changes in the gas surface density are insignificant.

Zoom In Zoom Out Reset image size

On the other hand, the dust surface density evolution is directly affected by the particle size owing to radial drift. As shown, in Figure 4, the large dust sizes lead to a pileup in the inner part of the disk. In contrast, the dust surface density for the small particles is nearly unchanged, because their drift velocity is much lower and the evolution time is too short to allow for a significant drift.

The azimuthally averaged outcomes of the RY32, RY1-42, and DY32 runs including the planet, after 400 orbits, can be seen in Figure 5. Only a slight difference among the three gas surface densities is visible. Thus, the impact of differently sized particles appears to be negligible for the gas evolution, since the overall shape and position of the curve are the same for all three simulations. The dust surface density, however, exhibits a completely different picture. In the case of the large fixed-size particles, the disk gets almost depleted, since mainly dust particles, which are trapped in the pressure bump behind the planet, remain in the disk. Particles from the outer edge of the disk migrate to this barrier, and the surface density peaks at this position within the disk. On the other hand, the small particles are well coupled to the gas, and as a result, they migrate slowly and are not trapped by the pressure bump. Comparing the two fixed-size simulations RY32 and RY1-42 (Figure 5) yields an enormous difference in the final dust surface density.

Zoom In Zoom Out Reset image size

Our main focus, however, is the comparison with the new subgrid model. In this respect the results presented so far suggest that, when a fixed-size ought to be used, small sizes (e.g., less than 10⁻¹ to 10⁻² cm) are a preferable choice relative to large sizes (e.g., more than 1 cm) because with the former one achieves a closer match with results obtained when the dust evolution model is employed. This arises from the fact that particles may fragment (or drift), which leads naturally to the production of smaller particles. There is only a small fraction of the inner disk that allows particles to grow up to 3 cm; hence, simulations with small sizes are more realistic than those using large sizes. Nevertheless, a more careful inspection leads to the conclusion that the simulation with small fixed-size particles (RY1-42) also differs appreciably from the simulation with dust size evolution. In the case of an evolving dust size, particles grow until they halt at the pressure-induced barrier, at which point the depletion of the inner part of the disk starts. Instead, in the corresponding fixed-size simulation (RY1-42) this does not happen. A possible objection to such inference could be that this particular comparison is made between runs employing large particles (with an initial size of 3 cm), which is unlikely to be a realistic starting point. To address the issue further, the same comparison was done with a 3 cm fixed dust size simulation RY32, a 10⁻⁴ cm fixed dust size simulation RY1-42, and a variable dust size simulation starting from 10⁻⁴ cm dust sized particles DY1-42; see Figure 6. In this case dust particles start growing from a realistic size and do not directly start with fragmentation. This, however, seems to favor the choice of a small grain size even more, most likely because of the same starting size of the particles.

Zoom In Zoom Out Reset image size

In any case, this result should be regarded with caution. Indeed, we compared the three simulations only after a relatively short timescale, 400 orbits, which corresponds to only a small fraction of the typical disk evolution timescale, which is of orders of Myr. Therefore, we expect that particles would grow further after 400 orbits, thus deviating progressively more from the simulation adopting small fixed-size grains. In order to investigate this further, we employ two additional simulations, RY1-12 and RY1-22, and compare again with the DY1-42 simulation. These simulations do not include dust evolution but employ particles with an intermediate size (10⁻² and 10⁻¹ cm, respectively); see Figure 7. Clearly these two new simulations match the dust evolution run even better, simply because at the time chosen for the comparison growth has proceeded to a stage such that the maximum dust size is closer to the size of grains in these two simulations. While this suggests that the match with a given fixed-size simulation will depend on the time at which the comparison is carried out, it is also arguable that, in fixed-size simulations, the choice of an intermediate dust size will in general yield a better match to the results of dust evolution simulations.

Zoom In Zoom Out Reset image size

In this section we compare the outcome of our dust evolution model against previous results obtained with the ATHENA code (Stone et al. 2008) using fixed dust size, presented by Zhu et al. (2014, hereafter ZH14). We choose to compare to this particular work because the ATHENA code also uses a Godunov-type finite-volume hydro method in polar coordinates, similar to RoSSBi. Dust, however, is treated with Lagrangian particles advected through the grid rather than with a second fluid as in our case. Another important difference is that there is no back-reaction in the ZH14 simulations. We first use the RoSSBi code with a fixed dust size to reproduce the ZH14 results, and afterward we rerun this with the dust evolution model on this particular simulation setup. We choose to reproduce the results of the simulation "M02D2" in ZH14. Hence, we construct a nearly identical disk setup with an embedded planet having a mass of 8 Earth masses (using the same thermal mass definition to initialize it as ZH14; see caption of Figure 8) and evolve it using an isothermal equation of state (see runs ATH and D-ATH in Tables 1 and 2). While the goal is to have the same resolution as ZH14 in our simulation, we use a logarithmically spaced grid in the radial coordinate, while ZH14 uses a linearly spaced grid. Therefore, we choose the radial resolution in order to account for this difference, constructing a grid with resolution 282 × 1024 (see Tables 1 and 2). The outcome of our fixed-size simulation aimed to reproduce the ZH14 results can be seen in Figure 8. The dust density profiles are similar in shape and magnitude in comparison to the original results. We argue that the agreement between our results and the ZH14 results is sufficient to investigate the impact of our dust evolution scheme. Residual differences are likely caused by the different grid spacing, which prevents us from comparing at truly identical resolution (we tested indeed that changing grid resolution has an effect on both the gas and dust surface density profiles, especially near the planet). The impact of the dust evolution method is displayed in Figure 8 as the model 282 (D) after 200 orbits. It shows a significant deviation from the fixed-size simulations, in line with our findings in Sections 3.1–3.3. Therefore, the comparison with ZH14 reinforces our general inference that including dust evolution has a major effect on the dust distribution in the disk, changing the spatial distribution and amplitude of overdense regions.

Zoom In Zoom Out Reset image size

Our numerical method is limited by several factors, which we briefly discuss here for clarity.

The original dust coagulation method was developed in the framework of an 1D model, whereas in this paper we use the method for 2D modeling. Therefore, the coefficients of the semianalytical formulae (see Equations (14) and (16)) might potentially change owing to this transition, but at the moment there is no possibility of running a 2D hydrodynamical simulation with the full dust coagulation calculation, using a grid code, that we could use to calibrate the method. Hence, we adhere to the value of the coefficients proposed by Birnstiel et al. (2012).

Another simplification that we already pointed out in Section 2.2 is that only one representative size of the dust fluid is used in each cell, not two sizes as proposed in the original algorithm (Birnstiel et al. 2012). This could change the resulting dust surface density up to a few percent. Nonetheless, this should not impact our conclusions, which are based on the relative difference between including and not including dust evolution.

The main challenge of our dust size treatment is that the dust size is not advected. The advection velocity depends on the local size, and thus the density evolution is impacted by dust growth and fragmentation, but we calculate the representative size in each grid cell and in each time step based on the local conditions at a given time only. The only exception is the initial growth stage, when particles grow gradually before they reach the fragmentation or drift barrier, and then their size depends on the evolution time. After this is completed, the particle size in a given cell changes based on local conditions such as the dust-to-gas ratio and sound speed (see Equations (14) and (16)), but it does not depend on the size that existed in the cell in a previous time step or on the sizes obtained in the neighboring cells. This is equivalent to assuming that dust coagulation does always have enough time to produce an equilibrium size distribution. This is true in the case of a homogeneous disk, but it may not be true when fast dynamical changes occur in the disk, such as the ones introduced by the planet. Using this model might lead to a discontinuous growth in particle size after each time step and hence affect the calculated surface densities. This is an obvious drawback that we aim to fix in our future work.

Additionally, our current code does not use any kind of sink particles, and consequently the particles are trapped behind the planet but not accreted. Introducing an accreting planet could potentially cause a decrease of the magnitude of overdensity formed behind the planet, but on the other hand, the increasing planetary mass would strengthen the pressure bump that is causing the overdensity.

The limitations we mentioned above contribute to the uncertainty on the dust density evolution we obtain in our models. However, we argue that even our approximate dust growth prescription yields a significant improvement over simulations with fixed dust size because it shows that the difference is both quantitative and qualitative, which potentially has important implications for the interpretation of observations. A similar result has been found for an SPH-based code; see Gonzalez et al. (2012, 2015a, 2015b).

By construction, our simulations are inviscid, and therefore no large-scale turbulence is explicitly modeled. However, we assume that there is some level of subgrid turbulence that triggers impact speeds between dust particles that limit their growth (see Equation (14)).

Since RoSSBi solves the hydrodynamical equations using the Godunov method, it also has intrinsically a very low numerical viscosity. The gas is in fact evolved in a very similar way to that with ATHENA, concerning the unsplit advection and second-order time accuracy. However, the RoSSBi code uses an exact Riemann solver rather than an HLLC-type (Stone et al. 2008), as well as a well-balanced scheme. In most published work codes either are inherently more viscous, such as FARGO (Masset 2000), which is a finite-difference scheme, or explicitly apply some form of viscosity to capture shocks, such as, again, in FARGO, but also in ZEUS (Stone & Norman 1992a, 1992b; Stone et al. 1992; Clarke 1996; Hayes et al. 2006), or in the SPH simulations of Gonzalez et al. (2015a, 2017). Viscosity would smooth out gaps and other features triggered by the presence of an embedded planets, generally reducing the sharpness of pressure bumps, therefore affecting the response of dust particles. This will have to be investigated in the future, but, once again, here we are focusing less on absolute effects and more on the relative differences between having and not having dust evolution accounted for.

On the other hand, our method includes the back-reaction from dust to gas, which some of the other models neglect. This effect becomes particularly important when the dust-to-gas ratio becomes high, as the gas disk structure is modified (Gonzalez et al. 2017; Kanagawa et al. 2017). The back-reaction should also lead to the well-known phenomenon of the streaming instability, which leads to a spontaneous clumping and formation of dusty filaments, which can subsequently become gravitationally unstable and collapse to planetesimals (Johansen et al. 2007; Kowalik et al. 2013; Simon et al. 2016). The models we present in this paper do not have enough resolution to resolve the streaming instability. Nevertheless, even a high-resolution 2D simulation might not capture such a phenomenon, since the streaming instability is mainly an effect of the xz-plane. However, the numerical method we developed may in fact help us in future work to revisit the streaming instability under more realistic conditions, as all of the previous models relied on treating dust with fixed size, despite that it is well known that the streaming instability is sensitive to dust size (Youdin & Johansen 2007; Bai & Stone 2010). Likewise, alternative models in which overdense dust clumps are produced by vortex-drag instabilities (Crnkovic-Rubsamen et al. 2015; Surville et al. 2016; Surville & Mayer 2018), which would later undergo streaming instability or collapse owing to their own self-gravity, will be affected by dust evolution, although the general character (but not the timescale) of such vortex-induced instabilities appears to be independent of dust size (Surville & Mayer 2018).

The main purpose of this work was to compare the dust surface density evolution with dust coagulation scheme to the fixed-size simulations. Because particle growth up to centimeter sizes is reasonably well understood and confirmed by laboratory experiments (e.g., Güttler et al. 2010), most of the previous simulations assumed a large grain size. As logical as this may sound, it leads, however, to a completely different evolution of the disk. The main problem encountered within this method is the fact that a protoplanetary disk tends to contain the largest particles in the inner part; thus, the drift timescales for particles in the outer disk are slower than in simulations with a fixed size. Consequently, while in simulations with fixed dust size a pileup in the inner part of the disk is not surprising, simulations with dust growth do not exhibit such an extremely overdense region. Including a planet in the simulation makes the difference even more pronounced. Our simple dust evolution model allows us to simulate multiple particle sizes within the dust fluid, hence also allowing small particles, which can then migrate through the planetary-induced pressure bump. In contrast to that in the fixed-size simulation, particles are held at the pressure bump in the disk. On the other hand, simulating a disk with small fixed-size dust particles also leads to problems. In this case the planet cannot stop particles from migrating; as a consequence, this approach yields different results than the dust growth simulation.

As already mentioned before, it is not sufficient to choose the maximum or minimum fixed size of dust particles to reproduce the evolution of a protoplanetary disk. According to our results, dust growth is important and cannot be neglected. Since the state-of-the-art grid code simulations do not allow particle growth, it is reasonable to choose an intermediate dust particle size, which yields outcomes that are closer to the evolution that includes dust growth. The results are shown in Figure 6, which compares the DY1-42 simulation, including dust growth prescription, to the fixed particle size simulations with 10⁻¹ and 10⁻² cm particles. In these last comparisons, the results of the fixed dust size simulations resemble more closely those of the dust evolution simulations relative to the results of the RY32 and R1-42 fixed-size simulations. Clearly, the choice of a fixed size in order to reproduce the results of the dust evolution algorithm will depend on the duration of the simulation. If a simulation starts with small growing particles and lasts only for a few thousand years, a fixed particle size that is only slightly higher than the starting size may be chosen. A very long simulation may require slightly larger particles, since the subgrid method would allow particles to grow to larger sizes as the simulation continues. Also, the inclusion of a planet affects the choice of the particle size. If a large planet is chosen, then particles should be smaller than those selected in the case of a small planet. This is justified by the fact that a large planet will induce a larger pressure bump, which will abruptly stop larger particles. On the contrary, small particles can pass through the induced pressure bump for a longer time.

A major trend of our results is that, in general, when our dust evolution model is included, sharp features in the dust component arising in disks with embedded massive planets tend to weaken significantly (see also Zhu et al. 2012). Such features, such as dust rings produced by dust pileups near gap edges, are a recurrent feature of 1D simulations that assume an axisymmetric background (e.g., Pinilla et al. 2012). We argue that these previous results have largely exaggerated the strength of dust rings and other features resulting from the perturbation of a massive planet and the formation of a gap. As a result, we caution about using observed dust rings in disks such as HL Tau to infer the presence of a planet, as the outcome is largely dependent on the role of dust evolution (Flock et al. 2015; Ragusa et al. 2017). Here we have just begun to address this important issue with a simple-minded subgrid model. The implications are of paramount importance for interpretation of upcoming surveys of T Tauri and debris disks carried out in various wavelengths with instruments such as ALMA and SPHERE (e.g., Gonzalez et al. 2012). As our simulations also show a distribution of dust sizes across the disk even after 400 orbits, synthetic observations will have to be carried out in multiple wavelengths with accurate Monte Carlo radiative transfer codes in post-processing in order to provide a useful test bed to interpret the upcoming observations, following on the footsteps of analogous work carried out using the RADMC3D tool for ALMA mocks of dust continuum emission from 3D spiral density waves and clumps in self-gravitating disks (Dipierro et al. 2015a; Mayer et al. 2016; Meru et al. 2017), as well as from 3D circumplanetary disks (Szulágyi et al. 2018). Indeed, while it is clear that with our dust evolution model sharp features and transitions are smoothed out, it is also expected that this and other effects depend on the specific dust size range considered and thus will be wavelength dependent.