VORTEX MIGRATION IN PROTOPLANETARY DISKS

In this section, we consider an isothermal disk with α = 3/2 and H₀ = 0.1r₀. Note that this is a disk with constant vortensity, and that, for numerical convenience, it is relatively thick for a protoplanetary disk. The dependence on H₀ is discussed in Section 7. We give the disk a circular velocity perturbation of magnitude c_s/2 over a scale H₀/2. The resulting vortensity and density perturbations after 10 orbits at r = 1 are shown in Figure 1.

Zoom In Zoom Out Reset image size

From Figure 1 we see that the vortensity perturbation remains localized, and that its radial extent roughly corresponds to the size of the initial perturbation. The density perturbation inside the vortex is modest (∼5%) compared to the vortensity perturbation (∼300%), which indicates that on scales ≪H, the disk response to the perturbation is approximately incompressible, as expected. For a vortex in equilibrium in a shear flow, there exists a relation between the vorticity perturbation and aspect ratio (Kida 1981; Lesur & Papaloizou 2009). In equilibrium, the relative perturbation in vorticity in a Keplerian disk should be (ω − ω₀)/ω₀ = −3(χ + 1)/(χ² − χ), where χ is the aspect ratio of the vortex. The strong vortex depicted in Figure 1 has an aspect ratio of ∼2.5, which would lead to an equilibrium relative vorticity perturbation of ∼−3. Noting that we can neglect density perturbations in the almost incompressible disk response, and can therefore equate vorticity and vortensity relative perturbations, this is in good agreement with the top panel of Figure 1. We note that this equilibrium is established for a range of vortex aspect ratios within a few orbits of the vortex. We therefore conclude that the vortices quickly adopt an equilibrium configuration, regardless of the initial perturbation.

The most notable density features start at a radial distance of approximately H from the vortex, where the flow relative to the vortex becomes supersonic and density waves can be excited. The resulting two-armed spiral wave pattern is reminiscent of an embedded planet (Bryden et al. 1999), with one inner wave propagating into the inner disk and one outer wave propagating into the outer disk. For embedded planets, it is well known that geometrical effects tend to favor the outer wave, which removes angular momentum from the planet and causes inward migration (Tanaka et al. 2002). An additional complication in the vortex case is that the vortex may be asymmetric with respect to r₀. We deal with vortex asymmetries in Section 5.4; here, we just note that for α = 3/2, the vortex is symmetric to a high degree. The resulting migration of the vortex is shown in Figure 2.

Zoom In Zoom Out Reset image size

Indeed the vortex migrates inward. Results are shown for three resolutions in Figure 2. The lowest resolution has 16 cells per scale height, radially, at r₀. In other words, the vortex is resolved by approximately 8 cells only, in the radial direction. In this case, the vortex weakens through numerical diffusion, which starts to slow down the migration rate after 10 orbits. Doubling the resolution pushes this time toward 20 orbits, and doubling it again gives a steady migration rate for at least 40 orbits. At early times, the migration rate is similar for all resolutions and we conclude that this migration rate is converged with respect to numerical resolution. It results in a migration timescale of $r_0/|\dot{r}| \approx 2000 \Omega _0^{-1}$ or 300 orbits at r = r₀. It is clear that for this disk, vortex migration is an important process to consider. Additional runs that included viscosity, with the kinematic viscosity ν parameterized using the α-prescription, ν = α_vc_sH, indicated that the numerical dissipation of these vortices roughly correspond to α_v = 10⁻⁴, 10⁻⁵ and 10⁻⁶, respectively, for the low-, medium-, and high-resolution cases shown in Figure 2.

Density waves propagate inward and outward away from the vortex. Provided the vortex is not too strong, we expect the waves to be in the linear regime close to the vortex, becoming nonlinear at larger distances where shocks are formed (e.g., Goodman & Rafikov 2001) and their amplitude reduces. While in the linear regime, the total rate of flow of angular momentum across a radial location that is associated with the waves, or equivalently the wave action, is conserved. When the nonlinear regime is entered and dissipation occurs, the wave action decreases. Angular momentum is advected inward in the radial direction at a rate given by

$$\begin{equation} A=2\pi r^2\langle \Sigma v_r v_{\varphi }\rangle, \end{equation} \tag{ 4 }$$

where the angle brackets denote an azimuthal average. Because there is no net mass flow associated with linear waves, the azimuthal velocity perturbation alone may be used in Equation (4) in this case. We remark that, by making use of the linearized equations governing the waves is possible to write Equation (4) in a different form. The waves have a pattern speed Ω_vort which is also the angular velocity of the vortex. Thus for the perturbations associated with them, we have

$$\begin{equation} \Omega _\mathrm{vort} \frac{\partial }{\partial \varphi }=\frac{\partial }{\partial t}. \end{equation} \tag{ 5 }$$

