An opening criterion for dust gaps in protoplanetary discs

TL;DR

This paper develops an analytical criterion for dust gap opening by low-mass planets in protoplanetary discs, validated with 3D simulations, aiding in interpreting high-resolution disc observations.

Contribution

It extends dust radial drift theory to include planetary and viscous forces, providing new criteria for dust gap formation and planet mass estimates.

Findings

Analytical grain size-dependent dust gap opening criterion.

Estimate of the dust gap's outer edge location.

Minimum Stokes number for planet-induced dust gaps.

Abstract

We aim to understand under which conditions a low mass planet can open a gap in viscous dusty protoplanetary discs. For this purpose, we extend the theory of dust radial drift to include the contribution from the tides of an embedded planet and from the gas viscous forces. From this formalism, we derive i) a grain size-dependent criterion for dust gap opening in discs, ii) an estimate of the location of the outer edge of the dust gap and iii) an estimate of the minimum Stokes number above which low-mass planets are able to carve gaps which appear only in the dust disc. These analytical estimates are particularly helpful to appraise the minimum mass of an hypothetical planet carving gaps in discs observed at long wavelengths and high resolution. We validate the theory against 3D SPH simulations of planet-disc interaction in a broad range of dusty protoplanetary discs. We find a remarkable agreement between the theoretical model and the numerical experiments.

Equations

$$\Lambda (r)=\mathrm{sgn}(r-r_{\mathrm{p}})f\frac{(\mathcal{G}{M}_{\mathrm{p}}{)}^2}{{\Omega }_{\mathrm{p}}^2}\frac{1}{{\Delta }^4},$$

$$r_{\mathrm{H}}\equiv r_{\mathrm{p}}{\left(\frac{M_{\mathrm{p}}}{3M_{\star }}\right)}^{1/3}.$$

$${\Delta }_{\mathrm{g}}=\max (\left|r-r_{\mathrm{p}}\right|,H,r_{\mathrm{H}}).$$

$${\Delta }_{\mathrm{d}}=\max (\left|r-r_{\mathrm{p}}\right|,r_{\mathrm{H}}).$$

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

$$\frac{\partial {\mathbf{v}}_{\mathrm{d}}}{\partial t}+({\mathbf{v}}_{\mathrm{d}}\cdot \nabla ){\mathbf{v}}_{\mathrm{d}}$$

$$K\equiv \frac{{\rho }_{\mathrm{d}}}{t_{\mathrm{s}}(1+ϵ)},$$

$$\mathrm{St}=\frac{{\rho }_{\mathrm{grain}}{s}_{\mathrm{grain}}{\Omega }_{\mathrm{k}}}{(1+ϵ){\rho }_{\mathrm{g}}{c}_{\mathrm{s}}}\sqrt{\frac{\pi \gamma }{8}},$$

$$\sigma_{r\theta}=\eta r\frac{\partial}{\partial r}\biggr{(}\frac{v_{\mathrm{g},\theta}}{r}\biggl{)},$$

$$\frac{v_{\mathrm{g},\theta }^2}{r}=\frac{v_{\mathrm{k}}^2}{r}+\frac{1}{{\rho }_{\mathrm{g}}}\frac{\partial P}{\partial r}+\mathcal{O}{\left(H/r\right)}^2.$$

$$\frac{\partial v_{\mathrm{g},R}}{\partial t}$$

$$\frac{\partial v_{\mathrm{g},\varphi }}{\partial t}$$

$$\frac{\partial v_{\mathrm{d},R}}{\partial t}$$

$$\frac{\partial v_{\mathrm{d},\varphi }}{\partial t}$$

$$v_{\mathrm{g},r}=$$

$$+\frac{2{\Lambda }_{\mathrm{d}}}{v_{\mathrm{k}}}\frac{ϵ}{(1+ϵ)(1+{\mathrm{St}}^2)}+\frac{2{\Lambda }_{\mathrm{g}}}{v_{\mathrm{k}}}\frac{1+{\mathrm{St}}^2(1+ϵ)}{(1+ϵ)(1+{\mathrm{St}}^2)},$$

$$v_{\mathrm{d},r}=$$

$$+\frac{2{\Lambda }_{\mathrm{d}}}{v_{\mathrm{k}}}\frac{ϵ\left(1+{\mathrm{St}}^2\right)+{\mathrm{St}}^2}{(1+ϵ)\left(1+{\mathrm{St}}^2\right)}+\frac{2{\Lambda }_{\mathrm{g}}}{v_{\mathrm{k}}(1+ϵ)(1+{\mathrm{St}}^2)},$$

$$v_{\mathrm{g},\theta }=$$

$$+\frac{ϵ({\Lambda }_{\mathrm{g}}-{\Lambda }_{\mathrm{d}})}{v_{\mathrm{k}}(1+ϵ)(\mathrm{St}+{\mathrm{St}}^{-1})},$$

$$v_{\mathrm{d},\theta }=$$

$$-\frac{{\Lambda }_{\mathrm{g}}-{\Lambda }_{\mathrm{d}}}{v_{\mathrm{k}}(1+ϵ)(\mathrm{St}+{\mathrm{St}}^{-1})},$$

$$\Delta v\equiv \frac{1}{{\rho }_{\mathrm{g}}{\Omega }_{\mathrm{k}}}\frac{\partial P}{\partial r},$$

$$v_{\mathrm{visc}}$$

$$v_{\mathrm{d},\theta }-v_{\mathrm{g},\theta }$$

$${\Lambda }_{\mathrm{g}\to \mathrm{d}}$$

$${\Lambda }_{\mathrm{d}\to \mathrm{g}}=r\frac{K}{{\rho }_{\mathrm{g}}}\left(v_{\mathrm{d},\theta }-v_{\mathrm{g},\theta }\right)=-ϵ{\Lambda }_{\mathrm{g}\to \mathrm{d}}.$$

$${\nu }_{\mathrm{d},\mathrm{eff}}\equiv \frac{1}{1+ϵ}\frac{1}{1+{\mathrm{St}}^2}\nu .$$

$$\left|-\frac{1}{{\rho }_{\mathrm{g}}{\Omega }_{\mathrm{k}}}\frac{\partial P}{\partial r}\right|$$

$$\left|\frac{1}{{\rho }_{\mathrm{g}}{\Omega }_{\mathrm{k}}}\nabla \cdot \sigma {|}_{\theta }\right|$$

$${\left|\frac{\Lambda }{v_{\mathrm{k}}}\right|}_{\max }\sim \frac{r_{\mathrm{H}}^6}{r^2{\Delta }^4}{v}_{\mathrm{k}}.$$

$${\left|\frac{{\Lambda }_{\mathrm{g}}}{v_{\mathrm{k}}}\right|}_{\max }$$

$${\left|\frac{{\Lambda }_{\mathrm{d}}}{v_{\mathrm{k}}}\right|}_{\max }$$

$$v_{\mathrm{d},r}{|}_{\Delta v}\sim \frac{\mathrm{St}}{1+{\mathrm{St}}^2}{\left(\frac{H}{r}\right)}^2{v}_{\mathrm{k}},$$

$$v_{\mathrm{d},r}{|}_{{\Lambda }_{\mathrm{d}}}\sim \frac{{\mathrm{St}}^2}{1+{\mathrm{St}}^2}{\left(\frac{r_{\mathrm{H}}}{H}\right)}^2{\left(\frac{H}{r}\right)}^2{v}_{\mathrm{k}}.$$

$$\frac{M_{\mathrm{p}}}{M_{\star }}\gtrsim {\alpha }^{1/2}{\left(\frac{H}{r_{\mathrm{p}}}\right)}^{5/2}.$$

$$\left(\frac{M_{\rm p}}{M_{\star}}\right)_{\mathrm{th}}\ {\raise-2.15277pt\hbox{$\buildrel>\over{\sim}$}}\ 3\biggl{(}\frac{H}{r_{\mathrm{p}}}\biggr{)}^{3}.$$

Full text

††pagerange: An opening criterion for dust gaps in protoplanetary discs–Acknowledgements††pubyear: 2016

1 Introduction

Dust rings and gaps-like structures have been recently revealed by high-resolution observations in both young and evolved protoplanetary discs (ALMA Partnership et al., 2015; Canovas et al., 2016; van der Plas et al., 2017; de Boer et al., 2016; Ginski et al., 2016; van Boekel et al., 2017; Andrews et al., 2016; Isella et al., 2016; Fedele et al., 2017). Various mechanisms have been proposed to explain the origin of these structures. A first category of models invokes discs that are dynamically young and in which planets have not yet formed. In those, rings may originate from self-induced dust pile-ups (Gonzalez et al., 2015; Gonzalez et al., 2017), zonal flows (Flock et al., 2015; Béthune et al., 2016), rapid pebble growth around condensation fronts (Zhang et al., 2015), aggregate sintering (Okuzumi et al., 2016), large scale instabilities due to dust settling (Lorén-Aguilar & Bate, 2016) or secular gravitational instabilties (Takahashi & Inutsuka, 2016).

The alternative and more natural explanation is to interpret the rings as an observational signature of embedded planets. The tides generated by planets of sufficient masses overcome the gap-closing contributions induced by the pressure gradient and the viscous spreading of the gas. As a result, the planet carves a gap in its vicinity by pushing material away from its orbit (e.g. Goldreich & Tremaine, 1979, 1980; Lin & Papaloizou, 1986, 1993; Rafikov, 2002a; Kley & Nelson, 2012; Baruteau et al., 2014). In the prototypal case of the disc around HL Tau, this explanation is supported by some observational features, such as the increase of the gap eccentricity at large orbital radii, as well as spectral index variations between dark and bright rings which suggests that dark rings are regions of low dust density (ALMA Partnership et al., 2015). So far, the structures observed in this disc (ALMA Partnership et al., 2015) and HD163296 (Isella et al., 2016) have been better reproduced by assuming the presence of planets (e.g. Dipierro et al., 2015; Dong et al., 2015; Jin et al., 2016; Isella et al., 2016). This scenario would be consistent with the increasing number of extrasolar planets detected (Laughlin & Lissauer, 2015), but requires Saturn mass planets within a few million years at most, challenging the scenario of planet formation through core accretion.

