Derive the drift and the gas response

Choose a rotating patch and a stopping time

Use local coordinates $x$ radially outward, $y$ along orbital rotation and $z$ vertically. The patch rotates with angular frequency $\Omega$ . Velocities are measured relative to the local Keplerian shear, so a freely orbiting particle at the reference radius has zero peculiar velocity. Denote gas velocity by $\mathbf u$ and dust velocity by $\mathbf v$ .

A fixed linear drag law makes dust acceleration $-(\mathbf v-\mathbf u)/t_{\rm stop}$ . Define $\mathrm{St}=\Omega t_{\rm stop}$ and local volume loading $\epsilon=\rho_{\rm d}/\rho_{\rm g}$ . The same Stokes number can correspond to different physical sizes in different gas conditions. Gas receives the opposite force per volume, so its drag acceleration has an extra factor $\epsilon$ .

Reading: Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1 · Benítez-Llambay, Krapp & Pessah 2019, FARGO3D multifluid tests, v3

Add the pressure force only to gas

For a uniform radial pressure gradient, write its outward acceleration as $2\eta v_{\rm K}\Omega$ . Positive $\eta$ corresponds to pressure decreasing outward and gas rotating more slowly than Keplerian motion. Neglect radial derivatives of the equilibrium velocities, vertical gravity and viscosity.

Steady radial and azimuthal force balances provide four equations. The factors two and one half combine Coriolis acceleration with the Keplerian background shear. Pressure appears in the gas equation only.

0=2v_y-\frac{v_x-u_x}{\mathrm{St}},\qquad0=-\frac{v_x}{2}-\frac{v_y-u_y}{\mathrm{St}}

Dust balances rotation and drag. Here time has been divided by the inverse angular frequency.

0=2u_y+2\eta v_{\rm K}+\frac{\epsilon(v_x-u_x)}{\mathrm{St}},\qquad0=-\frac{u_x}{2}+\frac{\epsilon(v_y-u_y)}{\mathrm{St}}

Gas also receives radial pressure support. Equal and opposite drag cancels when the equations are mass-weighted and added.

Reading: Youdin & Goodman 2005, Streaming Instabilities, arXiv v1 · Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1

Solve the coupled balance

Add the azimuthal gas equation to $\epsilon$ times the dust equation. Drag cancels, giving $u_x=-\epsilon v_x$ . The radial sum similarly gives $u_y+\epsilon v_y=-\eta v_{\rm K}$ .

The dust radial equation gives $v_y=(1+\epsilon)v_x/(2\mathrm{St})$ . Its azimuthal equation gives $u_y=v_y+\mathrm{St}v_x/2$ . Substitute both into the weighted radial relation. The coefficient of $v_x$ becomes $[(1+\epsilon)^2+\mathrm{St}^2]/(2\mathrm{St})$ , leading to the velocities below.

Q=(1+\epsilon)^2+\mathrm{St}^2,\qquad v_x=-\frac{2\mathrm{St}}Q\eta v_{\rm K},\qquad u_x=\frac{2\epsilon\mathrm{St}}Q\eta v_{\rm K}

Dust goes inward and gas outward for a positive pressure-support parameter. Their mass-weighted radial flow cancels in this uniform equilibrium.

v_y=-\frac{1+\epsilon}{Q}\eta v_{\rm K},\qquad u_y=-\frac{1+\epsilon+\mathrm{St}^2}{Q}\eta v_{\rm K}

These are small orbital-speed offsets. Add the Keplerian motion to obtain total orbital velocity.

Reading: Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1

Work a drift-time example

Take an illustrative orbit of radius $1.50\times10^{11}\,\mathrm m$ with $v_{\rm K}=3.00\times10^4\,\mathrm{m\,s^{-1}}$ and $\eta=10^{-3}$ . For negligible feedback and $\mathrm{St}=1$ , the radial drift speed is $-30.0\,\mathrm{m\,s^{-1}}$ . Dividing radius by speed gives a local drift-time scale of about 158 years. This is an order-of-magnitude time scale at that radius, not an integrated trajectory through a changing disk.

