© ESO, 2014

1. Introduction

Planet formation occurs in the gaseous and dusty protoplanetary disks that remain around young stars during their formation. The wide diversity of exoplanets detected in recent decades (e.g. Batalha et al. 2013) may have its origin in the different gas and dust distributions observed in protoplanetary disks (e.g. Mordasini et al. 2012). The understanding of how protoplanetary disks disperse and form planets is a rapidly evolving subject and many questions remain open (e.g. Alexander et al. 2014). Circumstellar disks with little or no excess emission at the near/mid-infrared, but with significant excess at longer wavelengths were first detected in 1989 (Strom et al. 1989), indicating disks with dust depleted inner regions. These inner cavities have more recently been directly imaged in the sub-millimetre regime (e.g. Brown et al. 2009).

These so-called transition disks are believed to be an intermediate step of the evolutionary phase between gas/dust-rich disks and debris disks, however, whether or not transition disks are the result of undergoing disk dispersal is still a subject of discussion (e.g. Espaillat et al. 2014). Different mechanisms have been proposed to explain the observed properties of transition disks such as protoevaporation (e.g. Alexander et al. 2006; Owen & Clarke 2012), interaction with one or multiple planets (e.g. Varnière et al. 2006; Dodson-Robinson & Salyk 2011), and dust growth (Dullemond & Dominik 2005; Birnstiel et al. 2012). Additional studies have been done assuming that different processes simultaneously take place in these disks (e.g. Rosotti et al. 2013), however, a full understanding of the physical properties of transition disks is far from complete.

Fig. 1 Sketch of the spatial segregation between small and large grains expected when two planets interact with the disk and lead to particle traps. Asymmetric structures may exist due to azimuthal trapping in vortices or/and disk eccentricity.

