Free Access
Volume 622, February 2019
Article Number A43
Number of page(s) 15
Section Interstellar and circumstellar matter
Published online 24 January 2019

© ESO 2019

1 Introduction

The process of star formation through the collapse of a pre-stellar core often leads to the creation of a protoplanetary disc made of gas and dust orbiting the proto-star (Williams & Cieza 2011). In the current paradigm of planet formation, protoplanetary discs are the birthplace of the first rocky bodies, also called planetesimals. These solids eventually become planets, either by core accretion or by gravitational instabilities, on timescales shorter than several Myr (Ribas et al. 2014). To date, more than 3500 confirmed exoplanets have been detected1, which clearly shows that planets are a common by-product of star formation. However, despite the large number of theoretical and experimental studies in the field, the process by which planets form from micron-sized dust grains remains elusive (Birnstiel et al. 2016, and references therein).

In fact, there are several barriers for planet formation such as the radial drift (Weidenschilling 1977; Laibe et al. 2012), the bouncing (Zsom et al. 2010), and the fragmentation (Blum & Wurm 2008) of solids during their evolution inside the disc. These physical processes severely compromise the survival and/or growth of such solids, and therefore prevent planet formation. Several mechanisms have been proposed to circumvent this problem: the presence of planetary traps (Paardekooper & Mellema 2004; Fouchet et al. 2007; Pinilla et al. 2012), streaming instabilities (Youdin & Goodman 2005; Jacquet et al. 2011; Johansen et al. 2012), anticyclonic vortices (Barge & Sommeria 1995; Meheut et al. 2012; Baruteau & Zhu 2016), dead zones (Kretke & Lin 2007; Dzyurkevich et al. 2010), meridional circulation (Fromang et al. 2011), photophoretic transport (Krauss & Wurm 2005; Cuello et al. 2016), radiation pressure effects (Vinković 2014), and self-induced dust traps (Gonzalez et al. 2017) among others. The common denominator of these processes is the concentration and/or growth of dust at a particular location of the disc, either by stopping the radial drift or by lowering the relative velocities of collision among solids. The former prevents solids from being accreted by the central star, while the latter allows particles to grow and eventually decouple from the gaseous phase.

In this work, we study the effect of shadow-triggered gaseous spirals on the dust distribution of a transition disc. In general, spiral wakes can be caused by different mechanisms: planet torques (Goldreich & Tremaine 1979), gravitational instabilities (Kratter & Lodato 2016, and references therein), stellar encounters (Pfalzner 2003), and illumination effects (Montesinos et al. 2016). Regardless of their origin, these perturbations cause a local symmetry breaking in the disc, which has profound consequences for structure evolution and planetesimal formation. In recent years, due to the spectacular increase in angular resolution and the development of extreme adaptive optics, it has been possible to reveal astonishing asymmetric features in several transition discs. These observations can be split into two kinds according to the wavelength: the near-infrared (NIR) or scattered-light observations, which trace the small micron-sized particles, and the (sub-)millimetre observations, which map the thermal emission of large (mm or cm) dust grains. In Table 1, we list some of the most interesting transition discs exhibiting visible spirals in the NIR and/or large asymmetric dust traps (e.g. HD 142527 and HD 135344B).

Elias 2-27 is a system of particular interest. Pérez et al. (2016) report two symmetric spirals in the dust thermal emission at 1.3 mm with ALMA. To our knowledge, this is the only detection of a grand-design spiral in the large grains distribution, i.e. close to the mid-plane. It remains unclear, however, whether the spirals and the shadows detected by Stolker et al. (2016) and by Benisty et al. (2017) are related to planet formation processes or not. For instance, in the protoplanetary disc around the F-type star HD 135344B, Stolker et al. (2016) report the detection of four shadows (A, B, C, and D) in the J band and at least three spirals arms. Surprisingly, shadow C is not detected in earlier observations (R, I, and Y bands), which seemsto indicate its transient nature. This could be explained by a local perturbation in the inner regions of the disc or by an accretionfunnel flow from the inner disc onto the star. In principle, an outer planet could excite a pair of symmetrical spirals, as suggestedby Dong et al. (2015), even though there is no evidence whatsoever of planetary gaps. Until now, the only works that have established a clear connection between shadows and spirals are Montesinos et al. (2016) and Montesinos & Cuello (2018). However, it is still unknown if this mechanism can explain the shadows in HD 135344B and how it affects the dust distribution in the disc.

After the detection of dusty spirals in Elias 2-27 by Pérez et al. (2016), many scenarios have been proposed and tested to explain the formation of such spirals. For instance, Tomida et al. (2017) report similar spirals in fairly massive three-dimensional discs (class I objects) formed via magnetohydrodynamic (MHD) effects during the phase of envelope collapse. In this case, the pair of grand-design spiral arms form due to gravitational instabilities in the disc. Interestingly, the spirals disappear after a few rotations and new spirals form due to growing instabilities in the gas. This spiral formation process occurs throughout the class 0 and class 1 phases. Meru et al. (2017) have also studied gravitationally unstable discs and show that it is possible to create spirals analogous to those of Elias 2-27 without considering MHD effects. In addition, the authors also consider the case of an internal (<300 au) and an external (>300 au) companion. The former configuration fails to produce spirals like the ones observed, while the latter is able to create grand-design spirals in the disc. However, the companion must be between 300 and 700 au away and have a mass below 13 Jupiter masses. Both scenarios, gravitational instability and the presence of an exterior companion, seem equally likely with the current observational constraints.

In this paper, we consider the formation of spirals triggered by shadows proposed by Montesinos et al. (2016), which form in the presence of light-obstructing material close to the star (Marino et al. 2015). This mechanism is due to the peculiar azimuthal acceleration caused by the shadows cast at the inner rim of the disc (see Fig. 3 in Montesinos et al. 2016). In this scenario, the spirals in the gaseous disc have lifetimes of the order of several thousand years and could in principle be detected in the NIR (H band, 1.6 μm). However, we are left with the following questions: can these shadow-triggered spirals efficiently trap dust? If yes, what are the observational signatures of this mechanism? It is not trivial to address these questions because the dynamics of dust strongly depends on particle size, intrinsic density, and gaseous local conditions (see Sect. 3.2).

The radial motion of dust particles embedded inside a gaseous sub-Keplerian disc, i.e. in the presence of a radial pressure gradient, has long been known and extensively studied (Weidenschilling 1977; Nakagawa et al. 1986; Laibe et al. 2012). More recently, the case of an azimuthal pressure gradient and the resulting dust dynamics was considered by Birnstiel et al. (2013). In both cases, solid particles tend to move towards the pressure maxima at velocities which depend on their Stokes number and the drag regime (see Sect. 2.2.1). So far, the dust behaviour in the presence of both radial and azimuthal pressure gradients has been modestly explored, and only for a handful of mechanisms: planetary wakes (Paardekooper & Mellema 2006; Ayliffe et al. 2012), disc instabilities (Rice et al. 2004; Dipierro et al. 2015; Tomida et al. 2017; Meru et al. 2017), and gravitational perturbations (Zsom et al. 2011; Dai et al. 2015). Thus, the dust dynamics in the presence of the shadow-triggered gaseous spirals proposed by Montesinos et al. (2016) remains to be studied.

The aim of this work is twofold. On the one hand, we analyse in detail the dust dynamics in the presence of shadow-triggered spirals, and on the other hand, we produce synthetic ALMA observations of the thermal emission of the dust. This allows us to address the two questions above about the dust trapping and its observability. The paper is structured as follows: in Sect. 2 we discuss the three main steps of our calculations; in Sect. 3 we report our results and the images obtained; in Sect. 4 we discuss our results; and in Sect. 5 we draw our conclusions.

