Disentangling Planets from Photoelectric Instability in Gas-rich Optically Thin Dusty Disks

We work in the thin disk approximation using the vertically averaged equations of hydrodynamics. Our model includes gas and dust; the gas is evolved in an Eulerian grid, whereas the dust is treated as Lagrangian particles. The equations solved are

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{g}}{\partial t}=-\left({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}\right){{\rm{\Sigma }}}_{g}-{{\rm{\Sigma }}}_{g}{\boldsymbol{\nabla }}\cdot {\boldsymbol{u}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}=-\left({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}\right){\boldsymbol{u}}-\displaystyle \frac{1}{{{\rm{\Sigma }}}_{g}}{\boldsymbol{\nabla }}P-{\boldsymbol{\nabla }}{\rm{\Phi }}-\displaystyle \frac{{{\rm{\Sigma }}}_{d}}{{{\rm{\Sigma }}}_{g}}\displaystyle \frac{\left({\boldsymbol{u}}-{\boldsymbol{v}}\right)}{{\tau }_{{\rm{f}}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial S}{\partial t}=-\left({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}\right)S-\displaystyle \frac{{c}_{V}}{T}\displaystyle \frac{\left(T-{T}_{p}\right)}{{\tau }_{{}_{T}}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\rm{\Phi }}=\displaystyle \sum _{i}^{n}\displaystyle \frac{{{GM}}_{i}}{\sqrt{{{R}_{i}}^{2}+{{b}_{i}}^{2}}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&P=\rho {{c}_{s}}^{2}/\gamma ,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{x}}}{{dt}}={\boldsymbol{v}},\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{v}}}{{dt}}=-{\boldsymbol{\nabla }}{\rm{\Phi }}-\displaystyle \frac{\left({\boldsymbol{v}}-{\boldsymbol{u}}\right)}{{\tau }_{{\rm{f}}}}.\end{eqnarray} \tag{ 7 }$$

In the equations above, Σ_g and Σ_d are the vertically integrated gas and bulk densities of dust, respectively. Their quotient defines the dust-to-gas ratio ε = Σ_d/Σ_g.⁵ The velocity of gas parcels is given by ${\boldsymbol{u}};$ P stands for the vertically integrated pressure, Φ is the gravitational potential, and ${\boldsymbol{v}}$ is the velocities of dust grains. The gravitational potential includes contributions from the star and embedded planets, treated as massive particles with an N-body evolution dynamically integrated with the hydrodynamics evolution. In Equation (4), G is the gravitational constant, M_i is the mass of particle i, and ${R}_{i}=| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{{p}_{i}}|$ is the distance relative to particle i, located at ${{\boldsymbol{r}}}_{{p}_{i}}$. The quantity b_i is the distance over which the gravity field of particle i is softened to prevent singularities. In Equation (3), S is the gas entropy, given by

$$\begin{eqnarray}&&S={c}_{{}_{V}}\left[\mathrm{ln}\left(\displaystyle \frac{P}{{P}_{0}}\right)-\gamma \mathrm{ln}\left(\displaystyle \frac{{{\rm{\Sigma }}}_{g}}{{{\rm{\Sigma }}}_{g0}}\right)\right],\end{eqnarray} \tag{ 8 }$$

where ${c}_{{}_{V}}$ is the specific heat at constant volume, $\gamma ={c}_{p}/{c}_{{}_{V}}$ is the adiabatic index, and c_P is the specific heat at constant pressure; P₀ and Σ_g0 are the reference pressure and gas density values, respectively. In Equation (3), T stands for the gas temperature,

$$\begin{eqnarray}&&T=\displaystyle \frac{{c}_{s}^{2}}{{c}_{p}\left(\gamma -1\right)},\end{eqnarray} \tag{ 9 }$$

where c_s is the sound speed; ${\tau }_{{}_{T}}$ is the thermal coupling time between gas and dust. The quantity T_p is the gas temperature set by equilibrium between the different heating and cooling processes, including photoelectric heating. Calculating it accurately necessitates a detailed radiative transfer that is beyond the scope of the present work. Instead, we use a simple prescription for the gas temperature set by photoelectric heating as used in Klahr & Lin (2005),

