Toward Detection of Exoplanetary Rings via Transit Photometry: Methodology and a Possible Candidate

Our simple model of a ringed planet adopted in this paper basically follows Ohta et al. (2009). The ring is circular, and has a constant optical depth τ everywhere between the inner and outer radii of R_in and R_out. We denote the radii of the star and planet by R_⋆ and R_p.

The configuration of the planet and ring during a transit is illustrated in Figure 1. The X-axis is approximately aligned with the projected orbit of the planet on the stellar disk, and the Z-axis is toward the observer. This completes the (X, Y, Z) coordinate frame centered at the origin of the ringed planet (left panel in Figure 1). The normal vector of the ring plane is characterized by the two angles θ and ϕ in a spherical coordinate (right panel in Figure 1).

Zoom In Zoom Out Reset image size

We also set up another coordinate system (x, y, z) centered at the origin of the star in such a way that the major and minor axes of the projected ring are defined to be parallel with x- and y-axes, respectively, with z-axis being toward the observer.

The ring is assumed to move along the planetary orbit with constant obliquity angles (θ, ϕ), and the planet is assumed to move on a Keplerian orbit around the star. The left panel in Figure 1 illustrates the transit of the ringed planet, whose impact parameter is b.

We assume a thin uniform ring with a constant optical depth τ for the light from the direction normal to the ring plane. Thus the fraction of the background stellar light transmitted through the inclined ring is given by $\exp (-\tau {(\sin \theta \cos \phi )}^{-1})$, and we define the shading parameter T as $1-\exp (-\tau {(\sin \theta \cos \phi )}^{-1})$. In our simple ring model, the value of T, instead of τ, fully specifies the effective optical transparency of the ring.

In summary, our simple ring model is characterized by five parameters; four (R_in, R_out, θ, ϕ) specify the geometry of the ring, the other is a shading parameter T. Instead of R_in and R_out, we use dimensionless parameters in fitting

$$\begin{eqnarray}&&{r}_{\mathrm{in}/{\rm{p}}}\equiv \displaystyle \frac{{R}_{\mathrm{in}}}{{R}_{{\rm{p}}}},\qquad {r}_{\mathrm{out}/\mathrm{in}}\equiv \displaystyle \frac{{R}_{\mathrm{out}}}{{R}_{\mathrm{in}}}.\end{eqnarray} \tag{ 1 }$$

The stellar intensity profile I(x, y) under the assumption of the quadratic limb darkening law is expressed in terms of two parameters u₁ and u₂:

$$\begin{eqnarray}\displaystyle \frac{I(x,y)}{{I}_{0}} & = & [1-{u}_{1}(1-\mu )-{u}_{2}{(1-\mu )}^{2}]\\ & & \times \,\left(\mu \equiv \sqrt{1-\displaystyle \frac{{x}^{2}+{y}^{2}}{{R}_{\star }^{2}}}\right),\end{eqnarray} \tag{ 2 }$$

where I₀ is the intensity at the center of the star. The physical conditions on the profile require the following complex constraints on u₁ and u₂:

$$\begin{eqnarray}&&{u}_{1}+{u}_{2}\lt 1,\ {u}_{1}\gt 0,\ {u}_{1}+2{u}_{2}\gt 0.\end{eqnarray} \tag{ 3 }$$

In this paper, we adopt ${q}_{1}={({u}_{1}+{u}_{2})}^{2}$ and ${q}_{2}={u}_{1}/(2({u}_{1}\ +{u}_{2}))$ instead of (u₁, u₂) following Kipping (2013). Then, Equations (2) and (3) are rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{I(x,y)}{{I}_{0}}=[1-2{q}_{2}\sqrt{{q}_{1}}(1-\mu )-\sqrt{{q}_{1}}(1-2{q}_{2}){(1-\mu )}^{2}],\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\mathrm{with}\ \ \ \ \ \ \ 0\lt {q}_{1}\lt 1,\ 0\lt {q}_{2}\lt 1.\end{eqnarray} \tag{ 5 }$$

In this parametrization, q₁ and q₂ vary independently between 0 and 1. This is useful in finding best-fit parameters (Kipping 2013). For reference, the Sun has q₁ = 0.49 and q₂ = 0.34 (u₁ = 0.47 and u₂ = 0.23) (Cox 2000).

Let D(x, y, t) be the blocked fraction of light coming from the location on the stellar disk (x, y). Due to the motion of the planet during a transit, D(x, y, t) is time-dependent and given as

$$\begin{eqnarray}&&D(x,y,t)\\ \quad =\,\left\{\begin{array}{ll}1 & :\,\mathrm{if}\ (x,y){\rm{}}\,{\rm{is}}\,{\rm{within}}\,{\rm{the}}\,{\rm{planetary}}\,{\rm{disk}}\\ T & :\,\mathrm{if}\ (x,y){\rm{}}\,{\rm{is}}\,{\rm{within}}\,{\rm{the}}\,{\rm{ring}}\,{\rm{disk,}}\,{\rm{but}}\,{\rm{out}}\,{\rm{of}}\,{\rm{the}}\\ & \quad {\rm{planetary}}\,{\rm{disk}}\\ 0 & :\,\mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 6 }$$

Then the normalized flux from the the system is given by

$$\begin{eqnarray}&&F(t)=1-\displaystyle \frac{{\int }_{\mathrm{stellar}\mathrm{disk}}\,I(x,y)D(x,y,t){dxdy}}{{I}_{\mathrm{all}}},\end{eqnarray} \tag{ 7 }$$

where the second term indicates the fraction of light blocked by a transiting ringed planet, and the total flux is

$$\begin{eqnarray}{I}_{\mathrm{all}} & = & {\displaystyle \int }_{\mathrm{stellar}\mathrm{disk}}I(x,y){dxdy}=\pi {I}_{0}{R}_{\star }^{2}\\ & & \times \,\left[1-\displaystyle \frac{2\sqrt{{q}_{1}}{q}_{2}}{3}-\displaystyle \frac{\sqrt{{q}_{1}}(1-2{q}_{2})}{6}\right].\end{eqnarray} \tag{ 8 }$$

We develop a reliable numerical integration method that solves the boundary lines of D(x, y, t) as described in Appendix A. Our method achieves the numerical error less than 10⁻⁷ in relative flux, and this is much smaller than a typical noise of the Kepler photometric data.

We briefly comment on three effects that we neglect in the analysis below; finite binning during exposure time, planetary precession, and forward scattering of the ring. While all of them are negligible for the Saturnian ringed planet with a long period, they may become important in other situations.

For the precise comparison of our light-curve predictions against the Kepler long cadence, we may have to take account of the finite exposure time (29.4 minutes) properly. In fact, the binning effect is shown to bias the transit parameter estimate in the case of short-period planets (Kipping 2010). For the long-period planets that we focus on here, however, the transit duration is sufficiently longer than the finite exposure time. Thus the binning effect is not important. In the case of the transit of Saturn in front of the Sun, for instance, the fractional difference of the relative flux is typically an order of 10⁻⁵ between models with and without the binning effect. This value is an order of magnitude smaller than the expected noise in the Kepler photometric data. Thus, we can safely neglect the binning effect in the present analysis.

The precession of a planetary spin would generate observable seasonal effects on the transit shape of a ringed planet (Carter & Winn 2010; Heising et al. 2015). Since our current target systems are extracted from those with a single transit, however, we can ignore the effect; the period of the precession is proportional to the square of the orbital period, and thus the precession effect during a transit is entirely negligible. Nevertheless, we note here that this could be an interesting probe of the dynamics of short-period ringed planetary systems that exhibit multiple transits.

In the present analysis, we consider the effect of light-blocking alone due to the ring during its transit. In reality, forward scattering (diffraction by the ring particles) may increase the flux of the background light. Let us consider light from the star to the observer through the ring particle with diameter d. First, light is emitted from the disk of the star, and arrives at the ring particles. The angular radius of the star viewed from the ring particles is about R_⋆/a, where a is the semimajor axis of the orbit, and R_⋆ is the stellar radius. Next, the light is diffracted by the ring particles, and the extent of the diffraction is described by the phase function (Barnes & Fortney 2004); the rough diffraction angle can be estimated from the first zero of the phase function θ ≃ 0.61λ/d, where λ is wavelength of light. In particular, the effect of the diffraction becomes significant when the viewing angle R_⋆/a is comparable to the diffraction angle. Let us define the critical particle size d_crit by equating R_⋆/a with θ = 0.61λ/d;