Table 1

Detections of spiral and asymmetric dust traps in transition discs.

2 Numerical method

The approach followed in this work is based on three main steps, which connect our hydrodynamical simulations with the observations in a consistent way:

  • 1.

    gaseous hydrodynamical simulations with shadows;

  • 2.

    dust evolution on top of the gaseous distribution;

  • 3.

    synthetic ALMA observations of the dust emission.

The first hydrodynamical step, described in Sect. 2.1, allows us to obtain the dynamical fields of thegas (velocity, temperature, density, acceleration) required for the dust evolution calculation. The equations of motion for the dust particles and their numerical integration are given in Sect. 2.2. Finally, based on the obtained dust distributions, we present the simulations of the thermal emission of the dust and emulate ALMA observations of our system in Sect. 2.3.

2.1 Hydrodynamical simulations with shadows

We model the evolution of a self-gravitating circumstellar gaseous disc by means of the public two-dimensional hydrodynamic code FARGO-ADSG2 (Masset 2000; Baruteau & Masset 2008). The code FARGO-ADSG solves the Navier-Stokes, continuity and energy equations on a staggered mesh in polar coordinates (r, θ). The simulations were performed with an improved version of this code, which includes physical mechanisms such as cooling, stellar heating, and shadows. This allows us to simulate the evolution of a transition disc in the presence of projected shadows, as in Montesinos et al. (2016). These are considered regions of the disc for which the incoming stellar irradiation is equal to 0. Below, we describe the disc model considered, the energy equation with shadows, and the initialisation of the hydrodynamical simulations. The evolution of the disc is followed for about 10 000 yr.

2.1.1 Disc model

We consider a gaseous disc orbiting a solar-like star of mass M = 1 M and luminosity L = 1 L. The disc extends from rin = 100 au to rout = 600 au. The total disc mass is set to 0.25 M and its initial density profile is given by (1)

This model corresponds to a marginally stable disc for which the Toomre parameter remains above 1 throughout the simulation (Appendix B). The computational domain in physical units extends from rin to rout au over nr = 512 logarithmically spaced radial cells. The grid samples 2π in azimuth with nθ = 1024 equally spaced sectors. We also ran simulations at higher resolution with nr = 512 and nθ = 2048 and obtained the same results, which means that the above-mentioned resolution is accurate enough for our model.

The turbulent viscosity of the gaseous disc is modelled by considering an α-disc (Shakura & Sunyaev 1973), with the α parameter set to 10−4. We assume a constant opacity equal to κ = 1 cm2g−1. It is worth noting that the opacity directly impacts on the cooling rate Q of Eq. (7) in Sect. 2.1.2. Interestingly, as shown in Montesinos et al. (2016), shadow-triggered spirals appear for high and low values of κ: 130 and 1 cm2 g−1, respectively.A more realistic opacity depending on both the temperature and the density (e.g. the Rosseland mean opacity; Semenov et al. 2003) would have a minor effect on the disc for the range of temperatures considered in this work. This justifies our choice of a constant opacity. The only heat source we consider is the stellar irradiation and we do not take into account the viscous heating, which is dominant only for the very inner regions of the disc (Appendix A). All our models include self-gravity, meaning that the gas feels the gravitational potential of the star and of the disc itself. Given the mass of our system, this term cannot be neglected. Moreover, it has important effects for the dust evolution, as discussed in Sect. 2.2.1.

2.1.2 Energy equation and shadow model

We implement the shadows in the disc in the same way as in Montesinos et al. (2016). If there are no shadows, the disc is irradiated by stellar light at all azimuths. The stellar irradiation heating per unit area is given by Frank et al. (2002) (2)

where β is the albedo, L the stellar luminosity, and H the disc scale height. We assume β = 0, which means that all the incoming radiation gets absorbed by the disc. Assuming local hydrostatic equilibrium, the scale height H is given by (3)

where cs is the sound speed and vk the Keplerian velocity. During the calculation, the value of H is continuously updated according to the evolution of the temperature of the disc. Also, the viscosity is given by ν = α cs H.

We further assume the presence of a non-precessing edge-on disc (i = 90°) inside the large cavity. This type of structure can be created by an inner planetary or stellar companion, as in Facchini et al. (2018), and blocks the stellar light. In this configuration, two static and diametrically opposed shadows are cast onto the outer disc, as shown in Fig. 1. In the case of prograde rotating shadows, with rotation periods smaller than 20 kyr, planetary-like spirals can appear at the corotation region with the shadows (Montesinos & Cuello 2018). This effect is neglected in this work. In other words, our assumption of static shadows is valid for shadows with precession periods larger than roughly 20 kyr. It is worth highlighting that this inner disc is not treated in the hydrodynamical simulations, only its radiative effect are considered. Thus, we modify the stellar irradiation rate to account for the shadows, (4)

where the function s(θ) is given by (5)

with A = 0.999 and σ = 0.05. This creates two shadows at θ = 0° and at θ = 90° with a width of about 29°, where the stellar heating rate is set to nearly zero. The width is motivated by the case of HD 142527, as reported by Marino et al. (2015), and is also comparable with HD 135344B’s shadow B, as reported by Stolker et al. (2016). By construction, the illuminated and shadowed regions are smoothly connected.

Our hydro code solves the non-stationary equation for the thermal energy density, e, given by (6)

where v is the gas velocity, P the gas pressure, the stellar heating rate per unit area (cf. Eq. (4)), and Q the radiative cooling of the disc. For our calculations, we use (7)

where σSB is the Stefan-Boltzmann constant and τ the optical depth given by . In order to close the system of equations, we consider an ideal equation of state, , with T the gas temperature and , where kB is the Boltzmann constant, μ the mean molecular weight of the gas, and mp the mass of the proton. Thus, the relationship between the thermal energy density and the temperature reads (8)

where γ is the adiabatic index. Finally, in order to avoid unrealistically low temperatures in the disc, we set the Cosmic Microwave Background (CMB) temperature (2.7 K) as the floor temperature in our simulations.

thumbnail Fig. 1

Sketch of the inner regions of the disc model considered in this work. The outer disc is face-on, while the non-precessing inner disc is edge-on. The light-obstructing material around the star casts shadows at two diametrically opposed regions of the outer disc, in dark blue.

Open with DEXTER

2.1.3 Initialisation

In a stationary regime, the energy received from the star equates the energy loss through the disc’s surface. In order to reach a steady-state regime as quickly as possible, the disc is initialised under the condition . Assuming vertical hydrostatic equilibrium, the gas temperature is related to the aspect-ratio of the disc as follows: (9)

Definingh = Hr = h0rf and equating Eq. (2) with Eq. (7), we obtain (10)

where f is the flaring index of the disc and is given by f = dln(Hr)∕dlnr. Equation (10) gives a recipe to initialise the disc’s aspect ratio as a function of the initial disc and star properties. Formally, this equation is not self-consistent because the aspect-ratio depends on its own derivative. However, the dependence on the flaring index is very weak (∝ f1∕7). Therefore, as a first-order approximation, we can assume that f is nearly constant and equal to f ≃ 1∕7.

2.2 Dust evolution calculation

The outputs of our hydro simulation are used as inputs for our dust code, which computes the evolution of an ensemble of N independent dust particles. Given that it is a post-processing calculation, the gaseous phase is not affected by the dust distribution. On the contrary, the dust is affected by the gas dynamics due to the aerodynamical drag, which strongly depends on the local disc conditions (P, T). It is useful to introduce the stopping time (ts) defined as the time it takes for the (Keplerian) dust to reach the same velocity as the (sub-Keplerian) gas. The dimensionless parameter called Stokes (denoted St) is also relevant for dust dynamics since it measures the dust coupling to the gas, i.e. the drag regime. In the Epstein regime, valid for particles with sizes smaller than 9∕4 of their mean free path, St reads as (11)

