TESTS OF THE PLANETARY HYPOTHESIS FOR PTFO 8-8695b

We conducted time-series photometric observations of fading events between 2012 and 2015. Below, in Section 2.2, we present ground-based observations of 26 different events. Figure 1 shows their distribution in time, and Figure 2 shows the light curves (including the seven highest-quality light curves presented previously by VE+12, for reference). In a few cases we observed the event through multiple broadband filters. The resulting multi-band light curves are shown in Figure 6. We also observed a candidate occultation at infrared wavelengths with one of the Magellan 6.5 m telescopes; those data are described in Section 2.3.1 and plotted in Figure 7. Finally, we analyzed the available Spitzer data, spanning a fading event as well as an expected occultation. Those data are described in Section 2.3.2 and shown in Figure 9. The dates of the Spitzer and Magellan observations are also indicated on Figure 1, along with the Keck spectroscopic observations described in Section 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We observed 13 fading events with the 1.2 m telescope at the Fred Lawrence Whipple Observatory (FLWO) on Mt. Hopkins, Arizona. The instrument, Keplercam, has a single 4096 × 4096 CCD with a 231 field of view. All the events were observed through an i' filter. For the events between 2014 January 9 and 18, we interleaved the i'-band exposures with g'-band exposures, although the g'-band data were only useful in two cases. Calibration was performed using standard IRAF¹⁵ procedures, including bias and flat-field corrections. The time stamps were placed on the BJD_TDB system using the code by Eastman et al. (2010). Circular aperture photometry was performed with the Interactive Data Language (IDL).

Another 13 events were observed with the 0.6 m TRAnsiting Planets and PlanetesImals Small Telescope (TRAPPIST), located at ESO's La Silla Observatory in Chile. This telescope is equipped with a thermoelectrically cooled 2048 × 2048 CCD with a 22' field of view (Gillon et al. 2011; Jehin et al. 2011). The observations were conducted with a custom "I + z" filter, which has transmittance >90% between 750 and 1100 nm. We refer the reader to Gillon et al. (2013) for descriptions of the procedures for observing and data reduction. The two events of 2012 December 14 and 15 were simultaneously observed with a Gunn r' filter, using EulerCam on the 1.2 m Euler-Swiss Telescope at the La Silla site, Chile. EulerCam uses a 2048 × 2048 CCD with a field of view of 147 on a side. For details on the instrument and data reduction procedures, please refer to Lendl et al. (2012).

A single event on 2014 January 19 was observed with the 6.5 m Magellan I (Baade) telescope at Las Campanas Observatory in Chile. The same event was observed simultaneously by TRAPPIST in the I + z band. With Magellan, we observed in the H band using FourStar, a 2048 ×2048 infrared array with a 109 square field of view. The data were reduced with IRAF and IDL procedures similar to those used on the FLWO data.

In all cases, the flux of PTFO 8-8695 was divided by the summed flux from several reference stars, leading to the light curves plotted in Figure 2. This figure also shows the seven light curves presented by VE+12 that cover the entire fading event; those observations were performed with the 1.2 m Palomar telescope and an R filter.

Outside of the fading events, the star varies gradually by ∼0.1 mag over several hours, in a manner consistent with its young age and late spectral type. Superimposed on those relatively gradual variations are periodic transit-like fading events lasting no more than about 2 hr. The depth and duration of the fading seems to vary from event to event. To derive the basic phenomenological parameters of the dimming events—depth, duration, and time of minimum light—we fitted a parameterized model describing both the gradual out-of-transit variations as well as the transit-like loss of light. We modeled the gradual variations as a polynomial function of time (2nd or 3rd order, depending on the event). The additional loss of light during the fading event, relative to the polynomial-corrected out-of-transit flux, was modeled as

$$\begin{eqnarray}&&{\rm{\Delta }}f(t)=\delta \mathrm{sech}\left[\displaystyle \frac{t-{t}_{0}}{w}\right]=\displaystyle \frac{2\delta }{{e}^{(t-{t}_{0})/w}+{e}^{-(t-{t}_{0})/w}},\end{eqnarray} \tag{ 1 }$$

where δ is the maximum fractional loss of light (the "transit depth"), w is the duration, and t₀ is the time of minimum light. We chose this model instead of a more physically motivated transit model (e.g., the model presented by Mandel & Agol 2002) because the asymmetries and other odd features in the light curves do not fit the standard models. Hence there is no advantage in fitting the physical model when a much simpler model can provide estimates of the basic transit parameters. One might be able to fit the data with a model based on transits of an oblate, oblique, precessing, gravity-darkened star (Barnes 2009), but such a model is far more demanding computationally. Our analytic model suffices to estimate the basic parameters of each event.

Figure 3 shows the light curves after dividing out the best-fitting polynomial functions. This gives a clearer view of the "transits" with most of the long-term trends removed. Table 1 gives the model parameters, as well as the value of ${\chi }_{{\rm{min}}}^{2}$ and the number of data points in each time series. In most cases, ${\chi }_{{\rm{min}}}^{2}$ is too large to be statistically acceptable, i.e., the simplified model of Equation (1) does not fit the data to within the photometric uncertainties. For this reason, Table 1 does not report the formal parameter uncertainties defined by the usual criterion Δχ² = 1. Rather, the reported parameter uncertainties have been enlarged by a factor of $\sqrt{{\chi }_{{\rm{min}}}^{2}/{N}_{{\rm{dof}}}}$ where N_dof is the number of degrees of freedom. These enlarged uncertainties were also adopted for our subsequent calculations.

Zoom In Zoom Out Reset image size

Table 1.  Best-fitting Parameters of Fading Events (see Equation (1))

UT Date	Filter	${\chi }_{{\rm{min}}}^{2}$	No. of	Time of Minimum Light, t0	Fractional Loss	Duration, w
			Data Points	[HJD−2455200]	of Light, δ	(days)
