# Gaps, rings, and non-axisymmetric structures in protoplanetary disks: Emission from large grains

J. Ruge, M. Flock, S. Wolf, N. Dzyurkevich, S. Fromang, Th. Henning, H.

Klahr, H. Meheut

► To cite this version:

J. Ruge, M. Flock, S. Wolf, N. Dzyurkevich, S. Fromang, et al. Gaps, rings, and non-axisymmetric structures in protoplanetary disks: Emission from large grains. Astronomy & Astrophysics - A&A, 2016, 590, pp.A17. 10.1051/0004-6361/201526616. hal-03366927

HAL Id: hal-03366927

https://hal.science/hal-03366927v1

Submitted on 13 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Gaps, rings, and non-axisymmetric structures in protoplanetary disks: Emission from large grains

J. P. Ruge<fnref n="1" />, M. Flock<fnref n="2" />,<fnref n="3" />, S. Wolf<fnref n="1" />, N. Dzyurkevich<fnref n="4" />, S. Fromang<fnref n="3" />, Th. Henning<fnref n="5" />, H. Klahr<fnref n="5" />, and H. Meheut<fnref n="3" />,<fnref n="6" />

Received 27 May 2015 / Accepted 15 March 2016

## ABSTRACT

**Aims.** Dust grains with sizes around (sub)mm are expected to couple only weakly to the gas motion in regions beyond 10 au of circumstellar disks. In this work, we investigate the influence of the spatial distribution of these grains on the (sub)mm appearance of magnetized protoplanetary disks.

**Methods.** We perform non-ideal global 3D magneto-hydrodynamic (MHD) stratified disk simulations, including particles of different sizes (50 µm to 1 cm), using a Lagrangian particle solver. Subsequently, we calculate the spatial dust temperature distribution, including the dynamically coupled submicron-sized dust grains, and derive ideal continuum re-emission maps of the disk through radiative transfer simulations. Finally, we investigate the feasibility of observing specific structures in the thermal re-emission maps with the Atacama Large Millimeter/submillimeter Array (ALMA).

**Results.** Depending on the level of turbulence, the radial pressure gradient of the gas, and the grain size, particles settle to the midplane and/or drift radially inward. The pressure bump close to the outer edge of the dead-zone leads to particle-trapping in ring structures. More specifically, vortices in the disk concentrate the dust and create an inhomogeneous distribution of solid material in the azimuthal direction. The large-scale disk perturbations are preserved in the (sub)mm re-emission maps. The observable structures are very similar to those expected from planet-disk interaction. Additionally, the larger dust particles increase the brightness contrast between the gap and ring structures. We find that rings, gaps, and the dust accumulation in the vortex could be traced with ALMA down to a scale of a few astronomical units in circumstellar disks located in nearby star-forming regions. Finally, we present a brief comparison of these structures with those recently found with ALMA in the young circumstellar disks of HL Tau and Oph IRS 48.

**Key words.** submillimeter: planetary systems – magnetohydrodynamics (MHD) – turbulence – planets and satellites: formation – accretion, accretion disks

## 1. Introduction

Recent observations of protoplanetary disks have revealed the existence of a wealth of structures with various properties, from dust asymmetries and concentrations in the outer regions of protoplanetary disks (e.g., Oph IRS 48: van der Marel et al. 2013; LkHα 330: Isella et al. 2013; SAO 206462 and SR 21: Pérez et al. 2014; or in HD 142527: Casassus et al. 2013; Fukagawa et al. 2013) to the amazing image of HL Tau (ESO press release 1436 and Partnership et al. 2015) obtained during the ALMA science verification phase, using its most extended configuration and revealing without ambiguity a multiple dust ring structure. These exciting observations have often been interpreted as the signature of the presence of a young planet in the disk. This is because the spatial distribution of dust is believed to trace the underlying gaseous structure of the disk.

Indeed, while small dust grains of submicron size are well coupled to the motion of the gas in circumstellar disks, this is

not the same for particles of (sub)mm size, which form from collisional growth owing to various underlying physical processes (e.g., Testi et al. 2014). For these larger particles, the coupling between the grains and the gas is reduced (Beckwith et al. 2000) and large-scale gas structures, such as radial pressure gradients or vortices, are able to locally concentrate these partially decoupled particles (Barge & Sommeria 1995; Tanga et al. 1996; Klahr & Henning 1997; Chavanis 2000; Johansen et al. 2004; Lyra et al. 2009a; Regály et al. 2012; Meheut et al. 2012; Pinilla et al. 2012; Birnstiel et al. 2013; Lyra & Lin 2013; Zhu & Stone 2014). Various mechanisms that create these large-scale structures have been proposed, such as variations in the accretion stress at the gas ionization boundaries (Varnière & Tagger 2006; Lyra et al. 2009b, 2015; Dzyurkevich et al. 2010; Faure et al. 2015; Flock et al. 2015), sublimation zones (Kretke & Lin 2007; Brauer et al. 2008; Dzyurkevich et al. 2013; Zhang et al. 2015), sintering zones (Okuzumi et al. 2016), planet-disk interactions (Paardekooper & Mellema 2004; Zhu et al. 2014), photo

<footnote n="1">¹ Universität zu Kiel, Institut für Theoretische Physik und Astrophysik, Leibnitzstr. 15, 24098 Kiel, Germany
e-mail: ruge@astrophysik.uni-kiel.de</footnote>

<footnote n="2">² Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109, USA
e-mail: mflock@caltech.edu</footnote>

<footnote n="3">³ CEA UMR AIM Irfu, SAP, CEA-CNRS-Univ. Paris Diderot, Centre de Saclay, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="4">⁴ Laboratoire de radioastronomie, UMR 8112 du CNRS, École normale supérieure et Observatoire de Paris, 24 rue Lhomond,
75231 Paris Cedex 05, France</footnote>

