Free Access
Issue
A&A
Volume 637, May 2020
Article Number A5
Number of page(s) 13
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201937198
Published online 01 May 2020

© ESO 2020

1 Introduction

Circumstellar disks are the cradles of planet formation. These disks, containing gas and dust particles, are observed to live for up to ~ 10 Myr before being dispersed (e.g., Williams & Cieza 2011). Planetary systems are assembled from the gas and dust in circumstellar disks. The evolutionary scenario of the disk can have a strong influence on the distribution of the dust content in the disk. Dust plays a critical role not only in the disk evolution and planet formation, but it also determines the observational appearance of the disks. One of the fundamental disk parameters that is important, not only for planet formation, but also for observations, is the dust-to-gas ratio in the disk (e.g., Johansen et al. 2007). Typically, the constant dust-to-gas ratio is used for disk mass derivation from continuum observations, giving the biggest source of uncertainty in disk mass estimates (Andrews & Williams 2007; Birnstiel et al. 2012).

However dust and gas are transported in the disk differently, thus leading to the order-of-magnitude local deviations of dust-to-gas ratio in the disk from the canonical 1:100 value (Vorobyov & Elbakyan 2019). Dust particle transport is an essential component of the pebble accretion model, where gas-giant cores are formed from the cm-sized pebbles that have been decoupled from gas (Ormel & Klahr 2010; Lambrechts & Johansen 2012). One of the main problems with planet formation theory is the radial migration of dust particles on timescales that are much shorter than their growth timescales (e.g., Weidenschilling 1977). The trapping of dust particles at the radial pressure bumps in the disk can prevent the dust particles from carrying out fast inward migration. (Whipple 1972; Haghighipour & Boss 2003). Radial pressure bumps could be formed in different parts of cicumstellar disks – at the snow line, at the dead zone inner edge, and at the disk inner edge (Drążkowska et al. 2013; Johansen et al. 2014; Charnoz et al. 2019).

Due to the small initial angular momentum or inefficient angular momentum transport in the disk, the resulting disk can become massive and prone to gravitational instability (Vorobyov & Basu 2005, 2006; Rice et al. 2010; Zhu et al. 2010). Recent observations have shown disks with spiral arms, indicating that they are possibly massive and gravitationally unstable (Pérez et al. 2016; Huang et al. 2018; Dullemond et al. 2018). Gravitational instability and fragmentation of disks are considered as one of the possible mechanisms of giant planet and brown dwarf formation (Boss 1998; Vorobyov & Basu 2010; Zhu et al. 2012; Stamatellos 2015). Gaseous clumps forming in the disk serve as a local pressure bumps. Dust particles accumulated inside the gaseous clumps grow rapidly and possibly form solid cores that could eventually become giant, icy, or terrestrial planets (Boley et al. 2010; Vorobyov 2013; Nayakshin 2017; Vorobyov & Elbakyan 2019).

In this paper, we study the evolution and dynamics of dust particles in the circumstellar disk and its surrounding envelope. We use FEOSAD numerical hydrodynamics code that allows us to study the formation and long-term evolution of circumstellar disks. Unlike other studies, where only the evolution of the disk with fixed size dust particles is considered, here we also take into account the dust growth. Moreover, the evolution of disk is considered inside the envelope, self-consistently; thus, the envelope serves as a reservoir of the small dust particles, which are constantly growing and migrating inwards, not allowing the dust in the disk to be exhausted. The studies of disk evolution starting from the pre-stellar core collapse phase are essential for understanding how the building blocks of planets (chondrules, calcium-aluminum-rich inclusions, etc.) are formed (Charnoz et al. 2015; Pignatale et al. 2018).

The paper is structured as follows: In Sect. 2, we introduce the numerical models used. Section 3 presents the main results obtained and in Sect. 4, we discuss the parameter space study. Our findings are summarized in Sect. 5. In Appendix A, we test the dependence of our conclusions on the choice of fragmentation velocity and turbulence strength.

2 Numerical model

In this section, we briefly describe the numerical hydrodynamic model used in this paper. The formation and evolution of a star and its circumstellar disk is studied using the FEOSAD code described in detail in Vorobyov et al. (2019). Here, we only briefly review its main features and equations.

