Effects of Ringed Structures and Dust Size Growth on Millimeter Observations of Protoplanetary Disks

TL;DR

This study uses hydrodynamical simulations to show that ringed structures and dust growth significantly influence millimeter observations of protoplanetary disks, explaining low spectral indices and brightness variations.

Contribution

It demonstrates how ringed structures and dust size evolution jointly affect millimeter flux and spectral indices, providing insights into disk observations and planet formation processes.

Findings

Ringed structures slow dust radial drift and promote growth.

High fragmentation velocity disks can have mm spectral index ~2.0.

Multiple rings can produce disks with spectral indices below 2.5.

Abstract

The growth of solids from sub-micron to millimeter and centimeter sizes is the early step toward the formation of planets inside protoplanetary disks (PPDs). However, such processes and their potential impact on the later stages of solid growth are still poorly understood. In this work, we test the hypothesis that most disks contain at least one ringed structure with a relative small radius. We have carried out a large family of 1D two-fluid (gas+dust) hydrodynamical simulations by evolving the gas and dust motion self-consistently while allowing dust size to evolve via coagulation and fragmentation. We investigate the joint effects of ringed structures and dust size growth on the overall sub-millimeter and millimeter (mm) flux and spectral index of PPDs. Ringed structures slow down the dust radial drift and speed up the dust growth. In particular, we find that those unresolved disks with a high fragmentation velocity ($\sim10\ {\rm m~s^{-1}}$) and a high dust surface density ($\sim10\ {\rm g\ cm^{-2}}$ in the ring) can have mm spectral indices as low as $\sim2.0$, consistent with mm observations of faint disks in nearby star forming regions. Furthermore, disks with more than one ringed structure can potentially reproduce brighter disks with spectral indices lower than $\sim2.5$. Future multi-wavelength high-resolution observations of these low spectral index sources can be used to test the existence of the ringed structures in the unresolved disks and differentiate the effects of dust size growth from optical depth.

Figures

Tables

Table 1. Table 1: Model parameters and simulation results

model	${\Sigma }_0$	$v_f$	${\alpha }_0$	$h_0$	$\gamma$	$r_{\mathrm{ring}}$	$w_{\mathrm{ring}}$	$d_{\mathrm{ring}}$	$F_{1mm}$	${\alpha }_{0.8-1.3\mathrm{mm}}$	${\alpha }_{1.0-3.0\mathrm{mm}}$	${\alpha }_{3.0-7.0\mathrm{mm}}$
	$\left(g{\mathrm{cm}}^{-2}\right)$	$\left(\mathrm{cm}{s}^{-1}\right)$				$\left(\mathrm{au}\right)$	$\left(\mathrm{au}\right)$		$\left(\mathrm{mJy}\right)$
m0v1	12.8	$3\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	4.79	2.92	3.01	2.41
m1v1	4.3	$3\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	1.52	3.00	3.00	2.30
m2v1	0.43	$3\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	0.11	2.67	2.62	2.12
m3v1	0.043	$3\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	0.005	2.41	2.36	2.08
m0v2	12.8	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	51.5	2.38	2.44	2.42
m1v2 (fiducial)	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	32.5	2.43	2.47	2.33
m2v2	0.43	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	15.7	2.24	2.35	4.10
m3v2	0.043	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	3.12	3.86	3.77	3.53
m0v3	12.8	$3\times {10}^2$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	69.7	2.69	2.67	2.71
m1v3	4.3	$3\times {10}^2$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	44.2	2.52	2.55	3.98
m2v3	0.43	$3\times {10}^2$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	9.32	3.52	3.53	3.51
m3v3	0.043	$3\times {10}^2$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	1.43	3.55	3.57	3.50
m0v4	12.8	$1\times {10}^2$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	69.9	2.80	2.91	3.50
m1v4	4.3	$1\times {10}^2$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	34.4	3.27	3.33	3.49
m2v4	0.43	$1\times {10}^2$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	7.41	3.41	3.45	3.50
m3v4	0.043	$1\times {10}^2$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	5.0	0.90	3.56	3.58	3.48
…	…	…	…	…	…	…	…	…	…	…	…	…
m1a1	4.3	$1\times {10}^3$	$5\times {10}^{-3}$	0.04	1.0	15.0	5.0	5.0	23.1	2.30	2.34	2.36
m1a2	4.3	$1\times {10}^3$	$2\times {10}^{-3}$	0.04	1.0	15.0	5.0	5.0	4.75	3.05	3.16	2.85
m1a3	4.3	$1\times {10}^3$	$5\times {10}^{-4}$	0.04	1.0	15.0	5.0	5.0	6.76	2.99	3.11	3.10
…	…	…	…	…	…	…	…	…	…	…	…	…
m1d1	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	2.0	29.8	2.82	2.97	3.40
m1d2	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	10.0	23.7	2.28	2.31	2.31
m1d3	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	5.0	20.0	16.6	2.41	2.48	2.80
…	…	…	…	…	…	…	…	…	…	…	…	…
m1r1	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	25.0	5.0	5.0	41.7	2.61	2.64	2.54
m1r2	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	50.0	5.0	5.0	52.8	2.95	3.01	3.28
m1r3	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	75.0	5.0	5.0	39.4	3.33	3.45	3.93
…	…	…	…	…	…	…	…	…	…	…	…	…
m1w1	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	2.5	5.0	27.5	2.84	2.89	2.33
m1w2	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	7.5	5.0	34.8	2.43	2.50	2.36
m1w3	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	15.0	10.0	5.0	32.5	2.43	2.47	2.32
…	…	…	…	…	…	…	…	…	…	…	…	…
m1h1	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.03	1.0	15.0	5.0	5.0	25.9	2.31	2.37	2.39
m1h2	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.06	1.0	15.0	5.0	5.0	35.9	2.58	2.50	2.65
…	…	…	…	…	…	…	…	…	…	…	…	…
m1b1	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	0.5	15.0	5.0	5.0	26.4	2.10	2.12	2.37
m1b2	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.5	15.0	5.0	5.0	45.2	2.28	2.35	2.52
…	…	…	…	…	…	…	…	…	…	…	…	…
nogap	4.3	$1\times {10}^3$	$1\times {10}^{-2}$	0.04	1.0	$-$	$-$	$-$	5.14	4.04	4.29	3.61

Equations

$${\Sigma }_{\mathrm{g}}(r)={\Sigma }_0{\left(\frac{r}{r_{\mathrm{c}}}\right)}^{-\gamma }\exp \left[-{\left(\frac{r}{r_{\mathrm{c}}}\right)}^{2-\gamma }\right],$$

$$\frac{c_{\mathrm{s}}}{v_{\mathrm{K}}}=h_0{\left(\frac{r}{r_0}\right)}^{0.25},$$

$$\frac{\partial {\Sigma }_{\mathrm{g}}}{\partial t}+\nabla \cdot ({\Sigma }_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}})=0,$$

$$\frac{\partial ({\Sigma }_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}})}{\partial t}+\nabla ({\mathbf{v}}_{\mathrm{g}}\cdot {\Sigma }_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}})+\nabla P=-{\Sigma }_{\mathrm{g}}\nabla {\Phi }_{\star }+{\Sigma }_{\mathrm{g}}{\mathbf{f}}_{\nu }-{\Sigma }_{\mathrm{d}}{\mathbf{f}}_{\mathrm{d}},$$

$${\mathbf{f}}_{\mathrm{d}}^i=\frac{{\Omega }_{\mathrm{k}}}{{\mathrm{St}}^i}({\mathbf{v}}_{\mathrm{g}}-{\mathbf{v}}_{\mathrm{d}}^i),$$

$${\mathrm{St}}^i=\frac{\pi {\rho }_{\mathrm{s}}{a}^i}{2{\Sigma }_{\mathrm{g}}}.$$

$$a_{\mathrm{St}=1}=\frac{2{\Sigma }_{\mathrm{g}}}{\pi {\rho }_{\mathrm{s}}}.$$