$$\begin{eqnarray}&&{T}_{p}={T}_{0}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{d}}{{{\rm{\Sigma }}}_{g0}}\right)}^{\beta },\end{eqnarray} \tag{ 10 }$$

where T₀ is a reference temperature. We keep β = 1 (see the Appendix of Klahr & Lin 2005 and Figure 1 of that work for the estimate of β ≈ 1 given by a balance between photoelectric heating of silicate grains and cooling via the C ii line). In the Lagrangian dust equations, ${{\boldsymbol{v}}}_{p}$ is the velocity of a dust grain, and ${{\boldsymbol{x}}}_{p}$ is its position. The last term in Equation (7) is the drag force by which gas and dust interact dynamically. It brings the velocities of the dust grains to that of the gas within the timescale ${\tau }_{{\rm{f}}}$, the friction time. The last term in Equation (2) is the backreaction of the drag force. We choose the friction time as proportional to the local orbital period

$$\begin{eqnarray}&&{\tau }_{{\rm{f}}}={{\tau }_{{\rm{f}}}}_{0}\left(\displaystyle \frac{{{\rm{\Omega }}}_{0}}{{\rm{\Omega }}}\right),\end{eqnarray} \tag{ 11 }$$

as in this case, the relative timescales (the PEI growth rate also scales with the orbital period) are similar at all radii.⁶ We keep ${{\tau }_{{\rm{f}}}}_{0}{{\rm{\Omega }}}_{0}=1$. The operator

$$\begin{eqnarray}&&\displaystyle \frac{D}{{Dt}}=\displaystyle \frac{\partial }{\partial t}+{\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}\end{eqnarray} \tag{ 12 }$$

represents the advective derivative.

Given that the goal of this work is to examine the competing roles of the PEI and gravitational perturbation by a planet, we ignore the effects of radiation pressure. Because the thermal time is expected to be very low, we use the instantaneous thermal coupling approximation (Lyra & Kuchner 2013). That means that instead of solving the energy equation, we equate T = T_p and update the sound speed accordingly. This change effectively amounts to choosing a new equation of state that depends on the dust density,

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{{\tau }_{{}_{T}}\to 0}P={c}_{{}_{V}}{T}_{0}(\gamma -1)\ \left(\displaystyle \frac{{{\rm{\Sigma }}}_{g}}{{{\rm{\Sigma }}}_{g0}}\right){{\rm{\Sigma }}}_{d}.\end{eqnarray} \tag{ 13 }$$

We solve the equations with the Pencil Code,⁷ which integrates the evolution equations with sixth-order spatial derivatives and a third-order Runge–Kutta time integrator. Sixth-order hyperdissipation terms are added to equations of motion to provide extra dissipation near the grid scale. They are needed because the high-order scheme of the Pencil Code has little overall numerical dissipation (McNally et al. 2012). The calculations are sped up by using the "fargo" orbital advection algorithm (Masset 2000; Lyra et al. 2017).

The simulations are performed on a uniform cylindrical grid (coordinates r, ϕ) with resolution (N_r, N_ϕ) = (528, 528), linear domain L_r = [0.4, 2.5], and full 2π coverage in azimuth. To represent the dust grains, 200,000 Lagrangian superparticles are scattered uniformly across the disk. The star and planet orbit around the barycenter of the system in circular orbits. Disk self-gravity is ignored. We ran simulations with three different planet masses. The mass ratios of the planet to the host star are q = 10⁻³, 10⁻⁴, and 10⁻⁵. For each planet, we simulated five different dust-to-gas ratios, from 10⁻⁴ to 1, logarithmically spaced. The runs from ε = 10⁻² to 1 are in the zone of highest growth rate for the PEI (Lyra & Kuchner 2013; Richert et al. 2018). We choose our units such that

$$\begin{eqnarray}&&{{GM}}_{\star }={{\rm{\Omega }}}_{0}={r}_{0}=1.\end{eqnarray} \tag{ 14 }$$

All quantities with subscript "0" are measured at r₀. For r₀ = 100 au and M_⋆ = 1M_⊙, our unit of time is ≈160 yr, and the unity of velocity is ≈3 km s⁻¹. We investigate disks of aspect ratios h ∈ {0.05, 0.1, 0.2}, where h = H/r, with H = c_s/Ω the disk scale height. For a solar mixture of molecular hydrogen and helium of mean molecular weight μ = 2.5 and γ = 7/5, these correspond to temperatures of T = 5, 20, and 80 K, respectively. We consider monodispersive particles of radius a_• = 1 μm. This choice, along with the choice of Stokes number unity,