Our simulations start from the gravitational collapse of a pre-stellar core of a certain mass, angular momentum, and dust-to-gas ratio and continue into the embedded phase of stellar evolution, during which the protostar and protostellar disk are formed. We introduce a “central smart cell” (CSC) at the coordinate origin with a radius of 1 au to avoid small time-steps and save computational time. The matter is accreted from the CSC on the star with two modes: regular and burst. During the regular accretion mode, it is assumed that the mass accreted from the CSC onto the star is a fraction ξ of the mass accreted from the disk to the CSC. For the CSC we use inflow-outflow boundary condition, in which the matter is allowed to flow freely from the disk to the CSC and vice versa. Thus, the mass accretion rate through the CSC-disk interface (disk) can be both positive (when the matter flows from the disk to the CSC) and negative (when the matter flows from the CSC to the disk). The mass accretion rate from the CSC onto the star during the regular accretion mode is calculated as M˙*,csc={ξM˙disk,if M˙disk>00,if M˙disk0 .\begin{equation*}\dot{M}_{\textrm{*,csc}} = \begin{cases} \xi \dot{M}_{\textrm{disk}}, & \mbox{if } \dot{M}_{\textrm{disk}}>0 \\ 0, & \mbox{if } \dot{M}_{\textrm{disk}}\leq0 \end{cases}. \end{equation*}(1)

Vorobyov et al. (2019) showed that the mass transport rate inside the CSC affects the evolution of protostellar disks. More specifically, a slow mass transport inside the CSC leads to the formation of more massive, warmer disks with dust particles reaching a decimeter in radius. In all our models, presented in this paper, we use ξ ≈ 1.0, meaning that the mass is transported with a similar efficiency both in the disk and CSC. We plan to consider the mass fluxes in the disk with the different mass transport rates inside the CSC in a follow-up study. More information on burst accretion mode and the model of central smart cell in general can be found in Vorobyov et al. (2019).

For the outer boundary of the computational domain, the free outflow boundary conditions are imposed so that the matter is allowed to flow out of the computational domain, but is prevented from flowing in. We use the polar coordinates (r, ϕ) on a two-dimensional numerical grid with 256 × 256 grid zones. The radial grid is logarithmically spaced, while the azimuthal grid is uniformly distributed. The simulations are terminated in the T Tauri phase of disk evolution when the age of the system reaches 1.0 Myr.

The evolution of the disk and the envelope are calculated taking into account the viscous and shock heating, irradiation from the central star and from the background, dust radiative cooling from the disk surface, momentum exchange between gas and dust, self-gravity of gaseous and dusty disks, and turbulent viscosity using the α-parametrization (Shakura & Sunyaev 1973). The equations of mass, momentum, and energy transport for the gas component are: Σgt+p(Σgvp)=0,\begin{equation*}\frac{{\partial \Sigma_{\textrm{g}} }}{{\partial t}} + \nabla_{\textrm{p}} \cdot ( \Sigma_{\textrm{g}} \bm{v}_{\textrm{p}} ) = 0, \end{equation*}(2) t(Σgvp)+[(Σgvpvp)]p=pP+Σggp+(Π)pΣd,grfp, \begin{align*}\frac{\partial}{\partial t} \left( \Sigma_{\textrm{g}} \bm{v}_{\textrm{p}} \right) + \left[\nabla \cdot \left(\Sigma_{\textrm{g}} \bm{v}_{\textrm{p}} \otimes \bm{v}_{\textrm{p}} \right)\right]_{\textrm{p}} =& - \nabla_{\textrm{p}} {\cal P} + \Sigma_{\textrm{g}} \, \bm{g}_{\textrm{p}} \nonumber\\ &\,+ (\nabla \cdot \mathbf{\Pi})_{\textrm{p}} - \Sigma_{\textrm{d,gr}} \bm{f}_{\textrm{p}}, \end{align*}(3) et+p(evp)=P(pvp)Λ+Γ+(v)pp:Πpp,\begin{equation*}\frac{\partial e}{\partial t} +\nabla_{\textrm{p}} \cdot \left( e \bm{v}_{\textrm{p}} \right) = -{\cal P} (\nabla_{\textrm{p}} \cdot \bm{v}_{\textrm{p}}) -\Lambda +\Gamma + \left(\nabla \bm{v}\right)_{\textrm{pp}^{\prime}}:\Pi_{\textrm{pp}^{\prime}}, \end{equation*}(4)

where the subscripts p and p′ refer to the planar components (r, ϕ) in polar coordinates; Σg is the gas surface density; e is the internal energy per surface area; P${\cal P}$ is the vertically integrated gas pressure calculated via the ideal equation of state as P=(γ1)e${\cal P}=(\gamma-1) e$ with γ = 7∕5; vp=vrr^+vϕϕ^$\bm{v}_{\textrm{p}}=v_r \hat{\bm r}+ v_{\phi} \hat{\boldsymbol \phi}$ is the gas velocity in the disk plane; and p=r^/r+ϕ^r1/ϕ$\nabla_p=\hat{\bm r} \partial / \partial r + \hat{\boldsymbol \phi} r^{-1} \partial / \partial \phi $ is the gradient along the planar coordinates of the disk. The term fp is the drag force per unit mass between dust and gas, describing the back-reaction of dust on gas according to the method described in Stoyanovskaya et al. (2018).

The gravitational acceleration in the disk plane, gp=grr^+gϕϕ^$\bm{g}_{p}=g_r \hat{\bm r} +g_{\phi} \hat{\boldsymbol \phi}$, takes into account self-gravity of the gaseous and dusty disk components found by solving the Poisson integral (see details in Vorobyov & Basu 2010) and the gravity of the central protostar when formed. Turbulent viscosity is taken into account via the viscous stress tensor Π, the expression for which can be found in Vorobyov & Basu (2010). The expressions for the cooling and heating rates Λ and Γ can be found in Vorobyov et al. (2018). We parametrize the turbulent viscosity in the disk using the α-prescription of Shakura & Sunyaev (1973) with constant α-parameter equal to 0.01, which is consistent with the angular momentum transport rates obtained from the three-dimensional MHD simulations of disks with winds (Suzuki & Inutsuka 2009). Models with lower values of the α-parameter will be considered in a follow-up study.

Dust in our model consists of two components: small (micron-sized) dust and grown dust with a minimum radius of 1.0 μm and a variable upper radius of ar. The dust size distribution has a fixed power law of p = −3.5 in the differential size distribution dn∕daap. All dust in the pre-stellar core is in the form of small dust particles, forming the initial dust mass reservoir. During the core collapse andthe disk evolution, the small dust from the reservoir gradually turns into the grown dust. Small dust particles are assumed to be coupled with the gas, while the dynamics of grown dust is controlled by friction with the gas and by the total gravitational potential of the star, and the gaseous and dusty components. At later stages, when the fragmentation becomes important, some mass of grown dust returns to small dust population. The continuity and momentum equations for small and grown dust are defined as: Σd,smt+p(Σd,smvp)=S(ar),\begin{equation*}\frac{{\partial \Sigma_{\textrm{d,sm}} }}{{\partial t}} + \nabla_{\textrm{p}} \cdot \left( \Sigma_{\textrm{d,sm}} \bm{v}_{\textrm{p}} \right) = - S(a_{\textrm{r}}), \end{equation*}(5) Σd,grt+p(Σd,grup)=S(ar),\begin{equation*}\frac{{\partial \Sigma_{\textrm{d,gr}} }}{{\partial t}} + \nabla_{\textrm{p}} \cdot \left( \Sigma_{\textrm{d,gr}} \, \bm{u}_{\textrm{p}} \right) = S(a_{\textrm{r}}), \end{equation*}(6) t(Σd,grup)+[(Σd,grupup)]p=Σd,grgp+Σd,grfp+S(ar)vp, \begin{align*}\frac{\partial}{\partial t} \left( \Sigma_{\textrm{d,gr}} \, \bm{u}_{\textrm{p}} \right) + [\nabla \cdot \left( \Sigma_{\textrm{d,gr}} \, \bm{u}_{\textrm{p}} \otimes \bm{u}_{\textrm{p}} \right)]_{\textrm{p}} =&\, \Sigma_{\textrm{d,gr}} \, \bm{g}_{\textrm{p}} \nonumber\\ &\ + \Sigma_{\textrm{d,gr}} \bm{f}_{\textrm{p}} + S(a_{\textrm{r}}) \bm{v}_{\textrm{p}}, \end{align*}(7)

where Σd,sm and Σd,gr are the surface densities of small and grown dust; up describes the planar components of the grown dust velocity; S(ar) is the rate of dust growth per unit surface area, the expression for which can be found in Vorobyov et al. (2019); and ar is the maximum radius of grown dust.

The properties of the central protostar (e.g., radius and photospheric luminosity) are calculated using the STELLAR evolution code (Yorke & Bodenheimer 2008; Hosokawa et al. 2013; Vorobyov et al. 2017; Elbakyan et al. 2019). The evolution of the central protostar and the circumstellar disk are connected self-consistently. The stellar mass grows according to the mass accretion rate from the disk. The radiative heating of the disk is also calculated self-consistently in accordance with the protostellar photospheric and accretion luminosities.

The initial surface density and angular momentum distributions in our models have the form of: Σg=r0Σg,0r2+r02,\begin{equation*} \Sigma_{\textrm{g}}=\frac{r_{0}\Sigma_{\textrm{g,0}}}{\sqrt{r^{2}+r_{0}^{2}}},\end{equation*}(8) Ω=2Ω0(r0r)2[1+(rr0)21],\begin{equation*} \Omega=2\Omega_{0}\left(\frac{r_{0}}{r}\right)^{2}\left[\sqrt{1+\left(\frac{r}{r_{0}}\right)^{2}}-1\right],\end{equation*}(9)

where Σg,0 and Ω0 are the angular velocity and gas surface density at the center of the core; r0=Acs2/πGΣg,0$r_{0}=\sqrt{A}c_{\mathrm{s}}^{2}/\pi G\Sigma_{\textrm{g,0}}$ is the radius of the central plateau, where cs is the initial isothermal sound speed in the core. Such a radial profile is typical of pre-stellar cores formed as a result of the slow expulsion of magnetic field due to ambipolar diffusion, with the angular momentum remaining constant during axially-symmetric core compression (Basu 1997). The positive density perturbation A is equal to 1.2, making the core unstable to collapse. We consider three numerical models with different total mass of the initial prestellar core. The model parameters are listed in Table 1. The initial dust-to-gas ratio in all models is 1:100. The initial cores of model L and M have ratio of the rotational to the gravitational energy β = 0.24%, while the core in model S has β = 0.07%. Such values are consistent with the observations of pre-stellar cores (Caselli et al. 2002). We use smaller value of β in model S to simulate more compact and less massive circumstellar disk. The initial gas temperature in collapsing cores is Tinit = 20 K.

Table 1

Model parameters.

3 Results

3.1 Global evolution

In this section, we present the global evolution of circumstellar disks in all three models. Figure 1 presents the gas surface density maps of the inner 1200 × 1200 au2 box for all three models with different initial prestellar core mass at distinct time instances. The computational domain has area 7000 × 7000 au2 (model L) and 3700 × 3700 au2 (model M and model S) including the accreting envelope, but here we show only the inner part with the disk. The time is calculated from the moment of star formation. It is clear from the top row that during its early evolution the disk in model L has a spiral structure and shows fragmentation, while in the disks of models M and S only spiral structure is visible. To check if the disks fulfill the Toomre gravitational instability and fragmentation criterion (Toomre 1964), in the upper-right corner of each panel we present the Toomre Q-parameter for all azimuthal grid points at a specific radial distance from the star. The Q-parameter for the mixture of the dust and gas is defined as: Q=c˜sΩπG(Σg+Σd,sm+Σd,gr),\begin{equation*} Q={\frac{\tilde{c}_{\textrm{s}} \Omega} {\pi G (\Sigma_{\textrm{g}}+\Sigma_{\textrm{d,sm}} + \Sigma_{\textrm{d,gr}})}},\end{equation*}(10)

where c˜s=cs/1+ζd2g$\tilde{c}_{\textrm{s}}=c_{\textrm{s}}/\sqrt{1+\zeta_{\textrm{d2g}}}$ is the modified sound speed (Vorobyov et al. 2018) in the presence of dust and ζd2g is the total dust-to-gas mass ratio.

The Q-parameter in all our models at t = 100 kyr shows values less than the threshold value of disk fragmentation Qfr = 1, meaning that the disks in all our models are characterized by gravitational instability. The spiral structures are clearly seen in all the models during the early disk evolution. The Q-parameter in all our models grows as the disks evolve. The growth of the Q-parameter is due to the decrease of the gas surface density for more than a order of magnitude during initial 900 kyr of disk evolution, while the temperature in the disk during the same time period deceases by only about factor of 3. At t = 250 kyr the Q-parameter in model L becomes greater than the Qfr = 1, while less than the threshold value for the spiral formation Qsp = 1.5. We note that precise Q stability value is still under debate and modern numerical simulations show that the disks become unstable and grow non-axisymmetric disturbances, as multi-armed spirals, for Q < 1.5 (e.g., Durisen et al. 2007; Kratter & Lodato 2016). The minimal value of the Q-parameter in the model L is higher than Qsp at t = 500 kyr, meaning that the disk becomes stable. In contrast, the disks in models M and S become stable much earlier and the Q-parameter in these models is higher than Qsp already at t = 250 kyr. The difference in duration of the gravitationally unstable phase of the disks is caused mainly by the difference in the disk masses of our models: the more massive the disk, the longer it will stay unstable. The maximum disk mass in model L, model M and model S are, respectively, 0.3, 0.2 and 0.1 M. The subsequent evolution of disks in all models after t = 500 kyr shows similar behaviour - the disks are gravitationally stable and nearly axisymmetric.

thumbnail Fig. 1

Gas surface density maps of the inner 1200 × 1200 au2 box at different time moments for all three models. The colorbar is shown in log scale. The insets show the Toomre Q-parameter for all azimuthal grid points at a specific radial distance from the star.

3.2 Dust growth

We further analyze the dust growth and the radial distribution of the maximum grain radius ar in the disk. Due to the aerodynamic drag with the gas, the dust particles lose angular momentum and drift towards the inner parts of the disk, thus experiencing fast radial drift. To grow to the planet embryo size, the dust particles need to overcome this so-called “drift barrier” (Weidenschilling 1977). On the other hand, due to the high relative velocities large dust particles are shuttered rather then stuck together, leading to the so-called “fragmentation barrier” (Blum & Wurm 2008; Brauer et al. 2008; Gonzalez et al. 2017). We do not consider the bouncing barrier in this study (Zsom et al. 2010).

Figure 2 shows the maximum dust radius at a given time and orbital distance for all azimuthal grid points in our models. The color of the dots presents the Stokes number (St) shown in the colorbar in log scale. The Stokes number, or the dimensionless grain stopping time, is a fundamental parameter controlling the dust dynamics and is defined as: St=ΩKρsarρgcs,\begin{equation*} \textrm{St}=\frac{\Omega_{\textrm{K}} \rho_{\textrm{s}} a_{\textrm{r}}}{\rho_{\textrm{g}} c_{\textrm{s}}},\end{equation*}(11)

where ΩK is the Keplerian angular velocity; ρs = 2.24 g cm−3 is the material density of dust grains; cs is the sound speed; ρg=Σg/2πHg$\rho_{\textrm{g}}=\Sigma_{\textrm{g}} / \sqrt{2\pi}H_{\textrm{g}}$ is the gas volume density; Hg is the height scale of the gas.

The thick black line shows the fragmentation barrier defined as: afrag=2Σgufrag23πρsαcs2,\begin{equation*}a_{\textrm{frag}} = \frac{2\Sigma_{\textrm{g}}u^2_{\textrm{frag}}}{3\pi \rho_{\textrm{s}} \alpha c_{\textrm{s}}^2}, \end{equation*}(12)

where ufrag = 30 m s−1 is a threshold value for the dust fragmentation velocity (e.g., Wada et al. 2009). Here we plot the azimuthally averaged value of fragmentation barrier. Given the uncertainties in the α-parameter and the fragmentation velocity ufrag, we discuss the expected outcomes from the changes of these model parameters in Appendix A.

Evidently, during less than 100 kyr of initial disk evolution, the main part of the initial micron sized dust in the disks of all our models is efficiently converted into the grown dust particles, reaching few cm in size at the radial distances ≲ 50 au. The growth of a dust particle ceases when its radius reaches the fragmentation limit defined with Eq. (12). The size of grown dust particles reaches the fragmentation barrier for the radial distances ≲ 20 au, while for the outer parts of the disk it never reaches the fragmentation barrier. As has been shown by Tsukamoto et al. (2017), the orbital radius for planetesimal formation by coagulation of fluffy dust particles is ≈ 20 au for a gravitationally unstable disk around a solar mass and the dust growth is regulated by the radial drift barrier on radial distance r ≳ 20 au. When dust is transported from regionswith larger afrag to regions with smaller afrag, we convert some of grown dust mass to small dust to simulate fragmentation. This means that dust growth and dynamics is regulated bythe fragmentation barrier in the inner parts of the disk, while in the outer parts of the disk it is regulated by the drift barrier (e.g., Birnstiel et al. 2016). During the subsequent disk evolution, the fragmentation barrier for all our models shows a well-defined peak at 20 au, decreasing towards the star. This means that the maximum dust radius in the disk is reached at several tens of au.

During the early evolution of the disk, the Stokes number of dust particles in the outer regions of the disks in our models shows high variability, reaching 0.3. Such a high Stokes numbers in the outer part of the disk are reached because the gas density is sufficiently low and the mean free path of gas particles is much higher than the size of dust particles. However, as the disk evolves, the Stokes number in the outer parts of the disk decreases. At the same time the maximum values (St ≈ 0.2) of the Stokes number are reached in the disk at few tens of au, where the dust-size radial distribution shows a peak. The Stokes number of the dust particles, similarly to the size of the particles, is decreasing towards the star for the radial distance of r < 20 au.

The envelope and outer parts of the disk, due to the long evolution time scale at large radial distances, play a role of feeding zone for the inner disk, providing not only gas but also small micron-sized dust particles, which grow and migrate inwards, thus compensating the rapid accretion of grown mm- to cm-sized particles onto the central star (Kornet et al. 2004; Garaud 2007). Grown dust particles with mm-cm sizes and relatively high Stokes numbers, also known as pebbles, that are decoupled from the gas, play an important role in formation of planetesimals and planetary cores. Johansen & Lacerda (2010) showed that protoplanets can accrete half of their mass from pebbles, thereby giving rise to a planet formation model called pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012; Ida et al. 2016; Johansen & Lambrechts 2017). Pebble accretion can reduce the planetary core growth timescales below the typical lifetime of the circumstellar disks even at 100 au orbital distances. Bitsch et al. (2015) showed that gas giants formed by pebble accretion are born between 10 and 50 au, and then migrate inward. Although the term “pebble” is widely used, it has no specific definition. Here we define pebbles as dust particles that have the radius greater than apeb = 0.5 mm. We are not using the Stokes number in the pebble definition, since it implies that the pebble size in the definition will change depending on the radial distance from the star, which in turn will bring a degeneracy in the pebble definition.

