PHOTO-REVERBERATION MAPPING OF A PROTOPLANETARY ACCRETION DISK AROUND A T TAURI STAR

The Spitzer observations were made in its warm mission with the Infrared Array Camera (IRAC; Fazio et al. 2004) under program 60109 for three periods of staring mode observations at 4.5 μm, on 2010 April 20, 22, and 24, respectively. Each session lasted about 8 hr, with frame time of 2 s and no dithering. We chose 4.5 μm over 3.6 μm with IRAC to maximize the fractional flux contribution from the disk.

Near-simultaneous ground-based observations were made with the FLAMINGOS on the 4 m Mayall telescope at Kitt Peak National Observatory, Arizona; the Spartan Infrared Camera on the 4 m SOAR telescope and the ANDICAM on the SMARTS 1.3 m telescope, both at Cerro Tololo Inter-American Observatory, Chile; and Camila on the Harold L. Johnson 1.5 m telescope of Observatorio Astronómico Nacional (OAN) at Sierra San Pedro Mártir, Mexico. We used the larger telescopes, Mayall/FLAMINGOS and SOAR/Spartan, in H band and the smaller telescopes, SMARTS 1.3 m/ANDICAM and OAN Johnson 1.5 m/Camila, in K band. The time coverage for each telescope is illustrated in Figure 1. The field of view of each participating telescope, including Spitzer, is illustrated in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A total of 27 sources that are in the 2MASS catalog (Skrutskie et al. 2006) were detected in the Spitzer/IRAC field of view. All these sources were covered by ground-based observations on at least one night. For all observations, we relied on the original time in the telescope systems as recorded in the image headers. To correct the light-travel time among participating telescopes at different locations, including Spitzer, the midpoint of each exposure is converted to the barycentric modified Julian date in dynamical time (BMJD_TDB) based on the celestial coordinates of the field.

2.2.1. Photometry

2.2.2. Corrections for IRAC Intrapixel Sensitivity Variations

Staring mode IRAC observations are known to be affected by intrapixel variations in detector response. Several approaches have been introduced to correct these effects, e.g., nonparametric sensitivity mapping (Ballard et al. 2010), polynomial fitting of the stellar centroid (Reach et al. 2005; Knutson et al. 2008), and power spectrum minimization (Cody & Hillenbrand 2011); see Spitzer Science Center (2015). None of the methods are applicable to saturated data, while YLW 16B was saturated on the BCD images on all three nights of our observations. To deal with the pixel phase effect and the saturation, we made photometry with two different apertures on the BCD images: a small aperture with a radius of 2.8 pixels, and a larger aperture of 5.6 pixels. Considering that an overexposed point-spread function (PSF) should first saturate the pixels around the centroid, we subtracted the small aperture flux from the large counterpart to obtain the brightness of the annular PSF wing, which is less prone to saturation and is fainter than the entire stellar flux by a constant factor only dependent on the instrumental PSF. Though with reduced S/N, an important benefit of this method is that the PSF wing contains more pixels with more even illumination and thus is much less affected by the intrapixel variations, so we needed to make no corrections for this effect.

On the first night, the opto-center of the YSO coincidentally lay close to and wandered around the corner of four adjacent pixels. We found the photometry values offset and noisier when the PSF centroid was on the northwest pixel, and we chose to reject the corresponding part of the data. On the second night, the photometry of the PSF wing worked well for all the data. Unfortunately, on the third night our approaches for detrending the pixel phase effect were not successful, likely owing to stronger saturation and the particular pixel location where the star fell on this night; we had to discard the 4.5 μm data. The final light curves, after corrections for the IRAC pixel phase effect were applied, are shown in Figure 3.

Zoom In Zoom Out Reset image size

As a general strategy to maximize the time resolution, the experiment design requested all participating telescopes, including Spitzer, to take new exposures immediately after the previous image was read out. The resultant mean cadence was 13 s for OAN Johnson 1.5 m/Camila, 17 s for SOAR/Spartan, 20 s for Mayall/FLAMINGOS and SMARTS 1.3 m/ANDICAM, and 3.4 s for Spitzer/IRAC. None of the instruments yielded evenly spaced time series, but the variations of the sampling frequency of each instrument were small (mostly less than 1 s with occasional omissions of a data point) and appeared to be random. These should have no effect on our science results.

To preserve all information in the original data, before computing each cross-correlation function (CCF), we used the faster of the average cadences in the two compared time series as the equidistant time step to interpolate both data sets, except for skipping the periods that one or both time series did not cover. Then, we computed the CCFs between the time series in H and K bands and at 4.5 μm. Each CCF within ±170 s (time for 100 Spitzer BCD frames) was fitted by a Gaussian to look for any time lag and by a quadratic polynomial for a sanity check. The time lags estimated this way are expected to be more robust (Zhang & Wu 2006) and do not necessarily coincide with the exact location of a local minimum, which is more prone to noise influence.

For CCFs involving 4.5 μm data, we achieved the most accurate time lag measurements on the second night. Because of the lower variation amplitudes and shorter overlaps of time coverage, on the first night the peaks of the CCFs were shallower and worse defined, resulting in less accurate time lag measurements. To validate the reality of the lag, we performed the following test. At each wavelength, we first fitted the light curve with a cubic spline function and subtracted it from the data, giving two signal streams: (1) the fitted variations without noise, and (2) the time-series photometric noise detrended for the source variations. We then added these two signals together with different artificial phase lags and extracted the CCFs from the artificial combined signal. The results suggested that the detected time lags between each pair of wavelengths have a 1σ uncertainty of ∼3.9 s and are free from noise aliasing.

