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

TL;DR

This study uses hydrodynamical models to distinguish between planetary influences and photoelectric instability in shaping structures within gas-rich, optically thin circumstellar disks, revealing conditions under which planetary effects can be identified.

Contribution

It demonstrates how PEI can obscure planetary signatures in disk structures and identifies criteria for detecting planets despite PEI effects.

Findings

PEI can mask planetary structures unless the planet's gap is sufficiently large.

High dust-to-gas ratios lead to dust accumulation and Rossby wave instability.

Planetary gaps wider than PEI wavelength reveal planetary presence.

Abstract

Structures in circumstellar disks such as gaps and rings are often attributed to planets. This connection has been difficult to show unequivocally, as other processes may also produce these features. Particularly, a photoelectric instability (PEI) has been proposed, operating in gas-rich optically thin disks, that generates structures predicted by planet--disk interactions. We examine the question of how to disentangle planetary effects on disk structure from the effects of the PEI. We use the Pencil Code to perform 2D global hydrodynamical models of the dynamics of gas and dust in a thin disk, with and without planetary perturbers. Photoelectric heating is modeled with an equation of state where pressure is proportional to dust surface density. The drag force on grains and its backreaction on the gas are included. Analyzing the situation without PEI, we find that gas--dust interactions alter the shape of the planetary gap from the dust-free case when the local dust-to-gas ratio $\varepsilon$ approaches unity. This result applies also to primordial disks, because dust drifting inwards accumulates at the edge of the planetary gap, and any initial dust-to-gas ratio eventually achieves $\varepsilon=1$ if the dust reservoir is sufficient. We find a result particular to high dust-to-gas ratio disks as well: as dust drifts inwards, the dust front becomes a sharp transition, and the backreaction triggers the Rossby wave instability. When PEI is included, we find that it obscures structures induced by planets unless the planet's mass is sufficiently large to carve a noticeable gap. Specifically, the instability generates arcs and rings of regular spacing: a planet is discernible when it carves a dust gap wider than the wavelength of the PEI.

Tables

Table 1. Table 1 Parameter space explored in our simulations: h ℎ h is the disk scale height of the disk, ε 𝜀 \varepsilon is the initial dust-to-gas ratio, and q 𝑞 q is the planet-to-star mass ratio. The letter identifier in the first column is shown in the upper left corner of each plot in this paper, for easier referral.

Run	h	$\epsilon$	$q$
Fiducial runs with no dust
i	0.05	$0$	${10}^{-3}$
ii	0.05	$0$	${10}^{-4}$
Runs with a single planet.
A	0.05	${10}^{-1}$	${10}^{-3}$
B	0.05	${10}^{-1}$	${10}^{-4}$
F	0.05	$1$	${10}^{-3}$
G	0.05	$1$	${10}^{-4}$
K	0.1	$1$	${10}^{-3}$
L	0.1	$1$	${10}^{-4}$
P	0.2	$1$	${10}^{-3}$
Q	0.2	$1$	${10}^{-4}$
AA	0.05	${10}^{-4}$	${10}^{-3}$
AB	0.05	${10}^{-4}$	${10}^{-4}$
AC	0.05	${10}^{-3}$	${10}^{-3}$
AD	0.05	${10}^{-3}$	${10}^{-4}$
AE	0.05	${10}^{-2}$	${10}^{-3}$
AF	0.05	${10}^{-2}$	${10}^{-4}$
Runs with photoelectric instability.
C	0.05	${10}^{-1}$	0
H	0.05	$1$	0
M	0.1	$1$	0
R	0.2	$1$	0
Runs with a planet and instability.
D	0.05	${10}^{-1}$	${10}^{-3}$
E	0.05	${10}^{-1}$	${10}^{-4}$
I	0.05	$1$	${10}^{-3}$
J	0.05	$1$	${10}^{-4}$
N	0.1	$1$	${10}^{-3}$
O	0.1	$1$	${10}^{-4}$
S	0.2	$1$	${10}^{-3}$
T	0.2	$1$	${10}^{-4}$
Runs with no gas backreaction.
BA	0.05	${10}^{-1}$	0
BB	0.05	${10}^{-2}$	0

Equations

$$\frac{\partial {\Sigma }_g}{\partial t}$$

$$\frac{\partial \mathbf{u}}{\partial t}$$

$$\frac{\partial S}{\partial t}$$

$$\Phi$$

$$P$$

$$\frac{d\mathbf{x}}{dt}$$

$$\frac{d\mathbf{v}}{dt}$$

$$S=c_{{}_V}\left[\ln \left(\frac{P}{P_0}\right)-\gamma \ln \left(\frac{{\Sigma }_g}{{\Sigma }_{g0}}\right)\right]$$

