© ESO 2020

1 Introduction

We have known for a long time that young stars harbor disks consisting of both dust and gas that accrete over time onto their host stars. In recent years with the advent of the Atacama Large Millimeter/sub-Millimeter Array (ALMA) and multiple high contrast optical telescopes (e.g., VLT/SPHERE, Gemini/GPI), further substructures has been found in these protoplanetary disks, such as rings, gaps, crescents, clumps, and spirals (see for a review, e.g., Andrews 2020). In this paper we focus on spirals. Spirals often share a similar morphology: two symmetric spiral arms on both sides of the disk, offset by 180° in an azimuthal direction and often accompanied by radial symmetric gaps close to the spirals. Spirals are found in scattered light (e.g., Fukagawa et al. 2004; Benisty et al. 2015; Akiyama et al. 2016), but also in the millimeter continuum (e.g., Pérez et al. 2016; Huang et al. 2018a; Rosotti et al. 2020) and more recently in CO gas emission (e.g., Kurtovic et al. 2018; Tang et al. 2017) and CO gas dynamics (Teague et al. 2019). This list is far from complete and many questions remain open. For example, in some objects spiral arms are found in gas tracers, but not in the dust continuum, while many spirals in the dust continuum or scattered light do not have a counterpart in the gas.

Some spirals can be explained as the result of the wake created by a massive planet or stellar companion, either located inside or outside the disk, which launches at Lindblad-resonances (see, e.g., Ogilvie & Lubow 2002). These phenomena are very interesting as they allow us to study the mass of the perturbing planet and the properties in the disk that set the shape and morphology of these spirals. The morphology and pitch angle of these spirals in the gas is well known from both theory (Rafikov 2002; Ogilvie & Lubow 2002) and numerical simulations (see, e.g., Dong et al. 2015a; Juhász & Rosotti 2018; Bae & Zhu 2018; Veronesi et al. 2019). However, some spirals can also be explained as being the result of gravitational instability, which happens for high disk masses, specifically when the Toomre Q parameter is lower than one (Rice et al. 2003; Cossins et al. 2010; Dipierro et al. 2014). While the number of spiral arms can vary in the self-gravitating case, for the most massive disks the result in surface density can look very similar to a planet-launched spiral, as in this case they also have two symmetric spiral arms (Dong et al. 2015b; Hall et al. 2018). This type of spiral has logarithmic pitch angles, while planet-driven spirals have linear pitch angles, but in most of the sources the spiral signal is too limited to discriminate between the two mechanisms directly, if a companion is not detected (Forgan et al. 2018). Elias 27 is currently the best candidate for having gravitational instability (GI) as the launching mechanism of its spirals (Meru et al. 2017). One way to discriminate between the two mechanisms can be the gas dynamics, as planet-driven spirals move at the same orbital velocity as the planet, while GI-driven spirals move along with the flow. As it is very time consuming to observe gas dynamics, at the moment of writing only one spiral arm has been observed in the velocity residuals (Teague et al. 2019).

Juhász & Rosotti (2018) show that the pitch angle of planetary spirals is higher in the NIR (near-infrared) than in the sub-millimeter dust emission, as the two tracers come from different layers in the disk at different sound speeds. This also serves as a method to test the planetary hypothesis, since the spirals launched by gravitational instability have no significant difference in pitch angle between observational tracers, as GI tends to erase or invert the temperature differences between the midplane and higher layers as a result of shocks in the midplane (e.g., Boss 2002). An often-made assumption is that the spirals in the ALMA continuum trace the same morphology as the spirals in the midplane gas. However, the fact that some spirals are visible in the dust but not in the gas or vice versa suggests that this assumption is potentially not always valid.

Both disks and mature planets have been observed in great detail, however the formation of planets is difficult to observe directly. Dust-growth models and planetesimal formation mechanisms often assume an initial grain size distribution, as we do not have high precision methods to determine the dust grain size in observations. However, we know that dust particles, which have no pressure of their own, experience radial drift as they feel a headwind from the gas pressure (see, e.g., Whipple 1972; Takeuchi & Lin 2002). This means that in general the dynamics of the dust are different from those of the gas. By analogy, we can expect changes in the morphology and amplitude of the spirals. In this work we focus on this problem, studying the differences between spirals in gas and dust and deriving semi-analytical relations that explain the differences. Surprisingly, up to now there has been no detailed study on the differences between gas and dust spirals. While Veronesi et al. (2019) performed numerical simulations of dust spirals for different Stokes numbers, they did not interpret them in the light of a semi-analytical model as we do in this work. Understanding the dynamics of planet-launched spirals can give us further insight into the planet formation process and can serve as a probe for the Stokes numbers in disks.

The paper is structured as follows. We first discuss the methodology in Sect. 2. We then present the results of the amplitude and morphology differences between gas and dust in Sect. 3, using numerical simulations and the interpretation of these results to derive a semi-analytical model that constructs the spiral in the dust from the deviations from the Keplerian velocity as a result of the spiral in the gas. We test our model on ALMA observations in Sect. 4. We discuss the limitations of our models, alternative explanations, and future prospects in the context of spiral arm detections in gas and dust in Sect. 5. We finally draw our conclusions in Sect. 6.

2 Methods

2.1 Numerical setup

We conducted multi-fluid simulations of the spiral structure formed by the gravitational potential of a planet. We evolved the dust and the gas at the same time using the FARGO-3D code (Benítez-Llambay & Masset 2016), using the Eulerian dust implementation described in Rosotti et al. (2016). The implementation uses the semi-implicit integrator introduced by Booth et al. (2015), but applied to a grid-based code rather than to a particle-based code. FARGO uses the ZEUS numerical algorithm in a co-rotating reference frame.

2.2 Initial conditions

When carrying out simulations, we run FARGO-3D in a 2D cylindrical geometry using an evenly spaced grid, ranging from 0.4 to 3 times the radius of the planet’s orbit, with a resolution of N_r × N_ϕ = 300 × 1000. The planet is kept on a circular orbit at r_pl = 1 and we do not allow it to migrate. We use a locally isothermal model with a constant flaring index, such that the aspect ratio in the disk is given by $$h(r)=\frac{c_{\text{s}}}{v_{\text{k}}}=\frac{H}{r}=h_{\text{pl}}{r}^f,$$(1)

where c_s is the sound speed in the medium,v_k is the Keplerian velocity, H is the pressure scale height, r is the radius in the disk, h_pl is the aspect ratio at the position of the planet, and f is the flaring index. We choose f = 0.25 and an aspect ratio at the planet position of h(r_pl = 1) = 0.05, which are typical numbers, but the value of these constants does not affect the results of this paper.

The gas surface density in the disk is assumed to follow a single exponent power law $$\Sigma (r)={\Sigma }_{\text{0}}{r}^{-p},$$(2)

where Σ₀ is the gas surface density at r_pl = 1 and p is set to 1. The value of Σ₀ is arbitrary as far as the dynamics are concerned and for this reason we show deviations in the surface density relative to the radial symmetric surface density, so that it can be scaled to all values of Σ₀. For the viscosity, we use the α prescription of Shakura & Sunyaev (1973) and fix α_ν = 10⁻². We will discuss the impact of a change in α_ν on our results in Sect. 5.1.