We found that the light curves of YLW 16B in H and K bands were consistently synchronized over all three nights with an insignificant average offset of −1.4 ± 2.4 s; the 4.5 μm signals lagged behind both H and K by 74.5 ± 3.2 s, accurately measured on the second night and consistent within 2σ errors with the less constraining data from the first night (Figure 4).

Zoom In Zoom Out Reset image size

To verify the detection of the central star and constrain its spectral type, we utilized the high-resolution spectrum taken with Keck/NIRSPEC on 2000 May 30 (R ∼ 18,000; Doppmann et al. 2005) to look for stellar absorption features in H and K bands. Although Doppmann et al. (2005) did not identify photospheric lines in the NIRSPEC spectra of YLW 16B, after careful inspection of the data we felt that several such lines might exist. To test their reality, the spectrum was cross-correlated with a BT-Settl synthetic model spectrum (Allard et al. 2012) from 1800 to 6000 K in increments of 200 K with solar metallicity and gravity $\mathrm{log}g\;=\;4.5$. Both observed and synthetic spectra were normalized. The NIRSPEC spectrum had higher spectral resolution than the BT-Settl template, so we boxcar-smoothed it over 3 pixels, which made the data less noisy and with comparable resolution to the model. We then identified the stellar temperature and wavelength shift that yielded the maximum cross-correlation. As shown in Figure 5, the largest cross-correlation appears to be at 4200 K in the spectral order around 2.1 μm. However, a closer inspection found inconsistency with other wavebands in both temperature and shift. Combining all three orders, the spectrum suggests the most likely effective temperature of 3600 ± 600 K, corresponding to a spectral type of ∼M2. Adopting a uniform spectral calibration provides solid identifications of the Na i absorption triplet at 2.21 μm, and the Mg i and Al i absorption lines at 2.09 μm at less significance.

Zoom In Zoom Out Reset image size

We experimented with line depths scaled to values from 5% to 25% to accommodate a corresponding range of star-to-disk flux ratios and possible emission filling-in of the absorption lines. However, we found that the cross-correlation produced no useful constraint to the line depth scale.

As a sanity check, we compared the findings with the lower-resolution spectra in H and K bands, observed by IRTF/SpeX weeks before and after the reverberation experiment on 2010 March 27, April 1, April 2, and May 7 (Faesi et al. 2012). Before cross-correlation, the spectra observed at all four epochs were combined for higher S/N. The synthetic BT-Settl templates were smoothed to the SpeX spectral resolution (R ∼ 1200). Both data sets were normalized. The cross-correlation, as shown in Figure 5, did not yield any obvious peak. However, the maximum also occurred at T_eff = 3600 K with consistent wavelength calibrations for both H and K bands, though at a weak significance (cross-correlation $\rho \sim 0.25$).

In addition, the combined SpeX spectrum also reveals several atomic emission lines that are typical for evaporating disks with disk winds (Figure 6). The most prominent are Br γ of H i at 2.17 μm, the Mg i lines at 1.49 and 1.50 μm, and the Na i lines at 2.21 μm. Interestingly, the 2.21 μm Na i lines were absorption features in the NIRSPEC spectrum, but intermittently appeared as emission lines with variable equivalent widths in the SpeX spectra. This suggests that emission filling of absorption lines plays a critical role in the detectability of some photospheric absorption features. Therefore, the line depth scale may not be a valid constraint to the fractional flux contribution in this case. The most reliable flux ratios for the star and the disk may have to be estimated by other means.

Zoom In Zoom Out Reset image size

We used a color–magnitude diagram to identify the mechanisms underlying the observed variations. Because the primary energy source of the system is the central protostar and the 4.5 μm photometry is most contaminated by the radiation of the disk, we focused on the H and K data, leaving us K versus H − K as the only option for a color–magnitude diagram. Because K and H − K are not independent, errors in K will affect the H − K color, producing a false negative slope. To reduce this effect to an insignificant level while retaining some time information, we averaged the H and K photometry every 204 s. The range is defined as 60 times the fastest average cadence in the experiment (3.4 s). Such averaged photometry was plotted in the color–magnitude diagram in Figure 7. The nominal errors were 0.002 mag in both K and H − K.

Zoom In Zoom Out Reset image size

In H and K bands, an important reference is the 2MASS calibration field photometry of YLW 16B (Parks et al. 2014). Over a time span of ∼900 days, we found that the YSO varied between 9.25 and 11.57 mag in K_s and between 2.55 and 3.20 in H − K_s. Tracking the variations in the color–magnitude diagram of K versus H − K, we found mixed behaviors of YLW 16B that can be classified into three types. The long-term variations of the YSO form a locus of positive slopes with moderate scatter over the time span of the 2MASS calibration field observations; the relative color–magnitude relation between the first two nights of our reverberation observations (for which we have comparable H and K photometry and good 4.5 μm data) had a negative slope, also with moderate scatter; finally, the hourly light curve for each night appears to have mixed behavior with no dominant trend.