$$T=\frac{c_s^2}{c_p\left(\gamma -1\right)}$$

$$T_p=T_0{\left(\frac{{\Sigma }_d}{{\Sigma }_{g0}}\right)}^{\beta }$$

$${\tau }_{\mathrm{f}}={\tau }_{\mathrm{f}}{}_0\left(\frac{{\Omega }_0}{\Omega }\right)$$

$$\frac{D}{Dt}=\frac{\partial }{\partial t}+\mathbf{u}\cdot \nabla$$

$$\underset{{\tau }_{{}_T}\to 0}{\lim }P=c_{{}_V}{T}_0(\gamma -1)\left(\frac{{\Sigma }_g}{{\Sigma }_{g0}}\right){\Sigma }_d$$

$$G{M}_{\star }={\Omega }_0=r_0=1$$

$$\mathrm{St}\equiv {\tau }_{\mathrm{f}}\Omega =1$$

$${\tau }_{\mathrm{f}}=\frac{a_{\bullet }{\rho }_{\bullet }}{\rho c_s},$$

$$\rho {|}_{\mathrm{St}=1}=\frac{a_{\bullet }{\rho }_{\bullet }}{c_s}$$

$$M_{\mathrm{disk}}=2\pi {\epsilon }_0{\Sigma }_0{r}_0^p{\int }_{r_{\mathrm{gap}}}^{r_{\mathrm{disk}}}{r}^{1-p}dr$$

$$M_{\mathrm{disk}}={\epsilon }_0{\Sigma }_0\pi \left(r_{\mathrm{disk}}^2-r_{\mathrm{gap}}^2\right).$$

$$M_{\mathrm{gap}}={\epsilon }_{\mathrm{gap}}{\Sigma }_{0}2\pi r_{\mathrm{gap}}\Delta r.$$

$$r_{\mathrm{disk}}={\left(\frac{{\epsilon }_{\mathrm{gap}}}{{\epsilon }_0}2r_{\mathrm{gap}}\Delta r+r_{\mathrm{gap}}^2\right)}^{1/2}$$

$$q_{\mathrm{th}}=3h^3$$

$$\Theta =\frac{t_{\mathrm{drift}}}{t_{\mathrm{orb}}}$$

$$v_r=-\frac{2\mathrm{St}\eta v_k}{1+{\mathrm{St}}^2},$$

$$\eta \equiv =-\frac{1}{2\rho {\Omega }^{2}r}\frac{\partial p}{\partial r}$$

$$x_s^2=1.68r_p^2\frac{q}{h}$$

$$\Theta$$

$$\approx$$

Full text

1. Introduction

Exoplanets form and dwell in disks of gas and dust, which show rings, spirals, and other patterns in optical (Grady et al., 2013; Konishi et al., 2016; Hornbeck et al., 2016), infrared (Currie & Kenyon, 2009; Benisty et al., 2015; Currie et al., 2014, 2016), and sub-millimeter images (Andrews et al., 2011; Dong et al., 2017; Muto et al., 2012; Dong et al., 2015; Dong & Dawson, 2016). Planets orbiting in disks perturb them gravitationally, creating similar patterns to those observed (Kley & Nelson, 2012) that have in turn been used to infer the presence of hidden planets and constrain their properties (e.g. Kuchner & Holman, 2003; Chiang et al., 2009; Dong et al., 2016; Dong & Fung, 2017; Dipierro & Laibe, 2017; Hord et al., 2017). A planet has been observed within the gap of the transition disk around PDS 70 (Müller et al., 2018; Keppler et al., 2018), a disk that indeed shows several of the features predicted by planet–disk interaction, validating the approach of using disk features as signposts for planets.

Yet, our understanding of disk processes and planet–disk interaction is not complete. The dynamics of planet–disk interaction are usually modeled in the dust-free or gas-free limits. The former (e.g. Kley & Nelson, 2012) is applicable to primordial disks, where gas dominates the dynamics, entraining the dust grains seen in optical and infrared images. The latter (e.g. Kuchner & Holman, 2003) is applicable to debris disks, where dust dominates the dynamics. But as disks evolve from gas-rich toward sparse, gas-free systems like our own solar system, they inevitably pass through a phase when they become optically thin, and the gas-to-dust mass ratio drops to the range of 0.1–1. In this phase, the disk physics becomes more complex, because dust grains and gas molecules both interact with stellar photons, and both gas and dust play an important role dynamically due to their similar total masses. (Klahr & Lin, 2005; Besla & Wu, 2007). This phase corresponds to young debris disks and to certain regions of every transitional and pre-transitional disk.