We use 40 dust species that have “St” logarithmically spaced between 10⁻⁴ and 1, where St is the Stokes number or dimensionless stopping time, defined as $$\text{St}=t_{\text{s}}{\Omega }_{\text{k}}=\frac{a{\rho }_{\text{d}}}{{\Sigma }_{\text{g}}},$$(3)

where t_s is the stopping time, Ω_k (r) is the Keplerian angular velocity at radius r in the disk, a is the size of the dust particles, and ρ_d is the bulk density of the dust, which we assume to be 3.6 g cm⁻³ (Love et al. 1994). The latter equality is true assuming that the particles are in the Epstein regime, which is adequate for protoplanetary disks (Birnstiel et al. 2010). We choose the maximum St to be well above the maximum grain size in grain-growth models (Birnstiel et al. 2012) and below the threshold where feedback from the dust on to the gas can be neglected (Rosotti et al. 2016). We choose the minimum St to be small enough such that the dust can be considered closely coupled to the gas as we will show in the next section, so for the purpose of this paper there is no need to take smaller St into account. The dust surface density is scaled with respect to the gas surface density using a constant dust-to-gas ratio of 100. Because we ignore the dust feedback on the gas, our simulations can be rescaled to any dust-to-gas ratio as long as it is small.

We modeled the spiral using three different planet masses: Earth-mass planet M_E = 3 × 10⁻⁶ M_*, super-Earth-mass planet M_SE = 10⁻⁵ M_* and Jupiter-mass planet M_J = 10⁻² M_*. The planet potential is delayed for two orbits and is gradually added to the model using a taper over 3 orbits. This prevents high frequency artifacts in the data and makes sure that the spiral sets quickly in time. We checked that the spiral morphology is in a steady state, which is the case after about 10 orbits, so the result after 15 orbits of integration is used in all further analyses.

2.3 Boundary conditions

For the gas simulations we used closed boundaries where the radial velocity is mirrored at the boundaries and the azimuthal velocity and surface density are scaled using the closest active cell. For the dust simulations we used open boundary conditions to avoid the dust piling up at the inner boundary due to radial drift of the larger dust grains. To make sure that the surface density stays constant over time and to prevent gradual depletion of the large dust grains, the radial velocity at the boundary is set to the radial drift velocity, as given in Takeuchi & Lin (2002), $$v_{r,\text{d}}=\frac{\eta v_{\text{k}}}{\text{St}+{\text{St}}^{-1}}\propto {\left(\frac{H}{r}\right)}^2{v}_{\text{k}}\propto r^{2f-1/2},$$(4)

where η = $-{\left(\frac{H}{r}\right)}^2\frac{\text{dlog}(P)}{\text{dlog}(r)}$ with P being the gas pressure.

In what follows we find it useful to define a linear pitch angle and we measure it using a linear fit assuming an Archimedean spiral ϕ = ar + b. The pitch angle of the spiral can then be written as β = tan⁻¹(a∕r). The spiral structure is analyzed at a radius of 1.7 r_pl throughout the paper, unless mentioned explicitly. This radius is well outside a potential gap opened by the planet, far from any possible boundary artifacts and in the regime where the spiral is approximately Archimedean.

Fig. 1 Synthetic FARGO-3D images of spiral arms created by the potential of a super-Earth-mass planet in the gas component (panels a–c) and in the large dust grain component with St = 1 (panels d–f). The images show the perturbation of the spiral normalized to the unperturbed disk. From left to right: normalized perturbation of the surface density (a and d), the perturbation in radial velocity with respect to the Keplerian velocity (b and e), and the normalized perturbation in azimuthal velocity (c and f).

3 Spiral modeling

We split theresults of our modeling up into three parts. In the first part we explore the structural changes of the spiral in different dust components and compare this to the spiral in the gas component. In the second part we derive semi-analytical relations to derive the density structure of the dust from perturbations in the Keplerian velocity of the gas component and vice versa. In the third part we combine the relations to explore the physical mechanisms that drive the changes in structure of the spiral in different dust components and show a model that accurately determines the spiral perturbation in the dust surface density from gas dynamics.

3.1 Comparison between spirals in gas and dust

In Fig. 1 we show the results of modeling spiral arms created by the potential of a super-Earth-mass planet using FARGO-3D for the gas and for the biggest dust component (St = 1). We show the amplitude of the spiral in surface density (δΣ), normalized to the surface density of the unperturbed disk (Σ₀), and the two components of the velocity perturbation, azimuthal velocity perturbation (δv_ϕ) and radial velocity perturbation (δv_r), normalized to the Keplerian velocity. In terms of the general structure, the spiral perturbation in surface density and the perturbation in radial velocity look largely similar, with the only significant difference being that the radial velocity perturbation changes signat the planet’s position. This is because the spiral moves with the orbital speed of the planet, which means that the gas overtakes the spiral in the inner disk, but the spiral overtakes the flow in the outer disk. A more detailed description of this can be found in Sect. 3.3. The perturbation in azimuthal velocity has a different structure and is almost an order of magnitude smaller than the perturbation in radial velocity. The structural difference between the spiral perturbation in the gas and the dust is small, but we see a clear difference in amplitude in all three components. The results are similar for planets with a different mass, but high mass planets induce a second spiral shifted 180° with respect to the first spiral.

In Fig. 2 we show the azimuthal cross section of the spiral at a radius of 1.7 r_pl in the three components of the perturbation for different St. The direction of the flow with respect to the spiral is shown from right to left. The shape of the spiral along the azimuth is largely similar in the surface density (panel a) and the radial velocity (panel b), but differs from the azimuthal velocity (panel c). The peak of the perturbation in the dust components with high Stokes numbers is shifted with respect to the peak of the perturbation in the gas, but the peak always lies on top of the curve of the perturbation in the gas (see Fig. 2). In the azimuthal velocity image we can qualitatively interpret this using the fact that components with high St are no longer closely coupled to the gas, which means that dust particles react in a delayed way to the velocity difference in the spiral: decelerating as long as the gas velocity is lower (from (1) to (2) in the right panel of Fig. 2), but never catching up before the gas speeds up as it leaves the spiral (from (2) to (3)). This results in a similar pattern in the radial velocity and surface density, which we will explain in more detail in the next section.

In Sect. 3.3 (see Eq. (14)) we will analyze further the different components of the velocity perturbation. For now, it is useful to consider the stream lines that a gas parcel would follow during an orbit in the vicinity of r = 1.7r_pl. We used a first order Eulerian integrator to determine the stream lines and we have plotted the results in the bottom panel of Fig. 2 with radial position as a function of position angle. We also introduce a new quantity, the line density perturbation. This is defined as the density of streamlines in the radial direction. The line density perturbation has the same shape and amplitude as the perturbation in surface density, showing that it is the critical quantity sculpting the spiral. The line density is related to the amount of “squeezing” in the radial direction, that is, to the radial velocity perturbation. A detailed interpretation of the line density and a toy model explaining this more in detail are presented in Appendix A.

In Fig. 3 we show the location of the dust maxima along the azimuthal direction for different St between a radius of 0.5 and 1.8 r_pl on top of a polar plot of the surface density perturbation in the gas. At a radius of 0.7 r_pl the maximum of the spiral changes side of the perturbation, causing an azimuthal jump of the maximum.