A criterion on the minimum mass required for a planet to open a gap in a gas disc was derived by Crida et al. (2006). They considered the balance between the gap-opening tidal torque and the gap-closing viscous torque, taking into account the non-local deposition of angular momentum by the density waves excited by planets (e.g. Goodman & Rafikov, 2001; Rafikov, 2002a). Recent investigations have refined this analysis, showing that less massive planets may open gaps as well (Dong et al., 2011; Duffell & MacFadyen, 2012, 2013; Duffell & Dong, 2015; Zhu et al., 2013), but more massive planets may not (Malik et al., 2015).

However, observations of protoplanetary discs probe mostly the dust content of the disc, not the gas, and the two phases are not necessarily coupled. Recently, numerical investigations have shown that gap opening is more effective in the dust than in the gas (Paardekooper & Mellema, 2004, 2006; Fouchet et al., 2007, 2010; Ayliffe et al., 2012; Gonzalez et al., 2012; Zhu et al., 2014; Picogna & Kley, 2015; Dipierro et al., 2015, 2016; Rosotti et al., 2016).

Dust is usually modelled as an almost collisionless fluid in protoplanetary discs. As such, the dust intrinsic pressure and viscosity are not effective at closing the gap. Moreover, the tidal torque is amplified by geometrical effects due to dust settling. In the vicinity of the planet, the dust dynamics depends strongly on the size of the grains. When the planet is massive enough to carve a gap in the gas, micron-sized grains couple to the gas viscous flow and enter the gap (Rice. et al., 2006), whereas drag produces a pile up of large grains into the pressure maxima at the gap edges (e.g. Pinilla et al., 2015). This radial size-sorting produces well defined features in the mm and scattered light emission, consistent with recent observations (e.g. Follette et al., 2013; Canovas et al., 2016).

Recently, the disc around TW Hydrae has been observed by ALMA and SPHERE, probing the dust continumm emission at $850\,\mu m$ and the scattered light in the H-band at $1.6\,\mu m$ respectively (Andrews et al., 2016; van Boekel et al., 2017). A comparison of SPHERE and ALMA images reveals that the gaps at $\sim 37$ and $\sim 43$ au observed with ALMA are absent in the SPHERE image (see Fig. 7 of van Boekel et al. 2017). This mismatch between the distribution of large dust grains probed by ALMA and of small dust grains (expected to be well mixed with the gas) probed by SPHERE indicates that large dust grains may be more susceptible to gap formation than gas.

Surprisingly, the minimum planet mass required to open a gap in a dusty disc has not been clearly identified theoretically yet. Recently, Rosotti et al. (2016) have found that a shallow gap can be carved in the dust for planets able to slightly affect the local gas structure without creating pressure maxima. Moreover, Dipierro et al. (2016) have shown numerically that even lower mass planets could open gaps in the dust only, without any perturbation in the radial pressure gradient at the planet location. In this case, the creation of gap results from the competition between the tidal torque and the drag torque outside the planet orbit since the planet is not able to affect the gas structure. The drag torque acting on dust is negative all through the disc, whereas the tidal torque exerted by the planet is positive outside its orbit and negative inside. The balance between these two torques outside the planet orbit determines if the planet is able to carve a gap in the dust or not. In this paper, we propose a theory to model this mechanism, with the aim of deriving a grain size-dependent criterion for dust gap opening by non-migrating planets in protoplanetary discs. To this purpose, we extend the formalism of dust drift introduced in Nakagawa et al. (1986) to include viscous forces and the disc-planet tidal interactions. We also infer the radial location of the outer edge of the dust gap and the minimum Stokes number above which low mass planet are able to carve gap only in the dust. The results of our analysis are thoroughly tested against 3D Smoothed Particle Hydrodynamics (SPH) gas and dust simulations of different disc models.

The paper is organised as follows: in Sect. 2 we describe the dust dynamics in disc hosting a non-migrating planet under the action of tidal and drag forces. In Sect. 3 we apply the formalism developed in Sect. 2 to derive a simple gap opening criterion for the dust in the regime where the planet does not significantly alter the gas disc. In Sect. 4, we perform a set of simulation to test our model. In Sect. 5 we discuss how our criterion can be used in practice to interpret observations and, finally, in Sect. 6 we summarise our findings and conclusions.

2 Dust dynamics in a viscous disc with a planet

2.1 Disc-planet tidal interaction

2.1.1 Gas disc

A planet transfers angular momentum to its surrounding disc through the excitation of spiral density waves at specific locations called Lindblad resonances. Goldreich & Tremaine (1979, 1980) derived an analytic expression for the tidal torque per unit mass $\Lambda$ (the excitation torque density) exerted by a fixed planet on an elementary ring of a pressureless disc:

[TABLE]

where $M_{\mathrm{p}}$ is the mass of the planet, $\mathcal{G}$ is the gravitational constant, $r_{\mathrm{p}}$ is the planet location, $\Omega_{\mathrm{p}}=(\mathcal{G}M_{\star}/r_{\mathrm{p}}^{3})^{1/2}$ is the Keplerian angular velocity at $r_{\mathrm{p}}$ where $M_{\star}$ is the mass of the central star, $\Delta\equiv\left|r-r_{\mathrm{p}}\right|$, and $f$ is a constant of order unity (e.g. $f\sim 0.4$ - Goldreich & Tremaine 1979, $f\sim 0.15$ - Lin & Papaloizou 1979, $f\sim 0.1$ - Rafikov & Petrovich 2012). While the key features of the tidal torque equation, such as the $\Delta^{4}$ dependence, have been widely accepted (Lin & Papaloizou, 1986; Bryden et al., 1999; Bate et al., 2003; Varnière et al., 2004; D’Angelo & Lubow, 2008), the proportionality coefficient $f$ is mostly inferred by analyzing the shape of the gap carved by the planet (e.g. Armitage & Natarajan, 2002). The $\mathrm{sgn}$ factor in front of the right-hand side of Eq. 1 shows that a planet tends to push material outside of its orbit. In a disc with non-zero pressure, the exchange of torque between the planet and the disc is a two-step mechanism.

The planet first excites density waves in the gas. The amount of initial torque density stored in these density waves is the one given by Eq. 1. Then, a fraction of this torque (the deposition torque) is transferred from the waves to the disc by virtue of damping processes such as viscosity or shocks (Takeuchi et al. 1996; Goodman & Rafikov 2001; Rafikov 2002a). This pressure-supported transport modifies also the effective location of the Lindblad resonance. An important consequence is that the deposition torque is essentially zero in the vicinity of the planet, since high order Lindblad resonances are shifted away from the planet by a typical length of the order of the scale height of the disc $H$ (Goldreich & Tremaine, 1980; Artymowicz, 1993). The linear theory (valid for low mass planets, e.g. Dong et al. 2011) shows that there is no Lindblad resonances in the region $\left|r-r_{\mathrm{p}}\right|\lesssim 2H/3$. This effect, called torque cut-off, prevents the tidal torque to diverge close to the planet orbit. Finally, no torque is deposited by Lindblad resonances in the co-orbital region of the planet, where the gravity of the planet dominates over the one of the star and particles experience horseshoe orbits (Bate et al., 2003; D’Angelo & Lubow, 2008). This region typically extends over a distance $r_{\mathrm{H}}$ from the planet’s orbit, where $r_{\mathrm{H}}$ denotes the Hill radius of the planet:

[TABLE]

To include the torque cut-off due to the pressure and the corotation region, the deposition toque is often calculated from the expression given by Eq. 1, but with the following modified $\Delta_{\rm g}$

[TABLE]

The prescription given by Eq. 3 has been extensively used in a variety of contexts, such as to model the disc-planet interaction (e.g. Bryden et al., 1999) or the evolution of supermassive black hole binaries embedded in a disc (e.g. Lodato et al., 2009).

In this work, we use Eq. 1 to treat the tidal torque density by a prescription with is analytically tractable. This expression assumes a few simplification that should be kept in mind. Firstly, Eq. 1 gives the tidal torque under the form of a smooth function. This regularity property originates from the fact that individual modes of the tidal potential interact with large regions of the disc rather than narrow rings centred over Lindblad resonances (Rafikov & Petrovich, 2012). Secondly, Eq. 1 neglects effects due to the non-linear propagation of the density waves excited by the planet, some details of the wave damping processes and the recently discovered negative torque correction. About the latter, in a disc of uniform surface density, the tidal torque density changes sign at a radial separation $\gtrsim 3H$ from the planet (Dong et al., 2011; Rafikov & Petrovich, 2012). However, we do not expect this correction to play a major role since we study processes developing in regions of size $\lesssim H$ around the planet (see Sect. 5.2), where the tidal torque exerted by the planet on the gas is always positive in the outer disc and negative in the inner disc. Moreover, we focus here on low mass planets embedded in viscous discs, where density waves are expected to be strongly damped close to the Lindblad resonsances. In this case, the exact role of the negative torque phenomenon remains an open question (Sect. 7 of Rafikov & Petrovich 2012). Hence, we choose to restrain our theory to a minimum but sufficient level of refinement.

2.1.2 Dust disc

From the discussion in Sect. 2.1.1, we expect the tidal torque to be more effective in the dust than in the gas for three reasons:

Dust is a pressureless fluid, where density waves cannot propagate far from the planet. The deposition torque equals therefore the excitation torque. 2. 2.

For the same reason, there is no torque cut-off in the dust, and the effective $\Delta_{\rm d}$ is

