THE K2-ESPRINT PROJECT. I. DISCOVERY OF THE DISINTEGRATING ROCKY PLANET K2-22b WITH A COMETARY HEAD AND LEADING TAIL

Our photometric pipeline follows the steps of similar efforts published to date (Vanderburg & Johnson 2014; Aigrain et al. 2015; Crossfield et al. 2015; Foreman-Mackey et al. 2015; Lund et al. 2015) that describe how to efficiently extract light curves from the calibrated pixel level data archived on MAST. The ingredients to generate the light curves, with our own choice described, are:

• 1.  
  Aperture selection: Our apertures have irregular shapes, which are based on the amount of light that a certain pixel receives above the background level. We selected this type of aperture to capture as much light as possible while reducing the number of pixels used, which in turn reduces the noise induced by a large background correction. Based on experiments we carried out on the engineering data release, selected pixels must be 30% higher than the background level (estimated from an outlier corrected median of all the pixels in the image) in 50% of the images in the case of a star brighter than Kepler magnitude K_p = 11.5. For stars fainter than magnitude K_p = 14, the selected pixels must be 4% higher than the background. For stars of intermediate magnitudes, a linear interpolation of these two thresholds is used. A simple algorithm groups contiguous pixels in different apertures, and the target star aperture is selected to be the one that contains the target star pixel position (obtained from the FITS headers). This type of aperture is similar to the ones used by Lund et al. (2015), and quite different from the circular apertures used in many of the other pipelines.
• 2.  
  Thruster event removal: As highlighted in Vanderburg & Johnson (2014), every 6 hr the telescope rolls to maintain the targets on the defined set of pixels that are downloaded, in what is known as a "thruster event." We recognized these events by calculating the centroid motion of a particularly well behaved star (EPIC 201918073), and selected those moments where the position of the star jumps much more than usual (in the case of these stars, 0.1 pixels in the x direction). The images obtained during thruster events are removed from the analysis.
• 3.  
  Data slicing: We split our dataset into eleven different segments chosen to have a length of approximately 7 days, but also to contain an integer number of telescope roll cycles. The first segment and the one after a large gap (in the middle of the dataset) are not used in the global search since they are poorly behaved in some cases, with systematic effects induced by thermal changes that appear after reorienting the telescope.
• 4.  
  Systematics removal: We calculate the centroid positions and obtain a fourth-order polynomial that describes the movement of the star in the x and y coordinates. This polynomial fit is used to determine a new set of coordinates, in which the star moves only along one direction. The fluxes are then decorrelated first against time and then against this moving coordinate, using a fourth-order polynomial in each case. This process is repeated three times, and the results are very robust against problems caused by the presence of low-frequency astrophysical sources of noise (see Vanderburg & Johnson 2014 for a more detailed description on how this process works).

This recipe was followed to generate the light curves of the 21,647 stars in Field-1. Among them, the light curve for EPIC 20163717 can be seen in the top panel of Figure 1. It is interesting to note that since our apertures depend on the amount of background light, which is increasing over the course of the observations, the number of pixels in the aperture decreases with time. This is a desired effect, since it tends to balance the increase of background light in the aperture by reducing the number of pixels, and therefore changes in the flux scatter are less severe.

Zoom In Zoom Out Reset image size

These light curves are generally analyzed using two different search algorithms: a more standard BLS routine (Kovács et al. 2002; Jenkins et al. 2010; Ofir 2014) to search for planets with orbital periods longer than 1 day, and a more specialized FFT pipeline used to detect planets with orbital periods shorter than 1–2 days (Sanchis-Ojeda et al. 2014). In this case, due to the short orbital period of the signal, we describe only the FFT search.

Our target star K2-22 was detected as part of our search for ultra-short period planets in the K2 Field-1 dataset using the FFT technique. The routine automatically identifies those objects for which a main frequency and at least one harmonic can be distinguished above the level noise in the FFT power spectrum (see Sanchis-Ojeda et al. 2014 for details). A total of 2628 objects were identified in that way, but a large fraction of them were caused by improper corrections of the 6-hr roll of the Kepler telescope. These false detections are easy to remove since their main frequency is always related to the fundamental roll frequency of 4.08 rolls per day, although this simple removal clearly affects the completeness of our search. A total of 390 objects were selected for visual inspection, and among them K2-22 was selected as the most promising ultra-short orbital period planet candidate.

The FFT by which this object was discovered is shown in Figure 2. Note the prominent peak at 2.62 cycles/day which is the base frequency corresponding to the 9.1457-hr period, as well as the next 8 higher harmonics which lie below the Nyquist limit. The overall slowly decaying Fourier amplitudes with harmonic number is characteristic of short-period planet transits.

Zoom In Zoom Out Reset image size

Our method for generating the light curves relies on a one-to-one relationship between the raw flux counts and the position of the star on the CCD chip. Any source of astrophysical variability could distort this relationship, and this is the case for both stellar activity induced signals and transits. A closer inspection of the raw light curve of K2-22 shows long-term trends that do not correlate with the centroid motion, and are likely due to the slow rotation of the host star. These trends are removed automatically as we fit a fourth-order polynomial in time to each of the 7-day segments, which effectively removes variability on scales longer than approximately 2 days (see upper panel of Figure 1).

In order to remove the effects of the transits, we first folded the original light curve given the period obtained from the FFT, after which we identified those orbital phases where no transit is expected. We then ran our photometric pipeline again, using the same aperture (see Figure 3) but only using the out-of-transit flux measurements to find the best fit polynomials to correct for both temporal and telescope motion variations. This process reduced the photometric scatter, and encouraged us to try different approaches to continue improving the light curve. We tried different combinations of polynomial orders and also different approaches to defining the apertures, but none of them improved the quality of the light curve (see top panel of Figure 1). After all the corrections, the final scatter per 30 minute cadence is 650 ppm, which is near the mean uncertainty obtained with our photometric pipeline for the typical K2 15th magnitude star.

Zoom In Zoom Out Reset image size

We also tried to produce a light curve in which other astrophysical signals would be preserved (e.g., starspot rotation). During the process of detrending each of the 11 segments of data, we saved the coefficients of the fourth-order polynomial in time, and used them to reconstruct the signal again after removing the centroid motion artifact. Since the aperture is individually defined for each segment, we had to adjust the mean flux level of each segment to create a continuous light curve. This astrophysically more accurate light curve is also shown in Figure 1, and exhibits a clear signal of starspots with a rotation that could either be 7–8 days or twice this value. The shorter quasi-periodicity would typically arise when the star has two active longitudes separated by 180° in longitude. We used an autocorrelation function to confirm this suspicion (see lower panel of Figure 1), and measured a rotation period of 15.3 days, following the techniques described in McQuillan et al. (2013).