$$\begin{eqnarray}&&{d}_{\mathrm{crit}}=0.61\,\displaystyle \frac{a\lambda }{{R}_{\star }}=0.63\,\mathrm{mm}\left(\displaystyle \frac{a/{R}_{\star }}{2060}\right)\left(\displaystyle \frac{\lambda }{500\,\mathrm{nm}}\right).\end{eqnarray} \tag{ 9 }$$

Barnes & Fortney (2004) discussed the effect of diffraction using d_crit. When d ≥ 10d_crit, the diffraction angle is small, and light just behind the ring particles is diffracted to the observer. In this case, the diffraction does not affect the direction of light, and we may express the extinction due to absorption with a single parameter T.

When d ≤ d_crit/10, the diffraction angle is large, and the ring particles diffract light to wider directions. Then, the amount of light that reaches the observer significantly decreases, and we may model the extinction in terms of T.

In both cases, d ≥ 10d_crit and d ≤ d_crit/10, the extinction can be modeled with a single parameter T. In the case of Saturn with the typical particle diameter d = 1 cm, for instance, d_crit ≃ 0.63 mm from Equation (9) satisfies d > 10d_crit, so our model can be used to calculate the light curves of Saturn observed far from the solar system.

We should note that when the typical size of particles satisfies ${d}_{\mathrm{crit}}/10\leqslant d\leqslant 10{d}_{\mathrm{crit}}$, the forward scattering induces the rise in the light curve before the ingress and after the egress, and this effect can become the key to identifying the signatures of the rings out of other physical signals. Incorporating the diffraction into the model, however, requires intensive computation, and this is beyond the scope of this paper.

The Kepler mission monitored more than 150,000 stars over four years, and identified about 8000 planet candidates as Kepler Objects of Interest (KOIs). In this paper, we focus on long-period planet candidates because icy ring particles, as observed around Saturn, are supposed to survive only at locations far from the host star. Considering that the temperature of the snow line is 170 K (Hayashi 1981), we choose 37 KOIs whose equilibrium temperatures are less than 200 K. In addition, we selected planet candidates reported by recent transit surveys; 41 candidates from a search by Wang et al. (2015) and 28 candidates from Uehara et al. (2016). In Table 1, the numbers of planetary candidates in three groups are listed with the number of transits.

We exclude several systems, which are not suited for our search. For KOI-5574.01 in KOIs and KIC 2158850 in Wang et al. (2015), we cannot find the transit signal among the noisy light curves. For KOI-959.01 in KOIs with P = 10 days and KIC 8540376 in Wang et al. (2015) with P = 31.8 days, we cannot neglect the binning effect due to the short transit duration. After removing these systems, 89 planet candidates are left in total for our search. Tables 1 summarizes the number of targets, and Figure 3 shows the overlapped objects among KOIs (Wang et al. 2015; and Uehara et al. 2016).

Zoom In Zoom Out Reset image size

Inevitably the signature of a possible ring around a planet is very tiny. Long-period planet candidates exhibit a small number of transits (Table 1), and the precision of the transit light curves is not improved much by folding the multiple events. Therefore, the search for a possible ring signature crucially relies on the quality of the few transiting light curves for individual systems.

According to the automated procedures described in Appendices B and C, we classify the long-period planet candidates into the following four categories.

(A) Insufficient S/N to constrain ring parameters: since the anomalous feature due to the ring is very subtle, one cannot constrain the ring parameters at all if the intrinsic light-curve variation of the host is too large to be explained by any ring model. Thus we exclude those systems that exhibit a noisy light curve out-of-transit. The exclusion criteria depend on the adopted ring model to some extent, but are determined largely by the threshold S/N that we set as S/N = 10. For definiteness, we consider four different ring models (Table 2), and the details of the procedure are described in Appendices B and C.

(B) Sufficient S/N and no significant anomaly: a fraction of the systems has a sufficiently good S/N and exhibits no significant anomaly. In such a case, we can put physically meaningful constraints on the possible ring parameters (Section 4).

(C) Too large anomaly for a ringed planet: in contrast to (B), some systems exhibit a large anomaly in the transiting light curve that exceeds the prediction in the adopted ring models. Nevertheless, different ring models may be able to explain the anomaly, and we still continue to search for ringed planets in this category (Section 5).

(D) Reasonable anomaly for a ringed planet: finally, a small number of systems with good S/Ns indeed exhibits a possible signature that could be explained by the ring model. We perform an additional analysis to test the validity of the ring hypothesis in a more quantitative fashion (Sections 5 and 6).

The above classification is done on the basis of observed anomalies, which are derived by fitting a planet model to light curves. The data are taken from the Mikulski Archive for Space Telescopes, and we use the Simple Aperture Photometry (SAP) data taken in the long-cadence mode (29.4 minutes). In fitting, we use only the first transit in the light curve for each candidate in deriving the observed anomaly for simplicity. After fitting the planet model to data, the long-period planet candidates are automatically classified into the above categories (A) ∼ (D). Table 2 summarizes the results of classification for four models. In a later section, we use the classification according to model I, which contains more candidates in categories (B) ∼ (D) than the other three. In fact, the choice of model I is partly reasonable because the distant planets potentially have tilted rings like Saturn because of low tidal force.

Because candidates in (A) have insufficient S/N for further analysis, we do not consider them in the following analysis. In Section 4, we obtain upper limits on R_out/R_p for candidates in (B). In Section 5, we first search for the ringed planets in categories (C) and (D) by visual inspection, and later examine the reliability of transits more quantitatively. In Section 6, we interpret the possible ringed planet candidate.

Figures 5 and 6 show the light curves of candidates in categories (C) and (D), respectively. Candidates in (C), where the observed anomaly exceeds the prediction of model I, may be consistent with other ringed planets in different configurations. Thus, we search for ringed planets not only in (D) but also in (C).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We extract ringed planet candidates by visual inspection of their light curves on the basis of the following properties expected for ringed planets.

• 1.  
  Duration of ingress and/or egress is long.
• 2.  
  Transit shape is asymmetric due to the non-zero ϕ.

As a result, we identify five systems KOI-771(D), KOI-1032(C), KOI-1192(D), KOI-3145(D), and KIC 10403228(D) as tentative ringed planets. For the other four candidates in (D), which show no visible ring-like feature in the light curves, we obtain the upper limits on R_out/R_p in the same method as in the previous section (Table 3). In total, we obtain the upper limits on R_out/R_p for 12 candidates, and 6 of them have R_out/R_p ≤ 1.5.

For six candidates in (C) with no ring-like features, we cannot set the upper limits of ring parameters, and we conclude that the signals are not due to rings, but are due to the temporal stellar activities.

We examine the reliability of transit signals for the five preliminary candidates. As a result, we find that four are false positives, and KIC 10403228 still passes all criteria. More specifically, we regard a target as a false positive if one of the following criteria is satisfied (Coughlin et al. 2016).

• Criterion 1:  
  The target object exhibits a significant secondary eclipse, which is expected for an eclipsing binary.

  • Results:  
    None of our candidates exhibit the secondary eclipse.
  • Criterion 2:  
    The signal originates from the other nearby stars or instrumental noise.
  • Results:  
    Inspecting Target Pixel Files, we found that the dips in the light curves of KOI-1032.01, KOI-1192.01, and KOI-3145 do not come from the target stars. Figure 7 shows an example of KOI-1192.01. The Community Follow-up Observing Program (CFOP) classifies KOI-1032.01 as a false positive (Uehara et al. 2016). Wang et al. (2015) and Uehara et al. (2016) also indicate that KOI-1192.01 and KOI-3145 are false positivities. Moreover, we find that the transit depths in the light curves of KOI 771.01 differ in many pixels, and the contaminations from the non-target stars are very strong. Wang et al. (2015) also pointed out that this system is false positive. For KIC 10403228, the transit depths differ in only two pixels, while it is constant in the other pixels, so we conclude that the signal is originated from the target star. The more detailed discussion of KIC 10403228 is presented in a later section.
  • Criterion 3:  
    The transit simultaneously occurs at different stars in different pixels. This indicates that the signal does not originate from the target but from the instrumental noise.
  • Results:  
    The transit events of KOI-1032.01 and KOI-1192.01 are located at the same time. This result is consistent with that of Criterion 2.
  • Criterion 4:  
    The shape of the light curve is inconsistent with that of a transiting object.
  • Results:  
    From Figures 5 and 6, all signals fit well to transit-like features.