$$\begin{eqnarray}&&\mathrm{St}\equiv {\tau }_{{\rm{f}}}{\rm{\Omega }}=1,\end{eqnarray} \tag{ 15 }$$

sets the unit of density.⁸ This is because the friction time in the Epstein regime is

$$\begin{eqnarray}&&{\tau }_{{\rm{f}}}=\displaystyle \frac{{a}_{\bullet }{\rho }_{\bullet }}{\rho {c}_{s}},\end{eqnarray} \tag{ 16 }$$

where ρ_• is the grain internal density. For icy grains, it should be ρ_• ≈ 1 g cm⁻³. Multiplying Equation (16) by Ω to obtain the Stokes number and setting to 1 yields, solving for density,

$$\begin{eqnarray}&&\rho \,{| }_{\mathrm{St}=1}=\displaystyle \frac{{a}_{\bullet }{\rho }_{\bullet }}{{c}_{s}},\end{eqnarray} \tag{ 17 }$$

which, for h = 0.1, results in ρ ≈ 10⁻⁸ g cm⁻³, with an associated column density of ${\rm{\Sigma }}=\sqrt{2\pi }\rho H\approx 2.5\times {10}^{-4}$ g cm⁻². We keep the gas and dust densities constant with radius, whereas the temperature decreases according to T = T₀r^−b. We usually use b = 1 corresponding to a negative global pressure gradient, causing the dust to slowly drift inward. Although a drift is present, the resolution is not high enough to capture the streaming instability (Youdin & Goodman 2005; Johansen et al. 2007; Lyra & Kuchner 2013). A subset of simulations uses b = 0.

The parameters of the simulations we ran are shown in Table 1. We explored the parameter space of planet mass, sound speed, and dust-to-gas ratio.

To investigate the effect of the PEI, we ran 36 simulations. We used two planets, analogous to a Neptune- and a Jupiter-mass planet. For each planet, we varied three disk temperatures, two dust-to-gas ratios, and three different cases for carving a gap in the disk. The cases are the effect of a single planet with no PEI, a disk with the instability on its own, and a combination of an embedded planet and the instability.

In order to study the effects of the PEI in planet–disk interactions, we need to first investigate how dust changes the shape of the gap. The vast majority of works on planet–disk interaction assumed either dust-free or gas-free scenarios. Even works with both gas and dust mainly explored the regime of low dust-to-gas ratios, appropriate for primordial disks, but where the dust is mostly passive. The regime appropriate for late-stage transition disks, or debris disks with gas, of a dust-to-gas ratio near unity has been severely unexplored. Here we seek to find at what point the dust backreaction generates a noticeable effect on the gap carving process. We explore the parameter space of dust-to-gas ratios ε = ρ_d/ρ_g and planet-to-star mass ratios q = M_p/M_⋆. Table 1 shows the explored values.

3.1.1. Dust Flattens the Gap Pressure Bump

3.1.2. Dependency on Extent of Mass Reservoir

Because of the global negative pressure gradient, the grains in the outer disk drift inward. At later times, the drift is also affected by the pressure gradient generated by the planetary gap. Notice the strong dependency of the drift speed on the dust-to-gas ratio, with the particles started at ε = 1 drifting slowest as they strongly affect the gas.

An immediate conclusion is that if the drift continues and all grains pile up at the gap edge, any initial dust-to-gas ratio will eventually achieve ε = 1. Indeed, the simulation that started with a global value of ε = 10⁻² (yellow line) is close to this threshold. Whether they can reach this value based on drift alone is determined by the extent of the disk, i.e., by the extent of the mass reservoir. To estimate the disk size needed to achieve ε = 1 at the gap, let us equate the disk mass in grains beyond the gap with the trapped mass. The former is