The top panel of Figure 1 shows the full data set from the K2 Field-1 observations covering an interval of ∼80 days. It is apparent that there are sharp dips in intensity whose depths are highly and erratically variable. In Figure 4 we can see a zoom in on two weeks of observations. The individual transits are now quite apparent, and the depth variations are dramatically evident.

Zoom In Zoom Out Reset image size

In order to analyze these variations, we first folded the light curve with the period obtained from the FFT, after removing any long-period signals. We did this by fitting for local linear trends using a total of 5 hr of observations right before and after each transit. We then fit the folded light curve with a simple idealized transit profile comprised of three straight-line segments: a flat bottom and sloping ingress and egress with the same slope magnitudes, hereafter referred to as the "three-segment (symmetric) model." This is very similar to the simpler "box model," but with non-zero ingress and egress times. We forced the fit to have a duration of at least 5 minutes on the lowest part of the transit, to make sure that the final fit did not look like a triangle. This fit gave a mean transit depth of 0.5%, and a total duration of 1 hr and 10 minutes. Given that the cadence of the observations is approximately 30 minutes, the real mean duration of the K2 transits must be close to 40 minutes, in agreement with what is expected of an USP planet orbiting an M-dwarf. This is confirmed by our follow-up ground-based observations with their better temporal resolution (see Section 4).

With the mean transit profile and a good estimate for the orbital ephemeris, we fit each of the 190 individual K2 transits with the same three-segment model, allowing only the transit time to vary. In the process, we evaluated the robustness of the detection of each of the individual transits. In some cases, the depth is so low that the transit cannot be detected, whereas on other occasions the transit occurs during a thruster event, so there are no data. In rare cases, the transit consists of only a single flux point, which does not allow for a clear determination of the transit time or depth. After removing all those cases, we were left with a reduced sample of 60 well detected transits. In addition to the best fit transit time, we estimated the formal uncertainty by finding the interval of times where the standard χ² function is within 1 of its minimum value, where a constant value of 650 ppm is used for the uncertainties in the flux. The transit times (including uncertainties) obtained in this manner are fit to a linear function, and the O – C ("observed minus calculated") residuals of that fit are displayed in Figure 5. A clear excess of scatter (above the formal statistical uncertainties) is detected, with a best-fit standard χ² of 480 for 60 transit times. Formal uncertainties, however, underestimate the true uncertainties when the model parameters are correlated, as could be the case here since we have removed a local linear trend for each transit which is known to be correlated with the transit time. We further examine this excess scatter below with simulations of the data train.

Zoom In Zoom Out Reset image size

We now use the new orbital ephemeris to fix the time of each transit, and repeat the process but now fitting for the transit depths. These transit depths are significantly different from one transit to the next, and the depths range from a maximum of 1.3%–0.27% which is close to the photometric detectability limit of our time series for individual flux points. We can also see in Figure 5 that the transit depths do not strongly correlate with the O – C timing residuals. These erratic transit depth variations are quite reminiscent of those exhibited by KIC 1255b (Rappaport et al. 2012; Croll et al. 2014), and cover much the same range in depths.

We checked for periodicities in the measured transit depths and timing residuals using a Lomb–Scargle periodogram, but no significant peak was detected with a false alarm probability smaller than 1%. We also attempted to constrain the true size of the planet by studying the shallowest K2 transits. We revisited the discarded transits, and selected the six shallowest cases in which the transit observations were complete (no thruster events). The mean depth of these transits is 0.14 ± 0.03%, which, given the radius of the star (see Section 5.2), translates into an upper bound of 2.5 ± 0.4 R_⊕ on the planet radius.

Even though the transit depth variations described above are later confirmed by ground-based observations, there is a potential concern that the variability could have been caused by the relatively comparable transit and sampling timescales. In this same regard, it is also possible that some or all of the excess variations in the O – C scatter (see Figure 5) over those expected from statistical fluctuations are due to the relatively short transit duration compared with the LC integration time. We have therefore carried out extensive numerical simulations of these effects.

The simulations of the depth and transit-time variations are based on a model that has a simple box transit profile which is integrated over the 30 minute LC time, and includes a sinusoid (and its first harmonic) to represent the rotating starspot activity. The amplitude of the sinusoids and the rotation period are fixed at 1% and 15 days, respectively. The duration of the model transit is fixed at 50 minute (see Section 4). The model is evaluated at the same times as the K2 observations, and only the same 60 transit windows are analyzed to ensure that the simulations represent a similar dataset to the one used in this paper. There are a total of four different numerical experiments, each one repeated many times, in which we either include starspots or not, and we either have a constant transit depth or depths drawn from a Gaussian distribution of a given variance, but always with a mean of 0.5%. In all cases a fixed orbital ephemeris is used. White noise of 700 ppm per LC sample is added. The simulated datasets are then processed with the same pipeline used as was used in this work for the K2 data.

The main conclusions from these simulations are: (i) The formal uncertainties in the "measured" transit depths (∼0.05%) are indeed underestimates of the actual uncertainties. The simulations using a constant depth have recovered depths with a scatter of 0.1%. This is likely the result of a combination of systematic effects induced by the short duration of the transits and the white noise terms. The simulations show that an actual scatter in the depths greater than 0.15% is easily recoverable. The scatter of our K2 depths is 0.2%, so we conclude that it is real. The ground based observations confirm this. (ii) The formal uncertainties also underestimate the true uncertainties of the transit times (defined as the scatter of measured times after removing a best linear trend). This effect is worse in the presence of starspots, but it does not depend on the transit depth variations. The effect can be large enough to explain all the scatter observed in the K2 timing analysis and therefore this scatter is not significant. We have added a systematic uncertainty in quadrature to the timings of 3 minutes, chosen to provide a best-fit reduced χ² of 1. (iii) Our treatment of stellar spots does not induce detectable depth variations (see Kawahara et al. 2013; Croll et al. 2015), as expected, mostly because we do not have the precision to detect them.

A fold of the activity corrected data (see Section 2.3) about the period we determined of 9.145704 hr is presented in Figure 6. The overall crudely triangular shaped transit profile is the result of a convolution of the intrinsic shape and the LC sampling time (see also Sections 4 and 7). The depth of the folded profile is 0.6%, and of course, this represents an average of the highly variable depths. The red curve is the mean out-of-transit normalized flux (averaged for orbital phases between $-0.5$ and $-0.25,$ and from 0.25 to 0.5). Note the clear positive "bump" in flux just after the transit egress and the smaller, but still marginally significant, bump just prior to ingress. These are significant at the 6-σ and 2-σ confidence limits, respectively.

Zoom In Zoom Out Reset image size

