© ESO, 2017

1. Introduction

The trapping of solid particles in gaseous vortices and its consequences on planet formation is an issue investigated today by many authors. Vortices in protoplanetary disks can be produced by various mechanisms such as inverse cascade in 2D geometries or specific 2D/3D instabilities like the Rossby wave or the baroclinic instability. Three-dimensional studies have confirmed that such vortices can form and survive in spite of the elliptical instability that predominantly affects the core of the vortices (Lesur & Papaloizou 2010; Richard et al. 2013). Ultimately, they can be depicted as quasi-2D structures with a nearly stationary evolution. In 2D Keplerian flows, vortices are known to capture very effectively the solid particles embedded in the gas (Barge & Sommeria 1995). The trapping mechanism is so efficient that particle concentration can be strong enough for significant back-reaction of the solid particles onto the gas. Dust feedback enables the growth of instabilities like Rayleigh-Taylor (Chang & Oishi 2010), Kelvin-Helmoltz (Johansen et al. 2006), or streaming instabilities (Youdin & Johansen 2007). In vortices, the confinement of the solid particles has recently been addressed by Fu et al. (2014), Raettig et al. (2015), Crnkovic-Rubsamen et al. (2015), Baruteau & Zhu (2016), Surville et al. (2016).

Here, our paper focuses on the global distribution of the solid particles in a disk hosting a persistent gaseous vortex but does not study the detailed complex evolution inside the vortex core which will be addressed elsewhere. Our goal is to explore how a giant vortex can impact the dust surface density of a protoplanetary disk at a global scale.

The most appropriate place for the formation and survival of these vortices are the optically thick regions of the disk where MRI is ineffective due to the weak ionization levels of the gas. For example, the boundary between the dead-zone and the turbulent inner regions has attracted a special interest since it can be the place of a density enhancement with the possible growth of the Rossby wave instability (Lovelace et al. 1999; Li et al. 2001; Varnière & Tagger 2006).

A study of these disk outer regions is of additional interest since the data can be confronted to circumstellar disk observations, particularly with the new generation of high-resolution instruments. Interferometric observations at submillimeter to centimeter wavelengths probe material in the disk with enough angular resolution to highlight asymmetries due to the presence of vortices. Images with a high signal-to-noise ratio in spectral line typically probe molecular gas in the disk upper layers because of the high optical depths, while imaging of the dust thermal emission allows the investigation of the spatial distribution of dust particles in the disk midplane, where vortices are expected to form. Strong azimuthal asymmetries in the dust thermal emission have been observed in a number of young disks at submillimeter and millimeter wavelengths, e.g. Oph IRS 48 (van der Marel et al. 2013, 2015), LkHα 330 (Isella et al. 2013), HD 142527 (Casassus et al. 2013), SAO 206462, SR 21 (Pérez et al. 2014). Such asymmetries are often attributed to dust trapped in large-scale vortices that are able to collect large amounts of dust grains with sizes of the order of ~1 mm. Another indirect way to detect vortices in disks is the fortuitous observation of eclipses of a young star by the trapped dust in a highly inclined disk (Barge & Viton 2003; Rostopchina-Shakhovskaja et al. 2012).

In this paper, numerical simulations are performed to follow the evolution of dust particles embedded in a gas disk that contains a vortex. The solid particles are aerodynamically coupled to the gas and the presence of the vortex is found to significantly change their global distribution around the star. Depending on their size, the particles are either captured and concentrated inside the vortex (Barge & Sommeria 1995; de la Fuente Marcos & Barge 2001; Chavanis 2000) or repelled from the vortex orbital radius where a gap forms. The carving of a gap in the dust distribution is a forgotten consequence of the dust-vortex interaction that we qualitatively present and discuss in this work. The results of the numerical simulations are presented in the form of surface density maps to make comparisons with the observations easier. They are also illustrated and corroborated by computing the trajectories of passive solid particles in a steady velocity field. Global simulations and trajectory computations were performed for a wide range of particle sizes.

The paper is presented as follows: Sect. 2 describes the disk-vortex model and the bi-fluid numerical simulations; Sect. 3 presents the main results of the numerical simulations and the density maps; Sect. 4 describes the dust vortex interactions and explains the different structures observed in the dust distribution; Sect. 5 presents simulated observations derived from the density and temperature maps using ALMA and Next Generation Very Large Arrays (ngVLA) simulators. The results of the simulated observations are discussed in Sect. 6, and the main results of this investigation are summarized in Sect. 7.

2. Models and numerical simulations

To study the evolution of the dust spatial distribution we performed 2D numerical simulations with a bi-fluid (gas+dust) code and use a simple model of a protoplanetary disk in the following framework:

• (i)

  the disk extends from 20 to 100 AU from thecentral star, a region large enough to provide predictions that canbe compared with the observations;

• (ii)

  the surface density and the temperature of the background structure of the disk, i.e. before a vortex is introduced, are given by simple power laws: Σ ∝ r^{− p} and T ∝ r^{− q}, respectively;