$$\begin{eqnarray}&&{M}_{\mathrm{disk}}=2\pi {\varepsilon }_{0}{{\rm{\Sigma }}}_{0}{r}_{0}^{p}{\int }_{{r}_{\mathrm{gap}}}^{{r}_{\mathrm{disk}}}{r}^{1-p}{dr},\end{eqnarray} \tag{ 18 }$$

where p ≥ 0 is the slope of the column density gradient, Σ = Σ₀ (r₀/r)^p. If p = 0, then

$$\begin{eqnarray}&&{M}_{\mathrm{disk}}={\varepsilon }_{0}{{\rm{\Sigma }}}_{0}\pi \left({r}_{\mathrm{disk}}^{2}-{r}_{\mathrm{gap}}^{2}\right);\end{eqnarray} \tag{ 19 }$$

if we consider the gap trap to be of width Δr, then the trapped mass is

$$\begin{eqnarray}&&{M}_{\mathrm{gap}}={\varepsilon }_{\mathrm{gap}}{{\rm{\Sigma }}}_{0}2\pi {r}_{\mathrm{gap}}{\rm{\Delta }}r.\end{eqnarray} \tag{ 20 }$$

Now we equate M_disk and M_gap and solve for r_disk,

$$\begin{eqnarray}&&{r}_{\mathrm{disk}}={\left(\displaystyle \frac{{\varepsilon }_{\mathrm{gap}}}{{\varepsilon }_{0}}2{r}_{\mathrm{gap}}{\rm{\Delta }}r+{r}_{\mathrm{gap}}^{2}\right)}^{1/2},\end{eqnarray} \tag{ 21 }$$

requiring ε_gap = 1, and for ε₀ = 10⁻², r_gap = 1.4, and Δr = H, then r_disk ∼ 5.5, i.e., over 500 au. A steep slope will make the demand for a larger disk even more pronounced, as the mass per ring would be less than in the idealized case we considered.

3.1.3. Differences in the Inner Disk for Different Planet Masses and Dust-to-gas Ratios

Another feature of Figures 2 and 3 is the reversal seen in the behavior of the inner disk with dust-to-gas ratio for the two planetary masses. In the outer disk, both planets show a depression in the gap edge as the dust-to-gas ratio increases. However, for the inner disk, while the Jupiter-mass planet presents the same behavior, the Neptune-mass planet does the opposite, with the inner gap wall increasing with dust-to-gas ratio.

Zoom In Zoom Out Reset image size

The difference in behavior seems to be associated with the particle drift. This effect is seen in Figure 4. The red lines show runs with ε = 1, and black lines show runs with ε = 0.1. Solid lines have b = 1 (b is the power law of the temperature, giving a global negative pressure gradient); the dashed lines correspond to b = 0 and thus no global negative pressure. Without the drift, the situations are similar, with ε = 1 (red dashed line) slightly more depressed than ε = 0.1 (black dashed line) because of the flattening effect already discussed. For the runs with the global drift, a pronounced difference is seen. As particles drift inward, to conserve angular momentum, they push the gas outward. For a higher dust-to-gas ratio (solid red line), more gas is pushed outward than for a lower dust-to-gas ratio (black solid line). This effect is not seen in the Jupiter-mass case because in this case, the planet carves a deep gap early enough. In the Neptune case, the shallow gap does not lead to a strong enough pressure gradient to virtually dominate the particle motion, as it does in the Jupiter case.

Zoom In Zoom Out Reset image size

In Figure 5, we show 2D plots of the gas density for the planet masses q = 10⁻³ (top panels) and 10⁻⁴ (bottom panels). Dust-to-gas ratios ε = 10⁻², 10⁻¹, and 1 are shown in the left, middle, and right panels, respectively. The effect of a nonzero dust-to-gas ratio is barely distinguishable for ε = 10⁻² and q = 10⁻⁴ (compare to Figure 1). For a dust-to-gas ratio of 0.1, the gas in the disk becomes more turbulent at the edges of the planet's orbit. The spirals that are caused by the planet in the dustless case are smooth and laminar, whereas the addition of dust to the necessary threshold causes them to become visibly turbulent at the ε = 0.1 case. For ε = 1, the wake and gap look even more turbulent, with the spiral clearly passing from laminar to turbulent at some distance from the planet in a pattern resembling cigarette smoke.

Zoom In Zoom Out Reset image size