These "bumps," which are not normal features of exoplanet transits, will be important for understanding the basic nature of the transits (see Section 9). We call these features the "pre-ingress bump" and the "post-egress bump." Even in the cases of KIC 1255b and KOI 2700b, the other two exoplanets which appear to have dusty tails, the main distinguishing feature of their transit profiles is a post-transit depression (Rappaport et al. 2012, 2014) which is attributed to a trailing comet-like dust tail. In addition, KIC 1255b exhibits a pre-ingress "bump" which has been attributed to forward scattering in the dust tail near the head of the dust cloud (see, e.g., Brogi et al. 2012; Rappaport et al. 2012; Budaj 2013; van Werkhoven et al. 2014). By contrast, the transit profile of K2-22b has no post-egress depression, and the most prominent bump comes after the egress, rather than before. These features will be crucial to the interpretation of the dust "tail" in this system.

Finally, even though we believe the transit is due to a dust tail, as a baseline reference model we attempted to fit a standard transit profile of a solid planet over a limb-darkened star (Mandel & Agol 2002) to our folded light curve. The high distortion of the transit light curve due to the 30 minute sampling precludes obtaining precise transit parameters, but we were able to constrain the scaled semimajor axis from the transit itself to be $d/{R}_{*}={4.2}_{-0.5}^{+0.15},$ with the uncertainties estimated from an MCMC analysis. In fact, from the stellar properties obetained in Section 5.2, and Kepler's 3rd law, we can make a better direct estimate of d/R_* = 3.3 ± 0.2, which is compatible with the transit fit. We ran a final transit model with a Gaussian prior on d/R_* with a mean value of 3.3 and standard deviation of 0.2, based on the inferred stellar density. This, in turn, allowed us to estimate a mean K2 total transit duration (first to fourth contact) of 46 ± 1 minutes, an impact parameter of b = 0.68 ± 0.06 (see Table 4) and a mean depth of ${({R}_{{\rm{p}}}^{\prime }/{R}_{*})}^{2}=0.55\%$, where ${R}_{{\rm{p}}}^{\prime }$ can be understood as the mean effective radius of the dust grains.

We have also used the Vanderburg & Johnson (2014) flux time series data set for K2-22 to check against the results of our own pipeline. With the application of a simple high-pass filter, this independently processed light curve looks nearly identical to the one shown in the upper panel of Figure 1. A fold of the data about the period we have determined yields a transit profile that is essentially the same as that shown in Figure 6, including the appearance of a convincing post-egress bump, and a somewhat less significant pre-eclipse bump. To the extent that our pipeline and that of Vanderburg & Johnson (2014) are independent, this is a satisfying test that our processing has introduced no artifacts into the lightcurve.

An early inspection of the SDSS images of the target star, K2-22, showed that it is quite cool and likely an M star. The image also indicated the presence of a faint companion at a distance ∼2'' in the south–west direction.

We observed the target star and its companion with Hyper-Suprime Cam (HSC, Miyazaki et al. 2012) on the Subaru 8.2-m telescope on 2015 January 25 (UT). The sky condition on that night was clear and photometric. The HSC is equipped with 116 fully-depleted-type 2048 × 4096 CCDs with a pixel scale of 017. We took images of 3 and 60 s exposures through g-, i-, and z-band filters (see Figure 9). The host star is saturated in i- and z-band with 60 s exposures, while the companion star is not clearly seen in the g-band image with 3 s exposure. We thus use the 60 s exposure images for g-band and 3 s exposure images for i- and z-band for the analysis presented here.

Zoom In Zoom Out Reset image size

The presence of a stellar companion can be problematic, particularly when the distance between the stars is smaller than the pixel size of Kepler. We measure the flux ratios of the companion to the host star in the HSC g-, i-, and z-band images (see Figure 9) by using the DoPHOT program (Schechter et al. 1993), which performs PSF-fitting photometry while de-blending stars. The background star is situated at a distance of 191, with a position angle of 134° west from north (see Table 3 for the coordinates). We obtain flux ratios of 4.0 ± 0.5%, 8.6 ± 0.5%, and 11.7 ± 0.5% for the g, i, and z bands, respectively, where the uncertainties have been increased to take into account systematic sources of noise.

Table 3.  Properties of the Host Star and Companion

Parameter (units)	Host Star	Companion
R.A. (J2000)	11:17:55.856	11:17:55.763
Decl. (J2000)	+02:37:06.79	+02:37:05.48
u-mag (SDSS)	19.07 ± 0.05	21.68 ± 0.17
g-mag (SDSS)	16.44 ± 0.05	19.61 ± 0.05
r-mag (SDSS)	15.01 ± 0.05	18.79 ± 0.05
i-mag (SDSS)	14.38 ± 0.05	17.34 ± 0.05
z-mag (SDSS)	14.05 ± 0.05	16.34 ± 0.05
J (2MASS)	12.74 ± 0.05	14.87 ± 0.05
H (2MASS)	12.09 ± 0.05	14.27 ± 0.05
KS (2MASS)	11.91 ± 0.05	13.93 ± 0.05
Teff (K)	3830 ± 100	3290 ± 120
$\mathrm{log}\;g$	4.65 ± 0.12	⋯
[Fe/H]	0.03 ± 0.08	0.06 ± 0.20
M* (R⊙)	0.60 ± 0.07	0.27 ± 0.05
R* (M⊙)	0.57 ± 0.06	0.30 ± 0.08
L* (L⊙)	${0.063}_{-0.007}^{+0.008}$	${0.010}_{-0.005}^{+0.007}$
Spec Type	M0V ±1	M4V ±1
Distance (pc)	225 ± 50	225 ± 50

Note. The magnitudes and colors are taken from a combination of the SDSS and 2MASS photometry. The stellar parameters are inferred from the combined analyses of the Keck-HIRES, IRTF-SpeX, and UH88-SNIFS spectra. See Section 5 for details.

Download table as:  ASCIITypeset image

The contribution of the background star points toward a cooler and fainter companion. The bright star has been characterized as part of the K2-TESS catalog, in which effective temperatures are obtained from the colors of the stars (Stassun et al. 2014). The temperature reported there for the brighter target star is ∼3700 K, which means that both stars are likely to be dwarfs stars (see Section 5.2). The colors used in that analysis could have been blended due to the proximity of the companion. We also obtained the ugriz SDSS magnitudes and the J 2MASS magnitude for both stars, confirming that the colors used in Stassun et al. (2014) were correct. We re-derived the flux ratios in g, i and z band (5.4%, 6.5% and 12% respectively) confirming the measurements obtained using the Subaru images. The minimum uncertainties used in the fits are 0.05 magnitudes to take into account systematic effects. Adopting the full line of sight extinction of A_V = 0.17 based on the expected distance to the objects, and a $\mathrm{log}\;g$ of 5 due to their red colors and dwarf nature (see Section 5.2), we obtained temperatures of 3800 ± 150 K for the bright star and 3150 ± 200 K for its fainter companion.

The magnitudes and colors of the two stars based on a combination of the SDSS, Subaru, and 2MASS photometry are summarized in Table 3.

5.2.1. NOT-FIES Spectrum

5.2.2. Keck-HIRES Spectra