Accordingly, the perturbed azimuthal component of the equation of motion gives

$$\begin{equation} (\Omega -\Omega _\mathrm{vort})\frac{\partial v_{\varphi }}{\partial \varphi }+\frac{\kappa ^2}{2\Omega }v_r = -\frac{c_s^2}{r\Sigma }\frac{\partial \Sigma ^{\prime }}{\partial \varphi }, \end{equation} \tag{ 6 }$$

where Σ' is the surface density perturbation. Noting that to linear order, the Lagrangian displacement ξ_r satisfies

$$\begin{equation} v_r =\left(\Omega - \Omega _\mathrm{vort}\right) \frac{\partial \xi _r }{\partial \varphi }, \end{equation} \tag{ 7 }$$

we see that

$$\begin{equation} \left(\Omega -\Omega _\mathrm{vort} \right)\left(v_{\varphi }+\frac{\kappa ^2}{2\Omega }\xi _r\right) = -\frac{c_s^2\Sigma ^{\prime }}{r\Sigma }. \end{equation} \tag{ 8 }$$

Using Equation (8) to substitute for v_φ in Equation (4) and making use of Equation (7), it is readily seen that we can rewrite Equation (4) in the form

$$\begin{equation} A= 2\pi \left\langle \frac{r c_s^2 \Sigma ^{\prime }v_r}{\Omega _\mathrm{vort}-\Omega }\right\rangle. \end{equation} \tag{ 9 }$$

We can now verify that this migration is due to the density waves by studying the behavior of wave action. The change in the angular momentum content of the fluid on account of the presence of the vortex is

$$\begin{equation} L_\mathrm{vort}=\int (\Sigma r^2\Omega -\Sigma _i r^2\Omega _i) {\it dS}, \end{equation} \tag{ 10 }$$

where subscripts i denote the initial conditions, and the integral is taken over the whole vortex surface area. We expect dL_vort/dt = ΔA, with ΔA the difference in wave action defined above evaluated for the outer and the inner wave.

We show the wave action for three resolutions in Figure 3. Equation (9) only has meaning in the wave region, at distances of at least 2H/3 from the vortex. Although our lowest resolution run clearly does not resolve the waves well enough to get the correct wave action on both sides, ΔA is very similar to the higher resolution runs. In Section 7, we make sure our resolution is at least as high as that of the dashed line in Figure 3. The similarity of the wave action for the highest resolution runs translates into a similar migration rate after 10 orbits (see Figure 2). We have measured the change in angular momentum of the vortex between 5 and 10 orbits and compared this to what is expected from the wave action. The results are summarized in Table 1.

Zoom In Zoom Out Reset image size

We point out that there is a large uncertainty in the measured ΔL. This is because L_vort can only be measured to an accuracy of ∼5% because of ambiguities in defining the vortex as a separate entity from the background flow, and ΔL results from a subtraction of two almost equal numbers (typically ΔL ≈ 0.1L_vort). We therefore estimate that there is an uncertainty of at least 50% in the measured values of ΔL, which makes the agreement between ΔAΔt and ΔL reasonable (see Table 1). Therefore, the inward migration of the vortex is consistent with the wave torque.

We show the time evolution of the wave action difference, or wave torque, in Figure 4, for the same three resolutions as in Figure 2. In all cases, a steady torque is set up after 3 orbits, at which time the vortex has reached its equilibrium with the background shear. After the equilibrium has been set up, the slow time evolution of the torque is due to numerical diffusion (the decline in the low resolution case) and slowly varying background disk properties due to the changing radial location of the vortex (the slow increase in the highest resolution cases).

Zoom In Zoom Out Reset image size

In order to understand the origin of the waves emitted by accretion disk vortices, let us introduce the local shearing box model (Hawley et al. 1995). In this model, we neglect curvature effects in Equations (1)–(2) and introduce a local Cartesian frame corotating with the disk at r₀ by defining x = r − r₀ and y = r₀φ. In this model, Equations (1) and (2) read

$$\begin{eqnarray} &&\frac{\partial \Sigma }{\partial t}+\bm {\nabla \cdot } \Sigma \bm {v}=0, \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \frac{\partial \bm {v}}{\partial t}+\bm {v\cdot \nabla v}&=& -c_s^2\bm {\nabla } \ln \Sigma -\bm {\nabla } (-q\Omega ^2x^2)-2\bm {\Omega \times v},\nonumber\\ \end{eqnarray} \tag{ 12 }$$