where ΩK is the Keplerian frequency, ρd the bulk density of the particle, a its radius, ρg the volumetric gas density, and cs the speed of the sound. There are three different regimes: for St ≪ 1 the dust is well-coupled to the gas; for St ~ 1 the dust particles are marginally coupled to gas and feel the strongest radial drift; and for St ≫1 the particles arecompletely decoupled from the gas. Given the density profile considered (Eq. (1)), for a fixed grain size a, the Stokes number is an increasing function of the distance to the star.

2.2.1 Equations of motion and numerical integrator

Our dust code is based on previous work by Paardekooper (2007) and Zsom et al. (2011). In addition to the aerodynamical drag, dust particles feel the gravitational potential from the star and from the gaseous disc. For a given dust particle of index i, we denote (ri, θi) and (vr,i, vθ,i) the polar coordinates and velocities, respectively. Therefore, the equations of motion read as

where Li is the specific angular momentum, and Fr,i and Fθ,i are respectively the radial and azimuthal drag forces for particle i. The gravitational potential (Φ) is defined as Φ = GM*R + Φd, where Φd is the disc’s gravitational potential, which is given by the Fourier transform of the surface density Σ of the disc. In order to avoid this computationally expensive calculation, we use the radial and azimuthal accelerations given by the hydro code. It is worth noting that this term increases the orbital velocity of the gas compared to non-self-gravitating discs. Thus, it is important to include it because otherwise the particles in the outer parts of the disc would orbit too slowly with respectto the gas. This would result in a constant and unrealistic tailwind, which would send dust particles outwards.

In general, Eqs. (12)–(15) also have inertial forces that arise because the coordinate frame of the hydro code is always centred on the star and not on the centre of mass. In the case of a binary system for instance, this introduces inertial forces that have to be taken into account. Given that our system is axisymmetric, the centre of mass is always at the centre of the star and we can therefore neglect these additional forces. The drag force on the dust particles is computed as (16)

where Δv is defined as the vectorial relative velocity between the dust, vdust, and the gas, vgas. The radial and azimuthal components of Fdrag are equal to Fr,i and Fθ,i in Eqs. (14) and (15), respectively. The expression of the Stokes number used in Eq. (16) is more general than that in Eq. (11), which is only valid in the Epstein regime. In Appendix C we report the method used to compute the Stokes number in Eq. (16) and the details of the numerical integrator. We use a standard first-order leapfrog integrator to solve Eqs. (12)–(15), similar to that used by Paardekooper (2007). In Appendix D we also detail the benchmark analysis of the dust evolution code, which ensures the accuracy of the numerical integrator.

It is worth highlighting that, in order to ensure an accurate integration of Eqs. (12)–(15), we perform interpolations between every two consecutive gas snapshots3. More specifically, we feed the dust code with the values of the required physical fields of the gas (i.e. Σ, T, P, vgas, and agas). By doing so, we consistently take into account the evolution of the gas disc to compute the resulting motion of the dust particles.

2.2.2 Dust simulation set-up

For each dust simulation, we use 100 000 dust particles with a fixed size between amin= 1 μm and amax = 10 cm. The particles areinitially distributed according to the initial gas surface density of Eq. (1). The bulk density of the solids, ρd, is set to 2.0 g cm−3 to model a mixture of pyroxene, refractory organics, and ices.

This distribution is a simplified initial condition since we are not considering a radial-dependent size distribution. This choice allows us to better follow the dust kinematics by directly comparing the different grain species considered. For instance, in a more evolved protoplanetary disc, one would expect to have a more compact disc of millimetric grains compared to the disc of micron-sized grains (Laibe et al. 2012). The implications of this assumption are discussed in Sect. 4.2.

2.3 Radiative transfer and synthetic ALMA observations

2.3.1 Radiative transfer with RADMC-3D

Since our goal is to connect our simulations with ALMA observations, we need to compute the thermal emission of the dust distribution at millimetre wavelengths. To perform radiative transfer calculations we use the RADMC-3D sofware4 v0.41 (Dullemond et al. 2012). However, RADMC-3D input is a continuous three-dimensional dust density distribution for each particle size in the radiative transfer model. Thus, we need a conversion from a discrete set of dust distributions of different sizes (Sect. 2.2) to a continuous dust distribution. This is done in three main steps:

1) We spread all dust particles over a two-dimensional surface density grid with an exponential function: (17)

Here, Σdust,i is the surface density that one particle of mass mi at location ri contributes to the overall dust surface density. The factor 1∕(πD2) conserves the mass of the particle, but the overall surface density has to be gauged later nevertheless. As kernel size, we use D = 100 au. This seems large compared to the size of our system, but as seen in Sect. 3.1, this is roughly the width of one spiral arm. Smaller values lead to noisy images.

2) We sort the particles into 12 particle size bins. Since we have particles between amin = 1 μm and amax = 10 cm, this leads to two bins per decade. We then sum the individual single-particle surface densities in 12 total dust surface densities according to their size, (18)

where aj,1 and aj,2 are respectively the lower and the upper limits of the size bin of index j.

3) Finally, these surface densities are vertically distributed with a Gaussian to get spatial mass densities according to (19)

with the dust scale height given by Birnstiel et al. (2010) as (20)

The combined mass of all our 100 000 dust particles is far too low for a typical disc. Thus, to obtain a dust-to-gas ratio of ~ 1%, we multiply our densities with a constant factor such that the total mass of the dusty disc is Mdust = 0.25 × 10−2 M.

We feed RADMC-3D with the dust temperature and density distribution obtained from the hydro and the dust evolution simulations, respectively. The radiative transfer is then performed through the Monte Carlo method of RADMC-3D based on Mie theory (Bohren & Huffman 1983). We use the optical constants taken from Dorschner et al. (1995) assuming the particles are spherical and consist of pure pyroxene (Mg_0.7Fe_0.3 SiO_3). It is worth highlighting that we modified RADMC-3D, such that no photons are sent from the star within two sectors of angle 15° at azimuth 0° and 180°, as shown in Fig. 1. This mimics the shadows that are cast by the inner disc as in Montesinos et al. (2016). The star has an effective temperature of Teff = 5700 K. The radiative transfer model is set up in spherical coordinates (R, θ, ϕ) with NR = 128, Nθ = 90, and Nϕ = 64 grid points in the radial, azimuthal, and polar directions, respectively. We use a number of photon packages equal to 1 million for our calculations. We run RADMC-3D simulations at a wavelengths of 1.6 μm and 1.3 mm in order to model the emission of micron- and millimetre-sized particles, respectively. In Fig. 7, we report the images obtained for our disc model after 10 000 yr of evolution.

2.3.2 Synthetic ALMA images with CASA

For direct comparison with recent ALMA observations of transition discs, we produce synthetic ALMA observations with the CASA package (McMullin et al. 2007). To do so, we scale the radiative transfer calculations of Sect. 2.3.1 to a distance of 140 pc from the Earth. The bandwidth is 8 GHz and the integration time is 5 h. We use the “alma.cycle15.cfg” array configuration. With this set-up we obtain a synthetic beam of 0.08 arcsec2. In Sect. 3.3 we present the synthetic ALMA images obtained. Finally, in Sect. 4.3 we discuss the detectability of the dusty spirals in our disc.

3 Results

In this work we consider two simulations: Sha-r100 and NoSha-r100, which correspond to the disc model described in Sect. 2.1.1 with and without shadows. In both cases, self-gravity (SG) is included. This allows us to disentangle shadowing effects from SG.