[TABLE]

Close to the planet, angular momentum can be deposited in the dust, but not in the gas. In detail, if $r_{\mathrm{H}}<H$, the torque exerted at the Lindblad resonances in the region between $r_{\mathrm{H}}$ and $H$ is effective in the dust, while in the gas is suppressed due to pressure effects. The ratio between the maximum torque density in the gas and in the dust is of order $\sim(r_{\mathrm{H}}/H)^{4}=(r_{\mathrm{p}}/H)^{4}(M_{\mathrm{p}}/3M_{\star})^{4/3}$. 3. 3.

For large grains, the thickness of the dust layer $H_{\mathrm{d}}$ is smaller than $H$, as a result from the competition between settling and turbulent stirring (Dubrulle et al., 1995). For dust layers with $H_{\mathrm{d}}\geq r_{\mathrm{H}}$ this results in an enhancement of the tidal torque due to local geometrical effects.

In absence of gas, the use of the smooth functional form given by Eq. 1 is not appropriate. Indeed, the tidal torque exerted on the pressureless dust phase concentrates at Lindblad resonances, where particles eccentricities are effectively excited (Ayliffe et al., 2012; Zhu et al., 2014). In this case, the planet-disc interaction is better described by Hill’s equations (Hill, 1878), additional effects related to dust pressure induced by mutual collisions and velocity fluctuations must be taken into account (Henon & Petit, 1986; Petit & Henon, 1987a, b; Rafikov, 2001). However, for a disc containing a low-mass planet, two arguments support the use of Eq. 1 for dust as well. Firstly, grains experience gas drag (see Sect. 2.2), a force which dominates the dynamics of the particles at the resonances as long as they are not too large (Fig. 18, top-center panel of Zhu et al. 2014). In this case, orbits of dust grains shrink around the planet and the ability for opening a dust gap is enhanced (Ayliffe et al., 2012; Zhu et al., 2014). Secondly, with a low-mass planet, $r_{\mathrm{H}}<H$, we will show that the width of the dust gap $\Delta_{\mathrm{gap}}$ satisfies $r_{\mathrm{H}}\lesssim\Delta_{\mathrm{gap}}\lesssim H$ (see Sect. 5.2). In this region, high-order Lindblad resonances of order $m$ are highly concentrated (r_{\mathrm{p}}/H\ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ m\ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ r_{\mathrm{p}}/r_{\mathrm{H}}) and degenerate into a continuum. This region is also sufficiently far away from low-order resonances (i.e.1:2, 2:3, 3:2 and 2:1), where eccentricity pumping is effective and can not be neglected (Zhu et al., 2014; Ayliffe et al., 2012). Thus, we use Eq. 1 to model the tidal torque density in the dust as well, and test this assumption with numerical simulations.

2.2 Equations of motion

The motion of dust particles in protoplanetary discs is affected by the aerodynamical interaction with the gas and vice-versa. The differential motion between the two phases gives rise to a drag force that damp this velocity difference. In discs hosting planets, an additional velocity difference can be produced due to the different tidal interaction of the planet with the gas and dust (Sect. 2.1).

We assume a thin, non-magnetic, non-self-graviting, dusty viscous and vertically isothermal protoplanetary disc hosting a non-migrating planet. We treat the dust phase as a continuous pressureless and viscousless fluid. The equations of motion for the gas and the dust are

[TABLE]

where the indices g and d refer to the gas and the dust phases, $v$ and $\rho$ denote the velocities and the densities, $\Phi$ and $\Phi_{\mathrm{p}}$ denote the gravitational potentials of the star and the planet. $P$ and $\sigma$ denote the pressure and the viscous tensor of the gas. $K$ denotes the drag coefficient, whose expression depends on the local values of the parameters of the grain and of the disc (e.g. Laibe & Price 2012). $K$ is related to the stopping time of the mixture $t_{\rm s}$ by

[TABLE]

where $\epsilon=\rho_{\mathrm{d}}/\rho_{\mathrm{g}}$ is the dust-to-gas density ratio. Instead of $K$, one uses generally the Stokes number $\mathrm{St}\equiv\Omega_{\mathrm{k}}t_{\mathrm{s}}$, defined as the ratio of the stopping time to the local dynamical timescale. In typical discs, the mean free path of the gas molecules is smaller than the dust particle size $s_{\mathrm{grain}}$. For this so-called Epstein regime, the Stokes number is given by (e.g. Price & Laibe, 2015)

[TABLE]

where $\rho_{\mathrm{grain}}$ is the intrinsic grain density, $c_{\mathrm{s}}$ is the sound speed and $\gamma$ is the adiabatic index.

To reduce Eqs. 5 – 6 to a system of equations that describes a steady steady solution for the gas and the dust, we follow the approach introduced in Nakagawa et al. (1986) (hereafter NSH86) and make two approximations. Firstly, we assume that the orbits have circularised after a transient regime occurring over a time $t_{\rm s}$ (Adachi et al., 1976). Secondly, the long term evolution of the gas surface density profiles due to viscous effect is neglected. To include the tides from the planets, we only consider the azimuthally averaged contribution of the tidal torque by replacing the source term $\mathbf{\nabla}\Phi_{\mathrm{p}}|_{\theta}$ by $\Lambda/r$, where $\Lambda$ is given by Eq. 1, adopting the prescription expressed in Eqs. 3 - 4 for each phases. For low mass planets, i.e. $r_{\mathrm{H}}<H$, the values of $\Lambda_{\rm g}$ and $\Lambda_{\rm d}$ are a priori different (the cut-off is at $H$ for the gas and at $r_{\rm H}$ for the dust). We assume classical shear viscosity and $\eta=\nu\rho_{\mathrm{g}}$ to denote the dynamical viscosity of the gas. This approach is commonly used to model angular momentum transport driven by turbulence generated e.g. by magneto-rotational or gravitational instabilities (Balbus & Hawley, 1991; Lodato & Rice, 2004; Rafikov, 2015). However, Rafikov (2017) has recently shown that accretion might not proceed viscously in protoplanetary discs, but may be driven non-diffusively by magnetohydrodynamic winds or spiral density waves (Rafikov, 2002a; Bai, 2016; Fung & Chiang, 2017). In this paper, we assume that accretion is mediated by a viscous-like mechanism (Lynden-Bell & Pringle, 1974). The only non-zero component of $\sigma$ for the axisymmetric sheared flow of the disc is

[TABLE]

where

[TABLE]

For low mass planets, the deviation of $\sigma_{r\theta}$ induced by the sub-Keplerian rotation of the unperturbed pressure profile is only of order $(H/R)^{2}\ll 1$ and can be neglected.

Under these assumptions, we perform a perturbative expansion relative to the Keplerian velocity $\mathbf{v_{\mathrm{k}}}=(0,r\Omega_{\mathrm{k}},0)$ in both phases, and look for stationary solutions assuming an axisymmetric disc. Note that from here, the notation $v_{\theta}$ will be used to refer to the perturbed azimuthal velocities for sake of simplicity. The equations of motion for the perturbed velocities can be expressed in polar coordinates as

[TABLE]

Eqs. 11 – 14 consist of a system of non-homogeneous differential equations of the form $\mathbf{X}^{\prime}(t)+\mathbf{A}\mathbf{X}(t)=\mathbf{B}$, where $\mathbf{A}$ and $\mathbf{B}$ are two constant matrices. Its steady state is $\mathbf{X}_{\rm stat}=\mathbf{A}^{-1}\mathbf{B}$. The system relaxes towards this stationary regime in the typical time $\min\left|\Re\left(\sigma_{\rm A}\right)\right|^{-1}=t_{\mathrm{s}}$, where $\sigma_{\rm A}$ are the eigenvalues of matrix $\mathbf{A}$ and $t_{\mathrm{s}}$ is the the stopping time of the mixture defined in in Eq. 8 in a dimensionless form.

2.2.1 Steady-state solution

The stationary solution of the linear system Eqs. 11 – 14 is

[TABLE]

where

[TABLE]

is the typical optimal drift velocity derived in NSH86, and

[TABLE]

is the viscous velocity derived by Lynden-Bell & Pringle (1974), since for low mass planet, surface density gradients develop over large scales (no gas gap). From Eqs. 17 – 18, the differential azimuthal velocity between the gas and dust is

[TABLE]

which is independent on the dust-to-gas ratio. The specific drag torque exerted by the gas phase on an elementary dust ring si therefore

[TABLE]

while the back reaction drag torque from the dust to the gas is

[TABLE]

2.2.2 Physical interpretation

In the limit of the dust grains being perfectly decoupled from the gas ($\mathrm{St}=+\infty$), particles orbit the star with Keplerian velocity ($v_{\mathrm{d},\theta}=0$) and are pushed outside of the planet orbit by the tidal torque at the constant velocity $v_{\mathrm{d},r}=2\Lambda_{\rm d}/v_{\rm k}$. The gas orbits the star at the sub-Keplerian velocity $\Delta v/2<0$, while the gas radial velocity is $v_{\mathrm{visc}}+2\Lambda_{\rm g}/v_{\rm k}$. When drag couples the two phases, the dust motion is dominated by the gas when $\epsilon\ll 1$ and $1/\left(1+\epsilon\right)\sim 1$. When $\epsilon\gg 1$ the gas motion is dominated by the dust. Both the gas and the dust strongly feel the other phase when the dust-to-gas ratio is of order unity (Gonzalez et al., 2017).