where Ω = Ω(r₀) is the Keplerian frequency at r₀, q = 3/2 for a Keplerian rotation profile and we have assumed an isothermal gas. In the following, we consider a vortex located at r = r₀ and perturbing the surrounding flow. The vortex core is a nonlinear solution to the above equations (see, e.g., Lesur & Papaloizou 2009), and we consider the linear perturbations produced by the vortex core at long distance. Since the vortex is, in first approximation, a quasi-steady structure, we consider stationary perturbations of the Keplerian flow v₀ = −qΩe_y introducing

$$\begin{eqnarray} &&\Sigma ^{\prime }(x,y)=\Sigma (x,y)-\Sigma _0, \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&\bm {u}(x,y)=\bm {v}-\bm {v}_0.\qquad \end{eqnarray} \tag{ 14 }$$

As a further simplification, we only consider the case of a constant background density profile Σ₀. To reduce the problem to a system of ordinary differential equations, we introduce a Fourier decomposition in the y direction:

$$\begin{eqnarray} &&\Sigma ^{\prime }=\tilde{\Sigma }^{\prime }(x)\exp ({\it iky}), \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} &&\bm {u}=\tilde{\bm {u}}(x)\exp ({\it iky}).\ \end{eqnarray} \tag{ 16 }$$

Using this decomposition in Equations (11) and (12) yields

$$\begin{eqnarray} &&-i\sigma (x)\tilde{\Sigma }^{\prime }+\Sigma _0\left (\frac{d \tilde{u}_x}{dx}\,+\,{\it ik}\tilde{u}_y\right)=0 \ \ \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} &&-i\sigma (x)\tilde{u}_x-2\Omega \tilde{u}_y+\frac{c_s^2}{\Sigma _0}\frac{d\tilde{\Sigma }^{\prime }}{dx}=0 \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} &&-i\sigma (x)\tilde{u}_y+(2-q)\Omega \tilde{u}_y+ic_s^2k\frac{\tilde{\Sigma }^{\prime }}{\Sigma _0}=0,\quad\ \ \ \end{eqnarray} \tag{ 19 }$$

where we have defined σ(x) = qΩx. After some algebra, one gets the following equation for the azimuthal velocity $\tilde{u}_y$:

$$\begin{equation} \frac{d^2\tilde{u}_y}{dx^2}+\Lambda ^2(x)\tilde{u}_y=0, \end{equation} \tag{ 20 }$$

where

$$\begin{equation} \Lambda ^2(x)=\frac{\sigma ^2(x)-2(2-q)\Omega ^2}{c_s^2}-k^2. \end{equation} \tag{ 21 }$$

This equation describes the propagation of a linear acoustic-inertial wave in the shearing box. As stated above, we consider a nonlinear vortex core located at x = 0 emitting a wave at large distance. Assuming that the radial extent of the vortex is very small compared to the disk thickness H = c_s/Ω, although this need not be the case for the azimuthal extent, the vortex core can be seen in first approximation as a boundary condition at x = 0 for the wave emission problem. Without loss of generality, we will consider the case of a wave emitted in the x > 0 region, the x < 0 solutions being deduced by symmetry arguments.

It is worth noting that Equation (20) is in fact the vorticity conservation equation in a compressible medium, with a constant background vorticity (2 − q)Ω. In a global disk however, background vorticity is not constant and density might also depend on radius. In this context, one would write a conservation equation for the vortensity ∇ × v/Σ which would be similar to Equation (20) if the background vortensity was constant. Therefore, the wave emission model we are discussing here should be understood as a local model for a constant vortensity disk.

The above equation admits solutions in the form of parabolic cylinder functions. However, to understand the physical origin of vortex waves, we present here Wentzel–Kramers–Brillouin (WKB) solutions. We first note that Equation (20) has a turning point at

$$\begin{equation} \nonumber x_s=\frac{\omega _{ai}}{q\Omega k}, \end{equation} \tag{ 22 }$$

where we have introduced $\omega _{ai}=\sqrt{k^2c_s^2+q(2-q)\Omega ^2}$, the natural frequency of acoustic-inertial waves. This particular location corresponds to the sonic line of the vortex. This is the line along which the flow speed is equal to the wave speed as given by the phase velocity. In the limit of large k this phase velocity becomes equal to the sound speed hence the term sonic line. For x < x_s ("subsonic" region), the wave has an exponential shape and it corresponds to the logarithmic tail of an incompressible vortex. Using first-order WKB theory in the limit x ≪ x_s, we get

$$\begin{eqnarray} \tilde{u}_y&=&A\exp \left (-\int _x^{x_s} [-\Lambda ^2(u)]^{1/2}\,du\right)\nonumber \\ &&+B\exp \left (\int _x^{x_s} [-\Lambda ^2(u)]^{1/2}\,du\right). \end{eqnarray} \tag{ 23 }$$

On the other hand, for x > x_s ("supersonic" region), the wave is essentially sinusoidal and can travel for very large distances. In the limit x ≫ x_s, we obtain

