FAST MODES AND DUSTY HORSESHOES IN TRANSITIONAL DISKS

Transitional protoplanetary disks possess inner cavities nearly devoid of dust (e.g., Espaillat et al. 2014). Outside these cavities, in some of the brightest disks, radio images reveal that dust is not axisymmetric, but clumps strongly to one side. At dust continuum wavelengths, surface brightness contrasts range from values on the order of a few (Brown et al. 2009; Tang et al. 2012; Isella et al. 2013; Pérez et al. 2014) to ∼30 (HD 142527; Ohashi 2008; Casassus et al. 2013; Fukagawa et al. 2013) to ∼130 (Oph IRS 48; van der Marel et al. 2013; Bruderer et al. 2014). These contrasts are lower bounds where observations are unresolved. Unlike dust, overdensities in CO gas are limited to factors of a few at most (van der Marel et al. 2013; Bruderer et al. 2014; Pérez et al. 2014b).

The lopsided dust concentrations—which may be giving us a first empirical glimpse into the process of planetesimal formation—cry out to be well understood. One proposed mechanism is dust trapping by gas vortices (e.g., Lyra & Lin 2013; Zhu & Stone 2014). The inner rims of transitional disks may be subject to the Rossby wave instability, which can spawn anticyclonic vortices (e.g., Lovelace & Romanova 2014). Whether vortices are radially and azimuthally large enough to match the imaging data is debatable.

We submit here a different explanation, one that follows from a simple fact: a lopsided disk moves the barycenter away from the host star. The gravitational potential therefore includes an "indirect" term (e.g., Murray & Dermott 1999). Our Letter solves for gas streamlines including the indirect potential, in the limit that gas does not feel its own gravity. We will discover in this limit that gas occupies a radially wide, horseshoe-shaped annulus that readily concentrates dust grains into patterns similar to those observed. Concentration is by gas drag, not dust self-gravity.

Our proposal is that the gas disk exhibits a global "fast" mode of azimuthal wavenumber m = 1, whose lopsidedness shifts the barycenter off the star. "Fast" means that the disk's pattern speed equals its mean Keplerian motion, not the apsidal precession frequency that characterizes m = 1 "slow" modes (for slow modes, see Tremaine 2001). Fast modes are more attractive than slow modes for shaping dust grain trajectories and explaining the radio observations, since the non-inertial forces arising from a slow mode's pattern speed are too small to compete with the substantial gas drag forces felt by dust.

We compute gas streamlines (test particle orbits) analytically in Section 2 and numerically in Section 3. How dust concentrates aerodynamically is described in Section 4. Limitations of our calculations—notably our neglect of gas self-gravity, which prevents us from computing fast modes self-consistently—are discussed in Section 5.

We consider analytically the dynamics of a test particle in the gravitational potential of a star with a prescribed motion. The star has mass M_* and executes a circular orbit of radius μ with fixed angular frequency Ω about the origin (star–disk barycenter). This orbit represents the star's "reflex motion" in response to the m = 1 component of the disk's potential. We work in a frame rotating with the star and centered on the origin (Figure 1). In our unit system, the gravitational constant times the stellar mass GM_* = 1 − μ, Ω = 1, and μ ≪ 1. Our problem is identical to the standard restricted three-body problem, except that the test particle does not feel the gravity of the secondary (read: disk). Because we assume the star's motion arises from the disk potential but neglect the latter when computing the test particle's motion, our calculations are not self-consistent; we proceed anyway in the hope that some of the grosser, qualitative features of our model will survive a more careful study.

Zoom In Zoom Out Reset image size

In Cartesian coordinates, the test particle obeys

$$\begin{eqnarray} \ddot{x} - 2\Omega \dot{y} & = & \frac{\partial U}{\partial x} \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \ddot{y} + 2\Omega \dot{x} & = & \frac{\partial U}{\partial y}\protect, \protect \end{eqnarray} \tag{ 2 }$$

where U is the celestial mechanician's centrifugal potential plus the potential of the star (offset from the origin by μ on the negative x-axis):

$$\begin{equation} U = \frac{\Omega ^2}{2} (x^2 + y^2) + \frac{1-\mu }{\sqrt{(x+\mu)^2 + y^2}}. \end{equation} \tag{ 3 }$$

Equilibrium points, where gravity and centrifugal forces balance, lie along the x-axis. Setting $\ddot{x} =\ddot{y} =\dot{x} = \dot{y} = y = 0$ in Equations (1) and (2), we find