KIC 10403228 is the single system that passes all the criteria. Thus, we move on to the detailed pixel-based analysis next.

Zoom In Zoom Out Reset image size

KIC 10403228 is considered to be an M dwarf and has a nearby star separated by about 3 arcsec (Rappaport et al. 2014). According to the data taken by the United Kingdom Infrared Telescope (UKIRT), the nearby star is located at (R.A., decl.) = (19^h24^m5425, +47°32'575) and its J-band flux is about one-fifth of KIC 10403228. Here we examine the possibility that the transit is associated with this nearby star rather than KIC 10403228.

Figure 8 shows the light curve and fractional depth of the transit event in each of the pixels around KIC 10403228. The small transit depths in pixels A and B suggest that the source of the transit is not the nearby star shown by a red filled star, because otherwise the transit depths should be larger in those pixels close to the nearby star. To confirm this fact in a more quantitative way, we also calculate the centroid offset using the pixel-level light curves. As a result, we find that the flux centroid moves toward the nearby star during the transit and that the displacement is comparable to the value expected from the observed transit depth (5%) and the flux ratio in J-band (5:1). The variation of the transit depth and the centroid displacement consistently indicate that the transit is not due to the nearby star. While we may be able to evaluate the contamination on the light curve from this nearby star more quantitatively, it does not change our conclusion in any case, and we do not perform the detailed analysis for simplicity.

Zoom In Zoom Out Reset image size

We note that the transit signal contains a clear short-period modulation (panel B in Figure 8). Since the modulation is not visible at panel C, it is most likely due to the nearby star. Actually, there is another long-period modulation with P ≃ 35 days in the light curve, which may come from the target star. If these periods are related to the stellar spins, the nearby star is a fast rotating star, and the target star is a slow rotator. Thus, we may ignore the effect of gravity darkening of the target star.

We fit various models with and without the ring to the light curve of KIC 10403228 by minimizing the value of χ² defined in Equation (29). In practice, we use  a ±3.09 day time window to trim 300 data points centered around T₀ = 744.773 days (BKJD [=BJD−2454833 days]). To remove the long-term flux variations in the light curve, we adopt the model in Equation (30) that is composed of a fourth-order polynomial and the transit model F(t) in Equation (7). The standard deviation σ is estimated to be 9.17 × 10⁻⁴ from the out-of-transit data. This value is about 1.3 times larger than the error recorded in the SAP data.

As the transit of KIC 10403228 is observed just once, we cannot infer the orbital period from the timing of the transit. However, we can infer it from Kepler's law. The depth and V-shape of the observed transit imply that the transiting object is relatively large and grazing. Thus, we approximate the total transit duration T_tra as

$$\begin{eqnarray}&&{T}_{\mathrm{tra}}\simeq P\left(\displaystyle \frac{2{R}_{\star }}{2\pi a}\right)\left(\displaystyle \frac{\sqrt{1-{e}^{2}}}{1+e\sin \omega }\right),\end{eqnarray} \tag{ 10 }$$

where the last factor is a correction term due to an eccentricity e with ω being the argument of periapse. From Kepler's law and Equation (10), one obtains

$$\begin{eqnarray}&&P\simeq 450\ \mathrm{years}\left(\displaystyle \frac{{\rho }_{\star }}{12.6\ {\rm{g}}\,{\mathrm{cm}}^{-3}}\right){\left(\displaystyle \frac{{T}_{\mathrm{tra}}}{2\mathrm{days}}\right)}^{3}{\left(\displaystyle \frac{1+e\sin \omega }{\sqrt{1-{e}^{2}}}\right)}^{3}.\end{eqnarray} \tag{ 11 }$$

We obtain P ≃ 450 years if we adopt e = 0 and T_tra = 2 days for the transit duration of KIC 10403228, and the stellar density ρ_⋆ = 12.6 ± 6.0 g cm⁻³ from Wang et al. (2015). The stellar density in Wang et al. (2015) is adopted from Dressing & Charbonneau (2013), who estimated the stellar properties by comparing the observed colors taken in 2MASS and SDSS with the Dartmouth model (Dotter et al. 2008).

Before fitting, we simply examine how often we expect to see a transit of a planet with P ≃ 450 years. Assuming that all the stars host planets with P ≃ 450 years, the expected number of transit detections is given by

$$\begin{eqnarray}&&{n}_{\mathrm{tra}}=0.045\left(\displaystyle \frac{{N}_{\mathrm{target}}}{{\rm{150,000}}}\right)\left(\displaystyle \frac{{t}_{\mathrm{obs},\mathrm{dur}}/P}{4\ \mathrm{years}/450\ \mathrm{years}}\right)\left(\displaystyle \frac{{R}_{\star }/a}{1/{\rm{25,000}}}\right),\end{eqnarray} \tag{ 12 }$$

where a/R_⋆ = 25,000 is the fiducial value estimated from Equation (10), t_obs,dur is an observational period, and N_target is the number of target stars. The adopted values of t_obs,dur and N_target are the typical values of Kepler. The frequency of planets with P = 450 years would be less than 1, so we may see n_tra as the optimistic upper limit of the expected value. This current value of n_tra = 0.045 is small, but not too unlikely. Apart from the tiny ring-like feature, the overall shape of the signal is clearly due to the transiting event, and it is very difficult to explain the feature from the stellar activities.

We would like to comment on the reliability of P ≃ 450 years. The key parameters are ρ_⋆ and the eccentricity in Equation (11). For example, if the system is a giant star rather than an M dwarf, the density and the period become smaller. In this sense, to specify the correct stellar density, we would need a follow-up observation. Moreover, the eccentricity can also change the estimated period in Equation (11). If e = 0.6, the period can be changed by the factor of (1/8.0)–8.0, and if e = 0.9, the factor of change is within (1/82.82)–82.82 (or 5 years < P < 34,000 years). Thus, the planet with a relatively short period and a large eccentricity can also explain the data. Although the period is uncertain, the different period does not change the fitting results, so we adopt P = 450 years for the fiducial value for the time being.

For fitting, we adopt P = 450 years, and q₁ and q₂ from the official catalog of Kepler. In summary, there are nine free parameters, t₀, R_p/R_⋆, b, a/R_⋆, and c_i (i = 0–4) for the model without the ring, and five additional parameters $\theta ,\phi ,{r}_{\mathrm{in}/{\rm{p}}},{r}_{\mathrm{out}/\mathrm{in}}$, and T for the model with ring. We set the initial values of c_i (i = 0–4) to those obtained from a polynomial curve fitting for the out-of-transit data.

First, we fit the planet-alone model to the data. The blue line in Figure 9 is the best-fit model without the ring. The best-fit parameters are listed in Table 4. The residuals from the fit clearly have some systematic features, and the planet-alone model fails to fully explain the light curve, in particular, around 745.8 days (BKJD) in Figure 9. Therefore, we attempt to interpret the data with the ringed planet model. After trying a lot of initial values for fitting, we finally find two solutions that give at least local minima of χ² in Equation (29). Figure 9 shows those two solutions in the red and green lines. The best-fit parameters are shown in Table 4. The geometrical configurations for both solutions are shown in Figure 10. Clearly, models with the ring significantly improve the fit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Best-fit Parameters for the Transit of KIC 10403228

	Single Planet	Ringed Planet (Solution 1)	Ringed Planet (Solution 2)
Fixed Parameters
P (years)	450	450	450
q1	0.6737	0.6737	0.6737
q2	0.0767	0.0767	0.0767
Variables
t0 (day)	0.065 ± 0.0021	0.017 ± 0.016	0.0386 ± 0.0052
Rp/R⋆	0.46 ± 9.09	0.27 ± 0.33	0.45 ± 0.05
b(=a cosi/R⋆)	1.16 ± 10.6	0.99 ± 0.44	1.14 ± 0.06
a/R⋆	24275.0 ± 31617.1	25394.5 ± 765.0	25743.9 ± 215.6
T	⋯	0.13 ± 0.07	1.0 (converged to upper bound)
θ(deg)	⋯	59.4 ± 45.0	12.3 ± 4.0
ϕ(deg)	⋯	52.3 ± 24.9	72.0 ± 4.5
Rin/Rp	⋯	1.003 ± 2.95	1.59 ± 0.42
Rout/Rin	⋯	2.89 ± 8.0	1.61 ± 0.59
c0	0.9972 ± 0.0002	1.00023 ± 0.00034	0.9999 ± 0.0003
c1	0.00066 ± 0.00008	0.0014 ± 0.0001	0.0015 ± 0.0001
c2	0.00067 ± 0.00011	−0.00047 ± 0.00015	−0.00030 ± 0.00013
c3	0.000158 ± 0.000013	0.000053 ± 0.000015	0.000038 ± 0.000015
c4	−0.000072 ± 0.000011	0.000021 ± 0.000014	0.000005 ± 0.000012
(Given R⋆ = 0.33 ± 0.05 R⊙)
Rp (RJa)	1.48b	0.88 ± 1.07	1.44 ± 0.27
${R}_{\mathrm{in}}\,({R}_{J}$ a)	⋯	0.89 ± 2.65	2.29 ± 0.72
${R}_{\mathrm{out}}\,({R}_{J}$ a)	⋯	2.56 ± 1.36	3.69 ± 0.23
(Statistical values)
${\chi }_{\mathrm{red}}^{2}$ (=χ2/dof)	2.70 (=784.0/290)	1.24 (=354.5/285)	1.23 (=349.1/286c)

Notes.

^aR_J is the radius of Jupiter. ^bR_p/R_⋆ = 0.46 and R_⋆ = 0.33 R_⊙ are assumed. ^cT is converged to the upper bound, so dof = 286 = 285 + 1.

Download table as:  ASCIITypeset image

In Table 4, values of R_p, R_in, and R_out are calculated on the assumption of R_⋆ = 0.33 ± 0.05 R_⊙ (Wang et al. 2015). It turns out that the resulting ratio of ring and planet radii is similar to that of Saturn: R_in ≃ 1.5R_p and R_out ≃ 2.0R_p.

We comment on the implication of the fitted model for KIC 10403228 in the following. The radiative equilibrium temperature of the ring particle is given by

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}\simeq 15.1\,{\rm{K}}{\left(\displaystyle \frac{{\rm{25,000}}}{a/{R}_{\star }}\right)}^{0.5}\left(\displaystyle \frac{{T}_{\star }}{3386\,{\rm{K}}}\right){\left(\displaystyle \frac{1-A}{1-0.5}\right)}^{0.25},\end{eqnarray} \tag{ 13 }$$