<footnote n="5">⁵ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="6">⁶ Laboratoire Lagrange, Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Bd de l'Observatoire, CS 34229,
06304 Nice Cedex 4, France</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Radial profile of the initial particles' Stokes number at the mid-plane. The green and the blue line mark the regime of the considered particle sizes.</figcaption></figure>

evaporation (Espaillat et al. 2012) or gravitational disk instabilities (Bertin & Lodato 2001; Vorobyov & Basu 2005). Further investigation confirmed their observability, e.g., for planet-disk interactions (Wolf & D'Angelo 2005; Gonzalez et al. 2012; Ruge et al. 2013), self gravitating disks (Cossins et al. 2010; Dipierro et al. 2014; Pohl et al. 2015; Takahashi et al. 2016) or vortices (Wolf & Klahr 2002; Pinilla et al. 2012). As we can see, the presence of a planet in the disk as the origin of the recently observed structure is only one possibility among others and must be examined critically. In addition, these disk systems are evolving on timescales of several millions of years with ongoing dust growth and reprocessing (Williams & Cieza 2011; Wolf et al. 2012) but the exact age of the observed system is difficult to estimate.

In this context, we have shown in Flock et al. (2015) that it is possible to explain the pronounced dust rings and local concentrations in magnetized disks without a planet. At the outer edge of the dead zone, the response of the magneto-rotational instability (MRI) to variations of the gas density supports the formation of a large-scale gap and jump structure in the gas density distribution. The surface density maximum triggers the Rossby wave instability (Lovelace et al. 1999; Lesur & Papaloizou 2009; Lin & Papaloizou 2010; Lin 2012a,b), which is a special form of Kelvin Helmholtz instability (Papaloizou & Pringle 1984), resulting in the formation of a vortex. In the present study, we further investigate the spatial and temporal distribution of large dust particles (radii ranging from 50 µm to 1 cm). We focus on the impact of these dust grains on the (sub)mm appearance of the circumstellar disk. For this purpose, we calculate the motion of each particle individually using a Lagrangian method and simulate the thermal re-emission maps based on the resulting dust density distribution with the continuum radiative transfer code MC3D (Wolf et al. 1999; Wolf 2003). The results are compared with our previous model, which include only small dust grains (Flock et al. 2015). Therefore, we make use of the same disk setup, which is characterized by a total disk mass of 0.042 $M_\odot$, a disk radius ranging from 20 au to 100 au, a dead zone from 20 au to 40 au, a central star with an effective temperature of 4000 K, a luminosity of 0.95 $L_\odot$, and a stellar mass of 0.5 $M_\odot$.

## 2. Modeling techniques

We use non-ideal global 3D MHD stratified simulations of a circumstellar disk model as described in our previous work (see model D2G_e-2 in Sect. 2 and Fig. 7 in Flock et al. 2015). The dust particles with perfect coupling to the gas motion have sizes of $a_{min} = 0.005$ µm ≤ $a$ ≤ 0.25 µm = $a_{max}$ and are called small dust particles in the present work. The grains are homogeneous, compact and spherical, and follow the size distribution

$$ \frac{dn(a)}{da} \propto a^{-3.5}. \qquad (1) $$

In addition to the distribution of the small dust, we consider spherical particles with the following radii: 50 µm, 90 µm, 0.2 mm, 0.3 mm, 0.5 mm, 0.9 mm, 1.7 mm, 3.1 mm, 5.5 mm, and 10 mm. We refer to these as large dust particles. We consider the external forces of gas drag and gravity. To calculate the grain trajectories, we solve the equation of motion in the spherical coordinates system ($r, \vartheta, \phi$):

$$ \frac{\partial v_r^{\text{par}}}{\partial t} = \frac{(v_r^{\text{gas}} - v_r^{\text{par}})}{\tau_t} + \frac{l_{\vartheta}^2}{r^3} + \frac{l_{\phi}^2}{r^3(\sin\vartheta)^2} - \frac{GM_{\star}}{r^2} \qquad (2) $$

$$ \frac{\partial l_{\vartheta}}{\partial t} = \frac{(v_{\vartheta}^{\text{gas}} \cdot r - l_{\vartheta})}{\tau_t} + \frac{l_{\phi}^2 \cdot \cos\vartheta}{r^2(\sin\vartheta)^3} \qquad (3) $$

$$ \frac{\partial l_{\phi}}{\partial t} = \frac{v_{\phi}^{\text{gas}} \cdot r \cdot \sin\vartheta - l_{\phi}}{\tau_t}, \qquad (4) $$

with the coupling time $\tau_t$, the gravitational constant G, the stellar mass $M_\star$ and the angular moments of the particle $l_\vartheta = v_\vartheta^{\text{par}} \cdot r$ and $l_\phi = v_\phi^{\text{par}} \cdot r \cdot \sin\vartheta$. The particle time integration is performed using a leapfrog integrator. We use outflow boundary conditions for the gas and the particles. In this work, we focus on the Epstein regime, which is valid when the particle size $a$ is smaller than 9/4 of the mean free path of the gas molecules (Weidenschilling 1977). In this case, the coupling time $\tau_t$ can be expressed as:

$$ \tau_t = \frac{\rho_d a}{\rho_g c_s}, \qquad (5) $$

with the solid density $\rho_d = 2.7$ g cm$^{-3}$, the gas density $\rho_g$, and the sound speed $c_s$. The Stokes number $St = \tau_t \cdot \Omega$ of the large dust particles is shown in Fig. 1. Initially, we place 50 000 particles in each size bin, randomly distributed in the simulation domain. The initial velocity of the particles is set to match the local Keplerian velocity. In Appendix B.3 we show that the effect of the number of particles and a different distribution of the particles onto the results remains insignificant.

Two snapshots of the simulation are considered in the following: one after 150 inner orbits of the disk, corresponding to 18 900 yr of disk evolution, and one after 600 inner orbits, corresponding to 75 600 yr. After 150 orbits of the inner disk, the (sub)mm dust particles have accumulated into ring structures (see Fig. 2, left) owing to a pressure jump at the outer dead-zone edge (see Fig. 3). In Fig. 3, we project the distribution of the larger particles in the 2D R-Z/R plane. We note that the larger particles are located along an azimuthal ring, while they slowly oscillate in the vertical direction. Another time snapshot could show a different distribution for $z > 0$ and $z < 0$. During the following hundreds of orbits, the Rossby wave instability generates a vortex that concentrates dust particles, depending on their size. We plot the spatial distribution of the
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Scatter plot of every tenth particle in four selected size bins of the disk in face-on orientation. **Left:** time step after 150 inner orbits (corresponding to 18 900 yr of evolution). **Right:** time step after 600 inner orbits (75 600 yr). The particle size is given in units of mm. After 150 inner orbits, the particles are concentrated into rings. After 600 inner orbits, the vortex has broken the ring structures and concentrated the larger particles.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Scatter plot of every tenth particle in four selected size bins of the disk in edge-on orientation for the time step after 150 inner orbits (top) and 600 inner orbits (bottom). The dead-zone outer edge is indicated in the plot. Particles are color-coded according to their size, as in Fig. 2. The black solid line shows the surface density profile of the disk after 150 inner orbits (18 900 yr), the dashed line shows this profile for the disk after 600 inner orbits (75 600 yr).</figcaption></figure>

grains that results at 600 orbits on the right panel of Fig. 2. In agreement with the findings of Meheut et al. (2012), we observe an increased concentration of larger particles towards the vortex center. We note that we can only observe the increase of the dust-to-gas mass ratio and not its saturated value

since we do not take the dust back-reaction onto the gas into consideration. Moreover, we do not resupply particles at the outer domain boundary, and the presented accumulation of solid material in the computational domain should thus be seen as a lower limit. During the full simulation run, we observe an increase in the dust-to-gas mass ratio by 20% when averaging over a disk scale height. However, we note that, on scales smaller than tenths of a disk scale height, the increase of the dust component can be larger. Only future very high-resolution dust and gas simulation, including the back-reaction of the dust on the gas and a consistent resupply of particles can investigate the effect of the solid component on the gas.

Inside the dead zone (for $r < 40$ au), all large particles remain in the midplane (see Fig. 3). Outside of the dead zone, the 50 µm particles are still efficiently mixed in the vertical direction up to about one gas scale height. In Fig. 3, we overplot the surface density of the gas at the two selected time steps. The gas surface density of the gap and jump decreases between the two time steps outward of 40 au (see Fig. 3). This decrease in gas density leads to reduced coupling between the particles to the gas, thereby favoring the maintenance of the ring over long timescales. The dynamical process is conducive to the formation of a dust ring by the evolution of the gas density and the subsequent decoupling of the dust have recently been demonstrated in a study by Meru et al. (2014).

Figure 4 shows, in addition to the gas surface density, the total dust surface density and the surface density for the reference dust components shown in Figs. 2 and 3. The plot shows that the dust surface density is mainly determined by the larger particles, which show efficient concentration at the pressure maximum. We also note that the density gap at 20 au is clearly induced by the outflow boundary since the surface density is truncated in a region between 20–25 au, which causes a local pressure maximum and thus the evolution of a vortex there. The actual inner dead-zone edge should be much closer in.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Surface density plot for the gas (thick black solid line), the total dust (dotted line) and the 4 representative dust components, 1.71 mm (blue), 527 µm (green), 162 µm (red), and 50 µm (black), for the time step after 150 inner orbits (top) and 600 inner orbits (bottom).</figcaption></figure>

**Radiative transfer (RT):** To calculate the temperature distribution and thermal re-emission maps of the disks, we use the RT code MC3D (Wolf et al. 1999; Wolf 2003). We focus on the continuum emission of the dust phase of the disk, assuming that the dust material is composed of 62.5% silicate and 37.5% graphite (optical data by Weingartner & Draine 2001). We calculate the wavelength-dependent absorption efficiencies for the different grain sizes with Mie theory, using the tool MieX (Wolf & Voshchinnikov 2004).

In this study, the dust phase consists of two components. The first component is the system of small dust particles that fully couple to the motion of the gas in the disk (see Sect. 2). The distribution of these particles follows the gas distribution. The second component of the dust phase is represented by large dust particles (see Sect. 2). As for the small particles, the large particles initially follow the continuous power-law size distribution given by Eq. (1). Each grain of size $a_g$ in the simulation represents a number of dust particles $\bar{N}_g(a_g)$. Appendix A presents a detailed description of the calculation of $\bar{N}_g(a_g)$.

### 3. Results and discussion

Our goal is to explore the impact of the large dust particles on the thermal re-emission of circumstellar disks. Our analysis is based on the simulated 441 µm, 871 µm, 1.3 mm, and 2.0 mm re-emission maps of our disks (assuming they are located at a distance of 140 pc from the observer) at the two selected time steps. In particular, we compare these maps to our previous models with small dust grains (Flock et al. 2015).

The additional large dust grains increase the total flux of the thermal disk re-emission at 1.3 mm by a factor of ~1.5. This is due to the absorption efficiency of the large grains that is several orders of magnitude higher in comparison to that of the small

grains at wavelengths >100 µm<fnref n="1" />. In scattered light at a wavelength of 2.2 µm (K band), the appearance of the disk is not influenced by the large dust particles at all, because they are all confined to less than one scale height (see Sect. 2 and Fig. 3), and this part of the disk is not traced by scattered light at this wavelength (Ruge et al. 2013).

The large dust grains are not homogeneously distributed in the disk. Their distribution is also time dependent (see Sect. 2 and Fig. 2). As a result, they are able to change the (sub)mm appearance of a circumstellar disk significantly within a few 10⁴ yr (one inner orbit ~126 yr), even though they only contribute a fraction ~0.07% of the total mass of the disk (compare left and right column of Fig. 5). This is in contrast with the results in Sect. 4.2 of Flock et al. (2015), where the re-emission structure of the disk – with only small grains – remains constant over the same period

The thermal emission maps in Fig. 5 show three dominating structures, which are due to the distribution of the (large) dust in the disk (see Fig. 2). At both time steps, the gap, which has already been discussed by Flock et al. (2015), is present in the re-emission maps. After 150 inner orbits (18 900 yr) the large dust particles have accumulated into rings. In the thermal emission maps, this leads to an enlarged brightness-contrast ratio between the innermost region of the gap and its edges. In comparison to the results obtained without the larger dust grains, the contrast is ~1.5 times higher (see Fig. 6). The gap appears deeper in this case. Similar results have been found for the local concentrations of larger, weakly coupled dust grains in the case of planet-disk interaction by Paardekooper & Mellema (2004), Rice et al. (2006), Fouchet et al. (2007, 2010), Lyra et al. (2009c), Ayliffe et al. (2012), Zhu et al. (2012).

After 600 inner orbits (75 600 yr) the disk shows a large vortex, which becomes visible owing to the concentration of larger particles (see Fig. 2). With the exception of the 50 µm particles, the vortex has destroyed the dust ring structures and concentrated the larger dust particles azimuthally. The size-dependent dynamics of the dust grains lead to different thermal emission maps at different wavelengths (see Fig. 5, right). Observations at longer wavelengths (≥1.3 mm) trace the dust concentration in the vortex very well, while the radiation of the smallest particles in our setup dominates at shorter wavelengths and is more smoothly distributed (see Fig. 5, right). A similar trend was seen in a recent observations of HL Tau at 7 mm (Carrasco-Gonzalez et al. 2016), which shows clear azimuthal variations inside the ring at 7 mm, compared to the smooth emission at 1.3 mm. Eventually, this offers the possibility of determining the size of the particles in the vortex by using spatially resolved maps of the spectral index (see Sect. 4).

As already mentioned in Sect. 2, the set of large dust particles samples a continuous grain-size distribution. Consequently, the distribution of dust grains of a certain size, that is not included in our setup, can be approximated by interpolating from the upper and lower size bins. We have checked that a higher sampling of the dust size bins has no effect onto the disk structures.

We note again that the vortex lifetime is around 40 local orbits until it is destroyed and reformed again on a timescale of 20 local orbits (Flock et al. 2015). We do not observe any migration of the vortex. We point out that we consider Ohmic resistivity as the sole source of magnetic field dissipation and we

<footnote n="1">¹ We note that the factor of 1.5 should be understood as a lower limit since we do not resupply the grains in the outer domain. Especially the particles (radii ≥3.1 mm) drift quickly radially inward and leave the simulation.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Selected re-emission maps for the disk after 150 inner orbits (18 900 yr, left column) and 600 (75 600 yr, right column). Each wavelength (441 µm, 871 µm, 1.3 mm and 2.0 mm) is presented in its own line. The square root of thermal re-emission flux per arcsec² is shown color-coded. At both time steps and for every wavelength a gap is visible in the maps. Depending on time and wavelength, the re-emission maps of the disk show different structures (rings and concentrations).</figcaption></figure>

neglect, for simplicity, the Hall effect and ambipolar diffusion. Recent results have shown their importance for MRI-driven turbulence in circumstellar disks (Bai & Stone 2013; Lesur et al. 2014; Béthune et al. 2016). They could potentially shift the MRI inactive zone radially outward or even trigger new instabilities. In addition, in this region the electron heating could become significant, which was recently shown by Mori & Okuzumi (2016).

## 4. Comparison to recent ALMA observations

Dust ring structures: the ring structures found in the density distribution of the disk after 150 inner orbits (18 900 yr) appear qualitatively similar to the structures recently observed in the disk of HL Tau (Partnership et al. 2015). Therefore, we explore how an equivalent (simulated) ALMA observation of our disk model would look like. We select the position of

<figure><img src="image_7.png" /><figcaption>Fig. 6. Comparison of radial brightness profiles of the disk after 150 inner orbits (18 900 yr) with (solid blue line) and without additional large dust particles (dashed red line) at a wavelength of 1.3 mm. The contrast ratio between the brightness minimum, which corresponds to the gap center, and the brightness maximum of the disk for radii >55 au is enlarged by a factor of ~1.5 through the additional large dust particles, in comparison to the case with only small grains. The total number of grains particles $N_{max}$ is the same for both the blue and the red curve (see Appendix A), which means that the number of small grains used to compute the red curve is reduced in comparison to the blue-lined case. This leads to a decreased flux at the gap center.</figcaption></figure>

HL Tau (RA: 4 h31min38.4s, Dec: +18°13′57.4″ J2000, distance: 140 pc; Kwon et al. 2011) and its inclination of 40° (Kwon et al. 2011). The thermal emission maps are re-calculated under the assumption of the stellar luminosity ($L_{\star} \sim 8.3 L_{\odot}$) and an effective temperature (4000 K) of HL Tau (Kwon et al. 2011). We focus on the wavelengths of 871 µm and 1.3 mm, and select very extended configurations of the ALMA array with maximum baselines of ~8 km and ~15 km. The observing time is chosen to be 6 h. The resulting spatial resolution is 49 mas × 39 mas at 871 µm and 37 mas × 32 mas at 1.3 mm. Thermal and phase noise are added (see the details on our ALMA setup in Flock et al. 2015). We emphasize that the goal of these simulations is to study qualitatively the extent to which the ring structure found in our simulations resembles the structure found in the HL Tau disk, while a quantitative fit to the observational results of HL Tau, which would require a detailed consideration of the envelope that has been investigated in RT models by Men'shchikov et al. (1999).

ALMA enables the detection of single dust rings with this configuration for an observing wavelength of 1.3 mm (see Fig. 7, top right). At a wavelength of 871 µm, the dust ring structure is less pronounced (see Fig. 7, top left). This is due to the emission of the smaller grains that are distributed homogeneously in the disk and still efficient at this wavelength (Fig. 5). These differences offer a possibility to identify dust concentrations as the origin of the ring structures. From the wavelength-dependence of the emission, one can constrain the size of the re-emitting particles. For this purpose, we calculate the spatially resolved map of the spectral index from both simulated ALMA observations (see Fig. 7, bottom left) and the profile of the spectral index in a selected cut through the map (see Fig. 7, bottom right). The spectral index is defined as

$$ \alpha_{Sp} = - \frac{\partial \log_{10} (F_{\lambda})}{\partial \log_{10} (\lambda)} \qquad (6) $$

At the position of the dust rings, the spectral index is reduced, which indicates the presence of larger particles. Within the gap in the disk, the spectral index is increased, which indicates that only smaller dust grains are present in this region. We note that the asymmetric feature, seen in the spectral index map, could
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** *Top:* simulated ALMA images of the disk after 150 inner orbits (18 900 yr) at a wavelength of 871 µm (left) and 1.3 µm (right). The position, the disk inclination, distance, and stellar luminosity of HL Tau were used for these simulation. At both wavelengths, a gap is visible. The gap in the images is emphasized and we note that its origin lies in the magnetic field in the disk, and not a planet. *Bottom left:* spatially resolved map of the spectral index calculated from the simulated observation on the top. The map is rotated clockwise by 40° in comparison to the maps at the top. *Bottom right:* profile of the spectral index along the cut, indicated by a dashed line on the bottom left figure.</figcaption></figure>

arise from the projection as the disk is inclined along that di-
rection. Finally, we conclude that similar (sub)mm appearances
of circumstellar disks can result from an inhomogeneous grain-
size distribution in the radial direction, as discussed above, with-
out the necessity to invoke the presence of a planet, as is usu-
ally done. We emphasize that the position of the gap and ring
structure is located at the dead-zone outer edge, which depends
mainly on the surface density amount and profile of gas and dust
(Dzyurkevich et al. 2013; Flock et al. 2015). For the parameters
used here, we expect this to happen at any location where the to-
tal surface density becomes around 15 g cm⁻², assuming a dust-
to-gas mass ratio of 0.01, since this value marks the threshold
for MRI activity. We note, again, that the ring state, featured in
Fig. 7, is present more often than just once in the simulation. As
we have pointed out, the lifetime of the vortex is around 40 local
orbits. After this time, the particles will again be only radially
concentrated. It is true that each individual ring state will look
different, especially if we do not resupply the particles inside
the domain. The vortex state is still more likely to be present
(2/3 vortex state vs. 1/3 ring state).

**Dust concentrations in vortices:** Due to the concentration of large dust particles in the vortex, the re-emission maps of the disk after 600 inner orbits (75 600 yr) appear similar to various recent ALMA observations of young circumstellar disks (see Sect. 1). To investigate this similarity, we explore what a simulated ALMA observation of this disk model would look like. We select the position, distance, and object inclination of Oph IRS 48 (RA: 16 h27min37.8s, Dec: -24°30′35.3″ J2000, distance: 120 pc, inclination: 50°; van der Marel et al. 2013). Other model parameters, such as stellar temperature and luminosity, or disk mass are not adapted to the parameters of Oph IRS 48. Our goal is to show the feasibility of observing the vortex in the disk. In comparison to the thermal emission maps in the right column of Fig. 5, we rotate the model counterclockwise by 90°. We select a wavelength of 1.3 mm, an observing time of 5 h, and a configuration of the ALMA array with a maximum baseline of 3.1 km, corresponding to a spatial resolution of 0.13″ × 0.11″.

The concentration of large dust particles in the vortex in-
creases the flux at the vortex position by a factor ~2.8 in the
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Simulated ALMA images of the disk after 600 inner orbits (75 600 yr) at a wavelength of 1.3 mm. The object's position, the disk inclination, and distance of Oph IRS 48 have been used. In comparison to Fig. 5, the model was rotated by 90° counterclockwise. We note that the inner edge of the disk represents the beginning of the computational domain.</figcaption></figure>

simulated ALMA observation, compared to the models without large dust grains. The peak flux at the vortex, compared to the background, is around 10, which is lower then the value 130 found by (Zhu & Stone 2014). However, we note that several factors influence this peak emission, especially the dust concentration and size of the vortex. For example, the planet-induced vortex by (Zhu & Stone 2014) appears to be larger in size. Finally, we show the feasibility of detecting a vortex through its ability to concentrate large dust grains in a circumstellar disk again, without the necessity of a planetary mass perturber.

## 5. Conclusion

In this work, we have followed the motion and concentration of large dust particles with different sizes in a typical protoplanetary disk using a non-ideal 3D MHD disk simulation. We have performed RT calculations to derive the temperature profile and thermal emission maps of the disk. Our goal was to investigate the impact of (sub)mm-sized dust grains on spatially resolved thermal emission maps of circumstellar disks and the feasibility of detecting selected structures with ALMA. In summary, we find:

- The larger dust particles (with radii ranging from 50 µm to 1 cm) are weakly coupled to the gas motion and accumulate in rings.

- Vortices in the disk can break up those rings and concentrate the dust particles azimuthally. This leads to an increased local dust-to-gas ratio in the vortex center.

- Owing to the high efficiency of emission of these large dust grains, the dust rings, as well as the dust concentration in vortices, are emphasized in the thermal emission maps of the disk. For example, the higher concentration of the larger dust grains enlarges the brightness contrast between the minimum in the gap and the peak outside by at least a factor of 1.5 at a wavelength of 1.3 mm, compared to the results of Flock et al. (2015).

- Depending on the wavelength, the gap, dust rings, and dust concentration in vortices can be observed with ALMA.

- The spatial distribution of dust grains of different size can also be identified by a decreased spectral index in spatially resolved spectral index maps calculated on the basis of a simulated ALMA observation.

- The ring structures, as well as the dust concentration in the vortex, are stable for at least 10⁴ yr. In addition, the disk alternates between a state with and without azimuthal dust concentration by a vortex, as already shown by Flock et al. (2015), which significantly increases the likelihood of observing the presented structures. The vortex state is more likely to observe (by 66%) compared to the ring state.

We conclude that inhomogeneous and non-axisymmetric distributions of larger dust grains have a major impact on the (sub)mm appearance of a circumstellar disk. Our work shows that these inhomogeneities can arise naturally as a result of internal dynamical processes associated with the presence of a weak magnetic field, without having to rely on a planet orbiting in the disk. More constraints will be required in the future to disentangle the origin of the observed (sub)mm structures in protoplanetary disks, including observations of the magnetic field (e.g., HL Tau Stephens et al. 2014) and/or of the disk kinematics (Simon et al. 2015).

**Acknowledgements.** We thank the anonymous referee for providing extensive and constructive comments, which significantly improved the paper. We acknowledge financial support from the German Research Foundation (J.P. Ruge: WO 857/10-1; H.H. Klahr: KL 14699-1). Mario Flock, Sebastien Fromang, and Heloise Meheut are supported by the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) / ERC Grant agreement nr. 258729. Parallel computations were performed on the IRFU COAST cluster located at CEA IRFU. This research was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration and with the support of the NASA Exoplanet Research program via grant 14XRP14_20153. Copyright 2016 California Institute of Technology. Government sponsorship acknowledged.