$$\begin{equation} \frac{\partial U}{\partial x} = x_{\rm eq} - \frac{(x_{\rm eq}+\mu)(1-\mu)}{|x_{\rm eq}+\mu |^3} = \protect 0. \protect \protect \protect \protect \protect \end{equation} \tag{ 4 }$$

To order μ, the two equilibrium points are located at

$$\begin{equation} x_{\rm eq} = 1 - \mu \,\,\, {\rm and} \,\, -1 - \mu /\protect 3. \protect \protect \protect \protect \protect \end{equation} \tag{ 5 }$$

Assuming that small displacements about these equilibria evolve with time t according to exp (λt), we linearize and combine Equations (1) and (2) to arrive at the relation for the eigenfrequency λ:

$$\begin{eqnarray} \lambda ^4 & + & \left[ 2 - \frac{1-\mu }{(x_{\rm eq} + \mu)^3} \right]\lambda ^2 \nonumber \\ & + & \left[ 1 + \frac{1-\mu }{(x_{\rm eq} + \mu)^3} -2 \frac{(1-\mu)^2}{(x_{\rm eq} + \mu)^6} \right] =\protect 0. \protect \protect \protect \protect \protect \end{eqnarray} \tag{ 6 }$$

The equilibrium point at x_eq < 0 is unstable (Re λ > 0). By contrast, the equilibrium point at x_eq = 1 − μ is stable for μ < 1/10 and is characterized by two frequencies of oscillation: one nearly equal to the mean motion Ω (and representing the usual epicyclic motion of an eccentric orbit), and another

$$\begin{equation} |\lambda | = \sqrt{3\mu } \, \Omega. \end{equation} \tag{ 7 }$$

This slower frequency describes the libration of the particle on trajectories that are shaped like horseshoes. To sketch these horseshoes, we switch to polar coordinates where the equations of motion read

$$\begin{eqnarray} &&\ddot{r} - r\dot{\theta }^2 -2\Omega r\dot{\theta } = \frac{\partial U}{\partial r} \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} &&r\ddot{\theta } + 2\dot{r}\dot{\theta } +2\Omega \dot{r} = \frac{1}{r}\frac{\partial U}{\partial \theta } \end{eqnarray} \tag{ 9 }$$

with

$$\begin{eqnarray} U & = & \frac{\Omega ^2}{2} r^2 + \frac{1-\mu }{\sqrt{r^2 + \mu ^2 + 2 r \mu \cos \theta }} \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} & \simeq & \frac{3}{2} + \frac{3}{2}\Delta ^2 - \mu (1 + \cos \theta). \end{eqnarray} \tag{ 11 }$$

In Equation (11), we have expanded U in small radial (but not small azimuthal) displacements Δ ≡ r − 1 about the stable point. We insert Equation (11) into Equations (8) and (9), keeping only leading-order terms and taking d/dt ≪ Ω to filter out fast epicycles; then,

$$\begin{eqnarray} &&{-}2\dot{\theta } = 3\Delta \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &&\ddot{\theta } = -6\mu \frac{\partial }{\partial \theta } \sin ^{2}(\theta /2). \end{eqnarray} \tag{ 13 }$$

Multiplying Equation (13) by $\dot{\theta }$, integrating over time, and substituting Equation (12), we obtain the following "shape function":

$$\begin{equation} \Delta ^2 = \frac{16}{3}\mu [B - \sin ^{2}(\theta /2)]. \end{equation} \tag{ 14 }$$

Here, B is a constant of integration that takes values between 0 (zero amplitude libration at θ = 0) and 1 (maximal libration). Figure 1 samples several B-values; evidently, the shape function (14) traces horseshoe-shaped trajectories enclosing the stable point. The horseshoes are radially widest at θ = 0:

$$\begin{equation} \max \Delta = \sqrt{\frac{16B\mu }{3}}. \end{equation} \tag{ 15 }$$

This width can be a large fraction of the radius; e.g., for B = 1 and μ = 0.01 (characteristic of a disk whose mass is ∼0.01 M_*), the full width is 2max Δ = 0.46.

Although the libration trajectories are shaped like horseshoes, the square-root scalings in Equations (7) and (15) for the libration frequency and width characterize tadpoles in the standard three-body problem. This kinship between our horseshoe orbits and conventional tadpoles is expected. In both types of motion, the turning points (which are located farthest from the positive x-axis) are effected by the offset host star, which exerts torques on the test particle that deflect it back toward the x-axis.