thumbnail Fig. 2

Radial distribution of the maximum radius of grown dust for all azimuthal grid points shown for consecutive times for all models. Color of dots shows the value of Stokes number for each azimuthal grid point. The colorbar is in log scale. Solid black line marks the dust fragmentation size afrag.

thumbnail Fig. 3

Azimuthally averaged absolute values of gas (blue curve), grown dust (red curve), and pebble (green curve) accretion (transport) rates versus radial distance from the star at distinct evolutionary times for all our models. Dashed part of the curve shows the outward migration, while the solid part – the inward migration.

3.3 Dust dynamics

In this section we study the dynamics of gas and dust in the disk. Much attention is paid to the dynamics of pebbles in the disk, which are essential for the planetesimal formation.

The mass flux in disk at the radial distance r is calculated as: M˙(r)=i=1n(Mi/Si)urirΔϕ,\begin{equation*}\dot{M}(r) = \sum_{i=1}^{n} (M^i/S^i) \bm{u}_{\textrm{r}}^i r \Delta\phi, \end{equation*}(13)

where Mi is the mass of matter confined in the ith computational cell with surface Si, uri$\bm{u}_{\textrm{r}}^i$ is the radial velocity of the matter in the ith computational cell at the radial distance r, and Δϕ is the spatial grid step in azimuthal direction. Summation is done for the all n azimuthal cells at each radial distance r.

Having the minimal radius of pebbles in each cell – apeb, we calculate the mass of pebbles Mpeb inside each computational cell as: Mpeb=apebamaxm (a)a3.5da=Mgr.dust(amaxapeb)amaxamin,\begin{equation*}M_{\textrm{peb}} = \int_{a_{\textrm{peb}}}^{a_{\textrm{max}}} m(a)a^{-3.5}\textrm{d}a = \frac{M_{\textrm{gr.dust}} (\sqrt{a_{\textrm{max}}} - \sqrt{a_{\textrm{peb}}})}{\sqrt{a_{\textrm{max}}} - \sqrt{a_{\textrm{min}}}}, \end{equation*}(14)

where amax and amin are, respectively, the maximum and the minimum radius of dust particles in the cell, and Mtot. dust is the total mass of dust in the cell.

It is interesting to see how the gas and dust are accreted on the disk from the envelope and later transported in the disk towards the central star. Figure 3 shows the absolute values of gas (blue curve), grown dust (red curve), and pebbles (green curve) mass accretion (transport) rates vs. radial distance from the central star for different time moments. Dashed part of the curves shows the outward migration, while the solid part shows the inward migration.