• (iii)

  the gas is assumed to be a mixture of hydrogen and helium with a mean molecular weight μ = 2.34m_H where m_H is the mass of a hydrogen atom;

• (iv)

  the dust is described as spherical particles of varying sizes (0.04 mm < s < 2 cm) and with a volume density of ρ_° = 0.8 g cm^-3, similar to the dust model considered by Tazzari et al. (2016);

• (v)

  the vortex is introduced at the beginning of the simulation in the form of an approximate solution of the incompressible Euler equations.

2.1. Disk parameters

The disk model we used discards the standard MMSN (Minimum Mass Solar Nebula) assumption and rests on recent observations of protoplanetary disks. We have chosen the primordial disk around the young star Oph IRS 48 as our reference case. This disk shows the strongest asymmetry in the azimuthal distribution of dust particles in a disk observed to date, which could be the signature of mm- and cm-sized dust trapped by a large-scale vortex (van der Marel et al. 2013, 2015; Lyra 2014; Zhu et al. 2014).

The required parameters for the disk model are then deduced from the observations of this peculiar object. The analysis of the ALMA observations of the dust continuum emission and CO isotopologues allowed Bruderer et al. (2014) to derive a model for the radial profile of the gas distribution in the regions out of the vortex. Although it is not our goal to find a model that reproduces in detail the observations for this disk, we use this case as a fiducial case to investigate the effects of a giant vortex on the global distribution of gas and dust in a real disk system.

Following Bruderer et al. (2014), the radial profile of the gas surface density can be described by a power-law function with index p = 1.1. We consider a gas surface density of 0.32 g cm^-2 at 60 AU from the star, which is the radial distance of the vortex structure as seen in the Oph IRS 48 disk (van der Marel et al. 2013). This is a factor of 10 greater than the value for the gas density considered in Bruderer et al. (2014) at the same location in the disk. We chose this larger value because subsequent radiative transfer modeling of disks has shown that masses derived from observations of CO isotopologues can increase by up to an order of magnitude when one properly treats the isotope-selective photodissociation of these molecules (Miotello et al. 2014, 2016).

On the other hand, the radial profile of the gas temperature is a power-law with an index of q = 0.56, with a value of 100 K at 60 AU (top right panel of Fig. 8 in Bruderer et al. 2014). The vertical scale height of the disk is derived from hydrostatic equilibrium at this temperature and the gravitational potential is for a 2 solar-mass star.

For the particles, different sizes are considered in the present work. We have chosen grain sizes that range from 0.04 mm to 6.4 mm to predict maps for the dust continuum emission at different wavelengths between λ = 0.45 mm and 0.9 cm, which can be probed at high angular resolution and sensitivity with the ALMA and future ngVLA interferometers. Although the observed emission realistically comes from a combination of grains spanning several orders of magnitude in sizes, the most efficient emitting grains have sizes s ≈ λ/ (2π | m |), where m is the refractive index of the emitting dust (m ≈ 2−3 for dust in protoplanetary disks, e.g. Rodmann et al. 2006).

The simulations performed in this paper allow us to explore the dependance of the dust-gas friction on the dust-vortex interaction and its consequence on the dust surface density and thermal emission. The initial distribution of the dust is assumed to trace the distribution of the gas with a dust-to-gas mass ratio equal to 0.001 across the whole disk (simulations performed with a mass ratio of 0.01 are also presented in Appendix A).

2.2. Vortex model

In this work we disregard the problem of vortex formation in a disk which has been addressed elsewhere (e.g. Richard et al. 2013; Barge et al. 2016). Hence, the vortex is introduced at the beginning of the simulations and, in a few rotations, adjusts to a quasi-steady vortex structure that lasts long enough to significantly perturb the motion of the dust particles. Two dimensional vortices are, indeed, quasi-steady gaseous structures that can survive over very long time periods of up to thousands of vortex rotations around the star (Surville & Barge 2015). Their lifetimes in numerical simulations are shorter when dust is included in the disk due to the gas drag, but the present simulations show that the vortices can survive long enough to significantly change the dust surface density and to imprint observable structures in the disks.

For the initial conditions we used a simple model of Rossby vortex developed by Surville & Barge (2015). The control parameters of the vortices were selected to get vortices in the “giant vortex class” in order to mimic the large and strong dust asymmetries observed in a number of disks (e.g. Casassus et al. 2013; van der Marel et al. 2013; Marino et al. 2015). In this model the gaseous vortex is centered on (r₀, θ₀) in polar coordinates and is described with a pseudo enthalpy ℋ(r,θ) approximated by a 2D Gaussian function parametrized by an amplitude ℋ₀ and two widths δr = r₀χ_r and δθ = χ_rχ_θ, where χ_r is the normalized radial extent of the vortex and χ_θ is its azimuthal aspect ratio. This function is maximum at the vortex center and smoothly connects to zero at the background level.