$$\frac{\partial {\Sigma }_{\mathrm{d}}^i}{\partial t}+\nabla \cdot ({\Sigma }_{\mathrm{d}}^i{\mathbf{v}}_{\mathrm{d}}^i)=\nabla \cdot \left({\Sigma }_{\mathrm{g}}{D}_{\mathrm{d}}^i\nabla \left(\frac{{\Sigma }_{\mathrm{d}}^i}{{\Sigma }_{\mathrm{g}}}\right)\right).$$

$$\frac{\partial {\Sigma }_{\mathrm{d}}^i{\mathbf{v}}_{\mathrm{d}}^i}{\partial t}+\nabla ({\mathbf{v}}_{\mathrm{d}}^i\cdot {\Sigma }_{\mathrm{d}}^i{\mathbf{v}}_{\mathrm{d}}^i)=-{\Sigma }_{\mathrm{d}}^i\nabla {\Phi }_{\star }+{\Sigma }_{\mathrm{d}}^i{\mathbf{f}}_{\mathrm{d}}^i.$$

$$v_{r,\mathrm{d}}=\frac{{\mathrm{St}}^{-1}{v}_{r,\mathrm{g}}-\eta v_{\mathrm{K}}}{{\mathrm{St}}^{-1}+\mathrm{St}},$$

$$v_{r,\mathrm{g}}\simeq -\frac{3}{{\Sigma }_{\mathrm{g}}\sqrt{r}}\frac{\partial }{\partial r}({\Sigma }_{\mathrm{g}}{\nu }_{\mathrm{g}}\sqrt{r}),$$

$$a_{\max }=\frac{4{\Sigma }_{\mathrm{g}}}{3\pi {\alpha }_{\mathrm{vis}}{\rho }_{\mathrm{s}}}\frac{v_{\mathrm{f}}^2}{c_{\mathrm{s}}^2},$$

$$I(r)=B_{\nu }\left(T_{\mathrm{d}}(r)\right)\left(1-e^{-{\tau }_{\nu }(r)}\right),$$

Full text

accretion, accretion disks — protoplanetary disks — planets and satellites: rings — submillimeter: planetary systems

††journal: ApJ††software: Astropy (Astropy Collaboration et al., 2013), LA-COMPASS (Li et al., 2005, 2009), Matplotlib (Hunter, 2007), Numpy (van der Walt et al., 2011), RADMC-3D (Dullemond et al., 2012)

1 Introduction

The discovery of thousands of exoplanets over the last couple of decades has clearly shown that the birth of planets is a very efficient process in nature (e.g., Burke et al., 2015). The commonly accepted theory is that planets form in young disks orbiting pre-main sequence (PMS) stars through the agglomeration of small $\mu$m dust particles into km-sized “planetesimals”, which are massive enough to gravitationally attract other solids in the disk (see review by Chiang & Youdin, 2010). Grain growth from sub-micrometer sizes, which are the typical sizes of dust in the interstellar medium (ISM), to millimeter and centimeter sizes is thus the first step toward the formation of planetesimals inside young circumstellar disks. However, the formation of planetesimals is still poorly understood due to i) the low sticking efficiency of grains larger than $\sim 1$ mm-cm, and ii) the mutual dynamical interaction between dust grains and the gas in the disk. In the case of a gas-rich disk with density and temperature both decreasing with distance from the central star, small solids radially drift inward because of the aerodynamic drag by the gas orbiting at sub-Keplerian speeds (Weidenschilling, 1977; Brauer et al., 2007). Models of the evolution of disk solids have predicted radial drift timescales which are too short ($100\sim 1000$ orbits) to form planetesimals (see review by Johansen et al., 2014). Understanding dust coagulation might be of fundamental importance to understand the formation of planetesimals.

Once dust particles grow from ISM-like sizes of $\mu$m to mm/cm sizes, they can be mainly traced by their continuum emission at sub-millimeter (sub-mm) and longer wavelengths. Therefore, sub-mm observations have the potential of setting strong constraints on the initial stages of the planet formation process. In the last two decades, several observations of protoplanetary disks (PPDs) at sub-mm/mm wavelengths have reported relatively shallow slopes with values of the mm spectral index $2.0\leq\alpha_{\rm 1-3mm}\leq 3.0$ ($F_{\nu}\propto\nu^{\alpha}$). These values are significantly lower than those measured for dust in the ISM, and are taken as observational evidence for dust growth to mm-sized grains in disks under the assumption that the dust emission is mostly optically thin (e.g., Testi et al. 2003; Natta et al. 2004; Rodmann et al. 2006; Ricci et al. 2010a; Pinilla et al. 2014; Andrews 2015; Ansdell et al. 2018). Ricci et al. (2012) explored the effect of local optically thick emission regions in the mm emission without mm/cm grains. They found that optically thick regions with relatively small filling factors can reproduce the mm spectral indices of disks without grain growth, although the dust over-densities required by the optically thick effect are generally larger than those predicated by the typical known physical processes proposed in the literature. These processes include vortex formation, streaming instabilities, disk viscosity transitions (e.g., Lovelace et al., 1999; Li et al., 2000, 2001; Klahr & Bodenheimer, 2003; Rice et al., 2004; Fromang & Nelson, 2005; Johansen & Youdin, 2007; Johansen et al., 2009; Boss, 2010; Regály et al., 2012).

Birnstiel et al. (2010b) have tested the grain growth model by millimeter observations of dust continuum, and they found that it can naturally reproduce the observed mm-slope if the radial drift of the dust is halted by an unknown mechanism. It is generally believed that some inhomogeneities in the gas density, either azimuthally symmetric, such as rings, or with azimuthal asymmetries, such as vortices, have to be invoked to slow down the dust radial drift because of the gas pressure structure (e.g., Pinilla et al., 2012, in the case of dust trapping by azimuthally symmetric rings). The dust can, therefore, be trapped effectively. Grains can collide frequently and stick together by van der Waals forces, which leads to the formation of large-size dust particles (e.g., Dominik & Tielens, 1997; Poppe et al., 2000; Blum & Wurm, 2008; Brauer et al., 2008; Birnstiel et al., 2010a). Such dust accumulation regions can appear as bright rings in PPDs observed from the dust continuum emission at millimeter (mm) wavelengths. There is now strong observational evidence for the existence of multiple rings (or rings sandwiched by gaps) in disks by high-resolution observations with the Atacama Large Millimeter Array (ALMA), e.g., HL Tau (ALMA Partnership et al., 2015), TW Hya (Andrews et al., 2016; Tsukagoshi et al., 2016), HD 163296 (Isella et al., 2016), AA Tau (Loomis et al., 2017), Elias 24 (Cieza et al., 2017; Cox et al., 2017; Dipierro et al., 2018), AS 209 (Fedele et al., 2018), GY 91 (Sheehan & Eisner, 2018), V1094 Sco (Ansdell et al., 2018; van Terwisga et al., 2018), MWC 758 (Boehler et al., 2018), more disks in a recent survey of young disks in Taurus (Long et al., 2018), 20 disks from the Disk Substructures at High Angular Resolution Project (DSHARP) observed very recently (Andrews et al., 2018; Dullemond et al., 2018; Huang et al., 2018b), and 16 disks from ALMA archival data (van der Marel et al., 2019). These high-resolution observations have also revealed that these narrow rings have relatively high optical depth at mm band (e.g., ALMA Partnership et al., 2015; Tsukagoshi et al., 2016; Andrews et al., 2018; Dullemond et al., 2018; Isella et al., 2016, 2018; Huang et al., 2018a, b). Dullemond et al. (2018) find that the narrow radial size of the rings strongly supports the dust trapping scenario. Such dust traps can also facilitate the grain growth.

Motivated by the fact that dust rings appears to be ubiquitous at least among the disk sample observed so far, we explore their consequences on the overall (sub-)mm disk emission properties. For the sake of simplicity, we use a non-monotonic disk viscosity profile (e.g., Figure 1) to produce a local gas surface density bump. We perform detailed two-fluids (gas+dust) hydrodynamics simulations, in which the dust coagulation and disk hydrodynamics evolution are treated together, and calculate the dust continuum emission at several mm bands with detailed radiative transfer. Our main goal is to interpret the low spectral indices observed for many young disks, and understand the link between the spectral slope and the dust dynamics in the disks. Based on the comparison between our simulations and the observed dust emission, we can constrain some fundamental parameters for the dust coagulation model and provide new insight into the mechanism of dust trapping and growth.