$$\begin{eqnarray} \tilde{u}_y&=& C\exp \left (i \int _x^{x_s} [\Lambda ^2(u)]^{1/2}\,du\right)\nonumber \\ &&+D\exp \left(-i \int _x^{x_s} [\Lambda ^2(u)]^{1/2}\,du\right), \end{eqnarray} \tag{ 24 }$$

where it can be shown that an outgoing wave (wave transporting angular momentum outward) is obtained when C = 0. This solution leads to the long "wakes" one observes in compressible simulations of vortices. In particular, in the limit x ≫ kc_s/qΩ, we find that the wavefronts of the supersonic region are given by

$$\begin{equation} y=y_0-\frac{q\Omega }{2 c_s}x^2, \end{equation} \tag{ 25 }$$

which produce a pattern similar to the wakes found in simulations (Figure 5).

Zoom In Zoom Out Reset image size

The connection of the supersonic and subsonic regions allows one to determine C and D knowing A and B. This can be done using asymptotic matching techniques based on Airy functions. Although the precise derivation of these solutions is beyond the scope of this paper, we find that |C| ∼ |D| ∼ max(|A|, |B|), i.e., the amplitude of the solution does not change significantly as one crosses the sonic line. However, it should be noted that the wave amplitude in the supersonic region is significantly reduced compared to its amplitude at x = 0 due to the exponential decay in the subsonic region. One can estimate this attenuation factor Γ with

$$\begin{eqnarray} \Gamma ^\pm =\frac{|\tilde{u}_y(\pm x_s)|}{|\tilde{u}_y(x_0)|}&\simeq &\exp \left(\mp \int _0^{\pm x_s} [-\Lambda ^2(u)]^{1/2}\,du\right)\nonumber \\ &=&\exp \left(-\frac{\pi \omega _{ai}^2}{4q\Omega kc_s}\right). \end{eqnarray} \tag{ 26 }$$

This attenuation factor has a simple physical interpretation: it relates the vortex amplitude to the emitted density wave amplitude at a particular wavenumber. One can show that the attenuation factor is minimum for $k=\sqrt{q(2-q)} H^{-1}$ where H is the disk thickness H ≡ c_s/Ω.

It is clear that in the shearing box, the wave emitted in the inner side and the wave emitted in the outer side have the same amplitude (Γ⁺ = Γ⁻). Because of this symmetry, one does not expect any migration. However, if one includes global effects such as curvature and the dependence of Ω as a function of radius, then an asymmetries might occur, resulting in vortex migration.

To extend our analysis to cylindrical geometry, we consider the case of a constant vortensity disk Σ ∝ r^−3/2 with a constant aspect ratio c_s ∝ r^−1/2. Furthermore, we assume the flow is barotropic and work in an inertial frame. We Fourier transform the velocity field in t and φ, defining $\bm {u}(t,\varphi,r)= \sum _m \int d\omega \exp [ i(\omega t-m\varphi)] \tilde{\bm {u}}(\omega,m,r)$. For a vortex corotating with the gas at r = r₀, we adopt ω = mΩ(r₀). We then perform the corresponding linearization of the equations to that described in the shearing box case. One eventually obtains

$$\begin{equation} \frac{1}{\Omega r}\frac{d}{dr}\Omega r \frac{d}{dr} r\tilde{u}_\varphi +\Lambda ^2(r) r\tilde{u}_\varphi =0, \end{equation} \tag{ 27 }$$

where

$$\begin{equation} \Lambda ^2(r)=\frac{\sigma ^2(r)-\Omega ^2}{c_s^2}-\frac{m^2}{r^2} \end{equation} \tag{ 28 }$$

and σ(r) = ω − mΩ(r). This equation is very similar to the shearing-sheet problem. As in the shearing-sheet case, we find two turning points r^±_s on both sides of the corotation radius r₀ defined as σ(r₀) = 0. However, it should be noted that |r⁺_s − r₀| < |r⁻_s − r₀|, i.e., the outer sonic line is closer to corotation than the inner sonic line. Indeed, imposing Λ(r^±_s) = 0 leads to

$$\begin{equation} r_s^{\pm }=r_0\left (1\pm \sqrt{\frac{1}{m^2}+\epsilon ^2}\right)^{2/3}, \end{equation} \tag{ 29 }$$

where is the constant aspect ratio = H/r. In the limit m → ∞ we have

$$\begin{equation} r_s^{\pm }=r_0\left(1\pm \frac{c_s}{r\Omega }\right)^{2/3}. \end{equation} \tag{ 30 }$$