3.1 Spiral formation and evolution

In Fig. 2 we show snapshots at four different times of the Sha-r100 simulation. Initially, the surface density of the disc is azimuthally symmetric by construction. In the temperature maps of the disc, the shadowed regions of the disc can be easily seen at azimuths 0° and 180°. The disc is colder in the shadowed regions because the stellar light is blocked by the inclined inner disc. The cold regions of the disc do not exactly overlap with the regions under the shadows. This effect is due to the anti-clockwise rotation of the disc and the time it takes for the material to reach thermal equilibrium. For instance, in the inner regions the gas crosses the shadows at high angular velocities, and thus remains hotter compared to the gas of the outer regions at the same azimuth. This also explains why the shadows’ interfaces are smoother than those defined by Eq. (5). This is also seen in the pressure map, where the shadowed regions have lower values compared to the illuminated regions. After 4000 yr, faint grand-design spirals appear in the pressure and temperaturedistributions, extending from 200 au to roughly 400 au. These asymmetries increase in contrast with time, as canbe seen for the additional snapshots (7500 and 10 000 yrs) where the spirals extend from 150 to 500 au.

In clear contrast, the NoSha-r100 simulation does not show grand-design spirals in the gas distribution for comparable timescales. In Fig. 3, we show the (r, θ) pressure maps for NoSha-r100 and Sha-r100 (left and middle panel, respectively). To easily identify the differences, we plot in the rightmost panel the normalised perturbation amplitude (PP0)∕P0. In the darker regions the pressure is 20 to 40% higher compared to the unperturbed case. As is shown in Sect. 3.2, these local pressure maxima are able to trap dust. On the contrary, brighter regions are characterised by lower pressures. We see for instance that the shadows create two bright lanes at θ = 0° and θ = 180°. Hence, the shadows trigger spirals in the disc, which then constitute regions of high pressure compared to the rest of the disc. Although we take into account self-gravity effects, the disc model considered in this study is gravitationally stable (Appendix B). It is worth noting, however, that for self-gravitating discs, spirals appear faster in the presence of shadows compared to the case without shadows (Montesinos et al. 2016). In other words, the shadows render the disc more prone to fragmentation. To make sure that the shadows are indeed causing the spirals in the disc, we performed a control run where we artificially turned off the self-gravity. In this case, the shadow-triggered spirals still form in the disc proving that self-gravity is a sub-dominant effect.

In our models the spiral pattern locally co-rotates with the disc, as in the case of self-gravity induced spirals. This contrasts with planetary-induced spirals, which follow the Keplerian frequency of the planet in solid rotation. Due to the quasi-Keplerian motion of our spirals, the dust and the pressure overdensity travel at almost the same angular speed, making it easier for the dust to concentrate in the local pressure maxima (Sect. 3.2).

We note as a side remark that the number of spirals arms in the outer disc is correlated to the number of shadows cast. By artificially setting one shadow or more than two in the simulations, we observe an equal number of spiral arms and shadows after a few thousand years. However, the structure and the evolution of each individual spiral arm remain the same. These spiral-shaped overdensities in the gaseous distribution constitute a “sweet spot” for dust trapping. The efficiency of this mechanism is detailed below.

thumbnail Fig. 2

Snapshots of the mid-plane gas temperature (top panels) and the vertically integrated pressure (bottom panels) at four different snapshots of the Sha-r100 simulation –from left to right panels: 2500, 5000, 7500, and 10 000 yr.

Open with DEXTER
thumbnail Fig. 3

(r, θ) maps of the vertically integrated pressure in NoSha-r100 (P0, left panel), in Sha-r100 (P, middle panel), and the normalised perturbation amplitude (PP0)∕P0 (right panel). All the snapshots are taken after 10 000 yr of evolution.

Open with DEXTER

3.2 Dust particles dynamics

Figure 4 shows the distribution of dust, after 7500 and 10 000 yr, obtained with the dust evolution code described in Sect. 2.2. We also show the pressure map to compare the location of the dust spirals with the local pressure maxima. We extract three particle sizes from the dust simulation as three different dust species: a ≈ 100 μm, a ≈ 1 mm, a ≈ 1 cm. Given the different sizes considered, each species has its own dynamical behaviour (Eq. (11)): the larger the particles, the larger their Stokes number. In the outer parts of the disc, the Stokes numbers are initially close to or below 1. In this case, the larger the particles, the faster the drift towards the pressure maxima. These particles are located inside the spirals and in the inner region of the disc. This can be seen in Fig. 4, where in the outer parts of the disc the centimetric particles drift efficiently towards the spirals. The centimetric particles below r ≈ 300 au, being far from the spirals, drift inwards and concentrate at the disc’s inner edge. Since the pressure (or equivalently the density) inside the spirals is higher than the local mean pressure, the Stokes number is expected to decrease. This effect can be seen in Fig. 4 for centimetric and millimetric particles. In fact, particles with Stokes numbers close to 1 are efficiently trapped inside the pressure maxima. Consequently, the regions between the spiral arms get depleted in dust with time. Hence, on longer timescales the dust distribution is expected to efficiently concentrate inside the shadow-induced spirals.

In order to study the characteristic behaviour of each species, we select individual but representative particles at a specific radial and azimuthal position of the disc at the end of the simulation. Then, we trace back their evolution in the disc to find their initial positions. We chose two sets of particles: one at r = 300 ± 5 au and another at r= 500 ± 5 au, both at θ = 0 ± 5°. It is worth noting that taking θ = 180 ± 5 we find similar results. Figures 5 and 6 show the distance to the star, the azimuth, the radial velocity, and the Stokes numbers as a function of time for particles of different sizes belonging to each set. We show three representative particles for each size bin: 100 μm in red, 1 mm in green, and 1 cm in blue. By selection, the particles plotted have the same final radial and azimuthal positions. However, these were different at t = 0 yr. In fact, in Fig. 5, we can see that the particles with small Stokes numbers (blue) remain at roughly the same radial position during the whole simulation, while those with St ≈ 1 (red) experience an outward movement. Particles with St ≈ 0.1 show the same behaviour with a less rapid outward drift. Similar but more pronounced effects are observed for particles at r = 500 au at the end of the simulation. This is due to the greater pressure gradients present in the disc at larger stellocentric distances (see Fig. 3). In addition, we see that some particles with St ≈ 1 beyond 500 au at t = 0 initially drift inwards and then concentrate at 500 au. This phenomenon is consistent with dust drifting towards5 the pressuremaximum located at the spiral position. This is in contrast with what is observed for an unperturbed disc (i.e. without shadows) where all the dust drifts inwards (see Appendix D). Moreover, the radial velocity is in agreement with the dynamical behaviour described. In fact, we observe that this velocity (typically negative) becomes slightly positive because of the presence of the pressure maxima. Also, it shows small oscillations due to the shadows’ crossing. Hence, the gaseous spirals act as efficient dust traps gathering material of varied sizes, starting at different radial positions.

thumbnail Fig. 4

Gas and dust distributions in the (r,θ)-plane after 7500 (left panels) and 10 000 (right) yr of evolution. Top row: vertically integrated pressure map of the disc. The rows below show the dust distribution coloured with the logarithm of the Stokes number – from top to bottom panels: 100 μm, 1 mm, 1 cm.Spirals in the gas and in the dust have the same morphology.

Open with DEXTER
thumbnail Fig. 5

From top to bottom panels: radial and azimuthal positions, radial velocity, and Stokes numbers as a function of time for three groups of particles, which end at r = 300± 5 au and θ = 0 ± 5. We plot 100 μm particles in blue, 1 mm particles in green, and 1 cm particles in red. Each group contains three particles of equal size.