We construct a model gas disk based on ideas developed by Paczynski (1977) and applied in many galactic contexts (e.g., Binney et al. 1991; Chang et al. 2007). Gas is assumed to occupy closed, non-intersecting test particle orbits. Crossing orbits are forbidden because they lead to shocks and energy dissipation. Gas streamlines approximate test particle orbits insofar as gas sound speeds are smaller than orbital velocities and gas self-gravity is negligible.

Working in the rotating frame of Section 2, we set μ = 0.01, M_* = M_☉, and $\Omega = \sqrt{\vphantom{A^A}\smash{{{GM}_{\odot }/a_0^3}}}$, where a₀ = 100 AU is the disk's characteristic radius. The star is positioned at x_* = −μa₀ = −1 AU and y_* = 0. Starting on the x-axis, we launch a test particle with an initial velocity $\dot{y} < 0$ ($\dot{x} = 0$) and integrate its trajectory until it re-crosses the x-axis. The final position and velocity are required to match initial values to within 1 part in 10⁴; the initial $\dot{y}$ is varied until these conditions are met. The integrations are performed using Python's explicit Runge–Kutta integrator of order 8(5, 3) (dop853) with an accuracy setting of 10⁻⁸. A total of 22 initial x-positions, distributed between the stable point at x = x_eq = 99.3 AU and x = x_max = 127.4 AU, are selected for the construction of 22 orbits/streamlines; see Figure 2. Their shapes match well the horseshoe orbits found analytically in Section 2. In particular, we confirm that Equations (7) and (15) for the libration frequency and radial width are obeyed.

Zoom In Zoom Out Reset image size

Each streamline is defined by N = 300 points marking the test particle's position recorded at equal intervals of time. To compute the gas surface density, we assign each streamline a mass ∝cos [(π/2)(x − x_eq)/(x_max − x_eq)], where x identifies the streamline's starting position; we assign more mass to inner horseshoes than outer ones. Lines are drawn between adjacent points on neighboring streamlines so that the entire space is tiled by quadrilaterals (see Figure 1). Each quadrilateral of area dA is assigned a surface density Σ = dM/dA, where dM = streamline mass/N.

For the most part this procedure yields surface densities that vary smoothly from quadrilateral to quadrilateral. However, we find a few spikes in density near the disk's edges (outermost horseshoes). These are numerical artifacts that arise from our particular choice of streamlines. Because our gravitational potential is only weakly non-axisymmetric, it admits a large variety of orbits that satisfy our tolerance criteria for being closed and non-intersecting. Our selection of streamlines is thus not unique. We chose streamlines that seemed to minimize the high density spikes, but were not able to eliminate them. Fortunately, the afflicted regions, which are restricted to disk edges, are tiny and easily masked out.

We smooth Σ by first applying the Python function scipy.ndimage.interpolation.map_coordinates, which sub-samples to a finer non-regular r–θ grid using a local order-1 spline. We then establish Σ on a regular r–θ grid with scipy.interpolate.griddata, which uses a local linear interpolation. The final step is to normalize the disk mass so that the barycenter remains at the origin. Figure 2 shows our final gas disk, of mass M_disk = 0.028 M_☉.

To zeroth order, the gas velocities are the (interpolated) velocities of the test particles used to construct the original 22 streamlines. We add a first-order correction to account for gas pressure gradients. With a fractional error of $\sqrt{\mu }$, pressure forces are balanced by Coriolis forces (the flow is geostrophic); the velocity corrections read

$$\begin{eqnarray} &&v_{r,1} = -\frac{c_{\rm s}^{2}}{2\Omega r }\frac{\partial \ln \Sigma }{\partial \theta } \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&v_{\theta,1} = \frac{c_{\rm s}^{2}}{2\Omega }\frac{\partial \ln \Sigma }{\partial r} \end{eqnarray} \tag{ 17 }$$

and are evaluated by applying Python's gradient function to Σ. Minor numerical artifacts in the first-order velocities at disk edges are smoothed using a combination of median and Gaussian filters. The sound speed c_s is fixed at 0.37 km s⁻¹, appropriate for gas at 40 K.