Fig. 2 Top, from left to right: azimuthal slab through the normalized perturbation in azimuthal velocity (a), radial velocity (b), and surface density (c) for Stokes numbers 10−4–1. The perturbation in the gas is indicated in black as a reference. Bottom: streamlines that a gas parcel would follow along an orbit at r = 1.7 rpl. The maximum radial deviation is marked with a black dot as a guide for the eye to the effect of the pitch angle on the flow of the particles. The perturbation in line density (see text) is indicated in red, sharing the same shape and amplitude as the perturbation in surface density (panel c).

Fig. 3 Polar plot of the perturbation in surface density over two orbits. The lines indicate the position of the maximum of the different dust components with a zoom of the region around r = 1.7 rpl that is used to determine the change in pitch angle.

3.2 Spiral structure as function of St

Morphology

In the left panel of Fig. 4 we show the change in pitch angle of the spiral as a function of the Stokes number, representingdifferent particle sizes of the dust, using the results of our simulations. The linear pitch angle is measured using a fit ϕ = ar + b of the position of the maxima along the position angle (PA) between 1.6 r_pl < r < 1.8 r_pl (see Fig. 3). There is some difference in the relative change of the pitch angle tan⁻¹dϕ∕(rdr) of the spiral for St > St_c (see next paragraph for a formal definition of St_c), although we find that even for St = 1, the change in dϕ/dr is only a few percent. This is less than the constraints we usually can put on the pitch angle due to the uncertainty in disk height and geometry, as well as observational errors. We find instead that there can be a significant azimuthal offset (20° for St = 1) between gas and dust spirals. The offset is higher for Jupiter-mass planets.

Fig. 4 Left panel: change in pitch angle with respect to the pitch angle of the spiral in the gas, fitted using the function ϕ = ar + b between r = 1.6,1.8 rpl as a function of the Stokes number. The solid lines represent the change in pitch angle relative to the pitch angle in the gas and the dotted lines show the constant azimuthal offset between the spiral in the gas and the dust in degrees. Right panel: amplitude of the spiral versus the Stokes number at r = 1.7 rpl, normalized by the amplitude of the spiral in the gas.

Amplitude

In the right panel of Fig. 4 we show the maximum amplitude of the spiral at a radius of 1.7 r_pl as a function of the Stokes number. The maximum spiral amplitude is for each planet mass normalized to the maximum amplitude of the spiral in the gas component for a better comparison. We find that the shape of the curve is almost identical for the three planet masses. We note, however, that the absolute amplitude of the spiral increases if the mass of the planet increases. The amplitude of the spiral in the surface density of the disk does not change for dust grains with St ≲ 0.05. This means that the amplitude of the spiral is constant for almost all dust particle sizes that are available in grain growth models. To interpret this result we introduce a naive critical Stokes number, determined by the point where the stopping time of the dust is of the same order as the crossing time through the spiral, $$t_{\text{cross}}=\left(\Delta \phi +{\Omega }_{\text{p}}{t}_{\text{cross}}\right)t_{\text{dyn}}=\frac{\Delta \phi t_{\text{dyn}}}{1-{\Omega }_{\text{p}}{t}_{\text{dyn}}},$$(5)

where Ω_p is the angular velocity of the planet and hence of the spiral, t_dyn is the dynamical time, and Δφ is the width of the spiral with respect to a whole revolution. The second term arises due to the fact that the spiral rotates slightly during crossing over. Using Eq. (3) the critical Stokes number can be estimated as $${\text{St}}_{\text{c}}=\frac{\Delta \phi }{(1-{\Omega }_{\text{pl}}{t}_{\text{dyn}})}.$$(6)

Estimating the width of the spiral using a Gaussian fit, we get St_c = 0.057, 0.061, and 0.063 for M_E, M_SE, and M_J respectively.The mean of these three numbers is shown in Fig. 4 as a vertical dotted line, which matches accurately with the turning point of the graph. The critical Stokes number depends on r, but for r far enough from the planet’s position, the shape of the curve and the critical Stokes number are very similar to the result shown (see Appendix C).

3.3 Semi-analytical description spiral shape

Using analytical relations we are able to determine the perturbation in surface density in the dust from the perturbation in azimuthal velocity in the gas alone. Three substeps are needed to interpret the results in Fig. 4 and understand the dynamics of dust while crossing the spiral:

• the relation between the perturbation in azimuthal velocity and radial velocity;

• how the perturbation in surface density arises from the perturbations in velocity;

• how the perturbations in the velocity of the dust are related to the perturbations in the velocity of the gas.

3.3.1 From δ v_ϕ to δ v_r

We can write the surface density as an azimuthal symmetric part and a small perturbation Σ = Σ₀ + δΣ and the corresponding velocity field as v = v₀ + δv, v_r = v_r,0 + δv_r, and v_ϕ = v_ϕ,0 + δv_ϕ, where v_r is the radial component of the velocity and v_ϕ is the azimuthal component of the velocity. Assuming that δv_ϕ ≪ v_ϕ,0 and using the fact that a perturbation in Keplerian velocity results in an immediate change of angular momentum, the effect of the perturbation in v_ϕ can be explained as spiraling in or outwards on circular orbits: $$v_{\varphi }\simeq r{\Omega }_{\text{k}}=\sqrt{\frac{GM}{r}}$$(7) $$v_r=\frac{\text{d}r}{\text{d}t}=\frac{\text{d}}{\text{d}t}\frac{GM}{v_{\varphi }^2}=GM\frac{\text{d}\varphi }{\text{d}t}\frac{\text{d}{v}_{\varphi }^{-2}}{\text{d}\varphi }=GM\frac{v_{\varphi }}{r}\left(\frac{-2}{v_{\varphi }^{3}r}\frac{\text{d}{v}_{\varphi }}{\text{d}\varphi }\right).$$(8)

Assuming v_r,0 = 0 and using the fact that ∂_ϕv_ϕ,0 = 0, the perturbation in v_r is related to the derivative of v_ϕ, $$\delta v_{\text{r}}=-2\text{sgn}(r-r_{\text{pl}}){\partial }_{\varphi }\delta v_{\varphi },$$(9)

where the sgn(r − r_pl) term arises from the fact that the spiral moves with the same speed as the planet, which means that outside the planet orbit the spiral is moving faster than the gas, but inside the planet orbit the gas will overtake the spiral, causing a sign flip in the derivative.

3.3.2 From δv to δΣ

When the density perturbation is small compared to the average density, we can determine a relation between the density perturbation and the perturbation in gas dynamics caused by the planet potential. We know that both Σ and Σ₀ have to obey the continuity equation, resulting in the following set of equations: $${\partial }_t{\Sigma }_{\text{0}}+v_{\text{0}}\cdot (\nabla {\Sigma }_{\text{0}})+{\Sigma }_{\text{0}}(\nabla \cdot v_{\text{0}})=0$$(10) $${\partial }_{\text{t}}({\Sigma }_{\text{0}}+\delta \Sigma )+(v_{\text{0}}+\delta v)\cdot (\nabla (\Sigma +\delta \Sigma ))+(\Sigma +\delta \Sigma )(\nabla \cdot (v_{\text{0}}+\delta v))=0.$$(11)