The paper is organized as follows. The numerical simulation method adopted in this work is presented in Section 2. We present the results in Section 3, and summarize the main results and discuss the physical implications in Section 4.

2 Model

We envision the following scenario: the disk gas surface density profile contains at least one ring in the inner disk region (which is unresolved for most disks currently). By performing 1D hydrodynamical simulations using LA-COMPASS (Li et al., 2005, 2009; Fu et al., 2014), we study the gas + dust dynamics spanning several million years, including dust coagulation and fragmentation self-consistently (Birnstiel et al., 2010a), and then produce the dust emission at different mm wavelengths by making use of the RADMC-3D package (Dullemond et al., 2012) to compare with observations.

2.1 Hydrodynamics and Model Assumption

Without loss of generality, we consider a disk around a PMS star with a mass of $0.6~{}M_{\sun}$ at a distance of 140 pc. The disk extends from 5 au to 240 au. Disk self-gravity in gas and dust is not included. To create dust rings in such disks, we use varying disk viscosity radial profiles to achieve such configurations. The gas viscosity is adopted from the $\alpha$-prescription $\nu_{\rm g}=\alpha_{\rm vis}c_{\rm s}h_{\rm g}$ with $h_{\rm g}$ being the gas disk height and $\alpha_{\rm vis}$ the dimensionless viscosity parameter (Shakura & Sunyaev, 1973). For simplicity, in our fiducial model, $\alpha_{\rm vis}(r)=\alpha_{0}=10^{-2}$ for the most region of the disk. We also assume a rectangular gap with a low viscosity $\alpha_{\rm vis}$ in its radial profile. The viscosity gap is centered at 15 au with a gap width $w_{\rm ring}=5\ {\rm au}$ and a value for the viscosity in the dip $\alpha_{\rm min}=\alpha_{0}/(1+d_{\rm ring})=1.7\times 10^{-3}$, which is shown in Figure 1. As we will see, the viscosity gap can produce a gas bump, which can then slow down the radial drift of dust particles in the disk, creating a dust ring.

The existence of a dust ring in the disk, though a hypothesis in this paper, is well motivated by current high-resolution observations mentioned above. Actually, any mechanism that can produce a similar gas bump, even without invoking a viscosity jump, could produce results similar to the ones discussed here. Both viscosity transition (e.g., Varnière & Tagger, 2006) and the presence of a planet (e.g., Bryden et al., 1999; Papaloizou et al., 2007; Dipierro et al., 2015; Jin et al., 2016; Dong et al., 2017; Liu et al., 2018; Zhang et al., 2018) could lead such a dust ring configuration. The viscosity transition can be characterized by a gap in the viscosity around a few au (e.g., see review by Armitage 2011), which is slightly smaller than the radius of our viscosity gap. However, the inner edge and the radial width of the dead zone are actually not well constrained and will depend in general on the ionization sources and the physical structure of the disk, as well as on the intensity of the magnetic field in the disk itself. Given these uncertainties, it makes sense to explore a broad range of values for the central radius and width of the viscosity jump, as will be done in this work. For the planet-disk interaction scenario, some kinematic evidence of protoplanets forming within ringed-structure disks have been found recently based on the detection of a localized deviation from Keplerian velocity in HD 163296 (Pinte et al., 2018; Teague et al., 2018), although other possibilities (e.g., zonal flow) for the inner ringed-structures cannot be fully ruled out yet.

The initial gas surface density profile $\Sigma_{\rm g}(r)$ is set as follows:

[TABLE]

where $r_{\rm c}=60~{}\rm au$ and $\gamma=1$ unless otherwise stated. These parameters represent a typical disk as seen from observations (Andrews et al., 2010). The normalization of gas surface density $\Sigma_{0}$ is treated as a parameter listed in Table 1. The locally isothermal sound speed $c_{\rm s}$ is chosen as

[TABLE]

where $v_{\rm K}(r)$ is the local Keplerian velocity, $r_{0}=5$ au unless otherwise stated. $h_{0}$ is listed in Table 1 with a typical value of $0.04$. This corresponds to a disk temperature profile as $T\propto r^{-0.5}$. Equation (2) also expresses the radial profile of the gas scale height $h_{\rm g}/r$. The gas and dust fluids are evolved following the conservation of mass, radial, and angular momentum equations.

The continuity equations of gas is

[TABLE]

where $\mathbf{v}_{\rm g}$ is the gas fluid velocity. The momentum equation for the gas is

[TABLE]

where $\bf{f_{\nu}}$ indicates the viscous force from the Shakura-Sunyaev disk (Shakura & Sunyaev, 1973). We adopt an equation of state $P=c_{\rm s}^{2}\Sigma_{\rm g}$ for the gas component, where $P$ is the vertically integrated gas pressure. $\Phi_{\star}$ is the gravitational potential of the central star, $\mathbf{f}_{\rm d}$ is the summation of $\mathbf{f}_{\rm d}^{i}$ over $i$ and $\mathbf{f}_{\rm d}^{i}$ is the drag force between the gas and dust species $i$, which is defined as

[TABLE]

where $\Omega_{\rm k}$ is the Keplerian angular velocity. ${\rm St^{i}}$ is the Stokes number of the dust species $i$. We can estimate this feedback term based on Takeuchi & Lin (2002). In the Epstein regime for the disk and dust parameters of interest here, the Stokes number of the particles with dust radius $a$ in the mid-plane of the disk is defined as

[TABLE]

The dust size corresponding to ${\rm St}=1$ is thus

[TABLE]

The Stokes number in the Epstein regime can then be re-written as ${\rm St}^{i}=a^{i}/a_{\rm St=1}$.

We treat the dust component as a pressure-less fluid. The dust feedback, i.e., drag forces between the gas and dust, are incorporated into the momentum equation for both the gas and dust (Fu et al., 2014). The continuity equation of dust species $i$ is:

[TABLE]

Here $D_{\rm d}$, which describes the dust diffusivity, is related with the Stokes number ${\rm St}$ as $D_{\rm d}=\nu_{\rm g}/(1+{\rm St}^{2})$ (Youdin & Lithwick, 2007), $\nu_{\rm g}$ is the gas viscosity. $\Sigma_{\rm d}^{i}$, $\mathbf{v}_{\rm d}^{i}$ are dust surface density and velocity for species $i$, respectively. The momentum equation for dust species $i$ after including the drag force from the gas is

[TABLE]

For an illustrative purpose, the radial velocity of the dust can be estimated as

[TABLE]

where $\eta=-\frac{c_{\rm s}^{2}}{v_{\rm K}^{2}}\frac{d\ln P}{d\ln r}$, and

[TABLE]

is the radial velocity of the gas. Here we omit the superscript of $i$ for different dust species. We can see that the dust radial velocity is composed of two terms, namely gas drag (the first term) and radial drift (the second one) in Equation (10). Then the drag force between dust and gas for each dust species can be obtained by combining Equations (5,10$-$11).

To model the dust size growth during disk evolution, we calculate the evolution of dust size distribution which is discretized logaritmically between $1\ \mu{\rm}$m and $100\ {\rm cm}$. We use 25 bins per size decade, which leads to 151 dust species with their size $a_{i}$. Such a dust size resolution can avoid an artificial dust size growth due to some numerical diffusion in the Smoluchowski algorithm (Smoluchowski, 1916), and ensure the convergence of dust growth rate (Ohtsuki et al., 1990).

The dust size evolution is computed within each spatial cell and is handled via an operator splitting approach between hydrodynamics and dust coagulation/fragmentation. Due to the high computational cost of solving dust coagulation, we implement a sub-stepping routine, where the dust coagulation module is executed every 50 time-steps of the hydrodynamical simulations, which still gives good calculation accuracy for the coagulation algorithm. The dust coagulation and fragmentation model follows Birnstiel et al. (2010a) with an explicit integration scheme by solving the Smoluchowski equation, which includes the radial drift and turbulent mixing as the sources of collision velocities. Turbulence is the major source of collision velocity. Since this velocity increases with Stokes number, we can derive the maximum size of grains that can grow before the impact velocity exceeds the fragmentation threshold velocity $v_{\rm f}$ as (Birnstiel et al., 2012; Pinilla et al., 2012)