Such complexity indicates that different variability mechanisms are at work on different timescales. We considered a number of possible variability mechanisms: (1) Variable extinction. Both interstellar and circumstellar extinction may vary with time (Carpenter et al. 2001; Cody et al. 2014). The interstellar extinction law has been well measured toward the Galactic plane and several star-forming regions (Indebetouw et al. 2005; Flaherty et al. 2007; Chapman et al. 2009). Thanks to ongoing dust growth, circumstellar dust around YSOs may have a different size distribution and optical properties from its interstellar counterpart (McClure et al. 2013), possibly leading to a different extinction law. However, because the Keplerian timescale of the inner disk is orders of magnitude shorter than the time needed for global disk evolution, circumstellar extinction is expected to vary quasi-periodically (Hillenbrand et al. 2013; Stauffer et al. 2015). (2) Starspots on the stellar photosphere (Guenther & Hessman 1993; Nguyen et al. 2009; Pozo Nuñez et al. 2015). Magnetic activity within the YSO can trigger cool starspots like sunspots, and the accretion mass flow from the disk through the magnetic field lines may produce accretion "hot spots" (Stauffer et al. 2014) near the magnetic poles on the photosphere of the protostar (Hartmann et al. 1994). Variability associated with starspots should be on the timescale of the rotation period of the protostar. Radiation from starspots is usually considered as blackbody, but with different temperatures from the rest of the photosphere. (3) The accreting gas along the magnetic field lines just above the photosphere of the protostar. Such material is supposed to be the precursor of accretion hot spots. Differences in spectral properties are that the accreting gas can be optically thin on its way to the star (free–free emission) and does not have to occupy or obscure any area of the stellar photosphere because of the edge-on line of sight, even if it is optically thick (free-bound or blackbody-like emission). (4) Variations of the observable thermally emitting area of the inner disk wall. The temperature of the inner disk wall is likely close to the dust sublimation temperature (Monnier & Millan-Gabet 2002; Muzerolle et al. 2003) and radiates like a blackbody. The timescale for the inner disk variability should be quasi-Keplerian.

We plotted and compared all the observations and theoretical models in the color–magnitude diagram in Figure 7(A). Each data point from the reverberation observations, averaged over 204 s, is individually plotted. The errors are smaller than the symbol sizes and are not drawn. Although these observations were not calibrated to the 2MASS standard, they are self-consistent and are valid for relative comparisons. The interstellar extinction vector is labeled with the 2MASS extinction law (Flaherty et al. 2007). The 2MASS data are divided into five bins of K_s magnitude, from 9.0 to 11.5 with a bin size of 0.5. The median H − K_s color of each bin is plotted, with the first and third quartiles indicated by horizontal bars.

The exact placement of the spotless protostar on the color–magnitude diagram is poorly constrained and subject to the uncertainties in extinction. In the diagrams, we place it arbitrarily for easier comparisons. The stellar photosphere is considered as a 3600 K blackbody. The solid lines of different colors represent the models of the variability mechanisms considered, as labeled in the legend. The cool and hot starspots are assumed to be 2000 and 8000 K, respectively. Three temperatures for the optically thick (blackbody-like) accreting gas are plotted, 6000, 10,000, and 14,000 K. For the cool and hot starspots and the optically thick accreting gas models, each line is drawn up to the point when the emitting areas of the added components are 40% of the area of a hemisphere of the stellar photosphere, no matter whether the photosphere is occupied (starspots) or not (accreting gas). The model line for variable disk-emitting area is plotted to the point at which the disk-emitting area is 12 times larger than a hemisphere of the stellar photosphere.

From Figure 7, it is clear that the overall positive slope (note the inverted Y-axis) of the long-term variations in the 2MASS data primarily reflects the variations of activities related to the protostar, like starspots or optically thick accreting gas, since those are the only models that have potentially matching positive slopes. The scatter in 2MASS H − K_s is similar to the amplitude of the short-term variations in the reverberation experiment, implying that changes similar to those seen in our monitoring have occurred frequently in the past.

The negative slope, as seen between the first two days of the reverberation experiment, appears to be caused by variation in the observable disk-emitting area. To verify this hypothesis, we referred to the time-resolved information of the IRTF/SpeX spectra at the four quasi-contemporary epochs. In the flux-calibrated data, shown in Figures 8(A) and (B), we noticed that the spectra on March 27 and May 7 had practically identical fluxes and spectral slopes in H band, with synthetic photometry agreeing to within 2%, probably better than the uncertainties. But the March 27 data were significantly brighter and redder in K, providing the basis to compare two different states of the inner disk. We subtracted the K-band spectrum on May 7 from that on March 27 and dereddened the difference with A_V = 30. We found that the differential spectrum, as shown in Figure 8(C), does not agree with optically thin free–free emission, but can be reasonably fitted with a blackbody. (The best-fit temperature is 1100 K, but this is only a very rough estimate because of the calibration uncertainties and limited wavelength coverage from which it is derived.) This reveals the spectrum "added" to K band on March 27 that made it brighter and redder than on May 7, and it confirms that a variable disk-emitting area is likely responsible for the negative slope in Figure 7.

Zoom In Zoom Out Reset image size

After associating the long-term and daily variations with their respective model components, now we focus on the hourly variations that led to the detected time lag. A zoom-in view of the color–magnitude diagram of the reverberation data is shown in Figure 7(B). Despite a significant day-to-day shift, the intraday observations appear to have mixed behavior with no dominant color trend. The squared linear correlation coefficients are ∼0.001 for the data of both nights.

The evolution between consecutive data points (separated by 204 s) is reversible back and forth in the diagram. The disk-emitting area could not possibly increase and then decrease on such short timescales, making the model incompatible with the hourly changes. Variable extinction may play a role over any timescale. But the fixed extinction vector alone cannot explain the random directions in the color–magnitude evolution. The most likely interpretation for the hourly variations is that they are caused by a changing continuum from the accreting gas (Guenther & Hessman 1993; Stauffer et al. 2014). The net effect on the output of the source is a varying combination of optically thick and optically thin emission, or that the optical depth was varying in an intermediate range. Variable extinction may still play a role, but even if it did, it does not introduce any time lag between different wavelengths. Therefore, the source of the variable part of the YLW 16B radiation should geometrically be in the magnetically confined accreting columns right above the photosphere of the protostar.