Neglecting quadratic terms in perturbation and assuming that v_r,0 = 0, we find $${\partial }_{\varphi }\delta \Sigma =-\frac{r{\Sigma }_{\text{0}}\delta v_{\text{r}}\left(\frac{{\partial }_r{\Sigma }_{\text{0}}}{{\Sigma }_{\text{0}}}+\frac{{\partial }_r\delta v_r}{\delta v_r}+\frac{1}{r}\right)+{\Sigma }_{\text{0}}{\partial }_{\varphi }(\delta v_{\varphi })}{v_{\varphi ,\text{0}}}.$$(12)

This can be simplified using Eqs. (2) and (9) to $${\partial }_{\varphi }\delta \Sigma =-\frac{r{\Sigma }_{\text{0}}\delta v_{\text{r}}}{v_{\varphi ,\text{0}}}\left(\frac{{\partial }_r\delta v_{\text{r}}}{\delta v_{\text{r}}}+\frac{1}{2r}-r^{-p}\right),$$(13)

where the first term is dominant in all cases, as we will show in Sect. 3.4.

The perturbation in surface density has the same structure as the perturbation in radial velocity, as we illustrate in Fig. 2. This can be interpreted using the fact that the spiral perturbation can be described as a 1D line in {r,ϕ}, as long as the azimuthal shape of the spiral does not change significantly in a radial cut through the spiral. The relation between taking the radial derivative and the azimuthal derivative is then given by the pitch angle of the spiral, d ϕ∕dr = a, with a the linear pitch angle. Using this assumption Eq. (13) reduces to a simple scalar relation between δΣ and δv_r : $$\delta \Sigma =-\text{sgn}(r-r_{\text{pl}}){\Sigma }_{\text{0}}r\frac{\delta v_{\text{r}}}{v_{\varphi ,\text{0}}}\frac{\text{d}\varphi {(r)}_{\text{sp}}}{\text{d}r}.$$(14)

This result can be interpreted using the bottom panel of Fig. 2, where we show stream lines at a few orbits close together that a gas parcel would follow during an orbit around the central object. Since the positional changes due to the spiral potential are tiny, the azimuthal shape of the spiral can be interpreted as being locally independent of the orbit radius. The only change between stream lines that are close together is the small azimuthal offset associated with the pitch angle of the spiral. This results in a vertical squeezing that is related to the amount of radial change, that is, the radial velocity. A simple sketch for a detailed explanation of this relation is given in Appendix A.

3.3.3 From gas to dust

The velocity perturbation of the dust can be derived from the dynamics of the gas by solving the azimuthal equation of motion for the dust. The only force acting on the dust is the drag force caused by the gas moving at slightly sub-Keplerian velocity due to the gas pressure. The equation of motion for the azimuthal velocity of the dust is hence given by $$v_{\varphi ,\text{d}}\frac{\text{d}{v}_{\varphi ,\text{d}}}{\text{d}\varphi }=-\frac{r}{t_{\text{fric}}}(v_{\varphi ,\text{d}}-v_{\varphi ,\text{g}}),$$(15)

where the subscript d is the dust component and g is the gas component. Rewriting this in terms of the Stokes number using Eq. (3), we can write the perturbation in the azimuthal velocity of the dust as the first order nonlinear ordinary differential equation $$\frac{\text{d}\delta v_{\varphi ,\text{d}}}{\text{d}\varphi }=\frac{r}{\text{St}{\Omega }_{\text{k}}}\left(\frac{v_{\varphi ,\text{g}}(\varphi )}{\delta v_{\varphi ,\text{d}}+{\Omega }_{\text{k}}r}-1\right).$$(16)

We approximated this differential equation with the standard Runge-Kutta method using the integrate.solve_ivp function in scipy (Virtanen et al. 2020). We integrate for three orbits in steps of 0.01 radians, using an initial value for the perturbation of 0 at all radii.

3.4 Combining everything

Combining the equations in Sect. 3.3 makes it possible to determine the shape and amplitude of the spiral in dust surface density directly from the azimuthal velocity perturbation in the gas. This model has three steps:

• 1.

  determine the perturbation in azimuthal velocity of the dust from the azimuthal velocity of the gas using Eq. (16);

• 2.

  determine the perturbation in radial velocity of the dust from the azimuthal velocity of the dust using Eq. (9);

• 3.

  determine the perturbation in surface density of the dust from the perturbation in radial velocity of the dust. This can be done using two different approaches:

• using the full perturbed continuity equation (Eq. (13));

• using the perturbed continuity equation and numerical derivation of the spiral position to determine the spiral pitch angle (dϕ/dr) as a function of the radius to relate the azimuthal derivative to the radial derivative (Eq. (14)).

The accuracy of the different substeps compared to the FARGO-3D output is given in Appendix B. In Fig. 5 we show the end results in the most extreme case, St = 1. In the top panel of the figure we show the normalized amplitude of the spiral throughout the whole disk as determined using the aforementioned procedure and on the bottom panel of the figure we show an azimuthal slice through all results as a comparison. The shape of the perturbation matches accurately with the predicted result from FARGO-3D for r > 1.5 r_pl. We find that the two approaches to determine the surface density from the perturbation in radial velocity give almost identical results, justifying the assumption that the perturbation in radial velocity sets the perturbation in surface density as explained in Sect. 3.3.2. In general, the surface density perturbation is more accurately recovered closer to the planet for smaller Stokes numbers. When Stokes numbers are low enough to consider thedust closely coupled to the gas, the surface density does not change with respect to the gas and the right hand side of Eq. (16) approaches zero, not adding any additional uncertainties to the final result.

Combining the equations in Sect. 3.3 gives a full understanding of the shape of the amplitude and pitch angle of the curves as a function of the Stokes number (Fig. 4). Using Fig. 2c we find that dust particles feel the change in the amount of drag when the azimuthal velocity of the gas changes (point (1) as indicated in the figure). Only very large dust particles with a high Stokes number have a reaction time that is longer than the crossing time through the spiral (Eq. (6)), so they accelerate and decelerate less than the gas. Since the dust particles leave the spiral before they are able to catch up (point (2) and (3) as indicated in Fig. 2), the amplitude of the spiral in azimuthal velocity decreases towards larger Stokes numbers. A smaller perturbation in azimuthal velocity with a damped acceleration and deceleration means a smaller perturbation in radial velocity (Eq. (9)), hence in surface density (Eq. (14)). A consequence of this mechanism is that the peak in the azimuthal velocity of the dust is always the intersection point with the azimuthal velocity of the gas, as the dust starts to decelerate after that point (see Fig. 2). This explains the small change in pitch angle of the spiral for larger St (Fig. 4). Since the spiral decreases in amplitude over radius and broadens in azimuthal shape, the position of the spiral in the dust is not changed by only a constant but the pitch angle changes as well.