As it is known in in absence of a planet ($\Lambda_{\rm g}=\Lambda_{\rm d}=0$), differential dynamics due to gas pressure and viscosity are communicated from one phase to the other via drag. The term proportional to $\Delta v$ in Eq. 16 implies that particles drift radially as a result of the residual differential orbital velocity between the two phases (Eq. 2.11 of NSH86). In the absence of any small scale pressure perturbation in the gas, $\Delta v<0$ in a typical disc since the inner regions are denser and warmer. This results in a radial inward drift of particles. Drift is most efficient for $\mathrm{St}\sim 1$. For $\mathrm{St}\gg 1$ (resp. $\mathrm{St}\ll 1$) particles maintain a fixed Keplerian (resp. sub-Keplerian) orbit. In Eq. 15, the term proportional to $\Delta v$ corresponds to the back-reaction term (Eq. 2.13 of NSH86): in discs with $\Delta v<0$, the gas is pushed back in the regions of lower pressure by the dust to conserve the global angular momentum (this motion becomes significant when $\epsilon$ approaches unity). Futhermore, the radial viscous gas flow is weighted by a factor which depends on $\epsilon$ and $\mathrm{St}$ (second term of the right-hand side of Eq. 15). For small particles $\mathrm{St}\ll 1$, the viscous flow is reduces by a factor $\left(1+\epsilon\right)^{-1}$. For large particles ($\mathrm{St}\gg 1$), $v_{\mathrm{g},r}=v_{\mathrm{visc}}$, i.e. dust does non affect the viscous gas motion regardless the value of $\epsilon$.

In Eq. 16, the term proportional to $v_{\mathrm{visc}}$ shows that additionally, drag makes very small grains stick to the gas. They are carried radially by the viscous flow. Hence, dust evolves “viscously” under the indirect effect of the viscous evolution of the gas, with a so-called “drag-induced dust viscosity”

[TABLE]

This contribution dominates over the pressure drift only when grains are tiny, i.e. $\mathrm{St}<\alpha$. Note that although Eq. 24 resembles to the dust diffusivity derived in Youdin & Lithwick (2007), the drag-induced dust viscosity does not describe the motion induced by turbulence over dust grains, but models how the viscous evolution of the gas affects the radial dynamics of grains. Importantly, when $\mathrm{St}\sim 1$ and $\epsilon\geq 1$, the gas dynamics is dominated by the back-reaction from the dust and not by the viscosity (Gonzalez et al., 2015; Taki et al., 2016).

In the general case, the differential motion between the gas and the dust induced by the different tides only affect the motion of smaller grains, i.e. \mathrm{St}\ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ 1. This concerns particularly the region $r_{\rm H}<\left|r-r_{\rm p}\right|<H$, where $\Lambda_{\rm d}\neq\Lambda_{\rm g}=0$. In this region, for a gas-dominated dynamics with $\epsilon\ll 1$, only large dust grains, i.e. \mathrm{St}\ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ 1, are pushed away from the planet orbit, since they experience a lower drag related to the different tidal torque between the two phases (see the third term in the right hand side of Eq. 22). In other words, for smaller grains, the motion induced by the tidal torque in this region is damped by the drag torque that tends to reduce the velocity difference with the unperturbed gas flow. As a result, small grains are forced to stick to the fixed gas. On the other hand, it can be noticed that the gas is not affected by the tides in this region if $\epsilon\ll 1$ (see the third and fourth terms in the right hand side of Eq. 15). Sufficiently far away from the planet, where $\left|r-r_{\rm p}\right|>H$ and thus $\Lambda_{\rm g}=\Lambda_{\rm d}$, there is no differential tidal torque, and the tidal barycentric velocity is spread over the two phases proportionally to the respective density of each phase.

2.2.3 Orders of magnitude

We compare the orders of magnitude of the different velocity terms in Eq. 16, related to the radial pressure gradient, the viscous and the tidal contribution. Until specified, we shall not restrict our analysis to the case of an unperturbed gas density profile. We obtain

[TABLE]

where $l$ denotes the typical length over which the gas surface density varies and where we have used the seminal Prandtl-like turbulent viscosity $\nu=\alpha c_{\rm s}H$ (Shakura & Sunyaev, 1973). The prefactor in Eq. 26 originates from the second derivative of the Keplerian deviation (Eq. 10). For an unperturbed gas profile, $l=r$ and the pre-factors in Eqs. 25 – 26 equal unity. The maximal tidal contributions, obtained at the cut-off locations, are of order

[TABLE]

For low mass planets, i.e. $r_{\mathrm{H}}<H$, the maximum tidal torque density for the gas and the dust are given by

[TABLE]

For the sake of clarity, we now assume $\epsilon\ll 1$ and limit our analysis to large grains which are most affected by the tidal torque without being decelerated by the drag torque arising due to the differential tidal torque between the two phases. For larger grains, i.e. \mathrm{St}\ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ 1, the terms in Eq. 16 proportional to $1/(1+\mathrm{St}^{2})$ are negligible. The remaining terms in Eq. 16 are given by

[TABLE]

[TABLE]

The term $\mathrm{St}^{2}/\left(1+\mathrm{St}^{2}\right)$ in Eq. 31 expresses that only large grains are entrained by the tides since small grains stick to the gas. For these grains the drag torque is dominated by the contributions from the pressure gradient and the tides.

3 Gap opening in dusty discs

3.1 The low mass planet regime

The formalism derived in Sect. 2 enables to study the gap opening process by a low massive planet embedded in dusty discs assuming that the pressure profile around the planet remains unperturbed. We now investigate under which condition this assumption is satisfied.

A gap is carved in the gas when the tidal torque overpowers the viscous torque. Assuming that the typical length over which the gas surface density varies by the tidal action of the planet is of order of $H$, we compare Eqs. 26 and 28 by replacing $l$ by $H$ to estimate the condition for gap opening in the gas. We find

[TABLE]

Although Eq. 32 provides an interesting scaling, more quantitative criteria are used in literature. The first one is based on the requirement that a strong shock forms within a scale height of the planet’s orbit (Lin & Papaloizou, 1993), giving

[TABLE]

according to which the planet Hill radius $r_{\mathrm{H}}$ must be greater than the vertical scale height of the disc $H$. However, recent 2D and 3D simulations of gas discs hosting planets have shown that planets with mass M_{\mathrm{p}}\ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ 0.2\,M_{\mathrm{p,th}} are able to create a pressure maximum outside the planetary orbit (Lambrechts et al., 2014; Rosotti et al., 2016). The second criterion is based on the requirement that gap opening should be faster than the viscous refilling of the gap, giving

[TABLE]

where $f$ is the constant in the tidal torque density formula (Eq. 1). Moreover, Eq. 34 is in agreement with the estimate derived in Eq. 32. If the pressure (resp. the viscous) force dominates the gap closing mechanism, we expect $M_{\mathrm{p,th}}$ to be larger (resp. smaller) than $M_{\mathrm{p,visc}}$. In this work, we consider that the minimum mass able to create a pressure maxima in the outer disc is given by the maximum of all the masses predicted by the previously mentioned criteria,

[TABLE]

However, Rosotti et al. (2016) have recently found that planet of masses slightly lower than the one given by Eq. 35 could create gap structures in the dust as well. In detail, if $M_{\mathrm{p}}\gtrsim 0.1\,M_{\mathrm{p,th}}$, the planet weakens the pressure gradient profile in its neighbourhood, which reduces the dust drift locally and leads to accumulation of particles. This traffic jam mechanism affects essentially marginally coupled particles ($\mathrm{St}\sim 1$) and lead to the formation of a dust gap. We consider this effect by assuming that the minimum planet mass able to affect the local gas profile is

[TABLE]

In Sect. 4, we confirm the validity of this condition. Eventually, it should be noted that all the criteria above assume planets remains on a fixed orbit. However, Malik et al. (2015) have shown that the migration may affect the ability of the planet to carve gaps, and the critical masses given by Eqs. 33 and 34 might be underestimated.

3.2 Dust gap width

We now focus our analysis on planets not able to affect the local pressure structure, i.e. $M_{\mathrm{p}}\lesssim M_{\mathrm{p,lim}}$, embedded in standard discs ($\partial P/\partial r<0$). The tidal interaction between the planet and the disc acts to carve the gap around the planetary orbit, whereas drag makes the grains drift inwards towards the central protostar. In particular, the flux of solids coming from the outside of the planet orbit tends to refill the dust depletion locally induced by the tides. In this case, the drag torque can be derived using the unperturbed pressure profile of the gas (Eq. 22).

We denote

[TABLE]

the exponent which characterises the steepness of the pressure profile of the disc. If we assume power-law profiles for the surface density ($\Sigma\propto r^{-p}$) and the temperature ($T\propto r^{-q}$), we have $\zeta=-(p+q/2+3/2)$. For typical disc ($\,p=1$, $q=1/2$, Andrews & Williams 2007; Williams & Best 2014), the power law exponent of the pressure profile is $\zeta\simeq-2.75<0$. With this notation, the azimuthally averaged radial dust velocity (Eq. 16) can be written as

[TABLE]

As an example, Fig. 1 shows the radial velocity of millimetre sized grains outside the orbit of the planet with different mass (lower than $M_{\mathrm{p,lim}}$) embedded in a typical disc model, assuming $f$ to be equal to the nominal value of $0.4$ introduced in Goldreich & Tremaine (1979). In this case, the planet is located at 40 au from the star, $M_{\mathrm{p,th}}\sim 0.4\,M_{\mathrm{J}}$, $\alpha=0.005$, $p=1$, $q=1/2$, $H/r=0.02$ at 1 au (corresponding to $H/r_{\mathrm{p}}\approx 0.05$ at $r_{\mathrm{p}}$) and $\epsilon=0.01$. The dust grains in the midplane have Stokes number $\mathrm{St}\sim 10$ at the planet location.

Fig. 1 shows that in the presence of a planet of very low-mass, the tidal torque is not strong enough to halt the inward radial flow triggered by the drag torque (e.g. the case \ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ 0.003M_{\rm J}, for which $v_{\mathrm{d},r}<0$). Larger planet masses empower stronger tidal torques in the disc, and the balance between the tidal and the drag torques leads to an outward drift (Fig. 1 shows regions where $v_{\mathrm{d},r}>0$ for $M_{\rm p}=0.01M_{\rm J}$ and $M_{\rm p}=0.03M_{\rm J}$).