[TABLE]

where $\rho_{\rm s}$ is the solid density of the dust particles. Note that this $a_{\rm max}$ is only valid for grains with ${\rm St}\leqslant 1$. This size limit for dust coagulation is also referred as the fragmentation barrier. The inverse dependence on viscosity parameter is because the turbulent velocity depends on $\alpha_{\rm vis}$.

We adopt the following procedure in simulations: Initially, only 1 $\mu$m sized dust particles are included in the disk, and the surface density distribution follows the radial profile of the gas with an initial radial-independent dust-to-gas mass ratio of $0.01$. The evolution for each dust species $a_{i}$ and the gas component is determined by the mass and momentum equations described above. We can self-consistently consider the hydrodynamics of both the gas and multiple dust species with the full dust coagulation and fragmentation. Previous studies by Pinilla et al. (2016) (see also Birnstiel et al. 2010a; Pinilla et al. 2012, 2014) have produced very interesting results on the effects of dust size growth during the disk evolution. In our current study, we treat the gas-dust dynamical interaction self-consistently, which goes beyond some of the previous work where the gas evolution is not considered (Pinilla et al., 2012, 2014). Birnstiel et al. (2010a) (see also Pinilla et al. 2016) has considered the viscous evolution of the gas, but without taking into account the full coupling (dust feedback) between the gas and dust.

Unless otherwise stated, we solve the 1D hydrodynamics equations with a logarithmically radial grid of $n_{r}=1024$. An outflow inner boundary condition and outer boundary condition are imposed on the dust and gas, which allow the gas/dust flow out and flow in from the boundary depending on their radial velocity.

2.2 Radiative Transfer

After having obtained the dust surface density distribution in disks, the dust continuum emission at different mm wavelengths was computed using the RADMC-3D package (Dullemond et al., 2012). A two-stage procedure is adopted for modeling the dust continuum radiative process. We first convert the 1D dust surface density produced from hydrodynamical simulation into a 3D distribution by assuming a dust scale height $h_{\rm d}(r)=0.1h_{\rm g}(r)$ with azimuthal symmetry as adopted in previous works (e.g., Isella et al., 2016; Liu et al., 2018)111For our fiducial model presented below, we have also tested a different dust scale height $h_{\rm d}=h_{\rm g}\min\left({1,\sqrt{\frac{\alpha_{\rm vis}}{\min(0.5,\rm St)(1+\rm St^{2})}}}\right)$ by considering dust vertical settling (Birnstiel et al., 2010a), which does not show a significant impact on the dust emission in terms of our interest. Specifically, it will lead to an increase of dust mm flux of $\sim 20\%$, and a slight change of the spectral index ($\sim 0.2$).. The vertical density structure is simply scaled as $\frac{\Sigma_{\rm d}}{\sqrt{2\pi}h_{\rm d}}\exp(-z^{2}/2h_{\rm d}^{2})$. We use 300 grids along the radial direction with a grid refinement around the bump and 40 uniform grids in the $\theta$ direction between $70^{\circ}-90^{\circ}$ with a mirror symmetry in the equatorial plane. This grid can recover all the dust mass from 1D hydrodynamics within a $5\%$ uncertainty. A larger grid number does not change the results.

We then ran RADMC-3D simulations to compute the dust temperature $T_{\rm d}(r,z)$ contributed by the stellar radiation. For simplicity, the dust temperature for all dust species is assumed to be the same, and controlled by the micro-sized dust particles.222Dust particles with different sizes should have different temperatures, with the smaller ones being slightly hotter due to a larger scale height. However, such effect should be minor since most of the emission comes from the region at the mid-plane. The dust opacity coefficient, which depends on grain sizes, chemical compositions, and shapes, is another important quantity for the modelling. We adopted the dust opacity as a function of wavelength and dust size from Isella et al. (2009) and Ricci et al. 2010a. Figure 2 presents the details of this model. We will also discuss how a different dust opacity can effect our results in the following sections (Semenov et al., 2003). The dust opacity is dominated by the diagonal region with $\lambda\simeq 2\pi a$ and $\kappa$ decreases with $a$.

With the dust temperature $T_{\rm d}$ obtained from RADMC-3D and the dust surface density $\Sigma_{\rm d}$ for each dust species derived from hydrodynamical simulations, we can calculate the continuum emission for each dust species by ray-tracing with RADMC-3D. The surface brightness of the continuum emission can be approximated as

[TABLE]

where the optical depth $\tau_{\nu}(r)=\kappa_{\nu}\Sigma_{\rm d}(r)/\cos\theta$, $B_{\nu}$ is the Planck function, $\theta$ is the disk inclination. We assume the disk inclination angle of $\theta=45^{\circ}$ and a position angle of $-30^{\circ}$ for all our models without losing generality. Using these quantities, we compute the corresponding spectral energy distribution and the radial profile of the surface brightness at several wavelengths between $0.89$ mm and $7.0$ cm by RADMC-3D.

In regions of the disk where the emission at a given wavelength is optically thick, the dust continuum flux and slope are fully determined by $B_{\nu}(T_{\rm d}(r))$, which results in a spectral slope of $2.0$ in the Rayleigh-Jeans limit. Conversely, in the optically thin regime the dust emission is related to the opacity coefficient as $F_{\nu}\propto\kappa_{\nu}B_{\nu}$. If the emission is in the Rayleigh-Jeans limit, then $F_{\nu}\propto\nu^{2+\beta}$, where $\beta$ is the slope of the dust opacity coefficient ($\kappa_{\nu}\propto\nu^{\beta}$). $\beta$ changes as dust grows bigger. $\beta$ can approach to zero once dust particle grows to size larger than a few mm, as seen in the lower panel of Figure 2. This suggests that dust growth can be an important factor for setting the spectral shape of the dust continuum in the optically thin regions of the disk (e.g., Wilner et al., 2000; Testi et al., 2001, 2003; Wilner et al., 2005; Draine, 2006).

3 Results

In this section, we describe our results for modeling the gas and dust dynamics for various disks, and discuss the comparison between the dust continuum observations of some young disks observed at mm wavelengths and the predictions of our models. Based on the comparisons with global spectral index, disk brightness with current observations, we aim at setting some constraints on the disk and coagulation parameter (e.g., disk viscosity, fragmentation velocity). The different behaviors of radial profiles for the spectral indices observed at different wavelengths in the future can also be used to break the model degeneracy in producing a low spectral index.

We first discuss the dust dynamics based on the coagulation and fragmentation model presented above. The main model parameters are listed in Table 1, otherwise their default values are adopted.

3.1 Fiducial Model

Here we start by discussing the model labeled as m1v2 in Table 1, which is characterized by $\Sigma_{0}=4.3~{}\rm{g~{}cm^{-2}}$ and fragmentation velocity $v_{\rm f}=10^{3}~{}\rm{cm~{}s^{-1}}$. The initial disk gas mass is $0.01\ M_{\odot}$.

After incorporating the low viscosity rectangular gap between 12.5$-$17.5 au as shown in Figure 1, we ran the numerical simulation for the two-fluid gas and dust evolution with dust coagulation and fragmentation for 2.1 Myr, which corresponds to 14700 orbits at 5 au. The dust settles down into an equilibrium state with a balance of coagulation and fragmentation after $\sim 2$ Myr.

The radial profile of the gas surface density is shown at different times in the upper panel of Figure 3. We generically separate this profile into the bump region and the smooth region in the following discussions. The gas bump is built up quickly in the viscosity gap region. In the case of mass conservation without disk wind, steady accretion implies that $r\Sigma_{\rm g}v_{r,\rm g}={\rm constant}$, and $v_{r,\rm g}\propto\nu_{\rm g}/r=\alpha_{\rm vis}c_{\rm s}h_{\rm g}/r$, so that $\Sigma_{\rm g}\propto\alpha_{\rm vis}^{-1}$. Therefore, a lower $\alpha_{\rm vis}$ in the gap region leads to a higher gas surface density. The minimum Toomre Q parameter ($\sim 10$) is reached in the gas bump region at the initial stage, justifying that the disk is gravitational stable during the whole evolution.