During the embedded phase, at t = 100 kyr, gas accretion from the envelope to the edge of the disk has a quite smooth character for all models, changing from g = few × 10−8 M yr−1 at the outer edge of the envelope to g = few × 10−6 M yr−1 at the outer edge of the disk. The accretion of grown dust inside the envelope also shows smooth accretion with accretion rates from g = few × 10−12 M yr−1 to g = few × 10−10 M yr−1. In the inner, r ≲ 20 au, part of the disk, the gas and the dust during the embedded phase show less variable inward migration, while in the outer parts of the disk (r ≳ 20 au) show highly variable transport rates not only for inward, but also outward migration. The transport rates in the outer region are changing by more than a order of magnitude. Such a behaviour is a result of gravitational instability in the disk during its initial embedded evolution. The transport rates of pebbles in the inner, r ≲ 100 au, part of the embedded disk are similar to the transport rates of grown dust particles, while having slightly lower values. Transport rates of pebbles in the disk show variations by more than two orders of magnitude. Such a behavior is a result of vigorous gravitational instability taking place in the young embedded disk. The temperature in the inner part of the disk is quite high during the early embedded phase, reaching 150 K (approximately the midplane temperature at the water ice line) at the radial distance ≈20 au, while during the later disk evolution 150 K is reached only at the radial distance ≈4 au. Similar results are found by Homma & Nakamoto (2018) showing that the ice line reaches about 12 au at 0.38 Myr and migrates to 3 au at 1 Myr.

At t = 500 kyr, when the embedded phase is already over, mass fluxes for both the gas and the grown dust at radial distances r ≳ 100 au become mainly positive, thus showing outward migration. Only the matter from very outer part of computational domain in model L shows inward migration. This is due to the still in-falling very outer parts of the relatively massive initial core. For the inner r ≲ 100 au of the disk both the gas and the grown dust show relatively smooth inward migration at about 10−7 Myr−1 and 10−9 Myr−1 rates, respectively. The pebbles also show smooth inward mass transport for r ≲ 100 au radial distances having transport rates slightly lower than the transport rates of grown dust particles, which are consistent with expected pebble fluxes in the disk (Lambrechts & Johansen 2014; Lambrechts et al. 2019). We note that the inner r ≲ 2 au part of the disk shows some variability connected with the boundary effects between the central smart cell and the active disk. At t = 900 kyr the gas and the grown dust in the outer r ≳ 100 au part of the disk show only outward migration because of the viscous spreading of the disk. The mass fluxes of gas, grown dust and pebbles at the inner r ≲ 100 au of the disk show qualitatively similar behaviour as at t = 500 kyr, but with slightly lower mass transport rates. Note that the pebbles are not forming in the disk on the radial distances r ≳ 100 au, which is consistent with the radial distance of pebble formation line obtained by Lambrechts & Johansen (2014).

The pebble mass fluxes in the disk, peb, plays an essential role in the formation of rocky Earth-like planets, hot super-Earths and gas giant planets (Lambrechts et al. 2019; Izidoro et al. 2019; Bitsch et al. 2019; Johansen et al. 2019; Lenz et al. 2019). An increase in total pebble mass flux in the disk by only a factor of 2 can cause the formation of super-Earth planets instead of Earth-like planets (Lambrechts et al. 2019). Thus, it is important to understand how the pebbles are transported in the disk. In Fig. 4 we show the pebble mass flux dependence on the gas mass flux. Color of the dots shows the evolutionary time in Myr. The pebble flux shows a tendency to decrease as the disk evolves. Clearly, the pebble flux in all our models shows a power dependence on gas flux. The dashed line in each panel shows the first order polynomial curve that is the best fit for each model. The gas vs. peb relation can be presented with a polynomial curve as: M˙peb[Myr1]=bM˙gasa[Myr1].\begin{equation*}\dot{M}_{\textrm{peb}}[M_{\odot} \textrm{yr}^{-1}] = {b} \dot{M}_{\textrm{gas}}^{a}[M_{\odot} \textrm{yr}^{-1}]. \end{equation*}(15)

The polynomial a and b coefficients for all models are presented in Table 2. The dotted lines show the ± 3σ deviation from the best-fit values. All the data points that lay between the two red dotted lines are in prediction interval with 99.87% probability. We overplot the peb = 0.01g dependence with the dash-dotted line for the reference.

In addition to the transport rates of matter in the disk, it is interesting to study the time dependence of accretion rates of matter from the disk onto the central star. The mass accretion rates could be estimated from the accretion luminosity and the stellar mass and radius (see Eq. (8) in Vorobyov et al. 2017). The accretion rates are found to correlate with the stellar mass (Calvet et al. 2004; Herczeg & Hillenbrand 2008; Manara et al. 2015) and stellar age (Hartmann et al. 1998; Manara et al. 2012; Venuti et al. 2014). In Fig. 5, we present the time dependence of mass accretion rates of gas M˙gas*$\dot{M}_{\textrm{gas}}^{*}$ and pebbles M˙peb*$\dot{M}_{\textrm{peb}}^{*}$ onto the star. The color of dots presents the age of the system in Myr. Evidently, a strong correlation between the gas and pebble mass fluxes exist at the early embedded phase of disk evolution, while during the later disk evolution the correlation becomes weaker. The reason for the weak correlation during the late disk evolution is the fact that as the dust particles grow to the pebble sizes and their Stokes number increases, they decouple from the gas. As a result, the pebble flux in the disk starts to decrease slower than the gas flux and becomes less correlated. We show the best-fit first order polynomial curve for each model with the dashed lines. The polynomial a and b coefficients for each best-fit curve are shown in Table 3.

Gas surface density distribution in circumstellar disks plays an important role in the process of planet formation. Figure 6 shows the radial distribution of azimuthally averaged dust-to-gas mass ratio (ζd2g), surface density of gas, grown dust, and pebbles for different time moments. The dust-to-gas ratio stays close to the canonical 1:100 value in the entire disk during the disk evolution in all our models, showing increase only in the inner 1 au near the central smart cell. This increase is a result of fast mass transport through the central smart cell leading to the gas depletion near it (Vorobyov et al. 2019). The radial profile of gas surface density during the early disk evolution is characterized by a slowly decreasing density in the central region (≲ 30 au) and a steeply declining tail at larger radii. On the other hand, during the early disk evolution the radial distribution of grown dust and the pebbles have no inner slow-decreasing region, showing steep power law distribution in the inner part of the disk and more steeply decreasing tail at the larger radii. As time progresses, the gas surface density distribution becomes smooth and less steep for the larger radii, while the grown dust and pebble distributions begin to obey a single power-law for the entire disk. The radial profile of pebbles shows r-dependency of Σpebr−1.2 in model L and M. In model S, the surface density of pebbles becomes proportional to r−1.5 after the embedded phase of evolution. Such a dependency is close to the MMSN dependency, which yield Σpebr−1.5. Similar results are obtained by Birnstiel et al. (2012) for the dust surface density profile in fragmentation limited regime with relatively large grains present. Our results also agree with the pebble surface densities obtained for constant mass flux and Stokes number for the pebbles (Lambrechts et al. 2019).

thumbnail Fig. 4

Relation between pebble mass flux peb and gas mass flux gas for all our models. Color of the dots presents the age of the system in Myr. The black dashed line shows the best-fit curve for each model. The red dotted lines show the ± 3σ deviation from the best-fit values. The dash-dotted line shows the peb = 0.01g dependence.

Table 2

a and b parameters for different threshold values of minimum dust size in the pebble definition.

thumbnail Fig. 5

Similar to Fig. 4, but for the mass accretion rates of gas and pebbles onto the central star.

Table 3

Similar to Table 2, but for the mass accretion rates of gas and pebbles onto the central star.

thumbnail Fig. 6

Azimuthally averaged dust-to-gas mass ratio (black), surface density of gas (blue curve), grown dust (red curve), and pebbles (green curve) versus radial distance from the star at distinct evolutionary times for all our models.

3.4 Comparison with observations

Here we compare surface densities obtained in our simulations with the surface densities of relatively massive disk in AS 209, HD 163296 and DoAr 25 systems. While neither of the models considered in this paper were specifically tuned for these systems, it is still interesting to check how our models perform against observationally derived disk parameters. Due to the angular momentum conservation and viscous spreading the disks become larger in time and less dense. Small dust particles that are coupled with the gas will follow the viscous evolution of disk and migrate outwards. Such a distribution of dust in the outer disk complicates the definition of outer edge of the disk, and disk physical and visible sizes for both dust and gas components can differ by a factor of a few (Facchini et al. 2017; Rosotti et al. 2019; Trapman et al. 2019). One of the main disk parameters used in planet formation theory is the disk surface density, which is converted into the total disk mass. However, different methods of disk mass determination are often inconsistent and can vary greatly (Bergin et al. 2013). Recently, Powell et al. (2019) applied the alternative method for disk mass determination suggested by Powell et al. (2017) to the multiwavelength observations of few young disks. The method determines the disk surface density without usage of tracer-to-H2 ratio or an assumed dust opacity model. The disk outer edge is inferred from the millimeter wavelength observations and then related to the maximum radial location of different particle sizes in the disk. The surface density profile of gaseous disk is determined based on the aerodynamic properties of the dust particles. The derived surface density profile is then used as benchmark to scale previously modeled surface density profiles derived from combined multiwavelength dust or CO emission observations.