Given the edge-on geometry of YLW 16B, the disk emission we see should mostly come from the far side of the inner rim (Lahuis et al. 2006; Bast et al. 2013). Therefore, to have the direct stellar light curve (from the near side) and responding disk light curve (from the far side) correlated, the source of the variability signals must be seen from both sides. The only satisfactory place for this source is near the limb on the star. Such visibility can be naturally explained if the variability signals are from the accreting gas, because the accretion funnels should hit the star near the magnetic poles, which are likely near or even aligned with the stellar rotation axis. These areas should always be near the limb of the edge-on star as seen from the observer.

Building a spectral energy distribution (SED) model that fits the available measurements of YLW 16B at all times is beyond the scope of this paper. Instead, with this analysis we aim to dissect the YSO to estimate the flux contribution of each emitting entity at our working wavelengths at the epoch of the reverberation experiment. This analysis will help us understand the observations of variations among the source components.

3.3.1. Color Blending

Considering the highly variable nature of the source, to deconvolve the colors of the individual source components, we should only adopt the data taken at the time of the reverberation experiment. However, using our photometry ignores the color terms for the filters and for the source SEDs, which are not possible to be corrected given the strong variability of the source. Therefore, the H- and K-band photometry of the reverberation observations is not precisely on the 2MASS standard scale. The differences in bandpass may introduce systematic errors, because the extinction law and colors of young stars are established based on the 2MASS definitions of JHK_s wavebands. To overcome the problem, we chose to use the 2MASS photometry, but only picked values that are representative of the state of the YLW 16B system at the time of the reverberation experiment.

For this purpose, we compared the averaged reverberation H and K magnitudes to their counterparts in the multiepoch 2MASS calibration field photometry (Parks et al. 2014) by computing their combined magnitude differences, defined as

$$\begin{eqnarray}&&{M}_{D}\;=\;\sqrt{{({H}_{\mathrm{reverb}}-{H}_{2\mathrm{MASS}})}^{2}+{({K}_{\mathrm{reverb}}-{K}_{s,2\mathrm{MASS}})}^{2}}.\end{eqnarray} \tag{ 1 }$$

The average magnitudes of the reverberation observations are H = 12.14, K = 9.29, and $[4.5]\;=\;5.75$. The smallest difference in all 2MASS entries is ${M}_{D,\mathrm{min}}$ = 0.132 mag, with the color of

$$\begin{eqnarray*}{(J-H)}^{\mathrm{obs}} & = & 4.122\pm 0.046\\ {(H-{K}_{s})}^{\mathrm{obs}} & = & 2.663\pm 0.014.\end{eqnarray*}$$

There are 21 sets of 2MASS measurements with ${M}_{D}\leqslant 3{M}_{D,\mathrm{min}}$. Averaging all of them yields essentially the same color indices within errors, which means that the color is insensitive to our selections. In the following calculations we used the single 2MASS measurement with the minimum M_D and assumed that it resembles what the reverberation photometry should be if calibrated to the 2MASS standard. In addition to reconciling the definitions of the wavebands, another important benefit of using the 2MASS data is the inclusion of J-band photometry, which is more sensitive to the extinction and stellar output and thus provides better constraints on the component deconvolution.

The observed color of the YSO is made up of three components, the radiation from the stellar photosphere, from the inner disk, and from the accreting gas, all behind the foreground extinction along the line of sight, yielding

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{F}_{J}}{{F}_{H}}\right)}^{\mathrm{obs}}=\displaystyle \frac{{10}^{-0.4{A}_{J}}({F}_{J}^{\ast }+{F}_{J}^{\mathrm{disk}}+{F}_{J}^{\mathrm{gas}})}{{10}^{-0.4{A}_{H}}({F}_{H}^{\ast }+{F}_{H}^{\mathrm{disk}}+{F}_{H}^{\mathrm{gas}})}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{F}_{H}}{{F}_{{K}_{s}}}\right)}^{\mathrm{obs}}=\displaystyle \frac{{10}^{-0.4{A}_{H}}({F}_{H}^{\ast }+{F}_{H}^{\mathrm{disk}}+{F}_{H}^{\mathrm{gas}})}{{10}^{-0.4{A}_{{K}_{s}}}({F}_{{K}_{s}}^{\ast }+{F}_{{K}_{s}}^{\mathrm{disk}}+{F}_{{K}_{s}}^{\mathrm{gas}})},\end{eqnarray} \tag{ 3 }$$

where F is flux density and A is the extinction in magnitudes. Subscripts and superscripts denote the wavebands and sources of radiation, respectively. Although the adopted set of 2MASS photometry is the closest in brightness to our reverberation observations, to be conservative, in the equations we did not connect it to the 4.5 μm reverberation photometry because of the nonsimultaneity. The averaged 4.5 μm photometry of the reverberation experiment is left for a consistency check of the JHK_s-based models.

We took the colors of young stellar photospheres from observations (Luhman et al. 2010). Given M2 as the approximate spectral type of YLW 16B, we adopted the intrinsic color of the photosphere as

$$\begin{eqnarray*}{(J-H)}^{\ast } & = & 0.70\\ {(H-{K}_{s})}^{\ast } & = & 0.20\\ {({K}_{s}-[4.5])}^{\ast } & = & 0.17.\end{eqnarray*}$$