We stress that although we consider debris disks and micron grains, this particular result depends only on the Stokes number and dust-to-gas ratio. It is therefore applicable to primordial protoplanetary disks if pebbles and boulders ever reach these high values of ε.

In Figure 6, we plot the dust (top) and gas (bottom) for a disk with H = 0.05 and an initial dust-to-gas ratio of ε = 10⁻¹. For each panel, the left columns are simulations containing an embedded planet and no photoelectric heating, the middle columns are simulations of a disk containing photoelectric heating but no embedded planet, and the right columns are disks containing both an embedded planet and photoelectric heating. The top and bottom rows for each panel represent a Jupiter-mass and a Neptune-mass planet, respectively.

Zoom In Zoom Out Reset image size

The top panels show the local dust density divided by the initial dust density on a logarithmic scale. The cavities carved in runs A and B manifest differently, with run A having a gap noticeably larger than that of run B. Run A also has an overall lower dust density throughout the disk compared to run B. These well-known results (Paardekooper & Mellema 2004) are expected as a higher-mass planet clears its orbit.

Run C demonstrates the effects of photoelectric heating on a disk. The structures created are similar to those of Lyra & Kuchner (2013), with rings and incomplete arcs.

Runs D and E correspond to a disk subject to photoelectric heating and containing a planet. Run E, which represents the Neptune analog, is nearly identical to the control run C, containing the same structures as in the pure instability case. We conclude that the addition of a Neptune-mass planet does not noticeably alter the shape or number of structures from the pure instability case. Run D, which contains a Jupiter-mass perturber, shows larger arcs and a discernible cavity. This is another novel result. When photoelectric heating is included, the PEI obscures structures induced by planets unless the planet's mass is sufficiently large to carve a noticeable gap. The shape of the structures caused by PEI are also changed in the Jupiter case compared to the Neptune case, likely due to the stronger eccentricity stirring at mean-motion resonances with higher planet mass (Murray & Dermott 1999; Zhu et al. 2014).

The bottom panel shows the same plots but for the gas (gas density divided by the initial gas density) on a linear scale. The most pronounced feature is the anticorrelation between the dust and the gas, such that a large concentration of dust coincides with gas depletion. Reiterating for the gas case, we see that a Jupiter-mass perturber carves a larger gap, as well as a more prominent spiral. The spiral launched by the Neptune-mass planet is almost lost in the midst of the structure generated by the PEI.

In Figure 7, we plot an averaged linear dust density as a function of radius from the star. The left panel coincides with a Jupiter-mass planet and the right with a Neptune-mass planet. In the left panel, the combined case (black line) creates larger concentrations of dust than the pure instability case (blue line). Additionally, the more massive planet (left) carves a considerable gap at the location of the planet (100 au). The width of the gap is on the order of 50 au. For the Neptune case, the blue and black lines coincide at various locations throughout the disk. This echoes our conclusion that a planet must carve a substantial gap in the dust for the patterns created by the planet to be differentiated from those associated with photoelectric heating. The shapes shown here in the dust can be used in tandem with actual observations in order to narrow down the possibility of an embedded planet in a disk where the gas is subject to photoelectric heating.

Zoom In Zoom Out Reset image size

In Figure 8, we plot the dust and gas for a disk with a value of h = 0.05 and dust-to-gas ratio of 1. For the dust, runs F and G show a gap due to the perturber's influence, as expected.

Zoom In Zoom Out Reset image size

The increase of ε to unity enhances the effect of the dust, as it is easier now for dust concentrations to affect or dominate the gas motions. The extent of the dust distribution is larger for the ε = 1 case than for ε = 0.1 because the dust dominated over the gas and stopped drifting. The lack of drift is clear from the inner disk, where the amount of dust in the ε = 1 case is substantial compared to the depleted inner disk in the ε = 0.1 case. Run H, which corresponds to the PEI case, shows various arcs, rings, and gaps. This is identical to the fiducial run of Lyra & Kuchner (2013; Figure 3 of that paper). Comparing to the PEI case in Figure 6, we see that the structures are significantly different, with mostly small-scale arcs in the ε = 0.1 case and no complete rings. This is similar to the backreaction-free models of Lyra & Kuchner (2013), who found that backreaction was needed to stabilize and maintain complete rings. In that paper, we examined only ε = 1 for the 2D runs, finding that in a run without the backreaction (Supplementary Figure 6 of that paper), the disk broke into several clumps. We explained it in terms of the tendency of the gas to counterrotate in and around regions of high pressure, which is, in turn, disrupted by the backreaction. For the ε = 0.1 case, the backreaction is weaker than in the ε = 1 case, so the ability of the PEI to structure the dust into rings is diminished. Run C, however, is more self-organized than the no-backreaction run in Lyra & Kuchner (2013), which appears strongly turbulent. The difference between these runs is not just the presence of backreaction but also the lower amount of dust that leads to less heating.