With this Gaussian approximation we found simple expressions for the velocity field of a gaseous vortex when considered in a frame of reference rotating at the mean angular velocity of the flow Ω_o$\begin{array}{ccc}\mathit{u}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{Ro}}\frac{\mathit{\theta }\mathrm{-}{\mathit{\theta }}_{\mathrm{0}}}{\mathit{\delta }{\mathit{\theta }}^{\mathrm{2}}}\frac{\mathrm{\mathscr{H}}}{\mathit{r}{\mathrm{\Omega }}_{\mathrm{K}}}\\ \mathit{v}& \mathrm{=}& \left[\mathrm{-}\mathit{r}\left(\frac{\mathit{r}\mathrm{-}{\mathit{r}}_{\mathrm{0}}}{\mathit{\delta }{\mathit{r}}^{\mathrm{2}}}\right)\mathrm{+}\frac{\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\gamma }\mathrm{)}\mathit{p}\mathrm{+}\mathit{q}}{\mathrm{2}\mathit{\gamma }}\right]\frac{\mathrm{\mathscr{H}}}{\mathit{r}{\mathrm{\Omega }}_{\mathrm{K}}}\end{array}$where Ro is the Rossby number at the vortex center and γ = 1.4 is the adiabatic index of the gas (the second term in the bracket is zero for an isentropic vortex). At the beginning of the simulations the vortex is characterized by Ro = 0.1, χ_r = 0.15 and χ_θ = 6.9. This structure is added on top of the equilibrium state of the disk, centered at a distance of 60 AU from the star. The evolution of gas and dust is then provided by a numerical integration of standard fluid equations (see, e.g., Inaba et al. 2005).

2.3. Numerical simulations