We also acquired five spectra of the host star with the Keck Telescope and HIRES spectrometer using the standard setup of the California Planet Search (CPS, Howard et al. 2010). Over the course of four nights, (2015 February 5–8) we observed under clear skies and average seeing (∼12). Exposure times of less than 5 minutes resulted in SNRs of 5–10 per pixel. Extensive scattered light and sky emission were unavoidable for exposures on 2015 February 5 due to the close proximity of the nearly full moon to the target star. All spectra were taken using the C2 decker, which is 087 wide by 14'' long, resulting in a spectral resolution of R = 60,000. The slit was oriented to minimize the amount of contamination from the companion star.

We estimated the effective temperature T_eff, surface gravity log g, iron abundance [Fe/H], and projected rotation velocity v sin i_⋆ of the host star from the co-added HIRES spectrum, which has a S/N of about 20 per pixel at 6000 Å. We used a modified version of the spectral analysis technique described in Gandolfi et al. (2008), which is based on the use of stellar templates to simultaneously derive spectral type, luminosity class, and interstellar reddening from flux-calibrated, low-resolution spectra. We modified the code to fit the co-added HIRES spectrum to a grid of templates of M dwarfs recorded with the SOPHIE spectrograph (Bouchy et al. 2008). We retrieved from the SOPHIE archive²⁸ the high-resolution (R = 75,000), high S/N (>80) spectra of about 50 bright red dwarfs encompassing the spectral range K5–M3 V. The stars were selected from the compilation of (Lépine et al. 2013). We downloaded spectra with no simultaneous thorium-argon observations, to avoid potential contamination from the calibration lamp.

The photospheric parameters of the template stars were homogeneously derived from the SOPHIE spectra using the procedure described in Maldonado et al. (2015), which relies on the ratios of pseudo-equivalent widths of different spectral features. Unfortunately, the low S/N prevented us from directly applying this technique to the co-added HIRES spectrum of the target.

Prior to the fitting procedure, the resolution of the template spectra was somewhat degraded to match that of the HIRES spectrograph (R = 60,000)—by convolving the SOPHIE spectra with a Gaussian function mimicking the difference between the two instrument profiles. A corrective radial velocity shift was estimated by cross-correlating the observed and template spectra. We restricted the spectral range over which the fit is performed to 5500–6800 Å and masked out the regions containing telluric lines. We selected the five best fitting templates and adopted the weighted means of their spectroscopic parameters as the final estimates for the target star. We found that the target star has an effective temperature of T_eff = 3780 ± 90 K, surface gravity of log g = 4.65 ± 0.12 (log₁₀ cm s⁻²), and iron abundance of [Fe/H] = 0.05 ± 0.08 dex. We also set an upper limit of 1.5 kms⁻¹ on the projected rotation velocity v sin i_⋆ by fitting the profile of several clean and unblended metal lines to the PHOENIX model spectrum (Husser et al. 2013) with the same parameters as the target star. Figure 10 shows the co-added HIRES spectrum in the spectal region around the Hα line, along with the best fitting SOPHIE template.

Zoom In Zoom Out Reset image size

5.2.3. IRTF-SpeX Spectrum

5.2.4. UH88 SNIFS Spectra

Spectro-photometric observations for 3 complete transits were obtained with OSIRIS on the GTC (see also Section 4). The GTC instrument OSIRIS consists of two CCD detectors with a FOV of 78 × 78 and a plate scale of 0127 per pixel. For our observations, we used the 2 × 2 binning mode, a readout speed of 200 kHz with a gain of 0.95 e-/ADU and a readout noise of 4.5 e-. We used OSIRIS in its long-slit spectroscopic mode, selecting the grism R1000R which covers the spectral range of 520–1040 nm with a resolution of R = 1122 at 751 nm. The observations and results we present here were taken using a custom built slit of 12'' in width, with the target and a comparison star both located in the slit. The use of a wider slit has the advantage of reducing the possible systematic effects that can be introduced by light losses due to changes in seeing and/or imperfect telescope tracking (Murgas et al. 2014).

During the first transit observed with the GTC (January 29) we took as a reference a very close comparison star to the east of the target. However, for the two last transits (February 13 and 14) we took a reference star with a similar brightness to K2-22 and located at a distance of 27 from the target. The position angle of the reference star with respect to the target was $-79\buildrel{\circ}\over{.} 2.$ The two stars were positioned equidistantly from the optical axis, close to the center of CCD#1, while CCD#2 was turned off to avoid crosstalk.

The basic data reduction of the GTC transits was performed using standard procedures. The bias and flat field images were produced using the Image Reduction and Analysis Facility (IRAF²⁹ ) and were used to correct the images before the extraction of the spectra. The extraction and wavelength calibration were made using a PyRAF³⁰ script written for GTC@OSIRIS long-slit data. This script automated some of the steps to produce the spectra such as: extraction of each spectrum, extraction of the corresponding calibration arc, and wavelength calibration (using the HgAr, Xe, Ne lamps provided for the observations). All spectra were aligned to the first spectrum of the series to correct for possible shifts in the pixel/wavelength solution during the observations caused by flexures of the instrument. Several apertures were tested during the reduction process, and the one that delivered the best results in terms of low scatter (measured in rms) in the points outside of transit for the white light curve was selected. The results presented here were obtained using apertures of 28, 40, and 44 pixels in width for the three transit observations, respectively. Final spectra were not corrected for instrumental response nor were they flux calibrated. Figure 11 shows the extracted spectrum of K2-22.

Zoom In Zoom Out Reset image size

The universal time of data acquisition was obtained using the recorded headers of the spectra indicating the opening and closing time of the shutter in order to compute the time of mid exposure. We then used the code written by Eastman et al. (2010)³¹ to compute the BJD time using the mid-exposure time for each of the spectra to produce the light curves analyzed here.

We constructed light curves over several spectral ranges to search for color dependencies of the transit shape and depth. In the case of the white light curve, the flux was integrated over almost the entire wavelength range of the observed spectra, between 570 and 900 nm, but avoiding the blue and red ends of the spectra where the SNR is lower. This step is particularly important, because observations at redder wavelengths suffer from fringing³² whereas using the data at bluer wavelengths not only did not improve the quality of the light curves, but increased the occurrence of outliers. Three narrower color light curves were created from the raw spectra for each star, by integrating all counts over the spectral ranges 570–680 nm, 680–790 nm, and 790–900 nm (see Figure 11). The final light curves used in the analyses were created by dividing the K2-22 light curves by the comparison star light curves. The scatter (std) for the relative white light curves can be found in Table 1.

The white light transit profiles of the three separate GTC transits were shown earlier in Figure 7. It is obvious that the depth of the transits change from one night to the next, confirming the rapidly varying nature of this object. Note that between transits 2 and 3, only 27.44 hr had passed and the transit depth nearly doubled.