The radiation of the inner disk should be close to a blackbody (Monnier & Millan-Gabet 2002; Muzerolle et al. 2003). Assuming a temperature of 1100 K (see Section 3.2 for the temperature estimate), the disk color should be

$$\begin{eqnarray*}{(J-H)}^{\mathrm{disk}} & = & 2.47\\ {(H-{K}_{s})}^{\mathrm{disk}} & = & 1.58\\ {({K}_{s}-[4.5])}^{\mathrm{disk}} & = & 2.51.\end{eqnarray*}$$

We find that the extinction law specifically determined in ρ Ophiuchi is not accurate enough for our purpose (Elias 1978; Kenyon et al. 1998). Fortunately, the infrared extinction law of ρ Ophiuchi is practically identical to that of other star-forming regions (Harris et al. 1978; Flaherty et al. 2007). Here we adopted the generic extinction law for 2MASS and Spitzer wavebands (Flaherty et al. 2007; Chapman et al. 2009), where

$$\begin{eqnarray}&&{A}_{J}=2.5{A}_{{K}_{s}}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{A}_{H}=1.55{A}_{{K}_{s}}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{A}_{[4.5]}=0.53{A}_{{K}_{s}}.\end{eqnarray} \tag{ 6 }$$

We assume S as the disk-emitting area in units of the area of the observable surface of the protostar. Then, we can construct all fluxes by introducing proper ratios of the fluxes between different wavebands, sources, and temperatures. If a flux ratio is ${\eta }_{{V}_{1}/{V}_{2}}^{C}$, where V₁ and V₂ are the varied conditions (e.g., wavelengths or temperatures) and C is the common condition, the unreddened flux at the ith waveband in a set of n band photometry can be written as

$$\begin{eqnarray}&&{F}_{{\lambda }_{i}}^{\ast }={E}_{{\lambda }_{n}}^{\mathrm{disk}}{\eta }_{{T}^{\ast }/{T}^{\mathrm{disk}}}^{{\lambda }_{n}}{\eta }_{{\lambda }_{i}/{\lambda }_{n}}^{\ast }\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{F}_{{\lambda }_{i}}^{\mathrm{disk}}={E}_{{\lambda }_{n}}^{\mathrm{disk}}S{\eta }_{{\lambda }_{i}/{\lambda }_{n}}^{\mathrm{disk}}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{F}_{{\lambda }_{i}}^{\mathrm{gas}}={F}_{{\lambda }_{n}}^{\ast }{\eta }_{\mathrm{gas}/\ast }^{{\lambda }_{n}}{\eta }_{{\lambda }_{i}/{\lambda }_{n}}^{\mathrm{gas}}\end{eqnarray} \tag{ 9 }$$

where ${E}_{{\lambda }_{n}}^{\mathrm{disk}}$ is the disk flux emitted per unit area at the longest (nth) wavelength. If the accreting gas is optically thin, its free–free emission should have a spectrum with

$$\begin{eqnarray}&&{\eta }_{{\lambda }_{i}/{\lambda }_{n}}^{\mathrm{gas}}\;=\;{\left(\displaystyle \frac{{\lambda }_{n}}{{\lambda }_{i}}\right)}^{-0.1}.\end{eqnarray} \tag{ 10 }$$

Then, Equations (2) and (3) become

$$\begin{eqnarray}{\left(\displaystyle \frac{{F}_{J}}{{F}_{H}}\right)}^{\mathrm{obs}}= & \displaystyle \frac{{10}^{-{A}_{{K}_{s}}}\left[{\eta }_{J/{K}_{s}}^{\ast }{\eta }_{3600\;{\rm{K}}/1100\;{\rm{K}}}^{{K}_{s}}+{\eta }_{J/{K}_{s}}^{\mathrm{disk}}S+{\eta }_{3600\;{\rm{K}}/1100\;{\rm{K}}}^{{K}_{s}}{\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}{\left(\displaystyle \frac{{\lambda }_{{K}_{s}}}{{\lambda }_{J}}\right)}^{-0.1}\right]}{{10}^{-0.62{A}_{{K}_{s}}}\left[{\eta }_{H/{K}_{s}}^{\ast }{\eta }_{3600\;{\rm{K}}/1100\;{\rm{K}}}^{{K}_{s}}+{\eta }_{H/{K}_{s}}^{\mathrm{disk}}S+{\eta }_{3600\;{\rm{K}}/1100\;{\rm{K}}}^{{K}_{s}}{\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}{\left(\displaystyle \frac{{\lambda }_{{K}_{s}}}{{\lambda }_{H}}\right)}^{-0.1}\right]}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}{\left(\displaystyle \frac{{F}_{H}}{{F}_{{K}_{s}}}\right)}^{\mathrm{obs}}= & \displaystyle \frac{{10}^{-0.62{A}_{{K}_{s}}}\left[{\eta }_{H/{K}_{s}}^{\ast }{\eta }_{3600\;{\rm{K}}/1100\;{\rm{K}}}^{{K}_{s}}+{\eta }_{H/{K}_{s}}^{\mathrm{disk}}S+{\eta }_{3600\;{\rm{K}}/1100\;{\rm{K}}}^{{K}_{s}}{\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}{\left(\displaystyle \frac{{\lambda }_{{K}_{s}}}{{\lambda }_{H}}\right)}^{-0.1}\right]}{{10}^{-0.4{A}_{{K}_{s}}}({\eta }_{3600\;{\rm{K}}/1100\;{\rm{K}}}^{{K}_{s}}+S+{\eta }_{3600\;{\rm{K}}/1100\;{\rm{K}}}^{{K}_{s}}{\eta }_{\mathrm{gas}/\ast }^{{K}_{s}})}.\end{eqnarray} \tag{ 12 }$$