The simulations are performed with a bi-fluid code where the solid component is considered a fluid without pressure. This approximation requires that the Stokes number and dust-to-gas mass ratio be smaller than one (Garaud & Lin 2004). In our problem the two conditions are satisfied except (i) when the particle size is larger than one centimeter or (ii) when located in the core of the vortex and close to its orbital radius where dust strongly accumulates. The code was previously developed to study the dynamical evolution of protoplanetary disks (Inaba et al. 2005; Surville & Barge 2015). It solves the inviscid continuity and Euler equations vertically integrated over the disk thickness. At equilibrium, the gas rotates at slightly less than the Keplerian velocity, due to the radial pressure gradient, with an angular velocity ${\mathrm{\Omega }}_{\mathrm{o}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{K}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\left(\mathrm{1}\mathrm{-}\mathit{\eta }\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}$(2)where $\mathit{\eta }\mathrm{=}\frac{\mathit{p}\mathrm{+}\mathit{q}}{\mathit{\gamma }}{\mathit{M}}_{\mathrm{A}}^{-2}\mathrm{\left(}\mathit{r}\mathrm{\right)}$ and M_A is the Mach number; the flow is also assumed compressible and adiabatic but non-homentropic. The solid particles are fully coupled to the gas by aerodynamical drag forces in the Epstein’s regime and characterized by a non-dimensional parameter or Stokes number $\mathit{St}\mathrm{=}\frac{\mathit{\pi }}{\mathrm{4}}\frac{\mathit{s}{\mathit{\rho }}_{\mathrm{◦}}}{\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{·}$(3)In a gas disk at equilibrium a solid particle drifts inward to the star due to the pressure gradient; its radial velocity is approximately given by ${\mathit{u}}_{\mathrm{dust}}\mathrm{=}\frac{\mathit{St}}{\mathrm{1}\mathrm{+}\mathit{S}{\mathit{t}}^{\mathrm{2}}}\mathit{\eta r}{\mathrm{\Omega }}_{\mathrm{K}}\mathit{,}$(4)which displays a maximum at St = 1 (e.g. Nakagawa et al. 1986; Inaba et al. 2005).

We start the simulations from a two-component disk model at equilibrium on top of which the vortex solution presented in the above model is added. The simulations are performed in a whole ring around the star, large enough to get significant changes in the radial distribution of the dust.

The code was previously used to study the formation of vortices by the Rossby-Wave or the baroclinic instability (Richard et al. 2013; Barge et al. 2016), and also to study the coupled evolution of gas and dust (Inaba & Barge 2006). This code, which uses the finite volume method and was first presented in Inaba et al. (2005) can maintain high accuracy over large integration periods (Surville & Barge 2015; Surville et al. 2016).

Test simulations show a three-step evolution: (i) the approximate vortex-solution adjusts in the flow to a vortex that numerically satisfies the fluid equations in ~10 vortex rotations; (ii) the formed vortex survives in the disk for hundreds of vortex rotations; (iii) the vortex is progressively destroyed by the back-reaction of the solid particles on a longer time scale.

The total duration of the integrations necessary for our simulations was determined by a compromise between the vortex lifetime and the time necessary for the vortex to capture and trap a significant amount of solid particles. It was then empirically set to 100 orbital periods of the vortex around the star.

The numerical resolution is (2048 × 4080) in the (r,θ) directions, sufficient to correctly describe particle captures by the vortex. In the radial direction, the resolution corresponds to 274 cells by scale height at 60 AU from the star. Simulations were performed for different values of the particle size.

3. Main results of the numerical simulations

The simulations provide the spatial distribution of the dust particles after 100 rotations of the vortex around the star. The results are discussed as a function of the particle size or the Stokes number estimated at the radial position of the vortex.

Fig. 1 Dust density in a disk containing a giant gaseous vortex that has completed 100 rotations around the star (units are in g cm-2). The vortex is initially located at 60 AU from the star; in each plot the solid particles have a single size: 0.046 mm, 0.092 mm, 0.14 mm, 0.32 mm, 0.74 mm, 1 mm, 2 mm, 6.4 mm, 2 cm, from left to right and from top to bottom.

3.1. Global distribution of the dust

Figure 1 shows the density map of the dust for different values of the particle size where we note a number of striking features.

• (1)

  Submillimeter-sized particles (St ≲ 0.063) concentrate in a regularlyshaped region centered on the vortex core (as commonlyexpected for passive particles).

• (2)

  Millimeter-sized particles (0.063 ≲ St ~ 1) concentrate in an irregular and inhomogeneous core region (containing dust knots), while a gap is carved just inside the vortex orbital radius (the gap depth increases with particle size).

• (3)

  Centimeter-sized particles (1 ≲ St) are strongly depleted inside the vortex orbital-radius, but concentrate in a narrow and knotted structure (arc or ring) where the dust-to-gas mass ratio can overcome unity (see Figs. 3 and 4).

In this last case we deliberately ignore some of our initial assumptions to explore weaker coupling with the gas avoiding complex developments or the use of a different code; even so, we note a piling up of the solid particles at the inner boundary of the simulated ring, due to the decrease in the drift velocity with the distance to the star. A similar effect occurs at the outer boundary for particle sizes of 0.14 mm and 0.32 mm but could be the result of a numerical artefact.

3.2. Evolution of the dust particles

Our simulations also provide a step-by-step description of the dynamical evolution of the dust particles. Three different cases have been determined following the values of the Stokes number.

• Light particles (St< 1). These particles are directly captured by thevortex and confined for a long time in its core where the dust-to-gas ratio can overcome unity allowing a significant back-reactionof the particles onto the gas. Figure 2 shows theconcentration of millimeter-sized particles in the core of thevortex with the formation of a high-density knot that becomesturbulent under instabilities in the dust-gas mixture. After about40 rotations the dust surface density undergoes a strong jump atthe orbital radius of the vortex and the gap that is formed presents alayered structure. Ridges correspond to a spiral lane of non-captured particles that are slowly drifting inward toward the star. The density contours of the gas are plotted on the same figure and indicate the position of the vortex (concentric elliptical paths) and the position of the density waves associated with the vortex (oblique rays at the top and bottom of the figure). The contours are noisy, but their large-scale deformations, after 100 rotations, reveal vortex stretching and turbulent motions due to the back-reaction of the particles onto the gas.

  Fig. 2 Dust-to-gas ratio for light particles of size s = 2 mm. The map is plotted in the (r,θ) plane and the evolution is shown after 20, 40 and 100 rotations of the vortex, from left to right and from top to bottom. Units of the x- and y-axis are in AU and radians, respectively.

• Optimal particles (St ~ 1). The capture efficiency is the greatest for particles with St ≃ 1. These particles rapidly concentrate in a small region of the vortex core (Fig. 3), which is turbulent due the back-reaction of the particles. In the inner part of the disk the gap is deeper than for lighter particles and also presents a layered structure. In this case, the back reaction of the particles onto the gas is strong enough to significantly change the vortex shape after 100 rotations.

  Fig. 3 Dust-to-gas ratio for nearly optimal particles after 20, 40 and 100 rotations of the vortex. The particle size is s = 6.4 mm. The layout is the same as in Fig. 2.

• Over-optimal particles (St ≳ 1). Particles in the centimeter range marginally violate our initial assumptions, but permit us to explore dust and gas evolution for weaker coupling. The vortex shape is less affected than for smaller particles (cf. Fig. 4). The particles focus on a very small region of the disk that looks like a knotted arc close to the orbital radius of the vortex. This structure evolves into a circular chain of high-density knots, each of which covers a surface on the order of 1 × 4 AU² or 25 × 40 numerical cells. Similar concentrations of particles were also observed in simulations performed by Surville et al. (2016) when St ~ 0.1. Inside this radius the surface density decreases dramatically, as in the previous case; a gap and its layered structure also form.

Fig. 4 Dust-to-gas ratio for a particle size of s = 2 cm after 20, 40, and 100 rotations of the vortex; layout as in Fig. 2.

4. Dust-vortex interaction

Our simulations show that the evolution of the solid particles and their distribution around the star strongly depends on the value of the Stokes number St. Particles are found either to be trapped and concentrated in the vortex or to form a high-density ring lining a strong density jump at the orbital-radius of the vortex. The capture by a vortex due to the local action of the Coriolis force is a well known mechanism. On the other hand, the carving of a gap with a strong density jump is a new output of the dust-vortex interaction that has to be explained. Here, we propose a number of carving mechanisms.

4.1. Capture and radial drift

In its motion around the star the vortex captures dust particles and concentrates them in its core. In a shearing sheet context particles are selected when their impact parameter is smaller than a threshold. There are also non-captured particles that can be divided into two classes: outward particles that were too far from the vortex orbital radius to be captured and inner particles that were too slow to catch the vortex.

In fact, due to the systematic inward drift, particles of the outward class clearly have a chance to be captured after a second turn of the vortex around the star; this obviously increases the capture efficiency of the vortex during its lifetime.

Our simulations show that particles of the inward class are not inevitably lost into the star but have a chance of being captured during additional turns of the vortex around the star. This mechanism concerns the particles that were not rapidly trapped by the vortex and are tracing a spiral wake in the disk.

Figure 5 shows the various steps of this mechanism with a dust lane that distorts gradually when approaching the vortex and that finally rolls up around the vortex core. This secondary trapping is observed throughout the vortex lifetime and systematically feeds the vortex with particles of the inner disk.

Fig. 5 Dust lane in the wake of the vortex that spirals around the star and is progressively captured by the vortex. The dust-to-gas ratio for particles with size s = 2 mm is plotted after 5, 7, 8, and 9 rotations of the vortex, from left to right and from top to bottom.

As a result, the systematic drift of the solid particles toward the star plays a major role in the carving of a gap in the dust distribution. On the one hand, it contributes to a strong depletion of the inner disk, due to the systematic drift toward the star and to the secondary trapping by the vortex. In contrast, in the outer disk, the depletion by the vortex can be balanced by the inward flux of solid particles that are crossing the outer boundary of the computational box. Indeed, in our simulations the outer boundary must be considered a permeable frontier for external particles since they are also drifting inward toward the star.

4.2. Collective drag

At a given stage, dust trapping becomes strong enough that collective effects and back-reaction on the gas come into play, and the dust concentration is dragged away and other concentrations form. They are in turn pulled away from the vortex and tend to form a knotted ring (see Fig. 4).

4.3. Trajectories of the dust particles

The dust surface density in the presence of a vortex can be also studied in terms of the dynamical evolution of dust particles. A simple approach to the problem is to use a vortex model as described in Sect. 2.2. In a frame of reference rotating at the angular velocity of the flow, the dynamical evolution of passive solid particles can be described by the motion equations: $\begin{array}{ccc}\begin{array}{c}\mathrm{¨}\\ \mathit{r}\end{array}\mathrm{-}\mathit{r}\mathit{\theta ̇}\mathrm{2}& \mathrm{=}& \mathrm{-}\frac{{\mathrm{\Omega }}_{\mathrm{K}}}{\mathit{St}}\mathrm{(}\mathit{ṙ}\mathrm{-}\mathit{u}\mathrm{)}\mathrm{+}\mathrm{2}{\mathrm{\Omega }}_{\mathrm{o}}\mathit{r}\mathit{\theta ̇}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{o}}^{\mathrm{2}}\mathit{r}\mathrm{-}\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{{\mathit{r}}^{\mathrm{2}}}\\ \mathit{r}\begin{array}{c}\mathrm{¨}\\ \mathit{\theta }\end{array}\mathrm{+}\mathrm{2}\mathit{ṙ}\mathit{\theta ̇}& \mathrm{=}& \mathrm{-}\frac{{\mathrm{\Omega }}_{\mathrm{K}}}{\mathit{St}}\mathrm{(}\mathit{r}\mathit{\theta ̇}\mathrm{-}\mathit{v}\mathrm{)}\mathrm{-}\mathrm{2}{\mathrm{\Omega }}_{\mathrm{o}}\mathit{ṙ}\mathit{,}\end{array}$where u and v are the two components of the gas velocity changes due to the presence in the flow of a persistent structure (vortex and/or spiral waves). These two motion equations are then solved numerically using a second-order Runge-Kutta method.

4.3.1. Vortex model

Particle trajectories were first computed in the vortex velocity-field described in Eqs. (1) (with no other contribution) and are plotted in Fig. 6. Two types of trajectories can be distinguished following the value of the impact parameter: vortex trapping (black line) or inward drift toward the star (green line). There is a steep change between the two. Testing various particle sizes, we also note that the transition between trapping and drift becomes sharper with increasing size; the threshold is located close the core orbital radius in accordance with the results of the numerical simulations.

Fig. 6 Trajectories of solid particles interacting with a gaseous vortex. The particle sizes are equal to 0.14 mm and 0.32 mm, from left to right.

The capture trajectories are similar to those computed with approximate vortex models (e.g. de la Fuente Marcos & Barge 2001).

4.3.2. With a spiral wave model

Vortices in Keplerian disks are always accompanied by spiral compressional waves that they themselves excite in the gas flow. These waves rotate around the star at the vortex angular-velocity and are also associated with velocity changes of the gas at the crossing of the spiral arms. Because they couple with the gas, particles have their motions perturbed at spiral-arm crossing and can exchange angular momentum with the waves, accelerating or decelerating depending on whether the orbit is located outside or inside the core orbital radius, respectively. The importance of this effect depends on the stopping time of the particles but also on the width, the amplitude and winding of the waves.

Fig. 7 Radial drift of solid particles interacting with spiral+vortex (time is in number of vortex rotations). Particle size is 0.32 mm and 6.4 mm, from top to bottom. The various curves are for different values of the impact parameter.

A detailed study of particle-wave interactions is beyond of the scope of the present work; our goal here is only to capture the main consequences of spiral waves on the dynamical evolution of dust particles coupled to the gas by a drag force. To this end we have taken advantage of the above computations by replacing the vortex model with a crude model of a spiral pattern. In our toy model the effect of the spiral wave reduces to perturbations of the gas velocity that we have taken in the form $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{u}\mathrm{=}{\mathit{ϵ}}_{\mathit{u}}\mathit{f}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathit{v}\mathrm{=}\mathrm{[}{\mathit{ϵ}}_{\mathit{v}}\mathrm{+}\mathit{r}\mathrm{(}{\mathrm{\Omega }}_{\mathrm{K}}\mathrm{-}{\mathrm{\Omega }}_{\mathrm{o}}\mathrm{\left)}\mathrm{\right]}\mathit{f}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}\end{array}& & \end{array}$(6)where ϵ_u and ϵ_v are the amplitude of the velocity perturbations in the radial and azimuthal directions, respectively, and f(θ) is an arbitrary cut-off function $\mathit{f}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\exp \left[\mathrm{-}{\left(\frac{\mathit{\theta }\mathrm{-}{\mathit{\theta }}_{\mathrm{s}}}{\mathit{\delta }{\mathit{\theta }}_{\mathrm{s}}}\right)}^{\mathrm{2}}\right]$(7)with θ_s ∝ (r/r₀−1)² to mimic the spiral shape of the perturbation.