The location of the outer edge of the gap $r_{\rm gap}$ can be estimated simply, by evaluating the distance to the planet where $v_{\mathrm{d},r}(r_{\rm gap})=0$, i.e. where the drift induced by the tides balances steadily the drift induced by the drag. We use $\Delta_{\mathrm{gap}}$ to denote the width of the dust gap outside the planetary orbit $r_{\rm gap}-r_{\rm p}$. To simplify the analysis, we assume that the temperature and surface density of the gas are uniform over the domain $\Delta_{\mathrm{gap}}$, i.e. the terms in Eq. 38 are evaluated at $r_{\mathrm{p}}$ with exception of $\Delta$ in $\Lambda_{\mathrm{d}}$ and $\Lambda_{\mathrm{g}}$ (see Eq. 1). Corrections of order $\Delta_{\mathrm{gap}}/r_{\mathrm{p}}$ are negligible for our analysis (see Sect. 4.3). The width of the dust gap is therefore given by

[TABLE]

where

[TABLE]

Note that the difference between the two expressions in Eq. 40 has no sensible effect on the value of $z$ in practice for larger grains. This suggests, as expected, that the different tidal torque experienced by the two phases does not cover the key role in the gap opening process. For $\mathrm{St}\gtrsim\alpha$, the second term of the numerator vanishes and for $\mathrm{St}\gtrsim 1+\epsilon$, Eq. 40 reduces to

[TABLE]

Eq. 39 is consistent with the fact that infinitely small grains at low dust densities follow the gas, in which no gap forms ($\Delta_{\mathrm{gap}}=0$ for $\mathrm{St}=0$ and $\epsilon=0$). For large grains, $\Delta_{\mathrm{gap}}\propto\mathrm{St}^{1/4}$, an expression only weakly sensitive to the exact value of $\mathrm{St}$. The location of the outer edge of the dust gap is a weak increasing function of the Stokes number, since it corresponds to a weaker drag torque, as found numerically by Dipierro et al. (2016) (see Sect. 4.3.3). Moreover, this analysis shows that, since large grains experience a lower drag related to the differential tidal torque between the two phases, the criterion is mostly based on the balance between the tidal torque $\Lambda_{\mathrm{d}}$ and the drag torque related to pressure forces. The differential tidal torque between the two phases mostly influence the dynamics of small grains, forcing them to follow the unperturbed gas flow. For large Stokes number, i.e. $\mathrm{St}>\alpha$, the increase of the dust-to-gas ratio at the midplane due to settling affect the motion of dust and gas, as described in Sect. 2.2. However, Eq. 39 shows that, since $\mathrm{St}\propto 1/\left(1+\epsilon\right)$, the gap width does not depend on the dust-to-gas ratio when large grains are considered, i.e. when Eq. 41 is valid. The eventual modification of the gas surface density profile due to the dust back reaction is not included in our model (see discussion in Sect. 5.2).

Finally, Eq. 39 shows that the distance between the planet and the outer edge of the gap is a decreasing function of the aspect ratio of the disc. Indeed the contribution of gas drag to gap closing goes as $\left(H/r_{\mathrm{p}}\right)^{2}$ (cf Eq. 38). Assuming power-law profiles for the surface density of the gas $\Sigma\propto r^{-p}$ and temperature $T\propto r^{-q}$, the location of outer gap edge scales as $\Delta_{\mathrm{gap}}/r_{\mathrm{p}}\propto r^{p+\left(q-1\right)/4}$. For a typical disc with $p=1$ and $q=1/2$, the outer gap edge is expected to increase with radius as $r^{0.87}$.

3.3 Gap opening criterion

3.3.1 Orders of magnitude

An order of magnitude for the minimum mass required for a planet to open a dust gap can be straightforwardly estimated from Sect. 2.2.3. By equating Eqs. 30 and 31 with the definition given in Eq. 2, we obtain

[TABLE]

Eq. 42 shows that the critical mass required to open a gap in the dust is lower for large Stokes numbers and large aspect ratios. More interestingly, the exponent $-3/2$ implies a sharp transition between small and large grains at $\mathrm{St}\simeq 1$. Fixing $H/r_{\rm p}$ and decreasing $\mathrm{St}$ in Eq. 42 shows that $M_{\rm p}$ increases efficiently, up to reach the critical mass required to open a gap in the gas as well, i.e. $M_{\rm p}/M_{\rm\star}\sim\left(H/r_{\rm p}\right)^{3}$ (Eq. 33, Lin & Papaloizou 1993). This suggests, as expected, that the condition \mathrm{St}\ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ 1 should be fulfilled for a gap to be carved in the dust only (see Sect. 3.4).

3.3.2 Necessary condition for dust gap opening

Eq. 39 provides a simple way to estimate the minimum planet mass $M_{\rm p}$ able to halt the inward radial drift induced by the drag. Noting that the minimum radius of the outer edge of the dust gap is the Hill radius, $M_{\rm p}$ is the planet mass for which the radial dust velocity is zero at $r=r_{\mathrm{p}}+r_{\mathrm{H}}$. Assuming that the temperature and surface density of the gas are uniform over the domain $\left|r-r_{\mathrm{p}}\right|\sim r_{\mathrm{H}}$ – an approximation of order $r_{\rm H}/r_{\rm p}\ll 1$ – we obtain

[TABLE]

with $z$ given by the first expression in Eq. 40 and with

[TABLE]

Eq. 43 provides the minimum mass for a planet to produce a density depletion in the dust outside the planet orbit. For the typical disc described above, the minimum planet mass has a value of $\sim 0.003\,M_{\mathrm{J}}$, as expected from the previous analysis about the radial velocity (see dotted line in Fig. 1). As expected, Eq. 43 is in agreement with the orders-of-magnitude estimate (see Eq. 42).

3.3.3 Sufficient condition for dust gap opening

Eq. 43 gives the minimum mass of a planet able to halt the radial inward drift of dust particles. However, for a gap to form, dust must not refill dust depletions as it explores different azimuths. Since refilling is a non-axisymmetric process, its effects is not picked up by the previous analysis based on averaged axisymmetric torques. An alternative approach based on timescales estimates is therefore developed hereafter to take dust refilling into account.

The gap opening timescale is the time required to evacuate all the dust contained between $r_{\mathrm{p}}$ and $r_{\mathrm{p}}+r_{\mathrm{H}}$. As previously mentioned, a gap of half width $r_{\mathrm{H}}$ is roughly the smallest gap that can be opened in the dust, since the Lindblad resonances are most effective at this distance from the planet. Focussing our analysis on the dust outside the planetary orbit, the angular momentum that must be removed to open the gap between $r_{\mathrm{p}}$ and $r_{\mathrm{p}}+r_{\mathrm{H}}$ is

[TABLE]

where $l$ denotes the specific angular momentum. The typical time to evacuate all the dust in this region is

[TABLE]

where $\left|\mathrm{d}J/\mathrm{d}t\right|$ is the one-sided total torque on the planet due to its interaction with dust outside the orbit. This total torque is the integral of the torque density given by Eq. 1 over the entire outer disc, i.e.

[TABLE]

The gap opening timescale is therefore given by

[TABLE]

Setting $\Lambda_{\rm g}=\Lambda_{\rm d}=0$ in Eq. 38 and assuming $v_{\mathrm{d},r}$ constant over the domain $\left|r-r_{\mathrm{p}}\right|\sim r_{\mathrm{H}}$, the closing time $t_{\rm close}=r_{\rm H}/v_{\mathrm{d},r}$ is

[TABLE]

The critical mass ratio $M_{\mathrm{p}}/M_{\star}$ above which a planet sustains its gap in a dust disc is obtained by equating the opening and the closing timescale, which gives

[TABLE]

with $z$ given by the second expression in Eq. 40. As already mentioned, the two expressions in Eq. 40 have the same value for our aims, i.e. for $\mathrm{St}\gtrsim\alpha$ and for $\mathrm{St}\gtrsim\epsilon$. Thus, Eq. 50 equals to $3^{3/2}\approx 5.2$ times the critical mass derived in Eq. 43. For the disc model described above the criterion gives a typical mass of $\sim 0.015\,M_{\mathrm{J}}$ to open a dust gap. For large grains, Eq. 50 reduces to

[TABLE]

This criterion provides a good estimator for the minimum mass to open a gap in dusty disc. As a remark, Eq. 43 originates from a balance of torques performed at steady-state (the gap is already opened), whereas Eq. 50 originates from a balance of torques performed in a transient regime (the gap is not opened yet). This explains why the two conditions are not rigorously identical and differ by a factor of order unity. Note that, since $\mathrm{St}\propto 1/\left(1+\epsilon\right)$, the gap opening criterion does not depend on the dust-to-gas ratio, assuming that the dust back reaction does not affect the local pressure profile. Therefore, our analysis shows that it is not necessary to estimate the local dust-to-gas ratio to derive the value of the minimum mass.

Fig. 2 displays the minimum mass required for a planet to open a gap in the dust as a function of the Stokes number. This limit is calculated from Eq. 50 in the disc model described above assuming $f=0.4$ (Goldreich & Tremaine, 1979). The red shaded area indicates the range of planet masses and Stokes numbers for which a gap is carved in the dust only. There is no ubiquitous lower mass for gap opening in the dust, since tides always overpower drag in the limit of very large and decoupled grains. The green area shows the domain for which the planet carves a gap in the gas as well, i.e. for $M_{\mathrm{p}}\geq M_{\mathrm{p,gap}}$. The small blue area indicates the range of masses $M_{\mathrm{p,lim}}\lesssim M_{\mathrm{p}}\lesssim M_{\mathrm{p,gap}}$ for which the local pressure profile is perturbed without creating a pressure maximum (Rosotti et al., 2016).