We now add dust particles to our gas disk. A dust particle behaves like a test particle except that it also feels a gas drag acceleration −(v_d − v_g)/t_stop, which damps dust–gas relative velocities over an aerodynamic stopping time t_stop. The stopping time scales inversely as the volumetric gas density ρ, which we compute from Σ by assuming the disk has a constant vertical thickness h = c_s/Ω. Our convention is to assign every particle a Stokes parameter St ≡ Ωt_stop evaluated at the position of peak ρ (i.e., at the stable point x = x_eq, y = 0); the stopping time at any other location then scales inversely as the local ρ (which deviates from the peak value by at most a factor of ∼3). A total of 2100 dust particles are laid down with an initial surface density matching that of the gas disk and with initial velocities equal to the zeroth-order gas velocities. Dust trajectories are computed using the same Runge–Kutta integrator as was used to obtain the gas streamlines. Gas densities and velocities are computed by local linear interpolation.

Figure 3 shows the evolution of dust surface density for St = 0.01, 0.1, and 1 (equivalent to grain sizes of ∼ 1 mm, 1 cm, and 10 cm for an assumed bulk density of 1 g cm⁻³). Because gas feels pressure forces while dust does not, dust and gas velocities generally differ. Dust particles feel tailwinds/headwinds that drag them toward gas pressure maxima (e.g., Chiang & Youdin 2010). For our disk, the global pressure maximum is located at the stable point, and indeed particles with St = 0.01 and 0.1 collect there. We estimate the concentration time as follows. In a gas disk with a radial pressure gradient, the radial speed of a particle is

$$\begin{equation} v_{{\rm dust},r} \sim \frac{c_{\rm s}^2}{\Omega r} \frac{\partial \ln \Sigma }{\partial \ln r} \frac{{\rm St}}{1+{\rm St}^2}\protect, \protect \end{equation} \tag{ 18 }$$

valid for any St (e.g., Nakagawa et al. 1986). For our problem, $\partial \ln \Sigma / \partial \ln r \sim \pm 1/\sqrt{16\mu /3}$ for streamlines to the left/right of the stable point (see Figure 3); we have here used the fact that our disk has a characteristic radial width $\sqrt{16\mu /3} \,a_0$ (Equation (15)). Thus, particles interior/exterior to the stable point drift radially outward/inward, traversing a series of ever-smaller horseshoes. The timescale to cross the radial width of the disk and collect onto the stable point is then

$$\begin{equation} t_{\rm conc} \sim \frac{ \sqrt{16\mu /3} \, a_0}{v_{{\rm dust},r}} \sim \frac{16}{3} \frac{\mu }{\Omega } \left(\frac{\Omega a_0}{c_{\rm s}} \right)^2 \frac{1+{\rm St}^2}{{\rm St}}\protect, \protect \protect \protect \end{equation} \tag{ 19 }$$