## References

Ayliffe, B. A., Laibe, G., Price, D. J., & Bate, M. R. 2012, MNRAS, 423, 1450

Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76

Barge, P., & Sommeria, J. 1995, A&A, 295, L1

Beckwith, S. V. W., Henning, T., & Nakagawa, Y. 2000, Protostars and Planets IV (Tucson: University of Arizona Press), 533

Bertin, G., & Lodato, G. 2001, A&A, 370, 342

Béthune, W., Lesur, G., & Ferreira, J. 2016, A&A, 589, A87

Birnstiel, T., Dullemond, C. P., & Pinilla, P. 2013, A&A, 550, L8

Brauer, F., Henning, T., & Dullemond, C. P. 2008, A&A, 487, L1

Carrasco-Gonzalez, C., Henning, T., Chandler, C. J., et al. 2016, ApJ, 821, L16

Casassus, S., van der Plas, G., M, S. P., et al. 2013, Nature, 493, 191

Chavanis, P. H. 2000, A&A, 356, 1089

Cossins, P., Lodato, G., & Testi, L. 2010, MNRAS, 407, 181

Dipierro, G., Lodato, G., Testi, L., & de Gregorio Monsalvo, I. 2014, MNRAS, 444, 1919

Dzyurkevich, N., Flock, M., Turner, N. J., Klahr, H., & Henning, T. 2010, A&A, 515, A70