where we fiducially adopt the Bond albedo of the ring particle A of 0.5. The stellar effective temperature T_⋆ = 3386 K of KIC 10403228 is taken from Wang et al. (2015). Since the equilibrium temperature expected from the model is much lower than the temperature 170 K at the snow line (Hayashi 1981), icy particles around the planet can survive against the radiation of the host star.

The best-fit values of θ = 594 for solution 1 imply a significantly tilted ring with respect to the orbital plane, and θ = 123 for solution 2 implies a slightly tilted ring. We examine the stability of those tilted rings on the basis of a simple tidal theory. Under the assumption that the ring axis is aligned with the planetary spin, the damping timescale of the ring axis is equal to that for the orbital and equatorial planes of the planet to be coplanar. This timescale is given by a tidal theory (e.g., Santos et al. 2015):

$$\begin{eqnarray}{\tau }_{\mathrm{tidal}} & = & \displaystyle \frac{{P}_{\mathrm{orb}}Q}{9\pi {k}_{2}}\displaystyle \frac{{\rho }_{p}}{{\rho }_{s}}{\left(\displaystyle \frac{a}{{R}_{\star }}\right)}^{3}\simeq \displaystyle \frac{{{GP}}_{\mathrm{orb}}^{3}Q}{27{\pi }^{2}{k}_{2}}{\rho }_{p}\\ & = & 6.94\times {10}^{16}\,\mathrm{years}{\left(\displaystyle \frac{{P}_{\mathrm{orb}}}{450\mathrm{years}}\right)}^{3}\\ & & \times \,\left(\displaystyle \frac{2.3\times {10}^{-4}}{{k}_{2}/Q}\right)\left(\displaystyle \frac{{\rho }_{p}}{0.70\,{\rm{g}}\,{\mathrm{cm}}^{-3}}\right),\end{eqnarray} \tag{ 14 }$$

where P_orb is an orbital period, Q is a dissipation factor, and k₂ is the second Love number. If we adopt k₂/Q = 2.3 × 10⁻⁴ (Lainey et al. 2012) and ρ_p = 0.70 g cm⁻³ (Cox 2000) of Saturn, the damping timescale is sufficiently long. Thus, the best-fit configurations are consistent with the spin damping theory even under the assumption that the equatorial plane of the planet is coplanar with the ring plane. Thus, the tilted rings of our best fits also imply the non-vanishing obliquity of the planet.

The ringed planet model is consistent with the data. However, the V-shape of the transit (Figure 9) is also a typical feature of eclipsing binaries, and the estimated period of ∼450 years may be too long to be detected in four years of Kepler's observation (Equation (12)). Therefore, we discuss other scenarios without a planetary ring. For this purpose, in the following, we present two possible hypotheses, which can also explain the data; a binary with a circumstellar disk and a hierarchical triple.

In this section, we pursue the possibility that the current transit is caused by an eclipsing binary with a circumstellar disk rather than a planetary ring. Actually, the fitting result in the previous section is also applicable to this binary scenario, so we may compare the plausibilities of the eclipsing binary and planet scenarios to test the circumstellar disk model. For this specific purpose, we use the public code VESPA (Validation of Exoplanet Signals using a Probabilistic Algorithm) (Morton 2012, 2015). With VESPA, we compare the likelihoods of the following four scenarios: "HEBs (Hierarchical Eclipsing Binaries)," "EBs (Eclipsing Binaries)," "BEBs (Background Eclipsing Binaries)," and "Planets" (Transiting Planets), adopting a variety of different periods.