Then, the motion equations (Eqs. (5)) are solved using the same method as in the above section. The results are reported in Fig. 7 for various values of the impact parameter.

The inward drift of the solid particles presents a sharp transition at the orbital radius of the vortex: the drift slows down outward while accelerates inward. The difference between the two drift rates could help to explain the formation of a plateau lined by a strong gap.

5. Simulation of the observations

In this section we discuss the observability of the structures predicted by the numerical simulations presented in the previous sections. We use the results of the numerical simulations for the different grain sizes to investigate the change in the appearance of the dust thermal emission at different wavelengths, as continuum observations at longer wavelengths are typically more sensitive to emission from larger grains (e.g. Pérez et al. 2012). For the smaller grains of our simulations (s< 1 mm) we simulate observations with ALMA, which provides the high angular resolution and sensitivity needed to trace the structures discussed in this work. In order to trace larger grains, which require observations at longer wavelengths, we simulate observations with the future ngVLA interferometer which will achieve high angular resolution and sensitivity at wavelengths longer than ~3 mm (Carilli et al. 2015), and will be particularly suitable for observations of large grains in protoplanetary disks (Isella et al. 2015). The surface brightness map of the dust thermal emission for our models was derived from the dust surface densities and temperatures provided by the numerical simulations. We set the absolute scale of each surface brightness map by requiring the total flux to be consistent with the measured sub-mm/mm spectral energy distribution of the Oph IRS 48 disk (van der Marel et al. 2015), which essentially fixes the values for the dust opacity, assumed constant throughout the disk, at each wavelength discussed here. The disk is seen with a face-on geometry.

