HIDING IN THE SHADOWS: SEARCHING FOR PLANETS IN PRE-TRANSITIONAL AND TRANSITIONAL DISKS

Transitional and pre-transitional disks are characterized by a flux deficit in the near-infrared wavelengths with respect to a classic protoplanetary disk, but otherwise have a similar spectral energy distribution (SED). The accepted explanation for the deficit is a lack of optically thick material in part or all of the inner accretion disk (Strom et al. 1989; Calvet et al. 2005). Pre-transitional disks have a gap in the accretion disk with hot optically thick material on the inside, while transitional disks have a complete hole and no detectable hot component to the SED (Espaillat et al. 2007).

The origin of the gap or hole in transitional objects is debated and many hypotheses have been suggested, such as particle growth (Dullemond & Dominik 2005), photoevaporation (Alexander & Armitage 2007), and magnetorotational instability (Chiang & Murray-Clay 2007). However, none of these processes can create the largest holes and gaps observed, such as those observed in GM Aur, DM Tau, and SAO 206462, which are tens of AU in radius (Calvet et al. 2005; Hughes et al. 2009; Mayama et al. 2012; Muto et al. 2012). This leads us to another very interesting hypothesis for the origin of the interruption in the protoplanetary accretion disk, namely, embedded planetary object(s). The largest holes and gaps require multiple embedded giant planets to explain their extreme size (Rice et al. 2006; Papaloizou et al. 2007; Dodson-Robinson & Salyk 2011) but many other less extreme transitional objects are consistent with one embedded, unobserved giant planet (Owen et al. 2011).

Transitional systems are all young (⩽10 Myr), thus, if there is an undetected giant planet in the system, there should also be a background planetesimal population. The orbits of the background population of planetesimals will be perturbed by the planet (Charnoz et al. 2001), increasing collision frequency in regions close to the embedded planet. Dust created will be visible as long as the planetesimals remain collisionally active. This second-generation dust will augment any primordial dust in both the outer gas-rich disk and the inner gas-poor hole by some fraction that is dependent on the survival time of second generation dust grains. In any disk where planet formation is in progress, second-generation dust will significantly enhance the dust grain population.

In this Letter we present a search for observational signatures of dust produced by planetesimals under the influence of a giant planet. We run two sets of N-body simulations of a giant planet embedded in a planetesimal population. In one simulation set the collision model (rubble) self-consistently treats debris production and in the other set (perfect merging) all planetesimal collisions are inelastic, with an analytical prescription for the size distribution of dust particles. We then use a radiative transfer code to create synthetic submillimeter images of the dust disk and identify features that change based on the planet orbit.

The numerical simulations presented here were calculated using the parallelized N-body gravity code PKDGRAV (Richardson et al. 2000; Stadel 2001; Leinhardt & Richardson 2005). PKDGRAV uses a hierarchical tree to calculate inter-particle gravity and a second order leapfrog integrator for time evolution. Each simulation took between three and six months to complete on 32 2.8 GHz processors.

2.1. Initial Conditions and Simulation Details

Most simulations begin with 10⁶ equal-mass planetesimals that interact with each other through gravitational forces and physical collisions. Ten simulations were performed in total, in simulations 0–8 × 10⁶ planetesimals are initially placed in a 3 AU wide ring from 0.8 to 3.8 AU around a 1 M_☉ potential. Simulation 9 has lower particle resolution (N = 10⁵), but a much wider initial disk extending from 0.8 to 10 AU. The embedded planet is placed at 2.8 AU, highlighted by a green circle in Figure 1. The total mass of planetesimals is 4.7 M_⊕ (Richardson et al. 2000), with a bulk density of 2 g cm⁻³ (Weidenschilling 1977). Planetesimals are initially distributed such that the surface density decreases with the semi-major axis, Σ(r)∝r^−1.5 (Kokubo & Ida 2002; Leinhardt & Richardson 2005). Eccentricities and inclinations were drawn from a Rayleigh distribution with dispersions 〈e²〉 = 2〈i²〉 = 0.007 (Ida & Makino 1992). Table 1 gives individual simulation details including the initial embedded planet mass, M_pl(0), the planetary growth time from M_pl(0) to one Jupiter mass (M_J), T_grow, the eccentricity of the planet, e_pl, and the final stop time in each simulation, T_final.

Zoom In Zoom Out Reset image size

Table 1. Numerical Simulation Initial Conditions

Simulation	Collision Model	Mpl	Tgrow	epl	Tfinal
		(0)	(yr)		(yr)