Over the past two decades, studies have shown that gas indeed is a common component of many debris disks (e.g. Chen & Jura, 2003; Redfield, 2007; France et al., 2007; Roberge & Weinberger, 2008; Moór et al., 2011; Welsh & Montgomery, 2013; Riviere-Marichalar et al., 2014; Moór et al., 2015; Lieman-Sifry et al., 2016; Marino et al., 2016; Hales et al., 2017; Marino et al., 2017). The bulk of the gas may be metals, or it may be molecular hydrogen, which is much more difficult to detect (Thi et al., 2001). The origin of the gas in these disks is a matter of continued debate. It is either leftover primordial gas, or it is produced from solid debris by photodesorption (Grigorieva et al., 2007) or collisions (Czechowski & Mann, 2007), processes that should occur in every debris disk at some level.

We have only just begun to understand which processes dominate in gas-rich optically thin disks, the nature of the hydromagnetic turbulence (Kral & Latter, 2016), the instabilities that can arise, and how these processes can affect the morphology of these disks. As a matter of fact, gap carving by a planet, a problem long benchmarked (de Val-Borro et al., 2006), has never been studied in the regime of dust-to-gas ratio close to unity. Gap carving is relatively well understood for gas-rich primordial disks, for which the dust back reaction can usually (but not always) be ignored as the dust-to-gas ratio is close to the interstellar value ($\varepsilon\approx 10^{-2}$).

In this case, gas is cleared away by planetary torques, and viscous torques tend to fill the gap (Kley & Nelson, 2012, and references therein). It is also relatively well understood for gas-free debris disks, as objects are cleared away by eccentricity pumping at resonances, leaving the planet neighborhood following the so-called Tisserand tails (e.g. Zhu et al., 2014). Yet, when both gas and dust are dominant, gap carving may deviate significantly from these idealized cases. This leads to the first question we pose: what happens to a disk gap when dust is added to the point that the gas drag backreaction changes the torques acting to close the gap?

A second question we pose, that depends on understanding the first, is related to the combined effect of the gravitational potential of a planet and a dynamical instability in gas-rich optically thin disks. The instability, suggested by Klahr & Lin (2005); Besla & Wu (2007), is affected by the combined operation of two dynamical processes. First, if the disk is optically thin, stellar UV photons eject electrons from dust grains, heating the gas (Kamp & van Zadelhoff, 2004; Jonkheid et al., 2004; Zagorovsky et al., 2010). Second, the hot gas provides a local pressure maximum that concentrates the dust grains (Whipple, 1972), establishing a positive feedback process. In just a few orbital periods, this photoelectric instability (PEI) can radically reshape the disk structure. In the first 2D and 3D models of this instability, Lyra & Kuchner (2013) found that the PEI would concentrate power in high wavenumbers if not for the effect of the gas drag backreaction, which organizes the dust in rings and arcs. They also found that the PEI leads to several patterns usually associated with planet–disk interactions, challenging any correlations drawn between dust structures and perturbing planets. Richert et al. (2018) incorporate radiation pressure into a PEI model, and show that while grains of different sizes feel very different amounts of different radiation pressure force, the PEI can still produce dramatic structures in such a more realistic disk model, depending on the parameters of the system. The morphology of the flow depends strongly on the gas density (which determines the gas drag), the dust density (which determines the photoelectric heating), and the dust-to-gas ratio (which determines the drag backreaction). For disks with modest level of gas $\approx 10^{-6}\ {\rm g\,cm^{-2}}$ and dust-to-gas ratio $\varepsilon=10^{-1}$, sharp concentric rings and arcs are seen, showing the anticorrelation between gas and dust seen in models of the PEI without radiation pressure (Lyra & Kuchner, 2013). For higher surface density $\approx 10^{-4}\ {\rm g\,cm^{-2}}$ and lower gas-to-dust ratio $\varepsilon=10^{-3}$, erratic spiral structure is seen throughout the disk, with transient small-scale structure in dynamical timescales, and no obvious tendency of the dust and gas to anticorrelate. At even higher gas densities, outward gas flows can entrain even larger grains that do not participate in the instability.

Given the variety of structures that the PEI can generate in a disk, which muddles an unequivocal association with planetary perturbations, a natural question to ask is: what are the expected patterns when both a planet and PEI are present in a disk? How do we disentangle them? This is the goal of the current study.

We structure this paper as follows. In Sect. 2 we describe the equations and numerical method used. In Sect. 3 we present our results, concluding in Sect. 5.

2. Method

2.1. The Equations

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

[TABLE]