Fig. 5 Top: azimuthal slice through the normalized density perturbation at a radius of 1.7rpl for the dust component with St = 1. The black line shows the time integrated FARGO-3D results, the blue line shows the results directly derived from the perturbation in azimuthal velocity of the gas using the perturbed continuity equation (Eq. (13)), the red line shows the result directly derived from the perturbation in the azimuthal velocity of the gas using a direct scaling between δvr and δΣ and the pitch angle of the spiral (Eq. (14)). The dotted line illustrates the density profile of the gas component for comparison. Bottom: normalized density perturbation in polar coordinates of the result for the St = 1 component. From left to right: time integrated FARGO-3D result, the results directly derived from the perturbation in the azimuthal velocity of the gas, and the result directly derived from the perturbation in azimuthal velocity of the gas, using a direct scaling between δvr and δΣ and the pitch angle of the spiral (Eq. (14)).

4 Testing on observational data

4.1 Source selection

Our models show that high signal to noise ratio (S/N) spiral structures in the gas dynamics, together with a spiral seen in the dust continuum surface density, can put constraints on grain size comparing the observed amplitude of the spiral in the dust with that in the gas. There are multiple candidate sources that are known to have spirals in the continuum or the gas. However, to give an accurate estimate of the surface density and to derive a velocity residual map, we need data with high spatial and velocity resolution. Out of the recent high-resolution Disk Substructures at High Angular Resolution Project (DSHARP) survey of 20 disks, 3 were found to harbor spirals: Elias 27, WaOph 6, and IM Lup (Huang et al. 2018a). Elias 27 harbors the highest contrast spiral arms of the survey, but gravitational instability is a convincing explanation for this source (Pohl et al. 2015; Meru et al. 2017; Forgan et al. 2018, although see also Hall et al. 2018), which makes it unusable for this purpose, as further discussed in Sect. 5.3. Furthermore, Elias 27 and WaOph 6 suffer from cloud and outflow contamination in the CO gas (Andrews et al. 2018), which makes it impossible to use the kinematic data published so far. The IM Lup disk has two high contrast spiral arms in the millimeter continuum, and no cloud contamination (Andrews et al. 2018), so is an ideal candidate to analyze. Unfortunately, no spiral signal is detected in the gas for the IM Lup disk, but in Sect. 4.3 we will show that we can make a prediction on what the spiral looks like in the two velocity components. Other sources are considered by surveying the literature. MWC 758 is a promising source (Dong et al. 2018; Boehler et al. 2018), but the observations are currently too limited to be used. HD 100453 has a spiral detected in scattered light, continuum, and CO gas (Benisty et al. 2017; Rosotti et al. 2020), but the spiral seen in the gas is beyond the size of the continuum disk and the inner region of the disk has a noisy velocity map, potentially due to a warp in the inner disk, which makes the source not useful in our analysis. TW Hya harbors the only spirals in gas dynamics, potentially in the azimuthal velocity (Teague et al. 2019). The spirals are detected in the CO gas emission as well, but no spiral signal is present inthe dust continuum images. In the next subsection, we will analyze the TW Hya disk further and show that the observed spiral signal in the dynamics is interesting for further study, but not useful to test our model.

4.2 TW Hya disk

TW Hya, located at a distance of 60.1 pc, hosts a well-studied protoplanetary disk for which a potential spiral in azimuthal velocity is suggested by Teague et al. (2019). However, using the exact same data set we found a different spiral than originally suggested inthe paper, which radically changes the interpretation of the velocity deviations. A detailed description of the observations, calibration, and further processing of the data can be found in the original papers (Huang et al. 2018b; Teague et al. 2019). In Fig. 6 we show the exact same data as in Fig. 1 of the original paper, but to a further extent, and derotated such that the semi-minor axis of the disk is vertical. The color bar that was used emphasizes the residual positive velocity, but the residual negative velocity can not be neglected in this case. Flipping the color bar, now emphasizing the blueshifted emission, we find a different spiral that has comparable morphology to one of the spirals found in the CO gas emission. The spiral is constant in amplitude along the azimuth and appears to have a kink at the position of the gap in the continuum emission at r = 1.65″ and PA = 160°. The shape of the spiral is consistent with the theoretical shape given by the wake equation of Rafikov (2002), shown as a dotted line in Fig. 6: $$\begin{array}{ll}\varphi (r)=\hfill & {\varphi }_{\text{pl}}-\frac{\text{sgn}(r-r_{\text{pl}})}{h_{\text{pl}}}{\left(\frac{r}{r_{\text{pl}}}\right)}^{1+\beta }\left\{\frac{1}{1+\beta }-\frac{1}{1-\alpha +\beta }{\left(\frac{r}{r_{\text{pl}}}\right)}^{-\alpha }\right\}\hfill \\ \hfill & +\frac{\text{sgn}(r-r_{\text{pl}})}{h_{\text{pl}}}\left(\frac{1}{1+\beta }-\frac{1}{1-\alpha +\beta }\right)\hfill \end{array}$$(17)

where α is the power exponent of the rotation angular frequency (Ω(r) ∝ r^−α), and β is the power exponent of the radial distribution of the sound speed (c_s ∝ r^−β). We set α to the typical value 1.5 and β to 0.3 as determined by Kama et al. (2016) to model the disk. Using a fit by eye, we find h_pl = 0.06 as a best fit, which is close to the value used in Kama et al. (2016) to model the disk ${(h/r)}_{\text{pl}}~$0.1. The fact that we see the whole spiral as blueshifted relative to the disk means that the spiral is a perturbation in vertical velocity, since both the radial velocity as well as the azimuthal velocity perturbations change sign along either the major or minor axis of the disk (see Teague et al. 2019, for a schematic of the different velocity components in a velocity residual map). These two clues provide potential evidence that a planet is carving out the gap, inducing a spiral that stirs the material up to larger heights in the disk. The spiral found in the original paper could be part of the sub-Keplerian rotation that we observe due to the fact that we look deeper inside the disk at the location of the gap cleared by the planet. The study conducted in this paper is 2D; 3D modeling would be required to study the relation between spiral patterns in the vertical component of the velocity and the in-plane radial and azimuthal velocity. This could give more insight into this particular case, but that is beyond the scope of this paper.

Fig. 6 Map of the velocity residual as in Teague et al. (2019) for the disk around TW Hya, but with a flipped color bar to highlight the blueshifted emission (now seen in red). The dotted line serves as a guide for the eye and is determined using the theoretical shape of a spiral launched by a planet (Eq. (17)). The solid line represents the radius of the potential planet’s orbit, which coincides with the blue- and redshifted circular lobes on either side that indicate sub-Keplerian rotation, probably due to a lower opacity inside the gap, so we probe deeper layers in the disk.

4.3 IM Lup disk

We can use the relations derived in Sect. 3 on observational data as a prediction of what spiral arms will look like in different velocity components. As we have shown in Sect. 4.1, the best target for this is the IM Lup disk as it has a high contrast spiral signal and no cloud contamination. IM Lup is a 0.5 Myr old K5 star located at a distance of 158 ± 1 pc (Alcalá et al. 2017; Gaia Collaboration 2018). The 1.25 mm continuum and ¹²CO J = 2–1 gas emission of this disk has been observed as part of the DSHARP survey. A detailed overview of the survey, including the observational setup, calibration, and imaging is provided in Andrews et al. (2018); a detailed description of the spirals in the millimeter continuum is provided in Huang et al. (2018a).