Fig. 8 Simulated ALMA images computed from the disk simulations illustrated in Fig. 1. Each map corresponds to a single particle size: 0.14 mm (wavelength of the simulated observations of 0.43 mm), 0.32 mm (λ = 0.87 mm), 0.74 mm (λ = 1.3 mm), and 1 mm (λ = 2.7 mm), from left to right and from top to bottom.

5.1. ALMA observations

We used the online ALMA Observation Support Tool¹ (version 4.0) to simulate realistic ALMA observations for a subset of our simulations. For our models with grain sizes of 0.14, 0.32, 0.74 mm, and 1 mm, we simulated ALMA observations at wavelengths of 0.43 (ALMA Band 9), 0.87 (Band 7), 1.3 mm (Band 6), 2.7 mm (Band 3), respectively. We considered array configurations that provide an angular resolution of about 80 milliarcsec, or about 10 AU at the distance of Oph IRS 48, for all the ALMA simulations. We found that this angular resolution offers a good trade-off between angular resolution and surface brightness sensitivity for the structures of interest for this paper. An observing time of 5 h was assumed for each observation. Precipitable water vapor (PWV) levels of 0.5, 1.3 and 1.8 mm were adopted for the simulated observations at wavelengths of 0.43, 0.87 and 1.3 mm, respectively. Figure 8 shows the maps returned by the ALMA Observation Support Tool after performing deconvolution with the CASA CLEAN algorithm (Briggs weighting with robust parameter of 0.5).