2010 Jan 5	R	193.4	228	1.8053 ± 0.0014	0.0442 ± 0.0052	0.0282 ± 0.0033
2010 Jan 6	R	216.5	244	2.6961 ± 0.0017	0.0789 ± 0.0208	0.0349 ± 0.0059
2010 Jan 9	R	217.0	191	5.8431 ± 0.0014	0.0979 ± 0.0142	0.0370 ± 0.0042
2010 Dec 9	R	169.9	270	339.9069 ± 0.0006	0.0498 ± 0.0023	0.0190 ± 0.0013
2010 Dec 10	R	248.4	298	340.8020 ± 0.0006	0.0425 ± 0.0011	0.0167 ± 0.0007
2010 Dec 13	R	191.1	271	343.9467 ± 0.0010	0.0558 ± 0.0038	0.0229 ± 0.0017
2010 Dec 14	R	473.2	288	344.8401 ± 0.0009	0.0412 ± 0.0018	0.0158 ± 0.0012
2012 Sep 4	I + z	194.0	209	974.8427 ± 0.0005	0.0510 ± 0.0027	0.0160 ± 0.0011
2012 Sep 13	I + z	388.8	265	983.8121 ± 0.0004	0.0521 ± 0.0014	0.0138 ± 0.0006
2012 Nov 23	I + z	240.1	243	1054.6614 ± 0.0009	0.0189 ± 0.0013	0.0128 ± 0.0014
2012 Dec 11	I + z	94.2	94	1072.5994 ± 0.0009	0.0366 ± 0.0024	0.0132 ± 0.0015
2012 Dec 14	I + z	154.0	176	1075.7379 ± 0.0006	0.0298 ± 0.0014	0.0098 ± 0.0008
2012 Dec 15	I + z	152.8	182	1076.6357 ± 0.0005	0.0338 ± 0.0029	0.0128 ± 0.0013
2013 Oct 12	i'	9.7	54	1377.9469 ± 0.0009	0.0324 ± 0.0024	0.0099 ± 0.0012
2013 Nov 23	I + z	102.7	93	1419.6524 ± 0.0022	0.0177 ± 0.0029	0.0114 ± 0.0030
2013 Dec 1	I + z	245.9	195	1427.7243 ± 0.0007	0.0316 ± 0.0026	0.0163 ± 0.0016
2013 Dec 11	I + z	190.3	233	1437.5864 ± 0.0007	0.0385 ± 0.0030	0.0153 ± 0.0016
2013 Dec 14	I + z	255.1	262	1440.7266 ± 0.0005	0.0344 ± 0.0016	0.0139 ± 0.0011
2014 Jan 9	i'	127.6	46	1466.7304 ± 0.0007	0.0385 ± 0.0016	0.0156 ± 0.0012
2014 Jan 14	i'	183.0	57	1471.6634 ± 0.0011	0.0337 ± 0.0025	0.0128 ± 0.0015
2014 Jan 17	i'	245.0	59	1474.8039 ± 0.0010	0.0377 ± 0.0042	0.0160 ± 0.0023
2014 Jan 18	i'	177.0	57	1475.7000 ± 0.0007	0.0375 ± 0.0014	0.0144 ± 0.0009
2014 Jan 19	I + z	253.6	240	1476.5974 ± 0.0007	0.0326 ± 0.0023	0.0167 ± 0.0016
2014 Jan 23	i'	331.9	136	1480.6289 ± 0.0004	0.0628 ± 0.0044	0.0218 ± 0.0013
2014 Jan 23	I + z	549.4	356	1480.6339 ± 0.0010	0.0295 ± 0.0018	0.0149 ± 0.0015
2014 Feb 5	i'	404.4	111	1493.6353 ± 0.0005	0.0349 ± 0.0013	0.0140 ± 0.0008
2014 Feb 9	i'	330.9	107	1497.6729 ± 0.0005	0.0402 ± 0.0013	0.0132 ± 0.0008
2014 Feb 13	i'	246.4	99	1501.7075 ± 0.0005	0.0451 ± 0.0022	0.0147 ± 0.0010
2014 Feb 22	i'	476.4	83	1510.6778 ± 0.0007	0.0463 ± 0.0027	0.0148 ± 0.0013
2014 Nov 17	i'	208.1	119	1778.7654 ± 0.0007	0.0304 ± 0.0020	0.0116 ± 0.0012
2014 Nov 29	i'	259.0	145	1790.8761 ± 0.0009	0.0185 ± 0.0014	0.0088 ± 0.0011
2014 Dec 27	i'	199.7	118	1818.6748 ± 0.0007	0.0222 ± 0.0012	0.0094 ± 0.0008
2015 Feb 27	I + z	200.0	206	1880.5567 ± 0.0009	0.0206 ± 0.0017	0.0087 ± 0.0012

Download table as:  ASCIITypeset image

Figure 4 shows the measured depths and durations. When the transits are deeper, they also tend to have longer durations; the measured depth and duration have a Pearson correlation coefficient of 0.84 (p < 10⁻⁵).

Zoom In Zoom Out Reset image size

The VE+12 light curves showed strong asymmetries in at least five out of the seven complete light curves. None of our new light curves show strong asymmetries, at least not as clearly as was seen by VE+12. We tried fitting a model in which the rates of brightness variation are not symmetric about the time of minimum light, by using a non-standard variant of the hyperbolic secant function (see, e.g., Ruan et al. 2000):

$$\begin{eqnarray}&&{\rm{\Delta }}f=\displaystyle \frac{2\delta }{{e}^{(t-{t}_{0})/{w}_{1}}+{e}^{-(t-{t}_{0})/{w}_{2}}}.\end{eqnarray} \tag{ 2 }$$

The asymmetric model does not seem to improve the quality of the fit to a significant degree. The number of cases for which the fitted asymmetry obeyed w₁ > w₂ (more prolonged "ingress") was nearly the same as the number of cases with w₂ > w₁, without any obvious pattern. Figure 5 shows the FLWO and TRAPPIST data as a function of $t-{t}_{0}$ (coverted to hours) after some averaging in time to increase the signal-to-noise ratio. In these averaged light curves there does seem to be a slight asymmetry, with a longer "ingress" than "egress" in both cases. The best-fitting asymmetric model is shown in Figure 5, and the parameters are given in Table 2.