Hence in this limit the sonic lines occur where the velocity relative to the vortex, that is associated with the background shear, is equal to the sound velocity. This corresponds to the large k limit in the local model. However, this does not occur for small m. In fact Equation (29) indicates that the distance of the sonic lines from the vortex increases as m decreases to the extent that when m = 1, the inner sonic line is out of the domain. This is because the phase velocity of the density waves increases as m decreases. WKB solutions very similar to the one presented for the shearing box case can be derived and one can compute the attenuation factor as

$$\begin{equation} \Gamma ^\pm =\frac{r_s^\pm |u_\varphi \big(r_s^\pm\big)|}{r_0|u_\varphi (r_0)|}\simeq \exp \left (\mp \int _{r_0}^{r_s^\pm } [-\Lambda ^2(u)]^{1/2}\,du\right). \end{equation} \tag{ 31 }$$

As it can be seen from this expression, the attenuation factor on the inner side (Γ⁻) is expected to be different from the attenuation factor on the outer side Γ⁺ because Λ is not symmetric on both sides of corotation. One can easily check that |Λ²(r₀ + δr)| < |Λ²(r₀ − δr)| when δr < r⁺_s − r₀. Moreover, since |r⁺_s − r₀| < |r⁻_s − r₀|, we conclude that Γ⁺ > Γ⁻, i.e., the outer wave always have an amplitude larger than the inner wave.

This can also be seen by making the transformation from r to ψ such that

$$\begin{equation} \frac{r}{r_0}=1+\sqrt{\epsilon ^2+\frac{1}{m^2}}\cos \psi, \end{equation} \tag{ 32 }$$

so that the domain (r⁺_s, r₀) is mapped to (0, π/2) and the domain (r⁻_s, r₀) is mapped to (π/2, π). The integrals

$$\begin{equation} \ln \Gamma ^\pm = \mp \int _{r_0}^{r_s^\pm } [-\Lambda ^2(u)]^{1/2}\,du \end{equation} \tag{ 33 }$$

transform to

$$\begin{equation} \ln \Gamma ^\pm = -\int _{-\pi (-1\pm 1)/4}^{-\pi (-3\pm 1)/4} \left(\frac{2m}{3\epsilon }\right) \frac{\big(\epsilon ^2+\frac{1}{m^2}\big)\sin ^2\psi }{ 1+ \sqrt{\epsilon ^2+\frac{1}{m^2}}\cos \psi }d\psi. \end{equation} \tag{ 34 }$$

On account of the differing integration ranges for the outward and inward propagating waves, the wave attenuation factor is readily seen to be larger for the inward propagating wave. Although the calculation can be done in full analytically, it is adequate for our purposes to regard and 1 m⁻¹ to be comparable small parameters and express the result correct up to order . Then

$$\begin{equation} \ln \Gamma ^\pm = - \frac{\pi }{6} \left(m\epsilon + \frac{1}{ m\epsilon }\right)\left(1\mp \frac{4}{3\pi }\sqrt{\epsilon ^2+\frac{1}{m^2}}\right). \end{equation} \tag{ 35 }$$

Thus the expected ratio of outward to inward wave flux arising from this asymmetry is

$$\begin{equation} \left(\frac{ \Gamma ^+}{\Gamma ^-}\right)^2 =\exp \left[- \frac{8 }{9}\sqrt{\frac{\epsilon }{m}} \left(m\epsilon + \frac{1}{ m\epsilon }\right)^{3/2}\right]. \end{equation} \tag{ 36 }$$

As we will see later, this explains the inward migration found in simulations with a similar profile. For example when = 0.1 and m = 5, Equation (36) predicts an asymmetry comparable to that seen in Figure 3. More generally, when m = 1/ we obtain a ratio ∼1–2.5 m⁻¹.

The above discussion on the local shearing box model assumed a constant vortensity and in this case the response to a symmetric vortex was found to be symmetric leading to equal inward and outward angular momentum flow rates and hence no net migration. However, if the constraint of constant vortensity is relaxed, an asymmetry in the response occurs in the subsonic region of the flow even when the vortex is symmetric. To see this, we consider a local shearing box model in which the background density has the form Σ = const × exp(−αx) with α being a constant. The vortensity of course varies as the inverse of Σ.

As we are interested in the subsonic region, we assume that the linear response satisfies the anelastic condition $\nabla \cdot (\Sigma \tilde{{\bf u}}) =0.$ After eliminating $\tilde{\Sigma }^{\prime }$ from Equations (17) and (18) and then using the anelastic condition to eliminate $\tilde{u_y},$ we obtain an equation for $Q=\tilde{v}_x\sqrt{\Sigma }$ in the form

$$\begin{equation} \frac{d^2 Q}{dx^2}= Q\left(k^2-\frac{\alpha }{3x}\right), \end{equation} \tag{ 37 }$$