We adopt JHK-magnitudes from 2MASS (J-mag =13.429 ± 0.028, H-mag = 12.793 ± 0.03, and K-mag =12.518 ± 0.027), (R.A., decl.) = (19^h24^m54413, +47°32'575), maxrad = 3.0 arcsec (angular radius of the simulated region), Kepmag = 16.064, and R_p/R_⋆ = 0.3. In reality, those observed colors might be contaminated by the nearby star discussed in Section 5.3, but we assume that the contamination is sufficiently small in the present analysis. Given these inputs, VESPA calculates the star populations and the probability distribution of transit shaped parameters for the above four scenarios. For our adopted set of input parameters, VESPA identifies the primary star as an M dwarf consistent with the classification of Dressing & Charbonneau (2013). We repeat the simulation 10 times with different initial random numbers according to the prescription of VESPA.

Figure 11 shows the relative probability of each scenario for different assumed periods. We define the relative probability as the product of the "prior" and "likelihood" computed by VESPA, multiplied by 1000 days/P. The last factor of 1000 days/P corrects for the probability that a long-period transit is observed in a given observing duration much shorter than the orbital period, which is not taken into account in the "prior" of VESPA. The plot shows the medians and the standard deviations of the probabilities computed from 10 sets of simulations. While the binary scenarios are more likely than the Planets scenario for the shortest and the longest periods investigated here, the Planets scenario is the most preferred in the intermediate region (10 years ≲ P ≲ 100 years). The result suggests that the planetary interpretation of the light curve is not so unlikely compared to the binary scenario, though there is a fair amount of probability that this is a false positive. Another important implication of Figure 11 is that the likelihood of orbital periods in the Planets scenario is much broader than what we intuitively thought before, and not sharply peaked around 450 years.

Zoom In Zoom Out Reset image size

While Figure 11 represents our final result from VESPA, we point out two additional factors that may be of importance for more detailed arguments.

First, the period distribution and the overall fraction of long-period planets and binaries have not been taken into account. Occurrence rate of giant planets around M dwarfs is given by Clanton & Gaudi (2016). They estimated the frequency of the planets with 10² M_⊕ < M_p < 10³ M_⊕ to be ${0.039}_{-0.025}^{+0.042}$ for 10³ days < P < 10⁴ days and ${0.013}_{-0.010}^{+0.025}$ for 10⁴ days < P < 10⁵ days. On the other hand, Janson et al. (2012) estimated the multiplicity distribution of the binaries in 3–227 au and found the overall occurrence rate of 0.27 ± 0.03 peaked around 10 au. These results imply that planets around M dwarfs are rarer than their stellar companions by one or two orders of magnitude. This difference in the overall frequency may further increase the relative plausibility of the EBs scenario compared to the Planets scenario.

Second, what also matters in reality is the frequency of planetary rings and circumstellar disks that produce the observed anomaly in addition to the transit signal. It is, however, far beyond our current knowledge to estimate these factors rigorously.

Given these difficulties, follow-up spectroscopy or high-resolution imaging would be more feasible to distinguish the EBs and Planets scenarios.

An eclipse due to a close binary (rather than a single star/planet) on a wide orbit around the primary M star is yet another possibility to explain the asymmetric and long transit-like signal observed for KIC 10403228. This is because the orbital motion of the occulting binary can produce the acceleration that modifies the in-eclipse velocity of the occulting object(s) relative to the primary. To test this possibility, we consider a hierarchical triple system consisting of a short-period binary ("inner" binary) orbiting around and eclipsing the primary M star on a wide orbit ("outer" binary). In the following, we only take into account the luminosity of the occulted star and ignore the flux from the smaller binary. We also assume that the orbits of both inner and outer binaries are Keplerian, and use the subscripts "in" and "out" to denote their parameters. A mass ratio of the inner binary is fixed to 1 for simplicity. In this model, the motion of the two components of the inner binary is specified by t_0,in (inferior conjunction of the inner binary), P_in, a_out/a_in, i_in, Ω_in (longitude of ascending node relative to that of the outer binary), in addition to the parameters for a single-planet model (now with the subscripts "out"). We fit all of these parameters, except for the stellar density ρ_⋆ = 13.0 gcm⁻³, q₁ = 0.6737, q₂ = 0.767, the time offset T₀ = 744.773 days in Equation (30), and the same baseline as obtained in Table 4 (solution 2), which are fixed. The mass ratio of the outer binary is related to P_in, P_out, and a_out/a_in as

$$\begin{eqnarray}&&q\equiv \displaystyle \frac{{M}_{\mathrm{in}}}{{M}_{\mathrm{out}}}=\displaystyle \frac{1}{{({a}_{\mathrm{out}}/{a}_{\mathrm{in}})}^{3}{({P}_{\mathrm{in}}/{P}_{\mathrm{out}})}^{2}-1},\end{eqnarray} \tag{ 15 }$$

where M_in is a total mass of the inner binary, and M_out is the mass of the primary star.

Figure 12 shows one of the best-fitting models with P_out = 1396.615 days, ${a}_{\mathrm{out}}/{a}_{\mathrm{in}}=54.53$, Ω_in = −0.00945, P_out = 10.96 days, b_out = 2.024, $\cos {i}_{2}=0.103$, t_0,out =1.965 × 10⁻³ days, and t_0,in = 1.965 × 10⁻⁵ days. In this solution, we find q = 0.0956, which leads to M_in = 30M_J. We also obtain χ²/dof = 379.3/292, which is comparable to the ringed-planet model. In this solution, the observed ∼2 day duration is reproduced despite the fact that the value of P_out is much shorter than required for the planetary scenario. This is made possible because the orbital motion of the inner binary cancels out the high orbital velocity of the outer binary. In addition to this particular solution, we find various other solutions with similar χ² values for a wide range of P_out. In general, the solutions with longer P_out are found to correspond to smaller q; for example, we find q ≃ 0.02 for P_out ≃ 30 years, and q ≃ 0.003 for P_out = 300 years. For M_⋆ = 0.3 M_⊙, these mass ratios q = (0.1, 0.02, 0.003) translate into M_in = (30M_J, 6M_J, M_J). Thus, in this scenario, the system can be composed of three low-mass stars or a star with a binary planet.

Zoom In Zoom Out Reset image size

The advantage of this scenario is that the observed long transit can be reproduced with much smaller P_out than in the ringed-planet model, which leads to far higher transit/eclipse probability. On the other hand, it is also true that the parameters need to be finely tuned to cancel the two orbital motions. While the degree of required fine tuning is crucial in comparing the evidence of this hypothesis with the planetary or stellar ring models, the evaluation of this factor is not trivial given the large parameter space. In addition, there still remain uncertainties in frequency of the hypothetical hierarchical triple (three low-mass stars or a star with a binary planet). Given these complexities, it is difficult to conclude whether or not this scenario is favored compared to the above two. Again, the follow-up observation will be effective for the further study.

So far, we have presented the three leading scenarios: the "planetary ring scenario," the "circumstellar disk scenario," and the "hierarchical triple scenario." However, there still remain other possibilities that may potentially account for the light curve of KIC 10403228. In this section, we examine these possibilities and show that they are unlikely to explain the data. Throughout this section, we basically assume that the transit is caused by a planet, but the results in this section are also applicable for the stellar eclipse.

6.4.1. Oblate Planet

6.4.2. Additional Transit Due to an Exomoon

In Section 6.3, we only consider the additional motion of an occulting object due to an accompanying object. However, a transit of the accompanying object (e.g., exomoon) itself is yet another possibility for the peculiar light curve of KIC 10043228. As shown below, this possibility is ruled out by the shape of the anomaly.

As shown in Figure 9, the anomaly in the light curve is significant only in the latter half. Motivated by this fact, we fit the light curve using the planet-alone model, masking the latter half of the transit and adopting the same baseline as obtained in Table 4 (solution 2); the difference between this model and the observed light curve would represent the anomalous contributions from anything other than the main transiting planet. The result in Figure 13 clearly shows that the anomaly consists of a short rise in the flux followed by a more significant dip. Such a feature is clearly inconsistent with the transit of an exomoon.

Zoom In Zoom Out Reset image size

6.4.3. Anomalies Specific to In-transit Data

6.4.4. Stellar Noise

The ring-like structure in the light curve shows up only for a short duration. Thus, the short-term stellar noise might mimic the ring-like anomaly just by chance. To discuss this possibility, we investigate the statistical property of the stellar activity of KIC 10403228. Specifically, we consider how frequently one encounters stellar noises comparable to the anomalous in-transit residuals. As will be shown, we find it difficult to reproduce the feature with stellar activities of KIC 10403228. In principle, we could check to see if the similar feature arises in stars other than KIC 10403228 more generally, but it is a separate question and does not answer if the signal for the particular star is due to that stellar activity. Therefore, we analyze the light curve of KIC 10403228 alone in this section.

To focus on the short-term noises, we remove the long-term variations by dividing the light curves into short segments and fitting each of them with polynomials. The more specific procedure is as follows. We exclude in-transit data as well as data around gaps in the light curve. From the remaining data, we pick up a segment of 6.18 day long light curve centered around a randomly chosen time and fit it with a quartic polynomial to remove the variation within the segment. In principle, one could use different functions (e.g., a spline function) or a different time window for detrending, but in any case the final results are insensitive to these choices. For consistency, we adopt the same baseline and time window as those used in Section 6.1.

We iterate "picking up a segment" and "detrending" procedures 1000 times and obtain 1000 segments of detrended light curves, whose centers are randomly distributed over the whole observing duration. We note that the total number of points in the detrended segments is 1000 × 300 = 3.0 × 10⁵, which is sufficiently large to sample all of the original data points (N = 10,000). By averaging the 1000 detrended light curves at each time, we obtain one light curve. This averaging operation suppresses the dependence on the choice of the central time of each segment. Figure 14 shows the resulting detrended light curve (bottom) along with the light curve before detrending (top).

Zoom In Zoom Out Reset image size

Now we move on to the comparison of the statistical property of stellar activities and the residuals of fit in Figure 9. Let us define F_data(t) as the flux ratio of the detrended light curve with respect to the mean. To investigate the short-term correlation of stellar activities, we divide the light curves into continuously brightening events (F_data(t) > 1) and fading events (F_data(t) < 1). Then, we compute the duration and amplitude (average of the deviation from the mean $| {F}_{\mathrm{data}}(t)-1|$) for each event. For comparison, we also calculate the duration and average relative flux for events in residuals in Figure 9. The left panel in Figure 15 is the scatter plot of the duration and average relative flux of events for three groups.

• (a)  
  All events out of the transit (the black data in Figure 14).
• (b)  
  Residuals of the ringed-planet fit (the red line in Figure 9).
• (c)  
  Residuals of the single-planet fit (the blue line in Figure 9).

The right panel in Figure 15 shows the distributions of durations for the three groups. In each duration bin, the vertical axis shows the total number of points in all events with that duration. The distribution of (a) is normalized to give the same total number of events as (b) and (c). The quoted error-bars are simply computed from Poisson statistics of the number of each event. Figure 15 shows that the distribution (b) is closer to (a) than (c). Thus, the ringed planet model is better than the planet model in terms of property of the correlated noise.

Zoom In Zoom Out Reset image size

So far, we have shown that the ring-like anomaly cannot be explained statistically. We further consider whether the stellar noise can mimic the light-curve shape itself. We examine this hypothesis by focusing on the most significant fading event in the out-of-transit data; see the left panel of Figure 15. The light curve of this event is shown in Figure 16. We would like to see if the combination of the planet model and this event can reproduce the ringed-planet-like feature. To do this, we appropriately embed the transit of the planet into the light curve around the fading event. Here, the parameters of the planet are the same as in Table 4. Then we fit the two models with and without a ring to those data, as shown in Figure 16(b). As a result, we find the difference in χ² of the two models to be 157.9, which is smaller than 434.9 obtained in Section 6.1 for solution 2. Thus, we conclude that it is difficult to reproduce the ring candidate by combining the stellar activities and the transit of the planet.

Zoom In Zoom Out Reset image size

6.4.5. Combination of the Above Mechanisms

Conditions for possible values of r_i are divided into the following three cases.

• (a)  
  Intersections of the edge of the planet (circle) and the edge of the ring (ellipse).
• (b)  
  Extreme points of the distance function from the center of the star to the edge of the planet (circle).
• (c)  
  Extreme points of the distance function from the center of the star to the edge of the ring (ellipse).

The number of r_i is at most eight for (a), two for (b), and two for (c). (a) and (b) are reduced to quadratic equations, which can be easily solved. The last case can be reduced to quartic equations. Here, we derive the quartic equations using the method of the Lagrange multiplier. Let the length of the major axis be 2R and that of the minor axis be $2R(1-f)$, where f is the oblateness. We set the center of the ellipse to be at (x_p, y_p). For (x, y) on the edge of the ellipse, we define the following function:

$$\begin{eqnarray}A(x,y,\lambda ) & = & {x}^{2}+{y}^{2}+\lambda \left[{\left(\displaystyle \frac{x-{x}_{{\rm{p}}}}{R}\right)}^{2}\right.\\ & & +\left.{\left(\displaystyle \frac{y-{y}_{{\rm{p}}}}{R(1-f)}\right)}^{2}-1\right]\,.\,f\end{eqnarray} \tag{ 17 }$$

From the condition, we need

$$\begin{eqnarray}&&\displaystyle \frac{\partial A}{\partial x}=2x+\displaystyle \frac{2(x-{x}_{p})\lambda }{{R}^{2}}=0\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial A}{\partial y}=2y+\displaystyle \frac{2(y-{y}_{p})\lambda }{{(1-f)}^{2}{R}^{2}}=0\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial A}{\partial \lambda }={\left(\displaystyle \frac{x}{R}\right)}^{2}+{\left(\displaystyle \frac{y}{R-{Rf}}\right)}^{2}-1=0\end{eqnarray} \tag{ 20 }$$