In the equations above, $\varSigma_{g}$ and $\varSigma_{d}$ are the vertically integrated gas density and bulk density of dust, respectively. Their quotient defines the dust-to-gas ratio $\varepsilon=\varSigma_{d}/\varSigma_{g}$. 111As we use the 2D infinitely thin approximation, the different scale heights of grains and gas is not resolved. The definition of dust-to-gas ratio used here is the ratio of their column densities; in 3D it should be the ratio of their volume densities, which, because of the different degree of stratification, is larger for loosely coupled grains in the midplane, by $\exp{[(\mathrm{St}+\delta)/\delta}]$ (Lyra & Lin, 2013). For well-sedimented grains, the 2D approach may underestimate the local dust-to-gas ratio by orders of magnitude. The velocity of gas parcels is given by ${\boldsymbol{u}}$ ; $P$ stands for the vertically integrated pressure, $\varPhi$ 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 a N-body evolution dynamically integrated with the hydrodynamics evolution. In Eq. (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 Eq. (3), $S$ is the gas entropy, given by

[TABLE]

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_{0}$ and $\varSigma_{g0}$ are reference pressure and gas density values, respectively. In Eq. (3) $T$ stands for the gas temperature

[TABLE]

where $c_{s}$ is the sound speed; $\tau_{{}_{T}}$ is the thermal coupling times 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 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)

[TABLE]

where $T_{0}$ is a reference temperature. We keep $\beta=1$ (see the appendix of Klahr & Lin 2005 and Fig 1 of that work for the estimate of $\beta\approx 1$ given by a balance between photoelectric heating off silicate grains and cooling via the CII line.) In the Lagrangian dust equations, ${\boldsymbol{v}}_{p}$ is the velocity of a dust grain, and ${\boldsymbol{x}}_{p}$ is their position. The last term in Eq. (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_{\mathrm{f}}$, the friction time. The last term in Eq. (2) is the backreaction of the drag force. We choose the friction time as proportional to the local orbital period

[TABLE]

as in this case the relative timescales (the PEI growth rate also scales with the orbital period) are similar at all radii222The drag time in the Epstein regime is (cf. Youdin, 2010) $\tau_{f}=a_{\bullet}\rho_{\bullet}/(c_{s}\rho)$, where $a_{\bullet}$ is the grain radius, $\rho_{m}$ is the grain’s material density, and $\rho$ the gas volume density. As such, Eq. (11) is equivalent to grains of identical size and material density scattered through a disk where $c_{s}\rho\propto\varOmega$. This would be the case for, e.g., $\rho\propto 1/r$ and $c_{s}\propto 1/\sqrt{r}$, which are not unreasonable choices for the density and temperature gradients.. We keep $\tau_{\mathrm{f}}{{}_{0}}\varOmega_{0}=1$. The operator

[TABLE]

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, and equate $T=T_{p}$ and update the sound speed is updated accordingly. This change effectively amounts to choosing a new equation of state that depends on the dust density

[TABLE]

We solve the equations with the Pencil Code333The code, including improvements done for the present work, is publicly available under a GNU open source license and can be downloaded at http://www.nordita.org/software/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).

2.2. Run Parameters and Units

The simulations are performed on a uniform cylindrical grid (coordinates $r,\phi$) with resolution $(N_{r},N_{\phi})=(528,528)$ and linear domain $L_{r}$=[0.4,2.5] and full 2$\pi$ coverage in azimuth. To represent the dust grains, 200,000 Lagrangian superparticles are scattered uniformly across the disk. The star and the 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 ratio of the planet to the host star are $q=10^{-3},10^{-4},10^{-5}$. For each planet we simulated 5 different dust-to-gas ratios, from $10^{-4}$ to 1, logarithmically spaced. The runs from $\varepsilon=10^{-2}$ to $\varepsilon=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

[TABLE]

All quantities with subscript “0” are measured at $r_{0}$. For $r_{0}=100\,AU$ and $M_{\star}=1M_{\odot}$, our unit of time is $\approx$ 160 yr, and the unity of velocity $\approx$ 3 km s*-1*. We investigate disks of aspect ratios $h\in\{0.05,0.1,0.2\}$, where $h=H/r$, with $H=c_{s}/\varOmega$ the disk scale height. For a solar mixture of molecular hydrogen and helium, of mean molecular weight $\mu=2.5$ and $\gamma=7/5$, these correspond to temperatures of $T$ = 5, 20, and 80 K, respectively. We consider monodispersive particles of radius $a_{\bullet}=1\mu$m. This choice, along with the choice of Stokes number unity

[TABLE]

sets the unit of density 444The Stokes number should in fact be dynamic, updated by density and temperature as given by Eq. (16). Setting it as a constant is a choice we make, to better isolate the physical phenomenon we want to study. Yet, notice that because the photoelectric heating introduces an anti-correlation between gas density and gas temperature, the product $c_{s}\rho$ varies a lot less than either $c_{s}$ or $\rho$, so $\mathrm{St}\approx{\rm const}$ is not a bad first approximation).. This is because the friction time in the Epstein regime is

