Streaming Torque in Dust-Gas Coupled Protoplanetary Disks

Previous studies by Squire & Hopkins (2018a, b) interpreted the YG streaming instability (Youdin & Goodman, 2005) as the resonance between epicyclic oscillations of gas and the drift velocity, termed as “Resonant Drag Instability” (RDI). Furthermore, they suggest that the drift velocity can resonate with various waves, including sound waves (Hopkins & Squire, 2018). RDI occurs when the dust drift velocity exceeds the sound speed ($W_{s}>c_{s}$), leading to the emergence of instabilities. In our scenario, where $W_{s}\ll c_{s}$, RDI is absent. However, the planet can still excite certain modes.

Through solving (9)-(14) with relaxation method, we can get as many sets of solutions as the number of $k_{y}$. For a specific case with $k_{y}=0.3$, Figure 1 illustrates the excitation of dusty and gaseous waves by a planet situated at the coordinate origin. And we select $f_{d}=0.01$, $\tau=1.0$ for a clear illustration. The perturbed densities and velocities are normalized by $\Sigma_{0}=qh_{p}^{-3}\Sigma_{p}$ and $V_{0}=qh_{p}^{-3}\Omega_{p}$, where $q$ represents the mass ratio of the planet to the host star ($M_{p}/M_{\star}$). In the figure, we notice that such a dust mass fraction has a minimal impact on gas, as depicted in subfigures (b), (d), and (f). The results resemble those calculated in a pure gas disk (e.g. Tanaka et al., 2002; Rafikov & Petrovich, 2012), except for the short-wavelength structure observed inside the planet. This short-wavelength structure arises due to the feedback of dust, leading us to refer to the gaseous waves as “quasi-density” waves.

Dusty perturbations exhibit distinct characteristics compared to gaseous perturbations. Subfigures (a), (c), and (e) illustrate that dusty perturbations display a short-wavelength structure inside the planet. It originates from the “quasi-drift” mode (Hopkins & Squire, 2018). In next subsection, we will solve its radial wavenumber and demonstrate that the appearance of short-wavelength structure is robust. Identifying these features necessitates a high spatial resolution. The lower right zoom-in panel in subfigures (a) shows that the radial wavelength near the inner boundary is about $2\times 10^{-2}\mathrm{H}$, which poses a significant challenge for resolution in hydrodynamic simulations. We use $2\times 10^{-4}\mathrm{H}$ resolution, corresponding to $10^{5}$ mesh points, to ensure that the refined structure is presented in our calculations. With such a high resolution, the solutions become converged and have no dependency on the radial domain extent or mesh points.

Away from the planet, dusty waves are dominated by “quasi-density” waves, except for the density perturbation inside the planet as shown in subfigure (a). This dominance is due to the slow damping of the quasi-drift mode and the small amplitude of the dusty quasi-density wave. The latter is depicted in the upper right zoom-in panel, where $\delta\Sigma_{d}\sim f_{d}\delta\Sigma_{g}$. As $f_{d}$ increases, the quasi-drift mode is damped faster, and quasi-density waves also come to dominate wavy behaviors further inside the planet. This will be further illustrated in the next subsection.

Hopkins & Squire (2018) derived the dispersion relation for dust-gas systems in the free-falling frame. The free drift mode exhibits a dispersion relation of $\omega=\boldsymbol{W}_{s}\cdot\boldsymbol{k}$. Similar results can be obtained in our framework. Employing WKB approximation, we numerically solve for the radial wavenumber of the quasi-drift mode. We obtained a fitting dispersion relation, $\tilde{\omega}=\left(1+\mu\right)\boldsymbol{W}_{s}\cdot\boldsymbol{k}$. The results are depicted in Figure 2, where the upper(lower) panel displays the real(imaginary) part of the radial wavenumber, $k_{x}$. In the upper row, the red solid line represents our numerical results. The blue dot-dash line indicates $k_{x}$ of the free drift mode, i.e. $k_{x}=\left(\tilde{\omega}_{d}-k_{y}W_{s,y}\right)/W_{s,x}$, where $\tilde{\omega}_{d}$ represents the Doppler-shifted frequency of dust, defined by $k_{y}r(\Omega_{p}-\Omega_{d})$. The green dot-dash line indicates $k_{x}$ of the quasi-drift mode, i.e. $k_{x}=\left(f_{g}\tilde{\omega}_{d}-k_{y}W_{s,y}\right)/W_{s,x}$. The quasi-drift mode aligns perfectly with our numerical calculations, indicating that these short-wavelength structures arise from this mode. Furthermore, $|k_{x}|$ increases linearly as $|x|$ increases. Specifically, at the inner boundary ($x=-10$), $k_{x}\sim 300$, a value consistent with the radial wavelength of $\sim 2\times 10^{-2}\mathrm{H}$ as depicted in Figure 1. Lastly, both the phase velocity and group velocity of the mode are $\sim\boldsymbol{W}_{s}$, indicating that the wave propagates inward in the $-x$ direction and ahead of the planet in the $+y$ direction.