To further investigate the nature of these changes in the GTC transit curves, in Figure 12 we plot the three color light curves for each transit separately. Transit 3 (bottom panel) is the deepest and shows a clearly increasing transit depth towards the blue. The SNR is lower in the other two sets of color light curves and the transits are shallower, and thus there is no clear trend in transit depth with wavelength. In order to interpret quantitatively the color dependence of the depths in the third GTC transit, we fit our standard three-segment transit profile model to the data from each of the nights, with three different transit depths to evaluate how the transit depth changes with color. We also fit each transit with five additional parameters, to construct a polynomial that reaches second order in time from mid-transit and it linearly depends on the mean-subtracted airmass and FWHM of the images of the host star. The total number of parameters is 21 for each observation, with a total 120, 132 and 114 data points in each of the three nights. After finding the best-fit model, we set the error of the flux measurements of each light curve to provide a best-fit standard χ² value equal to the number of degrees of freedom. We then took correlated noise into account computing the β factors for each light curve and each night (see Section 4), and multiplying the errors by their corresponding β factor. We use an MCMC routine to obtain posterior distributions for the transit depths, which are marginalized with respect to all the other 18 model parameters (see Table 5).

Zoom In Zoom Out Reset image size

From the marginalized posterior distributions we compute a quantity called the "Angström exponent," α, which is defined as $-d\mathrm{ln}\sigma /d\mathrm{ln}\lambda ,$ where σ is the effective extinction cross section. Here we take the transit depth (on each night) at wavelength, λ, to be proportional to the effective cross section at λ (under the assumption that the dust tail is optically thin). In order to use the transit depths, we first multiply each individual transit depth by (1 + D_i) to apply a dilution correction, where each D_i is the flux ratio between the faint companion and the host star evaluated at the flux weighted center of each selected band. These D_i values are obtained from our spectral energy distribution models (see Section 5.1), and have values of 4.4%, 6.6% and 8.9% with increasing wavelength, with uncertainties of 0.2%. We then find α = 0.83 ± 0.23, for the transit observed on February 15. For the other two transits, the values of α are 0.11 ± 0.37 and 0.13 ± 0.55, both slightly positive, but not statistically significantly different from zero. We made no corrections for the wavelength-dependent limb darkening properties, but utilized the quadratic limb darkening coefficients of Claret & Bloemen (2011) to estimate that the values of α would be lowered by between 0.15 and 0.02, for an equatorial transit (which is unlikely) and an impact parameter of b = 0.7, respectively. For higher impact parameters, up to 0.85 (see Table 4), the value of α could actually be raised by up to 0.15, and become even more significantly different from zero.

Table 4.  Properties of K2-22b

Parameter	Value
Orbital Perioda (days)	0.381078 ± 0.000001
Orbital Period (hr)	9.145872 ± 0.000024
Transit centera (BJD)	2456811.1208 ± 0.0006
${\dot{P}}_{\mathrm{orb}}/{P}_{\mathrm{orb}}$ (yr−1)a	≲3.5 × 10−7
Total transit durationa (minute)	48 ± 3
d/R*b	3.3 ± 0.2
Impact parameter, bc	0.68 ± 0.06
d (AU)	0.0088 ± 0.0008
Rp (R⊕)d	<2.5 ± 0.4
Mp (MJ)e	<1.4
θ*f	17° ± 15
${\dot{M}}_{\mathrm{dust}}$ (g s−1)g	≈2 × 1011
ℓ2 (leading tail only; units of R*)h	0.19–0.48
agrain (μm)h	0.3–1.0
Impact parameter, bh	0.42–0.78
ρ1/ρ2 (two-tail model, ℓ1 = ℓ2)h	<0.5

Notes.

^aDerived from the K2 and ground-based observations. ^bBased on the mass and radius of the host star given in Table 3 and Kepler's 3rd law. ^cDerived from the K2 observations using a standard Mandel & Agol (2002) fit to a hard-body transit. ^dBased on the shallowest K2 transits. ^e2-σ limit based on the Keck RVs. ^fEstimate of the half angle subtended by the star at the position of the planet. ^gFollowing the type of estimate made in Appendix D of Rappaport et al. (2014). ^h90% confidence limits based on the models of Section 9. The leading and trailing tails are assumed to have exponential scale lengths and maximum optical thicknesses, ℓ₂, ℓ₁, ρ₂, and ρ₁, respectively.

Download table as:  ASCIITypeset image

We interpret the value of the Angström exponents in terms of Mie scattering (for spherical dielectric dust particles) with a variety of different compositions. We computed the Angström exponent, α, over the range 630–840 nm for a set of power-law distributions for the dust particle sizes, with dN/da ∝ a^−Γ. In order to guarantee that the total cross section converges, we also need to specify a maximum grain size, a_max, as well as a minimum grain size which we fix at 0.01 μm. For the sake of specificity we adopted an illustrative dust composition of corundum, but the conclusions we draw are the same for a number of other common refractory materials, and in fact, for any material with real and imaginary indices of refraction of n ≃ 1.6 and 0.001 ≲ k ≲ 0.03. We plot the computed Angström exponent in Figure 13 as a function of a_max for five different power-law exponents, Γ. As can be seen from the figure, values of α in the range of ∼0–1, as indicated in Table 5, correspond to non-steep power-law indices of Γ ≃ 1–3 and maximum particle sizes of ∼0.4–0.7 μm. In turn, these correspond to "effective particle sizes" of 0.2–0.4 μm, where the effective radius is the average grain radius weighted by both the size distribution and the cross section, i.e.,

$$\begin{eqnarray}&&{a}_{\mathrm{eff}}={\displaystyle \int }_{0}^{{a}_{\mathrm{max}}}{a}^{1-{\rm{\Gamma }}}\;\sigma (a,\lambda )\;{da}\;/{\displaystyle \int }_{0}^{{a}_{\mathrm{max}}}{a}^{-{\rm{\Gamma }}}\;\sigma (a,\lambda )\;da\end{eqnarray} \tag{ 1 }$$

where σ(a, λ) is the wavelength dependent Mie extinction cross section for a particle of radius a.

Zoom In Zoom Out Reset image size

When taken together, all these observational results point clearly toward a planet in a 9-hr orbit that is disintegrating via the emission of dusty effluents. The lines of evidence pointing in this direction include: (1) erratically and highly variable transit depths (see Section 3); (2) transit profile shapes from ground-based observations that are likely variable (though with lesser confidence that the depth changes; see Sections 4 and 7); and (3) an average transit profile from the K2 data that exhibits clear evidence for a post-transit "bump" and also weaker evidence for a pre-transit "bump" (Section 3.2). The first of these is highly reminiscent of the disintegrating planet KIC 1255b (Rappaport et al. 2012), while the variable transit shapes are also detected in KIC 1255b from ground-based studies (R. Alonso et al. 2015, private communication; Bochinski et al. 2015). The transit profile in K2-22, with a post-transit "bump," is different from the transit profiles of KIC 1255b and KOI 2700b which show a post-transit depression as opposed to a post-transit bump. The first two of the above listed features point to obscuration by dusty effluents coming from a planet, while the third property needs to be explained in this same context. In the following sections we explore the significance and the interpretation of the transit profile.