Dzyurkevich, N., Turner, N. J., Henning, T., & Kley, W. 2013, ApJ, 765, 114

Espaillat, C., Ingleby, L., Hernández, J., et al. 2012, ApJ, 747, 103

Faure, J., Fromang, S., Latter, H., & Meheut, H. 2015, A&A, 573, A132

Flock, M., Ruge, J. P., Dzyurkevich, N., et al. 2015, A&A, 574, A68

Fouchet, L., Maddison, S. T., Gonzalez, J.-F., & Murray, J. R. 2007, A&A, 474, 1037

Fouchet, L., Gonzalez, J.-F., & Maddison, S. T. 2010, A&A, 518, A16

Fukagawa, M., Tsukagoshi, T., Momose, M., et al. 2013, PASJ, 65, L14

Gonzalez, J.-F., Pinte, C., Maddison, S. T., Ménard, F., & Fouchet, L. 2012, A&A, 547, A58

Isella, A., Pérez, L. M., Carpenter, J. M., et al. 2013, ApJ, 775, 30

Johansen, A., Andersen, A. C., & Brandenburg, A. 2004, A&A, 417, 361
Klahr, H. H., & Henning, T. 1997, Icarus, 128, 213

Kretke, K. A., & Lin, D. N. C. 2007, ApJ, 664, L55

