The interplay between radiation pressure and the photoelectric instability in optically thin disks of gas and dust
Previous theoretical works have shown that in optically thin disks, dust grains are photoelectrically stripped of electrons by starlight, heating nearby gas and possibly creating a dust clumping instability—the photoelectric instability (PeI)—that significantly alters global disk structure. In the current work, we use the Pencil Code to perform the first numerical models of the PeI that include stellar radiation pressure on dust grains in order to explore the parameter regime in which the instability operates. In models with gas surface densities greater than , we see a variety of dust structures, including sharp concentric rings and, for models with especially high gas surface densities (), non-axisymmetric arcs and clumps that represent dust surface density enhancements of factors of 5–20 depending on the run parameters. The gas distributions show various structures as well, including clumps and arcs formed from spiral arms. In models with lower gas surface densities, vortices and smooth spiral arms form in the gas distribution, but the dust is too weakly coupled to the gas to be significantly perturbed. In one high gas surface density model, we include a large, low-order gas viscosity, and, in agreement with previous radiation pressure-free models, find that it observably smooths the structures that form in the gas and dust, suggesting that resolved images of a given disk may be useful for deriving constraints on the effective viscosity of its gas. Broadly, our models show that radiation pressure does not preclude the formation of complex structure from the PeI, but the qualitative manifestation of the PeI depends strongly on the parameters of the system. The PeI may provide an explanation for unusual disk morphologies such as the moving blobs of the AU Mic disk, the asymmetric dust distribution of the 49 Ceti disk, and the rings and arcs found in the disk around HD 141569A.
Circumstellar disks play a key role in testing theories of planet formation and evolution, revealing the physical and chemical environment of planet-forming systems, including providing constraints on the properties of nascent planets. Resolved images of protoplanetary disks, transitional disks, and debris disks show a variety of complex morphologies, including cavities, gaps and rings (Debes et al., 2013; Wahhaj et al., 2014; Follette et al., 2015; Currie et al., 2015; van Boekel et al., 2017), as well crescent-shaped structures, arcs, and spiral arms (van der Marel et al., 2013; Grady et al., 2013; Biller et al., 2015; Perrot et al., 2016; Follette et al., 2017). These disk structures are frequently attributed to gravitational perturbation by unseen embedded planets (e.g., Kuchner & Holman, 2003; Nesvold & Kuchner, 2015; Richert et al., 2015; Dong & Dawson, 2016; Dipierro & Laibe, 2017; Dong & Fung, 2017; Dong et al., 2017).
The possibility of comparable masses of gas and dust in any given optically thin disk raises the possibility of hydrodynamical interactions that will give rise to features like gaps, rings, and clumps that are frequently attributed to gravitational perturbation by an unseen embedded planet (Klahr & Lin, 2005; Besla & Wu, 2007; Lyra & Kuchner, 2013). In optically thin disks, stellar far ultraviolet photons whose energies exceed the work function of the dust grains (a few eV; Besla & Wu, 2007) photoelectrically eject electrons which in turn heat nearby gas. Klahr & Lin (2005) and Besla & Wu (2007) propose that this leads to a clumping instability—the photoelectric instability (PeI)—wherein heating of the gas by dust grains creates a local pressure maximum, which then traps more dust, which further heats the gas, and so on. Klahr & Lin (2005) model the system in 1D, finding the instability, and extrapolate from the 1D results to suggest that in 2D the photoelectric instability will generate ring structures, similar to those observed in disks like the one around HR 4796A. Lyra & Kuchner (2013) model the system hydrodynamically with 2D global and 3D local simulations, and find that rings are not formed unless the backreaction of the drag force is considered. When that component is ignored, power concentrates in high azimuthal wavenumbers and only clumps are formed. When the action of the dust on the gas is considered, rings and incomplete arcs are seen to form in the dust distribution.
Other findings of Lyra & Kuchner (2013) are that 1) linear instability exists only for dust-to-gas ratio , with maximum growth rate at ; 2) non-linear instability is observed for ; 3) linear instability only exists if photoelectric heating is the dominant heating source; 4) the photoelectric instability supersedes the streaming instability when the conditions for both are present; and 5) the particular mode for which gas and dust velocity are equal, thus canceling the drag force and backreaction, executes free oscillations, which are seen as a small but finite eccentricity (0.03).
The photoelectric instability may provide an explanation for a number of observed systems with unusual morphologies. Scattered light images of the AU Mic disk, an edge-on system, reveal radially moving blobs not seen at longer wavelengths. The disk around HIP 73145 contains concentric rings in scattered light images (which reveal small grains), while larger grains, observed at ALMA wavelengths, are distributed more compactly around the central star (Feldt et al., 2017). For both of these systems, the differing behavior of small and large grains is not readily explained by a planetary perturber. The edge-on disk around 49 Ceti is known to be gas rich, though the total mass remains poorly constrained (Hughes et al., 2017); Hughes et al. (2017) identify asymmetric structure in the disk consistent with a warp or spiral arm, finding no such features in two resolved gas-poor disks. The HD 141569A transition disk contains rings and arclets of small grains (Perrot et al., 2016). The findings of Klahr & Lin (2005), Besla & Wu (2007), and Lyra & Kuchner (2013) also raise the question of whether the presence of gas plays a role in the formation of the sharp dust rings seen in scattered light images of disks such as those around HR 4796A (Milli et al., 2017) and Fomalhaut (Kalas et al., 2008).
As novel as they are, the hydrodynamical models of Lyra & Kuchner (2013) do not include radiation pressure from the central star on dust grains, which even around low-mass stars will put sufficiently small grains on highly eccentric orbits, in some cases blowing them out of the system. The ability of the photoelectric instability to explain the morphologies of optically thin disks depends vitally on whether it can operate in the presence of stellar radiation pressure on dust grains.
In this work, we conduct hydrodynamical simulations of optically thin disks that include both dust–gas photoelectric heating and radiation pressure on dust grains that span a range of sizes. In Section 2, we provide an analytical discussion of the role grain size with respect to the emergence of hydrodynamical instabilities. Equations solved and initial conditions are discussed in Section 3. Results are discussed in Section 4, while further conclusions and implications for future work are discussed in Section 5.
2. The role of grain size
A priori, it may be expected that radiation pressure from the central star on dust grains in an optically thin disk will inhibit the formation of the clumps, arcs, and rings seen in the models of Lyra & Kuchner (2013). For spherical grains, the radiation pressure strength , defined as the ratio of the radiation force to the gravitational force, depends on host star mass , host star luminosity , grain radius , and dust material density (Burns et al., 1979; Krivov, 2010), such that
When a dust grain is created on a Keplerian orbit, radiation pressure places it on an eccentric orbit where eccentricity (Burns et al., 1979; Strubbe & Chiang, 2006). A grain with receives a radiation force equal to half the gravitational force, causing it to become unbound (). The models produced in this paper will help to determine whether non-zero orbital eccentricities of dust grains affect the onset of clumping instabilities.
Clumping due to the photoelectric instability depends on aerodynamic drag. Lyra & Kuchner (2013) find that the instability is robust for this variable, meaning that grains of longer stopping time simply take longer to respond to the pressure maximum and concentrate. Yet, one can imagine that if other dynamical processes are modifying the state of the gas at timescales shorter than the stopping time of the grains, clumping by the photoelectric instability may be disrupted. In other words, although linear growth is present, the saturated state may be quite different for small and big grains. Given the stopping time , a nondimensional stopping time can be constructed, also known as Stokes number, , where is the Keplerian frequency. For a thin disk,
where is the surface density of the gas. Grain size therefore seems likely to play a dual role in the development of dust clumping instabilities in optically thin disks, which require grains small enough to be susceptible to aerodynamic drag, but large enough to remain bound to the star, preferably on a low-eccentricity orbit.
e show as a function of for a solar-type star for several values of . The gas surface density values shown are a few orders of magnitude below the densities where the optical thickness of the gas will impede both photoelectric stripping of dust and radiation pressure.
A priori, it would seem that a substantial level of gas is required in order to rapidly create dust clumping instabilities in the presence of radiation pressure, while larger grains may be able to generate such instabilities over longer timescales (Lyra & Kuchner, 2013). Yet, clumping by PeI in the presence of radiation pressure may require grains of Stokes number near unity. A shorter PeI onset timescale would make the PeI to operate in disks spanning a large range of dust production rates (which remain poorly constrained in observed systems).
In this work we perform two-dimensional global simulations of gas-bearing optically thin circumstellar disks using the Pencil Code, a high-order finite difference hydrodynamics code (Brandenburg & Dobler, 2002). Both the gas and dust are calculated in Cartesian coordinates111In the course of this study, we identified a bug in the Pencil Code that affected the calculation of azimuthal particle accelerations in polar coordinates. The bug has since been fixed, and does not seem to have adversely affected previous works using that part of the code.. The code evolves the gas according to the continuity equation,
and the equation of motion,
where and are the gas and dust surface densities, and and are the velocity and pressure of the gas. represents the Newtonian gravitational potential of the star. The dust drag term is defined below.
We model the dust using 400,000 Lagrangian superparticles of equal mass, except for two runs where the total dust mass is increased in part by doubling the number of superparticles to 800,000 (see § 3.2). Each superparticle contains subparticles of one physical radius between 0.1 and 10 , and all subparticles have a constant material density , the approximate material density of silicate grains. This yields a total dust mass of .
The overall grain size distribution follows the standard Dohnanyi (1969) power law. The use of superparticles of the same mass automatically contributes a dependence with (superparticles with smaller grains will contain more particles); grain sizes associated with superparticles follow a dependence in order to yield a size distribution overall. This scheme prevents numerical issues associated with consolidating large grains into a small numbers of very massive superparticles (which, especially when modeling photoelectric heating, can crash the code).
In all simulations, superparticles are inserted in a birth ring positioned at 100 AU from the central star; the birth ring is axisymmetric and has a radial Gaussian profile () in order to avoid an artificial sharp edge in the dust distribution. Superparticles are at first inserted gradually over the course of several orbits in order to avoid discontinuities and ensure that superparticles span a range of orbital phases, but then are inserted throughout the simulation so as to yield a constant number of superparticles, and therefore constant total mass of the disk. That is to say, when a superparticle crosses the inner (50 AU) or outer (800 AU) boundary of the disk, another superparticle is inserted in the birth ring with a new grain size chosen at random according to the distribution.
The dynamical equation for each dust superparticle with velocity depends on the effective gravitational potential and gas drag term ,
The expression for the drag acceleration is given by
where is the Keplerian orbital frequency at a given orbital radius, and gas–dust velocity differential . The low gas densities of debris disks imply very large mean free paths compared with dust grain radii (i.e., ), hence the drag coefficient is calculated for the Epstein regime as a function of sound speed such that
The effective gravitational potential term incorporates radiation pressure on the dust, such that
where for each superparticle, the ratio of the radiation pressure force to gravitational force . We calculate reference radiation pressure strength for a solar-type star according to the prescription of Burns et al. (1979), such that for a 1 grain with density 2 g cm. Our grain size range corresponds with radiation pressure strengths , which in turn corresponds with eccentricities ranging from near-circular to completely unbound orbits, where an unbound orbit () has , corresponding with grain size .
Given that the photoelectric heating time is small compared with the dynamical time (Besla & Wu, 2007; Lyra & Kuchner, 2013), we adopt the modified equation of state of Lyra & Kuchner (2013) that implements instantaneous heating of gas by dust. The gas pressure is assumed to be proportional to the dust density , such that
is a dimensionless parameter that sets the pressure contribution of photoelectric heating compared to the background temperature of the gas (Lyra & Kuchner, 2013), which is itself specified by the reference sound speed in code units, corresponding with a scale height of 0.05 which is assumed for both the gas and dust. The value of is calculated by interpolating particle positions onto the grid using a triangular-shaped cloud particle mesh scheme (Eastwood & Hockney, 1974). We assume that the gas is locally isothermal, which is appropriate for the short cooling times expected in debris disks (Lyra & Kuchner, 2013).
The high-order scheme used by the Pencil Code leads to little numerical dissipation, therefore we apply sixth-order hyperdissipation terms to the r.h.s. of Equations 4 and 5 to stabilize the density and velocity fields, respectively, at the grid scale (Lyra et al., 2008; McNally et al., 2012; Lyra et al., 2017).
3.2. Model parameters
We conduct seven simulations in total. For five models—runs A–E—we vary the initial gas surface density according to the values shown in
his range of surfaces densities corresponds with surface densities times lower than that of the Minimum Mass Solar Nebula at 100 AU (Weidenschilling, 1977). For the Pic disk, Brandeker et al. (2004) assume solar abundances to derive a gas surface density of , based on a scale height of 10 AU and mean molecular mass 2.5 amu. Our models have gas surface densities approximately 0.03 to 300 times that value. In run F, we duplicate our model with the highest value of (run E) but increase the total dust mass by a factor of 10; we achieve this by doubling the number of superparticles and increasing the mass per superparticle by a factor of five. In run G, we duplicate run F but add a large Laplacian viscosity to the r.h.s. of Eq. (5), corresponding with a Shakura–Sunyaev viscosity of 0.1 (Shakura & Sunyaev, 1973) at the reference orbital radius (100 AU). Lyra & Kuchner (2013) find that viscosity damps the PeI at high wavenumbers. This viscosity is presumably expected from the magnetorotational instability (MRI; Balbus & Hawley, 1991). Although the activity of the MRI in optically thin disks remains limitedly understood (Kral & Latter, 2016), we include this additional run to explore the effect of an eddy viscosity when photoelectric heating and radiation pressure are both operative.
The physical parameters that differentiate our models are summarized in Table 1, including the Stokes number of the smallest particle in each run, min() (calculated using Eq. 2 for ). Table 1 also provides the spatially averaged dust-to-gas ratio for after 400 orbits, denoted as , as well as the characteristic vertical geometric optical depth of any dust overdensities after 400 orbits, denoted as .
|Run||Total dust mass||min()|
|( )||()||(400 orbits)|
In each model, the gas is assumed to be initially uniformly distributed throughout the disk between 50 AU and 800 AU from the central star. The gas temperature is also uniform throughout the disk. This yields no global pressure gradient, allowing us to isolate the effects of photoelectric heating; specifically, it allows us to attribute any radial dust drift to the radiation pressure and PeI. Lyra & Kuchner (2013) show for the radiation pressure-free case that the photoelectric instability generates dust rings in the presence of a global pressure gradient and the streaming instability.
We run each model for 400 orbits, a sufficient amount of time to determine whether small-to-medium grains can trigger the formation of clumps or other features through the PeI. The largest grains in the lowest-gas runs are so poorly coupled that resolving the PeI growth timescale associated with them would require prohibitively long run times; in these low-gas runs, any participation of large grains in PeI-induced effects would presumably require some complex interplay between small and large grains. For instance, small grains, even unbound ones, could trigger gas overdensities that yield substantially shorter values of , leading to better coupling of larger grains.
To guide the interpretation of the simulations, we compute the analytical growth rates of the PeI as a function of Stokes number for the range considered. This is done by solving Eq. 26–29 of Lyra & Kuchner (2013). Without viscosity, the growth rates would grow unboundedly with wavenumber, eventually getting unphysically high at the grid scale. In reality the growth rates drop abruptly when the viscous range is approached. Because of this, we regularize the system with Laplacian viscosity. The result for , as in run G, is plotted in
he figure shows the maximum growth rate as a function of Stokes number and dust-to-gas ratio. The labels A–G indicate the minimum Stokes numbers and dust to gas ratios corresponding to simulations A–G (if runs A–F had Laplacian viscosity rather than hyperviscosity). Above dust-to-gas ratio unity no linear instability exists. The symmetry with respect to , seen in the nonlinear simulations of Lyra & Kuchner (2013) (Supplement Fig. 2 of that paper) is reproduced. The runs of this paper are labeled in the graph. The value chosen as representative of the run is at . In each simulation, the range of should reach an order of magnitude in each direction.
For runs A and B (Figs. 6 and 7), we see that a narrow, axisymmetric dust ring forms just outside the birth ring, with a central surface density corresponding with . In the gas, two vortices form just outside the birth ring. In run B, a gas gap appears, along with two additional vortices appear on the opposite side of the gap from the first two, matching them in azimuth. Each vortex orbits at sub-Keplerian speed, with an orbital frequency approximately 90% of that expected from Keplerian rotation.
In order to further investigate this behavior, we plot the gas surface density over time for run B in Figure 13. In the first few dozen orbits, a single vortex () emerges. In the next several dozen orbits, a gas gap forms, as well as two vortices () just outside it, positioned apart from each other in azimuth; the contemporaneity of the formation of the gap and the two inner vortices suggests that the Rossby wave instability (Lovelace et al., 1999) may be responsible. For the next several hundred orbits, two additional “matching” vortices appear just within the orbit of the inner gap edge, keeping pace with the outer vortices. After 520 orbits, the outer vortices have begun to migrate and the inner vortices are much less pronounced, and after 560 orbits only one inner–outer vortex pair remains.
In runs A–C (Figs. 6–8), we see that the increasing value of (thus reducing ) results in a narrow gas gap whose depth increases with . In run D (Fig. 9), we see a shallower, wider gap threaded by several spiral arm structures (not seen in the models of Lyra & Kuchner 2013).
In runs C and D (Figs. 8 and 9), compared with runs A and B, the gas distributions are less smooth and contain non-axisymmetric clumps and arcs. The dust distributions show concentric rings with factor of 10–20 dust enhancements, reminiscent of the models of Lyra & Kuchner (2013), though there a fewer rings in runs C and D, which have two and three main ring structures, respectively. The central optical depths of these rings are of order , with the exception of the outermost ring in run D (Fig. 9), which is about an order of magnitude fainter. In run C, the anticorrelation of gas and dust just outside the birth ring is reminiscent of the dust–gas anticorrelation seen in the models of Lyra & Kuchner (2013).
In run D, the dust rings are accompanied by more high-frequency structure, including some localized ripple patterns (a few AU in scale) strongly reminiscent of the tightly packed arcs and rings seen in the Lyra & Kuchner (2013) models. These rings are long-lived, having begun to form after only a few tens of orbits. Though the gas responds to the dust on dynamical timescales due to photoelectric heating, in runs A–D the Stokes numbers are high enough that the bound grains () respond to the gas only on much longer timescales; grains on low-eccentricity orbits see the azimuthally-averaged gas distribution over the course of many orbits, not the gas’s spiral structure. This asymmetrical coupling promotes the formation of dust rings. But note that some of the fine structure in the rings, like the bifurcation at the one o’clock position, seems to correspond with clumps in the gas.
No indications of the photoelectric instability emerge after 400 orbits for run E (Fig. 10). As shown in
he growth rate for run E is centered at the maximum value of . After 400 orbits, this low growth rate amounts to a mere 2% amplification. Runs D and C, though centered at the same low level of growth rate, and B, centered at even lower, reach to the left of the diagram and into regions of higher growth rate, as much as , (million-fold amplification in 220 orbits). Run A, though also centered at a low value of growth rate, and not reaching too deep into regions of high growth rate, may be non-linearly amplified due to the high dust-to-gas ratio.
For run F (Fig. 11), we find that increasing the total dust mass compared with run E yields the return of the photoelectric instability. The photoelectric instability radically transforms both the gas and dust distributions; the dust clumps and arcs seen in the right panel of Figure 11 correspond with , representing factor of 5–10 enhancement over the local dust surface density.
In Figure 15, we show the dust surface density every 50 orbits for 400 orbits. We find that the photoelectric instability sets in quickly, but has no obvious secular effect on the global structure of the disk. The behavior is dominated by transient clumps and arcs that appear and disappear on timescales of orbits to dozens of orbits. Qualitatively, the behavior of run F is quite different from that of run D, where the dust remains in mostly smooth rings. As seen in Figure 2, for run D, bound dust grains are weakly coupled with the gas, generating rings as discussed above. In run F, a large range of grain sizes are both bound and well-coupled to the gas. Setting and in Equation 3, we find that the threshold for having grains that are both bound and well-coupled occurs at a gas surface density of
The gas and dust distributions for run F (Fig. 11) are strikingly similar to the results of the Lyra & Kuchner (2013) model that excludes drag backreaction on the gas. In the absence of radiation pressure, where both the gas and dust move on roughly circular orbits, the expansion of gas due to high pressures induced by photoelectric heating will undergo Coriolis rotation. In models with backreaction, this rotation is opposed by the backreaction from the dust, and axisymmetry of the system is maintained. It is possible that when a dust grain is placed on a highly eccentric orbit by radiation pressure, it heats gas in a given region, creating a pressure maximum, but does not linger in a nearby circular orbit where it can stabilize the gas through drag backreaction.
Lyra & Kuchner (2013) find that in the presence of drag with backreaction and photoelectric heating, the gas and dust mutually displace each other, leading to alternating rings of gas and dust throughout the disk. It is apparent in Figures 6–8 that when a small amount of gas is present, dust displaces it, creating a gas gap. In order to test whether gas and dust anticorrelate when larger quantities of gas are present, in Figure 16 we plot the product of and (specifically, their mean-subtracted and standard deviation-normalized values) for run F after 400 orbits. We find that the dust and gas correlate (red regions) and anticorrelate (blue region), in roughly equal measure.
We find that the inclusion of a very strong low-order gas viscosity term in run G (Fig. 12) yields results are fairly similar to the no-viscosity case (run F; Fig. 11). It does, however, lead to the smoothing of some of the small-scale structure in both the gas and dust distributions seen in run F; dust surface density enhancement are of order factor of 5, somewhat smaller than seen in run F. The small, transient dust clumps seen in run F are less numerous and less elongated in run G, but are nonetheless present, also appearing and disappearing over orbital timescales. This result is consistent with the predictions and models of Lyra & Kuchner (2013), where viscosity suppresses the PeI at high wavenumbers, smoothing small-scale structure but not impeding the ability of the PeI to substantially reshape the disk.
In order to explore the roles played by different grain sizes in creating and sustaining the structures seen in the right-hand panels of Figures 6–12, we examine the radial distributions of dust in several bins of grain size, choosing runs A and F as representative cases. The upper, middle, and lower panels of Figure 17 show radial dust distributions for run A after 400 orbits, run F after 20 orbits, and run F after 400 orbits, respectively.
The radial dust distributions shown in the upper panel of Figure 17 confirm that it is large particles (i.e., small , corresponding with nearly circular orbits) that constitute the dust ring seen in run A (Fig. 6). This is consistent with the fact that these sharp rings take dozens of orbits to form, reflecting the large values of for large grains. The small values of place these large grains on circular orbit, keeping them close to the birth ring.
In the middle and lower panels of Figure 17, we identify two important grain size-dependent effects. The first, seen in the lower panel, is that even large grains (7–10 ) migrate outward from their low-eccentricity orbits in and around the birth ring, having now had enough time () to be entrained in PeI-induced flows. The second is that smaller grains are more efficiently entrained by the denser gas, leading ordinarily unbound particles () to remain in the disk. The middle panel shows the grain distribution for run F after only 20 orbits, at which point the photoelectric instability is just forming. Small-to-medium grains (0.1–3 ), as they are being blown outward by radiation pressure, accumulate outside the birth ring and trigger the photoelectric instability. Meanwhile, larger and therefore more weakly coupled grains migrate outward by gas drag over the course of orbits to tens of orbits. The dominance of large grains seen in the lower panel of Figure 17 suggests that while small, unbound grains help to trigger the photoelectric instability and participate in the resulting transient structures, many of them will ultimately be blown out, and the dust structure seen in Figure 11 eventually becomes dominated by larger, bound grains. This suggests that even in disks around more massive (say, A-type) stars, the extreme radiation pressure on dust grains (and subsequent lack of bound grains) does not necessarily inhibit the formation of the photoelectric instability. Note that Figures 3–9 show the dust surface density, which tends to be dominated by larger grains. But we found that plotting the optical depth (not shown) results in images that are nearly indistinguishable, despite emphasizing smaller grains.
We have produced hydrodynamical models of optically thin disks with gas and dust, simultaneously incorporating photoelectric heating and stellar radiation pressure on dust grains. We find the following:
The PEI growth rate is small for gas surface densities , but fast enough at higher gas surface densities that we see the PEI create dust density enhancements of up to a factor of 20 in our runs of just 400 orbits.
For a modest level of gas (; runs C and D), the photoelectric instability gives rise to axisymmetric dust rings over the course of dozens of orbits, as well as azimuthal structure in the gas. The dust and gas show the anticorrelation bevahior predicted by Lyra & Kuchner (2013).
For a higher level of gas (; run E), and similar dust-to-gas ratio , (runs F and G), the PeI emerges over the course of dozens of orbits, leading to erratic spiral structure throughout the disk (resembling the backreaction-free model of Lyra & Kuchner 2013), with small-scale structure appearing and disappearing over the course of orbits. In these models, there is no overall tendency of the dust and gas to anticorrelate.
The value of required to generate the PeI on short timescales in a given system, however, will depend on many parameters, including spectral type, total dust mass, initial gas profile, and grain size range. In some of our models with higher level of gas (; run E), but low effective dust-to-gas ratio , we found the PeI growth rate too slow to be captured by our simulations.
Previous models of debris disks with gas have lacked the physics to capture hydrodynamical instabilities like the PeI. Thébault & Augereau (2005) and Krivov et al. (2009) present models of debris disks with gas spanning a range of gas surface densities similar to the range used in the current work (though other key parameters vary between these works, such as stellar spectral type and the initial radial profile of the gas). In both cases, dust grains generated in a birth ring experience aerodynamic drag with a static gas cloud, precluding the emergence of the PeI, which requires dust–gas heating and drag backreaction. These one-dimensional models yield dust density distributions that decrease smoothly and monotonically with orbital radius, making them readily distinguishable from the clumps, arcs, and narrow rings seen in the models presented in Section 4.
Other models of optically thin disks reveal more complex morphologies that are not as readily distinguished from the results presented in this paper. For instance, the models of cataclysmic massive body collisions produced by Kral et al. (2015) produce non-axisymmetric structures that could be difficult to distinguish from the dust distributions shown in Figures 11 and 12, especially for disks in an edge-on viewing configuration. Nonetheless, the models of Kral et al. (2015) show smoother, more organized structure compared with the higher-frequency, more erratic structure seen in runs F and G (Figs. 11 and 12), and also do not show concentric, axisymmetric rings as in runs C and D. Comparisons of dust distributions at multiple wavelengths for a given disk could also help to disentangle these two effects, given that differential behavior by grain size should be greater for aerodynamic effects than gravitational ones.
Augereau & Papaloizou (2004) model a circumstellar debris-only disk with an external stellar perturber in order to study the origins of the spiral morphology of the HD 141569 disk. In general, their models produce smooth spiral structure in the circumprimary disk, however for perturbers with eccentric orbits, the structure produced in the disk is less smooth and somewhat reminiscent of the dust distributions seen in runs F and G (Figs. 11 and 12). Here again, multiwavelength image comparisons may be necessary to distinguish these models. Given that none of the numerical models so far produced of the PeI have given rise to smooth spiral arms of dust (Lyra & Kuchner 2013 and the current work), it would seem that perturbation by a massive companion is currently a more plausible explanation for such structures (though this may change as the PeI is modeled throughout a larger parameter space).
For a number of other observed disk morphologies, however, the PeI may provide a more plausible explanation than the presence of a massive perturber. The resemblance of the dust distributions for runs C and D (Figs. 8 and 9) to the concentric rings seen in the scattered light profile of the disk around HIP 73145 (spectral type A; Feldt et al., 2017) is particularly striking. The large grains are relatively compactly distributed, while small grains experience considerable radiation pressure. Though we do not model A-type stars in the current work, our models suggest that the photoelectric instability provides a promising explanation for such features, which can be triggered by unbound grains.
Also, the non-axisymmetric clumps and spiral arms seen in the gas distributions in several of our models suggest that interactions between gas and small-to-medium grains could underlie the asymmetric structures seen in the 49 Ceti and AU Mic edge-on disks (Hughes et al., 2017; Boccaletti et al., 2015), though explanations for the AU Mic moving blobs involving only rocky body collisions (no gas) have also been proposed (Sezestre et al., 2017; Chiang & Fung, 2017). We also underscore that so far the PeI is the best candidate for producing arcs, a feature that is not predicted from planet-disk interaction. That the PeI leads to arcs was predicted by Lyra & Kuchner (2013), before the discovery of these features in the disk around HD 141569A (Perrot et al., 2016).
The broad resemblance of our models to several observed systems notwithstanding, models that include more physical detail—magnetic fields, multiple gas species, and so on—and also explore a wider range of parameters (especially stellar spectral type) will help to confirm that the photoelectric instability can indeed provide a plausible explanation for these diverse and intriguing disk morphologies.
Future hydrodynamical models of optically thin disks should explore a number of physical processes not explored in the current work:
The magnetorotational instability (MRI) may operate efficiently in debris disks (Kral & Latter, 2016). The inclusion of a low-order viscous term in run G produced an observably different dust distribution from run F; the role of MRI-induced turbulence in redistributing gas, and subsequently dust, should be explored in detail.
Future studies should explore the competing roles of gas–gas and dust–gas photoelectric heating to determine the precise realm of disk parameter space in which dust–gas photoelectric heating is dynamically important. Kral et al. (2016) point out that in a carbon-rich disk like that around Pic, heating by photoelectrons released from carbon may overwhelm the heating produced by those released from dust. This in turn may inhibit the photoelectric instability by inhibiting the formation of pressure maxima.
In the current work, we model only a single gas species. However, as noted by Xie et al. (2013), different gas species can experience different values of . The presence of a modest radiation force on a certain gas species with, say, would significantly affect the mutual velocity of that species with the dust, potentially altering the effects of the PeI.
Future investigations should explore the role of Poynting–Robertson drag in the context of the photoelectric instability. The findings of Lyra & Kuchner (2013) suggest that even large, poorly coupled grains will eventually give rise to the photoelectric instability, but on very long timescales, inward drift due to PR drag could inhibit or affect this process.
In the current work, we have only modeled disks around solar-type stars. Around A-type stars, where many debris disks, including those with unusual morphologies, are found, the radiation force on dust grains will be considerably greater. This means that the bound grains will be larger (though the results presented in Fig. 17 suggest that even unbound grains—small and well-coupled—can trigger the PeI on their way out of the disk). Large and therefore weakly-coupled grains can generate and participate in the photoelectric instability, however the resulting structures will take longer to emerge (a likely explanation for the lack of dust clumping in runs A and B). It is possible for instance that the spatial and temporal frequencies of any non-axisymmetric structures resulting from gas–dust interactions involving large grains will be larger due to the larger values of .
Once the interactions of dust, gas, and radiation in optically thin disks are better understood, the signposts of embedded planets can be more accurately modeled. The substantial effect of drag backreaction and photoelectric heating seen in the planet-free models of the current work suggests that planet–disk interactions in optically thin disks may manifest themselves very differently than they do in existing, simpler models. The formation of gaps, rings, and other disk morphologies associated with planets may be inhibited, enhanced, or otherwise affected by the presence of gas.
- Augereau & Papaloizou (2004) Augereau, J. C. & Papaloizou, J. C. B. 2004, A&A, 414, 1153
- Balbus & Hawley (1991) Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214
- Besla & Wu (2007) Besla, G. & Wu, Y. 2007, ApJ, 655, 528
- Biller et al. (2015) Biller, B. A., Liu, M. C., Rice, K., Wahhaj, Z., Nielsen, E., Hayward, T., Kuchner, M. J., Close, L. M., Chun, M., Ftaclas, C., & Toomey, D. W. 2015, MNRAS, 450, 4446
- Boccaletti et al. (2015) Boccaletti, A., Thalmann, C., Lagrange, A.-M., Janson, M., Augereau, J.-C., Schneider, G., Milli, J., Grady, C., Debes, J., Langlois, M., Mouillet, D., Henning, T., Dominik, C., Maire, A.-L., Beuzit, J.-L., Carson, J., Dohlen, K., Engler, N., Feldt, M., Fusco, T., Ginski, C., Girard, J. H., Hines, D., Kasper, M., Mawet, D., Ménard, F., Meyer, M. R., Moutou, C., Olofsson, J., Rodigas, T., Sauvage, J.-F., Schlieder, J., Schmid, H. M., Turatto, M., Udry, S., Vakili, F., Vigan, A., Wahhaj, Z., & Wisniewski, J. 2015, Nature, 526, 230
- Brandeker et al. (2004) Brandeker, A., Liseau, R., Olofsson, G., & Fridlund, M. 2004, A&A, 413, 681
- Brandenburg & Dobler (2002) Brandenburg, A. & Dobler, W. 2002, Computer Physics Communications, 147, 471
- Burns et al. (1979) Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icarus, 40, 1
- Chiang & Fung (2017) Chiang, E. & Fung, J. 2017, ArXiv e-prints
- Currie et al. (2015) Currie, T., Lisse, C. M., Kuchner, M., Madhusudhan, N., Kenyon, S. J., Thalmann, C., Carson, J., & Debes, J. 2015, ApJ, 807, L7
- Debes et al. (2013) Debes, J. H., Jang-Condell, H., Weinberger, A. J., Roberge, A., & Schneider, G. 2013, ApJ, 771, 45
- Dipierro & Laibe (2017) Dipierro, G. & Laibe, G. 2017, MNRAS, 469, 1932
- Dohnanyi (1969) Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531
- Dong & Dawson (2016) Dong, R. & Dawson, R. 2016, ApJ, 825, 77
- Dong & Fung (2017) Dong, R. & Fung, J. 2017, ApJ, 835, 146
- Dong et al. (2017) Dong, R., Li, S., Chiang, E., & Li, H. 2017, ApJ, 843, 127
- Eastwood & Hockney (1974) Eastwood, J. W. & Hockney, R. W. 1974, Journal of Computational Physics, 16, 342
- Feldt et al. (2017) Feldt, M., Olofsson, J., Boccaletti, A., Maire, A. L., Milli, J., Vigan, A., Langlois, M., Henning, T., Moor, A., Bonnefoy, M., Wahhaj, Z., Desidera, S., Gratton, R., Kóspál, Á., Abraham, P., Menard, F., Chauvin, G., Lagrange, A. M., Mesa, D., Salter, G., Buenzli, E., Lannier, J., Perrot, C., Peretti, S., & Sissa, E. 2017, A&A, 601, A7
- Follette et al. (2015) Follette, K. B., Grady, C. A., Swearingen, J. R., Sitko, M. L., Champney, E. H., van der Marel, N., Takami, M., Kuchner, M. J., Close, L. M., Muto, T., Mayama, S., McElwain, M. W., Fukagawa, M., Maaskant, K., Min, M., Russell, R. W., Kudo, T., Kusakabe, N., Hashimoto, J., Abe, L., Akiyama, E., Brandner, W., Brandt, T. D., Carson, J., Currie, T., Egner, S. E., Feldt, M., Goto, M., Guyon, O., Hayano, Y., Hayashi, M., Hayashi, S., Henning, T., Hodapp, K., Ishii, M., Iye, M., Janson, M., Kandori, R., Knapp, G. R., Kuzuhara, M., Kwon, J., Matsuo, T., Miyama, S., Morino, J.-I., Moro-Martin, A., Nishimura, T., Pyo, T.-S., Serabyn, E., Suenaga, T., Suto, H., Suzuki, R., Takahashi, Y., Takato, N., Terada, H., Thalmann, C., Tomono, D., Turner, E. L., Watanabe, M., Wisniewski, J. P., Yamada, T., Takami, H., Usuda, T., & Tamura, M. 2015, ApJ, 798, 132
- Follette et al. (2017) Follette, K. B., Rameau, J., Dong, R., Pueyo, L., Close, L. M., Duchêne, G., Fung, J., Leonard, C., Macintosh, B., Males, J. R., Marois, C., Millar-Blanchaer, M. A., Morzinski, K. M., Mullen, W., Perrin, M., Spiro, E., Wang, J., Ammons, S. M., Bailey, V. P., Barman, T., Bulger, J., Chilcote, J., Cotten, T., De Rosa, R. J., Doyon, R., Fitzgerald, M. P., Goodsell, S. J., Graham, J. R., Greenbaum, A. Z., Hibon, P., Hung, L.-W., Ingraham, P., Kalas, P., Konopacky, Q., Larkin, J. E., Maire, J., Marchis, F., Metchev, S., Nielsen, E. L., Oppenheimer, R., Palmer, D., Patience, J., Poyneer, L., Rajan, A., Rantakyrö, F. T., Savransky, D., Schneider, A. C., Sivaramakrishnan, A., Song, I., Soummer, R., Thomas, S., Vega, D., Wallace, J. K., Ward-Duong, K., Wiktorowicz, S., & Wolff, S. 2017, AJ, 153, 264
- Grady et al. (2013) Grady, C. A., Muto, T., Hashimoto, J., Fukagawa, M., Currie, T., Biller, B., Thalmann, C., Sitko, M. L., Russell, R., Wisniewski, J., Dong, R., Kwon, J., Sai, S., Hornbeck, J., Schneider, G., Hines, D., Moro Martín, A., Feldt, M., Henning, T., Pott, J.-U., Bonnefoy, M., Bouwman, J., Lacour, S., Mueller, A., Juhász, A., Crida, A., Chauvin, G., Andrews, S., Wilner, D., Kraus, A., Dahm, S., Robitaille, T., Jang-Condell, H., Abe, L., Akiyama, E., Brandner, W., Brandt, T., Carson, J., Egner, S., Follette, K. B., Goto, M., Guyon, O., Hayano, Y., Hayashi, M., Hayashi, S., Hodapp, K., Ishii, M., Iye, M., Janson, M., Kandori, R., Knapp, G., Kudo, T., Kusakabe, N., Kuzuhara, M., Mayama, S., McElwain, M., Matsuo, T., Miyama, S., Morino, J.-I., Nishimura, T., Pyo, T.-S., Serabyn, G., Suto, H., Suzuki, R., Takami, M., Takato, N., Terada, H., Tomono, D., Turner, E., Watanabe, M., Yamada, T., Takami, H., Usuda, T., & Tamura, M. 2013, ApJ, 762, 48
- Hughes et al. (2017) Hughes, A. M., Lieman-Sifry, J., Flaherty, K. M., Daley, C. M., Roberge, A., Kóspál, Á., Moór, A., Kamp, I., Wilner, D. J., Andrews, S. M., Kastner, J. H., & Ábrahám, P. 2017, ApJ, 839, 86
- Kalas et al. (2008) Kalas, P., Graham, J. R., Chiang, E., Fitzgerald, M. P., Clampin, M., Kite, E. S., Stapelfeldt, K., Marois, C., & Krist, J. 2008, Science, 322, 1345
- Klahr & Lin (2005) Klahr, H. & Lin, D. N. C. 2005, ApJ, 632, 1113
- Kral & Latter (2016) Kral, Q. & Latter, H. 2016, MNRAS, 461, 1614
- Kral et al. (2015) Kral, Q., Thébault, P., Augereau, J.-C., Boccaletti, A., & Charnoz, S. 2015, A&A, 573, A39
- Kral et al. (2016) Kral, Q., Wyatt, M., Carswell, R. F., Pringle, J. E., Matrà, L., & Juhász, A. 2016, MNRAS, 461, 845
- Krivov (2010) Krivov, A. V. 2010, Research in Astronomy and Astrophysics, 10, 383
- Krivov et al. (2009) Krivov, A. V., Herrmann, F., Brandeker, A., & Thébault, P. 2009, A&A, 507, 1503
- Kuchner & Holman (2003) Kuchner, M. J. & Holman, M. J. 2003, ApJ, 588, 1110
- Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805
- Lyra et al. (2008) Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2008, A&A, 479, 883
- Lyra & Kuchner (2013) Lyra, W. & Kuchner, M. 2013, Nature, 499, 184
- Lyra et al. (2017) Lyra, W., McNally, C. P., Heinemann, T., & Masset, F. 2017, ArXiv e-prints
- McNally et al. (2012) McNally, C. P., Lyra, W., & Passy, J.-C. 2012, ApJS, 201, 18
- Milli et al. (2017) Milli, J., Vigan, A., Mouillet, D., Lagrange, A.-M., Augereau, J.-C., Pinte, C., Mawet, D., Schmid, H. M., Boccaletti, A., Matrà, L., Kral, Q., Ertel, S., Chauvin, G., Bazzon, A., Ménard, F., Beuzit, J.-L., Thalmann, C., Dominik, C., Feldt, M., Henning, T., Min, M., Girard, J. H., Galicher, R., Bonnefoy, M., Fusco, T., de Boer, J., Janson, M., Maire, A.-L., Mesa, D., Schlieder, J. E., & SPHERE Consortium. 2017, A&A, 599, A108
- Nesvold & Kuchner (2015) Nesvold, E. R. & Kuchner, M. J. 2015, ApJ, 798, 83
- Perrot et al. (2016) Perrot, C., Boccaletti, A., Pantin, E., Augereau, J.-C., Lagrange, A.-M., Galicher, R., Maire, A.-L., Mazoyer, J., Milli, J., Rousset, G., Gratton, R., Bonnefoy, M., Brandner, W., Buenzli, E., Langlois, M., Lannier, J., Mesa, D., Peretti, S., Salter, G., Sissa, E., Chauvin, G., Desidera, S., Feldt, M., Vigan, A., Di Folco, E., Dutrey, A., Péricaud, J., Baudoz, P., Benisty, M., De Boer, J., Garufi, A., Girard, J. H., Menard, F., Olofsson, J., Quanz, S. P., Mouillet, D., Christiaens, V., Casassus, S., Beuzit, J.-L., Blanchard, P., Carle, M., Fusco, T., Giro, E., Hubin, N., Maurel, D., Moeller-Nilsson, O., Sevin, A., & Weber, L. 2016, A&A, 590, L7
- Richert et al. (2015) Richert, A. J. W., Lyra, W., Boley, A., Mac Low, M.-M., & Turner, N. 2015, ApJ, 804, 95
- Sezestre et al. (2017) Sezestre, É., Augereau, J.-C., Boccaletti, A., & Thébault, P. 2017, ArXiv e-prints
- Shakura & Sunyaev (1973) Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337
- Strubbe & Chiang (2006) Strubbe, L. E. & Chiang, E. I. 2006, ApJ, 648, 652
- Thébault & Augereau (2005) Thébault, P. & Augereau, J.-C. 2005, A&A, 437, 141
- van Boekel et al. (2017) van Boekel, R., Henning, T., Menu, J., de Boer, J., Langlois, M., Müller, A., Avenhaus, H., Boccaletti, A., Schmid, H. M., Thalmann, C., Benisty, M., Dominik, C., Ginski, C., Girard, J. H., Gisler, D., Lobo Gomes, A., Menard, F., Min, M., Pavlov, A., Pohl, A., Quanz, S. P., Rabou, P., Roelfsema, R., Sauvage, J.-F., Teague, R., Wildi, F., & Zurlo, A. 2017, ApJ, 837, 132
- van der Marel et al. (2013) van der Marel, N., van Dishoeck, E. F., Bruderer, S., Birnstiel, T., Pinilla, P., Dullemond, C. P., van Kempen, T. A., Schmalzl, M., Brown, J. M., Herczeg, G. J., Mathews, G. S., & Geers, V. 2013, Science, 340, 1199
- Wahhaj et al. (2014) Wahhaj, Z., Liu, M. C., Biller, B. A., Nielsen, E. L., Hayward, T. L., Kuchner, M., Close, L. M., Chun, M., Ftaclas, C., & Toomey, D. W. 2014, A&A, 567, A34
- Weidenschilling (1977) Weidenschilling, S. J. 1977, Ap&SS, 51, 153
- Xie et al. (2013) Xie, J.-W., Brandeker, A., & Wu, Y. 2013, ApJ, 762, 114