[TABLE]

where $\rho_{\bullet}$ is the grain internal density. For icy grains, it should be $\rho_{\bullet}$ $\approx 1\,$g cm*-3*. Multiplying Eq. (16) by $\varOmega$ to obtain the Stokes number and setting to one, yields, solving for density

[TABLE]

which results, for $h=0.1$, in $\rho\approx 10^{-8}$ g cm*-3*, with associated column density of $\varSigma=\sqrt{2\pi}\rho H\approx 2.5\times 10^{-4}$ g cm*-2*. We keep the gas and dust densities constant with radius, whereas the temperature decreases according to $T=T_{0}r^{-b}$. We usually use $b=1$ corresponding to a negative global pressure gradient, causing the dust to slowly drift inwards. 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 use $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 photoelectric instability we ran 36 simulations. We used two planets, analogous to a Neptune and a Jupiter mass planet. For each planet we varied 3 disk temperatures, 2 dust-to-gas ratios, and 3 different cases for carving a gap in the disk. The cases are: the effect of a single planet with no photoelectric instability, a disk with the instability on its own, and a combination of an embedded planet and the instability.

3. Results

We first show the reproduction of the benchmarks of de Val-Borro et al. (2006) for embedded planets of mass ratio $q=10^{-3}$ (Jupiter) and $q=10^{-4}$ (Neptune) in a disk with no dust. In the upper plots of 1, we plot the gas densities. The Jupiter-mass planet carves out a deep gap in its orbit. Strong vortices are also seen in the inner and outer edges of the gap, spawned by the Rossby wave instability (Lovelace et al., 1999; de Val-Borro et al., 2007). The Neptune-mass case has the same qualitative structures, but much weaker in intensity, and the gap is shallower and narrower. In the bottom plots of 2, we plot the density divided by the initial density as a function of distance from the star. The usual gap shapes are reproduced. The runs with dust and photoelectric instability will be compared against these benchmark gap profiles.

3.1. The effect of dust

In order to study the effects of the photoelectric instability 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 dust explored mainly the regime of low dust-to-gas ratio, appropriate for primordial disks, but where the dust is mostly passive. The regime appropriate for late stage transition disks, or debris disk with gas, of dust-to-gas ratio near unity, has been severely unexplored. Here we seek to find at what point does the dust backreaction generate a noticeable effect on the gap carving process. We explore the parameter space of dust-to-gas ratios $\varepsilon=\rho_{d}/\rho_{g}$ and planet-to-star mass ratios $q=M_{p}/M_{\star}$. Table 1 shows the explored values.

3.1.1 Dust flattens the gap pressure bump

In 2, in the upper panel, we plot the gas density over the initial gas density as a function of distance from the star. The curves correspond to dust-to-gas ratios of 0 (gas-only), $10^{-4}$, $10^{-3}$, $10^{-2}$, $10^{-1}$, and 1, for $h=0.05$. We see that up to $\varepsilon=10^{-2}$ the shape of the gap does not deviate much from the gas-only case. However, for $\varepsilon=10^{-1}$ and more strongly for $\varepsilon=1$, a significant deviation is seen. This corresponds to the regime where the dust begins having a significant effect on the gas. In the lower panel we plot the local dust-to-gas density. The dashed line marks $\varepsilon=1$. We see that only the disks with dust-to-gas ratio of 1 or 0.1 cross this threshold. We conclude that only when the local value of $\varepsilon$ reaches unity the shape of the gap is affected, as one would indeed intuitively expect.

3.1.2 Dependency on extent of mass reservoir

Because of the global negative pressure gradient, the grains in the outer disk drift inwards. 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 $\varepsilon=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 $\varepsilon=1$. Indeed the simulation that started with a global value of $\varepsilon=10^{-2}$ (yellow line) is close to this threshold. If 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 $\varepsilon=1$ at the gap, let us equate the disk mass in grains beyond the gap with the trapped mass. The former is

[TABLE]

where $p\geq 0$ is the slope of the column density gradient, $\varSigma=\varSigma_{0}(r_{0}/r)^{p}$. If $p$=0 then

[TABLE]

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

[TABLE]

Now we equate $M_{\rm disk}$ and $M_{\rm gap}$ and solve for $r_{\rm disk}$

[TABLE]

requiring $\varepsilon_{\rm gap}=1$, and for ${\varepsilon_{0}}=10^{-2}$, $r_{\rm gap}=1.4$, and $\Delta r=H$, then $r_{\rm disk}\sim 5.5$, i.e., over 500 AU. A steep slope will make the demand of 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 ratio