Open with DEXTER
thumbnail Fig. 6

Same as Fig. 5, but for r = 500 ± 5 au.

Open with DEXTER

3.3 Synthetic observations

Figure 7 shows the radiative transfer calculation for Sha-r100 after 10 000 yr of evolution at λ1 = 1.6 μm and λ2 = 1.3 mm obtained following the steps in Sect. 2.3.1. As mentioned in Sect. 1, the former traces the small dust and the latter the millimetric grains. The value of λ2 was chosen in order to have the same wavelength as the observations reported by Pérez et al. (2016). We assume the system is 140 pc away from the Earth and has an inclination equal to 13°, i.e. almost face-on. Both imagesclearly show the two diametrically opposed shadows, which are expected by construction. In addition, faint spiral patterns are seen at λ1 = 1.6 μm, while these appear sharper and more radially concentrated at λ2 = 1.3 mm. The combined effect of radial drift and dust trapping inside the spirals explains the difference between thetwo images. We note, however, that the thermal emission from the disc inner rim at roughly 100 au is at least one order of magnitude brighter than the one emitted by the spiral arms. This has dramatic consequences for the detectability of such thermal emission.

In Fig. 8 we show the post-processed image with CASA, which emulates an ALMA observation (see Sect. 2.3.2). In particular, the presence of two dips at θ = 0° and θ = 180° shows that it is in principle possible to detect the shadows at millimetre wavelengths. This is in good agreement with the ALMA synthetic images obtained by Facchini et al. (2018), where an inclined inner disc casts shadows onto the outer disc. Unfortunately, in our case, due to the low signal-to-noise ratio (S/N) in the spirals it is not possible to detect any significant thermal emission from the dusty spiral arms. Even considering longer integration times and/or other antenna configurations, the dusty spirals remain below the detection threshold.

thumbnail Fig. 7

Dust thermal emission at 1.6 μm (left panel) and 1.3 mm (right panel) after 10 000 yr obtained with RADMC-3D. The source is placed at d = 140 pc and the disc is seen almost face-on (i = 13°).

Open with DEXTER
thumbnail Fig. 8

Synthetic observations with ALMA of Sha-r100 after 10 000 yr. The beam is indicated in the bottom left corner.

Open with DEXTER

4 Discussion

4.1 Grand-design dusty spirals

We considered a new mechanism able to produce grand-design dusty spirals in transition discs. Remarkably, there is no need to include a massive planetary or stellar companion, as in Dong et al. (2015). We only require two diametrically opposed shadows cast at the inner rim of a transition disc, as in Montesinos et al. (2016) and Facchini et al. (2018). The resulting gaseous grand-design spirals are able to trap dust of different sizes in the disc. Our calculations show that the dust particles with St ~ 1 are the onesthat gather most efficiently in the spirals. Interestingly, these particles move towards the pressure maxima from inwards and from outwards. This effect is not seen in the absence of shadows and the resulting spirals. Even though the morphology of our dusty spirals (see Fig. 7) is similar to those observed by Pérez et al. (2016), our synthetic ALMA image shows that these shadow-triggered spirals cannot be detected at 1.3 mm. Hence, shadows are a promising mechanism to produce dusty spirals in disc, but it is not possible to reproduce the observations of Elias 2-27. Moreover, Pérez et al. (2016) do not report any shadow detection. However, this mechanism should be considered for other transition discs exhibiting large cavities and shadows (e.g. HD 142527, HD 135344B, and HD 100453).

It is worth noting that we studied the evolution of shadow-triggered spirals for only 10 000 yr, a timescale too short to consider grain growth and fragmentation inside the spirals. Although our dust evolution code has a grain growth module, the typical grain growth timescales in the region of interest exceed the simulation time by at least one order of magnitude. That is the reason why we neglected this phenomenon, even though we expect particles to collide and grow inside the spirals. Unfortunately, simulating the disc for a prolonged time is problematic due to the large number of snapshots we need to store in the memory.

4.2 Dust trapping and mixing

In Sect. 3.2, we analysed the dust trapping by the shadow-triggered spirals in the gas. This phenomenon depends on the grain size and influences the dust dynamics in the disc dramatically. Figure 4 clearly shows that dust particles concentrate in the pressure maxima of the gas, which results in a dusty disc with spirals. The fact that particles are trapped in these overdensities prevents these solid bodies from being accreted by the central star. In the absence of spirals, only particles with St ~ 1 experience an extremely fast radial drift (Weidenschilling 1977). Eventually these particles are lost onto the star, which renders planetesimal formation rather problematic in the disc. Because dust particles can gather and grow efficiently inside the spirals, they constitute remarkable sweet spots for grain growth in the disc. This gathering effect could lead to streaming instabilities at this particular location (Youdin & Goodman 2005; Johansen et al. 2007; Jacquet et al. 2011), but this effect has not been taken into account in our simulations. Provided that strong dust clumping happens, the radial drift is slowed down in these overdense regions. Once the self-gravity of the solids becomes important, this mechanism can lead to the formation of planetesimals (Johansen et al. 2007, 2012).

In addition to this effect, we expect an efficient mixing of dust species inside the spirals. In Figs. 5 and 6 we observe a mixing of particles inside the spiral at ~ 300 au and at ~500, which initially started at different radial and azimuthal distances. In fact, particles with a different initial Stokes number drift at different radial velocities until they are trapped in the overdensity. This happens both in the inward and outward direction. Once this happens, there is a broad dust size distribution in the spiral. The fact that the trapped grains come from regions of the disc at different temperatures may result in different solid compositions (ices, volatile organics, refractory organics, silicates, iron, etc.). The mixing among different species is controlled by the numerator of Eq. (11), ρ and a, as discussed in Pignatale et al. (2017). In fact, a dense but small particle (ρd = 10 g cm−3, a = 0.1 mm) behaves in the same way as a lighter but larger particle (ρd = 1 g cm−3, a = 1 mm). This has dramatic effects on the resulting composition of solids in the disc, and also affects their further growth (Pignatale et al. 2017, and references therein).

However, caution is required when interpreting our results given that we are considering an initial size distribution that does not depend on the radius. If we consider a disc with a realistic size segregation, then the emission from the outer regions at λ = 1.3 mm would be fainter and that from the inner regions stronger. In fact, due to radial drift, we expect to have a more compact radial distribution of millimetric grains compared to micron-sized ones.

4.3 Detectability of shadow-triggered dusty spirals

In Fig. 8 we show a synthetic ALMA observation corresponding to Sha-r100. Although the shadows are clearly detected as dips in the emission at 1.3 mm, the thermal emission from the spirals is comparable to the noise level. Hence, even for longer integration times, their emission will remain below the detection limit. Unrealistically higher disc masses or dust-to-gas ratios could increase the thermal emission, thus rendering such detections feasible. If the canonical value of 1% of the latter is in fact too conservative, then shadow-triggered spirals could potentially be observable. The size of the cavity also plays an important role. As a matter of fact, discs with smaller cavities have higher fluxes because of the presence of hotter material closer to the star. This decreases the flux ratio between the spiral arms and the inner disc, which makes the detection even more challenging. Consequently, transition discs exhibiting shadows are ideal targets to look for spiral patterns in the thermal emission. It is worth noting nevertheless that considering a lower disc mass (e.g. 1% M), would produce even fainter spirals. Therefore, the disc model presented in this work should be seen as an optimistic scenario for this kind of spiral formation and detection.