A dust tail emanating from a planet, as inferred in the cases of KIC 1255b and KOI 2700b, trails the planet as is illustrated in Figure 6 of Rappaport et al. (2012; see also the middle panel of our Figure 15). Such dust tails are the way they would be seen in the reference frame of the planet, and note that the motion of the planet is implicitly in the opposite direction from the tail. In this case we would say that the tail "trails the planet." The reason for this is that the radiation pressure acting on the dust forces it into an eccentric orbit with its periastron located at the point where the particle was released. This orbit has a larger semimajor axis than that of the planet (see Appendix B of Rappaport et al. 2014). In turn, the larger orbit has a lower orbital frequency, and the particles appear to trail behind the planet in a comet-like tail.

In the case of a trailing dust tail, the ingress to the transit is sharp as the "comet head" moves onto the stellar disk (see Figure 6 in Rappaport et al. 2012). The trailing tail would lead to a depression upon egress as the tail slowly moves off of the stellar disk. In that case, the pre-transit "bump" is caused by forward scattering by the densest regions of the dust which have not yet reached the stellar disk. This is the case we believe we see in KIC 1255b (Brogi et al. 2012; Rappaport et al. 2012; Budaj 2013; van Werkhoven et al. 2014). If the tail of the planet is sufficiently short (i.e., compared to the radius of the host star) then there would be both a pre-transit and a post-transit "bump," the latter of which would dominate over the relatively shallow post-transit depression (see also Budaj 2013). In this case, the pre-transit bump would be somewhat larger due to the asymmetry in the direction of the tail.

A logical first guess as to how to produce a post-transit "bump" on the transit curve would be to reverse the direction of the comet-like tail. However, as we have seen, substantial radiation pressure forces inevitably lead to a trailing tail. Then, the question becomes how to produce a "leading dust tail" to the planet. One way to have dusty material lead the planet, i.e., moving faster, would be to have it overflow its Roche lobe (or, Hill sphere of influence) and fall in toward the host star. At first consideration this does not seem to work since a rocky planet in a 9-hr orbit will not be close to filling its Roche lobe. Rappaport et al. (2013) showed that the critical density for Roche-lobe overflow is largely a function of its orbital period, and is nearly independent of the properties of the host star. In particular, Equation (5) in Rappaport et al. (2013) suggests that critical density for a planet to be filling its Roche lobe is

$$\begin{eqnarray}&&{\rho }_{\mathrm{crit}}\simeq {\left(\displaystyle \frac{11.3\;\mathrm{hr}}{{P}_{\mathrm{orb}}}\right)}^{2}\;{{\rm{g}}\;\mathrm{cm}}^{-3}\simeq 1.5\;{{\rm{g}}\;\mathrm{cm}}^{-3}\end{eqnarray} \tag{ 2 }$$

where the right-hand value is for a 9-hr planet. It seems very unlikely that a planet with this low a mean density would be disintegrating via dusty effluents. Turning the problem around, we can ask what fraction of the Roche-lobe radius would be occupied by a planet with a mean density in the range of 5–8 g/cc, which might be more appropriate for a dust emitter (see, e.g., Rappaport et al. 2013). It is then evident that for densities which are ∼3–5 times higher than their critical densities, their radii are not even factors of ∼2 times smaller than their Roche lobes.

Therefore, we come to the conclusion that planets with rocky compositions in a 9-hr orbit are underfilling their Roche lobes by only a factor of ∼2. This is more than sufficient to prevent the planet from directly overflowing its Roche lobe. However, the potential difference between the planet's surface and the Roche lobe is about half the potential difference to infinity. Thus, if the mechanism that drives off the dust or the heavy metal vapors that condense into dust, e.g., via a Parker-type wind (Rappaport et al. 2012; Perez-Becker & Chiang 2013) imparts the full escape velocity of the planet or more, then the material can certainly reach the surface of its Roche lobe.

If the material reaching the Roche surface feels a substantial radiation pressure, it will still be blown back into a comet-like tail. By "substantial" we mean that β, the ratio of radiation pressure forces to gravitational forces, exceeds a certain value such as β ≳ 0.05. For sufficiently small values of β, however, particles initially directed toward the host star (which seems reasonable since the gas or dust emission should commence on the heated hemisphere) will largely fall toward the host star until the grains are pushed into orbit by coriolis forces. These particles will form a "leading tail" in the sense that the motion will be in front of the planet.

The values of β for the particles are computed according to

$$\begin{eqnarray}&&\beta \simeq \displaystyle \frac{{L}_{*}\sigma }{4\pi {{cGM}}_{*}\mu }\equiv \displaystyle \frac{3\;{L}_{*}\;\langle \tilde{\sigma }(a,\lambda )\rangle }{16\pi {{cGM}}_{*}\rho a}\end{eqnarray} \tag{ 3 }$$

where L_* and M_* are the luminosity and mass of the host star, a is the particle radius, σ is the particle cross section at wavelength λ, $\langle \tilde{\sigma }\rangle$ is the dimensionless cross section in units of πa² averaged over the stellar spectrum (see, e.g., Kimura et al. 2002), and μ and ρ are the mass and material density of the dust grains. Since, for low-mass main-sequence stars the luminosity scales roughly as ${M}_{*}^{4},$ we find that for a fixed material in the grains, β scales as

$$\begin{eqnarray}&&\beta \propto \displaystyle \frac{{M}_{*}^{3}\;\langle \tilde{\sigma }(a,\lambda )\rangle }{a}.\end{eqnarray} \tag{ 4 }$$

We evaluate the dimensionless spectrum-averaged cross section with a Mie scattering code (Bohren & Huffman 1983) for different assumed material indices of refraction (see Croll et al. 2014). Plots of β(a) for several different assumed compositions are shown in Figure 14. They are also compared against representative β(a) curves for KIC 1255b and KOI 2700b. The latter two systems have β(a) curves that are essentially a factor of 2 higher than for K2-22, largely due to the lower luminosity per unit mass of the latter host star.

Zoom In Zoom Out Reset image size

Expression (4) raises the interesting prospect that for low mass stars, the value of β will be sufficiently small that outflowing particles passing through their Roche lobes will be immune to radiation forces and act much in the same way as Roche-lobe overflowing material.