where k² = m² + α²/4. We remark that for assumed unit radius of the center of the box, k_y ≡ m is the azimuthal mode number and α is the power law index occurring in Σ ∝ r^−α for a global disk.

We are interested in solutions of Equation (37) that decay exponentially for large |x| and are matched as x → 0. However, we do dot require that derivatives match for x → 0. This is because we assume that there is a nonlinear vortex core there that is not modeled by Equation (37) which only describes the regions between the vortex core and the sonic lines.

Without loss of generality, we shall consider x > 0. The case x < 0 can be recovered by letting α → −α. An appropriate integral representation for the required solution follows from representations of "Coulomb Wave Functions" given in Abramowitz & Stegun (1972) and may be written as

$$\begin{eqnarray} Q &= &C_0\Gamma (q)(2k)^{-q}\exp (-kx) \nonumber \\ &&\times \int ^{\infty }_0 \left(x+\frac{t}{2k}\right)^{-q}t^{q} \exp (-t)dt, \end{eqnarray} \tag{ 38 }$$

where C₀ is an arbitrary constant and q = −α/(6k). The above form may be used to connect the form of the solution near x = 0 to that for large x giving

$$\begin{eqnarray} Q(x=0)&=& 1 \rightarrow Q \nonumber \\ &\sim & \Gamma (q+1)(2kx)^{-q}\exp (-kx) \ \ \ \ {\rm for} \ \ \ \ x \rightarrow \infty\nonumber\\ \end{eqnarray} \tag{ 39 }$$

which may also be expressed as

$$\begin{eqnarray} y(x=0)&=&1 \rightarrow \nonumber \\ y &\sim & \exp (-kx -q\ln (2kx) -\gamma q) \ \ \ \ {\rm for} \ \ \ \ x \rightarrow \infty.\nonumber\\ \end{eqnarray} \tag{ 40 }$$

Here, assuming that k is large, we have used Γ(q + 1) ∼ −γq, where γ is Euler's constant ∼0.577.

Recalling that q = −α/(6k), we find that the solution exceeds the pure exponential for large x > 0 when α is positive. That means the outgoing wave is favored in this case. Conversely, the ingoing wave is favored when α < 0 indicating the potential for inward migration in this case. But, noting that kx ∼ 1, this is found to lead only to a relatively weak modification of the wave propagation, except when |α| is large and the vortensity gradient is very strong.

However, it should be noted that the presence of a vortensity gradient leads to the modification of the vortex structure itself. In the above calculation, we have assumed that the vortex core was symmetrical, i.e., u_x(x = +) = u_x(x = −) when → 0. However, in the presence of a vortensity gradient, vortex solutions are not symmetric anymore. To show this effect, we have solved a set of anelastic equation in a shearing box including a radial vortensity gradient:

$$\begin{eqnarray} \bm {\nabla \cdot }(\Sigma _0(x)\bm {v})&=&0, \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} \frac{\partial \bm {v}}{\partial t}+\bm {v\cdot \nabla v}&=&-\bm {\nabla } \Pi -\bm {\nabla } (-q\Omega ^2x^2)-2\bm {\Omega \times v}. \end{eqnarray} \tag{ 42 }$$

In these equations, the vortensity gradient is imposed by the density profile Σ₀(x) and density/sound waves are filtered out. We have integrated numerically these equations assuming an exponential profile for the surface density Σ = Σ₀exp(−αx) using a spectral scheme. The initial conditions were chosen to be those appropriate to a Kida vortex (Kida 1981) of aspect ratio 8. This vortex then relaxes to a new quasi-stationary state corresponding to a vortex solution with a vortensity gradient. We show in Figure 6 the resulting vorticity distribution in the case α = 4 (the length unit being the vortex radial size). As it can be seen from this figure, the vortex core becomes asymmetric when a vortensity gradient is imposed. In particular, the perturbation extends for a greater distance on the less dense side and the vortex is "flattened" on the denser side. As a result, the less dense part of the vortex, being more extended, will excite stronger density waves than the inner side, resulting in a new source of wave asymmetry. It should be noted that this effect is purely hydrodynamic and is due to the fact that the perturber (the vortex) "feels" and responds to the surrounding flow. Therefore, no corresponding effect exists in the case of a planet embedded in a disk.

Zoom In Zoom Out Reset image size

We now summarize some important differences between vortex migration and Type I planet migration (see Table 2). Although both are driven by asymmetric density waves, the excitation mechanism as well as the sources of the asymmetry are fundamentally different. First of all, it should be pointed out that although it has very similar properties, in contrast to the planet migration problem, the asymmetry discussed in Section 5.3 does not involve any torque calculation. Instead, it appears as an attenuation in the linear propagation of density waves emitted by the vortex and it leads to a stronger outer wave. This is a major difference with the planet case in which asymmetric Lindblad torques acting at the sonic line are believed to drive planet migration. Furthermore, since we find the vortex always to move with the local gas velocity, there is no effect of the pressure buffer as discussed in the Appendix. We find that wave asymmetries only weakly depend on the background density profile, favoring the outer wave and therefore inward migration.