Interestingly, some recent techniques allow structure to be extracted from noisy data well below the 5σ threshold. For instance, in the context of galaxy observation, Akhlaghi & Ichikawa (2015) describe a method able to extract spiral arms at very low S/N. Contrary to other existing techniques, signal-related parameters (such as theimage point spread function or known object models) are irrelevant for this method. Thus, assuming that dusty spirals of this kind form in nature, some previous (sub-)millimetre detections may have missed them because of standard data reduction procedures. To our knowledge, this kind of observational signatures has not been observed with ALMA yet. However, based on our results, we are able to predict that, in the presence of shadows, the thermal emission of the dust should exhibit faint grand-design spirals and significant dips (see Fig. 7).

5 Conclusions

In this work, we have developed a novel numerical tool, compatible with FARGO-ADSG, to explore dust dynamics for hydrodynamical simulations in a consistent manner. Although we do not include the back-reaction of dust on gas, this does not affect our results given the dust-to-gas ratios considered here. The methodology detailed in Sect. 2 allows the evolution of the dust distribution and its resulting thermal emission to be computed. In this regard, this tool successfully bridges the gap between theory and observations. It is particularly well-suited to model some of the transitional discs exhibiting puzzling dust distributions (see Table 1).

With the detection of dusty spirals in Elias 2-27 (Pérez et al. 2016) and the mechanism described in Montesinos et al. (2016) in mind, we explored the trapping of dust in discs with gaseous spirals due to illumination effects. In our case, the shadows are provoked by an inclined inner disc which obstructs the stellar light at large radii. The main results of this work can be summarised as follows:

  • The gaseous spiral waves induced by shadows are able to accumulate and trap dust particles. This mechanism deeply affects the fate and the evolution of the solid material in the disc. We showed that for millimetre- to centimetre-sized particles it is possible to circumvent the radial-drift barrier.

  • In our simulations, the most efficient trapping occurs for solids with sizes ranging from 1 mm to 1 cm. However, this size range strongly depends on the disc model considered (Laibe et al. 2012).

  • The dust trapping that occurs inside the gaseous spirals results in an efficient mixing of dust particles of different sizes and densities. This has deep implications for planetesimal and comet formation at large stellocentric distances, since spirals gather material from different regions of the disc.

  • The obtained infrared (1.6 μm) and thermal (1.3 mm) emissions exhibit both shadows and spiral patterns in the disc. However, the brightness of the dusty spirals is too faint to explain the recent ALMA observations of Elias 2-27.

The methodology presented here can be applied to explore dust dynamics in transition discs with vortices, multiple planets, and moving shadows. By doing so, it is possible to produce robust observational signatures of different physical processes in action. This will hopefully help observers to interpret their current and future dust detections.


We thank the two anonymous referees who have reviewed this manuscript for the very constructive comments which allowed us to improve the quality of this work. N.C. acknowledges financial support provided by FONDECYT grant 3170680. N.C., M.M., and J.C. acknowledge financial support from Millenium Nucleus grant RC130007 (Chilean Ministry of Economy). M.M. acknowledges support from the Millennium Science Initiative (Chilean Ministry of Economy) and the CHINA-CONICYT fund, 4th call. S.M.S. gratefully acknowledges support through the PUC-HD Graduate Student Exchange Fellowship, which is part of the academic exchange programme between the Institute of Astrophysics of the Pontificia Universidad Católica (IA-PUC) and the Center for Astrophysics at the University of Heidelberg (ZAH), financed by the German Academic Exchange Service (DAAD). This work was partly carried out while J.C. was on sabbatical leave at MPE. J.C. and N.C. acknowledge the kind hospitality of MPE, and funding from the Max Planck Society through a “Partner Group” grant. J.C. acknowledges support from CONICYT-Chile through FONDECYT (1141175) and Basal (PFB0609) grants, and from the ICM (Iniciativa Científica Milenio) viathe Núcleo Milenio de Formación Planetaria grant. The Geryon/Geryon2 cluster housed at the Centro de Astro-Ingenieria UC was used for the calculations performed in this paper. The BASAL PFB-06 CATA, Anillo ACT-86, FONDEQUIP AIC-57, and QUIMAL 130008 provided funding for several improvements to the Geryon/Geryon2 cluster.

Appendix A Irradiation versus accretion heating

Here we aimto show that we can safely neglect the effect of the viscous heating in the disc. On the one hand, we consider the viscosity of the gaseous disc, which leads to the production of heat. The heating per gram of gas is proportional to the the viscosity coefficient ν and the square of the shear. In an accretion disc, the total heating per unit area, , reads (A.1)

Assuming a steady-state disc, the accretion rate can be written as = 3πΣν. Thus, we obtain (A.2)

On the other hand, the stellar heating directly depends on the amount of radiation received by the disc. The flux of stellar radiation at distance r from the star is given by . The irradiating flux is the projection of this flux onto the surface of the disc and includes a factor sin (ϕ), which can be approximated by ϕ for small grazing angles. Typically, for a flared disc ϕ ≈ 0.05. The disc receives radiation from the upper and the lower part; therefore, there is a factor of 2 as well. Finally, we obtain that the irradiation heating per unit area, is given by (A.3)

Hence, the criterion for which the irradiation heating dominates over the accretion heating reads as follows: (A.4)

Plugging in a reasonable value for the stellar accretion ( = 10−7M*), we obtain that r ≥ 11 au. This condition is fulfilled by our disc, so the viscous heating can be neglected.

Appendix B Disc temperature and Toomre parameter

In Figs. B.1 and B.2 we show the initial and final values of the disc temperatures and the Toomre parameter, respectively. Because the disc radiates energy during its evolution, the final temperatures are lower than the initial ones, roughly between 12 and 5 K. This is consistent with the expected temperatures between 100 and 600 au for an irradiated disc (see Eq. (9)). The Toomre parameter is defined as Q = csκfπGΣ, where κf is the epicyclic frequency. This dimensionless quantity compares the strength of thermal pressure and shear (numerator) to self-gravity (denominator). If Q > 1 (Q < 1), then the system is stable(unstable) to its own gravity. In our simulations, the value of the Toomre parameter throughout the entire evolution remains everywhere above 1. Regions with Q values slightly above 1, which are found in our model, can be considered as marginally stable (Kratter & Lodato 2016). However, even in that regime, the perturbation introduced by the shadows dominates over self-gravity effects, leading to spiral formation (see Montesinos et al. 2016, for a more thorough discussion).

thumbnail Fig. B.1

Initial and final radial temperature profiles for the Sha-r100 simulation. The range of temperature is consistent with Eq. (9).

Open with DEXTER
thumbnail Fig. B.2

Initial and final radial Toomre parameter profiles for the Sha-r100 simulation. Values above 1 correspond to gravitationally stable discs. This is the case for the disc model considered in this work.

Open with DEXTER

Appendix C Numerical integration of the equations of motion

As mentioned in Sect. 2.2, the expression of the Stokes number used in the code is different than in Eq. (11). In the dust code, we use the radial and the azimuthal gravitational accelerations given by FARGO-ADSG in order to compute the Stokes number. For a dust particle of index i, we have (C.1)

where Kn = λmfpa is the Knudsen number, λmfp the mean free path of the gas, and a the radius of the dust particle. The parameter ζdust = 3.3 g cm−3 is the bulk mass density of the dust particles, and ρgas the mid-plane mass density of the gas. The parameters gdrag and kdrag are the drag coefficients in the Epstein and in the Stokes drag regimes, respectively. They are given by (C.2)

and (C.3)

Re is the Reynolds number, which is defined as (C.4)

where m = |Δv|∕cs is the relative Mach number of the dust particles to the gas. This rather complex method for calculating the Stokes number has the advantage that it automatically interpolates between the different drag regimes and is therefore more stable.