In Figure 9, we check the effect of the dust backreaction for a simulation of ε = 0.1. The top panels show the gas, while the bottom panels show the dust; the left panels are runs without backreaction, and the right panels are runs with backreaction. We see that the backreaction indeed confines the dust into shapely and organized structures. The run without backreaction shows spiral features and accumulation of power in high azimuthal wavenumbers, albeit at a lesser degree than in the ε = 1 run in Lyra & Kuchner (2013). In Figure 10, we show the case for a much lower dust to gas ratio of 0.01. In the top panel for the dust BB, the dust does not arrange into structures despite the backreaction being present. This can be attributed to the lower quantity of dust present. The conclusion is that well-organized PEI structures should be expected only in disks with high dust-to-gas ratios, higher than primordial.

Zoom In Zoom Out Reset image size

Runs I and J show a disk with PEI and planets of q = 10⁻³ and 10⁻⁴, respectively. The more massive perturber clears out a noticeable gap in the dust. Comparing run J to run H, we see no noticeable differences. This echoes the conclusion that only planets that carve substantial gaps can be disentangled from the effect of PEI. For the gas in Figure 8, we see similar results to Figure 6. Runs I and J, which include a planet and PEI, show one important difference. The Lindblad spiral of the Neptune analog, which had already become weaker in run E (ε = 0.1), has all but disappeared in run J (ε = 1), washed out by the stronger effect of the PEI. The Jupiter-mass planet in run I is still visible, showing a gap and the Lindblad spiral.

Zoom In Zoom Out Reset image size

Figure 11 is the same as Figure 7, but now the dust-to-gas ratio is equal to unity. For the Jupiter-mass planet, the dust averages are shown on the left-hand side. The red and black lines both contain a planet and show a depletion of the dust around the location of the planet at r = 1. The blue line is for a simulation with only PEI, and the dust densities are unchanged where r = 1. The Neptune case is plotted on the right. Again, the red and black lines are for simulations with a planet. The pure instability case (blue line) coincides almost identically with the case including a planet and PEI, except for the peak near the planet, inside the gap. In this simulation, a case may be made that with sufficient accuracy, an observation would be able to distinguish between a PEI (blue line) and a planet + PEI (black line) case, judging by the lower dust peak at r = 1 in the planet + PEI simulation, so we do not categorically state that the case is hopeless for these lower-mass planets, only that significant effort from the observational side would be required to detect such a subtle difference. Our results help in the effort to remove the signatures from hydrodynamical effects in order to unambiguously detect planets in disks and may be applicable to disks such as HD 141569A and 49 Ceti, where the structures seen qualitatively match those produced by PEI.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because our results depend on the ability of the planet to carve a gap, the parameter that controls the observability is not the mass of the planet but the thermal mass. The thermal mass is that for which the Hill radius of the planet equals the scale height of the disk. That yields a mass ratio

$$\begin{eqnarray}&&{q}_{\mathrm{th}}=3{h}^{3}.\end{eqnarray} \tag{ 22 }$$

For h = 0.05, a Jupiter-mass planet is equivalent to $q=8{q}_{\mathrm{th}}/3\approx 2.7{q}_{\mathrm{th}}$, and the q = 10⁻⁴ planet is equivalent to ≈0.27q_th.

Increasing the aspect ratio h by a factor of 2 (to h = 0.1) increases the thermal mass by a factor of 8. The Jupiter-mass planet should become equivalent to q_th/3, roughly equivalent to the q = 10⁻⁴ planet in the h = 0.05 case.