0	merging	1 MJ	0	0.0	1000
1	merging	1 MJ	0	0.1	1000
2	merging	1 MJ	0	0.2	1000
3	merging	1 MJ	0	0.3	1000
4	merging	15 M⊕	1000	0.0	1000
5	merging	15 M⊕	1000	0.1	1000
6	merging	15 M⊕	1000	0.2	890
7	rubble	15 M⊕	1000	0.0	824
8	rubble	1 M⊕	1100	0.0	933
9a	merging	15 M⊕	460	0.0	1000

Note. ^aThe disk in this simulation extends from 0.8 to 10 AU.

Download table as:  ASCIITypeset image

Computational constraints placed a practical limit of 1000 yr on the simulations. We investigated two initial mass scenarios for the embedded planet: (1) Since the final stage of gasgiant planet growth is short (∼10³–10⁴ yr; Rice & Armitage 2003; Dodson-Robinson et al. 2008), the planetesimal disk would have little time to respond to the gravitational influence of the planet. Thus, in simulations 0–3 we began with a fully grown 1M_J planet; (2) to simulate more realistic scattering from a growing giant planet we also ran several simulations with an embedded planet that increased in mass exponentially from a terrestrial core to 1M_J over ∼10³ yr (simulations 4–8), and 460 yr (simulation 9). To grow the planet from a terrestrial-mass embryo to M_J,

$$\begin{equation} dm = \lambda M(t)\, dt \end{equation} \tag{ 1 }$$

was added at each time step, dt, to the embryo which had a mass of M(t) at time t, where λ = ln(M/M_pl)/T_grow is the growth constant (Dodson-Robinson et al. 2008).

The average eccentricity of Jupiter-mass exoplanets between 1.5 and 10 AU from their host star is 0.28.⁴ Thus, the eccentricity of the embedded planet was varied from 0.0 to 0.3. A more eccentric planet would create a more dynamically excited and collisionally active planetesimal disk, but this may result in a decrease in the planetesimal growth rate and loss of significant mass through collisional grinding. Therefore, it is initially unclear which orbital configuration produces more second-generation dust.

All simulations are in a gas-free disk for a number of reasons: (1) planetesimals in these simulations are all initially large (R ∼ 150 km) so aerodynamic gas drag over the simulation timescale is insignificant; (2) there is observational evidence that the inner hole of a transitional disk has low gas density (Pontoppidan et al. 2008; Brown et al. 2008; Ireland & Kraus 2008); (3) it allows us to isolate the effects of collisional dynamics; (4) in this work, dust is the instantaneous collisionally generated second generation dust, not primordial dust that is intimately mixed with the gas.

All simulations used a constant orbital (major) time-step of 0.01 yr, providing good temporal resolution. However, minor steps were used in resolving the details of a planetesimal collision in the rubble model.

2.2. Planetesimal Collision Models

2.3. Calculating the Dust Distribution

The rubble model includes debris production. However, the perfect merging collision model does not. In a previous work on debris disks, Booth et al. (2009) assumed that any collisional debris followed the mass. In his work each resolved planetesimal represented the massive end of a Dohnanyi collisional cascade (Dohnanyi 1969). However, in our simulations the planetesimal disk is not in a steady state; instead it is rapidly evolving as the young solar system develops. Therefore, a Dohnanyi collisional cascade is not applicable. Instead we assume an instantaneous size distribution n(r)dr = ζr^−3.5dr, where ζ is a normalizing factor, n(r) is the number of objects with a radius between r and r+dr, for the collision debris (Leinhardt & Stewart 2012). The purpose of our algorithm is to obtain a realistic instantaneous spatial distribution of the collisional dust. This will depend on the collisional environment of the planetesimals. Below we describe the algorithm:

Step 1. Determine the total mass crossing the orbit of each planetesimal, M_crossing. The orbits of two planetesimals intersect under the condition

$$\begin{equation} P_1 + P_2 \le \sqrt{A^2 + B^2}, \end{equation} \tag{ 2 }$$

where P is the orbital period, A = e₂P₁cos (ψ₂) − e₁P₂cos (ψ₁), B = e₂P₁sin (ψ₂) − e₁P₂sin (ψ₁), ψ is the phase of the orbit.

Step 2. Use the total mass and number of orbit crossing planetesimals (N_crossing) to calculate the average mass of a collider,

$$\begin{equation} m_{{\rm avg}}=M_{{\rm crossing}}/N_{{\rm crossing}}. \end{equation} \tag{ 3 }$$

Step 3. Find the number of collisions necessary to encounter one target mass worth of colliders,

$$\begin{equation} N_{{\rm coll}} = m_{{\rm targ}}/m_{{\rm avg}}, \end{equation} \tag{ 4 }$$

where m_targ is the mass of the target body.