From the lower row, we find that the quasi-drift mode is excited by the planet and propagates inward while undergoing rapid damping, particularly at locations slightly further away from the planet. This behavior is consistent with the findings of Hopkins & Squire (2018). In the absence of planets, the mode remains stable due to $W_{s}\ll c_{s}$. Consequently, its structures dominate the region near the interior of the planet, displaying asymmetry in the shearing sheet. When $f_{d}=0.01$, the damping rate is small, allowing the quasi-drift mode to dominate wavy behaviors far inside the planet. While $f_{d}>0.01$, the damping rate increases rapidly away from the planet, causing quasi-density waves to dominate wavy behaviors, as illustrated in Figure 4.

To visualize the waves in real space, we integrate density perturbations. Conclusively, the planetary wakes are depicted in Figure 3 with $f_{d}=0.01$. From left to right, we consider $\tau=0.1,1.0,3.0$, with every two panels displaying the dusty quasi-drift mode and gaseous quasi-density waves separately (note their ranges of $x$ axes differ). All perturbations in each panel are normalized by their respective maximum values. Comparing panels (b), (d), and (f), we find minimal changes in gaseous quasi-density waves. However, the dusty quasi-drift mode exhibits a short radial length-scale ($\sim 0.1\mathrm{H}$). These results align with our previous analysis indicating that the quasi-drift mode is rapidly damped as it propagates inward. As we talked before, both the phase velocity and group velocity of the mode are $\sim\boldsymbol{W}_{s}$. So, the quasi-drift mode remains ahead of the planet, exerting a positive torque on it, which is consistent with panels (a), (c) and (e). Due to the large amplitude of the mode, this effect is significant, termed as the “Streaming Torque (ST)”. When $\tau\sim 1$, the dusty wake exhibits a larger radial extent and relative density perturbation, indicating that the torque effect is most pronounced around $\tau\sim 1$, where the drift velocity, $\boldsymbol{W}_{s}$, reaches its maximum. Detailed calculations are presented in the subsequent section.

The morphologies of dust wave are consistent with Benítez-Llambay & Pessah (2018), who conducted hydrodynamic simulations for a dust-gas PPD. Their results present radial short wavelength structures inside the planet for as their figures 1 and 4 show. Dust surplus, deficit and multiple dust gaps inside the plane might imply the nonlinear evolution phase of dust drift mode.

For low-mass planets embedded within pure gas PPDs, the Lindblad torque (LT), $\Gamma_{\mathrm{L}}$, and corotation torque (CT), $\Gamma_{\mathrm{C}}$, exerted by the disk on the planet play crucial roles in driving type I migration (Ward, 1986; Goldreich & Tremaine, 1979; Korycansky & Pollack, 1993; Tanaka et al., 2002). These torques can be quantified by analyzing the AMF carried by density waves, given by:

Then LT and CT, are given by