Zoom In Zoom Out Reset image size

Another finding is that the loss of light is usually strongly chromatic, as shown in Figure 6. For those two cases in which we observed the same event in both the i' and g' bands, we found the loss of light to be 30%–50% larger in the g' band. For the single case in which we observed in both H and I + z, the loss of light was also ∼40% larger in the bluer band. We also observed two events simultaneously in the Gunn r' and I + z bands. On the first night the loss of light was 20%–30% larger in the bluer band. On the second night, the loss of light in r' was essentially the same as in I + z.

Zoom In Zoom Out Reset image size

The infrared emission from close-in giant planets is often strong enough that it is possible to detect the decrement in total flux when the planet is hidden by the star. The detection of such an occultation would be strong evidence supporting the planet hypothesis. We searched for the predicted occultations in the Magellan and Spitzer time-series infrared photometry, as discussed below.

2.3.1. Magellan Observations

The 6.5 m Magellan I (Baade) telescope and FourStar were also used to monitor PTFO 8-8695 for approximately 6 hr on 2014 January 21, spanning the predicted time of a planetary occultation. The prediction was based on the assumption that the planet's orbit is circular, i.e., that the occultations occur exactly halfway between transits. The observations were conducted in the H band and the data were processed in the manner described in Section 2.2. Figure 7 shows the resulting light curve. No occultation signal is evident, although some gradual variability is seen.

Zoom In Zoom Out Reset image size

We determined an upper limit on the relative brightness of the planetary dayside by fitting an occultation model to the data. The model had the same total duration and ingress/egress durations as the "transit" light curve observed with Magellan two days earlier. As with the transits, we fitted the gradual variations with a polynomial function of time (in this case, a cubic function). We used a Markov Chain Monte Carlo (MCMC) algorithm to calculate the posterior probability distribution for the parameters describing the cubic function as well as the loss of light during the occultation. The minimum χ² value was 2082.5, with 1876 degrees of freedom, indicating a statistically unacceptable fit; the cubic function is evidently not a completely satisfactory description of the observed flux variations. Rather than develop more elaborate models we simply inflated the parameter uncertainties by the factor $\sqrt{{\chi }_{{\rm{min}}}^{2}/{N}_{{\rm{dof}}}}.$ The resulting occultation depth was δ_occ = 0.00024 ± 0.00016, corresponding to a 3σ upper limit of δ_occ < 0.00072.

To decide if this upper limit rules out the planetary hypothesis, we need to know the expected occultation signal. Following the usual simplified model for transiting planets, the fractional loss of light during an occultation is