Kwon, W., Looney, L. W., & Mundy, L. G. 2011, ApJ, 741, 3

Lesur, G., & Papaloizou, J. C. B. 2009, A&A, 498, 1

Lesur, G., Kunz, M. W., & Fromang, S. 2014, A&A, 566, A56

Lin, M.-K. 2012a, ApJ, 754, 21

Lin, M.-K. 2012b, MNRAS, 426, 3211

Lin, M.-K., & Papaloizou, J. C. B. 2010, MNRAS, 405, 1473

Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805

Lyra, W., & Lin, M.-K. 2013, ApJ, 775, 17

Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2009a, A&A, 493, 1125

Lyra, W., Johansen, A., Zsom, A., Klahr, H., & Piskunov, N. 2009b, A&A, 497, 869

Lyra, W., Johansen, A., Zsom, A., Klahr, H., & Piskunov, N. 2009c, A&A, 497, 869

Lyra, W., Turner, N. J., & McNally, C. P. 2015, A&A, 574, A10

Meheut, H., Meliani, Z., Varniere, P., & Benz, W. 2012, A&A, 545, A134

Men'schikov, A. B., Henning, T., & Fischer, O. 1999, ApJ, 519, 257

Meru, F., Quanz, S. P., Reggiani, M., Baruteau, C., & Pineda, J. E. 2014, ApJ, submitted [arXiv:1411.5366]