After deprojection of the disk using a position angle of 144. °5 and an inclination of 47. °5 (Huang et al. 2018c), we estimate Σ₀ by taking the median value of the surface density along the azimuth and the spiral perturbation (δΣ) as the deviation from this value in the data. In Fig. 7 we show the spiral perturbation in the millimeter dust continuum surface density (left panel) with the derived perturbation in the radial velocity (middle panel) and azimuthal velocity (right panel), derived using Eqs. (9) and (14). The pitch angle of the spiral that is indicated in Fig. 7 and used in the model, is determined using a fit by eye through the spiral maximum in the surface density perturbation, assuming that the two parts of the spiral separated by a radial symmetric gap are part of one spiral arm. Using a parametrization of ϕ = ar + b we find ϕ = 590[°]r[′′] – 12[°] and ϕ = 590[°]r[′′] – 192[°], which is consistent with the results obtained in Huang et al. (2018a). There is evidence of a gap in the continuum emission at a radius of 0.70″ (Huang et al. 2018a) and a potential kink in the velocity channel maps at the same radius, with a potential planetary origin (Pinte et al. 2020), so we assume that r < r_pl in the region we analyze, adding an additional minus sign in Eq. (14).

The azimuthal velocity perturbation is two or three times smaller than the perturbation in radial velocity and is shifted ~30° with respect to the surface density perturbation. To compare our results with the measured dynamical perturbation, we used the bettermoments package described in Teague & Foreman-Mackey (2018) to determine the velocity map of the CO gas emission. The velocity map is determined using a fit of a quadratic curve to the three pixels closest to the peak intensity, providing a higher spectral precision than the velocity resolution of the data. The velocity perturbations can be determined from this map by subtracting a Keplerian rotation profile. To do so we fit a Keplerian model to the data using eddy (Teague 2019), where the line center is given by $$v_{\text{k}}=\sqrt{\frac{G{M}_{\ast }{r}^2}{{(r^2+z^2)}^{3/2}}},$$(18)

where z is the height in the disk using a power law profile with a power law cut off in the outer disk, $$z(r)=z_{\text{0}}{r}^{\psi }+z_{\text{1}}{r}^{\phi }.$$(19)

As the inclination and the mass of the star are degenerate with each other, we fix the inclination to the same value as found in the continuum. The resulting parameters of the fit can be found in Table 1. In order to compare our computational results with the observations, we project the velocity residuals to what we would be able to see with the known disk inclination. The radial velocity perturbation is projected as δv_r,proj = δv_rcos(PA + 90°)sin(i) and the azimuthal velocity perturbation as δv_ϕ,proj = δv_ϕcos(PA)sin(i) with PA the position angle in the disk and i the inclination. The results are illustrated in Fig. 8, where we show a projected estimate of the azimuthal velocity perturbation derived from the dust continuum emission using Eqs. (9) and (14) (left panel) and an estimate of the radial velocity perturbation using Eq. (14) (central panel). In the right panel of Fig. 8 we show theobserved velocity residual map. Our estimate in radial velocity is of the same order of magnitude as the residuals in the observed velocity map. However, the data contain “spokes” artifacts that are due to the velocity channelization of the data. With the current velocity resolution, it is therefore impossible to deduce any spiral structures, highlighting the need for data with higher velocity resolution.

Fig. 7 From left to right: observed surface density perturbation in the millimeter dust continuum of the IM Lup disk normalized to the radial symmetric surface density (a), the perturbation in radial velocity normalized to the Keplerian velocity, derived from the surface density perturbation using Eq. (14) (b), and the azimuthal velocity perturbation normalized to the Keplerian velocity derived from δvr,calc using Eq. (9) (c). The dashed lines represent a fit by eye through the maximum of the spiral in the surface density, which serves as a guide for the eye.

Fig. 8 Spirals in gas and dust for the IM Lup disk. From left to right: estimated perturbation in azimuthal velocity from the dust continuum using Eq. (9) and Eq. (14) (a), the estimated perturbation in radial velocity from the continuum using Eq. (14) (b), and the observed velocity residuals in the optically thick CO gas emission after subtracting a Keplerian model of the data (c).

5 Discussion

We have run FARGO-3D simulations and analyzed the spiral created by the wake of a planet in the different dust components. We find that dust grains with typical Stokes numbers (e.g., St = 0.1) such that they radially drift fairly quickly (see, e.g., Takeuchi & Lin 2002; Whipple 1972) are closely coupled to the gas in terms of the spiral morphology. This is due to the fact that the azimuthal width of the spiral is much smaller than a whole orbit, which means that the dust grains have to decouple more from the gas in order to react to the perturbation caused by the planet.

5.1 Importance of α_ν and r

Viscosity

Throughout the semi-analytical analysis we used a value of α_ν = 10⁻² for the viscosity to test and compare the derived relations. All relations used are independent of α_ν, although the shape of the spiral along the azimuth changes when using lower values for the viscosity; we chose a relatively high value for α_ν to smooth out high frequency oscillations along the position angle. Changing α_ν does not change the pitch angle of the spiral in the gas, but has an effect on the amplitude and width of the spiral as the perturbation diffuses out more efficiently with higher viscosity. This can slightly affect the pitch angle and amplitude of the spiral in the dust components. Once normalizing to the properties of the spiral in the gas (see Fig. 4), the shape of the pitch angle and spiral amplitude as a function of the Stokes number do not change significantly.

Radius

In order to make sure that no boundary effects play a role in our analysis, we analyzed the spiral at r = 1.7 r_pl and ran the code once with a much larger extent of 0.25 r_pl < r < 5 r_pl (N_ϕ × N_r = 1024 × 512), but apart from regions close to the boundary, this does not change the morphology of the spiral in the analyzed part of the disk (0.4 r_pl < r < 3 r_pl). We analyzed the spiral at a radius of 1.7 r_pl, far enough from the planet to avoid the effects of gap clearing and the direct gravitational effects of the planet. The model is unable to recover accurately the spiral in the dust close to planet, mainly due to nonlinearity of the pitch angle. The radial dependence of the model is small in the regime where the pitch angle is linear and the critical Stokes number changes by only a factor of a few over the radius (see Appendix C), which finds its main origin in the diffuse broadening of the spiral.

5.2 Small versus large dust grains

Rosotti et al. (2020) showed that the pitch angles of spirals in HD 100453 (~6°) are lower in the sub-millimeter dust continuum, which traces large dust grains, than in scattered light, which traces small dust grains (~16°). Using hydrodynamical simulations they show that this is a general property of spirals as a result of different sound speeds, due to the fact that for externally irradiated disks the midplane is colder than the upper layers of the disk. This mechanism can create differences in the morphology of observed spirals between gas and dust as well. One assumption they made in the modeling is that the spirals in the ALMA dust continuum image trace the same morphology as the spirals in the midplane gas. In a typical disk, the critical Stokes number when gas and dust spirals start to deviate, determined in Sect. 3.2, corresponds to centimeter sized dust grains, well above the maximum grain size that can be reached due to the bouncing and fragmentation barriers (e.g., Birnstiel et al. 2010; Windmark et al. 2012). This means that the assumption holds that the spirals in the ALMA dust continuum trace the same morphology as the midplane gas.