We reduce the above three equations to the following:

$$\begin{eqnarray}&&\displaystyle \frac{{\lambda }^{4}}{{(1-f)}^{4}{R}^{8}}+\displaystyle \frac{2{\lambda }^{3}}{{(1-f)}^{2}{R}^{4}}\left[\displaystyle \frac{1}{{R}^{2}}+\displaystyle \frac{1}{{(1-f)}^{2}{R}^{2}}\right]\\ &&\quad +{\lambda }^{2}\left[\displaystyle \frac{1}{{R}^{4}}+\displaystyle \frac{4}{{(1-f)}^{2}{R}^{4}}+\displaystyle \frac{1}{{(1-f)}^{4}{R}^{4}}\right.\\ &&\quad -\left.\displaystyle \frac{{x}_{{\rm{p}}}^{2}}{{(1-f)}^{4}{R}^{6}}-\displaystyle \frac{{y}_{{\rm{p}}}^{2}}{{(1-f)}^{2}{R}^{6}}\right]\\ &&\quad +\lambda \left[\displaystyle \frac{2}{{R}^{2}}+\displaystyle \frac{2}{{(1-f)}^{2}{R}^{2}}-\displaystyle \frac{{x}_{{\rm{p}}}^{2}+{y}_{{\rm{p}}}^{2}}{{(1-f)}^{2}{R}^{4}}\right]\\ &&\quad +1-\displaystyle \frac{{x}_{{\rm{p}}}^{2}}{{R}^{2}}-\displaystyle \frac{{y}_{{\rm{p}}}^{2}}{{(1-f)}^{2}{R}^{2}}=0.\end{eqnarray} \tag{ 21 }$$

In general, quartic equations are analytically solved, but we compute the solutions for the equation using a root-finding algorithm, because of the complexity of the analytic solution. x and y are calculated from derived λ as follows:

$$\begin{eqnarray}&&x=\displaystyle \frac{-{x}_{{\rm{p}}}}{1+(\lambda /{R}^{2})}+{x}_{{\rm{p}}},\ y=\displaystyle \frac{-{y}_{{\rm{p}}}}{1+(\lambda /{((1-f)R)}^{2})}+{y}_{{\rm{p}}}.\end{eqnarray} \tag{ 22 }$$

The number of solutions for (x, y) is at most four. We exclude the solutions including complex numbers and/or (x, y) not on the ellipse. Equation (22) gives the singular solutions when

$$\begin{eqnarray}&&\lambda =-{R}^{2},\ -{(1-f)}^{2}{R}^{2}.\end{eqnarray} \tag{ 23 }$$

Inserting the above values into Equations (18) or (19), we find x_p = 0 or y_p = 0. In this case, we cannot use Equation (22), but the conditions are reduced to the quadratic equations, which can be easily solved.

To derive θ_i,l(r), we calculate the intersections of a circle, centered at (0, 0), with the radius r, and a transiting object, composed of the circle (planet) and two ellipses (rings). The center of the ring system is (x_p, y_p). The intersections of two circles are easily computed and the number of the intersection points is two at most. Here, we derive the equations for intersection points of a circle and an ellipse. Let the radius of the circle be r. We select the same ellipse as before. For simplicity, we introduce the following parameters:

$$\begin{eqnarray}A & = & 1-{(1-f)}^{2},\ B=2{x}_{{\rm{p}}}{(1-f)}^{2},\\ C & = & {(1-f)}^{2}{R}^{2}-{r}^{2}-{(1-f)}^{2}{x}_{{\rm{p}}}^{2}-{y}_{{\rm{p}}}^{2},\ D=-2{y}_{{\rm{p}}}.\end{eqnarray} \tag{ 24 }$$

Then, an equation for x, the x-coordinate of intersections, is given by

$$\begin{eqnarray}&&{A}^{2}{x}^{4}+2{{ABx}}^{3}+(2{AC}+{B}^{2}+{D}^{2}){x}^{2}\\ &&\quad +2{BCx}+{C}^{2}-{D}^{2}{r}^{2}=0.\end{eqnarray} \tag{ 25 }$$

Equation (25) is a quartic equation, which is analytically solved. We solve this equation with the root-finding method in the same way as before. The number of the solutions for this equation is four at most. In total, there are up to 10 possible solutions for θ_i,l(r).

To test the precision and the computational time in our scheme, we simulate a transit of a Saturn-like planet with R_p/R_⋆ = 0.083667, R_in/p = 1.5, R_out/p = 2.0, θ = π/3, ϕ = π/3, and T = 1.0. We take P = 10759.3 days, a/R_⋆ = 2049.89, b = 0.5, q₁ = 0.49, and q₂ = 0.34 for orbital parameters and stellar parameters. For comparison, we prepare another integration scheme, which adopts pixel-by-pixel integration around the planetary center (e.g., Ohta et al. 2009).

First, we check the precision of the integration of our proposed method by comparing the precision of the pixel-by-pixel integration methods with 5000 × 5000 pixels. As a result, two methods are in agreement to the extent of 10⁻⁷. Thus, our proposed method achieves the numerical error less than 10⁻⁷, which is much smaller than the typical noise in the Kepler data 10⁻⁴.

Second, we check the computational time of our template. Our proposed method typically takes 3.0 ms for calculating one point and 200 s in fitting in Section 6.1. For comparison, we also check the computational time of the planetary transit using the PyTransit package (Parviainen 2015), and we find that it takes 0.3 ms to compute all the 300 data points and 0.3 s in fitting in Section 6.1.

Finally, we compare our method with the pixel-by-pixel integration. If we set the pixel sizes to satisfy the same computational time as that of our method, the precision of the integration becomes 10⁻⁵ in the fiducial configurations. This precision depends on the specific configuration; it becomes 10⁻⁴ if we adopt "R_p = 0.17 and b = 0.8" and 3 × 10⁻⁶ for "R_p = 0.042 and b = 0.3." In summary, when we need a high-precision model, one should use our proposed method, and, if not, one may use the pixel-by-pixel integration to save the amount of calculation.