We have found the condition that disentangles the carving of a gap from the structures created from the PEI. We take the radius of the dust gap and the frequency of the PEI to be ${r}_{{gap}}^{(d)}$ and λ_PEI, respectively. In order for the two processes to be differentiated, we must have ${r}_{{gap}}^{(d)}\gt {\lambda }_{\mathrm{PEI}}$. In this regime, the gap is large enough that the peaks caused by the instability are within the gap, as seen in Figure 11 for the Jupiter case. Comparing this to the Neptune case in the same figure, ${r}_{{gap}}^{(d)}\approx {\lambda }_{\mathrm{PEI}}$, showing no concrete difference between the PEI and combined case.

In the next subsections, we explore the role of temperature.

3.4.1. Disk with ε = 1 and h = 0.1

Figure 12 shows the dust and gas for a disk of aspect ratio h = 0.1. The top panel plots the dust and the bottom panel plots the gas, as in the previous cases. The Jupiter-mass planet has a thermal mass of 0.33, and the Neptune-mass planet has 0.03. We see that indeed, the Neptune analog does not carve a gap in the gas, staying deeply embedded. The Jupiter-mass planet carves a shallower gap, as expected from its thermal mass now being lower than 1.

The dust panel of run L shows a departure from axisymmetry and an enhanced dust drift compared to the Jupiter case for otherwise the same parameters (run K). Examining the gas plots below, we notice three prominent structures in run L, at r = 2. These structures resemble vortices, but if so, it is curious that they are not trapping dust. This should only happen if the vortices are triggered after the dust has drifted inward, so that there is no dust to be captured. Indeed, the outer boundary of the dust distribution, seen in the dust panel of run L, is inward of the radial location of the vortices, as seen in the gas panel of the same run.

It is suspicious that this phenomenon occurs near the boundary. Yet run K does not show the vortices, which rules out boundary conditions as the culprit. Instead, we conjecture that what we are seeing is RWI triggered by the dust drift in a high dust-to-gas ratio. As the dust drifts inward, the dust front constitutes a sharp discontinuity. At a high dust-to-gas ratio, the backreaction of the drag force has significant dynamical relevance, so the gas has a different rotational velocity immediately inward and outward of the dust front. This difference eventually becomes sufficient to trigger the RWI. This is a novel effect, RWI triggered at a dust front. While in this case, the dust front is artificial (it would not exist if there was a supply of dust from the outer disk), there are physical situations where a sharp dust discontinuity exists, as in, e.g., the birth ring model. In disks with a birth ring and significant gas, we thus predict that RWI should be triggered in the gas. If the vortices disturb the population of grains in the same way, the birth ring would turn eccentric, as the dust front in run L.

The Jupiter-mass planet does not show the same phenomenon because of the reduced dust drift, so the front is still located at the boundary region of the run, where the quantities are driven back to the initial conditions. For the h = 0.05 run, we conjecture that the same process would be present, but the dust has not drifted far enough inward yet for the dust front to be resolved out of the boundary region. The dust drift is proportional to the sub-Keplerian reduction parameter η, which is proportional to h². Increasing h by a factor of 2 increases the drifting time by a factor of 4, making the effect noticeable. Also, the resolution was probably too low to resolve the narrow band of the instability. Increasing h to 0.1 (Figure 13) doubled the width of the transition, allowing the modes to be resolved.

Zoom In Zoom Out Reset image size

In run M, which is a simulation with PEI and no planet, we see less arcs and more rings than in the h = 0.05 case; also, the structures are thicker with the increased value of H. The gas disk is empty in the center, owing to increased heating, according to Equation (10). The gas expands due to the higher temperature, thinning out. At the boundaries, it retains the high initial density value because of our boundary condition that drives the quantities back to their initial values.

The combined case of a planet and PEI, shown in runs N and O, shows rings and arcs as well. Run N shows a larger gap, with a lack of dust in the immediate area, apart from two concentrations at 12 and 5 o'clock. Run O, for a lower-mass planet, shows various arcs and a nondiscernible gap in the dust.

3.4.2. Disk with ε = 1 and h = 0.2

