New Spatially Resolved Mid-IR Observations of the Transitional Disk TW Hya and Tentative Evidence for a Self-Luminous Companion
We present spatially resolved observations of the canonical transition disk object TW Hya at 8.74 m, 11.66 m, and 18.30 m, obtained with the T-ReCS instrument on the Gemini telescope. These observations are a result of a novel observing mode at Gemini that enables speckle imaging. Using this technique, we image our target with short enough exposure times to achieve diffraction limited images. We use Fourier techniques to reduce our data, which allows high-precision calibration of the instrumental point spread function. Our observations span two epochs and we present evidence for temporal variability at 11.66 m in the disk of TW Hya. We show that previous models of TW Hya’s disk from the literature are incompatible with our observations, and construct a model to explain the discrepancies. We detect marginal asymmetry in our data, most significantly at the shortest wavelengths. To explain our data, we require a model that includes an optically thin inner disk extending from 0.02 to 3.9 AU, an optically thick ring representing the outer disk wall at 3.9 AU and extending to 4.6 AU, and a hotter-than-disk-equilibrium source of emission located at 3.5 AU.
Current planet formation theories are informed largely from observations of evolved planetary systems. Since the discovery of the first extrasolar planets (Mayor & Queloz, 1995) and the subsequent explosion of the field, theories of planet formation have grown and proliferated. Though the analysis of evolved systems is useful, perhaps the most effective way to test planetary formation theories is through the direct observation of young systems in the midst of planetary formation.
Objects thought to be in this stage, so-called transitional disks (Strom et al., 1989), are intermediate between embedded pre-main-sequence stars and evolved planetary systems. These objects exhibit evidence of a low-density, optically thin gap in the innermost regions of their primordial circumstellar disks. The presence of this gap was originally inferred from the near- and mid-infrared (mid-IR) properties of the spectral energy distributions (SEDs). The SEDs of these objects exhibit a deficit of flux at wavelengths less than about 10 m compared to the SEDs of classical T Tauri stars (CTTS), and a flux excess resembling CTT SEDs at wavelengths of 20-60 m. The deficit presumably occurs because the matter—and hence flux—in the inner regions of the disk is absent. The excess is caused by the outer disk, extending from the puffed up “edge” or “wall” at several AU to the extent of the disk at tens or hundreds of AU.
Though their geometries were originally inferred from the properties of their unresolved SEDs, transition disk objects have since been resolved at various wavelengths. Thalmann et al. (2010) analyze - and -band data (1.65 m and 2.2 m, respectively) and resolve the gap in the transition disk LkCa 15. They confirm a disk with a 46 AU truncation radius, inside of which is a largely evacuated gap. This was also confirmed by Andrews et al. (2011a) with imaging of LkCa 15 at 870 m. Several other transition disks have been resolved at millimeter wavelengths (Hughes et al., 2007; Isella et al., 2008; Brown et al., 2009; Isella et al., 2009, 2010; Andrews et al., 2011b).
The processes driving the formation of these gaps is not yet fully understood, though several hypotheses have been proposed (photoevaporation, grain growth, or a stellar or planetary companion). Many recent observations have been interpreted as evidence of a companion, supporting the hypothesis that young planets are important in the evolution and dissipation of the disk in these transition disk objects. Kraus & Ireland (2011) present the detection of a likely protoplanet located inside the known gap in LkCa 15. They use non-redundant aperture masking interferometry at three epochs to discover a faint companion located at 16-21 AU from the primary star. Huélamo et al. (2011) detect a source in the band (3.5 m) at a separation of 6.7 AU—within the disk gap—from the transition disk T Chamaeleontis. Eisner et al. (2009) present mid-IR observations of the transition disk SR 21. They spatially resolve the dust emission around this object, and suggest that the disk around SR 21 must be completely cleared within 10 AU. They propose a disk with a large inner hole and a warm companion near the outer edge of the cleared region.
TW Hydrae (TW Hya), a classical T Tauri Star, is one of the closest transition disks, located in the nearby (50 pc; Mamajek, 2005) TW Hya Association. Due in part to its proximity and age (10 Myr; Webb et al., 1999), TW Hya has become an intensely studied object. Calvet et al. (2002) first proposed that TW Hya had a developing gap located at 4 astronomical units (AU), inferred from various features of the SED. More recently, spatially resolved observations of TW Hya have offered additional insight. The outer radius of the disk has been resolved in millimeter wavelengths and determined to be in the range of AU (Wilner et al., 2000, 2003; Qi et al., 2004; Isella et al., 2009). Hughes et al. (2007) resolved the inner edge of the gap in the circumstellar disk at 7 mm, corroborating Calvet et al. (2002)’s determination of a hole at 4 AU. Published contemporaneously with Hughes et al. (2007), Ratzka et al. (2007) contradicted the claims of Calvet et al. (2002) in their analysis of interferometric observations of TW Hya in the mid-IR with the Very Large Telescope Interferometer (VLTI). They suggested that the disk gap occurs at a considerably smaller radius: AU. At yet shorter wavelengths, Eisner et al. (2006) measured the -band (2 m) visibility at a single baseline with the Keck Interferometer (KI), inferring an inner radius for the optically thin, evacuated region of AU.
Most recently, Akeson et al. (2011) presented near-infrared (near-IR) measurements from the Center for High Angular Resolution Astronomy (CHARA) array and the KI at various baselines. They model many of the past observations of TW Hya simultaneously with theirs. In combining the SED and spatially resolved observations from the literature (Eisner et al., 2006; Hughes et al., 2007; Ratzka et al., 2007), they review and extend past attempts at modeling the circumstellar disk, producing their own hybrid disk geometry. The Akeson et al. (2011) model is essentially the original Calvet et al. (2002) model with an added optically thick ring of emission at 0.5 AU, at roughly the disk equilibrium temperature (see Section 3).
Modeling of spatially unresolved data for TW Hya (e.g. SEDs or higher resolution spectra) can yield very different conclusions. Spatially resolved observations are necessary to properly constrain the disk geometry. It is important to note the ambiguity associated with modeling source geometries from the SED alone: Boss & Yorke (e.g., 1993, 1996) found that the interpretation of infrared (IR) flux deficits as central gaps does not offer a unique solution, and that opacity and geometrical effects produce degenerate solutions fitting the SED of an unresolved system. We too find that this is true, and that the choice of dust species and thus inner disk opacity can have a significant effect on the determination of disk geometry. Spatially resolved observations, ideally at multiple wavelengths, are needed to unambiguously constrain disk geometry.
In order to further constrain models of TW Hya’s circumstellar disk, and to resolve discrepant interpretations of observations, we present new speckle interferometric observations in the mid-IR. We observe TW Hya at three wavelengths (8.74 m, 11.66 m, and 18.30 m) and at two epochs (2007, 2009) using the Thermal-Region Camera Spectrograph (T-ReCS) instrument on the Gemini telescope. We resolve the disk at each wavelength, and construct a simple model to fit our observations. Combining our results with previous spatially resolved imaging and unresolved spectro-photometry, we constrain the properties of the inner disk (or lack thereof) in TW Hya.
Drawing from the literature, we first attempt to reproduce simple versions of the models first proposed by Calvet et al. (2002); Uchida et al. (2004), and Ratzka et al. (2007). In the literature we identify two general classes of models: “Calvet-like” models with an optically thick wall at 4 AU with a hole of optically thin material within (Calvet et al., 2002; Uchida et al., 2004; Hughes et al., 2007; Akeson et al., 2011; Gorti et al., 2011); and “Ratzka-like” models for which there exists a similar opacity hole, but located much closer in at 1 AU (Ratzka et al., 2007; Akeson et al., 2011). We approximately reproduce both model types, and compute synthetic mid-IR visibilities for comparison with our data; these are presented in Section 3.1. In Section 3.2, we describe a new model consistent with both the data from the literature and our new observations.
We present four observations of the transitional disk object TW Hya at three different wavelengths and two different epochs. A summary of these observations is presented in Table 1. The wavelengths and epochs of the observations are: 8.74 m (8.74) and 11.66 m (11.66a), obtained in 2007; and 11.66 m (11.66b) and 18.30 m (18.30), obtained in 2009. For each observation, we also took off-target nod pointings for infrared sky background subtraction. Each nod, or pointing, consists of several short exposures (“frames”) that are combined to yield a high resolution speckle image. It is the individual frames, statistically combined, that are ultimately used in our analysis. For each dataset, we observed a bright, point source calibrator in order to calibrate and remove instrumental and atmospheric effects.
We collected our data using a custom observing mode on the T-ReCS instrument at Gemini. We used short integration times (172 msec) to freeze the motion of the atmosphere and thus achieve high spatial resolution, enabled by the diffraction limit of the telescope. We employed a dither / nod pattern that moved observed objects around four positions on the detector. For the 18.30 m dataset, we employed standard chopping, and our frame integration times were a factor of 10 longer than at the other wavelengths. Total integration times are included in Table 1.
The fluxes for our SED comparisons were obtained from the Spitzer Space Telescope Infrared Spectrograph, from Uchida et al. (2004) (see their Figure 2). For our calculations of our best-fit model, we increase their stated errors of 10% to 30% to account for the larger number of points in the SED as compared to our resolved dataset, and to weight higher the value of the resolved data over the unresolved data in our fits. The errors displayed in the Figures are Uchida et al. (2004)’s 10%. We also use mid-IR interferometric measurements for our analysis. These data were obtained with the MIDI instrument at the VLT, and were presented in Ratzka et al. (2007).
|9 May 2007||II Hya||8.74||173||18||702||351||61|
|9 May 2007||II Hya||11.66||173||30||1170||1131||195|
|9,18 April 2009||HD 92036||11.66||181||100||3400||2505||454|
|22,23 May 2009||HD 92036||18.30||1813||76||228||228||413|
Note. – The data used in our analysis. We observed the target, TW Hya, at three different wavelengths and at two different epochs. We discard unusable frames, according to criteria in Section 2.2. Columns are: (1) date of observation; (2) name of calibrator star; (3) wavelength of observation in microns; (4) integration time of an individual frame in milliseconds (msec); (5) number of nods in the observation of the target; (6) number of individual frames in the observation of the target; (7) number of frames used for analysis, after removing flagged frames; and (8) total integration time for target using all usable frames, in seconds (sec).
2.2 Reduction and Analysis
2.2.1 Fourier Analysis
We reduce our data using Fourier analysis techniques. Each individual frame has an exposure time that alone is too short to provide a significant signal to noise ratio. However, a single long exposure would yield an image with insufficient angular resolution due to the effects of the Earth’s turbulent atmosphere. One cannot naïvely add together the short, individual frames, as the target centroid shifts on the sky due to atmospheric turbulence. We instead combine the power spectra of individual frames. This method of addition is independent of translations of the image centroid position. Furthermore, analyzing the data in Fourier space allows us to remove the instrumental point spread function (PSF) with a point-source calibrator with high accuracy. In Fourier space, we divide the power spectrum of the source by the power spectrum of a point source calibrator. This is considerably simpler than deconvolution in the image plane.
The first step in our data reduction procedure is to remove the sky background. We obtained our data in such a way that each series of frames (short exposure images) of a target is paired with a slightly offset series of frames of the same target; we call these two sets “adjacent nods”. This offset causes the target to appear on two different regions of our detector. We use these adjacent nod pairs to overcome the high sky background in the mid-IR by subtracting from each target its adjacent nod. As each nod contains many frames, the median value of all the frames in a particular nod is used for the background subtraction. The median value of an adjacent nod is subtracted from each individual short-exposure frame. For example, in a single nod pair, the median value for all the frames in Nod Type “A” is subtracted from each frame in Nod Type “B”, and likewise the median value of “B” frames is subtracted from frames in “A”. This process is repeated for each pair of nods.
Next, we flag unusable frames, without actually modifying any of the data used for analysis. We first attempt to locate a bright point source in the median image of a pointing: either the target—TW Hya—or a calibrator. To locate the position of the star in a pointing (or nod), we calculate the median image for that nod, subtract from the median image the mean value of all the pixels, and set negative pixel values to zero. We identify “hot pixels” as pixels with values more than three standard deviations larger than the mean pixel value of the resultant image. These “hot pixels” are assigned a value equal to the average value of their nearest neighbors using image convolution. This technique removes isolated “hot pixels” while leaving signal from a star relatively unaffected. The position of the star is then set to the location of the maximally valued pixel in the image. We classify the pointing as unusable if, after this process, there exist no pixels that are 5- deviations from the mean, or if the location of the point source determined by our algorithm is closer than the width of our subimages to the detector edge (see next paragraph). If we fail to detect a bright point source in a particular pointing, we flag all the frames in that pointing as unusable. This could happen due to poor seeing in these exposures, the presence of clouds, or some other effect. We note again that the steps described in this paragraph do not ultimately affect the data used for analysis, but only for data flagging. After a star is located, we return to the background-subtracted data described at the end of the previous paragraph.
After eliminating flagged data, we are left with frames containing usable point sources. For each frame, we cut out a “postage-stamp” size subimage from the raw telescope image, centered on the point source. These subimages are 64x64 pixels, or 5.765.76 arcseconds at T-ReCS’s plate scale, or 309309 AU at the distance of TW Hya. We are then left with order hundreds (ranging from to ) of short-exposure frames for target and calibrator for each dataset, cropped and centered at the star’s position.
For each frame, we subtract a residual sky value in addition to the initial IR sky background subtraction: the median of pixel values in the outer regions ( pixels or farther from the image center of each subimage). We then apply a Hanning window to the image by multiplying by the two-dimensional Hanning function:
where is the size (diameter) and is the center of the Hanning window. For our analysis, pixels (roughly arcseconds or AU). The Hanning window has some useful properties: it is unity at the center and falls smoothly to zero at the edges. This step removes spurious Fourier signals by smoothing sharp brightness transitions. We take the Fourier transform of the Hanning windowed image, and record the squared amplitude at each pixel (two dimensional power spectrum, or power image). We perform steps identical to those described above for an adjacent, equally sized, background-subtracted subimage of blank sky from the same parent image to produce a sky power image.
After computing the Fourier amplitudes, we use sigma-clipping to discard discrepant data. For each dataset, we calculate the mean and standard deviation of the squared visibility as a function of baseline for both target and calibrator. If any individual frame deviates 3- from the mean at any moderate, well-behaved baseline (between 0 m and 5 m), we exclude it from the sample. This step excludes a small fraction of our resulting frames: usually zero, but occasionally as high as 2.
After the reduction procedure outlined above, we are left with several hundred power images, for both target and calibrator. Additionally, we have power images for adjacent, blank sections of sky from parent images for each of these postage-stamp images. We first remove any bias due to detector artifacts by subtracting the sky power image from the object power image. We then subtract a residual bias: any non-zero value for the visibility amplitude at baselines well beyond our sensitivity (9 m and longer). We set negative values of the target’s power images to zero, and negative values of the calibrator to a small value (; to avoid eventual division-by-zero computational errors). We then normalize the target and calibrator power images by the maximum value in each respective set.
We calculate an azimuthal average of each power image by computing the mean of the pixels in an annulus corresponding to a particular baseline length. The power obtained in each annulus becomes the visibility amplitude at a corresponding baseline. We do this to increase the signal to noise and to make our plots more straight-forward to interpret. Furthermore, to rough approximation, it is an acceptable assumption that our target—a star with a nearly face-on circumstellar disk—is indeed symmetric about the azimuth. Calvet et al. (2002) also assume a face-on orientation, in agreement with estimates of TW Hya’s inclination in the literature (Krist et al., 2000; Qi et al., 2004; Pontoppidan et al., 2008). We discuss possible departures from our assumption of azimuthal symmetry in Section 3.3.
We use the azimuthally averaged visibility curves (one for each frame) to calculate a mean and standard deviation of the mean (SDOM, or standard error) for target and calibrator. We then use these derived values (, , , ) to calculate the calibrated visiblity () and its associated statistical uncertainty. The SDOM is calculated by dividing the usual standard deviation by the square root of the number of frames ().
The total statistical uncertainty associated with the calibrated visibility is given by the propagation of the two SDOMs. We calculate the propagated error using the usual relation,
for and for uncorrelated errors in and .
We also estimate a systematic error for our experimental setup, by comparing calibrator observations. For example, variation in seeing or sky background structure over timescales longer than our source-calibrator cadence will lead to systematic errors. Similar to the procedure described above, we average the azimuthal average of power images according to telescope pointing (i.e., we group calibrator frames of a single pointing together). We divide pointing-averaged calibrator visibility curves by the mean visibility curve of all the calibrators. We use the SDOM—not the standard deviation—of the several calibrator pointings as a systematic error: the magnitude at which we expect an uncertainty to exist for a given set of observations. For each baseline, the visibilities of the calibrators do appear to be normally distributed. We add this error in quadrature with the statistical uncertainties described above. The final error in our measurements is dominated by these estimated systematic errors. We present the data used for this method of systematic error estimation in Figure 1. This plot shows the visibility amplitude of each calibrator pointing, divided by the mean visibility amplitude of all the calibrators, as a function of baseline, for each dataset.
We expect the calibrator errors to appear correlated to some extent. That is, if a calibrator has a larger visibility amplitude at one baseline, it likely has a larger visibility at other baselines as well, as seeing variations are an important cause of the dispersion. One other possible cause for correlated errors is our application of a Hanning windown in the data reduction process. To check this effect, we perform our analysis without the Hanning window; we notice a slight decrease in apparent correlation. The size of the calculated error bars did not change appreciably, however, and we keep the Hanning window in our reduction process.
2.3 Size of Emitting Region in TW Hya
We estimate the size of the emitting region at each wavelength by fitting a ring of fixed thickness (width to radius ratio, 0.11) to each visibility dataset (see Equation 6). We performed a chi squared minimization where the size of the ring of emission was the only free parameter. Error bars (1-) were derived by finding the ring size value that corresponded to the minimum chi squared plus one. The results of our fit can be seen in Table 2 and a comparison of our data with the best fit ring visiblity curves can be found in Figure 2.
Note. – The results of a single wavelength ring fit to our observations. Error bars are derived by finding the value (in radius) on the chi squared contour that corresponds to one plus the minimum chi squared value. Columns are: (1) name of observation; (2) best fit radius with - error bars; (3) the ratio of ring width to inner radius, set to a constant for these illustrative purposes. We note the 2.5- difference in size between our two epochs of 11.66 m observations.
Since we have two epochs at 11.66 m, we are able to constrain variability over a timescale of roughly two years (720 days). Variability is observed between these two epochs at the 2.5- level (Figure 2; Table 2). For simplicity, we use both epochs of 11.66 m data in our modeling. We speculate as to the cause of this variability in Section 4.
We model the TW Hya system using largely geometric models. We approximate the central star as the Kurucz stellar atmosphere of a K7 star with radius , temperature K, surface gravity log (cm s) = 4.5, at a distance pc (Webb et al., 1999; van Leeuwen, 2007). We note, as has been mentioned previously in the literature (e.g., Sitko et al., 2000), that TW Hya is a known variable star, so we do not necessarily expect a Kurucz model to exactly model the emission of the star. Also, the shorter-wave flux density measurements from the literature were not taken contemporaneously with the mid-IR spectrum that comprises the bulk of the disk emission. While this Kurucz model does not perfectly reproduce the stellar flux at all wavelengths, and different analyses of TW Hya use slightly different stellar model parameters, it is sufficient for our modeling purposes.
Our models of TW Hya are composed of an optically thin disk (approximated as a series of concentric rings which follow a known temperature-radius relationship) and an optically thick ring of a single temperature representing the directly illuminated edge of the optically thick disk. For the optically thin component, we must choose and apply dust opacities. We discuss the details of this choice below.
3.1 Existing Models
We first attempt to reproduce the large cavity model of Calvet et al. (2002); Uchida et al. (2004). We use geometric model parameters and dust opacities from their work. We extract dust optical depth directly from Uchida et al. (2004)111To obtain the optical depths used by Uchida et al. (2004), we first record the reported flux density from their optically thin disk component (see their Figure 2, top panel). We then use these published flux densities, an assumption of constant surface density (as in Calvet et al. (2002)), and their stated disk geometry, to extract values for the dimensionless wavelength-dependent disk optical depth, , that they use., which presents results consistent with those of Calvet et al. (2002). This model consists of emission from three components: 1) a star, 2) an optically thin disk, extending from the dust sublimation radius at 0.02 AU to the directly illuminated disk wall at 3.3 AU, and 3) an optically thick ring at 3.3 AU, representing the puffed or flared wall of the outer disk edge. To approximate Calvet et al. (2002); Uchida et al. (2004)’s model of the stellar flux, we allow the flux density produced by the Kurucz stellar model to vary so that it matches the observed flux density at m; this is similar to their own strategy and the procedure of Sitko et al. (2000).
The flux from the optically thick disk wall is approximated by an unmodified blackbody. The flux density of an annulus of thickness is given by:
In the case of the optically thick ring, for all values of . For the optically thin disk, the temperature profile is given by
where is the temperature of the star, is the radius of the star, and is stellocentric radius. To calculate the flux from the disk, we use small concentric annuli. The temperature structure at every point in the optically thin disk behaves according to Equation 4, and the fluxes from each annulus are summed together to obtain the total disk flux. The temperature for the directly illuminated rim or edge of the optically thick disk is given by
We are motivated in the choice of these temperature profiles by the methods of Chiang & Goldreich (1997); Calvet et al. (2002); Eisner et al. (2006) (see Equation 1, Equations 1 and 3, and Equation 1, respectively).
All model components can be viewed as collections of single-temperature rings that obey the model temperature profile. Model fluxes are computed by summing Equation 3 over all annuli. The squared visibility for a ring (or annulus) is given by (see Eisner et al., 2003):
where is the - radius, is the wavelength, is the angular diameter of the object in radians, defines the ratio of radial thickness to size of the annulus (, where is the width of the annulus and is the inner radius of the annulus), is the Bessel function of the first kind of order one. The normalized visibility for the entire model is the flux-weighted average of the visibilities produced by each annulus.
From these model parameters (matched to those used in Calvet et al. (2002); Uchida et al. (2004)), we generate a synthetic SED and visibilities (Figure 3). The large bottom panel shows the SED, indicating the stellar flux density (thin solid line), the optically thin, low density inner disk (dotted-dashed line), the outer disk wall (dotted line), and the total flux (thick solid line). The flux data are also shown (blue circles). In the top half of the figure, we show our four separate observations: 8.74 (blue, top left), 11.66a (black, top right), 11.66b (black, bottom right), and 18.30 (red, bottom left). In each, the flux weighted visibilities of the different model components are indicated by different line types, as in the SED plot. The insets in the top right of each panel show zoom-ins of the data presented in this work, while the long baseline points at 45 m are from Ratzka et al. (2007). Vertical, dashed lines in the SED panel show the wavelengths of our observations at 8.74 m, 11.66 m, and 18.30 m.
Our estimate of Calvet et al. (2002)’s model reproduces the SED as well as done by the original authors. This is not unexpected, as we obtained the opacities, physical geometries, and temperatures directly from their work. More interesting are the mid-IR visibilities produced by this model, shown in the top half of the figure. This model reproduces well the 18.30 m data, a tracer of cooler portions of the disk, dominated by emission of the transition disk wall. Similarly, this model reproduces rather well both epochs of our 11.66 m data. One epoch (1166b) is reproduced better than the other (1166a), but we note that these two datasets are in fact inconsistent with each other due to temporal variability (see Section 2.3.1). This model is clearly inconsistent, however, with our very resolved 8.74 m observations, as seen in the top left panel of Figure 3.
3.1.1 Small Cavity Models
Though subsequent work has confirmed the validity of Calvet et al. (2002)’s model (e.g., Hughes et al., 2007), Ratzka et al. (2007) claimed that the evacuated cavity was much closer to the star, at 1 AU. Since Ratzka et al. (2007)’s model was based largely on mid-IR flux measurements, and presented resolved mid-IR data, we compare their work against our observations. Ratzka et al. (2007) used a more complicated disk modeling technique than did Calvet et al. (2002); we use our simple geometric models to distill out the most important elements. The main sources of the mid-IR emission in the model of Ratzka et al. (2007) were the optically thick, directly illuminated disk wall, and its optically thin atmosphere. We were thus motivated to produce emission at 1 AU from nearly coincident components, one optically thick and one optically thin.
Our methods to reproduce the Ratzka et al. (2007) model are similar to the process described previously in Section 3.1, with one exception. As the bulk of the optically thin emission from the Ratzka et al. (2007) model originates from the disk atmosphere, we use an optically thin ring at a single temperature (, ) instead of the optically thin inner disk described above (, ). The temperature of this ring is a model parameter and is free to vary. We also vary the surface density of the optically thin ring, . We vary these free parameters and minimize chi squared of the model compared to the mid-IR flux of TW Hya and the long-baseline (45 m) visibility points from Ratzka et al. (2007). As our motivation was to reproduce the model of Ratzka et al. (2007), we do not include our new visibility data in the minimization of chi squared.
Ratzka et al. (2007) used different optical depths than did Calvet et al. (2002); Uchida et al. (2004), and we obtained and used these in our modeling (T. Ratzka 2011, private communication). As in the modeling described above, we allow the flux density produced by the Kurucz stellar model to vary so that it matches the observed flux density at m; this is similar to their own strategy and the procedure of Sitko et al. (2000).
We performed a chi squared minimization and allowed the following parameters to vary freely to generate a model consistent with the Ratzka et al. (2007) geometry: Kurucz scaling, , , , , , . These properties and their best-fit values are listed in Table 3.
We present the calculated SED and visibility curves for our realization of this model in Figure 4. The large bottom panel shows the SED, indicating the stellar flux density (thin solid line), the optically thin, transitional disk wall atmosphere (dotted-dashed line), the transitional disk wall (dotted line), and the total flux (thick solid line). The flux density data are also shown (blue circles). In the top half of the figure, we show our four separate observations: 8.74 (blue, top left), 11.66a (black, top right), 11.66b (black, bottom right), and 18.30 (red, bottom left). In each, the flux weighted visibilities of the different model components are indicated by different line types, as in the SED plot. The insets in the top right of each panel show the data presented in this work, while the long baseline points at 45 m are obtained from Ratzka et al. (2007). Vertical, dashed lines in the SED panel show our wavelengths of observation.
Our realization of Ratzka et al. (2007)’s does not perfectly match the model SED presented in Ratzka et al. (2007) (particularly at m), but does a reasonably good job considering the complex disk properties (e.g., dust settling) we did not include. Our model visibilities do reproduce well those presented in Ratzka et al. (2007) at long baselines at 8.74 m and 11.66 m (see their Figure 8 for comparison). As can be seen in this figure, the Ratzka et al. (2007) data show the emission at 8.74 m to be more resolved than that at 11.66 m, though the model of Ratzka et al. (2007) does not reproduce this behavior. We see this trend of very-resolved 8.74 m emission in our own observations as well.
It is clear, however, that this small cavity model is inconsistent with nearly all our data. Though it is consistent with one epoch of our 11.66 m data (1166b), it is very inconsistent with all the other observations. It is unresolved at all baselines and all wavelengths relevant to this work. While we did not include an outer disk, this does not impact the observed discrepancy because the emission from these cool, outer regions is not relevant for the wavelength ranges we are interested in. This omission may have slightly worsened our reproduction of the small cavity model’s 18.30 m emission, but outer disk components at cooler temperatures will not cause the model to show an 8.74 m component as resolved as our data show.
We suspect that the difference in dust opacity choices is the principal reason that two distinct disk geometries—the 4 AU holes in Calvet et al. (2002) models, versus the 1 AU hole inferred by Ratzka et al. (2007)—can reproduce the same SED. In Figure 5, we compare the optical depths used by Ratzka et al. (2007) to those used by Calvet et al. (2002); Uchida et al. (2004). In order to examine the differences between the opacity choices in these two works, we have converted the optical depths from Ratzka et al. (2007) to a dimensionless optical depth with an assumption of surface density (informed from their work). We then scaled this value so that the mean values of optical depth for Calvet et al. (2002) and Ratzka et al. (2007) match, thus facilitating comparison. The Ratzka et al. (2007) optical depths show a relative paucity around m that explains why a smaller (1 AU versus 4 AU for a large cavity model) transitional disk rim is used in that model. We note that the optical depths used by Calvet et al. (2002); Uchida et al. (2004) were derived simultaneously with the other disk parameters, while those used by Ratzka et al. (2007) were not.
3.2 New Disk Model
Neither of the previous models considered (Calvet et al., 2002; Ratzka et al., 2007) adequately fit our new observations. The hybrid model of Akeson et al. (2011) also cannot explain our very resolved 8.74 m data. This is not surprising, as Akeson et al. (2011)’s model was designed to fit simultaneously the two model types proffered by Calvet et al. (2002); Ratzka et al. (2007). A new model is required. In order to produce more resolved emission at 8.74 m, without producing heavily resolved emission at larger wavelengths, we need a hot component at large stellocentric radius. We will show that this “hot” component must be considerably hotter than predicted by equilibrium disk models.
We adopt the optical depths from Uchida et al. (2004) (see Section 3.1). We assume a Kurucz stellar atmosphere with the scaling left as a free parameter. In addition to the components described in Section 3.1, we include an additional optically thick ring ( and ). We also vary the surface density of the optically thin inner disk, , but leave the dust sublimation radius used by Calvet et al. (2002); Uchida et al. (2004) unchanged, at 0.02 AU. We allowed the following parameters to vary freely (between physically reasonable limits) to generate our best-fit model: Kurucz Scaling, , , , , , , and . For definitions of these parameters, see Table 3.
After a thorough exploration of parameter space, we found a model that explains the existing data. The observables associated with this model are presented in Figure 6. In this figure we show the calculated SED and visibility curves for the model presented in this work. This model is a large cavity model, similar to Calvet et al. (2002) and consistent with the results of Hughes et al. (2007). The directly illuminated disk wall is represented by a ring of optically thick emission at 3.9 AU, and, as in Calvet et al. (2002), there is an optically thin disk of emission at low surface densities that extends from the dust sublimation radius to the outer disk wall ( AU). The crucial addition to this model is an additional optically thick ring of emission inside the disk wall, at 2.9 AU. In order to fit our 8.74 m data, this ring is at a hotter than equilibrium temperature (which would be 150 K, compared to 570 K for the hot ring). We also find that the best fit is achieved with an unscaled Kurucz stellar model, unlike the scalings of 2 required to fit the data using the models of Calvet et al. (2002); Uchida et al. (2004); Ratzka et al. (2007).
We note that with this model we produce a good fit to the SED, discrepant only in the 6-8 m region where possible stellar variability is a larger concern. Even so, the discrepancy in the SED fit is comparable in magnitude to that of the Uchida et al. (2004) deviation from observations in the 13 m region. We note that our fit to the long baseline mid-IR visibility amplitudes from Ratzka et al. (2007) is comparable to the model visibilities produced by Ratzka et al. (2007)’s own model. Most notably, however, we point out the simultaneous high quality fit of the long baseline mid-IR visibility amplitudes with the short baseline mid-IR visibility amplitudes presented in this work. Again, this was accomplished by forcing the ring to be hot enough to peak at shorter wavelengths—closer to 8.74 m—allowing the short baseline 8.74 m emission to be resolved, as we see in both our data and the VLTI data of Ratzka et al. (2007). Additionally, despite the increased degrees of freedom, the reduced chi squared value for our new disk model is 1% better than the value we calculate for the Calvet model. The reduced chi squared value for the Ratzka model is considerably larger. The value of the chi squared is largely driven by the unresolved epoch of our 11.66 m micron data, however. The less-resolved Calvet model fits this epoch slightly better. If we fit the average of the two 11 micron epochs, instead of attempting to fit both simultaneously, our new model has a reduced chi squared that is a factor of 3 lower.
3.3 Disk + Companion Model
The models described thus far are all symmetric about the azimuth. As described in Section 2.2.2, the reason we made this assumption was to decrease the uncertainty associated with our measurements of baseline-dependent visibility and to make our results more straight-forward to interpret. However, our solution of a hotter-than-equilibrium, geometrically thin ring of emission at a radius just inside the directly illuminated disk wall does not seem physical. That is, we cannot think of a physical process that could heat an annulus of matter far beyond the temperature achieved for a disk in equilibrium with stellar emission.
A self-luminous companion, on the other hand, can lead to such heating. Such an object would clearly not be azimuthally symmetric. The presence of a companion should thus manifest itself in our data through departures from azimuthal symmetry as a function of wavelength.
If all of our model components except the companion are azimuthally symmetric, we can predict relative degrees of asymmetry in the different wavelengths of observation. We expect the azimuthal asymmetry to be largest at 8.74 m; at 11.66 and 18.30 m, the flux contribution from the hot inner component is small compared to fluxes from other, symmetric model components (see Figure 6).
To test this, we create a model with a companion that has the same temperature and surface area as our added ring component and thus the disk plus companion model produces an SED identical to the model described in Section 3.2. For this reason we can utilize only our visibility data to constrain properties of the companion’s stellocentric radius and position angle (PA); the SED will not change as these parameters are varied.
We then calculate the degree of asymmetry generated by the presence of a companion in this model. Similarly, we measure the asymmetry in our reduced data. We follow the procedures outlined in Section 2.2 with one exception. Instead of performing an azimuthal average of the power image, we calculate the radial average as before, but include only small regions (of width 45 degrees) of azimuth separated by 23 degrees. That is, we calculate the radial average (power as a function of baseline) in different “wedges” of azimuth. We perform these steps for both a synthetic image of the disk plus companion model, and for our data.
Using the azimuth wedges for both synthetic model and reduced data, we compare the disk plus companion model to the data. We fix all the components of the model except the companion stellocentric radius and position angle. We then perform a best-fit chi squared minimzation over these two freely varying parameters (PA and stellocentric radius) for all epochs of our data (see Table 1). We find a best-fit companion radius of 3.5 AU, and a PA of 90 degrees. We note that as we have limited phase information using this technique, that we are unable to distinguish between PAs separated by 180 degrees. Thus, our best-fit PA is 90 degrees or 270 degrees. Position angle is measured counter-clockwise eastwards from north.
Since our data span two epochs separated by 720 days, we do not fit either epoch of the 11.66 m data perfectly. Further, if a companion indeed explains our data, then this object would undergo orbital evolution between our two epochs. For these reasons, we additionally fit PA separately for each epoch for a fixed radius: 3.5 AU, the best-fit companion radius for all the data. This exercise yields a best-fit PA of 90 degrees for our 2007 data (874 and 1166a), but the PA is not strongly constrained in the 2009 data (1166b and 1830). This is unsurprising, as the 2009 epoch does not contain short-wavelength observations in which the flux from the hot companion would be significant.
To more clearly illustrate asymmetry in our multi-epoch data and to ameliorate apparent discrepancies due to fitting our disk model to the combined 2007 and 2009 datasets, we perform a simple “de-trending” of the model visibility amplitudes. That is, we force the synthesized visibilities to match the shape of the visibilities obtained through our data reduction by applying a multiplicative scaling at each baseline. As visibilities greater than one are unphysical, we set detrended model values that exceed one to be equal to one. In Figure 7, we show the expected magnitude of the asymmetry and compare this to the magnitude of the asymmetry in our observations.
In this figure, the points with error bars show the degree of asymmetry revealed in our data reduction, using the wedge comparison technique described above. We choose the best-fit PA for each epoch for one wedge, and the best-fit PA plus 90 degrees for the second wedge. Though not highly statistically significant, this figure shows that the asymmetries are consistent with our expectation from an asymmetric source of emission at short wavelengths: that the asymmetry in the 8.74 m data is greatest, smaller at 11.66 m, and yet smaller at 18.30 m. The solid and dotted lines are the visibilities obtained by generating a synthetic companion model and directly extracting the normalized visibilities via Fourier techniques, described in the preceding paragraphs.
We present the azimuthally averaged visibilities from our disk plus companion model with best-fit values for companion stellocentric radius and PA in Figure 8, for comparison with our traditional disk models (Figures 3, 4, and 6). We present a synthetic image of this model in Figure 9.
We also examined the closure phases of our data. Closure phases—a measure of the complex phase that removes the phase ambiguity produced by the turbulent atmosphere—are used to discover and constrain asymmetric structure in observed objects (e.g., a companion). We found that the closure phases of our data are insensitive to structures on the small scales we are examining. The closure phases are sufficiently noisy to be largely insensitive to asymmetric structure more compact than 100 mas (our putative companion lies at 65 mas).
|Calvet / Uchida||Ratzka||New Disk Model||Disk + Companion Model|
|(16)||Position Angle (2007)||()||-||-||-||90|
|(17)||Position Angle (2009)||()||-||-||-||-|
Note. – The models used to generate SEDs and visibilities presented in this work. Columns separate models that are “Calvet-like”, “Ratzka-like”, or newly presented in this work. Rows describe the different components in the model, and are (a indicates the component is optically thin, otherwise it is optically thick): (1) the scaling factor applied to the Kurucz model, described in Section 3; (2) the dust sublimation radius, i.e., the inner radius of the optically thin disk; (3) the temperature at the dust sublimation radius; (4) the surface density of the inner, optically thin disk, compared to the surface density used by Calvet et al. (2002); (5) the radius of the directly illuminated optically thick wall; (6) the ratio of width to radius for the transition disk wall; (7) temperature of the optically thick wall; (8) radius of the optically thin wall atmosphere; (9) the ratio of width to radius for optically thin wall atmosphere; (10) temperature of the optically thin wall atmosphere; (11) the surface density of the optically thin wall atmosphere, compared to the surface density used by Calvet et al. (2002); (12) radius of the optically thick, hotter-than-equilibrium ring inside the transition disk wall. This component is only used for the model presented in this work; (13) the ratio of width to radius for the optically thick, hotter-than-equilibrium ring; (14) temperature of the optically thick, hotter-than-equilibrium ring; (15) opacities used in each model—“C” for Calvet, “R” for Ratzka. See Section 3.1.1 (16) The best-fit position angle of our companion, for our 2007 data (17) The best-fit position angle of our companion, for our 2009 data
The presence of a companion in the TW Hya system is consistent with theories of transition disk dissipation (e.g., Goldreich & Tremaine, 1982; Bryden et al., 1999; Rice et al., 2003). The fitted stellocentric radius of our putative companion is 3.5 AU, inside the optically thick disk wall at 3.9 AU. The Keplerian period for such a circular orbit would be 2860 days, assuming a mass for TW Hya of (Webb et al., 1999). The time between our two epochs is 720 days, and so during this time the change of the PA of the companion would be 91 degrees. Since we do not have an estimate of PA for our second epoch, we cannot test for this orbital evolution in our data. Though the 2009 data do not constrain a value for PA, the data are consistent with orbital motion of a companion. This is because the companion’s flux contribution at 11.66 m and 18.30 m is less significant; another epoch of 8.74 m observations would enable this exercise.
We have not quoted uncertainties on any of our fitted parameters. Our parameterizations for the reproduction of the Ratzka et al. (2007) model and our new disk model include 7 and 9 freely varying parameters, respectively. With parameter spaces of increasing dimensionality, model uniqueness and degeneracy of parameters becomes increasingly problematic. High dimensional chi-squared surfaces are complex, and may contain many local minima. To obtain our best-fit parameters, we employ grid-based chi-squared minimizations. We obtain a best-fit PA and stellocentric radius of the companion by simply substituting the hot ring in our new disk model ( 3.2) for the emission from an unresolved source. In doing this, we fix all of the many free parameters in our model (radius of optically thick disk rim, temperature of model components, etc.). Obtaining the uncertainties for PA and radius from a chi-squared surface assuming only one or two varying parameters instead of the actual, larger number, would not provide an accurate estimate of the error in these parameters. Since modeling the multi-dimensional chi-squared surface is computationally prohibitive using our methods, we do not report formal confidence intervals on our model parameters.
While the flux from the outer regions of the disk (at 11.66 m and longer wavelengths) is rather sensitive to changes in and (and has been modeled previously, Calvet et al. (2002, e.g.)), the quality of the fit to the data depend less strongly on and , due in part to the variable nature of TW Hya at the shorter wavelengths. We find that while changes of just 10s of Kelvins or tenths of an AU in the disk cavity wall can lead to an unacceptable fit, changes of a similar magnitude in the hot disk or companion cause a much smaller effect. In fact, a temperature for the hot companion as low as 430 K produces a lower-quality but still acceptable fit to our data: a deficit of flux at wavelengths shorter than 6-8 m can be explained by a larger stellar contribution. Though some model parameters are degenerate, we use this range of temperatures (430-500 K) to investigate potential observables predicted by planetary models.
Even if we assume our lowest acceptable value for , the fitted temperature of our companion is hotter than the disk equilibrium temperature at its stellocentric radius (150 K). Given our fitted surface area of emission, and the assumption that this object emits like a blackbody, the luminosity of this object is 3-6 . The age of the TW Hya system is 10 Myr (Webb et al., 1999), which puts an upper limit on the age of the putative companion.
These properties of our proposed companion are consistent with models of planetary formation. The temperature of our companion (of age Myr) is consistent with the predicted temperatures from planetary thermal evolution models presented in Fortney et al. (2008) for a planet beginning from the core accretion formation models of Hubickyj et al. (2005). The luminosity, however, is sufficiently large to require “hot start” models with arbitrary initial conditions; in this case, our companion’s luminosity is consistent a mass of 8-10 (Fortney et al., 2008). Indeed, the work of Marley et al. (2007) suggests that such a large luminosity is only possible with hot-start models, or if the companion is undergoing an accretion shock phase during the core-accretion formation process. Marley et al. (2007) state, however, that luminosities for the hot-start planetary models are highly uncertain and depend strongly on the initial entropy of the evolution tracks; this makes mass determination of very young planets difficult. We also note that the surface flux density of our companion as a function of wavelength is consistent with those of a 500-700 K model presented in Fortney et al. (2008), except for effects due to opacity, which we do not include (see their Figure 3).
The luminosity of our planet is feasible if it is undergoing an accretion shock phase, though the magnitude of an object undergoing this phase is highly uncertain, and this phase is short-lived (Marley et al., 2007). Another possible explanation for the high magnitude of the companion luminosity is that the companion is surrounded by and accreting circum-planetary material; an increase in the surface area of the emitting region would produce a larger luminosity.
Additionally, we have detected significant variability at 11.66 m in our multi-epoch observations. The origin of this variability is unclear. The fitted size (using a ring emitting at a single wavelength) for the more resolved epoch is AU, compared to the less resolved size of AU. At 11.66 m, emission is dominated by the optically thin transition disk cavity, at relatively small stellocentric radii; indeed, much of the emission at this wavelength comes from the hot inner edge of the optically thin disk, at the dust sublimation radius. The variability is thus likely due to changes in the properties of this inner region. Muzerolle et al. (2009) report “remarkable mid-IR variability” in the transition disk LRLL 31, in the same wavelength range as our observed variable sized emission and on timescales as short as one week. Flaherty & Muzerolle (2010) propose models to explain this sort of mid-IR variability, using non-axisymmetric perturbations in the disk (e.g., a warp with variable scale height, or spiral wave). Flaherty et al. (2011) expands on the earlier work of Muzerolle et al. (2009) with extensive observations of LRLL 31 over many epochs in the infrared. They find that the dust destruction radius stays relatively constant, that accretion and IR excess vary over timescales of weeks, and that changes in scale height and/or warping of the inner disk is likely responsible for the observed variation. They rule out accretion and stellar winds as a cause for the disk changes, and conclude that the most likely explanation for their observations is a companion or a dynamic interface between the stellar magnetic field and the disk. In particular, they describe how a companion, orbiting at an inclination relative to the disk, can drag dust from the disk midplane in a periodic fashion, producing variable IR emission. Given that the inner disk contributes most of the flux at 11.66 m in our model, variability in the disk itself—rather than variable flux from a companion—is the most likely explanation, as in Flaherty et al. (2011). As we only have one epoch of 8.74 m data, we cannot determine if the variability we observe persists at shorter wavelengths, though TW Hya has also been observed to be variable in the near-IR and optical as well (Eisner et al., 2010, and references therein). A multi-epoch, multi-wavelength study of TW Hya, as in Flaherty et al. (2011) for LRLL 31, would yield a greater understanding of the physical mechanisms responsible for this variability.
Recently, other companions have been detected or inferred around transition disk objects. Huélamo et al. (2011) presented evidence of a companion around T Chamaeleontis at mag and at a separation of 627 mas using sparse aperture masking and adaptive optics at the VLT. Eisner et al. (2009) present spatially resolved mid-IR observations from T-ReCS of SR 21, at a separation mas. They predict that in the band the companion has of the flux of the central star, and at . Kraus & Ireland (2011) detect a companion around LkCa 15 using non-redundant aperture masking interferometry at the Keck-II telescope. This object has projected separations of 72, 101, and 88 mas at three different epochs, with a contrast in the band of mag. Kraus & Ireland (2011)’s detection of a companion around LkCa 15 shows evidence of asymmetry and variability in stellocentric radius and flux (though they state that some of the asymmetrical smearing is indicative of poor data quality). They suggest that circumplanetary material is expected around the detected companion, as both planet and star are accreting protoplanetary disk mass.
We predict that in the band the companion around TW Hya will have 0.8-3.0% of the flux of the central star (5.3-3.6 mag), and in the band (4.7 m) the companion will have 2-7% of the flux. Our companion would lie at a separation of 65 mas from the central star. The contrast would be even more favorable at longer wavelengths; the hot companion’s contribution to the SED peaks at around 7 m. Due to its proximity then, TW Hya offers the possibility to detect a companion at a smaller stellocentric radius than, but at comparable angular separation to, these recent companion discovery publications. High resolution mid-IR observations (perhaps speckle interferometry or non-redundant aperture interferometry at the Large Binocular Telescope) should verify our predictions.
We note that Evans et al. (2012) have placed a lower limit on the -band contrast of a companion in a separation range of 40-80 mas around TW Hya at 5.28 mag, consistent with but on the faint end of our L-band contrast estimate of 5.3-3.6 mag. There are several reasons why our prediction for companion contrast at band could be too large, however. For example, if some of the emission from the companion is optically thin, the companion flux at band would be smaller (while maintaining the same flux at m), leading to a larger value for the contrast at while preserving the quality of our fit to the SED. Some optically thin emission would be evidence that this emission is partially due to circumplanetary matter, which would be more spatially extended to produce the same flux as a blackbody component.
We present new mid-IR, spatially resolved measurements of the transition disk object TW Hya, taken in a novel observing mode at the Gemini telescope using the T-ReCS instrument. We observed this object using speckle imaging and reduced the data using Fourier image analysis techniques. Our individual exposures were short enough to freeze the atmosphere’s turbulence, allowing for sub seeing-limited observations. In fact, due to high precision calibration of the PSF using our Fourier methods, we probe spatial scales at the diffraction limit of the Gemini telescope. At all our observational wavelengths, we resolve the science target, TW Hya.
We recreate and present simple models of TW Hya’s disk from the literature. We analyze the compatibility of these models with our new data, and show that existing models do not reproduce our very resolved 8.74 m emission. We create a new model that satisfactorily explains all the available data: our new resolved mid-IR measurements, long baseline interferometric measurements at mm (Hughes et al., 2007), the flux density in the mid-IR, and long baseline mid-IR visibility amplitudes presented by Ratzka et al. (2007). This model has a large, relatively empty cavity as shown in previous works (Calvet et al., 2002; Uchida et al., 2004; Hughes et al., 2007), but also includes a self luminous companion at stellocentric radius of 3.5 AU to explain our highly resolved 8.74 m data. We note also that the behavior of our data is similar to that of Ratzka et al. (2007)’s (8.74 m emission more resolved than 11.66 m emission), who obtain mid-IR observations at another epoch (2005) and at longer baselines (45 m).
In summary, we present new observations with 8.74 m emission more resolved than 11.66 m emission at baselines of 6 m. This unexpected result is consistent with other findings: Ratzka et al. (2007) show 8 m emission more resolved than 12 m emission at much longer baselines. If a companion is responsible for this emission, observations should show an asymmetry, most significant at the peak emission wavelengths of the companion. Our data are consistent with this expectation of asymmetry, but we cannot constrain the details of this asymmetry well. We provide estimates for separation and flux ratio of a putative companion, but these are uncertain due to large number of components in our model, and because the companion is not well-resolved. Evidence for a luminous source at a large stellocentric radius is fairly strong, but our constraints on its properties are weaker.
We measure the size of the emission in our observations at each wavelength, and detect temporal variability at the 2.5- level between two epochs at 11.66 m. We speculate that this variability reflects variable accretion through the inner disk, signatures of dynamical perturbations, and / or perhaps a variable brightness of a self-luminous companion.
The authors would like to thank Laird Close, Phil Hinz, George Rieke, and the referee for helpful suggestions and critiques. This work is based on observations obtained at the Gemini Observatory. TJA was supported by the National Science Foundation through a Graduate Research Fellowship. JAE acknowledges support from an Alfred P. Sloan Research Fellowship. Our Gemini Program IDs are GS-2008A-Q-18 and GS-2007A-Q-38.
- Akeson et al. (2011) Akeson, R. L., et al. 2011, The Astrophysical Journal, 728, 96
- Andrews et al. (2011a) Andrews, S. M., Rosenfeld, K. A., Wilner, D. J., & Bremer, M. 2011a, ApJ, 742, L5
- Andrews et al. (2011b) Andrews, S. M., Wilner, D. J., Espaillat, C., Hughes, A. M., Dullemond, C. P., McClure, M. K., Qi, C., & Brown, J. M. 2011b, ApJ, 732, 42
- Boss & Yorke (1993) Boss, A. P., & Yorke, H. W. 1993, The Astrophysical Journal, 411, L99
- Boss & Yorke (1996) —. 1996, The Astrophysical Journal, 469, 366
- Brown et al. (2009) Brown, J. M., Blake, G. A., Qi, C., Dullemond, C. P., Wilner, D. J., & Williams, J. P. 2009, ApJ, 704, 496
- Bryden et al. (1999) Bryden, G., Chen, X., Lin, D. N. C., Nelson, R. P., & Papaloizou, J. C. B. 1999, ApJ, 514, 344
- Calvet et al. (2002) Calvet, N., D’Alessio, P., Hartmann, L., Wilner, D., Walsh, A., & Sitko, M. 2002, ApJ, 568, 1008
- Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368
- Eisner et al. (2006) Eisner, J. A., Chiang, E. I., & Hillenbrand, L. A. 2006, The Astrophysical Journal, 637, L133
- Eisner et al. (2010) Eisner, J. A., Doppmann, G. W., Najita, J. R., McCarthy, D., Kulesa, C., Swift, B. J., & Teske, J. 2010, ApJ, 722, L28
- Eisner et al. (2003) Eisner, J. A., Lane, B. F., Akeson, R. L., Hillenbrand, L. A., & Sargent, A. I. 2003, ApJ, 588, 360
- Eisner et al. (2009) Eisner, J. A., Monnier, J. D., Tuthill, P., & Lacour, S. 2009, ApJ, 698, L169
- Evans et al. (2012) Evans, T. M., et al. 2012, ApJ, 744, 120
- Flaherty & Muzerolle (2010) Flaherty, K. M., & Muzerolle, J. 2010, ApJ, 719, 1733
- Flaherty et al. (2011) Flaherty, K. M., Muzerolle, J., Rieke, G., Gutermuth, R., Balog, Z., Herbst, W., Megeath, S. T., & Kun, M. 2011, ApJ, 732, 83
- Fortney et al. (2008) Fortney, J. J., Marley, M. S., Saumon, D., & Lodders, K. 2008, ApJ, 683, 1104
- Goldreich & Tremaine (1982) Goldreich, P., & Tremaine, S. 1982, ARA&A, 20, 249
- Gorti et al. (2011) Gorti, U., Hollenbach, D., Najita, J., & Pascucci, I. 2011, The Astrophysical Journal, 735, 90
- Hubickyj et al. (2005) Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2005, Icarus, 179, 415
- Huélamo et al. (2011) Huélamo, N., Lacour, S., Tuthill, P., Ireland, M., Kraus, A., & Chauvin, G. 2011, A&A, 528, L7
- Hughes et al. (2007) Hughes, A. M., Wilner, D. J., Calvet, N., D’Alessio, P., Claussen, M. J., & Hogerheijde, M. R. 2007, The Astrophysical Journal, 664, 536
- Isella et al. (2009) Isella, A., Carpenter, J. M., & Sargent, A. I. 2009, ApJ, 701, 260
- Isella et al. (2009) Isella, A., Carpenter, J. M., & Sargent, A. I. 2009, The Astrophysical Journal, 701, 260
- Isella et al. (2010) Isella, A., Carpenter, J. M., & Sargent, A. I. 2010, ApJ, 714, 1746
- Isella et al. (2008) Isella, A., Tatulli, E., Natta, A., & Testi, L. 2008, A&A, 483, L13
- Kraus & Ireland (2011) Kraus, A. L., & Ireland, M. J. 2011, ArXiv e-prints
- Krist et al. (2000) Krist, J. E., Stapelfeldt, K. R., Ménard, F., Padgett, D. L., & Burrows, C. J. 2000, ApJ, 538, 793
- Mamajek (2005) Mamajek, E. E. 2005, The Astrophysical Journal, 634, 1385
- Marley et al. (2007) Marley, M. S., Fortney, J. J., Hubickyj, O., Bodenheimer, P., & Lissauer, J. J. 2007, ApJ, 655, 541
- Mayor & Queloz (1995) Mayor, M., & Queloz, D. 1995, Nature, 378, 355
- Muzerolle et al. (2009) Muzerolle, J., et al. 2009, ApJ, 704, L15
- Pontoppidan et al. (2008) Pontoppidan, K. M., Blake, G. A., van Dishoeck, E. F., Smette, A., Ireland, M. J., & Brown, J. 2008, ApJ, 684, 1323
- Qi et al. (2004) Qi, C., et al. 2004, The Astrophysical Journal, 616, L11
- Ratzka et al. (2007) Ratzka, T., Leinert, C., Henning, T., Bouwman, J., Dullemond, C. P., & Jaffe, W. 2007, Astronomy and Astrophysics, 471, 173
- Rice et al. (2003) Rice, W. K. M., Wood, K., Armitage, P. J., Whitney, B. A., & Bjorkman, J. E. 2003, MNRAS, 342, 79
- Sitko et al. (2000) Sitko, M. L., Lynch, D. K., & Russell, R. W. 2000, The Astronomical Journal, 120, 2609
- Strom et al. (1989) Strom, K., Strom, S., Edwards, S., Cabrit, S., & Skrutskie, M. 1989, The Astronomical Journal, 97, 1451
- Thalmann et al. (2010) Thalmann, C., et al. 2010, ApJ, 718, L87
- Uchida et al. (2004) Uchida, K. I., et al. 2004, The Astrophysical Journal Supplement Series, 154, 439
- van Leeuwen (2007) van Leeuwen, F. 2007, Astronomy and Astrophysics, 474, 653
- Webb et al. (1999) Webb, R. A., Zuckerman, B., Platais, I., Patience, J., White, R. J., Schwartz, M. J., & McCarthy, C. 1999, The Astrophysical Journal, 512, L63
- Wilner et al. (2003) Wilner, D. J., Bourke, T. L., Wright, C. M., Jørgensen, J. K., van Dishoeck, E. F., & Wong, T. 2003, The Astrophysical Journal, 596, 597
- Wilner et al. (2000) Wilner, D. J., Ho, P. T. P., Kastner, J. H., & Rodríguez, L. F. 2000, The Astrophysical Journal, 534, L101