We also show the radial distribution of the total dust surface density $\Sigma_{\rm d}$, which sums up all the dust surface densities of each species $\sigma_{\rm d}$, in the middle panel of Figure 3. The dust accumulation in the gas bump region is clearly shown in the radial distribution of the dust surface density plot. We can see that the dust-to-gas ratio can increase from the initial value of 0.01 to 0.5 at the final stage in the dust trapping region, where the dust feedback could become important (Fu et al., 2014; Miranda et al., 2017).

A well-defined power-law distribution for the sizes of the dust particles in the dust trapping region is shown in the lower panel of Figure 3. The red dotted line shows the fitted power law with an index of $0.4$, which is equivalent to a power-law index of $q=3.6$ for the number density distribution $n(a)\propto a^{-q}$, close to the value of 3.5 for the dust size distribution in the ISM (Mathis et al., 1977).

The evolution of the total disk gas and dust masses is shown in Figure 4, where the masses within the “bump” region 12.5$-$17.5 au are plotted as well. As expected, the total disk mass has decreased significantly due to accretion during the 2.1 Myr evolution. The relative fraction of disk gas mass within the bump has increased. The total dust mass has reduced by about $30\%$ and, by 2.1 Myr, most of the dust mass is now residing within the bump region. The ring is essential to retain the dust mass from the rapid radial drift and loss through the disk boundary.

We further show the 2D contour plot of the dust surface density distribution for all dust species at all stellocentric radii throughout the disk in Figure 5. The fragmentation barrier, which was defined in Equation (12), is represented as the red dashed line. It can be seen that the bump in the fragmentation barrier that limits the dust size growth can be attributed to both the smaller viscosity and the larger gas surface density. The grain sizes corresponding to $a_{\rm St=1}$ (ref. to Equation (7)) are shown as the white line, which has the same slope as the gas surface density, i.e. a ringed structure around 15 au as in the upper panel of Figure 3. The fact that all the grain species in this model lie below the white line representing $a_{\rm St=1}$ indicates that ${\rm St}<1$ applies to all grains.

As expected, the dust spatial distribution shows significant trapping in the gas bump region as shown in Figure 5 (see also the middle panel of Figure 3). This is because the radial inward motion of the dust is significantly slowed down as a result of the positive pressure $(\Sigma_{\rm g}c_{\rm s}^{2})$ gradient in the inner edge of the gas bump. The large local concentration of dust grains retained in the gas bump can also facilitate the growth of dust particles from the initial 1 $\mu$m size to a few cm due to frequent collisions and sticking of small grains. Some particles above the fragmentation barrier shown in Figure 5 is due to particle diffusion. The larger particles drift faster, which leads to a sharper trap in the inner edge of the bump for the larger grains. Both dust trapping and size growth are important in determining the dust continuum emission.

Outside the gas bump region, the dust surface density decreases significantly with time due to its radial drift. In addition, the fragmentation barrier defined as in Equation (12) is significantly lower outside the bump region than inside the bump due to the radial viscosity profile we adopt, which leads to the average dust size (weighted by the dust surface density) being much smaller in the smooth region of the disk. This also explains an increase of dust size with a decreasing stellocentric radius $r$ because the dust surface density becomes increasingly higher towards the inner region of the disk which results in a higher fragmentation barrier.

We then calculate the dust temperature using RADMC-3D as described in Section 2.2. The radial profile of the temperature in the mid-plane of the disk is shown in Figure 6. The temperature profile is well consistent with the $r^{-0.5}$ scaling expected for a disk heated by the central star except for a bump around 15 au (Armitage, 2010). The dust temperature profile is also roughly in agreement with the profile in our hydrodynamical simulation based on Equation (2), which is shown as black dotted line in Figure 6. The temperature bump is associated with the dust bump feature in the same region, which results from the more efficient absorption (Bjorkman & Wood, 2001). A higher dust temperature in the bump could decrease the dust size in that region. But this is a small modification in dust size so it cannot significantly influence the dust emission, since the current maximum size is much bigger than $\sim 1{\rm mm}$. It is expected that the temperature in the mid-plane is lower than that at the disk surface in the outer disk region, while the temperature in the inner edge of the disk is close to the temperature at the disk surface determined by the stellar heating. This results in the steep slope for the mid-plane temperature profile in the inner disk. The total emission from the dust in this very inner region is, however, negligible due to the significant dust depletion as discussed above.

Due to the dust accumulation in the gas bump region, the local effective optical depth $\tau_{\rm eff}$ increases significantly. The radial profiles of the effective optical depth $\tau_{\rm eff}$ for all dust species, which are obtained by summing the optical depth over all dust species, are shown in the upper panel of Figure 7. We show the optical depth at two wavelengths, i.e., 1 and 7 mm, and at two epochs, close to the initial stage (0.04 Myr) and at the end of the simulation (2.1 Myr). At the wavelength of 1 mm, the dust-trapping region is optically thick, while other regions are still optically thin at the end of the simulation (2.1 Myr). At the longer wavelength of 7 mm, the emission becomes marginally optically thick even in the dust trapping region. Note that the optical depth in the ring is significant higher than that inferred from DSHARP samples (Dullemond et al., 2018). A lower gas surface $\Sigma_{0}$ can make a better agreement with these observations. Furthermore, the inclusion of dust scattering can result in an optically thick disk appearing as optically thin (Zhu et al., 2019). The dust optical depth in the trapping region does not evolve significantly with time, however, it decreases significantly outside the trapping region at mm wavelengths. This is mostly due to the dust radial drift effect.

After we obtained the dust size and radial distribution, its temperature, and optical depth for all dust species, we can obtain the radial profile of the dust continuum emission at different sub-mm and mm bands with a detailed radiative transfer calculation using RADMC-3D, which is shown in the lower panel of Figure 7. The flux at different bands can be obtained by integrating the dust continuum emission over the whole disk. We can then calculate the corresponding global spectral index and the 1mm dust continuum flux. The results for model m1v2 are shown in Table 1 and Figure 8. In Figure 8, we show the time variation of the spectral index calculated between three different sub-mm/mm wavelengths as a function of the 1 mm flux density. The 1 mm flux decreases significantly at the early stage of the evolution, and then show less temporal variations after $\sim$1.7 Myr. Such a decreasing flux is a consequence of the dust radial evolution as shown in the middle panel of Figure 3. In the early stages of the disk evolution, the region outside the trapping region contributes a comparable fraction to the total flux of the disk. The outer disk contribution is significantly reduced as time evolves to $\sim$1.7 Myr. This is also confirmed by the radial profile of the optical depth as shown in the upper panel of Figure 7. The emission outside the bump region could be initially optically thick for a massive disk (as is the case for model m1v2, see the upper panel of Figure 7), and it then becomes optically thin due to the dust radial drift. The spectral slope from the optically thin emission in the smooth region can be as large as about 3.5 since its value is determined by a population of small-sized dust, as shown in the lower panel of Figure 7. Therefore, the global spectral index could stay at a low value initially due to a large fraction of optically thick emission. It then increases with time rapidly when the emission begins to transfer to the optically thin regime, as long as the flux contribution from the smooth region is still comparable to that from the bump. As disk evolves further, the contribution from the smooth region becomes increasingly negligible, which means that the global spectral index is mainly determined by the emission from the bump, then the global spectral index decreases with a decreasing flux as shown in Figure 8.

In addition to the spectral index $\alpha_{\rm 1.0-3.0mm}$, we further calculate the index between $0.89-1.3$ mm, i.e., $\alpha_{\rm 0.89-1.3mm}$, which can be compared with the observations from the ALMA Lupus survey (Ansdell et al., 2018). It can be seen that the time evolution patterns of the spectral index in these two spectral windows are quite similar due to the closeness of the two windows.