which predicts that St ∼ 1 particles accumulate fastest. Our numerical experiments confirm this expectation and also verify the asymptotic scalings (t_conc∝1/St for St ≲ 0.1 and t_conc∝St for St ≳ 4; see Figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Surprisingly, not all particles concentrate at the global pressure maximum located on the x-axis (the gas disk symmetry axis). Figures 3 and 5 reveal that marginally aerodynamically coupled particles collect along an arc that extends off the x-axis in the direction of orbital motion. Starting from the x-axis at St ≲ 0.1, the position angles of final accumulation points advance with increasing St to a maximum of ∼45° at St ≈ 0.5; further increasing St causes the position angles to slide back down, returning to zero at St ≳ 10. The bottom right panel of Figure 3 illustrates that as long as the drag force is neither too strong nor too weak, a three-way force balance between the offset stellar gravity, the centrifugal force, and drag is possible, but only in one quadrant of the rotating frame. Particle size segregation off-axis offers an observational test of our model: the locations of peak intensity in dust emission maps may vary systematically with wavelength, to the extent that different observing wavelengths select for different size particles.

Zoom In Zoom Out Reset image size

A lopsided disk displaces its host star from the system barycenter. In the resultant "offset" stellar potential (one that has an indirect term), we have calculated that gas streamlines take the form of horseshoes. Our calculations are not self-consistent because on the one hand we have invoked the gas disk's mass to move the star, but on the other we have neglected the disk's gravity when computing gas streamlines. Our neglect of gas self-gravity is a severe approximation: the disk's gravity is as strong as the indirect forcing, since both scale as μ. Indeed, when we re-calculate a posteriori the potential U including the gas disk's contribution, we find that U's topology changes qualitatively. In particular, what was a global extremum at x = x_eq becomes a saddle point. Thus, the mode we have constructed in this Letter is not an equilibrium fast mode. Still, we might not be too far off the mark; reducing the disk mass by a factor of ∼2 (while keeping the stellar offset μ fixed) restores the original topology at x_eq. We are also encouraged by stellar dynamics calculations performed by Dury et al. (2008), who found that their "rotating lopsided mode" (i.e., a fast m = 1 mode) supports "banana" (i.e., horseshoe-shaped) orbits; see their Figure 14. Future work should incorporate disk self-gravity, not only to construct self-consistent fast modes, but also to search for instabilities that can create lopsided mass distributions.

One step toward more realistic gas equilibria is to replace our test-particle calculation with a hydrodynamical one that accounts for gas pressure. In this regard, we experimented with a two-dimensional inviscid hydro-simulation using the PLUTO code (Mignone et al. 2007), evolving an initially circular and uniform gas ring in an offset stellar potential (see the caption of Figure 2(b) for technical details). A few numerical difficulties complicated this first-cut experiment. The system did not settle into steady state; gas leaked outward and was lost off the grid; and our boundary conditions, chosen mainly for convenience, affected the propagation of spiral density waves. Nevertheless, ∼20–80 orbits into our hydro-simulation, we observed the gas disk conforming to the same horseshoe shapes obtained with our test-particle method—compare Figures 2(a) and (b)—and found consistent results for dust concentration. The qualitative agreement between our test-particle calculation and our hydro-simulation helps to validate the former, and motivates us to improve upon the latter by including gas self-gravity, viscosity, turbulence, and dust feedback.

Fast modes promise to better reproduce observed dust distributions with large azimuthal and radial extents (e.g., Casassus et al. 2013; van der Marel et al. 2013). In our demonstration model, particles with sizes of 0.1–1 mm have concentration times on the order of the disk age or longer, and should thus occupy horseshoe-shaped regions like those shown in the top right panel of Figure 3. In addition, larger particles, having Stokes parameters St ∼ 0.1–10 (particle sizes of 1–100 cm), naturally spread themselves in azimuth by ∼45°. Gas turbulence can also diffuse dust particles, keeping them on large horseshoe orbits. The Péclet number is the ratio of timescales for turbulent diffusion and concentration:

$$\begin{eqnarray} &&{\rm Pe} \equiv \frac{t_{\rm diff}}{t_{\rm conc}} \sim \frac{(\sqrt{16\mu /3} \, a_0)^2 / D}{t_{\rm conc}} \sim \frac{{\rm St}}{\alpha }. \end{eqnarray} \tag{ 20 }$$

Here, D ≈ D_g/(1 + St²) is the dust particle diffusivity (Equation (5) of Youdin & Lithwick 2007), D_g ∼ αc_sh is the gas diffusivity, and α < 1 is the Shakura–Sunyaev turbulence parameter. Small Pe ≲ 1—i.e., more spatially extended dust—obtains for St ≲ 0.01 (millimeter-sized or smaller particles) and α ∼ 0.01.

Our proposed lopsided gas mode is gravitational and global; hence, its radial width is not limited by the pressure scale height H. This limitation afflicts vortices, since the Keplerian shear across a vortex must remain subsonic (e.g., Lithwick 2009; Zhu & Stone 2014). Lyra & Lin (2013) acknowledge this difficulty in the last paragraph of their Section 6.5, noting that the radial half-width of the dust trap in Oph IRS 48 is 17 AU ≈3H, whereas their vortex theory (see their Equation (68)) predicts a maximum half-width of ∼H/2. Our horseshoe-shaped orbits are not constrained by H and may thus better reproduce the azimuthally and radially wide horseshoe-shaped emission seen in transitional disks like HD 142527 (Casassus et al. 2013). Because the gas mode is global, we speculate that it may be more resistant to the destructive effects of dust back-reaction than are vortices (Johansen et al. 2004; Meheut et al. 2012). Disk gravity underlies our mode but may undercut vortices; Lin & Papaloizou (2011a, 2011b) found that self-gravity shrinks vortices and frustrates their merging.

We have not identified the origin of the fast mode. It may be that, like the vortices of Zhu & Stone (2014), a planet is responsible. A passing star might also pull the disk to one side of the star and excite a long-lived fast mode (see Jalali & Tremaine 2012, who use a passing star to excite slow modes). However, for economy of hypothesis, one can do no better than look to the self-gravitational field of the gas disk itself (Dury et al. 2008).

We thank Robin Dong, Mir Abbas Jalali, Eve Lee, Chris Ormel, Roman Rafikov, and Scott Tremaine for useful discussions, and an anonymous referee for an encouraging report. E.C. acknowledges support from the Miller Institute and a NASA Origins grant.