In Equations (11) and (12), the flux ratios on the left are directly known from the observed color of the YSO; all η values on the right-hand side, except for ${\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}$, are known either from Planck's law or from the observation-based intrinsic color of the young stellar photosphere (Luhman et al. 2010). Therefore, the equations are model independent. Now we have two equations and three free parameters, ${A}_{{K}_{s}}$, S, and ${\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}$, the K_s-band scale factor of the emission of the accreting gas compared to the stellar photospheric flux. For any given value of ${\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}$, we can numerically solve for ${A}_{{K}_{s}}$ and S.

Taking the 2MASS photometry with the minimum M_D from the reverberation observations, we found that ${A}_{{K}_{s}}$ and S are not sensitive to ${\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}$ within a generous range. Although the accurate value of ${\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}$ at the time of the reverberation experiment is unknown, a lower limit is constrained by the hourly peak-to-peak amplitude of 20% observed in both H and K bands. At 4.5 μm, thanks to the dominance of the disk emission, the same variable component may only cause a peak-to-peak amplitude ≲10%, which would have been submerged in the noise (rms ∼ 8%) of the Spitzer/IRAC data. That being said, the amplitude of the disk response is much greater than that of the original signal. This would be naturally explained if the triggering signals are causing temperature (or chemical) changes in the disk (see Flaherty et al. 2014).

Accordingly, for ${\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}\;=\;0.25$, which would mean that the accretion excess can just account for the observed amplitude in K band, we obtained ${A}_{{K}_{s}}\;=\;3.32$ and S = 142; for ${\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}\;=\;1$, which would mean that the accretion excess is as strong as that of the stellar photospheric flux in K band, ${A}_{{K}_{s}}\;=\;3.46$ and S = 141. This degeneracy is because the flux density of the stellar photosphere of a young M2 star peaks in H band (Luhman et al. 2010), causing a relatively flat SED over the JHK_s wavebands, similar to that of the optically thin free–free emission assumed for the accreting gas. Extrapolating the color-blending equations to 4.5 μm, for ${\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}$ values between 0.25 and 1, the predicted magnitude is in the range of 5.76–5.78, in good agreement with the averaged reverberation photometry $[4.5]\;=\;5.75$. Larger values of ${\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}$ would mean that the accretion excess dominates over the stellar photospheric emission in K band, and would predict fainter [4.5] magnitude, which are not supported by the observations.

3.3.2. Implications for Disk Properties

In summary, the variable signals in H and K bands should be mostly from the emission of the accreting gas near the stellar photosphere, whereas the signals observed at 4.5 μm should be dominated by the disk response. This explains the observed simultaneity between the time series in H and K, as well as the lag at 4.5 μm.

Analysis of molecular spectral features suggests that YLW 16B is close to edge-on (Lahuis et al. 2006; Bast et al. 2013), in which case the observed disk continuum is mostly from the far side of the inner disk wall (Lahuis et al. 2006). Such a viewing angle requires the triggering signals from the protostar to travel an additional distance before reaching us. To first-order approximation, assuming that the inner disk radius is azimuthally uniform, the average additional travel distance of the lagged signals should be

$$\begin{eqnarray}&&\displaystyle \frac{{\displaystyle \int }_{0}^{\pi }{R}_{\mathrm{in}}{\mathrm{sin}}^{2}\theta d\theta }{{\displaystyle \int }_{0}^{\pi }\mathrm{sin}\theta d\theta }\;=\;\displaystyle \frac{\pi }{4}{R}_{\mathrm{in}}\end{eqnarray} \tag{ 13 }$$

where θ is the azimuthal angle and R_in is the radius of the inner disk rim, which, given the observed time lag between H/K and 4.5 μm in YLW 16B, should be 0.084 ± 0.004 au. The nominal uncertainty of 0.004 au only accounts for the propagated error from CCF centering and noise evalution and should be taken as a lower limit. In converting the time lag to the inner disk radius, we had to adopt several assumptions and simplifications, such as an isotropic radiation field of the protostar, azimuthal uniformity of the disk, negligible disk rim curvature, simplified geometry, etc., which may be significant error sources beyond the observational noise. However, these errors are very difficult, if at all possible, to quantitatively estimate with the available data. The total error of the radius of the inner disk rim is likely higher by a factor of several, on the order of 0.01 au.

The photo-reverberation measurement provides an independent verification of the current paradigm of the geometry and physics of the inner disk, which has been largely based on model-dependent interferometric results. Therefore, we will compare it with the expected radii from various truncation mechanisms.

To carry out this test, first we estimated the properties of the protostar and the accretion disk of YLW 16B at the time of the reverberation experiment. The total luminosity of a YSO is composed of the intrinsic stellar luminosity and the accretion luminosity, i.e., ${L}_{\mathrm{tot}}\;=\;{L}_{\ast }+{L}_{\mathrm{acc}}$. The Br γ emission is a well-established indicator of accretion (Muzerolle et al. 1998). We found that the equivalent width of this line is highly variable in the quasi-contemporary SpeX spectra (Figure 6(A)), ranging from −2.0 Å on March 27 to −4.8 Å on April 2. After converting the range to Brγ luminosity (${L}_{\mathrm{Br}\gamma }$) and compensating for extinction of ${A}_{{K}_{s}}\;=\;3.3$, we used the relation (Natta et al. 2006)