Incidentally, in a practical case of fitting with the Levenberg–Marquardt algorithm, our method is useful in a sense that gives the smooth value of χ². This is because the smoothness is needed to calculate the differential values for χ² in the LM method.

As we demonstrated in the main text, signatures of a ringed planet can be detected by searching for any deviation from the model light curve assuming a ringless planet. However, the deviation is often very tiny and comparable to the noise level, and so careful quantitative arguments are required to discuss the presence or absence of the ring in a given light curve. In the following, we present a procedure to evaluate the detectability of a ring based on the comparison between the residual from the "planet-alone" model fit and the noise level in the light curve.

Let us denote one light curve including a transit by I_i (i = 0, 1, ⋯, N_data), where N_data is the number of data points. We also define δ_i as the residual of fitting I_i with the planet-alone model. As a quantitative measure of this residual signal δ_i relative to the noise level, we introduce the following S/N:

$$\begin{eqnarray}{\rm{S}}/{\rm{N}} & = & \displaystyle \frac{{\displaystyle \sum }_{i}{\delta }_{i}^{2}}{{\sigma }^{2}}=\displaystyle \frac{{\displaystyle \sum }_{i}{\delta }_{i}^{2}}{{N}_{\mathrm{data}}}\ \displaystyle \frac{{N}_{\mathrm{data}}}{{\sigma }^{2}}\\ & = & {{\rm{\Delta }}}^{2}\displaystyle \frac{1}{{(\sigma /\sqrt{{N}_{\mathrm{data}}})}^{2}}\ \ \ \ \ \ \left({{\rm{\Delta }}}^{2}\equiv \displaystyle \frac{{\displaystyle \sum }_{i}{\delta }_{i}^{2}}{{N}_{\mathrm{data}}}\right).\end{eqnarray} \tag{ 26 }$$

In the last equality, we further define Δ² as the variance of the residual time series, and σ is evaluated as the standard deviation of the out-of-transit light curve. We use the subscript "obs" to specify the above quantities obtained by fitting the planet-alone model to the real observed data: ${\delta }_{i,\mathrm{obs}}$, S/N_obs, and ${{\rm{\Delta }}}_{\mathrm{obs}}^{2}.$

On the other hand, we can also compute the corresponding values of δ_i, S/N, and Δ, by fitting the simulated light curve of a ringed planet with the planet-alone model. We denote these values as ${\delta }_{i,\mathrm{sim}}^{2}(p)$, S/N_sim(p), and ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p)$, where p represents the set of parameters of the ringed-planet model. If these values are sufficiently large compared to the noise variance (see ${{\rm{\Delta }}}_{\mathrm{thr}}^{2}$ below), the signal of the ringed planet is distinguishable from the noise. In addition, comparing these theoretically expected residual levels with observed ones, we can relate the observed residuals to the parameters of the ringed model, even in the absence of clear anomalies.

To simplify the following arguments, we mainly use Δ² instead of S/N to evaluate the significance of the anomaly (see also Appendix B.3 for a detailed reason). Practically, conversion from one to the other is rather simple because the conversion factor $\sigma /\sqrt{{N}_{\mathrm{data}}}$ is well determined from the observed data alone; given a transit light curve, the transit duration T_dur and the bin size t_bin give the number of data points N_data = T_dur/t_bin, and the standard deviation σ can also be inferred from the out-of-transit flux.

For a given region of parameter space p, ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p)$ has the maximum value ${{\rm{\Delta }}}_{\max ,\mathrm{sim}}^{2}$. If ${{\rm{\Delta }}}_{\max ,\mathrm{sim}}^{2}$ is smaller than some threshold value ${{\rm{\Delta }}}_{\mathrm{thr}}^{2}$ determined by the noise level in the light curve, the ringed-planets with the corresponding value of p, even if they exist, cannot be detected in the system. Then, the comparison of ${{\rm{\Delta }}}_{\mathrm{obs}}^{2}$, ${{\rm{\Delta }}}_{\max ,\mathrm{sim}}^{2}$, and ${{\rm{\Delta }}}_{\mathrm{thr}}^{2}$ allows for classification into four categories schematically illustrated in Figure 18:

• (A)  
  ${{\rm{\Delta }}}_{\max ,\mathrm{sim}}^{2}\lt {{\rm{\Delta }}}_{\mathrm{thr}}^{2}$: The expected signal from the ring is so small compared to the noise level that we cannot discuss its detectability.
• (B)  
  ${{\rm{\Delta }}}_{\mathrm{obs}}^{2}\lt {{\rm{\Delta }}}_{\mathrm{thr}}^{2}\lt {{\rm{\Delta }}}_{\max ,\mathrm{sim}}^{2}$: Although the rings with ${{\rm{\Delta }}}_{\mathrm{thr}}^{2}\lt {{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p)$ could have been detected, no significant anomaly is observed (Δ_obs < Δ_thr) in reality. Thus, the parameter region that gives ${{\rm{\Delta }}}_{\mathrm{thr}}^{2}\lt {{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p)$ is excluded.
• (C)  
  ${{\rm{\Delta }}}_{\mathrm{thr}}^{2}\lt {{\rm{\Delta }}}_{\max ,\mathrm{sim}}^{2}\lt {{\rm{\Delta }}}_{\mathrm{obs}}^{2}$: A significant anomaly is detected, but its amplitude is too large to be explained by the ringed-planet model with the given range of p.
• (D)  
  ${{\rm{\Delta }}}_{\mathrm{thr}}^{2}\lt {{\rm{\Delta }}}_{\mathrm{obs}}^{2}\lt {{\rm{\Delta }}}_{\max ,\mathrm{sim}}^{2}$: A significant anomaly is detected, and its amplitude is compatible with the ring model. In this case, we may find the ring parameters consistent with the observed anomaly.

Zoom In Zoom Out Reset image size

The value of ${{\rm{\Delta }}}_{\mathrm{thr}}^{2}$ is arbitrary. In this paper, we choose ${{\rm{\Delta }}}_{\mathrm{thr}}^{2}$ so that it corresponds to S/N = 10 in Equation (26):

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{thr}}^{2}=\displaystyle \frac{10\,{\sigma }^{2}}{{N}_{\mathrm{data}}},\end{eqnarray} \tag{ 27 }$$

where σ and N_data are calculated from the observed data. The methods to calculate the other variances, ${{\rm{\Delta }}}_{\mathrm{obs}}^{2}$, ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p)$, and ${{\rm{\Delta }}}_{\max ,\mathrm{sim}}^{2}$ will be presented in the following subections.

Before proceeding further, let us consider the orbital period dependence of N_data = T_dur/t_bin in Equation (26). From Kepler's third law, ${T}_{\mathrm{dur}}\propto P\,({R}_{\star }/a)\propto {P}^{1/3}$. For the short-period planets, t_bin ∝ P because the number of folded transits is proportional to 1/P. Thus, the number of the data T_dur/t_bin is proportional to P^−2/3. This means that the detectability of rings (S/N) is higher for the shorter-period planets for a given value of Δ². This explains the strong constraints on the ring parameters obtained by Heising et al. (2015) for hot Jupiters.

B.2.1. Definition

The residual ${\delta }_{i,\mathrm{obs}}$ is obtained by fitting the planet-alone model to the data. If the ring does not exist, the value of S/N_obs in Equation (26), which is formally equivalent to the chi-squared, is expected to be close to the degree of freedom dof_obs. Contrary to this scenario, if the ring does not exist, S/N_sim(p) is equal to zero. This means that ${\rm{S}}/{{\rm{N}}}_{\mathrm{obs}}-{\rm{S}}/{{\rm{N}}}_{\mathrm{sim}}(p)\simeq {\mathrm{dof}}_{\mathrm{obs}}$ in the limit of the non-ring system. Thus, for comparison of ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p)$ and ${{\rm{\Delta }}}_{\mathrm{obs}}^{2}$, the value of (${\rm{S}}/{\rm{N}}-{\mathrm{dof}}_{\mathrm{obs}}$) serves as a good estimator of the observed anomaly rather than S/N. We thus slightly modify Equation (26) to define ${{\rm{\Delta }}}_{\mathrm{obs}}^{2}$ so that it corresponds to (${\rm{S}}/{\rm{N}}-{\mathrm{dof}}_{\mathrm{obs}}$):

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{\mathrm{obs}}^{2}=({\chi }^{2}-{\mathrm{dof}}_{\mathrm{obs}})\displaystyle \frac{1}{{(\sigma /\sqrt{{N}_{\mathrm{data}}})}^{2}},\end{eqnarray} \tag{ 28 }$$