As expected, the range of masses for which the planet is able to carve a gap in the dust only (red area) increases with increasing Stokes numbers, due to the reduced replenishment from the outer disc induced by the drag torque. For planets with masses inside the green area, the drag assists the gap opening in the dust, leading to an accumulation of dust particles at the pressure maximum and producing a well-defined dusty gap with a shape closely related to the Stokes number (e.g. Pinilla et al., 2015).

Eventually, the exact value of the minimum mass is related to the constant $f$ in front of the tidal torque density (Eq. 1). Eq. 51 shows that the critical mass to open a gap in the dust is sensible to the actual value of $f$ as it varies as $f^{-3/2}$ (see Eq. 44), as long as the hypothesis of the low-mass planet regime is satisfied. Measuring the critical mass provides therefore an effective way to measure $f$ in numerical simulations (see Sect. 4.3.1).

3.4 Critical Stokes number

An important feature mentioned in Sect. 2.2.3 appears clearly: low mass planets can carve dust gaps only if the Stokes number is above a critical value of order unity. The value of the critical Stokes number is obtained by comparing the minimum mass given by Eq. 51 to $M_{\mathrm{p,lim}}$, the minimum mass able to perturb the gas pressure profile. Since the value of the critical Stokes number is expected to be higher than unity, we assume $\mathrm{St}\gtrsim\alpha$ and $\mathrm{St}\gtrsim\epsilon$ and obtain

[TABLE]

For a typical disc with $\zeta=-2.75$ and $\epsilon=0.01$, $\mathrm{St_{crit}}\sim 5.2$ (see Fig. 2). Importantly, the value of the critical Stokes number does not depend on $H/r_{\mathrm{p}}$ since the gap opening criterion for the gas and the dust scales equally with the aspect ratio. Moreover, Eq. 37 shows that $\mathrm{St_{crit}}$ scales linearly with the power-law exponent of the pressure profile and is proportional to $\left(1+\epsilon\right)^{-1}$. Since $\mathrm{St}\propto\left(1+\epsilon\right)^{-1}$ as well, the value of the critical grain size does not depend on $\epsilon$.

4 Numerical simulations

4.1 Dust/gas simulations

We use the SPH code phantom to perform 3D global simulations of gas/dust discs containing an embedded protoplanet (Price et al., 2017). Importantly, every process involved in the physical problem (viscosity gravity and drag) is computed self-consistently. A calibrated non-zero viscosity is applied on each gas particle (Lodato & Price, 2010), to mimic the viscous transport of gas in disc described with an Prandtl-like model of turbulence (Shakura & Sunyaev, 1973). The star and the planet are treated as mobile point sources of mass. Their gravitational interactions with the gas/dust particles are computed using the sink particles approach (Bate, Bonnell & Price, 1995), which allows the bodies to migrate from their interactions with the disc. Dust and gas particles are accreted onto the sinks when two conditions are fulfilled: i) the SPH particle is found to be gravitationally bound the sink, and ii) the divergence of the velocity field at the location of the particle is negative. The gravitational acceleration between the $n_{\mathrm{sinks}}$ sink particles and the $i$-th gas/dust SPH particle is computed according to

[TABLE]

where $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$ denotes the differential location between the particles, $M_{j}$ is the mass of the $j$-the sink particle and $s_{j}$ is the usual softening parameter which prevents singularities at the sink locations. $s_{j}$ is also chosen to be the accretion radius of the sink particle.

The dust motion is computed using the two fluid algorithm described in Laibe & Price (2012). The drag force between a particle of one type and its neighbours of the other type is calculated in a pairwise manner, to conserve the linear and the angular momentum of both phases as well of the energy of the gas to machine precision. In particular, the drag from the dust onto the gas (sometimes referred as “back-reaction”) is included self-consistently. A specific double-hump drag kernel ensures the accuracy of the interpolation. To ensure better resolution within the gap, the smoothing length of the gas is used to compute drag terms. Outside of the gap, results do not depend from this choice. The algorithm has been extensively benchmarked on simple test problems including waves and shocks in dust and gas mixtures (Laibe & Price, 2011, 2012; Price & Laibe, 2015). We model spherical, compact and uncharged grains of constant sizes. The drag coefficient is computed consistently, based on the local values of the Knudsen number, as well as the Reynolds and the Mach number of the relative flow between the two phases (Kwok, 1975; Paardekooper & Mellema, 2006; Laibe & Price, 2012).

4.2 Initial conditions

The disc is setup in phantom by following the procedure outlined in Lodato & Price (2010). The system consists in a central star of mass $1.3\,M_{\odot}$ surrounded by a gaseous disc of $5\times 10^{5}$ SPH gas and $3\times 10^{5}$ SPH dust particles extending from $r_{\mathrm{in}}$ = 1 to $r_{\mathrm{out}}$ = 120 au. We model the initial surface density profiles of the discs using power-laws of the form $\Sigma(r)=\Sigma_{\rm in}(r/r_{\rm in})^{-p}$, where $\Sigma_{\rm in}$ is set such as the total gas mass contained between $r_{\rm in}$ and $r_{\rm out}$ is 0.0002 M*⊙* and a dust-to-gas ratio of 0.01. We adopt a power law exponent of the gas surface density profile $p=0.1$, and the aspect ratio of the disc is assumed to be $\sim 0.07$ at the planet location (corresponding to $0.04$ at $r_{\rm in}$). We assume a vertically isothermal equation of state $P=c_{\rm s}^{2}\rho_{\mathrm{g}}$ with $c_{\rm s}(r)=c_{{\rm s,in}}(r/r_{\rm in})^{-0.35}$. The exponent of the power-law profile of the pressure is therefore $\zeta=-1.95$. We set an SPH viscosity parameter $\alpha_{\rm AV}=0.1$ which ensures an effective Shakura & Sunyaev (1973) viscosity $\alpha_{\rm SS}\approx 0.004$. We study the evolution of the dust density resulting from the tides of one embedded planet located at a distance of 40 au from the central star. We perform a series of simulations varying the planet mass, the aspect ratio of the disc and the size of grains in order to test our criterion over a wide range of disc models.

4.3 Results

4.3.1 Planet mass

We simulate the evolution of 1 mm sized dust grains over 40 planetary orbits, which leaves enough time for the dust to settle and for the gap to form. The dust grains in our model have an initial Stokes number of $\mathrm{St}\sim 7$ at the disc midplane. Those large grains settle efficiently to the midplane of the disc in a stable dust layer with dust-to-gas scale height ratio of $\sim\sqrt{\alpha_{\rm SS}/\mathrm{St}}\sim 0.02$, consistent with the Dubrulle et al. (1995)’s model and SPH simulations of dusty discs (Laibe et al., 2008). A dust-to-gas ratio in density $\epsilon$ of $\sim 0.5$ is achieved in the midplane of the disc. However, as long as dust back-reaction remains weak enough to not affect the local pressure profile, our analysis shows that the minimum mass for dust gap opening and the location of the outer gap edge do not depend on the local dust-to-gas ratio. Therefore, we do not need to know the exact shape of the dust-to-gas ratio profile to determine the outer gap edge.

Our analytic criterion (Eq. 50) predicts for the minimum mass able to carve a gap in the dust disc to be $M_{\rm p}\simeq 0.052\,M_{\rm J}$ for grains with $\mathrm{St}>\mathrm{St_{crit}}=3.7$, assuming the nominal value of $f=0.4$ (Goldreich & Tremaine, 1979). Above this mass, we expect that planets of mass $\gtrsim M_{\rm p,lim}\simeq 0.13\,M_{\rm J}$ (Eq. 36) are expected to perturb the local pressure profile, weakening the gap closing effect induced by drag in the outer orbit (Rosotti et al., 2016). Moreover, we expect to see gaps only in the dust for planets up to $M_{\rm p,gap}\simeq 0.27\,M_{\rm J}$, according to our estimation of minimum mass to create a pressure maxima (Eq. 35). We therefore vary the planet mass in the range $0.01-0.25\,M_{\mathrm{J}}$ and look at the eventual structure of dust gap.

Fig. 3 shows rendered images of the dust surface density of the disc hosting a planet with different masses. A planet of mass $M_{\rm p}\lesssim 0.09M_{\mathrm{J}}$ causes a local depletion of grains in its close neighbourhood. This local void is permanently replenished by an incoming flux of drifting particles after the planet transit, leading to the formation of a non-axisymmetric depletion of dust. Fig. 4 shows that planets of mass $0.09M_{\mathrm{J}}\lesssim M_{\mathrm{p}}\lesssim 0.125M_{\mathrm{J}}$ are not able to disturb the local pressure profile. For $M_{\rm p}\gtrsim 0.13M_{\mathrm{J}}$, the pressure gradient is still always negative (bottom-right panel of Fig. 4), but the weakening of the local pressure profile produces a decrease in the radial drift velocity and therefore a deeper and wider gap compared to low-mass case. The surface density profiles of the dust shows that the width and the depth of the gap change sensitively at the transition between the two regimes (top left panel of Fig. 4). The minimum mass able to reshape the local pressure profile is $\sim 0.13\,M_{\mathrm{J}}$, in excellent agreement with the value of $M_{\rm p,lim}$ found by Rosotti et al. (2016). However, we obtain deeper dust gaps for $M_{\mathrm{p}}\gtrsim M_{\rm p,lim}$ compared to Rosotti et al. (2016). This discrepancy originates from a different choice for the smoothing length used to soften the tidal torque in the dust. They use a length of order the gas scale height, whereas we use a smaller length of order the dust scale height $H_{\mathrm{d}}$. Hence, we obtain a better estimate for the tidal torque in the dust in the region $\max(r_{\mathrm{H}},H_{\mathrm{d}})<\left|r-r_{\mathrm{p}}\right|<H$. In this region, tides dominate and shape the density distribution of large grains around low-mass planets, for which $r_{\mathrm{H}}<H$.