$$\begin{eqnarray}&&\mathrm{log}({L}_{\mathrm{acc}}/{L}_{\odot })\;=\;0.9[\mathrm{log}({L}_{\mathrm{Br}\gamma }/{L}_{\odot })+4]-0.7\end{eqnarray} \tag{ 14 }$$

to derive an accretion luminosity of YLW 16B between 0.31 and 0.63 ${L}_{\odot }$.

On the other hand, the stellar luminosity of YLW 16B is highly uncertain. The isochrones between stellar models are significantly divergent toward younger ages and lower masses (Baraffe et al. 2002; Hillenbrand et al. 2008). Together with the large error in the effective temperature, we could only estimate a vast range from 0.1 to 2 ${L}_{\odot }$. To obtain better constraints, we took an alternative approach by integrating the 2MASS photometry that is best matched with the reverberation averages with the minimum M_D. We found that the combined 2MASS JHK_s luminosity of the YSO is 0.56 ${L}_{\odot }$. Referring to the color-blending equations in Section 3.2, the viable range of ${\eta }_{\mathrm{gas}/\ast }^{{K}_{s}}$ from 0.25 to 1 means that the protostar, if integrated over the JHK_s wavebands, should contribute between 38% and 51% of the total flux. Meanwhile, we extrapolated the observed colors of young stellar photospheres (Luhman et al. 2010) with the Wien approximation (toward shorter wavelengths) and Rayleigh–Jeans law (toward longer wavelengths). For K6–M6 stars that correspond to the range of possible temperature of YLW 16B, we found that the combined JHK_s flux makes 21%–22% of the entire stellar luminosity (integrated from 0.1 to 5000 μm) and is insensitive to the spectral type. Therefore, the stellar luminosity of YLW 16B can be narrowed down to 0.96–1.35 ${L}_{\odot }$. Added to the range of accretion luminosity, this makes the total luminosity of the YSO between 1.27 and 1.98 ${L}_{\odot }$. Finally, we adopted ${M}_{\ast }\;=\;0.7$ ${M}_{\odot }$ for the stellar mass based on BT-Settl models, with the uncertainties allowing values from 0.5 to 1 ${M}_{\odot }$.

The accretion rate can be obtained by

$$\begin{eqnarray}&&\dot{M}\;=\;\displaystyle \frac{{L}_{\mathrm{acc}}{R}_{\ast }}{{{GM}}_{\ast }(1-{R}_{\ast }/{R}_{\mathrm{in}})}\end{eqnarray} \tag{ 15 }$$

where R_in = 0.084 au. For the range of L_acc from the quasi-contemporary SpeX observations, the accretion rate of YLW 16B varies between $5.4\times {10}^{-8}$ and $1.1\times {10}^{-7}$ ${M}_{\odot }$ yr⁻¹ around the time of the reverberation experiment, with an average of $\sim 8\times {10}^{-8}$ ${M}_{\odot }$ yr⁻¹.

The magnetic field of the protostar plays a critical role in the accretion process. Though the detailed topology of the magnetic field may be complex (Johns-Krull & Cauley 2014), an important and model-independent metric is the corotation radius, R_cor, at which the Keplerian orbital period of circumstellar material matches the rotation period of the star. The corotation radius is at the balance of star–disk interaction; for accretion onto the star to proceed, the magnetospheric truncation radius has to stay interior to R_cor (Bouvier et al. 2007). The corotation radius is given by

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{\mathrm{cor}}}{1\;\mathrm{au}}\;=\;1.957\times {10}^{-2}{\left(\displaystyle \frac{{M}_{\ast }}{{M}_{\odot }}\right)}^{1/3}{\left(\displaystyle \frac{{P}_{\ast }}{1\mathrm{day}}\right)}^{2/3}\end{eqnarray} \tag{ 16 }$$

where P_* is the stellar rotation period, which is unknown for YLW 16B. Nevertheless, given the color of $[3.6]-[8.0]\;=\;2.42$ (Gutermuth et al. 2009) and the correlation between rotation periods and IRAC excessess (Rebull et al. 2006), one stellar rotation of YLW 16B likely lasts between 6 and 10 days. Thus, the corotation radius should be in the range from 0.06 to 0.08 au. This range, if considered as the upper limit for the magnetospheric truncation radii, is marginally compatible with our photo-reverberation measurement. Therefore, magnetospheric truncations, which should occur at even smaller radii, are unlikely to be the determining mechanism for the size of the inner disk rim of YLW 16B.

Since the dust sublimation radius depends on the total luminosity of the central star, different disk models can be compared with the observations in a size–luminosity diagram (Millan-Gabet et al. 2007), as shown in Figure 9. Based on interferometric measurements, the size–luminosity diagram has been used to suggest that some accretion disks around high-mass Herbig Be stars are "classical,"¹⁴ i.e., optically thick but geometrically thin. In comparison, disks around intermediate-mass Herbig Ae stars can be well described by a directly heated and back-warmed "puffed-up" inner rim at 1500–2000 K (Eisner et al. 2004), resulting in a weak trend toward relatively increasing apparent sublimination temperatures. However, the trend does not necessarily extrapolate to the regime of solar-mass T Tauri stars, many of which appear to have larger inner disk cavities. Oversized inner disk holes around T Tauri stars have raised questions about the roles of possibly unrecognized physical processes, such as lower dust sublimation temperatures, magnetospheric pressure (Eisner et al. 2007), or viscosity heating (Muzerolle et al. 2004).

Zoom In Zoom Out Reset image size