In Fig. 7, we present the azimuthally averaged surface density distributions of gas, grown and small dust at t = 946 kyr for all our models. The radial distributions show qualitatively same behavior for all our models. The gas surface density slowly decline in the inner few tens of au region and smoothly becomes more steep at larger radii. The radial distribution of small dust surface density for all models follows the shape of gas surface density distribution, except the very inner 2 au of the disk. Fast mass transport through the central smart cell leads to gas and small dust depletion in the inner 2 au of the disk (Vorobyov et al. 2019). Thus, the inner 2 au of the disks is mainly populated by grown dust particles. The radial distribution of grown dust particles shows different behavior, having a close to a power-law distribution for the entire disk. On the radial distance of about 1000 au the radial distribution of grown dust shows steep decrease in all our models. All the grown dust particles on such high radial distance has drifted inwards. In the top panel of Fig. 7, we overplot the disk surface densities of three1 young stellar objects derived by Powell et al. (2019). The ages and disk masses for the observed objects, along with the disk masses of our models at t = 946 kyr are presented in Table 4. It is assumed that the disk mass is equal to the total mass of gas and dust confined inside the radial distances at which the radial distribution of grown dust shows steep decrease. We compare the observed objects only with our model L because model M and S have much lower masses compared to the disk masses of observed objects. The surface densities of observed objects show quite good agreement with gas surface density of model L for radial distances 10 ≲ r ≲ 300 au. For smaller radial distances (r ≲ 10 au), disk surface densities for DoAr 25 and HD 163296 (CO) show quite good agreement, while for AS 209 and HD 163296 (dust), they overestimate our model results.

It is important to note that dust surface density distribution is not necessarily mirrored by its thermal emission. The radiation intensity of an isothermal slab of dust with temperature T and absorption coefficient κν can be calculated as: Iν=Bν(T)(1exp(κνΣd)).\begin{equation*} I_{\nu}=B_{\nu}(T)(1-\exp\left(-\kappa_{\nu}\Sigma_{\textrm{d}}\right)). \end{equation*}(16)

Here we neglect the scattering, to role of which is actively discussed at the moment (Zhu et al. 2019; Liu 2019; Carrasco-González et al. 2019). To calculate the absorption coefficient κν we use Mie theory for a mix of compact carbonaceous and silicate dust as in Pavlyuchenkov et al. (2019), assuming their power-law distribution as in our dust evolution model (Vorobyov et al. 2018). Dust opacity κν experiences a distinct maximum near amaxλ∕2π (Birnstiel et al. 2018) that is reflected as a knee in the intensity radial distribution at the location where grains attain this characteristic size.

Vertical dashed lines in each panel of Fig. 7 mark the radial distance beyond which all the dust particles in the model have radii of ar < 0.85∕2π mm, meaning that these particles will not effectively contribute to the continuum emission at wavelengths λ > 0.85 mm. The vertical lines for all models lay on the radial distance 168 ± 10 au from the central star, which is consistent with the observational estimates of visible sizes of protoplanetary disks in dust continuum (Ansdell et al. 2018).

Figure 8 presents the radial distributions of model L synthetic intensities in ALMA Band 3 (2.80 mm), Band 6 (1.30 mm), and Band 7 (0.85 mm). Images are non-convolved, but we present intensities in units per 34 mas beam. Stars mark the radial distance at which the grown grain size is equal to λ∕2π and the distribution of intensity becomes steeper. We define these radial distances as the outer dust radius for each band, as usually done for the observational data. Figure 8 shows that the disk emission outside the outer radius traces the density distribution well as the disk periphery is relatively isothermal due to the external heating by interstellar radiation field and the dust opacity coefficient, κν, also does not vary due to the small size of the grains. The inner disk emission increases much faster towards the star than dust surface density as temperature also increases inwards. We compare the outer dust radii obtained from reproduced intensities of model L with the ones for AS 209, DoAr 25 and HD 163296, taken from Powell et al. (2019). The results are shown in Table 5. The outer dust radius of model L is in quite good agreement with the outer dust radius of DoAr 25, while the ones for As 209 and HD 163296 are overestimated.

thumbnail Fig. 7

Radial dependence of the azimuthally averaged surface densities of gas (blue curve), grown dust (red curve), and small dust (green curve) at evolutionary time t = 946 kyr for all ourmodels. The dashed vertical lines mark the radial distance at which the radius of dust particles ar = 0.85∕2π mm. Thin curves in top panel show the gas surface densities of three observed disks derived by Powell et al. (2019).

Table 4

Disk masses of the observed objects obtained from Powell et al. (2019) along with the disk masses of our models.

thumbnail Fig. 8

Synthetic intensities in ALMA Band 3 (2.80 mm), Band 6 (1.30 mm), and Band 7 (0.85 mm) and total dust surface density for model L at t = 964 kyr. The stars mark the radial distance at which dust size becomes equal to λ∕2π.

Table 5

Outer dust radii of the three observed disks for different ALMA bands obtained from Powell et al. (2019) and outer dust radius of our model L obtained from the reproduced intensities.

4 Discussions

4.1 Parameter space study

In order to check how our results for pebble fluxes depend on the threshold value of minimum dust particle radius in pebble definition (apeb), we perform a parameter space analysis varying the value of apeb. In Table 2 we present the a and b polynomial coefficients for three different apeb values in pebble definition that are used in Eq. (15) for all our models. The last row in the table shows the averaged values for all the models. Clearly, the power-law index for the models with different apeb changes slightly, meaning that the correlation between transport rates of gas and pebbles is weakly dependent on the minimum size of pebbles in range of 0.5−2 mm.

We also perform parameter space study to check how the mass accretion rates of gas and pebbles onto the star depend on the apeb. The a and b polynomial coefficients of best-fit curves for different apeb values are shown in Table 3. Unlike the transport rates in the disk, the polynomial coefficients show slight increase for mass accretion rates onto the star. This result is expected because the dust particles are growing during their inward migration in the disk, thus more grown particles with larger sizes are accreted onto the star.

4.2 Model caveats

In the current model, we consider only two dust populations of small and grown dust grains with the assumption of a power-law size distribution, which is yet rather simplistic. Multi-grain size simulations of dust growth in hydrodynamically evolving disks have recently become available (Drążkowska et al. 2019), however full-disk simulations with good spatial resolution are very challenging because of high CPU costs. Snow lines play important role in the evolution of the dust particles in protoplanetary disks (e.g., Gárate et al. 2020; Musiolik & Wurm 2019; Ros et al. 2019; Vericel & Gonzalez 2020; Ziampras et al. 2020). In future studies, we plan to improve our model by introducing H2 O, CO, and CO2 snow lines and adopting dust fragmentation velocity (ufrag) depending on the position of dust particle relative to the snow line. In addition, other important factors such as bouncing barrier (Zsom et al. 2010) and grain charging (Okuzumi et al. 2011; Akimkin 2015; Steinpilz et al. 2019) are planned to be added to the model. It was shown by Vorobyov & Elbakyan (2018) that the numerical resolution becomes important in the study of the migration of dense gaseous clumps. However, we expect that the change in numerical resolution will not dramatically influence the global disk evolution in our model. The model reproduces similar results for gas and dust components in a test problem with a numerical resolutions varying by factor of 2 (Vorobyov et al. 2018).

4.3 Planet formation pathways

The classical planet formation models assume that the dust particles grow up to become kilometer-sized planetesimals, which, in turn, grow into protoplanets colliding and sticking with each other (Safronov 1972). Recently, as an alternative to the classical model, the promising streaming instability (Youdin & Goodman 2005; Johansen et al. 2015; Carrera et al. 2017) and the pebble accretion scenarios have been proposed (Lambrechts & Johansen 2012; Johansen & Lambrechts 2017). Tanaka & Tsukamoto (2019) showed that the pebble accretion could possibly lead to the early planet formation in the Class 0/I disks. Our models showed that only inner r ≲ 2 au part of the disk shows dust-to-gas ratios higher than the canonical 1:100 value, which are needed for the streaming instability to take place. Thus, the combination of streaming instability and the pebble accretion could possibly form protoplanetary embryos in the inner, r ≲ 2 au, part of the disk. We also predict that the dust-to-gas ratio in the models with lower α-parameter will become higher than the canonical value for the radial distances up to 10 au, thus leading to the protoplanetary embryos formation not only in the inner few au of the disk. In our models with α = 0.01 the viscous torques are strong in the entire disk, while in the models with the lower α-parameter the viscous torques are weak and the matter is transported by the gravitational torques. However, both gravitational and viscous torques are weak in the inner warm region of the disk with low α-parameter, leading to the formation of pressure maximum and accumulation of the matter.