We use a symplectic, first-order leapfrog method to integrate Eqs. (12)–(15). Equations containing the Stokes number are integrated implicitly. First, Eqs. (12) and (13) are integrated using only half of a time step Δt. Omitting the particle index i, the integration from the kth to the (k + 1)th time step is performed as follows:

Then, Eqs. (14) and (15) are integrated with a full time step and using the new values and ,

where facc,r and facc,θ are the gravitational accelerations in r and θ directions due to the gravitational potential of the disc, extracted from FARGO-ADSG.

Finally, Eqs. (12) and (13) are integrated again with half of a time step using Lk+1 and vr,k+1 in the following way:

To ensure the accuracy of this method, every quantity is linearly interpolated onto the location of the dust particle. To be stable, the particles perform several time steps between two FARGO-ADSG outputs. Between two outputs, the gas quantities are linearly interpolated in time. To be as precise as possible, two consecutive FARGO-ADSG snapshots should not lie too far apart. As a rule of thumb, ten or more snapshots per orbit at the inner boundary of the system are enough for this numerical scheme.

Furthermore, a random walk term can be added to simulate the turbulent diffusion of the dust. If activated, after every time step, the particles are

displaced bya length lturb in a randomdirection. The distance of displacement is given by Youdin & Lithwick (2007) as (C.11)

This option is turned off in all the simulations presented here.

Appendix D Benchmark of the dust evolution code

To prove that the equations of motions are integrated accurately, we performed simulations of several test scenarios. In all cases, we kept the gas parameters constant with a surface density of and a temperature profile of . We consider a disc orbiting a 1 M star. For anaxisymmetric disc, the radial evolution of a given dust particle of index i is governed by (D.1)

where uP is the maximum particle drift velocity caused by the sub-Keplerian gas (Birnstiel et al. 2016). In a fixed gas disc, the particle’s radial evolution is solely dependent on the particle’s Stokes number.

To check whether the code is stable and correct for a broad variety of Stokes numbers, we performed four test scenarios. In three cases we kept the Stokes number of the particles constant at 10−100, 10100, and 1. In the first case the particles are completely coupled to the gas and are therefore being accreted with the radial gas velocity vgas,r. In the second case the particles are completely decoupled from the gas, are on Keplerian orbits, and do not change their radial position. In the third case with St = 1 the particles experience the maximum radial-drift velocity of . We also performed a fourth test case where we did not fix the Stokes number at all. Here, the Stokes number will decrease with time as the particles drift to regions with higher surface densities.

Figure D.1 shows the radial evolution of selected test particles in each of the four scenarios with blue lines. We also directlyintegrated Eq. (D.1) for all cases and overplotted the expected particle tracks with broad grey lines. The top left shows a particle with a very small Stokes number. It is perfectly coupled to the gas and is accreted with the radial gas velocity. The particle on the top right has a very high Stokes number and is therefore on a Keplerian orbit. Its radial position does not change at all. The particle on the bottom left has a Stokes number of unity and is accreted with the maximum drift velocity. Due to our choice of the surface density and temperature profiles, this drift velocity is constant with radial distance from the star. The gas velocity can be neglected here. The particle on the bottom right has a variable Stokes number that decreases with decreasing distance from the star due to higher surface densities in the inner disc. Therefore, the particle’s radial velocity decreases slightly over time. The results from the particle simulation and the direct integration match perfectly. We therefore conclude that the particle code is stable for all Stokes numbers and returns the correct results.

thumbnail Fig. D.1

Test cases of the particle simulation. Comparison of the radial evolution for test particles in four different scenarios. The results obtained with our dust code are shown in blue, while those obtained through the direct integration of Eq. (D.1) are in grey. From top to bottom and from left to right panels: St = 10−100, St = 10100, St = 1, and variableStokes numbers. The gas disc does not evolve in any of these test cases.