Mori, S., & Okuzumi, S. 2016, ApJ, 817, 52

Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610

Okuzumi, S., Momose, M., Sirono, S.-i., Kobayashi, H., & Tanaka, H. 2016, ApJ, 821, 82

Paardekooper, S.-J., & Mellema, G. 2004, A&A, 425, L9

Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721

Partnership, A., Brogan, C. L., Perez, L. M., et al. 2015, ApJ, 808, L3

Pérez, L. M., Isella, A., Carpenter, J. M., & Chandler, C. J. 2014, ApJ, 783, L13

Pinilla, P., Birnstiel, T., Ricci, L., et al. 2012, A&A, 538, A114

Pohl, A., Pinilla, P., Benisty, M., et al. 2015, MNRAS, 453, 1768

Regály, Z., Juhász, A., Sándor, Z., & Dullemond, C. P. 2012, MNRAS, 419, 1701

Rice, W. K. M., Armitage, P. J., Wood, K., & Lodato, G. 2006, MNRAS, 373, 1619

Ruge, J. P., Wolf, S., Uribe, A. L., & Klahr, H. H. 2013, A&A, 549, A97

Simon, J. B., Hughes, A. M., Flaherty, K. M., Bai, X.-N., & Armitage, P. J. 2015, ApJ, 808, 180