5 Conclusions

In this paper, we study the long-term evolution of accreting circumstellar disks with their surrounding envelopes using the numerical hydrodynamics code FEOSAD. The joint dynamics of gas and dust (including dust growth) is calculated for about 950 kyr of disk evolution. Three models with different disk masses are considered. The main focus of the study is the evolution and dynamics of grown dust particles, in particular the pebbles.

Our main results can be summarized in the following:

  • The dust in the disk grows to the sizes of few cm in the inner r ≲ 100 au during less than 100 kyr after the disk formation. Because of a constant supply of small dust from the envelope, the grown few cm sized particles exist in the disk on the radial distances r ≲ 100 au for more than 900 kyr.

  • The main part of the grown dust particles in the inner r ≲ 100 au of the disk have pebble sizes. A strong correlation exists between the gas and pebble fluxes in the disks, showing almost universal power-law dependency for disks with different masses. The power-law dependency is more steep for the early embedded phase of disk evolution, while becomes more declivous as the disk evolves.

  • The radial distribution of pebbles in the disks shows almost universal power-law distribution, close to the MMSN distribution slope of Σ ∝ r−1.5.

  • The gas surface density of model L shows quite good agreement with the gas surface densities of observed disks obtained from Powell et al. (2019). The outer dust radius of model L obtained from the reproduced intensities shows a good agreement with the outer dust radius of DoAr 25 obtained from Powell et al. (2019).

Acknowledgements

We thank the anonymous referee for a insightful report, which helped to improve this paper. Research was financially supported by the Ministry of Science and Higher Education of the Russian Federation (State assignment in the field of scientific activity, Southern Federal University, 2020). V.G.E. acknowledges the Swedish Institute for a visitor grant allowing to conduct research at Lund University. A.J. was supported by the Swedish Research Council (grant 2018-04867), the Knut and Alice Wallenberg Foundation (grant 2012.0150) and the European Research Council (ERC Consolidator Grant 724687-PLANETESYS). M.L. was supported by the Knut and Alice Wallenberg Foundation (grant 2012.0150). V.A. was supported by RFBR grant 18-52-52006.

Appendix A The effect of fragmentation velocity and midplane turbulence on the dust evolution

In this appendix, we study how the change in fragmentation velocity ufrag and the α-parameter in the definition of fragmentation barrier (Eq. (12)) affects the evolution and the dynamics of dust particles in model L. We define the modified model L as model L2.

The threshold values of fragmentation velocity ufrag for different dust particle compositions and sizes have been heavily studied by number of authors either experimentally or numerically (Blum & Wurm 2008; Teiser & Wurm 2009; Wada et al. 2009, 2013; Zsom et al. 2010; Meru et al. 2013; Yamamoto et al. 2014; Gundlach & Blum 2015; Bukhari Syed et al. 2017), giving a wide range of possible values for the velocities. This range varies between 1 and 30 m s−1 for the silicate particles, and between10 and 80 m s−1 for icy-grains (Charnoz et al. 2019). In model L2 we set ufrag = 1 m s−1.

In our models, the same α-parameter is used for the description of disk viscosity and the strength of midplane turbulence that affects the fragmentation of dust particles. In model L2, we consider the viscous parameter αv = 10−2 and the turbulence strength αt = 10−4. The choice of αt < αv if justified by the models of protoplanetary disks where the midplane is found to be less active (e.g., Turner et al. 2014). Okuzumi & Hirose (2011) found that αt is for about an order of magnitude lower than the αv, while based on the observations of the disk around HL Tau, Pinte et al. (2016) found that αt is of the order of a few 10−4. In model L2 in Eq. (12) we set α = αt = 10−4, which is also consistent with the values of αt used in other protoplanetary disk models (e.g., Carrera et al. 2017; Drążkowska & Dullemond 2018).

The left column in Fig. A.1 shows the radial distribution of maximum grown dust radius for all azimuthal grid points at distinct time moments. The Stokes number is shown with the color. Solid black line marks the dust fragmentation size. Unlike model L, the maximum radius of dust particles in model L2 reaches only few mm. However, the evolution of dust in model L2 is qualitatively similar to the one in model L. The dust evolution is regulated with the drift for the radial distance r ≳ 20 au, while with the fragmentation for the radial distance r ≲ 20 au. As a result, a well-defined peak in the maximum dust radius is developed at ≈ 20 au. The fragmentation radius decreases closer to the star, reaching few × 10−2 cm. After the end of embedded phase, the fragmentation radius at r = 1 au becomes lower than 0.5 mm (the minimum pebble size), meaning that no pebbles (dust particles greater than 0.5 mm in radius) are accreted onto the star. As a result of lower maximum dust radii, the Stokes number of dust particles in model L2 is about an order of magnitude lower than its counterpart in model L. The Stokes number of dust particles in the inner 10 au of the disk is of order 10−3, meaning that the dust particles are well coupled with the gas.

thumbnail Fig. A.1

Left column: radial distribution of the maximum radius of grown dust for all azimuthal grid points in model L2. Color of dots shows the value of Stokes number for each azimuthal grid point. Solid black line shows the dust fragmentation size afrag. Middle column: azimuthally averaged dust-to-gas mass ratio (black), surface density of gas (blue curve), grown dust (red curve),and pebbles (green curve) vs. radial distance from the star. The dashed line shows the power-law slope ∝ r−1.5. Right column: azimuthally averaged absolute values of gas (blue curve), grown dust (red curve), and pebble (green curve) accretion (transport) rates vs. radial distance from the star. Dashed part of the curve shows the outward migration, solid part – the inward migration.

The middle column in Fig. A.1 shows the radial distribution of azimuthally averaged dust-to-gas mass ratio (black), surface density of gas (blue curve), grown dust (red curve), and pebbles (green curve) at distinct evolutionary times. The power-law slope ∝ r−1.5 is shown with the dashed line. The evolution of gas surface density in model L2 is similar to the one in model L, while the evolution of grown dust and pebble surface density differs. It is easy to notice that unlike model L, in model L2 the grown dust and pebble surface density for the radial distance r ≲ 20 au does not follow the r−1.5 slope and becomes more shallow. This is a result of lower fragmentation radius for the inner part of the disk of model L2 compared to the model L. We note that the surface density of pebbles sharply drops at radial distance r ≈ 1.1 au because at this distance the fragmentation radius becomes lower than the minimum radius of pebbles. The azimuthally averaged dust-to-gas ratio (ζd2g) in model L2 stays equal to the initial 10−2 value for the entire disk at all evolutionary time moments.

The right column in Fig. A.1 shows the radial distribution of azimuthally averaged absolute values of gas (blue curve), grown dust (red curve), and pebble (green curve) accretion (transport) rates at distinct evolutionarytimes. Dashed line shows the outward migration, while the solid line shows the inward migration. At the early embedded phase the transport rates of gas, grown dust, and pebbles in model L2 similar to model L show high variability with inward and outward migration. The values of transport rates in model L2 are of the same order as in model L, differing only by factor of few. At the evolutionary times after the end of embedded phase, transport rates of gas, grown dust and pebbles show relatively smooth inward migration at the radial distance r ≲ 100 au with values similar to the ones in model L. In contrast to model L, the transport rate of pebbles rapidly decreases at r ≈ 1.1 au.

We also check how the pebble and gas mass fluxes in the disk are affected by the change of the parameters in the fragmentation radius definition. The top panel of Fig. A.2 presents the dependence of the pebble mass flux peb on the gas mass flux gas in the disk. The color of the dots represent the age of the system. Similarly to model L, there is a strong correlation between the pebble and gas fluxes in model L2. The polynomial coefficients a and b of the best-fit curve for gas vs. peb relation are presented in the first row of Table A.1. The polynomial coefficients in model L2 are close to the ones in model L, but slightly higher. The bottom panel of Fig. A.2 shows the time dependent relation between the mass accretion rates of gas M˙gas*$\dot{M}_{\textrm{gas}}^{*}$ and pebbles M˙peb*$\dot{M}_{\textrm{peb}}^{*}$ onto the central star in model L2. We note that the lower values of M˙peb*$\dot{M}_{\textrm{peb}}^*$ are missing due to the fact that the pebbles are absent at the radial distance r ≲ 1.1 au after the end of embedded phase. The polynomial coefficients a and b of the best-fit curve for M˙gas*$\dot{M}_{\textrm{gas}}^{*}$ vs. M˙peb*$\dot{M}_{\textrm{peb}}^{*}$ relation are presented in the second row of Table A.1. The coefficients in model L2 are slightly higher than in model L forthe pebbles sizes of 0.5 and 1 mm, while slightly lower for the pebble size of 2 mm.

thumbnail Fig. A.2