As discussed in Sect. 3.3.3, we use the value of the critical mass obtained by our simulations to determine the value of the constant $f$ in front of Eq. 51. We obtain

[TABLE]

in agreement with theoretical estimates (see Sect. 2.1). The low uncertainty of $0.01$ is due to the sensitivity of the critical mass with respect to $f$ (see Eq. 44). Hereafter, we adopt this sole value of $f$ to compare the results of all numerical simulation to our theoretical model. Using Eq. 54, the critical Stokes number above which we expect to observe a gap only in the dust is $\mathrm{St_{crit}}=5.3$.

To estimate the location of the gap outer edge from the dust surface density profile, we follow the approach of Dong & Fung (2017), which appears to work better for shallow gaps where the density does not drop below an empirical threshold. The gap outer edge $r_{\mathrm{gap}}$ is defined as the location outside the planet orbit where the dust surface density $\Sigma_{\mathrm{d}}(r_{\mathrm{gap}})$ reaches the geometric mean between its minimum value in the gap $\Sigma_{\mathrm{d}}(r_{\mathrm{min}})$ and its unperturbed value at the same location $\Sigma_{\mathrm{d,0}}(r_{\mathrm{min}})$, i.e.

[TABLE]

Fig. 5 shows that the position of the gap outer edge estimated from all the simulations is consistent with the value predicted by Eq. 39 for the entire range of masses where our analysis is valid, i.e. $0.09\,M_{\mathrm{J}}\lesssim M_{\rm p}\lesssim M_{\rm p,lim}$. For more massive planets, shallower pressure profiles reduce the drag torque outside the planet location, which translates the outer edge of the gap further away from the planet compared to our analytical predictions. In this case, our model underestimates the value found in numerical simulations by $\sim 20\%$. We also evaluate the location of the gap outer edge from the zero dust-velocity condition of Eq. 38, using consistent radial profiles of surface density and temperature across the gap obtained from numerical simulations (dashed line in Fig. 5). We note that the analytical approximation given by Eq. 39 works fairly well – the discrepancy being of order a few per cent – due to the low values of $\Delta_{\mathrm{gap}}/r_{\mathrm{p}}$. More massive planets create denser dust gaps outer edges. Since resolution in SPH increases with density, errors on the location of the gap decrease with the mass of the planet.

It is worth remarking that the limiting case where the outer gap edge is equal to the Hill radius (Fig. 5, dotted line) is obtained for a planet of mass $\sim 0.017M_{\mathrm{J}}$, where our simulations does not show a gap-like structure. Finally, in the low-mass planet regime, no gap forms in the gas. We test this condition by running a simulation where the gravity is switched off in the dust, to discriminate the mechanism at the origin of the dust gap opening. We find that the gap is indeed carved by the mechanism explained in this paper, as expected (see lower right panel in Fig. 1 of Dipierro et al. 2016).

4.3.2 Aspect ratio of the disc

Eq. 50 shows that dust gap opening depends on the ratio $H/r_{\mathrm{p}}$ since radial drift is triggered by the radial pressure gradient of the gas. We test the criterion derived above by performing a simulation with a disc satisfying $H/r_{\mathrm{p}}=0.035$ at the planet location, a value two times lower than the one previously adopted. From Eq. 50, the minimum mass of the planet able to carve a dust gap in this disc model is expected to be $0.011\,M_{\mathrm{J}}$. For planets more massive than $M_{\mathrm{p,lim}}=0.017\,M_{\mathrm{J}}$, the local pressure profile is perturbed (Eq. 35). We perform simulations using the higher adequate resolution to ensure the viscosity is the same as for the other models. We verify that in this simulation suit, the local pressure profile remains monotonic (see Fig. 7). Fig. 6 shows that a planet of mass \ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ 0.012\,M_{\mathrm{J}} ($3\,M_{\oplus}$) is able to carve a gap in the dust, a result consistent with our predictions. Fig. 8 shows that the location of the outer edge of the gap predicted by our model reproduces well the results of simulations. From Eq. 51, the minimum planet mass scales as $\left(H/r_{\mathrm{p}}\right)^{-3}$ and is therefore eight times less massive than the one in our model of reference. Moreover, Eq. 39 indicates that the outer edge of the gap carved by the planet with the minimum mass expressed in Eq. 51 scales with $\left(H/r_{\mathrm{p}}\right)^{2}$. Thus, the gap outer edge of the minimum mass planet is four times lower that the one in the reference model. This can be verified by comparing Figs. 5 and 8.

4.3.3 Grain sizes

When the size of dust grains increases, the drag torque weakens and gap closing caused by drag becomes less efficient (see the first term in Eq. 22). Hence, planets of very low masses can open dust gaps as long as locally, the grains are large enough. Moreover, we expect the outer edge to be further away from the planet for larger Stokes number (Sect. 3.3.3). To test these predictions, we perform a series of simulations using the disc model studied in Sect. 4.2, a planet of mass $0.1\,M_{\mathrm{J}}$, and varying the initial size of the grains. The simulations are evolved over 100 planetary orbits for the grains to relax in a steady state outside of the planet orbit. Fig. 9 and the left panel of Fig. 10 show that dust gaps of large grains present the asymmetric W-shape evidenced by Ayliffe et al. (2012) and Dipierro et al. (2016) around the orbit of the planet. In this region, the drag torque is too weak to prevent the formation of a large and stable population of dust grains in the corotation region. The left panel of Fig. 10 compares the location of the gap outer edge obtained in SPH simulations and the one derived from our analytic model. The agreement between the theory and numerical results is very good (5-10 %). The moderate errors between our theoretical estimate and numerical simulations are due to the peculiar shape of the dust gap.

4.4 Summary

We have considered a disc hosting a low-mass planet which does not disturb the local pressure profile of the gas. We obtained two analytic criteria for the minimum mass of the planet required to i) stop the inflow of dust particles (axisymmetric mechanism) and ii) ensure that drift can not refill the inner regions of the disc in dust (non-axisymmetric mechanism). These two criteria represent the necessary and sufficient conditions for dust gap opening, respectively. The exact value of the minimum masses predicted by these conditions depends on the proportionality constant in front of Eq. 1 (planet migration is neglected). By combining our various numerical tests, our final sufficient condition for gap opening in dusty discs is given by

[TABLE]

for grains with

[TABLE]

We predict the outer edge of the dust gap to be located at a distance $\Delta_{\mathrm{gap}}$ from the planet, where

[TABLE]

Planets with masses larger than the limit given by Eq. 56 but lower than $\sim M_{\mathrm{p,lim}}$ (Eq. 36) shall carve a deep gap in the dust without affecting the gas structure. More massive planet in the range $M_{\mathrm{p,lim}}\lesssim M_{\mathrm{p}}\lesssim M_{\mathrm{p,gap}}$ (Eq. 35) are expected to slightly perturb the local pressure profile, leading to the formation of a dust gap due to the combined action of tidal torque and the weakening of the drag. More massive planet, $M_{\mathrm{p}}\gtrsim M_{\mathrm{p,gap}}$, carve a gap both in the gas and dust phase (Lambrechts et al., 2014; Rosotti et al., 2016). Fig. 11 shows that numerical simulations corroborate the different dust gap opening regimes predicted in Sect. 3 for low-mass planets.

5 Using the criterion

5.1 Interpreting observations of gaps

5.1.1 Gap detectability

Relating the morphology of dust gaps to the properties of the planet and the disc gives insights about the planet formation process. Ideally, multiple-wavelengths observations should be combined to infer density distributions of grains experiencing different aerodynamical regimes. Scattered light emission at optical and near-infrared frequencies trace small dust grains ($\simeq 0.1-10~{}\mu$m) at the surface of the disc, where stellar photons are absorbed or scattered (Watson et al., 2007). The scattering emission intensity probes the gas structure at the surface of the disc, since these grains are efficiently coupled with the gas (\mathrm{St}\ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ \alpha). Emission at (sub)-millimetre wavelengths probes surface density of large grains in the mid-plane of the disc ($\simeq 0.1-10$ mm), since discs are usually optically thin at these wavelengths in the vertical direction (Dullemond et al., 2007; Williams & Cieza, 2011). A narrow beam is required to resolve the gap, together with a large signal-to-noise ratio to discriminate its weak emission. The gap depth can then be extrapolated, assuming that the weak emission in the gap is solely due to a low dust surface density.

5.1.2 Estimating the Stokes number

Eq. 56 involves the Stokes number of the grains. $\mathrm{St}$ can be estimated directly if the gas surface density is known. Unfortunately, hydrogen density is a quantity which is not directly measurable in a disc. Gas masses are therefore usually estimated by processing the (sub-)millimetre continuum or line measurements of CO isotopologues. This requires to model fractions of isotopologues abundances, dust grains opacities and local gas-to-dust mass ratios (Williams & Cieza, 2011; Williams & Best, 2014; Miotello et al., 2014). An estimate of the value of the Stokes number $\mathrm{St}$ can also be inferred indirectly via the ratio $\mathrm{St}/\alpha$ coming from the thickness of the dust layer (e.g. Dubrulle et al., 1995), assuming a fixed value for the turbulence parameter $\alpha$. Similarly, assuming a fixed value of the Stokes number, it is possible to infer the level of gas turbulence by analyzing the dust settling. As an example, Pinte et al. (2016) measured a dust scale height of $\sim 1$ au at $r\sim 100$ au in the disc around HL Tau. Assuming coupled grains (e.g. $\mathrm{St}\simeq 0.01$), this thickness implies $\alpha\simeq 10^{-4}$. However, i) this value is $1-2$ orders of magnitude smaller than the value consistent with typical accretion rates of protostars, ii) such low $\mathrm{St}$ implies typical Minimum Mass Solar Nebula gas density, which would make the planet migrate very quickly onto the star, and iii) when planets interact with discs of such low viscosity, vortices develop via the Rossby wave instabilities, trap grains and produce non-axisymmetric structures detectable by ALMA (Lyra & Lin, 2013), not detected in HL Tau. Assuming $\mathrm{St}\simeq 1$, the value of the measured thickness implies $\alpha\simeq 10^{-2}$, more consistent with the expected value of accretion rate of protostar.