Stephens, I. W., Looney, L. W., Kwon, W., et al. 2014, Nature, 514, 597

Takahashi, S. Z., Tsukamoto, Y., & Inutsuka, S.-i. 2016, MNRAS, 458, 3597

Tanga, P., Babiano, A., Dubrulle, B., & Provenzale, A. 1996, Icarus, 121, 158

Testi, L., Birnstiel, T., Ricci, L., et al. 2014, Protostars and Planets VI (Tucson: University of Arizona Press), 339

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199

Varnière, P., & Tagger, M. 2006, A&A, 446, L13

Vorobyov, E. I., & Basu, S. 2005, ApJ, 633, L137

Weidenschilling, S. J. 1977, MNRAS, 180, 57

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67

Wolf, S. 2003, Comput. Phys. Comm., 150, 99

Wolf, S., & D'Angelo, G. 2005, ApJ, 619, 1114

Wolf, S., & Klahr, H. 2002, ApJ, 578, L79

Wolf, S., & Voshchinnikov, N. V. 2004, Comput. Phys. Comm., 162, 113

Wolf, S., Henning, T., & Stecklum, B. 1999, A&A, 349, 839

Wolf, S., Malbet, F., Alexander, R., et al. 2012, A&ARv, 20, 52

Zhang, K., Blake, G. A., & Bergin, E. A. 2015, ApJ, 806, L7

Zhu, Z., & Stone, J. M. 2014, ApJ, 795, 53

Zhu, Z., Stone, J. M., & Rafikov, R. R. 2012, ApJ, 758, L42

Zhu, Z., Stone, J. M., Rafikov, R. R., & Bai, X.-N. 2014, ApJ, 785, 122
## Appendix A: Calculating $\tilde{N}_g$

In the case of large dust grains, each of the initial “test” particles of a certain radius $a_g$ represents a certain number of dust particles $\tilde{N}_g(a_g)$. We assume that the initial dust mass $M_d$ amounts to 1% of the total gas mass $M_{gas}$ and is the sum of the mass of the small dust $M_{sd}$ and the mass of the larger dust $M_{ld}$. In the radiative transfer calculation, the dust properties of small dust are discretized into 1000 bin sizes, therefore the mass $M_{sd}$ can be expressed in the following way:

$$M_{sd} = \frac{4}{3} \pi \rho_{md} \sum_{i=1}^{1000} a_i^3 \cdot N(a_i), \qquad (A.1)$$

where $\rho_{md} = 2.7 \text{ g cm}^{-3}$ is the density of the dust material and $N(a_i)$ is the absolute number of particles with radius $a_i$. Following the grain size distribution (see Eq. (1)), $N(a_i)$ can be written as a fraction of a reference number of particles. Without loss of generality, we select the total amount of particles $N_{max}$ with radius $a_{max}$:

$$N(a_i) = N_{max} \cdot \left( \frac{a_i}{a_{max}} \right)^{-2.5}. \qquad (A.2)$$

The mass of the large dust grains can be described in the same way:

$$M_{ld} = \frac{4}{3} \pi \rho_{md} \sum_{g=1}^{10} a_g^3 \cdot N(a_g). \qquad (A.3)$$

Because of Eq. (A.2), we finally get