Note that the opacity we have adopted (Isella et al., 2009; Ricci et al., 2010a) is higher by a factor of a few ($\sim 5$) than that adopted in other works (e.g., Semenov et al., 2003; Birnstiel et al., 2018). To test the sensitivity of our spectral results on the choice of dust opacity, we adopt a norm silicates (NRM) dust opacity model333http://www2.mpia-hd.mpg.de/homes/henning/Dust_opacities/Opacities

/RI/new_ri.html. It turns out that the global spectral index around 1 mm only changes only slightly and the final disk flux $F_{\rm 1mm}$ decreases by a factor of $\sim 2$. The decrease of the total flux is mainly because the optical depth in the trapping region becomes lower. Such modifications do not significantly change our main conclusions. In the following discussion, we will only adopt the opacity from Isella et al. (2009; Ricci et al. 2010a) unless otherwise stated.

3.2 Dependence on Parameters

In this section, we study how our main results depend on the model parameters we have used. The main parameters we will explore include the fragmentation velocity $v_{\rm f}$, the initial gas surface density $\Sigma_{0}$, disk global viscosity $\alpha_{0}$, viscosity gap (ring) width $w_{\rm ring}$, gap (ring) depth $d_{\rm ring}$, gap (ring) location $r_{\rm ring}$, disk scale height $h_{0}$, and disk gas profile index $\gamma$, which are listed in Table 1. In the following, we evolve all the disk models for 2.1 Myr when the dust/gas dynamics has reached the equilibrium state.

3.2.1 Fragmentation Velocity

We first study the effect of fragmentation velocity $v_{\rm f}$ of the dust coagulation model. We ran several models with $v_{\rm f}$ in the range of $1\times 10^{2}$ to $3\times 10^{3}\ {\rm cm\ s^{-1}}$ while fixing all other parameters. The final global spectral index and 1 mm flux for models with varying fragmentation velocity are shown in the upper left panel of Figure 9.

In comparison with the case of $v_{\rm f}=1\times 10^{3}\ {\rm cm\ s^{-1}}$ in Figure 8, the spectral index around 1 mm $\alpha_{\rm 1.0-3.0mm}$ is insensitive to $v_{\rm f}$ until $v_{\rm f}$ decreases to $1\times 10^{2}\ {\rm cm\ s^{-1}}$. The spectral index $\alpha_{\rm 1.0-3.0mm}$ for a very small $v_{\rm f}=1\times 10^{2}\ {\rm cm\ s^{-1}}$ can reach a large value of $\sim 3.5$. The two effects contribute to this large index. As $v_{\rm f}$ becomes very small (i.e., lowering the fragmentation barrier), it limits the maximum grain size obtained from coagulation. With the decrease of $v_{\rm f}$ by one order of magnitude, the maximum dust size decreases by two orders of magnitude, giving the largest grain sizes around only $\sim 2$ mm. In addition, without the large dust grains, the dust emission contributed by the small dust becomes optically thin. Such a large population of small-sized grains with optically thin emission thus produces emission with a high spectral index.

When the fragmentation velocity becomes even larger ($v_{\rm f}=3\times 10^{3}\ {\rm cm\ s^{-1}}$), most of the dust mass now resides in the largest dust grains, reducing the overall surface density in smaller grains. Because opacity is typically dominated by the small grains, the effective optical depth can get significantly reduced, i.e., becoming mostly optically thin. As a result, the total mm flux, which is dominated by the dust bump region, also decreases significantly when $v_{\rm f}$ becomes large enough, as shown in the upper left panel of Figure 9.

The spectral index $\alpha_{\rm 0.8-1.3mm}$ for different $v_{\rm f}$ is also shown in Figure 9, which is similar to $\alpha_{\rm 1.0-3.0mm}$. The sub-mm/mm fluxes for most of the faintest disks in the Lupus sample can be reproduced by models with the high fragmentation velocity $\sim 30\ {\rm m\ s^{-1}}$. However, it is difficult to use models with high fragmentation velocity to explain the disks with low spectral indices as discussed below.

3.2.2 Disk Mass

Another important disk parameter is the disk mass, which, in our models, is determined by the gas surface density $\Sigma_{0}$. In addition to the fiducial model, we chose three other values for $\Sigma_{0}$ as listed in Table 1. The most massive disk corresponding to $\Sigma_{0}=12.8\ {\rm g\ cm^{-2}}$ (model labels including m0) has a total disk mass of $0.03\ M_{\odot}$. We also find that, for this model, the minimum Toomre Q parameter in the gas bump region can approach 3 during its evolution. The dependence of the total disk flux $F_{\rm 1mm}$ on $\Sigma_{0}$ is rather straightforward. An increasing $\Sigma_{0}$ will result in a higher $F_{\rm 1mm}$ as shown in the upper right panel of Figure 9 because more dust can be available initially. But the increasing flux is not linearly related with $\Sigma_{0}$, especially when the emission is dominated by the optically thick regions. We find that the disk flux is very likely saturated at a flux level of $\sim 100$ mJy, as shown in Table 1 (and also the upper right panel of Figure 9). In this optically thick regime, the global spectral index is insensitive to $\Sigma_{0}$ since this is determined solely by the spectral slope of the Planck function, and does not depend on the size of the emitting dust.

When the initial dust surface density is decreased by two orders of magnitude through the decrease of $\Sigma_{0}$, the dust surface density becomes low globally, and most of the dust mass is retained in the small particles. As a result, the whole disk, even the dust trapping region, becomes optically thin. At the same time, due to the lower fragmentation barrier $a_{\rm max}$, which linearly depends on the gas surface density $\Sigma_{\rm g}$, the maximum dust size also decreases with the decreasing $\Sigma_{0}$. In this case, not only the dust continuum flux decreases with $\Sigma_{0}$, but also the global spectral index $\alpha_{\rm 1.0-3.0mm}$ increases to values reflecting the dust opacity index $\beta$ of grains with sizes much smaller than $\sim 1$ mm.

We have also made a comprehensive parameter study for the combined effects of different $v_{\rm f}$ and $\Sigma_{0}$, and list their results in Table 1 and Figure 10.

We find several trends that are described below. For a fixed $v_{\rm f}$, the dependence of $\alpha_{\rm mm}$ on $\Sigma_{0}$ is as follows. When $v_{\rm f}$ is as high as $3\times 10^{3}~{}{\rm cm\ s^{-1}}$, $\alpha_{\rm mm}$ is insensitive to $\Sigma_{0}$ as the fragmentation barrier $a_{\rm max}$ is already larger than 1 cm even for the lowest values of $\Sigma_{0}$ explored in this study. However, the global spectral index is usually larger than $2.5$, except for extremely low mass disk (i.e., model m3v1). Such a low mass disk is below the detection limit of current observations. While for the case of low $v_{\rm f}$, $\alpha_{\rm mm}$ (1 mm flux $F_{\rm 1mm}$) is anti-correlated (correlated) with $\Sigma_{0}$ due to the optical depth effect. Based on our parameter survey for different $\Sigma_{0}$ and $v_{\rm f}$ listed in Table 1, the tentative anti-correlation between $\alpha_{\rm mm}$ and $F_{\rm 1mm}$ found by (Ansdell et al., 2018) for some faint sources, as well as the low spectral index for some slightly bright disks from different surveys (Ricci et al., 2012; Ansdell et al., 2018), is consistent with a fragmentation velocity in the range of $3\times 10^{2}\sim 1\times 10^{3}~{}{\rm cm\ s^{-1}}$. This conclusion still holds when we use the opacity from Semenov et al. (2003). As there are still some uncertainties for the dust opacity and disk parameters (the viscosity profile and the associated gas bump structures as we will discuss later), we should point out that such a constrain on the fragmentation velocity could be not so stringent.

To summarize, both a high fragmentation velocity $v_{\rm f}$ and high dust mass can produce a global spectral index $\alpha_{\rm mm}$ smaller than 2.5, which are consistent with the mm observations for some young PPDs (Ricci et al., 2012; Ansdell et al., 2018), as shown in Figure 9. Such a model degeneracy can be broken by checking the variation of the spectral index radial profile with different wavelengths, which will be discussed in Section 3.3. When both the fragmentation velocity and dust mass are below a critical value such that the fragmentation barrier $a_{\rm max}<1\ {\rm mm}$ and the effective optical depth $\tau_{\rm eff}\ll 1$, the global spectral index $\alpha_{\rm mm}>3.0$.