5.1.3 Gap in the gas

The mass of the hypothetical planet and the local properties of the disc can be estimated when observing a gap in optical or NIR scattered light emission, given a degeneracy over the ratio $M_{\rm p}/\alpha$ (Fung et al., 2014; Kanagawa et al., 2016; Rosotti et al., 2016). If the signal-to-noise ratio is high enough, the depth of the gap (Fung et al., 2014; Kanagawa et al., 2015) or its shape (Kanagawa et al., 2015, 2016) can also be used to the same purpose. These methods suffer large systematic errors, mostly due to uncertainties on the local disc geometries which heavily affect the surface brightness around the gap (Jang-Condell & Turner, 2012). Scattering emission may additionally reveal spiral structures, whose morphologies may be related to the mass of the planet and the aspect ratio of the disc (Zhu et al., 2015; Juhász et al., 2015; Dong et al., 2015).

An additional method to probe gaps in the gas is the detection of line emissions of CO isotopologues, such as 12CO, 13CO and C18O. Those might be optically thin at the corresponding wavelengths and trace the gas down to the disc midplane (Williams & Best, 2014; Miotello et al., 2014). Isella et al. (2016) claim evidence of a decrement in the density of CO isotopologues within the middle and outer continuum gaps in the disc around HD163296. Yet, the decreased emission of the CO molecular lines might also be produced by a reduced density of large grains. In detail, since photodissociation by ultraviolet radiation (UV) is the primary process that regulates the abundance of gas phase CO in the emitting layer of discs, a reduced dust density around the planet location might induce a less efficient absorption of the UV photons. The decreased shielding of the CO molecules by dust lead the UV photons to penetrate into the disc and become optically thick at higher column densities, i.e. closer to the midplane. The higher efficiency of UV photodissociation at the planet orbit with respect to the adjacent regions might therefore produce a decreased emission of CO molecular lines, that can be misinterpreted as a real gas gap. However, since small grains are much more efficient in absorbing UV radiation, a gap in only large dust grain might not affect remarkably the shielding of CO molecules since larger grains have less opacity in the UV and do not shield CO strongly (Visser et al., 2009). Connecting the variation of CO isotopologues emission lines with real gas density variations still remains an open question.

5.1.4 Applying the criterion

A good use of the criterion starts with two preliminary remarks. Firstly, Fig. 2 shows that M_{\mathrm{p}}\ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ M_{\mathrm{p,gap}} is the only condition required for planets to open gaps in both the gas and the dust and this, independently of the Stokes number. Thus, if a gap is detected in NIR scattering and thermal-mm emission, no information can be extracted on $\mathrm{St}$ with the criterion derived in this study. Indeed, grains of all sizes tend to drift towards the pressure maximum at the outer edge of the gas gap, as long as $\mathrm{St}>\alpha$. Further analysis of the morphological details of the gap should be conducted to infer the properties of the system. Secondly, the absence of any gap in both gas and dust does not necessarily reflect the absence of any gravitational body in the disc. Grains may replenish the orbit of the planet as they drift inwards if $\mathrm{St}\lesssim\mathrm{St_{crit}}$. However, if any additional detection limit is given for the maximum planet mass a disc can embed, Eq. 57 provides a condition on the minimum Stokes number compatible with the eventual existence of the planet.

We now focus on the non-trivial case, i.e. a gap detected only at millimetre wavelengths. This work shows that any low mass planet can create this structure as long as the Stokes number is large enough. As an example, for $\mathrm{St}\simeq 20$, the necessary mass of the planet required to open the dust gap is $\simeq 100$ lower than the one required to open a gap in the gas due to the factor $\mathrm{St}^{3/2}$ in Eq. 56. We therefore expect observations of dust only gaps to be more frequent in the outer disc, where grains of a given sizes have Stokes numbers much larger than unity. At least, $\mathrm{St}\gtrsim 1$, which constrains the maximum local density of the gas. Eq. 58 gives the expression of the distance between the planet mass and the outer edge of the dust gap. Even if the Stokes number is only roughly approximated, the weak sensitivity brought by the factor $\mathrm{St}^{1/4}$ allows to determine the planet mass relatively precisely (assuming that the aspect ratio of the disc is known). The absence of any gap in the gas additionally implies that $M_{\mathrm{p}}<M_{\mathrm{p,lim}}$. Combined with the roughly known value of $M_{\rm p}$, this condition provides a minimum value for $\alpha$ in the disc. Comparing this value and the degree of dust settling may give a way to infer if the seminal diffusive description of turbulence in disc is relevant or not. Interestingly, if several dust only gaps are detected in the same disc, the degeneracy over the constant in front of the criterion can be broken, helping to determine the masses of the planets more precisely.

5.2 Limitations

In this study, we have restricted our analysis to planets on fixed orbits. Migration may strongly affect the ability of a planet to carve gaps in the gas phase (Malik et al., 2015) and can appreciably change the density structure around the planet (Rafikov, 2002b). On the contrary, Fig. 12 shows that the dust gap opening time scale obtained from Eq. 48 is much shorter than the migration time scale, estimated from the differential Lindblad torque derived in Tanaka et al. (2002) in our reference disc model with an embedded planet of mass $0.1M_{\mathrm{J}}$. Hence, planet migration is not expected to affect the ability of a planet to carve gaps in our mechanism. Note that the ratio between both timescales depends indirectly on the Stokes number through the local gas density.

We have also assumed dust grains of constant sizes. It is known that dust coagulation or fragmentation may strongly affect the dust dynamics (e.g. Laibe et al., 2008; Birnstiel et al., 2010). However, we can safely neglect the grains size evolution over the small time required to open the dust gap. Finally, we have assumed that back-reaction is weak enough and does not affect the gas surface density significantly as grains are repelled outside of the planet orbit. However, large grains with $\mathrm{St}>1$ are not able to affect significantly the gas structure (see Eq. 15). Whether back-reaction can trigger the formation of gaps in the gas would be worth investigating, but remains beyond the scope of this study. Hence, the criterion has not been proven to work at dust-to-gas ratio larger than unity.

In Sect. 2.1 we rationalised the use of Eq. 1 to model the tidal torque on dust by advocating that gaps open in a region where high-order Lindblad resonances are highly concentrated and degenerate into a continuum, i.e. $r_{\mathrm{H}}\lesssim\Delta_{\mathrm{gap}}\lesssim H$. We test this assumption by comparing the value of $\Delta_{\mathrm{gap}}$ predicted by our theory with the scale height of the disc, for planet masses within the range $\left[M_{\mathrm{p,crit}},M_{\mathrm{p,lim}}\right]$ and various Stokes number. We first note that

[TABLE]

In both cases, the outer edge of the gap increases linearly with the aspect ratio of the disc. Fig. 13 shows the range of theoretical locations of the gap outer edge as a function of the Stokes numbers. For $M_{\mathrm{p}}=M_{\mathrm{p,lim}}$, the gap outer edge is smaller than the local gas scale height if

[TABLE]

For a typical disc with $H/r_{\mathrm{p}}=0.05$ and $\zeta=-2.75$, the Stokes number above which the gap is larger than the disc scale height is $\sim 55$. Hence, the location of the outer gap edge is smaller the local gas scale height for a large range of disc model parameters, which supports our initial assumption.

6 Conclusion

We derived an analytic criterion that predicts the minimum mass required for a planet to open a gap in the dust phase of a viscous protoplanetary disc in the case where the planet do not perturb the local pressure profile of the disc. In this regime, a gap opens in the dust if the tidal torque overpowers the drag torque outside the planet orbit. We generalised the approach of NSH86 to include the disc-planet tidal interaction and the viscous forces in the equations of motion. Gas and dust velocities in steady state were obtained analytically (Eq. 15-18). From there, assuming that the planet is not able to affect the local pressure structure, we derived a relation between the minimum mass required to open a gap in the dust, and the key parameters of the dust motion: the Stokes number, the aspect ratio of the disc and the dust-to-gas ratio (Eq. 51). We benchmarked the value of the scaling constant in front of the tidal torque density formula using 3D dust-and-gas SPH simulations of various discs.

Our final opening criterion for dust gaps is given by Eq. 56. We found that low mass planets are able to carve dust gaps when the Stokes number $\mathrm{St}\geq\mathrm{St_{crit}}\simeq 1$ (Eq. 57). We also derived an analytic formula for the radial extension of the outer dust gap edge (Eq. 58). The criterion and the location of the gap outer edge estimated by our analysis have been tested through 3D SPH simulations of a variety of dusty disc models with an embedded planet. The numerical results appear consistent with our analysis.

Observations in near-infrared scattering or millimetre thermal dust emission might reveal gaps, in both phases or in the dust only. This criterion can be used to constrain the mass of the planet embedded in the disc, or the Stokes number of the grains.

Acknowledgements

GD thanks Monash for CPU time on NeCTaR. We wish to thank Giuseppe Lodato, Daniel J. Price and Richard Alexander for fruitful discussions and the referee for an insightful report of the manuscript. GD acknowledge funding via PRIN MIUR 2010-2011 ‘Chemical and dynamical evolution of the Milky Way and Local Group Galaxies’. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 681601). GL is funded by European Research Council FP7 advanced grant ECOGAL. We used SPLASH (Price, 2007).