5.3 Spiral origin

Spiral perturbations in protoplanetary disks can be created in two different ways. In this paper we focused on spirals that are created by Lindblad resonances in the disk, generated by planets. This differs significantly from spirals caused by the gravitational instability of massive disks, as these instabilities are thought to move with the flow. This means that self-gravitating spirals can trap dust (e.g., Rice et al. 2004, 2006; Booth & Clarke 2016) and their amplitude increases with the Stokes number. Instead, in the planet case, spirals move at the same angular velocity as the planet throughout the disk; spiral crossing is short (compared to the orbital timescale) and there is no particle trapping. In this case, the amplitude of the spiral decreases with the Stokes number, as we have shown in this paper. In principle, this difference can be used in observations to discriminate between the two origins using continuum multi-wavelength observations.

In the IM Lup system there is evidence of a gap in the continuum emission at a radius of 0.70″ or 111 au (Huang et al. 2018a) and a potential kink in the velocity channel maps at the same radius with a potential planetary origin (Pinte et al. 2020), which makes the assumption that the spiral is planet-launched reasonable. The IM Lup disk mass (M_disk∕M_*) is considered too small to cause gravitational instability (e.g., Kratter & Lodato 2016), although there is some debate on the origin of the spiral (Huang et al. 2018a). If the spirals we observe in the continuum of IM Lup are generated by gravitational instability, our model would not work, and the estimate of the projected dynamics derived from the surface density would be incorrect. In addition, to differentiate between grains of different sizes, data with a better velocity resolution can be used to end this debate by comparing the velocity perturbations from Keplerian rotation with our estimates based on the planetary origin.

5.4 Limitations and future prospects

In this paper we did not consider the impact of dust evolution on the local dust density structures. Dust grain coagulation and fragmentation can potentially have a strong impact on the local dust density structures as well (see, e.g., Drążkowska et al. 2019). However, in contrast to, for example, gaps and rings, the spiral structure is established in only a few orbits, while grain growth is a significantly slower process: measured in orbits, the collisional timescale between two grains is of order of the gas-to-dust ratio (e.g., Birnstiel et al. 2012). It is thus a reasonable assumption to take the grain properties as fixed as we have done in this paper. This does not exclude the possibility that, on secular timescales, grain growth will modify the grain properties and the spiral will readjust accordingly, following the relations we have derived in this work. The same argument holds for another effect we neglected in this paper, namely planet migration. The spiral reaches steady state in a few orbits, which is fast compared to the migration timescale (hundreds or more orbits, e.g., Baruteau et al. 2014). This justifies the assumption made in Sect. 3 that the planet moves on a fixed circular orbit.

Another effect we neglected in this paper is the feedback of the dust onto the gas. This would not be justified close to the pressure maximum created by the planet, where the dust accumulates and dust-to-gas ratio becomes high. However, in our analysis we focused on the spiral structure further from the planet, where there is no accumulation effect. It should be noted that recently Dipierro et al. (2018) pointed out that the effect of dust feedback could be non-negligible even in smooth parts of the disk, away from pressure maxima, with the magnitude of the effect depending on the value of α_ν (though see, e.g., Rosotti et al. 2019 for calculations including the effect of dust feedback where this is not seen). If this is confirmed infuture studies, the effect of dust feedback on planetary spirals would need to be re-assessed. Our paper provides the starting point to build a theoretical description including dust feedback.

Since the ¹²CO line is optically thick, it is useful to observe the dynamics of rarer isotopologues to study the in-plane motions in the spiral, as the in-plane perturbation fades out higher up in the disk. However, studying the dynamics of, for example, C¹⁸ O does take a considerable amount of observing time. Understanding velocity perturbations in the vertical direction can be useful to get a better understanding of the mechanism that drives the spiral in the gas, for example in TW Hya. Three-dimensional simulations and analysis needs to be done in order to account for vertical motions in our model, but that is beyond the scope of this paper.

For now, it is not possible to determine the gas surface density perturbation of observations based on the continuum and vice versa, since the velocity of the spiral with respect to the disk is required (i.e., the position of the planet). Since we show in Fig. 4 that the difference in spiral morphology and amplitude changes only slightly for millimeter emission, we used the approximation that the spiral is unaffected in the continuum emission in order to test our model and make a prediction for IM Lup. If in the future spirals in both continuum and gas components are observed, this property might be useful to determine the location of the planet that launches the spiral.

6 Conclusions

In this paper, we systematically studied the variations in the properties of planet-induced spirals in gas and dust caused by the dynamics. We can conclude the following:

• The morphology and amplitude of a planet-driven spiral only changes for very large dust grains. The shape of the curve amplitude as function of the Stokes number is almost independent of planet mass, making it suitable as a tracer of dust grain size if a spiral signal in both gas and continuum is detected.

• A planet-launched spiral in surface density finds its origin in a perturbation in radial velocity, which can in turn be described as the immediate result of a perturbation in azimuthal velocity loosing or gaining angular momentum, hence pushing matter in or out.

• The azimuthal offset and amplitude change between spirals in gas and dust can be explained by the fact that dust grains feel the drag of the gas, but do not react fast enough, causing the curve to flatten.

• The above two points can be combined in a model that accurately estimates the morphology and amplitude in surface density of the spiral in dust with St ≤ 1, using only the perturbation in the azimuthal velocity of the gas.

• We find a new, convincing spiral in the dynamics of the gas in the TW Hya disk with a kink at the position of the gap in the continuum that matches well with the theoretical shape of a spiral wake. This spiral extends along both the major and minor axis of the disk, which means that we see a spiral pattern of gas that is launched from the disk in a vertical direction, probably driven by a planet’s gravitational potential.

• We prove that our model is able to handle observational data and can be used to make an estimate of how spirals will look in in-plane gas dynamics. However, data with better velocity resolution are needed to be able to compare the continuum spirals with spirals in the gas dynamics.

Acknowledgements

We thank the referee for their valuable comments. We also thank Ewine van Dishoeck for useful discussions as well as a careful reading of our manuscript and many constructive comments, Myriam Benisty for invaluable discussions and Richard Teague for an effortless providing of the reduced data of the TW Hya disk. G.R. acknowledges support from the Netherlands Organization for Scientific Research (NWO, program number 016.Veni.192.233).

Appendix A Interpretation of line density and Eq. (14)

Figure A.1 shows a toy model of an infinitesimally small part of the stream lines that particles will follow during an orbit around the central object, separated from each other by an infinitesimally small change in orbit radius, dr. The azimuthal shape of the density perturbation for these two stream lines is the same and the pitch angle of the spiral causes a small azimuthal shift between the stream line profiles, δϕ. This azimuthal shift results in a change in radial distance between the stream lines δr at certain position angles.

Fig. A.1 Toy model of an infinitesimally small part of the stream lines a particle will follow during its orbit around the central object. The median distance between the stream lines is dr, δr is the radial perturbation, δϕ is the azimuthal shift between the stream lines caused by the pitch angle of the spiral (d ϕspir ∕dr), and θ is the angle of the stream line from a circular orbit.