$$0.01 M_{gas} = M_d = M_{sd} + M_{ld} = \frac{4}{3} \pi \rho_{md} \frac{N_{max}}{a_{max}} \left( \sum_{i=1}^{1000} a_i^{0.5} + \sum_{g=1}^{10} a_g^{0.5} \right). \qquad (A.4)$$

Through reordering of Eq. (A.4), $N_{max}$ – and therefore every $N(a_i)$ and $N(a_g)$ – can be calculated. Finally, $\tilde{N}_g$ is given by

$$\tilde{N}_g = \frac{N(a_g)}{50000}. \qquad (A.5)$$

To calculate the dust surface density we simply integrated over a radial bin size $\Delta r$ and then calculate the given dust mass per unit area:

$$\Sigma_{\text{dust}} = \frac{\int_{\Delta r} M_d dr}{\int_{\Delta r} \int_{\phi} r dr d\phi}. \qquad (A.6)$$

## Appendix B: Particle method validation

In this section, we performed global hydrodynamical simulations to validate our method. For the test models, we use the setup by Nelson et al. (2013) which describes the hydrostatic disk solution that is dependent on the two parameters $p$ and $q$.

$$\rho = \rho_0 \left(\frac{R}{R_0}\right)^p \cdot e^{\frac{GM}{c_s^2} \left[\frac{1}{r} - \frac{1}{R}\right]} \qquad (B.1)$$

$$v_{\phi} = R\Omega_{K} \cdot \left( 1.0 + (p+q) \cdot \frac{H^2}{R} + q \cdot (1.0 - \sin\theta) \right)^{0.5}, \qquad (B.2)$$

with the cylindrical radius $R = r \sin(\theta)$, the spherical radius $r$, the orbital frequency $\Omega_{K}$, the sound speed $c_s = H \cdot \Omega_{K}$, and the

<figure><img src="image_10.png" /><figcaption>Fig. B.1. Vertical motion of a particle with Stokes number of $St = 300$ (green solid line), a model of a damped oscillator with orbital frequency period (dotted line) and the amplitude damping $z(t)$ (blue solid line) plotted over time.</figcaption></figure>

scale height $H = H_0(\frac{R}{R_0})^{(q+3)/2}$. For both tests, we use a global disk model with a radial extent of 20–40 au, $\Delta\theta = 0.72$ radian and $\Delta\phi = 0.4$ radian. The resolution is set to $(96 \times 96 \times 48)$ in $(r, \theta, \phi)$ and $H_0 = 0.1$ at $R_0 = 40$ au. In the following two test runs, we check the vertical and radial integration and drag regimes individually.

### B.1. Vertical motion

The first disk model uses the parameter set of $p = q = 0$ so that the radial pressure gradient vanishes. In this case, the particles do not drift and we can focus on the vertical motion. In this test, we inject a large particle (Stokes number of 300) at two scale heights, located at 30 au. The vertical motion of large particles can be described by a damped oscillation, with oscillation frequency $\omega = \Omega_{K}$ and the amplitude damping $z(t)$

$$z(t) = z_0 \cdot e^{-t/(2\tau_r)}, \qquad (B.3)$$

with the stopping time $\tau_r$. Figure B.1 shows the particle evolution over height. The particle oscillates with orbital rotation frequency and its motion fits a damped oscillation very well.

### B.2. Radial motion

In the second test, we inject particles of different sizes at the mid-plane at $R = 35$ au to study the radial drift. For this setup, we use $p = -1.0$ and $q = -0.5$. We set the particles’ initial velocity to Keplerian rotation. The particle drift is plotted in Fig. B.2. The particles with Stokes number $St = 1$ drift radially inwards. The maximum drift velocity $v_{max}$ is proportional to the velocity difference $v_K - v_{\phi}$ and is, in our case, around 33 meters per second. The drift velocity of different particles sizes can be expressed using

$$v_{\text{dust}} = \frac{-2v_{\text{max}}}{\text{St} + \frac{1}{\text{St}}}. \qquad (B.4)$$

We average the particle velocity between 2 and 8 local orbits. Figure B.2 shows the particle velocity, overplotted with Eq. (B.4). The velocities match the prediction of Eq. (B.4) very well.
<figure><img src="image_11.png" /><figcaption>**Fig. B.2.** Averaged radial particle velocity (black crosses) with Eq. (B.4) overplotted (blue solid line).</figcaption></figure>

## B.3. Benchmark: The radiative transfer

Following the MHD simulation presented in Appendix 2, we performed an additional global simulation with $5 \times 10^6$ particles. For the setup, we use the identical same initial conditions for the MHD setup as the previous model. To distribute the larger particles, we use a $1/R$ profile instead of a uniform distribution. We use the concept, which has been already defined in Appendix A, to incorporate these particles into the radiative transfer simulation. All further parameters remain constant. The gas evolution, and thus the distribution of the small particles, is the same as in the previous models. The larger particles quickly redistribute and are concentrated at the same position. We perform the radiation transfer calculation at the same time output after 150 inner orbits.

The results of both radiative transfer calculations at a wavelength of 1.3 mm are presented in Fig. B.3. Both simulations show a very similar emission map with differences in the order of several percent, see Fig. B.3 bottom. From this, we conclude that the influence of the number of large particles and the initial distribution remains small. This comparison strengthens our conclusion since the structures resulting from the positions of the larger dust particles (see Sect. 3) are robust.

<figure><img src="image_12.png" /><figcaption>**Fig. B.3.** Selected re-emission map at a wavelength of 1.3 mm of a disk model including $5 \times 10^5$ (top), $5 \times 10^6$ (middle) initial larger dust particles and the relative intensity deviations between those two models (bottom). The differences show deviations in the order of several percent.</figcaption></figure>