5.2. ngVLA observations

For observations at longer wavelengths tracing larger grains, we considered the ngVLA interferometer, which is currently being designed as a future large radio facility operating at wavelengths between about 3 mm and 30 cm (Carilli et al. 2015). We simulated ngVLA observations at wavelengths of 3.7 mm (frequency of 80 GHz) and 1.0 cm (30 GHz) for our models with grain sizes of 6.4 mm and 2.0 cm, respectively.

For the ngVLA simulations we used the SIMOBSERVE task in CASA adopting an array configuration made of 300 antennas with longest baselines of 300 km and an inner core to recover the emission from the largest scales of the disk. Thermal noise was added using the CASA tool SETNOISE and was scaled to provide a final rms noise of about 0.11 and 0.040 μJy beam^-1 in the simulated images at 3.7 mm and 1.0 cm, respectively, and for the highest angular resolution at these frequencies, i.e. 5 and 12 mas. These noise values would correspond to a total integration time of about 100 hours with the current ngVLA specifications. Deconvolution was performed using CLEAN with Briggs weighting with a robust parameter of 0.5. The obtained images are shown in Fig. 9².

Fig. 9 Simulated ngVLA images computed from the disk simulations illustrated in Fig. 1. Each map corresponds to a single particle size: 6.4 mm (λ = 3.7 mm) and 2.0 cm (λ = 1.0 cm), in top and bottom rows, respectively. The images in the left panels were obtained after tapering the interferometric visibilities with a 0.1 arcsec beam; the images in the right panels fully exploit the maximum angular resolution expected for the ngVLA at these wavelengths (see text).

Fig. 10 Dust density and ALMA images simulated in a disk containing a vortex with an aspect ratio of 12. Top: density maps corresponding to a single population of solid particles with size equal to 0.046 mm and 0.74 mm, from left to right, respectively. Bottom, simulated ALMA images at a wavelength of 0.45 mm and 1.3 mm, from left to right, respectively.

6. Discussion of the simulated observations

We studied the evolution of solid particles in a protoplanetary disk that hosts a persistent gaseous vortex, and listed the various structures formed in the global dust distribution. One of our main goals was to identify the possible observational signatures of these structures to better assess the presence of a gaseous vortex in the disk. To this end we performed numerical simulations of a protoplanetary disk using ρ Oph IRS 48 as a reference, as this disk is suspected to host a large-scale gaseous vortex.

We found that, depending on their size (or Stokes number), the solid particles are either captured and confined in the vortex forming a high density region in the core where dust-gas instabilities can grow or form high density dust blobs that align in a knotted structure (irregular arc or ring) close to the orbital radius of the vortex core.

Inside this radius the dust surface density decreases dramatically forming a gap with a sharp edge.

Simulated images produced with the ALMA and the future ngVLA simulators show that a number of these large-scale structures could be observed.

The most striking features in the simulated ALMA images are (i) the high-density region in the vortex core that is clearly observable at all the wavelengths investigated here; (ii) the depletion of solid particles around the vortex and the piling up at the inner and outer boundaries of the disk that are particularly visible at λ = 0.43 and 0.87 mm, respectively; and (iii) the sharp density jump at the orbital radius of the vortex that is distinguishable mostly at λ = 1.3 mm, whereas at λ = 2.7 mm only the bright dust emission from particles trapped by the vortex and from the inner disk is detected, due to the lower ALMA sensitivity to dust emission at longer wavelengths (Fig. 8).

At each of the investigated ALMA observing bands, the peak of the surface brightness map is always found at the location of the dust trapped by the vortex. This indicates that only observations with enough sensitivity and angular resolution can highlight the presence of substructures that are due to the interaction between the vortex and the disk.

Furthermore, the azimuthal variation of the location of trapped particles with different sizes is clearly seen from the maps at different wavelengths. Figure 10 shows the density maps and the simulated ALMA images obtained in the case of a vortex with an aspect ratio two times larger than in the standard case studied above. The azimuthal asymmetries produced in the simulated images clearly have a larger azimuthal extent that could better mimic true observations.

At longer wavelengths, the high-resolution and sensitivity observations achievable with the ngVLA have the potential to probe the high-density regions in the vortex as well as the knotted ring structures expected for pebbles at λ = 3.7 mm and λ = 1 cm (Fig. 9).

7. Summary