Any clearing mechanism that creates a gas surface density depletion leads to a decoupling of small (micron-sized particles) and large grains (millimetre-sized particles; Rice et al. 2006; Pinilla et al. 2012; Gonzalez et al. 2012; Zhu et al. 2012, which may result in dissimilar observed structures at different wavelengths. Observations of transition disks by the Strategic Explorations of Exoplanets and Disks with Subaru survey (SEEDS) show that for some cases, a cavity is not detected at near-infrared scattered light emission, contrary to millimetre observations (Dong et al. 2012). Molecular gas and polycyclic aromatic hydrocarbon (PAH) emission have also been detected inside millimetre dust cavities (Geers et al. 2007; Salyk et al. 2009; Pontoppidan et al. 2008; Maaskant et al. 2013). This “missing cavities” problem could be the result of considering the evolution of dust when a massive planet is interacting with a disk (e.g. Rice et al. 2006; de Juan Ovelar et al. 2013). In these models, the radial distribution of micron-sized particles is expected to be similar to the gas distribution, while a large radial separation between the gas and dust is expected for millimetre particles (Pinilla et al. 2012). This interesting radial disentanglement has already been observed for different transition disks, such as HD 135344B (Garufi et al. 2013), HD 100546 (Ardila et al. 2007; Walsh et al. 2014), and SR 21 (Follette et al. 2013; Pérez et al. 2014). Moreover, observations show interesting structures of transition disks, such as spiral arms at scattered light (e.g. Muto et al. 2012; Grady et al. 2013), eccentric gaps (e.g. Thalmann et al. 2014), and strong azimuthal asymmetries at the millimetre emission (Casassus et al. 2013; van der Marel et al. 2013; Fukagawa et al. 2013), which may also result from planet disk-interaction(s).

The knowledge of gas and dust density distributions at different wavelengths may hint at the physical properties of potential embedded planet(s) in the disk. In this paper, we aim to study the observational signatures in the case of two planets interacting with a disk, by analysing the resulting images and spectral energy distributions (SEDs) obtained by combining hydrodynamical, dust evolution, and radiative transfer numerical simulations. The main questions to address are illustrated in Fig. 1:

• Considering different disk and planet properties, how efficient is dust trapping, and where are the dust traps located?

• What is the resulting dust size distribution across the disk? Can it lead to spatial segregation between micron- and mm-grains?

• What are the conditions to induce an asymmetry in the disk? Is it due to vortex formation or disk eccentricity?

Thanks to the wealth of available multi-wavelength observations, the disk around SR 21 is an excellent candidate to apply these models. The star SR 21 is located in the Ophiuchus star-forming region, whose disk was identified as a transition disk by Brown et al. (2007). They reported Infrared Spectrograph (IRS) data of the Spitzer Space Telescope, and inferred a 0.45–18 AU gap by modelling the SED. Observations of CO ro-vibrational transition with CRIRES on the Very Large Telescope (VLT) indicates a narrow ring of gas at ~7 AU distance from the star (Pontoppidan et al. 2008), with a width of ~1 AU. Follette et al. (2013) present scattered light imaging obtained with Subaru HiCIAO camera in H band (1.6 μm), finding that large (~mm) and small (~μm) particles are decoupled. They conclude that the polarised intensity radial profile of their observations is decreasing steeply outwards (r^-3), with no break at the location of the cavity radii observed at millimetre wavelength (~36 AU, Andrews et al. 2011). The sub-millimetre observations at 350 GHz (850 μm) by Brown et al. (2009) and Andrews et al. (2011) suggest an asymmetric feature, which was confirmed with ALMA observations at 690 GHz (450 μm) (Pérez et al. 2014) with an emission contrast of around two. To model this asymmetry, Pérez et al. (2014) fit a 2D Gaussian at ~46 AU from the central star. In the following, we will consider the stellar parameters of SR 21, and aim to reproduce the main characteristics of these observations.

The structure of this paper is: in Sect. 2, hydrodynamical, dust evolution, and radiative transfer numerical models are described, with the corresponding set-up. The results of the simulations are shown in Sect. 3. We present synthetic images and SEDs in Sect. 4. Sections 5 and 6 are the discussion and main conclusions of this work.

2. Models and setups

In this section, we describe the codes and setups used for the hydrodynamical, dust evolution, and radiative transfer simulations. Our model parameters are chosen to reproduce transition disks with large millimetre cavities and a narrow ring of material inside the cavity. For this, we assume two non-migrating planets embedded in the disk.

2.1. Hydrodynamical simulations

We explore the evolution of the gas surface density of a disk hosting two non-migrating planets, by considering the gravitational effect of the planets onto the disk. For these simulations, we used the 2D version of the fast advection hydrodynamical code FARGO (Masset 2000). In this code, the length and mass units are scaled by the planet location and the mass of the central star, respectively. In this case, we use FARGO to solve the Navier-Stokes and continuity equations for a flaring disk (H ≡ h/r = c_s/v_k, with c_s the sound speed and v_k the Keplerian velocity) interacting with two planets. The aspect ratio of the disk is increasing with radius as H = h₀r^f, with a flaring index f of 0.25 and an aspect ratio h₀ of 0.05 at the position of the inner planet.

For vertically isothermal laminar disks, a ~1 M_Jup mass planet can clear a gap (e.g. Lin & Papaloizou 1993). For these simulations, we consider the planet-to-central-stellar mass ratio for the inner planet (P1) equal to q_P1 = 4 × 10^-4, which corresponds, for example, to a 1 M_Jup planet around a 2.5 M_⊙ star. For the outer planet (P2), the mass ratio is q_P2 = [ 2 × 10^-3,4 × 10^-3,6 × 10^-3 ], which for a stellar mass of 2.5 M_⊙ implies planet masses of 5,10 and 15 M_Jup, respectively. These masses are taken because large separations between the distribution of gas and millimetre dust particles are expected for such massive planets embedded in the disk (Pinilla et al. 2012) and long-lived asymmetries due to eccentricity or vortex formation may exist at the edge of the gap (e.g. Ataiee et al. 2013). The dimensionless radial grid was taken to be logarithmically spaced from 0.1 × r_P1 to 24 × r_P1, and the resolution was n_r × n_φ = 512 × 587. We used non-reflecting conditions for the inner boundary and two values for the disk viscosity were taken for each case α_turb = [ 10^-3,10^-2 ]. The gas surface density is initially a power-law function, Σ = Σ₀ × (r/r_P1)^-1, and Σ₀ is taken such that the disk mass is 2.5 × 10^-3M_⋆.

The two planets are considered to be in fixed circular orbits, such that r_P2/r_P1 = 3.5, and they do not feel the gravity of the disk nor each other. We adopt the position of the inner planet to be at 5 AU, hence the outer at 17.5 AU. The motivation for these locations are based on the SR 21 features (although additional tests were done varying the location of these planets, finding similar conclusions). When a planet opens a gap in a disk, the outer gap edge is expected to be at most 5r_H for the gas and ~7−10r_H for mm-sized particles (Pinilla et al. 2012), where the Hill radius r_H is r_H = r_p(q/ 3)^{1 / 3} and q the planet-stellar mass ratio. Thus, under the assumptions of our simulations, a gas ring at the outer edge of the gap carved by the inner planet is expected at ~0.26 × r_P1 from the planet (i.e. ~6.3 AU from the star). In the outer regions, when particles are trapped at the edge of the second gap, the peak of the mm-emission is expected to be at ~[ 0.61,0.77,0.88 ] × r_P2 for q_P2 = [ 2 × 10^-3,4 × 10^-3,6 × 10^-3 ], respectively i.e. ~[ 28,31,33 ] AU from the star.

Fig. 2 Left panels: 2D gas surface density after 1000 orbits of the inner planet, when two planets with rP2/rP1 = 3.5 are interacting with the disk, with a planet-to-stellar mass ratio of qP1 = 4 × 10-4 for the inner planet and qP2 = 2 × 10-3 (left), qP2 = 4 × 10-3 (middle), and qP2 = 6 × 10-3 (right), and two different values of disk viscosity αturb = 10-2 (top) and αturb = 10-3 (bottom). Right panels: azimuthally averaged gas surface density from the corresponding hydrodynamical simulations and time-averaged over 200 orbits (from 1000 to 1200 orbits).

2.2. Dust evolution models

We adopted the dust evolution model explained in Birnstiel et al. (2010). This is a numerical 1D grid-based model for the global evolution of the dust surface density as introduced by Brauer et al. (2008). In this case, we consider that the dust is transported because of drag forces and mixing by turbulence. A grid of 180 particle sizes is taken, covering grain sizes from ~1 μm to ~200 cm. For the particle growth, sticking, erosion, and destructive collisional outcomes are considered. The threshold for fragmentation occurs when particles reach velocities of v_f = [ 10,30 ] m s^-1. Results from numerical simulations of collision between particles outside the snow line, where grains may have an ice mantle, indicate that the sticking efficiency is high and that fragmentation velocities may be as high as ~50 m s^-1 (Wada et al. 2009). We study the influence of planets on the dust evolution by considering the gas surface density and velocity from the azimuthally averaged values from the hydrodynamical simulations time-averaged over 200 inner orbits (between ~1000–1200 orbits of the inner planet) as the initial conditions for the dust evolution. Dust radial drift velocities are proportional to the pressure gradient, which depends on the gas surface density, and the drag velocities depend on the gas velocities. These two components of the dust radial velocity and turbulent motion also depend on the coupling of the dust particles to the gas i.e. the dust particle size and gas surface density. By taking the azimuthally averaged values, we do not consider the disk eccentricity and its effect on the final dust density distribution. In our simulations, we consider that the planets are not accreting gas or dust. For all of the cases, the disk mass and the initial gas-to-dust ratio are assumed to be M_disk ~ 0.006 M_⊙ (as Andrews et al. 2011, for SR 21) and 100, respectively.

2.3. Radiative transfer

To generate the SEDs and images at different wavelengths, we used the 2D Monte-Carlo radiative transfer code MCMax (Min et al. 2009). The vertical structure is computed considering the particle size distribution from the dust evolution models, and assuming settling and vertical turbulent mixing consistently with α_turb = [ 10^-3,10^-2 ] as for the hydrodynamical and dust evolution models (see Mulders & Dominik 2012, Eqs. (2) and (3)). For the grain composition, we considered a mixture of silicates (~58%), iron sulphide (~18%), and carbonaceous (~24%) grains as in de Juan Ovelar et al. (2013). For the inclination and distance to the source, we assume i = 15° and d = 135 pc. All parameters are summarised in Table 1.

3. Results of numerical simulations

In this section, we present the main results from the hydrodynamical simulations, together with the dust evolution models.

3.1. Planet-disk interactions

The left panels of Fig. 2 show 2D gas surface density after 1000 orbits of the inner planet (hereafter used as the time unit for the hydrodynamical results) when two planets interact with the disk. When the planet-to-stellar mass ratio of the outer planet increases, the corresponding gap is wider and deeper because of the higher angular momentum transfer, as shown by several authors (see Kley & Nelson 2012, for an extensive review). The effect of lowering the viscous torque, by assuming a lower value for α_turb, also has an effect on the gap shape (e.g. Crida et al. 2006). This is particularly true in the outer edge of the second gap, where a vortex is formed as in the case of a massive single planet interacting with a disk of moderate viscosity (e.g. Ataiee et al. 2013; Fu et al. 2014). The right panels of Fig. 2 illustrate the azimuthally averaged gas surface density from the corresponding hydrodynamical simulations, time averaged over 200 orbits (between 1000–1200 orbits). This explicitly shows the effect of the planet mass and viscosity on the gap shape. These gas surface density profiles are the initial conditions for the dust models and we assume that these time-averaged profiles remain constant during the dust evolution.

Fig. 3 Azimuthally averaged vortensity ⟨ η2D ⟩ φ (blue lines) and disk eccentricity (red lines) versus the disk radial extension for the hydrodynamical simulations of Fig. 2 after 1000 (dotted lines) and 2000 (solid lines) orbits. The vertical lines correspond to the positions of the two planets.

Fig. 4 Dust surface density for micron- (upper) and millimetre-sized particles (bottom), assuming two planets interacting with the disk at 5 AU and 17.5 AU and after 1 Myr of dust evolution. Two different thresholds are considered for the fragmentation velocity of particles vf = 10 m s-1 (left) and vf = 30 m s-1 (right). The vertical lines correspond to the positions of the two planets.

The interaction of a planet(s) with a disk induces gas asymmetries, which may be the result of disk eccentricity (e.g. Papaloizou et al. 2001; D’Angelo et al. 2006; Kley & Dirksen 2006; Hosseinbor et al. 2007; Lubow 2010), vortex formation due to Rossby wave instability (RWI; Li et al. e.g. 2000, 2005; Lin e.g. 2012, 2014; Meheut et al. e.g. 2012), or a mix of both. The Rossby wave instability develops at the local minima of potential vorticity or vortensity (η), triggering anticyclonic vortices. The vorticity is the tendency of a fluid to locally rotate around a central point, thus for a fluid with u velocity, vorticity is defined as ∇ × u, and the vortensity is the vorticity weighted by the density. If the minimum of η is deep, the RWI is strong and therefore the resulting vortex is also strong. In 2D, vortensity η_2D is defined as the ratio of the vertical component of the vorticity i.e. ẑ·∇ × u to surface density (e.g. Lovelace et al. 1999; Lin 2014).

Figure 3 illustrates the azimuthally averaged vortensity and disk eccentricity (calculated as in Kley & Dirksen 2006) versus the disk radial extension for the hydrodynamical simulations of Fig. 2, after 1000 and 2000 orbits. Both vortensity and eccentricity increase with higher planet-to-stellar mass ratio as well as for lower disk viscosity (α_turb). Nonetheless, the vortensity is almost flat at the edges of these gaps for any case of α_turb = 10^-2, implying no long-lived anticyclonic vortices. For α_turb = 10^-3, there are local vortensity minima at the outer edge of the gap carved by the second planet when q_P2 = 4 × 10^-3 and q_P2 = 6 × 10^-3, implying potential RWI that may trigger the formation of anticyclonic vortices, which can survive over thousands of orbits (Fig. 2). For these two cases, the azimuthal asymmetries are a combination of eccentricity and vortex formation. After 2000 orbits, the maximum values of the eccentricity for q_P2 = [ 2 × 10^-3,4 × 10^-3,6 × 10^-3 ] are e_max ~ [ 0.05,0.08,0.11 ] for α_turb = 10^-2 and e_max ~ [ 0.12,0.18,0.20 ] for α_turb = 10^-3.

3.2. Dust density distribution

In Fig. 4, we show the dust surface density for micron- and millimetre-sized particles after 1 Myr of dust evolution, assuming the gas surface density from the hydrodynamical simulations as the initial condition. Two different thresholds are considered for the fragmentation velocity of particles v_f = 10 m s^-1 and v_f = 30 m s^-1. The distribution of small and large grains is strongly affected by the shape of the gap i.e. the initial gas surface density, the turbulent parameter α_turb, and the fragmentation velocity v_f.

Because of the enhancement of the gas surface density at the outer edges of the gaps, there are pressure maxima and therefore preferential regions for particles to drift. Nevertheless, turbulence can push particles out from these pressure bumps. The balance between the turbulent mixing and the strength of the positive pressure gradient at these locations determines if grains can be trapped or not. Indeed, in the regions where the radial drift is reduced, the relative velocities between particles are dominated by turbulent motion, and the maximum particle size before they fragment (a_max) is parametrised as in Birnstiel et al. (2009)${\mathit{a}}_{\max }\mathrm{\propto }\frac{{\mathrm{\Sigma }}_{\mathrm{g}}}{{\mathit{\alpha }}_{\mathrm{turb}}{\mathit{\rho }}_{\mathrm{s}}}\frac{{\mathit{v}}_{\mathrm{f}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathit{,}$(1)where Σ_g is the gas surface density and ρ_s the volume density of a grain (of the order of 1g cm^-3). Thus, if the maximum particle size in pressure bumps results in bodies that are more affected by the turbulent mixing than radial drift, particle trapping would not happen.

Fig. 5 Dust density distribution at 1 Myr of evolution when two planets are interacting with a disk at 5 AU (1 MJup) and at 17.5 AU (5 MJup) with αturb = 10-2 (top), vf = 10 m s-1 (left) and αturb = 10-3 (bottom), vf = 30 m s-1 (right). The solid line corresponds to the particle size for which particles experience the highest radial drift, which is directly proportional to Σg (Eq. (1)).

The simulations with α_turb = 10^-2 and v_f = 10 m s^-1, show how turbulence prevails over the potential trapping of particles and there is no dust accumulation at the outer edge of the gap carved by the second planet, independent of M_P2 (Fig. 4). In this case, the increment of the dust surface density for small particles at the location of the planets happens because, in the outer regions, large grains fragment and the small dust co-moves with gas radial velocities without being trapped. The inward motion of dust particles from the outer edge through the gap creates a pile-up of small particles close to the planet location because of fragmentation, where the maximum size of particles (a_max) is lower than outside the gap (because of the lower gas surface density, see Eq. (1)). When the fragmentation velocity is higher, particles can reach larger sizes (Eq. (1)), leading to trapping of dust particles at pressure maxima, except for the case of M_P2 = 5 M_Jup, where the positive pressure gradient is not high enough to stop the rapid inward drift of particles. For these high fragmentation velocities (v_f = 30 m s^-1) and viscosity (α_turb = 10^-2), micron-sized particles are distributed all over the disk, with an enhancement at the location of the outer trap due to the continous fragmentation of large grains.

In contrast, in the case of α_turb = 10^-3, there is effective trapping of millimetre-grain particles in both gap outer edges, and independent of the fragmentation velocities. When v_f = 10 m s^-1, the constant fragmentation, by turbulent motion of large grains enriches the disk of smaller grains. In particular, this fragmentation occurs at the location of the pressure traps, which also leads to an increment of the surface density of micron-sized particles (see Andrews et al. 2014). When α_turb = 10^-3 and v_f = 30 m s^-1, fragmentation is unlikely to happen, and there are fewer micron-sized particles in the disk. For these cases, the dust surface density of mm-grains (Fig. 4) is dominated by the very large grains i.e. a> 10 mm.

The location and shape of the corresponding millimetre rings depend on the planet mass and viscosity. For a more massive outer planet, the mm dust particles accumulate further away (Pinilla et al. 2012). The width is affected by the values taken for α_turb, leading to a wider ring when turbulent mixing is higher. The amount of small and large grains in the inner disk (r< 5 AU), also depends on the turbulence strength and potential trapping. For instance, by comparing the cases of v_f = 10 m s^-1 and the different values of α_turb, when trapping occurs (α_turb = 10^-3), there is a much smaller amount of micron-sized particles outside the pressure traps (the same happens for the same α_turb, but with higher fragmentation velocity).

Figure 5 illustrates in more detail the dust distribution of all grain sizes and how dust grains are trapped or not for each condition, for the particular case of q_P2 = 2 × 10^-3 (i.e. 5 M_Jup around a 2.5 M_⊙ star). When fragmentation does not occur (as in the case of v_f = 30 m s^-1 and α_turb = 10^-3), particles continue growing to the largest size considered in the simulations (200 cm), leading to an optically thin disk at all wavelengths.

Fig. 6 Spectral energy distribution resulting from the radiative transfer modelling, assuming grain size distribution after 1 Myr of evolution when two planets interact with the disk, and assuming different values for the disk viscosity αturb and fragmentation velocity of particles vf. The disk data are for SR 21 and are taken from Brown et al. (2007); Follette et al. (2013) and PACS (DIGIT PI: N. Evans).

From the performed simulations, three cases allow radial trapping of particles at the outer edges of the carved gaps, a smooth radial distribution of micron-sized particles i.e. no total filtration of particles, and an effective fragmentation of dust grains i.e. no continuous growth to planetesimal-sized objects. These cases are: α_turb = 10^-3, v_f = 10 m s^-1 and M_P2 = 5 M_Jup, and α_turb = 10^-2, v_f = 30 m s^-1 and M_P2 = { 10,15 } M_Jup. Tables 2 and 3 summarise the main findings of the simulations done for this paper. To ensure the trends described in this paper, we repeated some of the simulations considering higher disk mass (by a factor of 2), larger distance between planets (r_P2/r_P1 = 4.5), and smaller initial grain size (0.1 μm), and found similar results.

In our models, we assume that the planets do not migrate and that they do not accrete gas or dust. Migration can create resonances, affect the eccentricity, and change the location of the pressure traps, which can affect the final distributions of dust. For instance, the inner planet partially opens a gap (Fig. 2), and therefore it can experience a rapid migration, inhibiting the particle trapping at the outer edge of its gap. The masses considered for the outer planet are high enough for a slow migration, allowing particles to be trapped at the outer edges of the corresponding gaps. These particle concentrations would move while the planet is migrating. When planets do not migrate, resonances between the planets are not expected because their separation is large. In this case, the maximum value found for the eccentricity between the two planets is ~0.15 (Fig. 3). This eccentricity is not high enough to overlap the gaps, implying that for all cases there is a long-lived gas ring between planets that, depending on different parameters (Tables 2 and 3), may or may not trap particles. Moreover, accretion of gas and dust onto the planet may also change the gas and dust distributions close to the planet orbits (e.g. Owen 2014).

4. Comparison with observations

In this section, we present the SEDs and synthetic sub-millimetre and polarised infrared images computed from the models.

4.1. Spectral energy distributions

Figure 6 shows the SEDs from radiative transfer simulations for all the cases introduced in Sect. 3, when two planets interact with the disk, different values for the disk viscosity α_turb, and fragmentation velocity of particles v_f. We consider the grain size distribution from the dust evolution models after 1 Myr of evolution. For comparison, in Fig. 6, we include data from Follette et al. (2013), the PACS (DIGIT data) and IRS spectra (Brown et al. 2007) for the case of SR 21 disk.

Fig. 7 Synthetic H-band HiCIAO (left panels) and ALMA Band 7 (850 μm; right panels) observations for the models of an outer planet MP2 = 5 MJup, and each value of the fragmentation velocity vf and turbulence parameter αturb. The white dashed line in the left panels indicates the inner of HiCIAO images masked in the Follette et al. (2013) observations of SR 21. In the right panels, the size of the ALMA Cycle 2 beam is indicated by a stripped circle.

With a fragmentation velocity of v_f = 10 m s^-1 and α_turb = 10^-2, there is no effective trapping of particles in the outer pressure maxima for any M_P2 (Fig. 4) and no dust depleted cavity, thus there is no characteristic deficit of emission at near/mid-infrared as seen in transition disks. In addition, the emission at far-infrared and mm-wavelength is under-predicted. When α_turb = 10^-3 and trapping of particles is possible, the mid-infrared emission decreases compared to the cases of α_turb = 10^-2 because of the cleared cavity. In the cases of α_turb = 10^-3, the planet mass slightly influences the SEDs. For instance, when M_P2 = 5 M_Jup, where there is less dust filtration at the corresponding gap (and therefore a larger amount of micron-sized particles at these locations, but less in the inner disks), there is lower near-infrared emission, but larger mid-infrared flux. Between M_P2 = 10 M_Jup and M_P2 = 15 M_Jup, there is almost no difference in the SEDs profiles because dust filtering is similar.

When v_f = 30 m s^-1, α_turb = 10^-2, and M_P2 = 5 M_Jup, there is no trapping and not enough dust to produce mid- and far-infrared emission (Fig. 4). Models with M_P2 = 10 M_Jup and M_P2 = 15 M_Jup allow many large grains in the outer part of the disk and therefore there is a higher far-infrared and millimetre emission. However, there are also a large amount of dust particles within the cavity, so no strong deficiency of the mid-infrared emission is generated. For these cases, the fact of having higher turbulence helps to increase the scale-height of the disk and the amount of dust in the upper layers, implying a larger amount of reprocessed stellar light. In the case of no fragmentation (v_f = 30 m s^-1 and α_turb = 10^-3), particles grow to very large sizes, causing the disk to be optically thin at most wavelengths and producing a debris-disk type SED.

The best models found in Sect. 3, which allow trapping of millimetre particles, effective fragmentation, and continuous distribution of micron-sized particles, are suitable for reproducing the characteristic SEDs of transition disks. The models of α_turb = 10^-3, v_f = 10 m s^-1, and M_P2 = 5 M_Jup may under-predict the typical far-infrared and mm-emission. Assuming a more massive disk or/and higher dust-to-gas ratio could increase the emission at mm-wavelengths. Some of the simulations were repeated with a higher disk mass (by a factor of 2), without significantly changing the resulting trends of the SEDs. In this case, the lack of far-infrared flux may be a consequence of a small disk scale height at the outer disk rim, where h(r) could be larger than in the inner disk (e.g. Maaskant et al. 2013). On the other hand, the models of α_turb = 10^-2, v_f = 30 m s^-1, and M_P2 = { 10,15 } M_Jup reproduce a better fitting of the far-infrared and mm-emission, however, the dip at 10 μm is less evident because of the homogeneous distribution of grains in the inner part of the disk (Fig. 4).

4.2. Polarised near-infrared emission and millimetre maps

In Fig. 7 we show synthetic HiCIAO H-band polarised intensity and ALMA Band 7 (850 μm) intensity images, for the models of an outer planet M_P2 = 5 M_Jup and each value of the fragmentation velocity v_f and turbulence parameter α_turb, assuming an inclination angle of i = 15° and a distance of 135 pc. The HiCIAO images are obtained convolving full resolution 1.6 μm images with a measured HiCIAO point spread function (PSF) in H-band taken from the ACORNS-ADI SEEDS Data Reduction Pipeline software (Brandt et al. 2013). The ALMA images are obtained with the CASA simulator (version 4.1.0) assuming Cycle 2 capabilities and one hour of total observation time (antenna configuration that gives a beam of ~0.15″ × 0.14″¹). These images trace opposite ends of the dust particle size distribution i.e. polarised intensity images at short wavelengths trace small (1–10 μm) dust grains, while sub-millimetre observations trace large (1–10 mm) grains.

Polarised intensity observations of the cases of M_P2 = 5 M_Jup, v_f = { 10,30 } m s^-1 and α_turb = 10^-2, in which there is no trapping of dust particles and an enhancement of the dust surface density of small particles close to the position of the planets, show ring structures at those locations. Note that the full resolution images are azimuthally symmetric and therefore the speckle features seen in the H-band images are due to the shape of the HiCIAO PSF and do not correspond to morphological features of the disk. The images show an asymmetry in brightness between the far (top) and near (low) regions of the disk which is a result of the forward scattering of starlight by the small dust particles and the fact that we are looking at an inclined disk (i.e. a face-on disk would not show such asymmetry in brightness). The ALMA images for these two cases do not show any gap because the disk is depleted of large particles at the outer parts of the disk.

When there is trapping, and there is still a population of small-dust particles distributed over the disk (case of M_P2 = 5 M_Jupv_f = 10 m s^-1 and α_turb = 10^-3), the intensity should be distributed over the whole disk with a smooth depletion between the two planets (Fig. 4). The PSF features mask this region, however, it is still possible to detect a ring of small particles that accumulate at the location of the outer particle trap. The ALMA image in this case shows a slight cavity between the two planets, but the two traps are too close for the gap to be resolved.

For α_turb = 10^-3 and v_f = 30 m s^-1, the HiCIAO image reflects a lack of small particles, when fragmentation and trapping are unlikely. The ALMA image shows a very faint ring of millimetre particles created at the outer edge of the gap, where grains are continuously growing to the allowed maximum size (Fig. 5). Note that the flux in both HiCIAO and ALMA images for all simulations with v_f = 30 m s^-1 has been multiplied by a factor of 1.5 and 20, respectively, in order to show the features, but these features are therefore difficult to detect.

In the cases of Fig. 7, to reproduce a disentanglement between large and small grains, the most feasible case is when v_f = 10 m s^-1, α_turb = 10^-3, and M_P2 = 5 M_Jup, where micron-sized particles are located over most of the disk (r> 5 AU), and mm-dust particles are trapped in narrow rings at ~7 AU and ~27 AU. In this case, ALMA in Cycle 2 is not capable of resolving the gap between the two particle traps, however, a future configuration of 50 antennas can resolve this gap.

5. Discussion

In this section, we focus the discussion on comparing the previous results to the case of the disk of SR 21. For this target, we consider that the observed gas-ring detected with CRIRES at 7 AU is the result of the inner planet opening a gap, which is located at 5 AU. The mm-cavity at ~36 AU arises from the interaction with the additional outer massive planet at ~15–20 AU.

5.1. The asymmetry

To fit the SR 21 ALMA observations, Pérez et al. (2014) fit a Gaussian profile in the radial and azimuthal direction following the vortex morphology introduced by Lyra & Lin (2013) and compare it to a Gaussian radial ring (see also Birnstiel et al. 2013). The best-fit model corresponds to a vortex, located at ~46 AU from the central star (~10 AU further than the fitting of a ring), centre at ~178° (measured from east to north), with a radial and an azimuthal width of 14.4 AU and 40.4 AU, respectively. The best fit corresponds to a vortex with aspect ratio (χ) of 2.8. Assuming particles of a = 1 mm size, with ρ_s = 1 g cm^-3, and Σ_g = 10 g cm^-2, their turbulent velocities within the vortex are about 16% of the sound speed. These turbulent velocities v_turb (${\mathit{v}}_{\mathrm{turb}}\mathrm{\propto }\sqrt{{\mathit{\alpha }}_{\mathrm{turb}}}{\mathit{c}}_{\mathrm{s}}$, see Ormel & Cuzzi 2007) correspond to a diffusion parameter α_turb of around ~2 × 10^-2 with temperatures of ~30–35 K at the vortex location. Although magneto-rotational instability simulations suggest values of α_turb to be in the range of 10^-3−10^-2 (e.g. Johansen & Klahr 2005), steady-state vortices are unlikely to form and be long-lived at the outer edge of the gap if these values for α_turb are assumed for the whole disk (Fig. 3). However, this high turbulence might be only at the vortex location and created by the vortex itself. On the other hand, their best fit for the vortex width (σ_r,V ~ 14.4 AU) seems to be too wide as compared to the scale height of the disk (~2 AU, considering a temperature of T ~ 35 K at the peak of the emission ~40 AU). To sustain long-lived vortices in protoplanetary disks, the gas velocities must remain sub-sonic and the radial width of a vortex therefore cannot be much larger than the scale height of the disk (see e.g. Barranco & Marcus 2005). Because of the ALMA beam size of ~27 AU, the vortex width is still not resolved and we cannot disprove the vortex-scenario.

Fig. 8 Synthetic H-band HiCIAO (left panel) and ALMA Band 7 (850 μm) (right panel) observations for the case of an outer planet of MP2 = 10 MJup, vf = 30 m s-1 and αturb = 10-2. White dashed circle indicates the inner of HiCIAO images masked in the Follette et al. (2013) observations of SR 21. Stripped circle in right panel indicates the size of the ALMA Cycle 2 beam.

An alternative scenario for this asymmetry is disk eccentricity. Eccentric gaps have been observed for some transition disks. For example, using NIR high contrast imaging, Thalmann et al. (2014) report an eccentric gap (e ≈ 0.3) for the LkCa 15 disk. Contrary to an anticyclonic vortex, over-densities in an eccentric disk do not trap particles. Because of the low gas velocities at the aphelion compared to the perihelion, gas, and dust densities are expected to be higher at the aphelion, creating asymmetries that could be also detected with high resolution and sensitivity observations of CO line profiles and continuum (Regály et al. 2011, 2014; Hsieh & Gu 2012; Ataiee et al. 2013). The contrast of the gas and dust densities are expected to be similar, which is very different from the case of a vortex where the contrast of the millimetre-sized dust density can be much higher than in the gas (Birnstiel et al. 2013). For example, with a gas surface contrast of ~4 due to the vortex (case of α_turb = 10^-3 and q_P2 = 6 × 10^-3) dust contrasts of more than 100 are expected.

From the hydrodynamical simulations in Sect. 3, both eccentricity and vortex formation are expected in different cases (Fig. 3). Eccentricity dominates any possible asymmetry when α_turb = 10^-2, with a maximum value of e ~ 0.1 for q_P2 = 6 × 10^-3 or M_P2 = 15 M_Jup and creates a slight gas surface density contrast in the azimuthal direction. However, spiral arms created by the outer planet may strengthen those asymmetries (Fig. 2). For α_turb = 10^-3, potential asymmetries result from the combination of vortex formation and disk eccentricity, besides the case of q_P2 = 2 × 10^-3 or M_P2 = 5 M_Jup for which the formed vortex dissipates at longer time-scales.

5.2. Spatial segregation of small and large grains

According to the dust distribution shown in Fig. 4 and SED fitting (Fig. 6), the best models to reproduce the SR 21 features are v_f = 30 m s^-1 and α_turb = 10^-2 for outer planet masses of M_P2 = { 10,15 } M_Jup. Figure 8 shows the HiCIAO and ALMA observations for one of these cases: the outer planet of mass M_P2 = 10 M_Jup where we see that the distribution of small particles is uniform, while the ALMA observations only show the ring caused by the outer dust trap, in reasonable good agreement with observations of SR 21.

In the synthetic images of Fig. 8, we do not reproduce the potential asymmetry due to eccentricity, since we only conduct dust evolution models in the radial direction by azimuthally averaging the results from hydrodynamical simulations. As explained in the previous section, the asymmetry in flux between the upper and lower regions of the disk in the HiCIAO images is caused by the strong forward scattering of the small dust particles due to the inclination of the disk (i = 15°). Although low inclination, it is enough to make this effect noticeable. This is an effect not detected by the Follette et al. (2013) observations which can be attributed to the fact that, although we have used realistic composition and size distribution for the dust particles, these are in our simulations smooth particles whereas real dust particles are expected to have rough surfaces. Particles with rough surfaces have a flatter phase function at high viewing angles (90 ± 15 in this case) and as a consequence the forward scattering of real particles is expected to be lower than that of the smooth particles used for these simulations (Min et al., in prep.). This would cause the asymmetry in flux to be reduced in real observations.

6. Summary

Transition disks reveal intriguing structures such as spiral arms, asymmetric cavities, and spatial segregation between small and large grains. In this paper, by combining hydrodynamical and dust evolution models, we explore the gas and dust density distributions when two planets are embedded in a disk. After computing the radiative transfer, we predict the observational signatures in disks with various disk viscosity (α_turb), fragmentation velocity threshold (v_f), and planet mass. In particular, we investigate the possibility of having a smooth distribution of micron size particles over the entire disk and a large cavity in millimetre grains, as recently observed in transition disks, such as e.g. SR 21 (Follette et al. 2013; Pérez et al. 2014), HD 169142 (Osorio et al. 2014), and MWC 758. (Andrews et al. 2011; Grady et al. 2013). We also look for possible asymmetries in the gas distribution.

Our framework is the following: we fix a 1 M_Jup planet in the inner disk at r_P1 around a 2.5 M_⊙ star and a more massive planet further out in the disk at 3.5 × r_P1, α_turb = [ 10^-3,10^-2 ], v_f = [ 10,30 ] m s^-1, and M_P2 = [ 5,10,15 ] M_Jup. Our conclusions are the following (see also Tables 2 and 3):

• We confirm that for specific values of α_turb and v_f, trapping at pressure maxima is possible as a result of planets embedded in the disk (Pinilla et al. 2012). Although trapping is very efficient for α_turb ~ 10^-3, for higher turbulence (α_turb = 10^-2), dust accumulation only occurs for massive planets (M_P2 = [ 10,15 ] M_Jup) and low fragmentation efficiency (i.e. high threshold velocities for destructive collisions). The sharpness of the ring-like structure seen at millimetre wavelengths strongly depends on α_turb (Fig. 4).

• When trapping is efficient, a smooth distribution of micron-sized particles over the entire disk can only be observed if the combination of planet mass and turbulence is such that small grains are not fully filtered out. For α_turb ~ 10^-3, an outer planet more massive than 5 M_Jup will lead to micron-grain filtration. On the other hand, with α_turb ~ 10^-2, the disk is very diffusive, and micron-sized particles are distributed throughout the disk.

• Depending on the disk viscosity and planet masses, asymmetries in the gas may exist at the outer edges of the gaps due to disk eccentricity or vortex formation. Eccentricities are induced by massive planets (Kley & Dirksen 2006), while long-lived vortices are created by massive planets in a low viscosity disk (Fig. 3). The simulations presented in this paper show that the cases in which particle trapping occurs and small particles are spread over the entire disk are also likely to show gas asymmetries caused by disk eccentricity. In this case, it is expected that the azimuthal dust distribution of all particles is identical to gas distribution.

• In the particular case of SR 21 transition disk, the observed cavity at millimetre wavelength (Andrews et al. 2011) can result from a planet located around 17–18 AU. In this case an inner planet of 1 M_Jup at 5 AU can create a ring of gas around 6–7 AU, which is in agreement with CRIRES observations of the gas (Pontoppidan et al. 2008), where dust particles can also accumulate and lead to an NIR excess. To obtain micron-sized particles in the entire disk, the outer planet should not be too massive (M_P2 ≈ 5 M_Jup), and α_turb = 10^-3, v_f = 10 m s^-1. If the disk viscosity is higher, the mass of the outer planet must increase to have effective trapping (M_P2 = [ 10,15 ] M_Jup), and therefore a ring-like emission at mm-wavelengths. In such cases, the fragmentation velocities should be velocities need to be significantly higher (v_f = 30 m s^-1 than for bare silicate grains (few m s^-1). These fragmentation velocities are expected for icy grains (Blum & Wurm 2008; Wada et al. 2009). For any of these cases, the resulting asymmetry is likely caused by eccentricity and it has low contrast of around ~1.5–2 in the gas, which may be enhanced by spiral arms.

Acknowledgments

The authors are thankful to H. Meheut for fruitful discussions and D. Fedele for providing the DIGIT data of SR 21. P. P. is supported by Koninklijke Nederlandse Akademie van Wetenschappen (KNAW) professor prize to Ewine van Dishoeck. T. B. acknowledges support from NASA Origins of Solar Systems grant NNX12AJ04G. Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA).

References

All Tables

All Figures

	Fig. 1 Sketch of the spatial segregation between small and large grains expected when two planets interact with the disk and lead to particle traps. Asymmetric structures may exist due to azimuthal trapping in vortices or/and disk eccentricity.
In the text

Fig. 2 Left panels: 2D gas surface density after 1000 orbits of the inner planet, when two planets with rP2/rP1 = 3.5 are interacting with the disk, with a planet-to-stellar mass ratio of qP1 = 4 × 10-4 for the inner planet and qP2 = 2 × 10-3 (left), qP2 = 4 × 10-3 (middle), and qP2 = 6 × 10-3 (right), and two different values of disk viscosity αturb = 10-2 (top) and αturb = 10-3 (bottom). Right panels: azimuthally averaged gas surface density from the corresponding hydrodynamical simulations and time-averaged over 200 orbits (from 1000 to 1200 orbits).

In the text

	Fig. 3 Azimuthally averaged vortensity ⟨ η2D ⟩ φ (blue lines) and disk eccentricity (red lines) versus the disk radial extension for the hydrodynamical simulations of Fig. 2 after 1000 (dotted lines) and 2000 (solid lines) orbits. The vertical lines correspond to the positions of the two planets.
In the text

	Fig. 4 Dust surface density for micron- (upper) and millimetre-sized particles (bottom), assuming two planets interacting with the disk at 5 AU and 17.5 AU and after 1 Myr of dust evolution. Two different thresholds are considered for the fragmentation velocity of particles vf = 10 m s-1 (left) and vf = 30 m s-1 (right). The vertical lines correspond to the positions of the two planets.
In the text

	Fig. 5 Dust density distribution at 1 Myr of evolution when two planets are interacting with a disk at 5 AU (1 MJup) and at 17.5 AU (5 MJup) with αturb = 10-2 (top), vf = 10 m s-1 (left) and αturb = 10-3 (bottom), vf = 30 m s-1 (right). The solid line corresponds to the particle size for which particles experience the highest radial drift, which is directly proportional to Σg (Eq. (1)).
In the text

Fig. 6 Spectral energy distribution resulting from the radiative transfer modelling, assuming grain size distribution after 1 Myr of evolution when two planets interact with the disk, and assuming different values for the disk viscosity αturb and fragmentation velocity of particles vf. The disk data are for SR 21 and are taken from Brown et al. (2007); Follette et al. (2013) and PACS (DIGIT PI: N. Evans).

In the text

Fig. 7 Synthetic H-band HiCIAO (left panels) and ALMA Band 7 (850 μm; right panels) observations for the models of an outer planet MP2 = 5 MJup, and each value of the fragmentation velocity vf and turbulence parameter αturb. The white dashed line in the left panels indicates the inner of HiCIAO images masked in the Follette et al. (2013) observations of SR 21. In the right panels, the size of the ALMA Cycle 2 beam is indicated by a stripped circle.

In the text

	Fig. 8 Synthetic H-band HiCIAO (left panel) and ALMA Band 7 (850 μm) (right panel) observations for the case of an outer planet of MP2 = 10 MJup, vf = 30 m s-1 and αturb = 10-2. White dashed circle indicates the inner of HiCIAO images masked in the Follette et al. (2013) observations of SR 21. Stripped circle in right panel indicates the size of the ALMA Cycle 2 beam.
In the text