Now take $\epsilon=1$ at the same stopping time. Then $Q=5$ , dust moves inward at $12.0\,\mathrm{m\,s^{-1}}$ and gas outward at the same speed. Their equal densities make the net radial mass flow zero. The dust’s radial drift is reduced because gas responds to the drag.

Let a prepared pressure peak evolve

A finite trap has spatial gradients, so the local steady formula alone cannot evolve it. Instead use continuity $\partial_t\rho+\partial_x(\rho v_x)=0$ for each phase, transport momentum, apply gas pressure $P=c_s^2\rho_{\rm g}$ and exchange drag. In the radial experiment, length is in $H=c_s/\Omega$ and velocity in $c_s$ . These differ from the $\eta r$ and $\eta v_{\rm K}$ scales of the later linear-wave experiment.

Initially the gas density is proportional to $\exp[A\cos(2\pi(x-x_{\rm peak})/L)]$ on a periodic domain $L=8H$ , normalized to mean one. Dust begins uniform, and local coupled drift balances supply initial velocities. This is a prepared initial condition, not an exact global steady solution: gradients subsequently alter both phases. The gas peak is allowed to deform.

The actual frozen results below use 128 cells, an eight-inverse-angular-frequency warm-up and the same saved gas/dust fields that generate the scene. The table reports the concentration reached by the evolved gas and dust.

Actual synthetic finite-volume results after eight inverse angular frequencies:128periodic cells, gas sound speed one, domain eight gas scale heights, prepared log-density contrast0.6, initially uniform dust. Source freeze-radial.ts saves both density and momentum fields and resolution/time-step comparisons.
Case	Stokes number	Initial mean loading	Dust peak / mean	Largest gas radial speed · sound speed
Default trap	$0.300$	$0.300$	$1.74$	$0.00508$
Tight coupling	$0.0300$	$0.300$	$1.07$	$0.000986$
Rapid drift	$1.00$	$0.300$	$2.97$	$0.00467$
Strong feedback	$0.300$	$3.00$	$1.15$	$0.0121$

Reading: Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1 · Benítez-Llambay, Krapp & Pessah 2019, FARGO3D multifluid tests, v3

Separate physical feedback from numerical smoothing

The finite-volume solver exchanges every face flux between neighboring cells, including the periodic boundary. Each phase therefore keeps its mass. During a local drag step the center-of-mass velocity stays fixed while the relative velocity decays as $\exp[-(1+\epsilon)\Delta t/t_{\rm stop}]$ . For drag alone, relative kinetic energy decreases.

\frac{dE_{\rm rel}}{dt}=-\frac{\rho_{\rm d}}{t_{\rm stop}}|\mathbf v-\mathbf u|^2,\qquad E_{\rm rel}=\frac{\rho_{\rm g}\rho_{\rm d}}{2(\rho_{\rm g}+\rho_{\rm d})}|\mathbf v-\mathbf u|^2

Drag dissipates relative energy while conserving mixture momentum. The evolving shear and gas pressure can supply motion, and the isothermal model removes the heat.

Reading: Benítez-Llambay, Krapp & Pessah 2019, FARGO3D multifluid tests, v3 · Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1

Check the scales you have resolved

The radial solver uses first-order spatial Rusanov fluxes, second-order SSP Runge–Kutta transport and symmetric exact drag/epicycle source steps. Its spatial numerical diffusion smooths sharp structures. At the default time, 128 cells differ from a 1024-cell reference by 2.92 percent in relative L1 dust density, averaged over the whole periodic domain; the peak is about 6.1 percent lower. Those are measured resolution differences, not uncertainty bars on a real disk.

A pressureless fluid also assigns one dust velocity to each point. It cannot represent crossing particle streams or dense multistream caustics. The finite early experiment is useful for trapping and gas feedback; later nonlinear concentrations need a particle or more complete kinetic treatment. Nothing in the radial equations supplies the vertical mode needed for linear streaming instability. The next chapter introduces that missing direction.