However, unlike a planet, a vortex directly feels the disk and its associated density profile, which leads to an asymmetric vortex core as discussed in Section 5.4 in the case of a vortensity gradient. The resulting asymmetric wave emission can either counteract or reinforce inward migration due to geometrical effects as discussed above.

We also point out that since there is no gravitational interaction between disk and vortex, there are no horseshoe orbits close to corotation. In the planet case, the horseshoe drag (Ward 1991) is an additional source of angular momentum exchange (Paardekooper & Papaloizou 2009), that can be strong enough to counteract the wave torque in non-isothermal disks (Paardekooper et al. 2010). In the vortex case, we only have the wave torque.

We have measured the migration speed in isothermal disks with H₀ = 0.1r₀ for various background density gradients; the results are depicted in Figure 7. The migration speed is approximately proportional to −α, with outward migration happening for surface density profiles that increase outward. As indicated above a constant surface density disk shows almost no vortex migration. For such a disk, the geometrical effects favoring inward migration (see Section 5) are almost exactly balanced by the asymmetry of the vortex, which favors outward migration. The vortex asymmetries make vortex migration more strongly dependent on α than is the case for linear Type I planetary migration (Tanaka et al. 2002) that requires α ≲ −2.5 for migration reversal.

The asymmetry in the vortex can be seen, for example, in the Fourier components of the radial velocity. In Figure 8, we show the m = 10 component for a constant vortensity disk (black curves) and a constant surface density disk (gray lines). The approximate size of the vortex is indicated by the vertical dotted line. The constant vortensity disk gives rise to a vortex that is almost symmetric for small |x|, as expected. The constant surface density disk shows a stronger response in the inner disk (dashed curve), which favors outward migration and works against the geometrical effects of Section 5. From Figure 7 we see that for this disk, these competing effects almost cancel each other, resulting in almost no migration for α = 0. We comment that for a disk with a temperature gradient (through a spatially varying, but constant in time, sound speed), this neutral state corresponds to a disk with constant pressure.

Zoom In Zoom Out Reset image size

Until now, we have focused on a disk with H₀ = 0.1r₀. Since we need to be able to resolve the vortex, which is of size H₀/2 typically, it is numerically convenient to work with relatively thick disks. However, we expect the angular momentum transport, and therefore the migration speed of the vortex, to depend on H₀. From Equation (9), we expect vortices of size ∼H₀ and velocity perturbations ∼c_s, that A ∝ H³₀. Since the most important contributions to the angular momentum flux come from m ∼ r₀/H₀, we see from Equation (36) that the ratio of wave actions outside the orbital radius of the vortex and inside the orbital radius of the vortex is A_o/A_i ≈ 1 − CH₀/r₀, with C being a numerical constant of order unity. This means that the action difference ΔA = A_i − A_o ∝ H⁴₀. This equals the rate of angular momentum transfer between vortex and disk and therefore governs vortex migration.

In the top panel of Figure 9, we show the measured wave action inside the orbit of the vortex (circles), outside the orbit of the vortex (triangles), and their difference for different values of H₀. The size of the initial perturbation is a fixed fraction (1/2) of H₀, and the magnitude of the velocity perturbation is a fixed fraction of the sound speed (3/4). With these initial conditions, we expect above scalings to hold approximately. We have overplotted fits through the numerical results ∝H³₀ (for A_i and A_o) and ∝H⁴₀ (for ΔA), showing that indeed these scalings hold for a wide range of H₀, even in cases where the disk can no longer be called "thin."

Zoom In Zoom Out Reset image size

The angular momentum of the vortex scales as H²₀, according to Equation (46). Combined with an angular momentum transfer rate ∝H⁴₀ as derived above, we derive a migration timescale ∝H⁻²₀, or a radial velocity of the vortex ∝H²₀. This is the relation we explore in the bottom panel of Figure 9, where we show the magnitude of the migration velocity of the vortex for different values of H₀. The solid curve is a fit through the numerical data ∝H²₀, showing that again the scaling derived shows good agreement with the numerical simulations. Note that unlike for the case of planet migration, vortex migration is faster in thicker disks. This is because the size of the vortex scale with H₀, so that for all H₀ we get a similar perturbation. This would amount to considering planets of higher mass in thicker disks.

Even though vortex migration slows down for thinner disks, it is still very substantial at H₀ = 0.05r₀, a value that is though to be appropriate for protoplanetary disks. From the measured radial velocity, we derive a migration timescale of r₀/|v_r| = 1700 orbits. This is very short, but we point out that these results were obtained for a strong vortex with a size comparable to H₀. Below, we discuss how migration depends on the size of the vortex.