We next carried out a large number of simulations of dust particles ejected from a planet and initially moving in various directions with different speeds. For most of our simulations, we took the direction of the particles to be uniformly distributed within a cone of 30° radius and centered in the direction of the host star. The velocities were somewhat arbitrarily taken to be one half the escape speed from the surface of the planet's Roche lobe to infinity in the absence of the host star, i.e., $\sqrt{{{GM}}_{{\rm{p}}}/2{R}_{L}}$ where M_p is the planet's mass and R_L is its Roche-lobe radius. This is the minimum that can be contemplated for a Parker-wind type outflow. For each dust grain, a radius, a, is chosen from an assumed particle size distribution via Monte Carlo methods. We considered a power-law differential size distribution for the dust grains with a slope of $-2,$ as illustrative, with the maximum and minimum sizes in this distribution taken to be 1 μm and 1/20 μm, respectively. The luminosity of the host star is taken to be that of a main-sequence M-K star of 0.6 M_⊙ (see Table 3). For purposes of the numerical calculations of the particle dynamics only, we used a simple analytic fit to the plots shown in Figure 14 of the form:

$$\begin{eqnarray}&&\beta \propto \displaystyle \frac{\langle \tilde{\sigma }(a,\lambda )\rangle }{a}\propto \displaystyle \frac{{c}_{1}+{c}_{2}{a}^{3}}{1+{c}_{3}{a}^{4}}\end{eqnarray} \tag{ 5 }$$

where the c's are constants to be fit, and which depend on L_*, M_*, and the dust composition.

The results of our particle dynamics simulations are shown in Figure 15. The top panel represents the trajectories of all particles regardless of their size or corresponding value of β. The dust density is assumed to decay exponentially in time, with a time constant of 10⁴ s. Note that there is a forward moving "leading tail" in addition to a bifurcated trailing tail. The inner of the two trailing tails arises from particles with small values of β that are launched with a substantial velocity component in the forward direction. In the bottom panel of Figure 15 we see the two tails from the perspective of an observer viewing an equatorial transit; there are far more particles in the leading tail. For comparison, we show in the middle panel a case where the computed value of β was arbitrarily multiplied by a factor of 4; all other parameters remained the same. Note that the result is a purely trailing dust tail as one would have in a solar system comet. This latter exercise ensures that all the particles have a β that exceeds a critical value of ∼0.02, guaranteeing that the particles will go into a trailing tail. The bottom line is that dust-emitting planets around luminous stars should have predominantly trailing tails while the reverse is true for low-luminosity stars.

Zoom In Zoom Out Reset image size

In all of these calculations, we have ignored possible ram pressure forces on the dust grains due to a stellar wind from the host star. For a justification of why this may be a good approximation, see Appendix A of Rappaport et al. (2014).

9.3.1. Idealized Transit Models

The transit profiles expected for a planet with a trailing comet-like tail are already sufficiently complicated compared to a conventional hard-body planet transit. They include such issues as the attenuation caused by the (unknown) absorption profile, typically characterized by an exponential scale length³³ ; possible forward scattering which is influenced by the grain size and composition (which may be changing along the tail as different materials sublimate at different rates); convolution of the scattering phase function with the finite angular profile of the host star; and another convolution with the density profile of the tail. When adding the possibility of both a trailing and leading tail, as we have postulated in this work, the situation quickly becomes even more complicated. In an attempt to keep the numbers of free parameters to a minimum, we have adopted the following seven-parameter model: an exponential attenuation profile of the tail in both directions, each with its own scale length, ℓ; a relative density between the two tails, ρ₁/ρ₂; the size of the dust grains, a; and an overall normalization factor that yields the correct mean transit depth. In addition to these parameters associated with the dust tails, there is also an impact parameter for the transiting planet, and the time of orbital phase zero. We assume that the hard body of the planet itself does not make a significant contribution to the transit profile.

