Ring Morphology with Dust Coagulation in Protoplanetary Disks

Tidal interactions between the embedded planets and their surrounding protoplanetary disks are often postulated to produce the observed complex dust substructures, including rings, gaps, and asymmetries. In this Letter, we explore the consequences of dust coagulation on the dust dynamics and ring morphology. Coagulation of dust grains leads to dust size growth that, under typical disk conditions, produces faster radial drifts, potentially threatening the dust ring formation. Utilizing 2D hydrodynamical simulations of protoplanetary disks that include a full treatment of dust coagulation, we find that if the planet does not open a gap quickly enough, the formation of an inner ring is impeded due to dust coagulation and subsequent radial drift. Furthermore, we find that a “buildup” of submillimeter-sized grains often appears in the dust emission at the outer edge of the dust disk.


Introduction
The high-resolution observations of protoplanetary disks (PPDs) by the Atacama Large Millimeter/submillimeter Array (ALMA) have revealed some of the most stunning structures yet (e.g., Andrews et al. 2018;Huang et al. 2018a;Isella et al. 2018;Macías et al. 2019). Many disks exhibit a system of rings and gaps in the dust emission. Several mechanisms have been proposed for creating rings of dust in PPDs; one of the most commonly accepted scenarios is that embedded protoplanets interact tidally with the disk, opening up a gap and creating pressure bumps that effectively trap dust into rings. Recently, Dong et al. (2017Dong et al. ( , 2018 found that a single super-Earth-sized planet may open up several rings and gaps within its own orbit. However, in order to trap dust within the planet's orbit, the planet must compete with dust loss from radial drift in the inner disk. In a typical protoplanetary disk, the gas drag on dust grains tends to cause dust to drift toward the star (Whipple 1972;Weidenschilling 1977;Brauer et al. 2008;Birnstiel et al. 2012). This can be formulated in terms of the Stokes number of the grain , where a is the dust grain size, r p its internal density, and S g the gas surface density. The radial drift timescale is given by where P is the vertically integrated gas pressure, » S c g s 2 , and | | g = d ln P d ln r . For typical disks, using g = 1.41, = c v s K 0.05, and = St 1 (which gives the largest drift velocity), the drift timescale is about 90 orbits.
Due to the computational expense as well as physics uncertainties, dust coagulation is often neglected and singlespecies or multi-species approaches are used with prespecified, fixed distributions (Fu et al. 2014;Miranda et al. 2016;Dong et al. 2017;Ricci et al. 2018;Zhang et al. 2018). However, as dust grains grow in size, their Stokes numbers increase according to St and drift timescale shortens according to Equation (1). Simulations that neglect coagulation may therefore underestimate the effects of radial drift, especially if the initial dust size is small. Unless there is a feature blocking the radial drift, large dust grains are typically lost as they travel inward and accrete onto the star. Hence, including dust coagulation may change the outcome of ring formation dramatically.
Turbulent velocities are postulated to be the main source of collisions for dust grains, and these velocities also grow with the Stokes number. Once turbulent velocities exceed the fragmentation velocities of the grains, large grain growth will halt as the particles fragment into smaller pieces. For Stokes numbers 1, we have the following maximum size to which dust particles grow, beyond which the dust size distribution falls exponentially Pinilla et al. 2012;Li et al. 2019): where a vis is the viscosity parameter for coagulation, v f is the fragmentation velocity for the dust, and c s is the sound speed in the gas. For our studies here, we use v f =10 m s −1 . This implies thatã 8 cm max , if we use S = -20 g cm . This fragmentation limit on particle size is strict, if the initial dust size is much smaller than a max .
There is much uncertainty about the coagulation timescale in PPDs (Brauer et al. 2007), as the typical micron-sized dust grains from the ISM gradually grow in size inside the PPD. We adopt the approach of Birnstiel et al. (2012), which derives the particle growth timescale as approximately , where  0 is the initial dust-to-gas ratio in the disk (Ormel & Cuzzi 2007;Youdin & Lithwick 2007;Brauer et al. 2008). The coagulation timescale to the maximum fragmentation size a max , then, is given by which is about 180 orbits if we use =  0.01 0 , = a 1 0 μm, and a max from Equation (2).
By considering when t t » drift coag , Birnstiel et al. (2012) derives the upper limit on dust grain size due to losses from radial drift to be where S d is the dust surface density, and the prefactor is tuned by simulations. Putting in the typical values, we getã 5.5 cm drift , which is slightly smaller than a max .
When a planet is embedded in a PPD, it tries to open a gap in the disk. In order to estimate the gap-opening timescale, we adopt the equation from Lin & Papaloizou (1986), which neglects the effects of viscous stress: In this equation, H 0 is the aspect ratio of the disk at the planet's radius and M p is the mass of the planet. The viscous timescale is typically much longer than the gap-opening timescale, especially for the nearly inviscid disks in this work, and so for the purposes of our study it is safe to ignore viscous stress. orbits, similar to t gap . Hence, we expect a significant dust fraction will coagulate and drift inward from the planet's radius before it opens a gap. (Note that the combined effect of coagulation and drift is more complicated than adding these two timescales.) In the inner regions of a PPD, the dust coagulation and subsequent loss due to radial drift competes with the gapforming process of an embedded planet. Due to coagulation, dust grains in the inner disk grow to their fragmentation-limited size a max over a time t coag . After reaching size~a max , the large grains drift inward over a time t drift if no pressure trap has formed. Hence, if t t t > + gap coag drift , we do not expect a ring of dust to form inside the planet's orbit, and vice versa.
In this Letter, we explore how the gap-opening process (which is determined by M M p * and H 0 ) affects the ring formation whenever full dust coagulation is included in the disk's evolution. In Section 2, we will discuss our simulation methods and setups. In Section 3, we present our main results, followed by a summary and discussion in Section 4.

Methods
We perform 2D hydrodynamic simulations coupled with a full-coagulation treatment of the dust size evolution in a model PPD. We then post-process synthetic dust continuum images to analyze the observational signatures of dust coagulation.

Hydrodynamic Model
Hydrodynamic simulations of the PPDs were performed with LA-COMPASS code (Li et al. 2005(Li et al. , 2009. We choose an isothermal sound speed profile of where H is the scale height of the gas and r p is the planet's radius, which we set to 30au. At this radius, 1000 orbits corresponds to approximately 150,000 yr. This profile corresponds to a locally isothermal temperature profile that scales as µ -T r 1 2 . We note that Miranda & Rafikov (2019) recently showed that an adiabatic equation of state with exponent g = 1.001 can have an effect on the gap-opening process for planets. For all of our simulations, the star has a mass of 1M e , and the viscosity parameter a =´-5 10 5 using the α prescription. Our simulation domain extends from 12 to 480 au, with a grid of resolution of´=f n n 1024 1024 r . We examined the numerical convergence of our simulations using higher-resolution (´2048,´3072) single-species runs. We found the inner ring formation outcome to be the same for the single-species runs below, while the structure and peak density of the outer ring could change with the resolution. The disk is initialized with an axisymmetric gas profile of The dust is set with an initial dust-to-gas ratio of 1%. These parameters lead to a total disk mass of around 38M J . The planet masses vary between 10 Å M and 50 Å M , and H 0 varies between 0.03 and 0.07. We choose M p and H 0 combinations that lie in the vicinity of the critical region where t t t » + gap coag drift . We grow the planet over a period of 10 orbits at the beginning of the simulation.

Coagulation Model
The gas and dust dynamics equations are as in Fu et al. (2014), and we include the dust feedback reaction on the gas. The coagulation scheme is described in Li et al. (2019) and Drążkowska et al. (2019), using the model of Birnstiel et al. (2010). In order to simulate dust dynamics and coagulation, we use 151 species of dust logarithmically spaced between 1 μm up to 1 m sized boulders, which gives us 25 species per size decade. The dust distribution is initialized at a size of 1 mm. LA-COMPASS includes a fluid for each species of dust, and to determine the evolution of the dust sizes we explicitly integrate the Smoluchowski equation in each spatial cell (Smoluchowski 1916). Turbulence and radial drift are considered as sources of relative velocities between the dust particles. The turbulence for dust is governed by a separate viscosity parameter a vis , for which we choose a value of 10 −3 . We note that a a ¹ vis in our runs. This is because a vis represents the dust turbulent velocities at the midplane, which may differ from the global viscous evolution of the disk characterized by α.  (Brauer et al. 2008;Birnstiel et al. 2010).
We handle the dust coagulation in each cell using an operator splitting approach alongside the gas and dust hydrodynamics. To save computing time, we implement a substepping routine and calculate the coagulation outcomes every 50 time steps.

Radiative Transfer
To create synthetic dust emission images of our disk models, we utilize the radiative transfer code RADMC-3D of Dullemond et al. (2012). The detailed procedure is the same as described in Huang et al. (2018b) and Li et al. (2019). For RADMC-3D, we reduce the grid size in the r and f dimensions to´q n n r =´f n 300 40 100. The algorithm used to downsample and extrapolate the dust density recovered over 94% of the total dust mass for all cases. In order to calculate the dust temperature for emission, we utilized the thermal photon Monte Carlo capabilities of RADMC-3D with 500 million photons. The grain-sizedependent dust opacity is adopted from Isella et al. (2009) and Ricci et al. (2010). The star is assumed to have a blackbody temperature of 4200 K.
For all of the disk models, we set the observation distance to 100 pc and simulated the dust emission at zero inclination at 1000 orbits. The star is assumed to have a blackbody temperature of 4200K and a mass of 1M  . We imaged the disk models at 890 (ALMA Band 7), 1330 (ALMA Band 6), and 3000 μm (ALMA Band 3). The emission map is convolved with a Gaussian beam with FWHM 35 mas in order to simulate observation of the disk, which is comparable to the resolution achieved with the recent DSHARP survey .

Results
In order to illustrate the effects of planet mass and aspect ratio on the formation of rings interior to the planet's orbit, we will study three representative models in more detail with a combination of parameters = M 24, 40 p Å M and = H 0.04, 0 0.05. Additionally, for each of them we ran a separate singlespecies simulation in order to compare. We will refer to them as M24H005C/S, M24H004C/S, and M40H005C/S, respectively. (The letters C and S refer to coagulation and single species, respectively.) We use the density-weighted average dust size in the inner ring from the coagulation runs as the characteristic sizes for single-species runs. The dust size is 1.34, 1.26, 2.23 cm for M24H005S, M40H005S, and M24H004S, respectively. Many other runs are also made to map out the parameter space. First, we will use the model M24H005C to verify the analytical formulas for the characteristic scales given in the introduction.

Characteristic Scales
In Figure 1, we have the radial dust size distributions taken at 250 orbits and 1000 orbits (the last frame of our simulation). Along with the size distributions, we have plotted a max and a drift . At 250 orbits, the dust in the inner regions of the disk is approaching the fragmentation-limited size a max , while dust in the outer regions is still submillimeter. At 1000 orbits, the particles have reached a max in the rings, and they fit the shape of a drift well in the gaps. The edge at a drift is less sharp than a max since the radial drift growth limit is not as strict as the fragmentation limit. This illustrates that in our simulations the coagulation, fragmentation, and radial drift are all behaving close to expectation.

Full-coagulation Runs
In Figure 2, we have plotted the dust density for each of the models with full coagulation. In M24H005C, t t t > + gap coag drift , and we see that the innermost ring is very weak. However, increasing the planet mass or decreasing H 0 lowers t gap , and we see that both M40H005C and M24H004C form strong inner rings. In M24H005C, an inner pressure maximum forms in the gas, but the dust coagulates and drifts beforehand.
The density-averaged size of the dust within the outer ring at 37au of M24H005C is given in the bottom right panel of Figure 2. For comparison, we have also plotted the average size for a region just outside the ring at 40au. Inside the ring, we see an enhancement of the dust size due to the fact that the dust density is higher in this region. Size enhancement within forming rings boosts radial drift and contributes to the diminished inner ring in M24H005C. There are also variations in the average size along the f-direction, which we expect to be time-dependent at this stage of disk evolution.

Single-species Runs
The dust profiles for the single-species runs are presented in Figure 3. For all three runs, the ring within the planet's radius either fails to form completely or is greatly diminished. This is due to the fact that the dust is large from the start of the simulation and can immediately drift quickly. For these reasons, the singlespecies dust emission plots in Figure 3 exhibit a much darker inner disk than M24H005C. One ring forms outside of the planet, and only a portion of the "horseshoe" dust region at the planet's radius is present. We also performed these single-species runs with a small dust size of 0.02 cm, and in each case three rings formed. Because of the small size and slow drift, the rings were wider and the inner disk was brighter.
For comparison, M24H005C clearly exhibits four distinct rings. Three immediately surround the planet's orbit, and a fourth faint ring is located outside of these planetary gap rings. The faint, outermost ring is from a buildup of submillimetersized dust particles that have grown from the μm-sized grains in the outer disk but have not yet reached larger sizes where they drift to the planet very quickly. As dust grains grow past 1 mm near this ring, they drift inward and are quickly lost. This creates the gap inside the "buildup" ring, since radial drift speed is nonlinear. This feature is not present in M40H005S due to the fact that the dust size is uniformly large near ∼1 cm.

Varying Aspect Ratio and Planet Mass
Using the same method of comparing peak dust density values between the inner and outer rings as above, we have compiled the ratios for all of our simulations after 1000 orbits into Figure 4.
We find good agreement between analytical predictions and our simulation results. We used a max as the characteristic size for the critical region where t t t » + gap coag drift . Inside and close to this line, we see a gradient between the two regimes, where the rings partially exist but are weaker than the region in the bottom right of the parameter space. If the planet gap opens in approximately the same amount of time as the dust coagulates and drifts inward, we still expect some dust to be trapped in the inner ring.

Synthetic Dust Emission
The dust emission maps for our coagulation runs have been plotted in Figure 5, for wavelengths of 890, 1330, and 3000 μm. In all three wavelengths, the innermost ring of M24H005C is much fainter than that of M24H004C or Decreasing the disk aspect ratio similarly creates a much stronger inner ring. (d) Density-averaged dust size for M24H005C in the outer ring ( = r 37 au) and outside the outer ring ( = r 40 au). We see dust size enhancement inside the ring where the dust density is higher along with variation of the average density in the f-direction.

M40H005C.
The "buildup" outermost dust ring is present in all three wavelengths, but it is relatively strongest in 890 μm. The emission maps in Figure 5 all exhibit asymmetries in the outer ring adjacent to the planet, due to the asymmetrical dust profile in this region. These asymmetries are most apparent in the longer wavelengths, corresponding to larger dust particles that are less coupled to the gas. In the overdense regions of the outer ring, elevated dust densities could precipitate more efficient particle growth.

Discussion and Summary
In this Letter, we have conducted hydrodynamic and dust coagulation simulations of PPDs spanning planet masses from 10 to 50 Å M and disk aspect ratios between 0.03 and 0.07 in order to test our analytical predictions for the formation of a dust ring interior to the planet's orbit. The results of our simulations matched our predictions well, and we found that dust coagulation and subsequent radial drift impedes the inner dust ring from forming fully for high aspect ratios and low planet masses, when t t t > + gap coag drift . The dust distributions in simulations that included coagulation differed from equivalent single-species runs in several major ways: First, the inner dust ring was saved from radial drift for runs with high planet masses and low aspect ratios. Second, there was an additional "buildup" ring outside of the rings immediately surrounding the planet gap. Third, the ring exterior to the planet's orbit tends to be asymmetric.
One possible explanation of the asymmetrical outer dust ring could be that the planet is triggering the Rossby wave instability (Lovelace et al. 1999;Li et al. 2000Li et al. , 2001. The asymmetries in the dust emission were strongest in the longest wavelength. Vortices would be populated by larger particles because dust grains with high Stokes numbers move more quickly toward the pressure maxima and the elevated gas densities there allow particles to grow to larger sizes. The varying dust density around the outer ring changes the size distribution of the dust, which affects the spectral index and opacity of the region. These effects are exacerbated if the dust grows to larger sizes and concentrates into vortices more effectively. There are several physical uncertainties in our coagulation model, such as the fragmentation velocity and collision outcomes. The size a max is particularly sensitive to the fragmentation velocity, with µ a v max frag 2 from Equation (2). However, changes to the fragmentation velocity would shift the boundary line for the parameter space in Figure 4 but maintains its general behavior. Additionally, our assumption that the initial dust size is uniform at 1 μm throughout the disk is simplistic. In the future, more sophisticated initial dust size distributions can be studied along with different planets masses and disk parameters. Given these uncertainties, our study is meant to highlight the interplay between dust growth, radial drift, and ring formation in an approximate manner.
Our simulations neglected disk self-gravity (DSG), which can play an important role in the gap-opening process of a disk. In particular, for some initial gas surface density values,  including DSG may boost the gap-opening process and decrease t gap (Zhang et al. 2014). This would effectively move the critical band in Figure 4 upward. In preliminary LA-COMPASS simulations that include DSG, we do observe this effect, as t gap decreases and a stronger inner ring forms for a run equivalent to M24H005C. This Letter provides simulation evidence that dust coagulation can impact the dust ring formation, brightness, and asymmetries. They also provide valuable inferences about the planet mass, disk aspect ratio, and coagulation dynamics in a PPD system. These results motivate additional detailed studies on the physics of dust coagulation in a more realistic PPD environment.