As an independent confirmation, the reverberation inner radius of YLW 16B is larger than the model prediction for a geometrically thin "classical" disk, but is consistent with a "puffed-up" inner disk rim truncated at ∼1500 K in the presence of backwarming, a natural outcome of the bulk disk behind the inner wall. This resembles the disk sizes around higher-mass Herbig Ae stars (Millan-Gabet et al. 2007; Dullemond & Monnier 2010) and does not require any additional mechanism. The result is qualitatively consistent with recent work that shows the interferometric radius of T Tauri star disks consistent with dust sublimation radii if scattered light is considered (Anthonioz et al. 2015).

Harries (2011) built a radiative transfer model based on a flared disk geometry and computed the expected time lags over a range of wavelengths. By coincidence, the starting point was a hypothetical star with properties similar to those of YLW 16B: ${T}_{\mathrm{eff}}\sim 4000$ K, ${L}_{\ast }\sim 1$ ${L}_{\odot }$, and inclination $i\;=\;60^\circ$. In his model, the radius of the inner disk rim is assumed as 0.093 au. The response signal at 4.5 μm is a combination of thermal and scattered emission from a range of stellocentric radii in the disk. The model predicts a time lag at 4.5 μm of ∼180 s, about twice the round-trip light-travel time between the protostar and the disk inner rim. If this relation holds for YLW 16B, given the possible range of the total luminosity, the disk inner rim would have a radius of 0.037 au with temperatures of 2200–2400 K (Dullemond & Monnier 2010), probably inside the corotation radius and much hotter than the typical dust sublimation temperature of ∼1500 K. If this were the case, the dust in the inner disk would have to be either highly refractory or short-lived before sublimation. The disk would be truncated by the magnetosphere very close to the protostar and flared at all stellocentric distances as assumed in the Harries (2011) model. Otherwise, the disk may be self-shadowed with only a limited radial span exposed to stellar insolation, in which case the Harries (2011) model is not applicable, and the observed time lag at 4.5 μm simply corresponds to the light-travel time to the effective radius of the inner disk rim.

To identify the disk geometry of YLW 16B, we noticed that the SED of the system suggests some degree of self-shadowing (Figure 8(D); also see the detailed modeling in Lahuis et al. 2006). Here we additionally test whether the widths of the CCFs can be used to constrain the radial extent of the inner disk (see Isella & Natta 2005; Tannirkulam et al. 2007). A considerable radial span of the inner disk that is exposed to the starlight should lead to a superposition of various time lags and thus broaden the widths of the CCFs.

A problem specific to our case of YLW 16B is that each pair of wavebands has different overlaps of time coverage, while the lengths of the overlaps and the waveform of the light curve therein could be determining factors in yielding different widths of the CCFs. Therefore, we carried out the test by computing the CCFs based on the second night's data between 12131.0 and 19179.2 s in BMJD_TDB, where we have the longest uninterrupted time coverage overlap in all wavebands (Figure 3). The standard deviations of the Gaussian fits of the CCFs over this time are fairly close: 520 s for $H\star K$, 540 s for $H\star [4.5]$, and 510 s for $K\star [4.5]$. The similarity of all three CCF widths suggests the lack of a significant "broadening effect" at 4.5 μm, as would be expected from a fully flaring disk (Dullemond & Dominik 2004; Harries 2011), and instead favors the geometry of a self-shadowed disk.

For a more quantitative evaluation, we created a set of artificially smoothed light curves in K by computing the unweighted moving average of the original K-band time series over a window of varying length of time and computed their CCFs with the original H-band data. The purpose of using both H and K is to mimic a time series with superposed time lags and independent noise and assess the sensitivity of the CCF width to the degree of superposition, because the variations in the original H and K time series are both expected to be from the accreting gas with no time lag.

As a result, given the limited waveform and common time coverage, we found the CCF width fairly insensitive to time lag superposition. If rounded to the nearest 10 s, the standard deviation of the $H\star \;\mathrm{smoothed}\_K$ CCF does not reach 530 s (the next possible value broader than that of the original $H\star \;K$) until K is averaged over a window $\gt 400$ s. To make the width of $H\star \;\mathrm{smoothed}\_K$ greater than 580 s (i.e., 3σ broader than those of the original CCFs), the smoothing window has to be >950 s, which is the upper limit of any "broadening effect" at 4.5 μm. If the thermal reprocessing time of the disk is negligible, the size of the smoothing window translates to a 3σ upper limit of ∼1 au for the starlight-exposed radial extent of the disk inner edge. This limit is much smaller than the extent of the whole disk, as indicated by the very red WISE color $W3-W4\;=\;2.83$, corroborating a self-shadowed geometry over intermediate disk radii. However, this does not eliminate the possibility of optically thin gas and/or small amounts of refractory dust inside the rim (e.g., Tannirkulam et al. 2008; Benisty et al. 2010; Fischer et al. 2011).

As an explanation for larger interferometric disk hole sizes, it has been suggested that scattered light may inflate the interferometric measurements of the sizes of the inner disk cavities around T Tauri stars (Pinte et al. 2008). However, the reverberation results for YLW 16B are incompatible with such scattered light models, which rely on the assumption of a flaring disk. The models require a large radial extent of the disk to contribute to the scattering (Pinte et al. 2008), fundamentally incompatible to a self-shadowed disk geometry. Such diversity in disk geometry might also explain why the "composite disk model," which considers both thermal and scattered emission, provides better fits to the interferometric data of some T Tauri stars but is ineffective in other cases (Anthonioz et al. 2015).