Step 4. Scale the mass of a full size distribution of debris, m_smash, by the number of target masses encountered per orbit,

$$\begin{equation} \alpha _{{\rm coll}} = \frac{P}{N_{{\rm coll}} t_{{\rm coll}}}, \end{equation} \tag{ 5 }$$

where t_coll = 1/(χvn_avg) is the collisional timescale of the target, $\chi = \pi (R+r_{{\rm avg}})^2 (1 + v_{{\rm esc}}^2/\sigma _v^2)$ is the collisional cross-section of the target, R is the target radius, r_avg is the radius of m_avg, v_esc the escape velocity from an object with a mass equal to the combined mass, m_targ + m_avg, σ_v is the velocity dispersion of the orbit crossing planetesimals, v is the speed of the target, and n_avg is the number density of colliders.

Step 5. Scale the mass of debris produced from the total disruption of the target, m_smash, by α_coll and calculate the mass of a single dust pixel,

$$\begin{equation} m_{{\rm dust\,pixel}} = \alpha _{{\rm coll}} m_{{\rm smash}}/ N_{{\rm dust\,pixel}}, \end{equation} \tag{ 6 }$$

where N_dust pixel is the number of dust units to be spread around the target's orbit (in this work N_dust pixel = 1000).

Step 6. Spread the dust around the target orbit in θ such that each dust pixel is spaced equally in time. Therefore,

$$\begin{equation} d\theta = 2 \pi a b/(N_{{\rm dust\,pixel}} r^2), \end{equation} \tag{ 7 }$$

where a is the semi-major axis of the orbit, b is the semi-minor axis of the orbit, and r is the distance from the stellar focus to the dust pixel.

The fixed-mass embedded planet simulations (simulations 0–3) are qualitatively similar to the corresponding growing-planet simulations (simulations 4–9) in all major respects (i.e., planetesimal, surface density, and dust distributions). However, the large impulse from the 1M_J planet at the start of the simulations reduces the disk particle number drastically, giving poor collision statistics. Therefore, only simulations with a growing planet (simulations 4–9) are discussed here.

Figure 1 shows an example of the temporal evolution of the planetesimal surface density (in this case for simulation 4). The structures in the surface density evolution are qualitatively similar in both the equivalent rubble simulation (simulation 7) and the fixed-planet mass simulation (simulation 0). All of the e = 0.0 simulations result in a clearly defined gap around the orbit of the planet and additional partial gaps at low-order mean-motion resonances. In addition there is a horseshoe ring of planetesimals that orbit with the planet in 1:1 resonance (visible from 330 yr). These obvious structures in the planetesimal disk become less clear with increasing planet eccentricity.

The collisionally active zones of the planetesimal disk are highlighted in Figure 2. Significant planetesimal growth is generally confined to the inner planetesimal disk with more modest growth in the outer disk and the 1:1 resonance. In the low-eccentricity simulations, regions of growth are stimulated by low-eccentricity pileups such as those at the 1:2 and 2:3 inner resonances. However, it is clear from the extended planetesimal disk simulation (simulation 9; Figure 2(d)) that planetesimal growth is enhanced out to at least the outer 2:1 resonance. Unlike the smaller disk simulations (simulations 0–8), there are enough planetesimals at low eccentricity and larger semi-major axis to collisionally interact with the highly eccentric planetesimals that have been excited by the planet (those planetesimals on eccentric spurs—sharp arc-like features), increasing dust production both interior and exterior to the planet orbit. In the high-eccentricity simulations, dynamical stirring increases the collision rate over a broad range in the semi-major axis as in Figure 2(b). All simulations show eccentric spurs which increase the velocities and probabilities of collisions and are the main mechanism by which planetary stirring enhances collisional dust production.

Zoom In Zoom Out Reset image size

Figure 3(a) shows the surface density in dust if dust mass is tied directly to the planetesimal mass as is assumed in debris disks (Booth et al. 2009). Figures 3(b)–(e) use our new method described in Section 2.3 which should be more realistic. If we assume that these disks are disruptive and in a quasi-steady state of collisional grinding, then the inner disk will dominate the dust surface density. However, the reason that the inner disk is so bright in Figure 3(a) is because the planetesimals at a small semi-major axis are not significantly perturbed by the embedded planet and have grown to large masses. In the Booth et al. (2009) method, each planetesimal represents an entire collisional cascade, which generates significant dust mass. In reality, the inner region of the simulations is dynamically quiet and produces only a small amount of collisional debris because most collisions result in accretion events not disruption events. Thus, for the rest of this paper we will discuss only the instantaneously generated dust mass calculated with our algorithm in the context of young transitional objects.