Finally, we consider a disk with a scale height of h = 0.2, which is on the hotter side for primordial disks but fairly typical for the scale heights of debris disks with gas, including β Pic. In the top panel of Figure 14, we plot the dust structures, and in the bottom panel, we plot the gas structures. For the Jupiter case (top left), we see a gap carved in the dust. However, for the Neptune case, we only see a small ring inward of the planet and nothing outward. Clearly, the dust drift is more intense for the Neptune case. The reason for this is twofold. First, as seen in the gas plots (bottom panel), the Jupiter-mass planet, though having a thermal mass of 0.04, carves a very shallow gap. The Neptune-mass planet has a thermal mass of 0.004 and is unable to open a gap. The high density outside the Jupiter-mass planet traps the drifting grains. In the Neptune-mass planet case, they stream past the planet orbit and are trapped at the density maximum inward of it.

Zoom In Zoom Out Reset image size

It is curious that although there is no gas gap, the grains are able to stream past the planet, bypassing the dust filtering that planets usually do. This seems to happen because the drifting time is faster than a planetary orbit. We can compare the time it takes a grain to drift past the corotational region of the planet to the orbital period of the planet,

$$\begin{eqnarray}&&{\rm{\Theta }}=\displaystyle \frac{{t}_{\mathrm{drift}}}{{t}_{\mathrm{orb}}}.\end{eqnarray} \tag{ 23 }$$

If Θ > 1, the drifting is slow and the planet is able to scatter the particle. If, on the other hand, Θ < 1, the planet is too slow and the grain can drift past it. Consider the drift velocity of a grain (Youdin 2010),

$$\begin{eqnarray}&&{v}_{r}=-\displaystyle \frac{2\mathrm{St}\,\eta \,{v}_{k}}{1+{\mathrm{St}}^{2}},\end{eqnarray} \tag{ 24 }$$

where v_k = Ω_k r is the Keplerian velocity and

$$\begin{eqnarray}&&\eta \equiv =-\displaystyle \frac{1}{2\rho {{\rm{\Omega }}}^{2}r}\displaystyle \frac{\partial p}{\partial r}\end{eqnarray} \tag{ 25 }$$

is the sub-Keplerian parameter. The time it takes to drift is ${v}_{r}/2{x}_{s}$, where x_s is the half-width of the corotational region, given by Paardekooper & Papaloizou (2009),

$$\begin{eqnarray}&&{x}_{s}^{2}=1.68{r}_{p}^{2}\displaystyle \frac{q}{h},\end{eqnarray} \tag{ 26 }$$

and r_p is the planet's location. Solving for Θ,

$$\begin{eqnarray}&&{\rm{\Theta }}=\displaystyle \frac{{\rm{\Omega }}{x}_{s}}{\pi {v}_{r}}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\approx 0.4\ {q}^{1/2}{h}^{-5/2}{\chi }^{-1}{\mathrm{St}}^{-1}\left(1+{\mathrm{St}}^{2}\right),\end{eqnarray} \tag{ 28 }$$

where χ = α + β with $\alpha =-\partial \mathrm{ln}\rho /\partial \mathrm{ln}r$ and $\beta =-\partial \mathrm{ln}T/\partial \mathrm{ln}r$. For q = 10⁻⁴, h = 0.2, and $\mathrm{St}=1$ and considering α + β = 1, we get Θ = 0.5. For the Jupiter-mass planet, Θ = 1.4. Thus, it is reasonable to expect that in the Neptune-mass case, the $\mathrm{St}=1$ grains can drift past the corotational region in this hot disk but that the Jupiter-mass planet, with a larger corotational region, will be able to scatter the grains and carve a deep dust gap.⁹ In the simulation with PEI (run R), we see a much denser ring of dust in the center of the domain. The combined case of instability plus planet (runs S and T) are shown in the plots on the right-hand side. For the massive planet (run S), the main differences are that the grains, instead of being scattered roughly uniformly as in run P, now concentrate in arcs. In addition, an arc coinciding with the location of the Lagrange point L5 is prominent. The lower-mass planet (run T) seems indistinguishable from the pure instability case (run R).

The behavior of the gas (bottom panels) resembles the situation with h = 0.1. The gas density is much lower throughout most of the disk due to the enhanced photoelectric heating being driven back to the initial value at the boundaries. A spiral is barely discernible for the Jupiter-mass planet. No clear distinction is present for the lower-mass planet.