We illustrate in Figure 16 what some of the possible transit profiles might look like with just two of these parameters varying. First, we consider only a trailing dust tail, and we vary the exponential scale length (ℓ, in units of the host star's radius). These transit profiles are shown in the top panel of Figure 16. The tails' exponential scale length, ℓ, is color coded to help in distinguishing the different curves. No forward scattering is included in this case. Note how the shorter ℓ is, the smaller is the post-transit depression. The latter becomes quite noticeable when ℓ for the tail becomes greater than about 30% of the radius of the star. In the middle panel of Figure 16 the contribution from forward scattering by relatively large particles (∼ 0.5 μm) has been included. In all cases there is a noticeable pre-transit bump due to the forward scattering, but only for the case of short dust tails (i.e., small ℓ) is there a perceptible post-transit bump; it is typically suppressed by the post-transit depression as the dust tail is in the process of moving off the stellar disk. The grain sizes assumed here are sufficiently large so that 1.2λ/4a, the characteristic Mie scattering angle, is comparable with the angular size of the host star as seen from the planet. In that case, the effective angular scattering pattern from a single grain is somewhat larger than that of the angular size of the host star (i.e., ∼17°). In the bottom panel of Figure 16 the transit profiles are for the case where there is both a trailing and a leading dust tail. The exponential lengths of both dust tails have been taken to be the same (i.e., ℓ₁ = ℓ₂), but the leading tail is taken to have twice the optical thickness everywhere as the trailing tail. For large ℓ there are no pre- or post-transit bumps, as they are pulled below the discernible level by the corresponding depressions caused by the tails. By contrast, for the case of short tail lengths, i.e., ℓ ≲ 0.2 R_*, there is both a pre- and a post-transit bump, with the latter being slightly larger.

Zoom In Zoom Out Reset image size

9.3.2. Model Fit to the K2-22b Transit

Finally, we have attempted to fit a simple two-tail dust model to the observed transit profile of K2-22b. Once the model and its free parameters were selected, we computed the scattering pattern for the dust particles (taken to be of a single fixed size throughout both tails) via a Mie scattering code (Bohren & Huffman 1983). Then the scattering pattern was convolved numerically via a 2D integral with the radiation profile coming from the host star (including limb-darkening). In turn, this net effective scattering profile was convolved in 1D with the assumed exponential tails to produce the scattering pattern as a function of orbital phase. The attenuation of the beam from the host star was simply computed by placing a double-sided exponential profile in front of a limb-darkened stellar disk in different longitudinal locations (i.e., as a function of orbital phase) and at different vertical locations to get a best estimate of the impact parameter b (in units of the stellar radius). The sum of the absorption and forward scattering contributions was then added and convolved with the Kepler long cadence integration time.

We first explored a wide range of parameter space, especially focusing on the physically interesting parameters: b, the impact parameter, ℓ₁, ℓ₂, ρ₁/ρ₂, the exponential scale lengths and ratios of optical thickness of the two tails, respectively, and a, the grain size. We also considered models with a power-law distribution of grain sizes, but these produced no better fits than using a single grain size. To simplify the search, we fixed ℓ₁ = ℓ₂ under the assumption that the dust sublimation time and the speed away from the planet is similar in the two tails. However, we allowed ρ₁/ρ₂ to be a free parameter since the amount of dust in the two tails could be quite different. From this broad search, we were able to draw several interesting conclusions. (1) The best fits, by only a small margin, are found for the single (leading) tail models, i.e., no trailing tail is needed. (2) For single-leading tail models, the impact parameter is constrained to be 0.42 < b < 0.78 and 0.19 < ℓ₂ < 0.48 R_* (both 90% confidence limits). (3) For two-sided tail models the allowed ratio of optical thicknesses is constrained to be ρ₁/ρ₂ < 0.5 (assuming ℓ₁ = ℓ₂). (4) The best-fitting single grain-size models have 0.3 < a < 0.5 μm. The final parameter search and error estimation were done with an MCMC routine. These results are summarized in Table 4.

The overall best-fitting model is for b = 0.65 and ℓ₂ = 0.32 R_*. The results are shown in Figure 17. The black histogram is the same average transit profile as is shown in Figure 6. The red curve is the overall fit to the transit profile, including the convolution with the LC integration time. The overall model transit profile is comprised of two components: the direct attenuation of the flux from the host star (green curve) and the forward scattering from the dust tail back into the beam directed toward the observer (blue curve). Note how the forward scattering produces small bumps both before and after the transits, with the larger bump following the egress.

Zoom In Zoom Out Reset image size

The dust sublimation timescale depends on the equilibrium temperature (T_eq), the mineral composition, and the size of the dust grains (see, e.g., Kimura et al. 2002; Appendix C of Rappaport et al. 2014). T_eq for the dust grains in K2-22 is necessarily somewhat uncertain, ranging between ∼1500 and 2100 K, depending on exactly how it is computed (to within unknown factors of order unity). However, if we take T_eq = 1700 K as a typical grain equilibrium temperature, then the lifetime of a 1 μm grain composed of Fe, SiO, or fayalite is a matter of a few tens of seconds. On the other hand, minerals such as enstatite, forsterite, quartz, and corundum, would have sublimation lifetimes of 1/2, 3, 15, and 30 hr, respectively. SiC and graphite could last for a very much longer time than any of these. (All the dust parameters for the above estimates, including the heat of sublimation, are taken from the conveniently tabulated values of van Lieshout et al. 2014). Thus, there is no one "expected" dust sublimation timescale that we can anticipate. If on the other hand, we believe that the dust tail is only a few degrees of the orbit in length (see Table 4), and that β ≈ 0.05, then the "inferred" lifetime is ∼1 hr (see Equations (4) and (6) of Rappaport et al. 2014 for an explanation), and therefore the dust composition might resemble enstatite or forserite. However, one should not draw too many conclusions from this since the actual tail length depends on a combination of numerous (uncertain) parameters such as the particle sizes (and indeed the size distribution), the parameter β, the actual equilibrium temperature, and so forth. Nonetheless, we can say that minerals such as enstatite or forsterite do seem consistent with the observations.

We note that there is an interesting possibility for the dust streaming from the planet to interact with dust that has already been orbiting for a substantial while—provided that the sublimation lifetime will allow the dust to survive for a sufficiently long time. The time required for orbiting dust in a trailing dust tail (i.e., "outer-track orbits" with 0.05 ≳ β ≳ 0.02) to meet up again with the planet is approximately P_orb/(2 β) (see Equation (4) of Rappaport et al. 2014); for K2-22 this amounts to 10–25 P_orb (or ∼ 4–10 days). For particles with β ≲ 0.02 the orbits likely lie inside the planet's orbit ("inner-track orbits") with radii between ∼94% and 98% of d. This leads to synodic orbital periods for those dust particles of between ∼10 and 30 P_orb, fairly similar to the outer-track orbits. Thus, any orbiting dust that might collide with newly emitted dust must last for 4–10 days. The dust grains on the inner-track orbits have speeds in the rest frame of the planet that range from 0 to 6 km s⁻¹ (prograde), while the corresponding range of outer-track orbits is 0–25 km s⁻¹ (retrograde). Thus, the inner track particles would catch up to, and collide with, newly emitted dust grains with only a few km s⁻¹ relative velocity. By contrast, the outer track particles could collide with newly emitted dust at relatively high speeds of more than 10 km s⁻¹. At such speeds, a collision would deposit a mean energy of ∼1 eV per atomic mass unit within the dust grain. This seems likely to completely destroy most dust particles that happen to collide on the outer track.

Additionally, any orbiting dust grains which are moving quickly may catch up, and collide, with other more slowly orbiting dust. Again, the relative speeds of any colliding dust particles on the inner track would likely be of order a km s⁻¹ and would not necessarily destroy the dust, whereas the reverse is true for the outer-track dust grains.

We can make a crude estimate of the collision frequency of a dust grain that is orbiting the host star. Take the mean grain number density along the observer's line of sight during a typical transit to be n₀, and the collision cross section to be σ. In that case, the collision mean free path is ℓ = 1/(n₀ σ). If we take the radial thickness of the dust tail that is causing the transits to be Δr, then Δr ≃ 1/(n₀ σ) if both the optical cross section for extinction and the collision cross section are approximately equal (i.e., the geometric cross section appropriate for larger particles), and if the optical depth of the dust cloud during transits is of order unity. (The latter is necessary since the dust tail is thin in the vertical direction and covers only a small fraction of the disk of the host star (see middle panel of Figure 15). Using the above expression, we can estimate the mean collision time, τ_coll, for particles orbiting in the dust disk

$$\begin{eqnarray}&&{\tau }_{\mathrm{coll}}\approx \displaystyle \frac{1}{2\pi }\displaystyle \frac{{\rm{\Delta }}r}{d}\displaystyle \frac{{v}_{\mathrm{orb}}}{{v}_{\mathrm{rel}}}\displaystyle \frac{{n}_{0}}{\langle n\rangle }\;{P}_{\mathrm{orb}}\end{eqnarray} \tag{ 6 }$$

where v_orb and v_rel are the mean orbital speed in inertial space and the mean relative speeds of the particles in a dust disk, and $\langle n\rangle$ is the mean density in dust grains around the azimuth of the dust disk. These factors are extremely uncertain, but as an illustrative example, if we take Δr/d ≈ 0.05, v_orb/v_rel ≈ 10, and ${n}_{0}/\langle n\rangle \approx 100,$ then τ_coll is of order several times P_orb.

What is the fate of dust particles that collide but are not destroyed in the collision? There could be locations in the dust disk where the dust is slowed down by the collisions and the density builds up. In principle, these regions of enhanced density could be approximately fixed in the rotating frame of the planet and host star. However, since we see only a single transit-like feature during an orbital period, we conclude that such pile-up of dust is not significant anywhere, with the possible exception of near the planet where the dust density is the highest.

A careful examination of what the effects are of dust-dust collisions in this system or, for that matter, KIC 1255b and KOI 2700b are much beyond the scope of this paper, and we leave that study to another work.