3.2.3 Viscosity

The viscosity is expected to influence the gas bump properties and modify the dust dynamics as well. We simulate several models with different viscosity parameters $\alpha_{0}$. The spectral index $\alpha_{\rm mm}$ at different bands and 1 mm flux are listed in Table 1 and are shown in Figure 9.

According to Equation (12), the maximum dust size increases as viscosity decreases. This affects the optical depth significantly. Indeed, as shown in Figure 11 where the run is similar to our fiducial run except that $\alpha_{0}=5\times 10^{-3}$ is a factor of 2 lower, the fragmentation barrier in the whole disk $a_{\rm max}$ is larger. This leads to two consequences. One is that $F_{\rm 1mm}$ decreases with decreasing $\alpha_{0}$ because the dust in the whole disk can coagulate to a larger size thanks to the larger fragmentation barrier, which results in the emission shifting from the optically thick to thin regime for most of the disk region. The other is that $\alpha_{\rm mm}$ becomes increasingly larger. It turns out that the dust bump region is still optically thick and it gives a spectral index around $2.0$, but the flux aside from the bump region can contribute a considerable fraction of the total emission. This flux has a spectral index close to $2.5\sim 3.5$ because of the optically thin emission. The global spectral index in the whole disk then gives a value of about $3.0$.

In order to produce brighter disks, a relatively high viscosity parameter $\alpha_{0}$ is thus required. Note that the increase of global spectral index $\alpha_{\rm mm}$ is mainly attributed to the global increase of the fragmentation barrier $a_{\rm max}$. In order to produce a relatively low spectral index with a low disk viscosity, we need to adopt a lower fragmentation velocity $v_{\rm f}$ to limit the dust sizes.

3.2.4 Viscosity Gap Depth, Width and Location

We now study how viscosity gap depth affects the emission properties. We used three $d_{\rm ring}=2.0,10.0,20.0$ with fixed $\alpha_{0}$, which are shown as the models m1d1, m1d2, and m1d3 in Table 1 and Figure 9.

Overall, $\alpha_{\rm mm}$ is insensitive to $d_{\rm ring}$. This is because the emission in the trapping region can either be optically thick or dominated by larger dust particles, both of which lead to a shallow spectral shape, while the emission outside the bump region does not vary significantly. Furthermore, we find that even a slightly shallower gap for the viscosity profile, i.e., $d_{\rm ring}=2.0$, can still trap enough dust and result in dust size growing. The value of global spectral index $\alpha_{\rm 0.8-1.3mm}$ increases slightly, but is still smaller than 3.0. In order to reproduce a spectral index lower than 2.5, the gap depth for the gas viscosity, or equivalently the gas density contrast, should be larger than $\sim 3$.

In the situation where the gas bump is produced via a mechanism different from a viscosity gap, such a gap should not be used in the coagulation/fragmentation calculation. This is because the real viscosity should be smooth in the whole disk, and our viscosity gap is only used for generating the bump. This has been tested for our fiducial model. We find that the maximum dust size reached in the bump region decreases by a factor of 6, which is $\sim 2\ {\rm mm}$. The mm flux and the spectral index are 42.2 mJy and 2.46, respectively, which are actually close to our fiducial model results.

We have also adjusted the gap width by using different $w_{\rm ring}$, $2.5$ au, $7.5$ au, and $10.0$ au, of the viscosity profile. Results are listed in Table 1 and shown in Figure 9. We find that the final spectral slopes and $F_{\rm 1mm}$ are insensitive to changes of $w_{\rm ring}$ between $2.5~{}{\rm au}$ and $10~{}{\rm au}$.

The dependence of disk emissions on the ring location is investigated by placing the viscosity gap at $r_{\rm ring}=$ $25~{}{\rm au}$, $50~{}{\rm au}$, and $75$ au. The results are shown in Figure 9. In general, the 1 mm flux increases slightly and the spectral slope increases when $r_{\rm ring}$ moves outward. This is because the dust mass interior to the dust bump region, $r\lesssim r_{\rm ring}-w_{\rm ring}/2$, gets larger. This region has a relatively steep spectral slope because it is optically thin emission from small dust grains. The emission from the dust trapping region can still be optically thick. Therefore, as $r_{\rm ring}$ increases, both the global spectral index and the total disk flux increase. As shown in Figure 9, the global spectral indices from these models give somewhat higher values than the observations.

3.2.5 Gas Scale Height and Surface Density Profile

We further consider two other disk parameters, gas scale height $h_{0}$ and gas surface density profile $\gamma$. The gas scale height can influence the gas temperature (pressure) profile. We show one example for the variation of the mid-plane dust temperature with the disk scale height in Figure 6. The temperature is slightly higher because of a higher sound speed. The $\gamma$ parameter can change the total amount of dust. The models labeled as m1h1(2) and m1b1(2) with the corresponding results are shown in Table 1. For all these models, the spectral index and mm flux variations are only modest. This is because the dust trapping region is mostly optically thick, and in the optically thick regime the emission depends on temperature, not on the grain size distribution. The small variation of the flux and spectral index is due to the small variation in the dust temperature in the dust trapping region as shown in Figure 6.

3.2.6 No Viscosity Gap – No Ring

In order to quantify the role of the viscosity gap/ring, we study the case with $\alpha_{\rm vis}=\alpha_{0}=10^{-2}$ throughout the whole disk, and all the other disk parameters the same as in our fiducial model m1v2. At the end of the 2.1 Myr evolution, the total dust mass is only $1.3\times 10^{-6}M_{\odot}$, or $\leq 1\%$ of the original dust is retained in the disk and is two orders of magnitude smaller than for model m1v2, as shown in the left panel of Figure 12. In addition, low dust surface density severely limits the ability of dust particles growing to larger sizes, as shown in the right panel of Figure 12. In fact, almost all the bigger grains will drift radially inward due to their larger Stokes number.

We then repeat the same radiative transfer calculations to produce the dust continuum emission for this disk model. The obtained values for the spectral index for the emission within 100 au from the star are quite large, reaching values of $\approx 4.0$ due to an inefficient size growth and an optically thin emission for the dust. As expected, the flux at 1 mm is only 5.0 mJy due to the low total dust mass. It is difficult to decrease the global spectral index down to the observed values by changing fragmentation velocity and/or initial dust mass because the disk emission is dominated by sub-mm sized particles in the outer region of the disk.

3.3 Differentiating Effects of Optical Depth and Dust Size

Both optically thick emission and efficient dust growth can be responsible for the low values of the global spectral indices, and one key question is how to discriminate between these two scenarios. We resort to the spatially resolved radial profile of the disk emission at multiple wavelengths to break such degeneracy.

To quantify our method, we choose two models with different values of fragmentation velocity $v_{\rm f}$ to compare with our fiducial model m1v2. We convolve the image produced from RADMC-3D with a Gaussian beam of the FWHM $0.003^{\prime\prime}\times 0.003^{\prime\prime}$, close to the full resolution of our RADMC-3D grid, which can recover all the substructures produced in the hydrodynamical model. Although such resolution is beyond the current limit of ALMA observations, our goal is to diagnose the physical processes that contribute to different spectral shapes of the dust emission.

The radial profile of flux and spectral index at mm wavelengths for model m1v1 (high $v_{\rm f}$) are shown in the upper panels of Figure 13. Most of the disk is optically thin for both 1 mm and 7 mm, which is responsible for the relatively high spectral index throughout the disk. There exists a narrow region within the dust ring where it is marginally optically thick and the spectral index dips to $\sim 2.0$ at that location. This is caused by the biggest particles which dominate the emission have grown to be much larger than $\sim 7/2\pi\ {\rm mm}\simeq 1.0\ {\rm mm}$. The sharp peak near the inner edge of the dust ring in the optical depth profile is contributed to the large particles which are sharply trapped there. The overall radial profile of spectral index is somewhat similar to model m1v2 case except that the gap in the viscosity profile becomes narrower, especially at the longer wavelength, as shown in the upper panel of Figure 13.