Relation between pebble mass flux peb and gas mass flux gas in the entiredisk (top panel) and onto the central star (bottom panel) for model L2. Color of dots presents the age of the system in Myr. The black dashed line shows the best-fit curve for the data in each panel. The red dotted lines show the ± 3σ deviation from the best-fit values. The dash-dotted line shows the M˙peb*=0.01M˙g*$\dot{M}^*_{\textrm{peb}}\,{=}\,0.01\dot{M}^*_{\textrm{g}}$ dependence.

Table A.1

Parameters a and b for different threshold values of minimum dust size in the pebble definition.

References

  1. Akimkin, V. V. 2015, Astron. Rep., 59, 747 [NASA ADS] [CrossRef] [Google Scholar]
  2. Andrews, S. M., & Williams, J. P. 2007, ApJ, 671, 1800 [CrossRef] [Google Scholar]
  3. Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21 [NASA ADS] [CrossRef] [Google Scholar]
  4. Basu, S. 1997, ApJ, 485, 240 [CrossRef] [Google Scholar]
  5. Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Nature, 493, 644 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  6. Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Birnstiel, T., Fang, M., & Johansen, A. 2016, Space Sci. Rev., 205, 41 [NASA ADS] [CrossRef] [Google Scholar]
  8. Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bitsch, B., Izidoro, A., Johansen, A., et al. 2019, A&A, 623, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Boley, A. C., Hayfield, T., Mayer, L., & Durisen, R. H. 2010, Icarus, 207, 509 [NASA ADS] [CrossRef] [Google Scholar]
  13. Boss, A. P. 1998, ApJ, 503, 923 [NASA ADS] [CrossRef] [Google Scholar]
  14. Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Bukhari Syed, M., Blum, J., Wahlberg Jansson, K., & Johansen, A. 2017, ApJ, 834, 145 [NASA ADS] [CrossRef] [Google Scholar]
  16. Calvet, N., Muzerolle, J., Briceño, C., et al. 2004, AJ, 128, 1294 [NASA ADS] [CrossRef] [Google Scholar]
  17. Carrasco-González, C., Sierra, A., Flock, M., et al. 2019, ApJ, 883, 71 [NASA ADS] [CrossRef] [Google Scholar]
  18. Carrera, D., Gorti, U., Johansen, A., & Davies, M. B. 2017, ApJ, 839, 16 [NASA ADS] [CrossRef] [Google Scholar]
  19. Caselli, P., Benson, P. J., Myers, P. C., & Tafalla, M. 2002, ApJ, 572, 238 [NASA ADS] [CrossRef] [Google Scholar]
  20. Charnoz, S., Aléon, J., Chaumard, N., Baillié, K., & Taillifet, E. 2015, Icarus, 252, 440 [NASA ADS] [CrossRef] [Google Scholar]
  21. Charnoz, S., Pignatale, F. C., Hyodo, R., et al. 2019, A&A, 627, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Drążkowska, J., & Dullemond, C. P. 2018, A&A, 614, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Drążkowska, J., Windmark, F., & Dullemond, C. P. 2013, A&A, 556, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Drążkowska, J., Li, S., Birnstiel, T., Stammler, S. M., & Li, H. 2019, ApJ, 885, 91 [NASA ADS] [CrossRef] [Google Scholar]
  25. Dullemond, C. P., Birnstiel, T., Huang, J., et al. 2018, ApJ, 869, L46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Durisen, R. H., Boss, A. P., Mayer, L., et al. 2007, Protostars and Planets V (Tucson, AZ: University of Arizona Press), 607 [Google Scholar]
  27. Elbakyan, V. G., Vorobyov, E. I., Rab, C., et al. 2019, MNRAS, 484, 146 [NASA ADS] [CrossRef] [Google Scholar]
  28. Facchini, S., Birnstiel, T., Bruderer, S., & van Dishoeck, E. F. 2017, A&A, 605, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Gárate, M., Birnstiel, T., Drazkowska, J., & Stammler, S. M. 2020, A&A, 635, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Garaud, P. 2007, ApJ, 671, 2091 [NASA ADS] [CrossRef] [Google Scholar]
  31. Gonzalez, J. F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984 [NASA ADS] [Google Scholar]
  32. Gundlach, B., & Blum, J. 2015, Icarus, 257, 126 [NASA ADS] [CrossRef] [Google Scholar]
  33. Haghighipour, N., & Boss, A. P. 2003, ApJ, 583, 996 [NASA ADS] [CrossRef] [Google Scholar]
  34. Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385 [NASA ADS] [CrossRef] [Google Scholar]
  35. Herczeg, G. J., & Hillenbrand, L. A. 2008, ApJ, 681, 594 [NASA ADS] [CrossRef] [Google Scholar]
  36. Homma, K., & Nakamoto, T. 2018, ApJ, 868, 118 [NASA ADS] [CrossRef] [Google Scholar]
  37. Hosokawa, T., Yorke, H. W., Inayoshi, K., Omukai, K., & Yoshida, N. 2013, ApJ, 778, 178 [NASA ADS] [CrossRef] [Google Scholar]
  38. Huang, J., Andrews, S. M., Pérez, L. M., et al. 2018, ApJ, 869, L43 [NASA ADS] [CrossRef] [Google Scholar]
  39. Ida, S., Guillot, T., & Morbidelli, A. 2016, A&A, 591, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Izidoro, A., Bitsch, B., Raymond, S. N., et al. 2019, A&A, submitted, [arXiv:1902.08772] [Google Scholar]
  41. Johansen, A., & Lacerda, P. 2010, MNRAS, 404, 475 [NASA ADS] [Google Scholar]
  42. Johansen, A., & Lambrechts, M. 2017, Ann. Rev. Earth Planet. Sci., 45, 359 [NASA ADS] [CrossRef] [Google Scholar]
  43. Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022 [Google Scholar]
  44. Johansen, A., Blum, J., Tanaka, H., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson, AZ: University of Arizona Press), 547 [Google Scholar]
  45. Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, Sci. Adv., 1, 1500109 [NASA ADS] [CrossRef] [Google Scholar]
  46. Johansen, A., Ida, S., & Brasser, R. 2019, A&A, 622, A202 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Kornet, K., Różyczka, M., & Stepinski, T. F. 2004, A&A, 417, 151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271 [NASA ADS] [CrossRef] [Google Scholar]
  49. Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [Google Scholar]
  50. Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Lambrechts, M., Morbidelli, A., Jacobson, S. A., et al. 2019, A&A, 627, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Lenz, C. T., Klahr, H., & Birnstiel, T. 2019, ApJ, 874, 36 [CrossRef] [Google Scholar]
  53. Liu, H. B. 2019, ApJ, 877, L22 [NASA ADS] [CrossRef] [Google Scholar]
  54. Manara, C. F., Robberto, M., Da Rio, N., et al. 2012, ApJ, 755, 154 [NASA ADS] [CrossRef] [Google Scholar]
  55. Manara, C. F., Testi, L., Natta, A., & Alcalá, J. M. 2015, A&A, 579, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Meru, F., Geretshauser, R. J., Schäfer, C., Speith, R., & Kley, W. 2013, MNRAS, 435, 2371 [NASA ADS] [CrossRef] [Google Scholar]
  57. Musiolik, G., & Wurm, G. 2019, ApJ, 873, 58 [NASA ADS] [CrossRef] [Google Scholar]
  58. Nayakshin, S. 2017, PASA, 34, e002 [NASA ADS] [CrossRef] [Google Scholar]
  59. Okuzumi, S., & Hirose, S. 2011, ApJ, 742, 65 [NASA ADS] [CrossRef] [Google Scholar]
  60. Okuzumi, S., Tanaka, H., Takeuchi, T., & Sakagami, M.-a. 2011, ApJ, 731, 95 [NASA ADS] [CrossRef] [Google Scholar]
  61. Ormel, C. W., & Klahr, H. H. 2010, A&A, 520, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Pavlyuchenkov, Y., Akimkin, V., Wiebe, D., & Vorobyov, E. 2019, MNRAS, 486, 3907 [NASA ADS] [CrossRef] [Google Scholar]
  63. Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519 [NASA ADS] [CrossRef] [Google Scholar]
  64. Pignatale, F. C., Charnoz, S., Chaussidon, M., & Jacquet, E. 2018, ApJ, 867, L23 [NASA ADS] [CrossRef] [Google Scholar]
  65. Pinte, C., Dent, W. R. F., Ménard, F., et al. 2016, ApJ, 816, 25 [NASA ADS] [CrossRef] [Google Scholar]
  66. Powell, D., Murray-Clay, R., & Schlichting, H. E. 2017, ApJ, 840, 93 [NASA ADS] [CrossRef] [Google Scholar]
  67. Powell, D., Murray-Clay, R., Pérez, L. M., Schlichting, H. E., & Rosenthal, M. 2019, ApJ, 878, 116 [NASA ADS] [CrossRef] [Google Scholar]
  68. Rice, W. K. M., Mayo, J. H., & Armitage, P. J. 2010, MNRAS, 402, 1740 [NASA ADS] [CrossRef] [Google Scholar]
  69. Ros, K., Johansen, A., Riipinen, I., & Schlesinger, D. 2019, A&A, 629, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Rosotti, G. P., Tazzari, M., Booth, R. A., et al. 2019, MNRAS, 486, 4829 [NASA ADS] [CrossRef] [Google Scholar]
  71. Safronov, V. S. 1972, Evolution of the Protoplanetary Cloud and Formation of the Earth and Planets (Israel: Program for Scientific Translations) [Google Scholar]
  72. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  73. Stamatellos, D. 2015, ApJ, 810, L11 [NASA ADS] [CrossRef] [Google Scholar]
  74. Steinpilz, T., Joeris, K., Jungmann, F., et al. 2019, Nat. Phys., 1 [Google Scholar]
  75. Stoyanovskaya, O. P., Vorobyov, E. I., & Snytnikov, V. N. 2018, Astron. Rep., 62, 455 [NASA ADS] [CrossRef] [Google Scholar]
  76. Suzuki, T. K., & Inutsuka, S.-i. 2009, ApJ, 691, L49 [NASA ADS] [CrossRef] [Google Scholar]
  77. Tanaka, Y. A., & Tsukamoto, Y. 2019, MNRAS, 484, 1574 [NASA ADS] [CrossRef] [Google Scholar]
  78. Teiser, J., & Wurm, G. 2009, A&A, 505, 351 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]
  80. Trapman, L., Facchini, S., Hogerheijde, M. R., van Dishoeck, E. F., & Bruderer, S. 2019, A&A, 629, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  81. Tsukamoto, Y., Okuzumi, S., & Kataoka, A. 2017, ApJ, 838, 151 [NASA ADS] [CrossRef] [Google Scholar]
  82. Turner, N. J., Fromang, S., Gammie, C., et al. 2014, Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning (Tucson, AZ: University of Arizona Press), 411 [Google Scholar]
  83. Venuti, L., Bouvier, J., Flaccomio, E., et al. 2014, A&A, 570, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Vericel, A., & Gonzalez, J.-F. 2020, MNRAS, 492, 210 [NASA ADS] [CrossRef] [Google Scholar]
  85. Vorobyov, E. I. 2013, A&A, 552, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Vorobyov, E. I., & Basu, S. 2005, ApJ, 633, L137 [NASA ADS] [CrossRef] [Google Scholar]
  87. Vorobyov, E. I., & Basu, S. 2006, ApJ, 650, 956 [NASA ADS] [CrossRef] [Google Scholar]
  88. Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896 [NASA ADS] [CrossRef] [Google Scholar]
  89. Vorobyov, E. I., & Elbakyan, V. G. 2018, A&A, 618, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  90. Vorobyov, E. I., & Elbakyan, V. G. 2019, A&A, 631, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Vorobyov, E. I., Elbakyan, V., Hosokawa, T., et al. 2017, A&A, 605, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Vorobyov, E. I., Akimkin, V., Stoyanovskaya, O., Pavlyuchenkov, Y., & Liu, H. B. 2018, A&A, 614, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  93. Vorobyov, E. I., Skliarevskii, A. M., Elbakyan, V. G., et al. 2019, A&A, 627, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  94. Wada, K., Tanaka, H., Okuzumi, S., et al. 2013, A&A, 559, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  95. Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2009, ApJ, 702, 1490 [NASA ADS] [CrossRef] [Google Scholar]
  96. Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
  97. Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius (New York: Wiley Interscience Division), 211 [Google Scholar]
  98. Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67 [NASA ADS] [CrossRef] [Google Scholar]
  99. Yamamoto, T., Kadono, T., & Wada, K. 2014, ApJ, 783, L36 [NASA ADS] [CrossRef] [Google Scholar]
  100. Yorke, H. W., & Bodenheimer, P. 2008, ASP Conf. Ser., 387, 189 [NASA ADS] [Google Scholar]
  101. Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
  102. Zhu, Z., Hartmann, L., & Gammie, C. 2010, ApJ, 713, 1143 [NASA ADS] [CrossRef] [Google Scholar]
  103. Zhu, Z., Hartmann, L., Nelson, R. P., & Gammie, C. F. 2012, ApJ, 746, 110 [NASA ADS] [CrossRef] [Google Scholar]
  104. Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJ, 877, L18 [NASA ADS] [CrossRef] [Google Scholar]
  105. Ziampras, A., Ataiee, S., Kley, W., Dullemond, C. P., & Baruteau, C. 2020, A&A, 633, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