where

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i}{\left(\displaystyle \frac{{\delta }_{i,\mathrm{obs}}}{\sigma }\right)}^{2}.\end{eqnarray} \tag{ 29 }$$

The residual ${\delta }_{i,\mathrm{obs}}$ is defined for the best-fit planet-alone model obtained by minimizing χ² as described in Appendix B.2.2 below. The value of χ² is computed using the data just around the transit (within 0.6 T_dur from the transit center) so that the value is not strongly affected by the out-of-transit data. We assume ${\mathrm{dof}}_{\mathrm{obs}}={N}_{\mathrm{data}}\mbox{--}{N}_{\mathrm{para}}\mbox{--}1$, where N_para is the number of fitted parameters.

B.2.2. Detail of Fitting

In fitting, we minimize χ² using the Levenberg–Marquardt algorithm by implementing cmpfit (Markwardt 2009). The adopted model M(t) is composed of a fourth-order polynomial and a transit model F(t):

$$\begin{eqnarray}M(t) & = & F(t)[{c}_{0}+{c}_{1}(t-{T}_{0})+{c}_{2}{(t-{T}_{0})}^{2}\\ & & +{c}_{3}{(t-{T}_{0})}^{3}+{c}_{4}{(t-{T}_{0})}^{4}],\end{eqnarray} \tag{ 30 }$$

where c_i are coefficients of polynomials, and T₀ is a time offset. The polynomials are used to remove the long-term flux variations in the light curve. The transit model F(t) is implemented by the PyTransit package (Parviainen 2015). PyTransit generates the light curves based on the model of Mandel & Agol (2002) with the quadratic limb darkening law.

The above model M(t) includes 12 parameters, ${t}_{0},{R}_{{\rm{p}}}/{R}_{\star },b,a/{R}_{\star },P,{q}_{1},{q}_{2}$, and c_i (i = 0 ∼ 4). For KOIs, the initial values of a/R_⋆, R_p/R_⋆, and b for fitting are taken from the KOI catalog. The initial values of the limb darkening parameters are taken from the Kepler Input Catalog. For a single transit event, where we cannot estimate the orbital period from the transit interval, we choose P, instead of a/R_⋆, as a fitting parameter and estimate a/R_⋆ from P using Kepler's third law and the mean stellar density given in the catalog.

In fitting, we remove outliers iteratively to correctly evaluate χ². We first fit all the data with the model M(t), and flag the points that deviate more than 5σ from the best model. We then refit only the non-flagged data using the same model, and update the flags of all the original data points, including the ones classified as outliers before, on the basis of the new best model and the same 5σ criterion. We iterate this procedure until the flagged data are converged. While this process gives a more robust evaluation of χ², it may also erase the signature of the ringed planet; thus we visually check all the light curves in any case not to miss the real ringed planets.

The noise variance σ² is estimated for each transit light curve by fitting the out-of-transit light curve with a fourth-order polynomial, and calculating the variance of the residuals. Flare-like events are excluded from the estimation of the noise variance.

Since the parameter space p for a ringed planet is very vast, we wish to reduce the volume we need to search with simulations as much as possible. First, we show that ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p)$ does not depend on P and a/R_⋆ with other parameters fixed including limb darkening parameters q₁, q₂, the transit impact parameter b, planet-to-star radius ratio R_p/R_⋆, inner and outer ring radii relative to the planetary radius r_in/p and r_out/p, the direction of the ring (θ, ϕ), and a shading parameter T. This property becomes apparent by rewriting ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p)$ into the following integral form approximately, assuming that the sampling rate (t_bin) is sufficiently small compared to the duration T_dur:

$$\begin{eqnarray}{{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p) & \simeq & \displaystyle \frac{{\displaystyle \int }_{-{T}_{\mathrm{dur}}/2}^{{T}_{\mathrm{dur}}/2}{\delta }_{\mathrm{sim}}^{2}(t,p){dt}}{{\displaystyle \int }_{-{T}_{\mathrm{dur}}/2}^{{T}_{\mathrm{dur}}/2}{dt}}=\displaystyle \frac{{\displaystyle \int }_{-1/2}^{1/2}{\delta }_{\mathrm{sim}}^{2}({T}_{\mathrm{dur}}t^{\prime} ,p){dt}^{\prime} }{{\displaystyle \int }_{-1/2}^{1/2}{dt}^{\prime} }\\ & = & {\displaystyle \int }_{-1/2}^{1/2}{\bar{\delta }}_{\mathrm{sim}}^{2}(t^{\prime} ,p){dt}^{\prime} ,\end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray}&&{\bar{\delta }}_{\mathrm{sim}}(t,p)\equiv {\delta }_{\mathrm{sim}}({T}_{\mathrm{dur}}t,p)\end{eqnarray} \tag{ 32 }$$

and the origin of time is shifted to the transit center. Assuming that the values of q₁, q₂, b, R_p/R_⋆, r_out/p, r_in/p, θ, ϕ, and T are fixed, $\bar{\delta }(t,p)$ defined above does not depend on T_dur explicitly. Therefore, ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p)$ given by Equation (31) does not depend on the timescale of the transit T_dur, which is determined by P and a/R_⋆, and we do not need to simulate the dependence of ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p)$ on these two parameters.

To constrain the parameter space further, we use the observed transit depth. Here we also assume that the values of q₁, q₂, b, T, r_in/p, and the ring direction are fixed and that R_p/R_⋆ and r_out/p are the only free parameters. Then, the constraint on the observed transit depth leaves only one degree of freedom, specified by contours in the R_p/R_⋆-r_out/p plane; henceforth, we rewrite ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}(p)$ as ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}({r}_{\mathrm{out}/{\rm{p}}})$ to explicitly show this dependence.

To compute the relation ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}({r}_{\mathrm{out}/{\rm{p}}})$ for a given transit depth, we first calculate the value of ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}$ and the transit depths for a sufficient number of points in the (r_out/p, R_p/R_⋆) plane. The necessary number of points depends on the fiducial model, and, in our simulation, we prepare about two hundred points for each model in Table 2. For any r_out/p, the observed transit depth uniquely translates into R_p/R_⋆ by the interpolation in the R_p/R_⋆-transit depth plane, because the transit depth is a monotonically increasing function of the R_p/R_⋆. Thus, the given value of r_out/p is uniquely related to ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}$ given the transit depth. By repeating this procedure for many different values of r_out/p, we can compute the relation ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}({r}_{\mathrm{out}/{\rm{p}}})$. We note that once a sufficient number of interpolated lines are prepared, one transit depth determines the relation ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}({r}_{\mathrm{out}/{\rm{p}}})$ without an additional calculation.

Figure 19 shows ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}({r}_{\mathrm{out}/{\rm{p}}})$ curves created in this way, for 4 × 4 = 16 different sets of impact parameters, ring directions, and transit depths. The four sets of p adopted here (model I ∼ model IV) are summarized in Table 2, and four transit depths are chosen to be 0.001, 0.005, 0.01, and 0.05. We fix T = 1 and r_in/p = 1 in all of these simulations.

Zoom In Zoom Out Reset image size

Here we simulate ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}({r}_{\mathrm{out}/{\rm{p}}})$ only for $1\leqslant {r}_{\mathrm{out}/{\rm{p}}}\leqslant {r}_{\mathrm{eq}}$, where r_eq is the value of r_out/p for which the minor axis of the sky-projected outer ring is equal to the planetary radius, computed for each model. This is because the value of ${{\rm{\Delta }}}_{\mathrm{sim}}^{2}({r}_{\mathrm{out}/{\rm{p}}})$ shows no R_p/R_⋆ dependence beyond r_eq, when T = 1 and r_in/p = 1 are adopted; if this is the case, the planetary disk is within the outer disk and the transit depth is solely determined by the latter.

In this paper, we only use the observed constraint on the transit depth. However, this is just for simplicity and we can certainly take into account the constraints on other parameters including b, q₁, and q₂ from the morphology of the observed transit light curve (e.g., egress and ingress durations). Such constraints further restrict the ring models that could be consistent with the observed light curve and thus help more elaborate discussions on the ring parameters, which we leave to future works.