The result of Eq. (14) can be checked using this simple geometrical model. Defining a line density as the number of stream lines in an arbitrary radial interval, we can relate the normalized perturbation in line density (− δr/dr, shown in red in Fig. 2) to the normalized perturbation in surface density (δΣ∕Σ₀): $$\frac{\delta \Sigma }{{\Sigma }_{\text{0}}}=-\frac{\delta r}{\text{d}r}.$$(A.1)

We can find the normalized line density perturbation (δr/dr) by using the fact that the pitch angle of the spiral fixes dr/δϕ = dϕ_spir∕dr and that the angle θ, as indicated in Fig. A.1, is given by tan(θ) = δr∕δϕ. Since δr can be written in terms of the radial velocity as δr = v_rdt and δϕ can be written in terms of the angular velocity as δϕ = Ω_kdt, we can write this as tan(θ) = δr∕δϕ = v_r∕Ω_k = rv_r/v_ϕ. Dividing these two quantities gives $$\delta \Sigma =-{\Sigma }_{\text{0}}r\frac{\delta v_r}{v_{\varphi ,\text{0}}}\frac{\text{d}\varphi {(r)}_{\text{sp}}}{\text{d}r},$$(A.2)

which is the same equation as Eq. (14).

Appendix B Result of the different substeps

Fig. B.1 Results of the different substeps as discussed in Sect. 3.4. Top panel: perturbation in radial velocity as determined by FARGO-3D and determined using Eq. (9) with the azimuthal velocity perturbation as a reference. Central panel: perturbation in surface density as determined by FARGO-3D and determined from the perturbation in velocity using Eqs. (14) and (13). Bottom panel: result of the numerical integration to derive the azimuthal velocity perturbation of the dust from the azimuthal velocity of the gas, together with the FARGO-3D results and the azimuthal velocity perturbation as comparison.

Appendix C St_c as function of r

Fig. C.1 Critical Stokes number as function of radius. The changes in Stc are minor, apart from positions near the planet or boundary. The semi-analytical approach follows the numerical approach accurately, which means that the differences are mostly caused by a change in spiral width.

References

All Tables

All Figures

Fig. 1 Synthetic FARGO-3D images of spiral arms created by the potential of a super-Earth-mass planet in the gas component (panels a–c) and in the large dust grain component with St = 1 (panels d–f). The images show the perturbation of the spiral normalized to the unperturbed disk. From left to right: normalized perturbation of the surface density (a and d), the perturbation in radial velocity with respect to the Keplerian velocity (b and e), and the normalized perturbation in azimuthal velocity (c and f).

In the text

Fig. 2 Top, from left to right: azimuthal slab through the normalized perturbation in azimuthal velocity (a), radial velocity (b), and surface density (c) for Stokes numbers 10−4–1. The perturbation in the gas is indicated in black as a reference. Bottom: streamlines that a gas parcel would follow along an orbit at r = 1.7 rpl. The maximum radial deviation is marked with a black dot as a guide for the eye to the effect of the pitch angle on the flow of the particles. The perturbation in line density (see text) is indicated in red, sharing the same shape and amplitude as the perturbation in surface density (panel c).

In the text

	Fig. 3 Polar plot of the perturbation in surface density over two orbits. The lines indicate the position of the maximum of the different dust components with a zoom of the region around r = 1.7 rpl that is used to determine the change in pitch angle.
In the text

Fig. 4 Left panel: change in pitch angle with respect to the pitch angle of the spiral in the gas, fitted using the function ϕ = ar + b between r = 1.6,1.8 rpl as a function of the Stokes number. The solid lines represent the change in pitch angle relative to the pitch angle in the gas and the dotted lines show the constant azimuthal offset between the spiral in the gas and the dust in degrees. Right panel: amplitude of the spiral versus the Stokes number at r = 1.7 rpl, normalized by the amplitude of the spiral in the gas.

In the text

Fig. 5 Top: azimuthal slice through the normalized density perturbation at a radius of 1.7rpl for the dust component with St = 1. The black line shows the time integrated FARGO-3D results, the blue line shows the results directly derived from the perturbation in azimuthal velocity of the gas using the perturbed continuity equation (Eq. (13)), the red line shows the result directly derived from the perturbation in the azimuthal velocity of the gas using a direct scaling between δvr and δΣ and the pitch angle of the spiral (Eq. (14)). The dotted line illustrates the density profile of the gas component for comparison. Bottom: normalized density perturbation in polar coordinates of the result for the St = 1 component. From left to right: time integrated FARGO-3D result, the results directly derived from the perturbation in the azimuthal velocity of the gas, and the result directly derived from the perturbation in azimuthal velocity of the gas, using a direct scaling between δvr and δΣ and the pitch angle of the spiral (Eq. (14)).

In the text

Fig. 6 Map of the velocity residual as in Teague et al. (2019) for the disk around TW Hya, but with a flipped color bar to highlight the blueshifted emission (now seen in red). The dotted line serves as a guide for the eye and is determined using the theoretical shape of a spiral launched by a planet (Eq. (17)). The solid line represents the radius of the potential planet’s orbit, which coincides with the blue- and redshifted circular lobes on either side that indicate sub-Keplerian rotation, probably due to a lower opacity inside the gap, so we probe deeper layers in the disk.

In the text

Fig. 7 From left to right: observed surface density perturbation in the millimeter dust continuum of the IM Lup disk normalized to the radial symmetric surface density (a), the perturbation in radial velocity normalized to the Keplerian velocity, derived from the surface density perturbation using Eq. (14) (b), and the azimuthal velocity perturbation normalized to the Keplerian velocity derived from δvr,calc using Eq. (9) (c). The dashed lines represent a fit by eye through the maximum of the spiral in the surface density, which serves as a guide for the eye.

In the text

	Fig. 8 Spirals in gas and dust for the IM Lup disk. From left to right: estimated perturbation in azimuthal velocity from the dust continuum using Eq. (9) and Eq. (14) (a), the estimated perturbation in radial velocity from the continuum using Eq. (14) (b), and the observed velocity residuals in the optically thick CO gas emission after subtracting a Keplerian model of the data (c).
In the text

	Fig. A.1 Toy model of an infinitesimally small part of the stream lines a particle will follow during its orbit around the central object. The median distance between the stream lines is dr, δr is the radial perturbation, δϕ is the azimuthal shift between the stream lines caused by the pitch angle of the spiral (d ϕspir ∕dr), and θ is the angle of the stream line from a circular orbit.
In the text

Fig. B.1 Results of the different substeps as discussed in Sect. 3.4. Top panel: perturbation in radial velocity as determined by FARGO-3D and determined using Eq. (9) with the azimuthal velocity perturbation as a reference. Central panel: perturbation in surface density as determined by FARGO-3D and determined from the perturbation in velocity using Eqs. (14) and (13). Bottom panel: result of the numerical integration to derive the azimuthal velocity perturbation of the dust from the azimuthal velocity of the gas, together with the FARGO-3D results and the azimuthal velocity perturbation as comparison.

In the text

	Fig. C.1 Critical Stokes number as function of radius. The changes in Stc are minor, apart from positions near the planet or boundary. The semi-analytical approach follows the numerical approach accurately, which means that the differences are mostly caused by a change in spiral width.
In the text