We have shown the results of 2D simulations obtained with our bi-fluid code to calculate the spatial distribution of solid particles with different sizes as expected in a disk hosting a giant gas vortex. Our simulations show that a giant vortex not only captures dust grains with Stokes number St< 1, but can also significantly affect the distribution of larger grains (St ≳ 1) by carving a gap and by producing a narrow ring-like distribution of pebbles due to dust concentrations escaping from the vortex core. We used the results of our code to calculate the expected surface brightness map at different sub-mm to cm wavelengths, tracing emission from solids with different sizes. We then use ALMA and ngVLA simulators to test the observability of the substructures predicted by our code using high angular resolution and sensitivity observations in dust continuum with these sub-mm and radio interferometers. Most of the substructures predicted in this work could be observed for a disk with properties similar to those derived for the ρ Oph IRS 48 disk. However, this is not the case with the currently available observations of this object.

Acknowledgments

Computations were performed using HPC resources from GENCI [TGCC and CINES] (Grant – x2016047407) and also at LAM on the mib MPI cluster maintained by CESAM. We would also like to thank the anonymous referee for the useful comments.

References

Appendix A: Stronger dust-to-gas ratio

Simulations starting with a dust-to-gas ratio equals to 0.01 (i.e., the cosmic abundance) were also performed in the same context as in Sect. 2. To put all the heavy material of a disk into a single population of solid particles is, of course, an unrealistic assumption that was considered only as a limiting case. Moreover, this is not a favorable assumption for maintaining a vortex against friction drag and dust-gas instabilities. Figure A.1 shows an example of a simulation where the particle size is equal to 6.4 mm. The evolution of the dust distribution is similar to that observed in Fig. 4 with shorter time scales and different numerical values. We estimated that it does not make sense to simulate observations in this case.

Fig. A.1 Dust-to-gas ratio for nearly optimal particles after 20, 40 and 100 rotations of the vortex. The particle size is s = 6.4 mm. The layout is the same as in Fig. 2.

All Figures

	Fig. 1 Dust density in a disk containing a giant gaseous vortex that has completed 100 rotations around the star (units are in g cm-2). The vortex is initially located at 60 AU from the star; in each plot the solid particles have a single size: 0.046 mm, 0.092 mm, 0.14 mm, 0.32 mm, 0.74 mm, 1 mm, 2 mm, 6.4 mm, 2 cm, from left to right and from top to bottom.
In the text

	Fig. 2 Dust-to-gas ratio for light particles of size s = 2 mm. The map is plotted in the (r,θ) plane and the evolution is shown after 20, 40 and 100 rotations of the vortex, from left to right and from top to bottom. Units of the x- and y-axis are in AU and radians, respectively.
In the text

	Fig. 3 Dust-to-gas ratio for nearly optimal particles after 20, 40 and 100 rotations of the vortex. The particle size is s = 6.4 mm. The layout is the same as in Fig. 2.
In the text

	Fig. 4 Dust-to-gas ratio for a particle size of s = 2 cm after 20, 40, and 100 rotations of the vortex; layout as in Fig. 2.
In the text

	Fig. 5 Dust lane in the wake of the vortex that spirals around the star and is progressively captured by the vortex. The dust-to-gas ratio for particles with size s = 2 mm is plotted after 5, 7, 8, and 9 rotations of the vortex, from left to right and from top to bottom.
In the text

	Fig. 6 Trajectories of solid particles interacting with a gaseous vortex. The particle sizes are equal to 0.14 mm and 0.32 mm, from left to right.
In the text

	Fig. 7 Radial drift of solid particles interacting with spiral+vortex (time is in number of vortex rotations). Particle size is 0.32 mm and 6.4 mm, from top to bottom. The various curves are for different values of the impact parameter.
In the text

	Fig. 8 Simulated ALMA images computed from the disk simulations illustrated in Fig. 1. Each map corresponds to a single particle size: 0.14 mm (wavelength of the simulated observations of 0.43 mm), 0.32 mm (λ = 0.87 mm), 0.74 mm (λ = 1.3 mm), and 1 mm (λ = 2.7 mm), from left to right and from top to bottom.
In the text

Fig. 9 Simulated ngVLA images computed from the disk simulations illustrated in Fig. 1. Each map corresponds to a single particle size: 6.4 mm (λ = 3.7 mm) and 2.0 cm (λ = 1.0 cm), in top and bottom rows, respectively. The images in the left panels were obtained after tapering the interferometric visibilities with a 0.1 arcsec beam; the images in the right panels fully exploit the maximum angular resolution expected for the ngVLA at these wavelengths (see text).

In the text

	Fig. 10 Dust density and ALMA images simulated in a disk containing a vortex with an aspect ratio of 12. Top: density maps corresponding to a single population of solid particles with size equal to 0.046 mm and 0.74 mm, from left to right, respectively. Bottom, simulated ALMA images at a wavelength of 0.45 mm and 1.3 mm, from left to right, respectively.
In the text

	Fig. A.1 Dust-to-gas ratio for nearly optimal particles after 20, 40 and 100 rotations of the vortex. The particle size is s = 6.4 mm. The layout is the same as in Fig. 2.
In the text