Another feature of 1 and 2 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, the inner gap wall increasing with dust-to-gas ratio.

The difference in behavior seems to be associated with the particle drift. This effect is seen in 4. The red lines show runs with $\varepsilon=1$, and black lines show runs with $\varepsilon=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 $\varepsilon=1$ (red dashed line) slightly more depressed than $\varepsilon=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 particle drift inward, to conserve angular momentum they push the gas outwards. For higher dust-to-gas ratio (solid red line), more gas is pushed outwards than for 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.

In 5 we show 2D plots of the gas density for the planet masses $q=10^{-3}$ (upper panels) and $q=10^{-4}$ (lower panels). Dust-to-gas ratios $\varepsilon=10^{-2}$, $\varepsilon=10^{-1}$, and $\varepsilon=1$ are shown in the left and right panels, respectively. The effect of a non-zero dust-to-gas ratio is barely distinguishable for $\varepsilon=10^{-2}$ and $q=10^{-4}$ (compare to 1). For dust-to-gas ratio 0.1, the gas in the disk becomes more turbulent at 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 $\varepsilon=0.1$ case. For $\varepsilon=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.

We stress that although we consider debris disks and $\mu m$ 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 $\varepsilon$.

3.2. Disk with $\varepsilon=10^{-1}$ and $h=0.05$

In 6, we plot the dust (top) and gas (bottom) for a disk with $H=0.05$ and an initial dust-to-gas ratio of $\varepsilon=10^{-1}$. 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 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.

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 discernable cavity. This is another novel result. When photoelectric heating is included, the photoelectric instability 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 photoelectric instability are also changed in the Jupiter case compared to the Neptune, 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 photoelectric instability.

In Figure 6, we plot an averaged linear dust density as a function of radius from the star. The left panel coincides with a Jupiter-mass and the right to a Neptune-mass planet. On the left panel the combined case in solid black, creates larger concentrations of dust than the pure instability case in solid blue. 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.

3.3. Disk with $\varepsilon=1$ and $h=0.05$

In 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.

The increase of $\varepsilon$ 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 $\varepsilon=1$ case than for $\varepsilon=0.1$ because the dust dominated over the gas and stopped drifting. The lack of drift is clear from the inner disk, where we amount of dust in the $\varepsilon=1$ is substantial compared to the depleted inner disk in the $\varepsilon=0.1$ case. Run H, which corresponds to the photoelectric instability 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 6, we see that the structures are significantly different, with mostly small-scale arcs in the $\varepsilon=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 $\varepsilon=1$ for the 2D runs, finding that in a run without the backreaction (Suppl. Fig 6 of that paper), the disk broke into several clumps. We explained it in terms of the tendency of the gas to counter-rotate in and around regions of high pressure, which is turn is disrupted by the backreaction. For the $\varepsilon=0.1$ case, the backreaction is weaker than in the $\varepsilon=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), that appears strongly turbulent. The differences between these runs is not just the presence of backreaction, but also the lower amount of dust, that leads to less heating.

In 9, we check the effect of the dust backreaction for a simulation of $\varepsilon=0.1$. The upper panels show the gas while the lower panels show the dust; the left panels are runs without backreaction, the right panels 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 $\varepsilon=1$ run in Lyra & Kuchner (2013). We expect that for $\varepsilon=10^{-2}$ the backreaction will be less important and thus the structures will be more turbulent. The conclusion is that well-organized PEI structures should be expected only in disks with high dust-to-gas ratio, higher than primordial.

Runs I and J show a disk with PEI and a planet, of $q=10^{-3}$ and $q=10^{-4}$, 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 8, we see similar results to 6. Runs I & 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 ($\varepsilon=0.1$), has all but disappeared in run J ($\varepsilon=1$), washed out by the stronger effect of the PEI. The Jupiter-mass planet in run I is still visible in the gap showing a gap and the Lindblad spiral.

11 is the same as 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 dashed red line and the solid black line both contain a planet, and show a depletion of the dust around the location of the planet at $r=1$. The solid blue line is for a simulation with only PEI, and dust densities are unchanged where $r=1$. The Neptune case is plotted on the right. Again, the red dashed line and solid black line are for simulations with a planet. The pure instability case, in the solid blue line, coincides almost identically to 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 subtle a difference. Our results help on 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 produces by PEI.

3.4. Dependence on thermal mass

Because our results depends 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

[TABLE]

For $h=0.05$, a Jupiter-mass planet is equivalent to $q=8q_{\rm th}/3\approx 2.7q_{\rm th}$, and the $q=10^{-4}$ planet to $\approx 0.27q_{\rm th}$.