For a lower fragmentation velocity model m1v3, the emission can be optically thick at 1 mm, but optically thin at 7 mm band in the dust trapping region, as shown in the lower left panel of Figure 13. This is also contrary to the case of m1v2 shown in the upper panel of Figure 7, where the emission in the dust trapping region is optically thick both at 1 mm and 7 mm. We show the corresponding radial profile of the spectral index and flux in the lower right panel of Figure 13. The lowest spectral index between $0.8-1.3$ mm is about 2.0, but at longer wavelengths, e.g. at $3.0-7.0$ mm, it reaches much higher values in the same region. This reason is as follows. First, this is due to the optically thin emission at longer wavelength as shown in the lower left panel of Figure 13. Second, in the optically thin regime, the grains cannot grow to a size larger than $1.0\ {\rm mm}$ because of the lower fragmentation velocity.

As a comparison to our fiducial model with $v_{\rm f}=1\times 10^{3}\ {\rm cm\ s^{-1}}$ (model m1v2), the spectral index at the dust trapping region can reach 2.0 both at short waveband ($0.8-1.3$ mm) and long waveband ($3.0-7.0$ mm), which is a consequence of the large optical depth.

As shown in these cases, the different radial behavior of the spectral index at multiple mm wavelengths can be used to break the model degeneracy as well as set some constraints on the dust size444These constraints will be affected by the uncertainties in the dust opacity or the unknown porosity of the particles.. This can be tested by future high resolution observations of ALMA and ngVLA at longer wavelengths.

4 Summary and Discussions

In this work, we present detailed 1D hydrodynamical simulations of PPDs to study the two-fluid (gas+dust) co-evolution with a state-of-the-art dust evolution code (Li et al., 2005, 2009; Fu et al., 2014; Birnstiel et al., 2010a) with dust size growth via coagulation and fragmentation (Birnstiel et al., 2010a). This allows us to self-consistently obtain the quasi-steady state gas and dust distributions including dust feedback, dust radial drift, and coagulation/fragmentation balance. The comparisons with observations are then made using the RADMC-3D code for dust continuum radiative transfer.

By assuming a viscosity gap in the inner region of the disk ($\sim 15$ au), a gas bump as well as a dust trap are produced. With a systematic parameter study for different disk gas and dust parameters and coagulation velocities, we find:

• •

A gas and dust ring at the inner disk can produce a lower global spectral index.

• •

High fragmentation velocities $v_{\rm f}$ ($1000~{}{\rm cm~{}s^{-1}}\lesssim v_{\rm f}\lesssim 3\times 10^{3}~{}{\rm cm~{}s^{-1}}$) can facilitate efficient dust size growth (Figure 5), which produces a dust mm spectral index $\alpha_{\rm mm}$ close to 2.0, even for optically thin emission. For lower values of $v_{\rm f}$ (300 to $1000~{}{\rm cm~{}s^{-1}}$), a high dust surface density in the dust trapping region leads to optically thick emission, which then results in a low $\alpha_{\rm mm}$ as well (Figure 8, Figure 9).

• •

The high flux and low spectral index sources from a few surveys (Ricci et al., 2012; Ansdell et al., 2018) can be explained by a fragmentation velocity $v_{\rm f}$ in the range $300-1000~{}{\rm cm~{}s^{-1}}$. The tentatively negative relation between $\alpha_{\rm mm}$ and $F_{\rm 1mm}$ shown from Ansdell et al. (2018) for the faint PPD sources can also be interpreted by a relatively low $v_{\rm f}$ ($\lesssim 1000~{}{\rm cm~{}s^{-1}}$) with a variation of dust mass surface density. This is simply because the optical depth becomes smaller as the dust mass, and mm flux, decrease (Figure 9).

• •

The spectral index at mm wavelengths $\alpha_{\rm mm}$ is not very sensitive to the gas bump width $w_{\rm ring}$, gap depth $d_{\rm ring}$, disk scale height $h_{0}$, and disk profile index $\gamma$. A very low global disk viscosity $\alpha_{0}\lesssim 10^{-3}$ cannot reproduce the observed low global spectral index.

• •

The one-ringed disk models can reproduce the observed properties of the faint part of the observed disks ($\lesssim 100$ mJy) reasonably well, but they cannot reproduce the brightest disks in the observational surveys (Figure 9). We expect that multiple rings could be responsible for these bright disks as discussed later.

• •

While both large size grains and optically thick emission contribute to a small spectral index $\alpha_{\rm mm}\lesssim 2.5$, the future high resolution observations, e.g., by ALMA and ngVLA, can distinguish these two effects based on their spatial distribution properties (Figure 7, Figure 13).

In this work, we adopt a viscosity parameter $\alpha_{\rm 0}=10^{-2}$ for most of our disks. This value is a factor of a few higher than some observational constrained values based on the vertical thinness of some disks with rings and gaps, e.g., HL Tau (Pinte et al. 2016; $\alpha_{\rm vis}\sim 3\times 10^{-4}$) and HD 163296 (Flaherty et al. 2017; $\alpha_{\rm vis}\lesssim 3\times 10^{-3}$). On the one hand, such a global low $\alpha_{\rm vis}$ could reflect the low turbulence in the ringed/gap regions. Therefore, it indicates that some regions in the disk could be indeed magnetorotational instability (MRI) inactive compared with the theoretical expectation of MRI. On the other hand, for a lower viscosity $\alpha_{\rm vis}$ in a very thin disk, we find that an even lower fragmentation velocity ($\lesssim 3\ {\rm m\ s^{-1}}$) is required to produce a low spectral index ($\alpha_{\rm mm}\lesssim 2.5$). One issue we should point out is that we keep the viscosity profile $\alpha_{\rm vis}(r)$ independent on dust evolution. This is contrary to the expectation of the MRI dead zone scenario. Dzyurkevich et al. (2013) has shown that $\alpha_{\rm vis}$ depends on ionization fraction, which depends on dust properties.

We have ignored the dust scattering effect in the radiative transfer. Recently, Zhu et al. (2019) (see also Liu 2019) found that dust scattering can reduce the emission from an optically thick region. It leads to the apparent interpretation that bright disks tend to be harder. It can also modify the spectral slope at mm band depending the variation of the dust albedo with wavelength. However, the modification of the spectral shape is not significant (close to the expectation of 2.0) if the dust size is much larger than $\lambda/2\pi$, which is the case in the dust trapping region in most of our models. This issue will be addressed with detailed radiative transfer calculations in the future.

We have also explored the effects of multiple rings. A very low global spectral index (smaller than 2.5) can be obtained with a reasonable combination of disk parameters and fragmentation velocity ($\sim 10\ {\rm m\ s^{-1}}$), and more importantly, the disk becomes systematically brighter. Multiple rings can increase the fraction of optically thick emission regions, thus significantly increase the disk total flux, while still keeping low values for the global spectral indices around mm wavelength. These are consistent with recent ALMA surveys for 12 disks in Taurus star-forming region (Long et al., 2018) and 18 single-disk systems in the DSHARP program (Andrews et al., 2018; Huang et al., 2018b), where disks with multiple sub-structures tend to be brighter, although an observational bias might exist.

The ringed structures we inferred could be the location to trigger further planetesimal formation. If the ringed structure is associated with the MRI dead zone, the less turbulence in the gas bump would result in lower collision speeds for particles, which would be more suitable for further grain growth. Such a gas bump can also slow down the rapid radial drift to allow more time for planetesimal formation. In addition, the concentration of particles in high densities can trigger planetesimal formation via streaming instability (Youdin & Goodman, 2005) and the gravitational collapse (e.g., see review by Pinilla & Youdin, 2017). The planetesimal formation process can remove some mm-size dust particles, and hence reduce the mm flux of the disk, which makes our model even harder to explain the brighter disks.

We would like to thank Ruobing Dong and Pinghui Huang for helpful discussions, and the referee for useful comments on the paper. YPL, HL and SL gratefully acknowledge the support by LANL/CSES and NASA/ATP. T.B. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 714769 and funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Ref no. FOR 2634/1. FY acknowledges the support by Natural Science Foundation of China under grant 11661161012. Simulations of this work were performed with LANL Institutional Computing.