Zoom In Zoom Out Reset image size

Simulations with a circular embedded planet (Figures 3(b) and (c)) have a bright inner ring, a dark gap with bright 1:1 resonance, and a bright outer ring. However, the simulations with an eccentric planet (e > 0.0; Figures 3(d) and (e)) have less pronounced 1:1 resonances and an incomplete gap. The fixed-mass simulations 0–3 share these qualitative traits. The murky gap in the eccentric simulations is at odds with the predictions from fluid simulations that predict a cleared eccentric gap (Marzari et al. 2010; Moeckel & Armitage 2012). However, the short timescale of our simulations may be the cause of this discrepancy.

Flux density plots in the bottom row of Figure 3 were created with RADMC3D (Dullemond 2012), a Monte Carlo radiative transfer code, at a wavelength of λ = 850 μm. The total observed flux depends strongly on the dust mass, which in turn depends on the assumed disk mass. We do not make any attempt to predict flux but instead concentrate on the qualitative features expected in a disk that is initially optically thick. Small dust was modeled as two dust species, amorphous carbon and silicates, both species having radii of 0.1 μm and 0.631 μm, and relative abundances of 0.2 and 0.8, respectively. Larger dust from 1 μm to 1000 μm is modeled using simple Mie scattering spheres with three size bins per decade. This range was chosen as dust of all sizes should be created in a collision. The smallest size is approximately defined by the dust blow-out radius, and the largest is when emission in the sub-mm range is no longer significant. Features highlighted in the surface density plots are still present. In addition pericenter glow is visible in Figures 3(h) and (i) along with an azimuthal asymmetry in the flux distribution of the inner disk.

Figure 4 shows an illustrative example of the observational outcomes of our work. Assuming a dust mass of 2.8 × 10⁻⁷M_☉ (the mass of the debris produced in simulation 7) and the dust composition outlined above, the maximum signature is ∼200 mJy when observed face on at 10 pc. As this is derived from the surface density, this depends on: the size range of the dust, the extrapolation model from planetesimals to dust, and the particle number. However, the disk gap would be visible up to a beam size of ∼0.5 AU (005) and the asymmetry of the disk when e > 0.0 up to a beam size of ∼2 AU (02).

Zoom In Zoom Out Reset image size

In this work we present results from ten numerical simulations of collisional dust production from dynamical stirring by a young embedded planet. The simulations that we present use two different planetesimal collision models, perfect merging and rubble. Most simulations use a small 3 AU wide planetesimal disk and 10⁶ planetesimals, but we have also completed one lower resolution large disk simulation with a 10.2 AU wide disk and 10⁵ planetesimals. A significant advance of this work is the use of a more realistic dust production model, increasing the accuracy of the spatial distribution of instantaneous collisionally produced dust and the resulting observables.

The series of simulations found that dust production is enhanced by the presence of a planetary companion. Dynamical stirring from the planet perturbs some planetesimals onto highly eccentric orbits (Figure 2) which increases collisional activity from the inner 1:3 resonance to the outer 2:1 resonance producing two bright regions interior and exterior to the gap (Figure 3). The eccentric simulations show that a clean gap is not opened in contrast to the fluid case (Marzari et al. 2010; Moeckel & Armitage 2012), but the disk asymmetry can be easily observed (Figure 4).

The next steps will be to test and improve this model in numerous ways, observing the effect of multiple planets, improving the treatment of dust dynamics, and taking account of the dust lifetime.

The algorithm in Section 2.3 neglects the effects of gas and primordial dust on the dust dynamics. This work addresses the extreme case of no gas, but as the density of gas in the inner regions of transitional and pre-transitional disks is not well constrained (Pontoppidan et al. 2008; Brown et al. 2008; Ireland & Kraus 2008), other scenarios should be investigated.

The rubble model is computationally intensive in a collisionally active disk. In the future EDACM (Leinhardt & Stewart 2012) will be used. EDACM has the same fragmentation and debris tracking capabilities whilst using analytical expressions to determine outcomes, reducing computational load.

Finally, the exact magnitude and lifetime of the collisional dust cannot be determined from these simulations. Further investigation into this area will have to treat the dynamics and lifetime of collisional dust carefully.

J.D. is grateful for support from NERC, Bristol School of Physics, and McDonald Observatory. Z.M.L. is supported by a STFC Advanced Fellowship. S.D.R. is supported by National Science Foundation CAREER award AST-1055910. N.A.T. is supported by the Leverhulme Trust. Our thanks to Andy Young for help with data visualization. This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol.⁵ We thank an anonymous referee for thoughtful comments that greatly improved the manuscript.