For a two-fluid system, distinct AMFs arise for dusty and gaseous waves. Figure 4 illustrates the density perturbations and AMFs for both components with $\tau=1.0$, $k_{y}=0.3$, and varying $f_{d}$. And the normalization is $F_{0}=q^{2}h_{p}^{-3}\Sigma_{p}r_{p}^{4}\Omega_{p}^{2}$. The first row, consistent with Figure 1, demonstrates minimal dusty AMF due to a low $f_{d}$. Gaseous AMF tends towards a constant at large distances from the planet. However, when $f_{d}$ increases, as shown in subsequent rows, both dusty and gaseous AMFs diminish progressively, indicative of damped “quasi-density” waves.It means we cannot get the disk torque on planets through (21)-(25). Actually, Equation (22) does not apply in dissipation flows (Balbus, 2003). In the system, dust and gas interacts through friction. Then, part of kinetic energy is dissipated into heating. But we use isothermal state for gas. The cooling process would exert a thermal damping on propagating waves. Similar behaviors appear in Miranda & Rafikov (2020) because of thermal damping. However, our current investigation focuses on low $f_{d}$ and the quasi-drift mode, which exhibits substantial torques. So, we do not go into detail about the damped AMFs. Notably, the quasi-drift mode experiences heightened damping as $f_{d}$ increases, which is consistent with the discussion in section 3.2. Ultimately, we refrain from utilizing Equations (21) to (25) for torque calculations. This decision stems from the inherent limitation of this method in capturing ST effects accurately.

We can calculate the torque by integrating planetary force, which reads

In Fourier space,

We use Equation (27) with $k_{y}$ ranging from $0.01$ to $15$, logarithmically discrete to 320 points. The convergence for total torques is confirmed with high spatial resolution.

The results are presented in Figure 5 for $f_{d}=0.01$. The $x$-axis represents different $\tau$ values. The solid dots represent our calculations, linked by three curves. The blue curve indicates torques from dusty modes, labeled as “Dusty Torque”, while the green curve represents those from gaseous modes, labeled as “Gaseous Torque”. The black curve depicts the total torque. Here, $h_{p}\Gamma_{0}=q^{2}h_{p}^{-2}\Sigma_{p}r_{p}^{4}\Omega_{p}^{2}$. We set $h_{p}=0.03$ for comparison with Tanaka et al. (2002).

The torque analysis in Figure 5 reveals several key insights. Dusty torques predominantly exhibit positive and dominate total torques in most scenarios. Conversely, gaseous torques appear negative, driving inward migration for planets solely.

Dusty torques reach extreme points when $\tau\sim 1$, corresponding to the peak drift velocity. The dominance of dusty torques is attributed to the ST, which surpasses contributions from CT and LT for dust. Notably, when $\tau$ is sufficiently low ($\lesssim 10^{-1}$), dusty torques become negative due to strong gas-dust coupling, causing dust behavior to mimic that of gas.

The negative nature of gaseous torques aligns with their role in facilitating inward migration, a phenomenon extensively studied in prior literature (e.g. Ward, 1986, 1997; Korycansky & Pollack, 1993; Tanaka et al., 2002). If we consider $\tau$ to be very low or large, the gaseous torque must converge to that in a pure gas PPD. Tanaka et al. (2002) employed the “Modified Local Approximation”, which incorporates the curvature effect, surface density gradient, and pressure gradient compared to the shearing sheet. In this paper, only the effect of the pressure gradient is considered, resulting in a change of the rotation curve from Keplerian frequency. This alteration affects the locations of Lindblad resonance (LR) points, leading to a net LT. Our calculations yield $\Gamma_{g}\simeq-2.0h_{p}\Gamma_{0}$, which is consistent with the value obtained by Tanaka et al. (2002) of $\Gamma_{g}\simeq-2.3h_{p}\Gamma_{0}$, attributed solely to the pressure gradient. When $\tau\sim 1$, the gaseous radial velocity amplifies CT, leading to maximum negative gaseous torque. This phenomenon is reminiscent of type III migration, originating from planetary migration (Masset & Papaloizou, 2003; Ogilvie & Lubow, 2006; Paardekooper, 2014).

It is then easy to estimate the planetary migration timescale. Using some typical parameters, Tanaka et al. (2002) found it would be typically $10^{6}$ years with inward migration, which is comparable to or shorter than the lifetime of PPDs. Tanaka et al. (2002) gives the (outward) migration timescale

where $L_{p}$ is the angular momentum of the planet $M_{p}(GM_{\star}r_{p})^{1/2}$. With $f_{d}=0.01$, Figure 5 allows $\Gamma$ ranges from $-3h_{p}\Gamma_{0}$ to $30h_{p}\Gamma_{0}$. Using the same parameters as those used in Tanaka et al. (2002), the final result would be $-10^{6}$ to $10^{5}$ years (the negative sign means inward migration), depending on $\tau$. This illustrates how the ST can dominate both the CT and LT, consequently driving outward planetary migration.