Vortex migration is driven by perturbations at the sonic line, located roughly at a distance 2H₀/3 away from the center of the vortex. Away from the core, perturbations due to the vortex decay exponentially in the subsonic region (see Equation (26)). We therefore expect smaller vortices to show less migration. Note that a smaller radial vortex size s can be compensated for in principle, as far as migration is concerned, by a stronger initial velocity perturbation. Unfortunately, the regime s ≪ H₀ is difficult to study numerically because of computational constraints. We therefore restrict ourselves to s ⩾ H₀/10 to study the effect of vortex size on migration.

We show numerical results for H₀ = 0.1r₀ in Figure 10. Note that the radial size s quoted in the figure is actually the size of the initial perturbation rather than the actual vortex size. We have always found the two to agree well. We see that indeed smaller vortices migrate slower, with a factor 50 between a vortex of size H₀/10 and a vortex of size H₀/2. For s > H₀/2, migration is roughly constant, but for the larger cases the vortex extends beyond the sonic line, a situation that is difficult to maintain.

Zoom In Zoom Out Reset image size

It is difficult to provide simple scalings in this case, since for most of our results the vortex core occupies a substantial fraction of H₀. We have checked that similar results are obtained for different values of H₀. Therefore, even a vortex of size H₀/10 with H₀ = 0.05r₀ would have a migration timescale of ∼30, 000 orbits only.

As mentioned above, a stronger velocity perturbation will induce stronger wave emission and therefore faster migration. For a given radial vortex size, a stronger velocity perturbation will lead to a vortex that is more circular, while weaker vortices have a larger aspect ratio χ. To study the influence of the strength of the vortex, we have varied the initial velocity perturbation at a fixed radial size of the vortex. The aspect ratio of the resulting vortex was measured by considering contours of constant vortensity.

The resulting wave action and migration speeds are displayed in Figure 11. As expected, stronger vortices migrate faster. This is not only due to a larger perturbation amplitude in the velocity but also due to a change in the perturbation spectrum. Since the attenuation factor peaks at m ∼ r₀/H₀, the perturbation amplitudes around this value of m play a major role in determining the migration speed. A vortex of aspect ratio unity and radial size H₀ gives rise to large amplitudes at m ∼ r₀/H₀ and therefore to strong migration. Increasing the aspect ratio shifts the peak of the vortex spectrum to lower values of m, which are less effective in inducing angular momentum transfer. This is what is observed in Figure 11. Migration slows down by an order of magnitude when the vortex weakens from χ = 2 to χ = 3.5.

Zoom In Zoom Out Reset image size

The strong dependence of migration rates on the background surface density profile opens up the possibility of halting vortex migration at special locations in the disk, where the local gradient of surface density is small. This is explored using a density profile as depicted in the top panel of Figure 12. The corresponding profile of α = −dlog Σ/dlog r is shown in the bottom panel of Figure 12. From the discussion above we expect a vortex to migrate inward outside r = r₀, where α>0, and outward inside r = r₀, where α < 0. In this configuration, r = r₀ is an equilibrium radius for migrating vortices.

Zoom In Zoom Out Reset image size

We have put in a strong vortex, with initial velocity perturbation amplitude of 3c_s/4 over a scale H₀, at r/r₀ = 1.5 in a disk with H₀ = 0.1r₀. With such a strong vortex, we make sure it does not dissipate before reaching the equilibrium radius. In Figure 13, we show the time evolution of the orbital radius of the vortex (bottom panel), the strength of the vortex S_v (defined by the minimum in the relative vortensity perturbation, top panel), and the local value of α (middle panel). From the top panel we see that the vortex gets weaker with respect to the background flow, partly because of numerical diffusion, partly because of the migration itself. Since the vortex, being a negative perturbation in vortensity, is moving into a region of lower vortensity, the relative perturbation is expected to go down, as observed in Figure 13. Results for two different resolutions are shown in Figure 13; the evolution of S_v that appears similar in both cases is due to migration, while differences are most likely caused by numerical diffusion.

Zoom In Zoom Out Reset image size

From the bottom panel of Figure 13, we see that migration starts off fast, slows down over time, and almost stalls at r = r₀. This slowing down of migration is completely due to the decrease in α as depicted in the middle panel of Figure 13. At r/r₀ = 1.5, α ≈ 4, giving rise to fast inward migration. As the vortex approaches the surface density maximum, α approaches zero and migration slows down. Since we obtain essentially the same result for both resolutions, we conclude that the halting of inward migration is due to r = r₀ being an equilibrium radius rather than due to the vortex dissipating. This therefore provides an example of vortex trapping at a surface density bump.