1

The surface density for HD 163296 is derived using two different techniques (dust and CO lines) and here we show the both results.

All Tables

Table 1

Model parameters.

Table 2

a and b parameters for different threshold values of minimum dust size in the pebble definition.

Table 3

Similar to Table 2, but for the mass accretion rates of gas and pebbles onto the central star.

Table 4

Disk masses of the observed objects obtained from Powell et al. (2019) along with the disk masses of our models.

Table 5

Outer dust radii of the three observed disks for different ALMA bands obtained from Powell et al. (2019) and outer dust radius of our model L obtained from the reproduced intensities.

Table A.1

Parameters a and b for different threshold values of minimum dust size in the pebble definition.

All Figures

thumbnail Fig. 1

Gas surface density maps of the inner 1200 × 1200 au2 box at different time moments for all three models. The colorbar is shown in log scale. The insets show the Toomre Q-parameter for all azimuthal grid points at a specific radial distance from the star.

In the text
thumbnail Fig. 2

Radial distribution of the maximum radius of grown dust for all azimuthal grid points shown for consecutive times for all models. Color of dots shows the value of Stokes number for each azimuthal grid point. The colorbar is in log scale. Solid black line marks the dust fragmentation size afrag.

In the text
thumbnail Fig. 3

Azimuthally averaged absolute values of gas (blue curve), grown dust (red curve), and pebble (green curve) accretion (transport) rates versus radial distance from the star at distinct evolutionary times for all our models. Dashed part of the curve shows the outward migration, while the solid part – the inward migration.

In the text
thumbnail Fig. 4

Relation between pebble mass flux peb and gas mass flux gas for all our models. Color of the dots presents the age of the system in Myr. The black dashed line shows the best-fit curve for each model. The red dotted lines show the ± 3σ deviation from the best-fit values. The dash-dotted line shows the peb = 0.01g dependence.

In the text
thumbnail Fig. 5

Similar to Fig. 4, but for the mass accretion rates of gas and pebbles onto the central star.

In the text
thumbnail Fig. 6

Azimuthally averaged dust-to-gas mass ratio (black), surface density of gas (blue curve), grown dust (red curve), and pebbles (green curve) versus radial distance from the star at distinct evolutionary times for all our models.

In the text
thumbnail Fig. 7

Radial dependence of the azimuthally averaged surface densities of gas (blue curve), grown dust (red curve), and small dust (green curve) at evolutionary time t = 946 kyr for all ourmodels. The dashed vertical lines mark the radial distance at which the radius of dust particles ar = 0.85∕2π mm. Thin curves in top panel show the gas surface densities of three observed disks derived by Powell et al. (2019).

In the text
thumbnail Fig. 8

Synthetic intensities in ALMA Band 3 (2.80 mm), Band 6 (1.30 mm), and Band 7 (0.85 mm) and total dust surface density for model L at t = 964 kyr. The stars mark the radial distance at which dust size becomes equal to λ∕2π.

In the text
thumbnail Fig. A.1

Left column: radial distribution of the maximum radius of grown dust for all azimuthal grid points in model L2. Color of dots shows the value of Stokes number for each azimuthal grid point. Solid black line shows the dust fragmentation size afrag. Middle column: azimuthally averaged dust-to-gas mass ratio (black), surface density of gas (blue curve), grown dust (red curve),and pebbles (green curve) vs. radial distance from the star. The dashed line shows the power-law slope ∝ r−1.5. Right column: azimuthally averaged absolute values of gas (blue curve), grown dust (red curve), and pebble (green curve) accretion (transport) rates vs. radial distance from the star. Dashed part of the curve shows the outward migration, solid part – the inward migration.

In the text
thumbnail Fig. A.2

Relation between pebble mass flux peb and gas mass flux gas in the entiredisk (top panel) and onto the central star (bottom panel) for model L2. Color of dots presents the age of the system in Myr. The black dashed line shows the best-fit curve for the data in each panel. The red dotted lines show the ± 3σ deviation from the best-fit values. The dash-dotted line shows the M˙peb*=0.01M˙g*$\dot{M}^*_{\textrm{peb}}\,{=}\,0.01\dot{M}^*_{\textrm{g}}$ dependence.

In the text

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

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

Initial download of the metrics may take a while.