Open with DEXTER


  1. Akhlaghi, M., & Ichikawa, T. 2015, ApJS, 220, 1 [NASA ADS] [CrossRef] [Google Scholar]
  2. Ayliffe, B. A., Laibe, G., Price, D. J., & Bate, M. R. 2012, MNRAS, 423, 1450 [NASA ADS] [CrossRef] [Google Scholar]
  3. Barge, P., & Sommeria, J. 1995, A&A, 295, L1 [NASA ADS] [Google Scholar]
  4. Baruteau, C., & Masset, F. 2008, ApJ, 678, 483 [NASA ADS] [CrossRef] [Google Scholar]
  5. Baruteau, C., & Zhu, Z. 2016, MNRAS, 458, 3927 [NASA ADS] [CrossRef] [Google Scholar]
  6. Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Benisty, M., Stolker, T., Pohl, A., et al. 2017, A&A, 597, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  9. Birnstiel, T., Dullemond, C. P., & Pinilla, P. 2013, A&A, 550, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Birnstiel, T., Fang, M., & Johansen, A. 2016, Space Sci. Rev., 205, 41 [NASA ADS] [CrossRef] [Google Scholar]
  11. Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Bohren, C. F., & Huffman, D. R. 1983, Absorption and Scattering of Light by Small Particles (Weinheim: Wiley-VCH Verlag GmbH) [Google Scholar]
  13. Casassus, S., Wright, C. M., Marino, S., et al. 2015, ApJ, 812, 126 [NASA ADS] [CrossRef] [Google Scholar]
  14. Christiaens, V., Casassus, S., Perez, S., van der Plas, G., & Ménard, F. 2014, ApJ, 785, L12 [NASA ADS] [CrossRef] [Google Scholar]
  15. Cuello, N., Gonzalez, J.-F., & Pignatale, F. C. 2016, MNRAS, 458, 2140 [NASA ADS] [CrossRef] [Google Scholar]
  16. Dai, F., Facchini, S., Clarke, C. J., & Haworth, T. J. 2015, MNRAS, 449, 1996 [NASA ADS] [CrossRef] [Google Scholar]
  17. Dipierro, G., Pinilla, P., Lodato, G., & Testi, L. 2015, MNRAS, 451, 974 [NASA ADS] [CrossRef] [Google Scholar]
  18. Dong, R., Zhu, Z., Rafikov, R. R., & Stone, J. M. 2015, ApJ, 809, L5 [NASA ADS] [CrossRef] [Google Scholar]
  19. Dorschner, J., Begemann, B., Henning, T., Jaeger, C., & Mutschke, H. 1995, A&A, 300, 503 [Google Scholar]
  20. Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, Astrophysics Source Code Library [record ascl:1202.015] [Google Scholar]
  21. Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Facchini, S., Juhász, A., & Lodato, G. 2018, MNRAS, 473, 4459 [NASA ADS] [CrossRef] [Google Scholar]
  23. Follette, K. B., Grady, C. A., Swearingen, J. R., et al. 2015, ApJ, 798, 132 [NASA ADS] [CrossRef] [Google Scholar]
  24. Fouchet, L., Maddison, S. T., Gonzalez, J.-F., & Murray, J. R. 2007, A&A, 474, 1037 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics, 3rd edn (Cambridge: Cambridge University Press), 398 [Google Scholar]
  26. Fromang, S., Lyra, W., & Masset, F. 2011, A&A, 534, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Garufi, A., Quanz, S. P., Schmid, H. M., et al. 2016, A&A, 588, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  29. Gonzalez, J.-F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984 [NASA ADS] [Google Scholar]
  30. Grady, C. A., Muto, T., Hashimoto, J., et al. 2013, ApJ, 762, 48 [NASA ADS] [CrossRef] [Google Scholar]
  31. Jacquet, E., Balbus, S., & Latter, H. 2011, MNRAS, 415, 3591 [NASA ADS] [CrossRef] [Google Scholar]
  32. Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  33. Johansen, A., Youdin, A. N., & Lithwick, Y. 2012, A&A, 537, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Kratter, K., & Lodato, G. 2016, ARA&A, 54, 271 [NASA ADS] [CrossRef] [Google Scholar]
  35. Krauss, O., & Wurm, G. 2005, ApJ, 630, 1088 [NASA ADS] [CrossRef] [Google Scholar]
  36. Kretke, K. A., & Lin, D. N. C. 2007, ApJ, 664, L55 [NASA ADS] [CrossRef] [Google Scholar]
  37. Laibe, G., Gonzalez, J.-F., & Maddison, S. T. 2012, A&A, 537, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Marino, S., Perez, S., & Casassus, S. 2015, ApJ, 798, L44 [NASA ADS] [CrossRef] [Google Scholar]
  39. Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
  40. McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Data Analysis Software and Systems XVI, eds. R. A. Shaw, F. Hill, & D. J. Bell, ASP Conf. Ser., 376, 127 [NASA ADS] [Google Scholar]
  41. Meheut, H., Meliani, Z., Varniere, P., & Benz, W. 2012, A&A, 545, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Meru, F., Juhász, A., Ilee, J. D., et al. 2017, ApJ, 839, L24 [NASA ADS] [CrossRef] [Google Scholar]
  43. Montesinos, M., & Cuello, N. 2018, MNRAS, 475, L35 [NASA ADS] [CrossRef] [Google Scholar]
  44. Montesinos, M., Perez, S., Casassus, S., et al. 2016, ApJ, 823, L8 [NASA ADS] [CrossRef] [Google Scholar]
  45. Muto, T., Grady, C. A., Hashimoto, J., et al. 2012, ApJ, 748, L22 [NASA ADS] [CrossRef] [Google Scholar]
  46. Nakagawa, Y., Sekiya, M., & Hayashi, C. 1986, Icarus, 67, 375 [NASA ADS] [CrossRef] [Google Scholar]
  47. Paardekooper, S.-J. 2007, A&A, 462, 355 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Paardekooper, S.-J., & Mellema, G. 2004, A&A, 425, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Paardekooper, S.-J., & Mellema, G. 2006, A&A, 453, 1129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Pérez, L. M., Carpenter, J. M., Andrews, S. M., et al. 2016, Science, 353, 1519 [NASA ADS] [CrossRef] [Google Scholar]
  51. Pfalzner, S. 2003, ApJ, 592, 986 [NASA ADS] [CrossRef] [Google Scholar]
  52. Pignatale, F. C., Gonzalez, J.-F., Cuello, N., Bourdon, B., & Fitoussi, C. 2017, MNRAS, 469, 237 [NASA ADS] [CrossRef] [Google Scholar]
  53. Pinilla, P., Benisty, M., & Birnstiel, T. 2012, A&A, 545, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Ribas, Á., Merín, B., Bouy, H., & Maud, L. T. 2014, A&A, 561, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Rice, W. K. M., Lodato, G., Pringle, J. E., Armitage, P. J., & Bonnell, I. A. 2004, MNRAS, 355, 543 [NASA ADS] [CrossRef] [Google Scholar]
  56. Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A, 410, 611 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
  58. Stolker, T., Dominik, C., Avenhaus, H., et al. 2016, A&A, 595, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  59. Tomida, K., Machida, M. N., Hosokawa, T., Sakurai, Y., & Lin, C. H. 2017, ApJ, 835, L11 [NASA ADS] [CrossRef] [Google Scholar]
  60. van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  61. van der Marel, N., Cazzoletti, P., Pinilla, P., & Garufi, A. 2016, ApJ, 832, 178 [NASA ADS] [CrossRef] [Google Scholar]
  62. Vinković, D. 2014, A&A, 566, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Wagner, K., Apai, D., Kasper, M., & Robberto, M. 2015, ApJ, 813, L2 [NASA ADS] [CrossRef] [Google Scholar]
  64. Weidenschilling, S. J. 1977, MNRAS, 180, 57 [NASA ADS] [CrossRef] [Google Scholar]
  65. Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67 [NASA ADS] [CrossRef] [Google Scholar]
  66. Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [NASA ADS] [CrossRef] [Google Scholar]
  67. Youdin, A. N., & Lithwick, Y. 2007, Icarus, 192, 588 [NASA ADS] [CrossRef] [Google Scholar]
  68. Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Zsom, A., Sándor, Z., & Dullemond, C. P. 2011, A&A, 527, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]


As a rule of thumb, we require at least ten hydro snapshots per orbit at the inner boundary.


From both larger and smaller radii.

All Tables

Table 1

Detections of spiral and asymmetric dust traps in transition discs.

All Figures

thumbnail Fig. 1

Sketch of the inner regions of the disc model considered in this work. The outer disc is face-on, while the non-precessing inner disc is edge-on. The light-obstructing material around the star casts shadows at two diametrically opposed regions of the outer disc, in dark blue.

Open with DEXTER
In the text
thumbnail Fig. 2

Snapshots of the mid-plane gas temperature (top panels) and the vertically integrated pressure (bottom panels) at four different snapshots of the Sha-r100 simulation –from left to right panels: 2500, 5000, 7500, and 10 000 yr.

Open with DEXTER
In the text
thumbnail Fig. 3

(r, θ) maps of the vertically integrated pressure in NoSha-r100 (P0, left panel), in Sha-r100 (P, middle panel), and the normalised perturbation amplitude (PP0)∕P0 (right panel). All the snapshots are taken after 10 000 yr of evolution.

Open with DEXTER
In the text
thumbnail Fig. 4

Gas and dust distributions in the (r,θ)-plane after 7500 (left panels) and 10 000 (right) yr of evolution. Top row: vertically integrated pressure map of the disc. The rows below show the dust distribution coloured with the logarithm of the Stokes number – from top to bottom panels: 100 μm, 1 mm, 1 cm.Spirals in the gas and in the dust have the same morphology.

Open with DEXTER
In the text
thumbnail Fig. 5

From top to bottom panels: radial and azimuthal positions, radial velocity, and Stokes numbers as a function of time for three groups of particles, which end at r = 300± 5 au and θ = 0 ± 5. We plot 100 μm particles in blue, 1 mm particles in green, and 1 cm particles in red. Each group contains three particles of equal size.

Open with DEXTER
In the text
thumbnail Fig. 6

Same as Fig. 5, but for r = 500 ± 5 au.

Open with DEXTER
In the text
thumbnail Fig. 7

Dust thermal emission at 1.6 μm (left panel) and 1.3 mm (right panel) after 10 000 yr obtained with RADMC-3D. The source is placed at d = 140 pc and the disc is seen almost face-on (i = 13°).

Open with DEXTER
In the text
thumbnail Fig. 8

Synthetic observations with ALMA of Sha-r100 after 10 000 yr. The beam is indicated in the bottom left corner.

Open with DEXTER
In the text
thumbnail Fig. B.1

Initial and final radial temperature profiles for the Sha-r100 simulation. The range of temperature is consistent with Eq. (9).

Open with DEXTER
In the text
thumbnail Fig. B.2

Initial and final radial Toomre parameter profiles for the Sha-r100 simulation. Values above 1 correspond to gravitationally stable discs. This is the case for the disc model considered in this work.

Open with DEXTER
In the text
thumbnail Fig. D.1

Test cases of the particle simulation. Comparison of the radial evolution for test particles in four different scenarios. The results obtained with our dust code are shown in blue, while those obtained through the direct integration of Eq. (D.1) are in grey. From top to bottom and from left to right panels: St = 10−100, St = 10100, St = 1, and variableStokes numbers. The gas disc does not evolve in any of these test cases.

Open with DEXTER
In the text

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

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

Initial download of the metrics may take a while.