$$\begin{eqnarray*}&&{\delta }_{{\rm{occ}}}={A}_{g}{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{2}+{\left(\displaystyle \frac{{R}_{p}}{{R}_{\star }}\right)}^{2}\displaystyle \frac{{\displaystyle \int }_{{\lambda }_{1}}^{{\lambda }_{2}}{B}_{\lambda }({T}_{p})d\lambda }{{\displaystyle \int }_{{\lambda }_{1}}^{{\lambda }_{2}}{B}_{\lambda }({T}_{\star })d\lambda }.\end{eqnarray*}$$

The first term is due to reflected starlight, in which A_g is the geometric albedo, and R_p and R_⋆ are the planetary and stellar radii. The second term is due to the planet's thermal emission, in which the observing bandpass extends from λ₁ to λ₂ (1.48–1.76 μm, for H band), B_λ(T) is the Planck function, T_⋆ is the star's effective temperature, and T_p is the planet's dayside effective temperature. The latter is calculated from the condition of radiative equilibrium, giving

$$\begin{eqnarray}&&{T}_{p}={T}_{\star }{\left(\displaystyle \frac{{R}_{\star }}{a}\right)}^{1/2}{\left(\displaystyle \frac{1-A}{f}\right)}^{1/4},\end{eqnarray} \tag{ 3 }$$

where A is the Bond albedo and f is a dimensionless number depending on the manner of radiation. If the entire surface radiates isotropically as a blackbody, then f = 4. If instead the dayside radiates uniformly and the nightside radiation can be neglected, then f = 2. Furthermore, if the angular dependence of the planet's radiation is assumed to follow Lambert's law, then A_g = 2A/3.

In this case it is difficult to establish the key parameters ${({R}_{p}/{R}_{\star })}^{2}$ and R_⋆/a, because of the changing morphology of the fading events with time and wavelength. Barnes et al. (2013) found ${({R}_{p}/{R}_{\star })}^{2}\approx 0.027$ and R_⋆/a ≈ 0.58 using a model incorporating the effects of gravity darkening and orbital precession. Using those parameters, we calculate the expected value of ${\delta }_{{\rm{occ}}}$ and plot it as a function of the Bond albedo in Figure 8, for both f = 4 and f = 2. The expected δ_occ ranges from a minimum value of 0.0028 for f = 4 and A = 0, to a maximum value of 0.0061 for f = 2 and A = 1. Such large occultation depths are ruled out by our Magellan observations.

Zoom In Zoom Out Reset image size

2.3.2. Spitzer Observations

The Spitzer Space Telescope monitored PTFO 8-8695 on 2012 April 28 for about 12 hr, slightly longer than a full photometric period. The data consist of 1369 full-array images from the IRAC detector operating at 4.5 μm with an integration time of 26.8 s. The data were calibrated by the Spitzer pipeline version S19.1.0. These observations were carried out in non-cryogenic mode under program No. 80257 (PI: Stauffer), and are publicly available on the Spitzer Heritage Archive database.¹⁶

We converted the data from the Spitzer units of specific intensity (MJy sr⁻¹) into photon counts, and then performed IRAF aperture photometry on each subarray image. Best results were obtained with an aperture radius of 2.5 pixels and a background annulus extending from 11 to 15.5 pixels from the center of the point-spread function (PSF). The center of the PSF was measured by fitting a two-dimensional Gaussian function to each image. At this stage, 30 discrepant fluxes were discarded by applying a 5σ median clipping algorithm.

Next we needed to remove the apparent flux variations associated with motion of the image on the detector, the main source of systematic effects in time-series photometry with the IRAC InSb arrays (Knutson et al. 2008). This effect is caused by the combination of (i) the coarse sampling of the PSF, (ii) the significant inhomogeneity of the pixels, and (iii) fluctuations in the telescope pointing. To mitigate this effect we chose the Bi-Linearly Interpolated Sub-pixel Sensitivity (BLISS) mapping method presented by Stevenson et al. (2012). In this method, the flux data themselves are used to constrain a model for the subpixel sensitivity variations. In our implementation of the method, the detector area probed by the PSF center was divided into a 13 × 13 grid. With this degree of sampling, the PSF center visited each grid point at least 10 times throughout the course of the observations.

Figure 9 shows the Spitzer light curve after BLISS correction. The variability can be described as the combination of quasi-sinusoidal variation with a period of ≈0.5 day, and the transit-like dip in brightness at the expected time, with an amplitude of 0.5% and a duration of approximately 1.4 hr. No occultation is seen at the expected time (0.224 days after the transit). Figure 10 gives a better view of the transit-like event, and the data surrounding the predicted time of occultation. In this figure, only the data within 0.1 days of each event are shown, and the data have been rectified by fitting a quadratic function of time to data outside of the event and then dividing by the best-fitting function.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Just as with the Magellan light curve, we determined an upper limit on the occultation loss of light by fitting a parameterized model to the data. The model included a quadratic function of time to describe the out-of-occultation variations. The occultation model was required to have the same durations between first, second, third and fourth contacts as observed earlier with Magellan. The loss of light δ_occ was a free parameter. We used an MCMC algorithm to calculate the posterior probability distribution of δ_occ and the parameters of the quadratic function. The minimum χ² value was 565.9, with 535 degrees of freedom. The result for the occultation depth was ${\delta }_{{\rm{occ}}}=-0.0008\pm 0.0006,$ i.e., the best-fitting value corresponds to a brightness increase rather than a loss of light. This corresponds to a 3σ upper limit of δ_occ < 0.0010. Again, as illustrated in Figure 8, the upper bound on δ_occ given by Spitzer is smaller than the occultation depth implied by the parameters of the Barnes et al. (2013) model.

Careful inspection of Figure 9 shows a candidate flux dip of centered around a time coordinate of 0.425, with an amplitude of ≈0.3%. One might be tempted to attribute this dip to the occultation of a planet on an eccentric orbit, for which the occultation need not be halfway between transits. However, the statistical significance of this dip is dubious, and the required value of the eccentricity would be e > 0.35, using Equation (33) of Winn (2010). Such a high eccentricity would be unprecedented and unexpected for a short-period planet. In general, giant planets with periods shorter than three days have nearly circular orbits, a fact that is attributed to the gradual action of tidal dissipation. Given the youth of the star, it is possible that there has not yet been sufficient time for orbital circularization; however, a higher eccentricity and a potentially smaller pericenter distance would also put the planet in even more danger of violating the Roche limit.

Our non-detections of occultation signals at both 1.7 and 4.5 μm bands rule out the existence of a planet that radiates like a blackbody in these two bands. We have not pursued more realistic models for the planetary emission spectrum, given that the atmospheric composition is unconstrained, but seems unlikely that atmospheric absorption features would suppress the planetary flux in both bands to such a degree that it would be undetectable in our data.

Perhaps the most important finding of all the photometric observations is that the fading events are not strictly periodic. The top panel of Figure 11 shows the residuals after subtracting the best-fitting linear function of epoch from the measured times of minimum light,

$$\begin{eqnarray}&&{t}_{n}={t}_{0}+{Pn},\end{eqnarray} \tag{ 4 }$$

for which ${\chi }_{{\rm{min}}}^{2}=15573$ with 35 degrees of freedom. This poor fit is the result of the large scatter (≈15 minute) of the residuals within each season, and the even larger deviation (≈1.3 hr) of the most recent season's residuals relative to the earlier data. The pattern of residuals suggests that the period was nearly constant up until the 2014/5 observing season, when the fading events began occurring earlier than expected. This apparent change in period or phase can be readily checked by gathering additional data over the next few seasons. The best-fitting parameters of the linear ephemeris are

$$\begin{eqnarray*}{t}_{0} & = & 2455201.832\pm 0.007\;{\rm{days}},\\ P & = & 0.448391\pm 0.000003\;{\rm{days}}.\end{eqnarray*}$$

In these expressions the uncertainties have been scaled up by a factor of $\sqrt{{\chi }_{{\rm{min}}}^{2}/{N}_{{\rm{dof}}}}$ to account for the statistically poor fit.

Zoom In Zoom Out Reset image size

We also tried fitting a quadratic function of epoch,

$$\begin{eqnarray}&&{t}_{n}={t}_{0}+{P}_{0}n+\displaystyle \frac{1}{2}\displaystyle \frac{\;{dP}}{\;{dn}}{n}^{2},\end{eqnarray} \tag{ 5 }$$

for which ${\chi }_{{\rm{min}}}^{2}=4980$ with 34 degrees of freedom. After enlarging the parameter uncertainties as described above, the best-fitting parameters are

$$\begin{eqnarray*}{t}_{0} & = & 2455201.790\pm 0.006\;{\rm{days}},\\ {P}_{0} & = & 0.448438\pm 0.000006\;{\rm{days}},\\ \;{dP}/\;{dn} & = & (-2.09\pm 0.25)\times {10}^{-8}\;{{\rm{days\; epoch}}}^{-1}.\end{eqnarray*}$$

The bottom panel of Figure 11 shows the residuals between the observed and calculated times. The implied fractional change in period per epoch, calculated as $\frac{1}{P}\frac{\;{dP}}{\;{dn}},$ is equal to $-4.66\times {10}^{-8}.$ If this period change were to continue steadily, the period would shrink to zero after ${P}_{0}/\dot{P}\sim {10}^{4}$ years.

The primary purpose of the spectroscopic observations was to seek evidence for the RM effect, the spectroscopic anomaly that is seen during a planetary transit due to stellar rotation. During a transit, a planet blocks different portions of the rotating stellar photosphere, leaving a particular imprint on the rotationally broadened stellar absorption lines. The exact shape and time development of the spectral deformations depend on the transit parameters, and in particular on the angle between the stellar rotation axis and the orbital axis as projected on the sky plane. For PTFO 8-8695, the planet hypothesis requires a large misalignment between these angles. We attempted to detect the RM effect and test that prediction.

After the initial data reduction we corrected for the blaze function by using calibration lamp exposures to estimate the blaze function for each order, and then fitting a linear function of wavelength to remove the residual variations and normalize the continuum to unity. Following the barycentric correction, all of the out-of-transit exposures were co-added to create a single spectrum with a higher signal-to-noise ratio. This was used for a final differential normalization, wherein the summed spectrum was subtracted from each observed spectrum, and a fourth-order polynomial was fitted to the residuals in each order. These polynomials were subsequently subtracted from the corresponding spectrum. This was done to minimize the potential influence of any time variations in the blaze function throughout the night. We verified that the details of this normalization process did not have a significant influence on the following analysis.

For each spectrum, we calculated the cross-correlation function (CCF) with reference to a synthetic spectrum. The synthetic spectrum was obtained from the PHOENIX database (Husser et al. 2013), for a star with T_eff = 3500 K, $\mathrm{log}g=3.5$ and solar metallicity. We selected the appropriate wavelength ranges for creating the CCF via a visual inspection. We needed to locate areas for which the normalization seemed reliable, and where there were at least a few well-defined absorption lines. As had been reported by VE+12, there are only a few regions between 5000 and 7000 Å suitable for this work.

We calculated the mean CCF based on all of the out-of-transit observations. Then we subtracted the mean CCF from each individual CCF. When ordered in time, the resulting "differential CCFs" should display the shadow of the transiting object in velocity space. The deformation due to a transiting object would be seen as a dark line. The slope of this line in the velocity-time plane would depend on the projected obliquity. For example, in the case of good spin–orbit alignment, there would be a deficit of blue light (negative radial velocities) in the first half of the transit, followed by a deficit of redshifted light during the second half. As can be seen in Figure 13, no such signal—neither aligned nor misaligned—can be discerned.

Zoom In Zoom Out Reset image size

Could we have detected the spectroscopic transit of the hypothetical planet, given the quality of our data? To answer this question, we simulated the RM signal of a transiting planet with ${({R}_{p}/{R}_{\star })}^{2}=0.026,$ the approximate transit depth measured by TRAPPIST during the events of 2013 December 11 and 14, closely bracketing the event observed with Keck. We assumed v sin i_⋆ = 100 km s⁻¹, a value consistent with the line broadening seen in our Keck spectra, and adopted a macroturbulent velocity of 15 km s⁻¹. Then we injected RM signals into the data, for various choices of λ, the sky-projected stellar obliquity.

Some representative examples of the simulated RM effect are shown in Figure 14. The left panel shows the simulated RM signal of a well-aligned planet (λ = 0°). Such a signal would easily have been detectable with the data at hand. The same is true for λ = 45°, shown in the middle panel. However, for λ = 90°, the signal would have been more difficult to detect. This is because in this case the RM signal is nearly stationary in velocity, as shown in the right panel, making it more difficult to separate from the noise in the velocity–time plane. Based on visual inspection of figures similar to Figure 14, we conclude that we can rule out any trajectory except for those within about 15^◦ of perpendicularity (λ = 90° or 270^◦). Given the non-Gaussian and correlated nature of the noise, it is difficult to make a firmer statistical statement.

Zoom In Zoom Out Reset image size

A secondary goal was to seek changes in the sky-projected rotation rate (v sin i_⋆) that would be expected if the star's rotation axis is precessing around the total angular momentum of the system. If PTFO 8-8695 does consist of a star and a planet whose rotation axes and orbits precess around the common angular momentum, then $v\mathrm{sin}{i}_{\star }$ should change with time. We searched for such a change between the two epochs for which Keck/HIRES data have been obtained. Five spectra were obtained in 2011 April and presented by VE+12. Another epoch is represented by our 2013 December data.

We derived CCFs from the 2011 data in exactly the same way as for the 2013 data. We compared the CCFs with theoretical absorption lines taking into account uniform rotation, limb darkening, macroturbulence, and gravity darkening. We adopted the quadratic limb darkening parameters from the tables of Claret et al. (2012) for $\mathrm{log}g=3.5$ and T_eff = 3500 K, solar metallicity, and the Johnson V band. We adopted a macroturbulent velocity of 15 km s⁻¹. To model the gravity darkening, we assumed the same effective temperature, a rotation period of 0.448 days, a stellar mass of 0.4 M_⊙, a stellar radius of 1.4 R_⊙, a stellar inclination angle of 90°, and a gravity-darkening exponent of $\beta =0.25.$ We neglected any oblateness of the stellar photosphere. As an example, Figure 15 shows the CCF for one of the pre-transit observations from 2013 December, along with the best-fitting model.

Zoom In Zoom Out Reset image size

The results for v sin i_⋆ are 103.6 ± 0.3 km s⁻¹ in 2011 and 104.1 ± 0.7 km s⁻¹ in 2013. The quoted uncertainties are based on the scatter between the different observations for each epoch, and do not include any additional systematic uncertainties due to the limitations of the model (such as uncertainties in the treatment of limb darkening and gravity darkening, or the neglect of differential rotation and oblateness). Therefore the relative variation in v sin i_⋆ is bounded to less than a percent, although the absolute value is probably uncertain by at least 10%.

To obtain a better idea about the absolute value of $v\mathrm{sin}{i}_{\star }$ and its uncertainty, we tried fitting individual absorption lines rather than the CCF. Specifically, we fitted seven apparently isolated lines between 5300 and 7700 Å. The standard deviation in the v sin i_⋆ measurement form these seven absorption lines varies between 7 and 15 km s⁻¹ for the different observations, which suggests that an uncertainty of 10 km s⁻¹ in $v\mathrm{sin}{i}_{\star }$ should be a reasonable estimate. This leads to our final estimate of 103 ± 10 km s⁻¹.

Our result for v sin i_⋆ is higher than the value of 80.6 ± 8.1 km s⁻¹ reported by VE+12, but we do not think that this necessarily (or even likely) implies that the projected rotation rate is varying in time. This is because it is difficult to compare the results directly, given that VE+12 used a completely different instrument and analysis technique. Our internal comparison is much more sensitive, since it is between two Keck/HIRES spectra obtained at different times and analyzed in exactly the same way.

Our optical spectra also reveal strong emission lines from the hydrogen Balmer series as well as the Ca ii H & K transitions. The top panel of Figure 16 shows the median Hα line profile, based on all the "out-of-transit" spectra observed on 2013 December 12 (at least an hour before or after the time of minimum light). The line is very broad. Most of the emission is confined to velocities ≲100 km s⁻¹, consistent with the star's rotation rate, but the velocity profile extends to at least 300 km s⁻¹, particularly on the blue side. This is suggestive of at least a low level of ongoing accretion. Material that falls onto the star from large distances could attain the free-fall velocity $\sqrt{2{{GM}}_{\star }/{R}_{\star }}\approx 330$ km s⁻¹, given the nominal parameters M_⋆ = 0.4 M_⊙ and R_⋆ = 1.4 R_⊙. The equivalent width of the Hα line is 8.7 Å, placing the star near the traditional borderline between the categories of "weak-lined" and "classical" T Tauri stars. For simplicity it is often said that the classical stars are actively accreting while the weak-lined stars are not accreting, although in reality there seems to be no sharp distinction between these categories, and PTFO 8-8695 presents an intermediate case.

Zoom In Zoom Out Reset image size

Figure 16 also shows the time sequence of observed changes in the Hα line profile. Specifically, for each of the 22 spectra, we plotted the residuals between the observed line profile and the median "out-of-transit" line profile. Also indicated are the times of minimum light, as well as "ingress" (one hour prior) and "egress" (one hour afterward). The red dashed line indicates a slope of 37 km s⁻¹ hr⁻¹. This is the expected radial acceleration of any feature attached to the stellar photosphere, which would move from −100 to +100 km s⁻¹ over the course of P_rot/2 = 0.224 days.

Evidently the line profile varied in a complex pattern on a timescale of minutes. Several excess-emission features do seem to be rotating along with the star; for example, a pattern of positive residuals appears in the fourth-to-final spectrum (time coordinate t = 1.015) at velocity +50 km s⁻¹ and shifted redward at the expected rate throughout the final three observations. A similar pattern—perhaps originating from the same feature—is seen starting at minimum light at velocity −50 km s⁻¹ and proceeding redward until about an hour after minimum light. These particular components of the emission line seems likely to be caused by active regions on the stellar surface.

During the fading event, the residuals show relative absorption at a redshifted velocity of 25–100 km s⁻¹. The absorption seemed to disappear at around the same time as the end of the fading event. The onset of the absorption feature was at least 2 hr before minimum light, which is at least an hour before what seems to be "ingress" of the fading event. Thus, the transient redshifted absorption does not seem to be exactly coincident with the fading event, although it does at least suggest that the fading event was associated with hot infalling material in front of the star.

The Hγ line profiles (not shown here) tell a similar story but with a lower signal-to-noise ratio. The Hβ line was not observed, given the spectral format. Figure 17 shows the median Ca ii H line profile, along with the time series of deviations from the median. In this case the fractional variations were even stronger and seemingly faster; there is no straightforward narrative to the sequence of residuals.

Zoom In Zoom Out Reset image size

To seek independent evidence for ongoing accretion, we used the available broadband photometry to construct the spectral energy distribution of PTFO 8-8695 and search for any infrared excess. Figure 18 shows the result, based on a query of the VizieR website,¹⁷ which gave measurements in the VRJHK bands as well as the WISE W1–W4 bands. (The W4 observation gave an upper limit.) We corrected for dust extinction with the dust map from the NASA/IPAC website,¹⁸ and fitted the results to a grid of zero-metallicity stellar atmosphere models from the library of Castelli & Kurucz (2004). The best-fitting stellar parameters were T_eff = 3500 ± 120 K and log g = 4.0 ± 0.9. The effective temperature is in agreement with the previously reported value of 3470 K (Briceño et al. 2005). The apparent lack of an infrared excess out to 10 μm is characteristic of weak-lined T Tauri stars. This lack of evidence for the existence of an accretion disk within 1 AU does not necessarily rule out accretion, but does suggest that any accretion is relatively weak.

Zoom In Zoom Out Reset image size

The seemingly erratic variations in the depth and duration of the fading events, along with the slight asymmetries in the phase-folded light curves (Figure 5), bring to mind the case of KIC 12557548 (Rappaport et al. 2012). This object was identified in data from the Kepler spacecraft. It exhibits transit-like flux dips with a very short period (15.7 hr) and duration (1.5 hr), an erratically varying depth ranging from 0%–1.3%, and an egress of longer duration than ingress. Rappaport et al. (2012) interpreted the dips as transits by a dusty tail being emitted by a small rocky planet. Two other similar cases have since been identified (Rappaport et al. 2014; Sanchis-Ojeda et al. 2015). In at least one case the transit depth has been shown to be smaller at infrared wavelengths, as expected for extinction by small dust grains.

PTFO 8-8695 shares much of the phenomenology that has just been described. The main difference is that the other systems are not rapidly rotating young stars. They appear to be slowly rotating main-sequence stars, and are not even close to being synchronized with the transit period. Furthermore it is not clear whether this hypothesis could be reconciled with the apparent change in period that was seen in the most recent season of observations.

It is tempting to invoke tidal dissipation as a mechanism for gradually shrinking the orbit and shortening the period, but this is implausible for a low-mass rocky body. Conceivably, orbital decay could be a consequence of the dust emission. If the dust acquires additional specific angular momentum from radiation pressure while leaving the system, it would be driven into a higher orbit. There it would pull back on the planet and potentially extract angular momentum from the planetary orbit. However, the magnitude of this effective drag is difficult to estimate from first principles, particularly because the dust may represent only a modest fraction of the total mass loss, and the gas need not behave the same way as dust as it leaves the system. In any case the lifetime of ∼10⁴ years implied by the observations (see Section 2.4) is uncomfortably short. It would require a special coincidence to observe such a short-lived phase of evolution.

The synchronization of the stellar rotation period and the period of the fading events, along with the changing depth and duration of the fading events, raises the suspicion that the periodic dips are caused by starspots being carried around by rotation. The star is expected to be heavily spotted, given its youth. Moreover, the photometric variations produced by stellar activity are expected to be weaker in the infrared than in the optical, consistent with our observations. Gradual changes in the spot pattern could be invoked to explain the changes in depth, duration, and timing of the fading events.

VE+12 have already pointed out the main weakness of this hypothesis. Flux variations caused by starspots have a natural timescale of half the rotation period, the interval over which a patch of the stellar surface is typically visible to the observer. In contrast, the fading events of PTFO 8-8695 last only ≈15% of the rotation period. One can arrange for short-duration dips by locating the starspots near one of the rotation poles, and tipping the star such that the spots are only on the visible hemisphere over a narrow range of longitudes. Indeed, T Tauri stars are well known for displaying stable spot patterns near the magnetic poles (e.g., Granzer et al. 2000). However, in such a circumstance the spot is highly foreshortened and limb-darkened, and it is difficult to produce variations of several percent. It is also difficult to produce the apparently sharp features that have been occasionally seen in the light curves, such as the apparent point of "fourth contact" in our Magellan H-band time series (see the bottom panel of Figure 6). It might be possible to save this hypothesis by invoking multiple spots and in a complex and stable pattern, such that the summation of their photometric variations is coincidentally narrow in time. This model, though, would be rather contrived.

As we have just seen, the basic difficulty of any model in which the variations are due to features on the stellar photosphere is the relatively short duration of the fading events. Eclipses by orbiting bodies avoid this problem because the time spent in front of the star scales as R_⋆/a, and the orbital distance a can be adjusted to match the observations. This is the basis of the planet hypothesis, for which Barnes et al. (2013) found ${R}_{\star }/a=0.5$–0.6. However the eclipsing body need not be a planet. Perhaps it is a feature in the innermost, corotating portion of the stellar magnetosphere.

At an age of a few million years, low-mass T Tauri stars are often still actively accreting material from their circumstellar disks. PTFO 8-8695 lacks any evidence of an accretion disk, yet, as noted in Section 3.3, the strength and breadth of the Hα line profile places the star in between the traditional categories of "classical" and "weak-lined" T Tauri stars and suggests that PTFO 8-8695 may be weakly accreting. Therefore, we will consider the implications of a small amount of dust that would not be detected as an IR excess, but would still be able to produce the flux dips. The accretion process is thought to proceed as follows (see, e.g., Bouvier et al. 2007, p. 479, for a review). First, matter spirals inward through a thin disk, until it reaches an orbital distance of a few stellar radii, at which point the disk is disrupted by the star's magnetosphere. Then, within the magnetosphere, the material is magnetically funneled into narrow tubes or columns, extending from the inner edge of the accretion disk onto the star's magnetic poles. The material falls freely along these columns and crashes onto the star, producing shock waves and a luminous hotspot. In many models, the stellar rotation rate becomes synchronized with the Keplerian orbital velocity at the inner radius of the magnetosphere (the corotation radius), a phenomenon known as disk locking.

This picture contains several elements that could naturally lead to flux variations with a period equal to the stellar rotation period. For example if the accretion disk is warped or has other non-axisymmetric variations in thickness near the innermost, corotating portion of the disk, then the star may be periodically eclipsed by these irregularities. This is thought to be the basic explanation for the quasiperiodic eclipses of AA Tau (Bouvier et al. 1999, 2003). Alternatively, periodic eclipses could be produced by stable patterns or concentrations of dust in the accretion flow. This may explain the observed variability of the "short-duration dippers" recently identified by Cody et al. (2014) and Stauffer et al. (2015). Those authors studied a number of T Tauri stars in the young cluster NGC 2264 that exhibit short-duration fading events. They found that the flux dips are quasi-periodic and exhibit changes in depth and shape from epoch to epoch over a period of several years. The dips appear shallower in the infrared than in the optical, and the light curves have rounded minima rather than being flat-bottomed (Stauffer et al. 2015). The flux dip periods are usually equal to the rotation periods of the stars, and some dips have been observed to persist for years. All these properties have been attributed to extinction by infalling material from the innermost portion of the accretion disk (McGinnis et al. 2015). And, all these properties are consistent with our observations of PTFO 8-8695.

In some respects, though, PTFO 8-8695 is different from the rest of the dippers. Its rotation period of 0.45 days is shorter than the 3–10 day periods of most of the stars observed by Cody et al. (2014) and Stauffer et al. (2015), or the 8 day period of AA Tau. The duration of the fading events is also relatively shorter at 15% of the period, compared to the more typical value of ∼30%. It is also unclear whether disk warps or dust concentrations could produce the apparently sharp features and flat bottoms observed in a few of our light curves. Furthermore, the dippers all have spectral energy distributions with a detectable infrared excess, but the currently available data for PTFO 8-8695 show no evidence for any infrared excess (Figure 18).

Furthermore, it is questionable whether dust can exist in solid form so close to the star, with an orbital distance less than 2 R_⋆. Assuming that stellar radiation is the dominant mechanism of heating, the dust sublimation radius R_s is given by (Monnier & Millan-Gabet 2002)

$$\begin{eqnarray}&&{R}_{s}=\displaystyle \frac{1}{2}\sqrt{{Q}_{R}}{\left(\displaystyle \frac{{T}_{\star }}{{T}_{s}}\right)}^{2}{R}_{\star },\end{eqnarray} \tag{ 6 }$$

where T_s is the dust sublimation temperature (≈1500 K), and ${Q}_{R}={Q}_{{\rm{abs}}}({T}_{\star })/{Q}_{{\rm{abs}}}({T}_{s})$ is the ratio of the dust absorption efficiencies for radiation at the color temperatures of the incident and reemitted fields, respectively. Adopting the customary value of T_s = 1500 K, and assuming Q_R to be of order unity (as expected in this case for silicate grains), this formula gives ${R}_{s}=2.7{R}_{\star }.$ Therefore the hypothesized dust with R_⋆/a = 0.5–0.6, corresponding to a = 1.7–2 R_⋆, would likely be vaporized.

Hotspots are another aspect of the magnetospheric accretion model that has previously been invoked to explain some of the quasiperiodic variations of T Tauri stars (Herbst et al. 1994). We advance here a related hypothesis for the case of PTFO 8-8695: perhaps the fading events represent occultations of one of the hotspots that is produced by ongoing low-level accretion. In this scenario, the accreting material is funnelled onto a hotspot near one of the star's magnetic poles, which is displaced from the star's rotation pole. Furthermore, the star's rotation pole is tipped toward the observer such that the hotspot is on the visible hemisphere for ≈85% of each rotation period. When the hotspot is hidden from view, we observe a fading event.

This is similar to the starspot hypothesis (Section 4.2) in that the photometric modulations are the result of the rotation of a photospheric feature of contrasting intensity. However, replacing the dark starspot with a luminous hotspot could potentially solve some of the problems faced by the starspot model. Hotspots have been inferred to have effective temperatures ranging up to 10⁴ K, and thereby present much higher contrast than starspots with the surrounding photosphere. Furthermore, hotspots are probably not confined to a vertically thin surface layer of the photosphere, and as such they are not subject to the effects of limb-darkening or foreshortening. Together these factors may make it easier to produce abrupt modulations of a few percent in the total light even for a small feature near the stellar limb.

If this model is correct, then the mass accretion rate can be estimated from the loss of light during fading events, which is observed to be of order 5%. This requires the accretion luminosity to be

$$\begin{eqnarray*}&&{L}_{{\rm{acc}}}=\displaystyle \frac{{{GM}}_{\star }\dot{M}}{R}\sim 0.05{L}_{{\rm{bol}}}.\end{eqnarray*}$$

Using the nominal values M_⋆ = 0.4 M_⊙, R = 1.4 R_⊙, and L_bol = 0.25 L_⊙, we find $\dot{M}\sim {10}^{-9}$ M_⊙ yr⁻¹, at the low end of the range of mass accretion rates that has been inferred for accreting T Tauri stars (10⁻⁹–10⁻⁷ M_⊙ yr⁻¹). This seems reasonable: a relatively low accretion rate is also consistent with the relatively weak Hα equivalent width of 8.7 Å and the absence of a detectable infrared excess.

One concern with this model is that in the systems where hotspots have been previously invoked to explain photometric variability, the pattern of variability is not as consistent or long-lasting as is seen in PTFO 8-8695 (Herbst et al. 1994). Oftentimes the hotspot variability shows no periodicity or, at best, short-lived periodicity, sometimes with period changes of 20% or more. However, those previous results pertain to classical T Tauri stars with much higher inferred accretion rates; perhaps we are seeing the different and more stable behavior of a hotspot in a more weakly accreting system.

We have discussed five hypotheses for the fading events of PTFO 8-8695. The precessing giant-planet model has failed several key tests, the most serious of which are probably the nondetection of the planetary occultations, and the apparent shift in the phase of the fading events in the most recent season of observations. The planet hypothesis also struggles to explain the observed coincidence of the rotation and orbital periods. This same problem afflicts the hypothesis of the dust-emitting rocky planet.

The other three models share the virtue of a natural explanation for the equality of the rotation period and the period of the fading events. However, the starspot model has difficulty reproducing the observed duration and occasionally sharp ingress/egress of the fading events. The other two models invoke the presence of an accretion disk, for which the evidence is ambiguous or negative: the Hα line profile does extend to higher velocities than can be explained by stellar rotation, but the equivalent width is relatively low and there is no detectable infrared excess.

The occulted-hotspot model seems quantitatively promising, as it is consistent with a low accretion rate of ∼10⁻⁹ M_⊙ yr⁻¹. There is no deterministic theory for the expected photometric variations due to magnetically funnelled accretion, making it difficult to achieve a firm confirmation of this hypothesis. Nevertheless all our observations seem at least consistent with this picture. At present, we consider this hypothesis to be the best explanation for PTFO 8-8695. To come to a firmer conclusion will probably require more photometric and spectroscopic observations, seeking changes in the timing and behavior of the fading events, variations in the Hα line profile, and more sensitive searches for any infrared (or ultraviolet) excess or other indicators of low-level accretion.

At the outset of this project, and throughout most of this paper, we have been mainly concerned with the status of the planetary hypothesis for this intriguing planetary candidate. In fact this object may turn out to be useful for understanding magnetospheric accretion, due to a fortuitous geometry, thereby joining the ranks of such systems as AA Tau (Bouvier et al. 1999) and KH 15D (Hamilton et al. 2012).