Increasing the aspect ratio $h$ by a factor 2 (to $h=0.1$) increases the thermal mass by a factor 8. The Jupiter-mass planet should become equivalent to $q_{\rm th}/3$, roughly equivalent to the $q=10^{-4}$ 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 photoelectric instability. We take the radius of the dust gap and the frequency of the photoelectric instability to be $r_{gap}^{(d)}$ and $\lambda_{\rm PEI}$, respectively. In order for the two processes to be differentiated, we must have that $r_{gap}^{(d)}>\lambda_{\rm PEI}$. In this regime, the gap is large enough that peaks caused by the instability are within the gap, as is seen in 13 for the Jupiter case. Comparing this to the Neptune case in the same figure $r_{gap}^{(d)}\approx\lambda_{\rm 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 $\varepsilon=1$ and $H=0.1$

12 shows the dust and gas for a disk of aspect ratio $h=0.1$. The top panel plots the dust and the bottom panels the gas, as in the previous cases. The Jupiter mass planet has thermal mass $0.33$, and the Neptune mass planet $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 departure from axisymmetry and an enhanced dust drift compared to the Jupiter case for elsewise 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 inwards, 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 inwards 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 culprit. Instead, we conjecture that what we are seeing is Rossby Wave Instability triggered by the dust drift in a high dust-to-gas ratio. As the dust drifts inwards, the dust front constitutes a sharp discontinuity. At high dust-to-gas ratio, the backreaction of the drag force has significant dynamical relevance, so the gas has a different rotational velocity immediately inwards and immediately 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 inwards yet for the dust front to be resolved out of the boundary region. The dust drift is proportional to the sub-Keplerian reduction parameter $\eta$, which is proportional to $h^{2}$. Increasing h by a factor 2 increasing the drifting time by a factor 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 doubled the width of the transition, allowing the modes to be resolved.

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, the structures are thicker also with the increased value of $H$. The gas disk is empty in the center, owing to an increased heating, according to Eq. (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 non-discernable gap in the dust.

3.4.2 Disk with $\varepsilon=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 $\beta$ Pic. We plot in 14, on the top panel, the dust structures, and in the lower panel the gas structures. For the Jupiter case (upper 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 outwards. Clearly the dust drift is more intense for the Neptune case. The reason for this is twofold. Firstly, as seen in the gas plots (lower panel), the Jupiter-mass planet, though having a thermal mass of 0.04, carves a very shallow gap. The Neptune mass planet has thermal mass 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 inwards of it.

It is curious that although there is not a 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 planet to the orbital period of the planet

[TABLE]

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

[TABLE]

where $v_{k}=\varOmega_{k}r$ is the Keplerian velocity and

[TABLE]

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

[TABLE]

and $r_{p}$ the planet’s location. Solving for $\Theta$,

[TABLE]

where $\chi=\alpha+\beta$ with $\alpha=-\partial\ln\rho/\partial\ln r$ and $\beta=-\partial\ln T/\partial\ln r$. For $q=10^{-4}$, $h=0.2$, $\mathrm{St}=1$ and considering $\alpha+\beta=1$, we get $\Theta=0.5$. For the Jupiter-mass planet, $\Theta=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 555 Built-in in this estimate is our assumption that the Stokes number is constant. In fact, if a planet carves a deep gas gap, the Stokes number should go to infinity (as there is not gas, there is no drag). As such, this analysis applies only to partially depleted gaps (low and intermediate-mass planets).. In the simulation with photoelectric instability (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 in 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 concentrates 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 (lower plots) resembles the situation with $h=0.1$. The gas density is much lowered 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 Jupiter-mass planet. No clear distinction is present for the lower mass planet.

4. Caveats

The presented models are admittedly simplified. In this section, we state what we consider the main limitations of our calculations to be.

The most stringent limitation of the models is the 2D approximation. The twodimensionality imposes that the vertical direction is degenerate, so we do not resolve the vertical scale height. As a result the backreaction of the drag force is largely underestimated. In 3D, the midplane grain layer would be far denser due to sedimentation.

Another caveat is that we kept the Stokes number constant. In reality, the Stokes number in the Epstein regime should depend on density and temperature, according to Eq. (11). The effect of a self-consistent Stokes number in the PEI has not been explored, but it is not expected to be significant – as dust concentrates they increase the gas temperature, which also decreases the density through gas expansion. The former decreases $\mathrm{St}$, the latter increases $\mathrm{St}$, so both combined mitigate the effect on $\mathrm{St}$.

5. Conclusions

In this paper we have considered for the first time disk-planet interaction in the presence of photoelectric instability. We also study for the first time planet–disk interaction for disks of dust-to-gas ratio around unity, where gas and dust have similar inertia and none dominate over the other.

We find that the gaps carved in disks with dust and gas deviate from the gas-only case when the local dust-to-gas ratio approaches unity. The gap walls, because they are dust traps, are the places where the deviation occurs. At global dust-to-gas ratio $\varepsilon=10^{-2}$, the effect starts to become noticeable as a slightly reduction in the height of the gap wall. For $\varepsilon=10^{-1}$ and $\varepsilon=1$, the gap wall is significantly depressed as the local dust enhancement goes well above unity. This result does not depend on photoelectric instability and thus applies to optically thick (primordial and gas-rich transition) disks as well.

Our simulations show that when photoelectric instability is included, the features produced by low-mass planets are not conspicuous enough to be distinguishable from the features produced by the instability alone. The signature of gap carving and a conspicuous spiral are the tell-tale signatures of planets, as is usual for conventional transition disks. The most clear sign of a planet is a radial clearing in the PEI-induced structures, wider than the typical separation between them. That is, we can distinguish a planet amid of photoelectric instability structures if the planet carves a dust gap larger than the wavelength of the photoelectric instability. Quantitatively, we can state that the criterion to distinguish a planet in a disk with photoelectric instability is $\Delta r^{(d)}_{\rm gap}>\lambda_{\rm PEI}$, where $\Delta r^{(d)}_{\rm gap}$ is the radial width of the dust gap and $\lambda_{\rm PEI}$ the wavelength of the instability. The quantity $\Delta r^{(d)}_{\rm gap}$ is the simpler one to quantify. The boundary of the dust gap is located at the edge of the gas gap, that functions as a dust trap. As for $\lambda_{\rm PEI}$, linear instability analysis predicts growth even for $k\rightarrow\infty$ (Lyra & Kuchner, 2013). Regularizing with viscosity produces the most unstable wavelength of the order of the scale of viscous dissipation, which is very small in primordial and transitional disks, but of

order AU inlate transition and debris disks. Yet, a 3D run in (Lyra & Kuchner, 2013) shows ring spacing of $\approx 0.1H$ for a simulation with grid resolution $dx=H/256\approx 0.004$, i.e., $\lambda_{\rm PEI}/dx\approx 25$, and thus far from the viscous range. A specific prediction for $\lambda_{\rm PEI}$ is still lacking in the literature, a situation that is compounded by the fact that the PEI also exists in the nonlinear regime.

Our simulations also show that the PEI leads to well-defined sharp structures only for high dust-to-gas ratios, because of the pivotal role of the drag force backreaction in confining the dust. Simulations without backreaction, or with too little dust, $\varepsilon\approx 10^{-1}$, easily break into spiral structures and turbulence. We conclude that well-organized PEI structures should be expected only in disks with dust-to-gas ratio higher than primordial.

Finally, incomplete arcs are not produced in the simulations with planets alone: they require the photoelectric instability. However, we did not consider situations with self-gravity, which is thought to be a component in a possible mechanism to maintain incomplete arcs in the ring system of Neptune (Namouni & Porco, 2002; Tsui, 2018), albeit unlikely in the case of debris disks, of very low mass, to produce a noticeable effect.

The models of Richert et al. (2018) show that radiation pressure can have a significant effect on the global dust distributions of optically thin disks. The models we present in Section 3 are most straightforwardly applicable to optically thin disks around low-mass stars (spectral types K and M) and white dwarfs, whose low luminosities imply a minimal amount of radiation pressure even on small, well-coupled ($\tau_{\mathrm{f}}\varOmega=1$) grains. On the other hand, for more luminous stars, the greater radiation pressure will likely cause the disk morphologies to deviate substantially from the models presented in this paper. Nonetheless, our models and those of Richert et al. (2018) demonstrate that radiation pressure, massive planets, and the photoelectric instability can all play an important role in shaping optically thin disks, suggesting that future models should more extensively explore the combination of these and other effects. In particular, the MRI is another potentially important process in such systems (Kral & Latter, 2016) whose interplay with these other processes remains poorly explored.

W.L. acknowledges support of Space Telescope Science Institute through grant HST-AR-14572, and the NASA Exoplanet Research Program through grant 16-XRP16_2-0065. M.K. acknowledges support provided by NASA through a grant from the Space Telescope Science Institute (HST Cycle 21 AR13257.01), which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. The simulations presented in this paper utilized the Stampede2 cluster of the Texas Advanced Computing Center (TACC) at The University of Texas at Austin, through XSEDE grant TG-AS140014. We are indebted to the referee, Dr Sébastien Charnoz, for thoughtful comments that helped improve the manuscript.