Systematic Search for Rings around Kepler Planet Candidates: Constraints on Ring Size and Occurrence Rate

If a transiting planet has a ring, the ring signature should be imprinted equally in each light curve at different transit epochs. Since the ring parameters, in particular the orientation angles of the ring, are not supposed to vary for the timescale of T_obs, the S/N of the signature should increase by stacking all the light curves properly. To produce such precise phase-folded light curves requires an accurate determination of both the transit center and baseline of each light curve at different transit epochs.

We use the short-cadence presearch data-conditioned simple aperture photometry (PDC-SAP) fluxes of the target objects from the Mikulski Archive for Space Telescopes (MAST). We adopt the transit model F(t) implemented by the Pytransit package (Parviainen 2015) for transiting ringless planets, which generates the light curve based on the model of Mandel & Agol (2002) with the quadratic limb-darkening law.

We first apply the ringless model separately to each transit by varying the transit centers and baseline functions alone. Here we take the fourth-order polynomials as the baseline functions, and we retrieve the transit duration and other parameters of the transiting planets from the MAST pipeline with the help of the Python interface kplr (http://dan.iel.fm/kplr/). We extract the light curve during the epoch of ±2 times the transit duration with respect to each transit center for the subsequent analysis.

After fitting, we exclude outliers exceeding 5σ amplitude in the flux so as to determine the baseline of the light curve accurately. We repeat the fitting procedure and removal of outliers until no outliers are left. Then we visually check each transit in order to exclude inappropriate transits that may be strongly affected by instrumental systematics.

Several transits exhibit large transit-timing variations, which our pipeline cannot automatically deal with. In such cases, we appropriately choose the initial transit centers before fitting so as to correctly identify the transits. Finally, we obtain the best baseline using the out-of-transit (outside ±0.6× transit duration around the transit center) data alone and normalize the light curve with the fitted baseline. Our fit to the transit model light curve is performed with the public code mpfit (Markwardt 2009) that is based on the Levenberg–Marquardt (LM) algorithm.

We stack the obtained normalized light curve at each transit and make the phase-folded light curve. We derive the transit duration by applying the ringless model to the phase-folded light curve. With the updated transit duration, we repeat the above procedure to obtain the final phase-folded light curve.

We extract the phase-folded light curve during an epoch within ±1 transit duration around the transit center. To finish the fitting procedures in realistic time, the light curve is divided into 500 bins with an equal time interval. Here we require one bin to accommodate at least 10 points to guarantee the appropriate binning. So, for systems with a number of the phase-folded data less than 5000, we choose the bin width for one bin to have 10 data points.

Finally, we have phase-folded light curves for 168 planets, which are analyzed for the ring search in the next subsection.

Our search for ring signatures is based on the comparison between the separate best solutions for a planet with and without a ring for all our targets.

In order to find the best solution in the seven-parameter space for a ringless planet model, we randomly generate 1000 different initial sets of parameters from the homogeneous distribution in a finite range. Then, we use the LM method to find the local minima starting from each of the initial values, and we choose the best solution among the solutions. In fitting, we use the binned data that are produced in Section 3.1. We confirm that generally 100 initial sets of parameters are sufficient to find the minimum for our purpose.

Finally, we calculate the χ² value,

$$\begin{eqnarray}&&{\chi }_{\mathrm{ringless}}^{2}=\displaystyle \sum _{{\rm{i}}}{([d({t}_{{\rm{i}}})-m({t}_{{\rm{i}}})]/{\rm{\Delta }}d({t}_{{\rm{i}}}))}^{2},\end{eqnarray} \tag{ 3 }$$

from the binned data. Here $d({t}_{{\rm{i}}})$, $m({t}_{{\rm{i}}})$, and ${\rm{\Delta }}d({t}_{{\rm{i}}})$ are the observed flux, the expected flux of the model, and the uncertainty in observed flux at t = t_i, respectively. We assume ${\rm{\Delta }}d({t}_{{\rm{i}}})$ to be a standard deviation of the normalized flux of each light curve estimated from its out-of-transit epoch.

The same procedure is performed for a ringed-planet model. In this case, we have 12 free parameters: t₀, b, R_p, r_out/in, ${r}_{\mathrm{in}/{\rm{p}}}$, θ, ϕ, T, $a/{R}_{\star }$, c, q₁, and q₂. We calculate the χ² value ${\chi }_{\mathrm{ring}}^{2}$, which has a definition similar to ${\chi }_{\mathrm{ringless}}^{2}$.

One fit of the ringed model takes about a few minutes in a laptop, and the fits to the entire data set were carried out with PC clusters at the Center for Computational Astrophysics (CfCA) at the National Astronomical Observatory, Japan.

Our next procedure is to create a list of tentative ringed-planet candidates from the comparison between the best-fit values for the two models, ${\chi }_{\mathrm{ringless},\min }^{2}$ and ${\chi }_{\mathrm{ring},\min }^{2}$. Specifically for this purpose, we adopt an F-test with F statics (e.g., Lissauer et al. 2011) and define

$$\begin{eqnarray}&&{F}_{\mathrm{obs}}=\displaystyle \frac{({\chi }_{\mathrm{ringless},\min }^{2}-{\chi }_{\mathrm{ring},\min }^{2})/({N}_{\mathrm{ring}}-{N}_{\mathrm{ringless}})}{{\chi }_{\mathrm{ring},\min }^{2}/({N}_{\mathrm{bin}}-{N}_{\mathrm{ring}}-1)},\end{eqnarray} \tag{ 4 }$$

where N_bin is the number of in-transit bins of the phase-folded light curve (typically 500), and N_ring = 12 and N_ringless = 7 are the number of free parameters in the planetary models with and without a ring, respectively.

The numerator of the right-hand side of Equation (4) corresponds to the improvement in χ² of the ring model divided by the number of the additional degrees of freedom characterizing a ring. The denominator is the χ² per degree of freedom for the ringed model. Thus, F_obs represents a measure of relative improvement of the fit by introducing the ring. The large F_obs prefers the ringed-planet model. Note, however, that F_obs is defined simply through the ratio of the minimum values of χ² for the two models. Therefore, it has nothing to do with the goodness of the fit for either model, which needs to be checked separately.

According to the F-test, the measure of the null hypothesis that our ringed model does not improve the fit relative to the ringless model is given by the p-value, defined as

$$\begin{eqnarray}&&p=1-{\int }_{0}^{{F}_{\mathrm{obs}}}F(f| {N}_{\mathrm{bin}}-{N}_{\mathrm{ring}}-1,{N}_{\mathrm{ring}}-{N}_{\mathrm{ringless}}){df},\end{eqnarray} \tag{ 5 }$$

where $F(f| {N}_{\mathrm{bin}}-{N}_{\mathrm{ring}}-1,{N}_{\mathrm{ring}}-{N}_{\mathrm{ringless}})$ is the F-distribution with the degrees of freedom $({N}_{\mathrm{bin}}-{N}_{\mathrm{ring}}-1$, ${N}_{\mathrm{ring}}-{N}_{\mathrm{ringles}})$.

A large value of F_obs, or a small value of p, statistically disfavors the current null assumption, i.e., the ringed model better fits the data than the ringless model. In this paper, we adopt the condition of p < 0.05 for the rejection of the null hypothesis. For those tentative candidates of ringed planets, we attempt to understand the origins of anomalies by examining individual light curves and statistics (e.g., ${\chi }_{\mathrm{ring},\min }^{2})$ further.

We also test the robustness of possible ring signatures by dividing the multiple transits into those at even and odd transit numbers, creating the phase-folded light curves separately, and computing the p-values $({p}_{\mathrm{even}},{p}_{\mathrm{odd}})$. Unlike the other analyses, we use the nonbinned data here in order to evade the additional uncertainties in the light curves due to the extra binning step, especially for systems with a low number of data. For the calculation of (p_even, p_odd), we approximate the best-fit model of the binned data as that of the nonbinned data, and then we calculate F_obs in Equation (4) for nonbinned data. If rings mainly account for signals in light curves, we expect p_even to be close to p_odd because of the consistency of the signals.

Finally, we comment on the validity of applying the F-test to our ring search. The F-test needs to satisfy two conditions (e.g., Protassov et al. 2002). One is that the two models are nested in a sense that the more complicated model reduces to the simpler one if the additional parameters in the former model are removed. This is trivially satisfied in the present case. The other condition is that the simpler model should not be located at the edge of the parameter space of the more complicated model. Strictly speaking, this condition may not hold because our ring model reduces to the ringless model in the limit of ${R}_{\mathrm{out}}\to {R}_{{\rm{p}}}$. Nevertheless, the F-test gives us a practically useful criterion, and we decide to use it in selecting tentative candidates for further analysis.

Even for planetary systems without any detectable signatures of a ring, we may constrain the property of a possible ring within the observational detection limit. To proceed realistically, we need to reduce the number of free parameters characterizing the ring. Thus, we fix the inner radius of the ring as R_in = R_p and set the opacity of the ring as T = 1 just for simplicity. Furthermore, we focus on two cases for the orientation angles of the ring, as we describe in the next subsections. Thus, we are left with a single parameter, the outer radius of the ring R_out. In practice, we place upper limits on the ratio R_out/R_p from the fit to the light curves.

3.4.1. Aligned with the Planetary Orbit

Under the strong tidal interaction with the star, the ring becomes aligned to the orbital plane of the planet. Indeed, Brown et al. (2001) gave the upper limit on the ring size of a hot Jupiter, HD 209458 b, as 1.7R_p assuming the alignment.

In a similar manner, we place upper limits on the ring size assuming the tidal alignment. The tidal alignment leads to the orientation of $\theta =\arcsin (b/(a/{R}_{\star }))$ and ϕ = 0. In addition, the small value of θ enhances the effective optical depth viewed from the observer relative to that from the top view. Thus, we assume T = 1 even though rings can be very thin, like Jupiter's rings.

In summary, we fix ϕ = 0°, $\theta =\arcsin ({{bR}}_{\star }/a)$, T = 1, and R_in = R_p for fitting. Assuming these conditions, we fit the ringed model to the data using at least 100 sets of randomly chosen initial parameters, and we pick up the best solution among the local optimum solutions.

After obtaining the best solutions with fixed values of R_out/R_p, we define the 3σ limit ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp}}$, where

$$\begin{eqnarray}{\rm{\Delta }}{\chi }^{2}({R}_{\mathrm{out}}/{R}_{{\rm{p}}}) & \equiv & ({\chi }_{\mathrm{ring},\min }^{2}({R}_{\mathrm{out}}/{R}_{{\rm{p}}})\\ & & -{\chi }_{\mathrm{ringless},\min }^{2})/({\chi }_{\mathrm{ringless},\min }^{2}/\mathrm{dof})\end{eqnarray} \tag{ 6 }$$

becomes 9. In practice, we compute ${\rm{\Delta }}{\chi }^{2}({R}_{\mathrm{out}}/{R}_{{\rm{p}}})$ at 11 values of R_out/R_p: 1.1, 1.3, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, and 10.0. Then, we interpolate them to find ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp}}$.

Our procedure for setting the upper limit is illustrated in Figure 2 for KOI-97.01. In this example, the interpolated curve crosses the ${\rm{\Delta }}{\chi }^{2}=9$ threshold at R_out/R_p = 1.55. Thus, we obtain ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Aligned}}=1.55$ for KOI-97.01. Figure 3 plots three corresponding fitting curves with R_out/R_p = 1.5, 2.0, and 2.5 along with the curve of the ringless model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If Δχ² < 9 for R_out/R_p = 10.0, we do not place upper limits ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp}}$. These cases are marked as ⋯ in Tables 1–3 below.

Table 1.  Parameters and Statistics of 14 Systems with p < 0.003

KOI	Kepler	Porba (days)	tdamp (Gyr)	${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Aligned}}$	${({R}_{{\rm{p}}}/{R}_{\star })}_{\mathrm{ringless}}$	p	(${\chi }_{\mathrm{ringless},\min }^{2}$, ${\chi }_{\mathrm{ring},\min }^{2}$, Nbin)	(S/N)	Commentb
2.01	2 b	2.20	2.60e–05	1.11	0.07755 ± 2e–05	5.02e–07	(534.22, 494.94, 500)	4357.03	GD
3.01	3 b	4.89	8.76e–04	2.75	0.05886 ± 3e–05	2.60e–09	(2011.23, 1820.95, 500)	2403.33	Spot
13.01	13 b	1.76	7.92e–06	1.63	0.064683 ± 4e–06	<1e–10	(8830.73, 5904.78, 500)	6359.84	GD
63.01	63 b	9.43	5.18e–03	2.44	0.06481 ± 4e–05	<1e–10	(1222.04, 1078.89, 500)	732.29	Spot
70.02	20 b	3.70	1.23e–03	2.47	0.01799 ± 9e–05	1.34e–03	(671.05, 644.41, 500)	128.93	Bad fold
102.01	⋯	1.74	7.01e–05	5.51	0.02810 ± 6e–05	9.78e–05	(554.11, 525.66, 500)	432.74	Bad fold
676.01	210 c	7.97	5.35e–03	7.20	0.0520 ± 5e–04	2.00e–10	(693.78, 620.69, 500)	302.61	Spot
868.01	⋯	236.00	2.95e+01	6.65	0.144 ± 1e–03	7.88e–06	(202.04, 171.63, 203)	182.23	Others
971.01	⋯	0.53	9.00e–10	8.13	0.1 ± 1e+00	3.54e–04	(319.16, 304.58, 500)	257.01	FP and others
1416.01	840 b	2.50	3.55e–05	1.74	0.1459 ± 2e–04	5.14e–04	(579.31, 553.80, 500)	919.15	FP and spot
1539.01	⋯	2.82	1.20e–05	1.10	0.2568 ± 2e–04	<1e–10	(982.31, 791.58, 500)	1364.17	FP and spot
1714.01	⋯	2.74	5.29e–06	1.10	0.17618 ± 2e–05	<1e–10	(6978.85, 2755.22, 500)	688.33	FP and spot
1729.01	⋯	5.20	2.88e–04	1.80	0.1764 ± 3e–04	<1e–10	(688.30, 610.78, 500)	816.17	FP and spot
3794.01	1520 b	0.65	3.27e–07	⋯	0.101 ± 3e–03	2.37e–06	(678.22, 632.70, 500)	265.86	Evap

Notes. ^aValues from Kepler Object of Interest (KOI) Catalog Q1–Q17 DR 25 (https://exoplanetarchive.ipac.caltech.edu/). ^bFP = possible false positive (https://exofop.ipac.caltech.edu/cfop.php); GD = gravity darkening (Appendix B.1); Evap = evaporating planet (Appendix B.2); Spot = spot crossing (Appendix B.3); Bad fold = incorrect data folding (Appendix B.4); Small = nonsignificant signal (Appendix B.5); Others = Appendix B.6.

Download table as:  ASCIITypeset image

The alignment condition is determined by the tidal dissipation function Q_p and the Love number k_p, the planet/star mass ratio, the dimensional moment of the inertia of the planet C, the orbital period P_orb, and the normalized semimajor axis a/R_⋆. As discussed in Appendix A, 154 out of the 168 planetary systems are supposed to become aligned within a timescale of 1 Gyr, if we adopt the fiducial values Q_p = 10^6.5, k_p = 1.5, and C = 0.25 and the mass–radius relation (Equation (8) in Weiss et al. 2013). We compute the upper limit on R_out/R_p for all 168 systems in any case, even if their alignment timescale is long.

3.4.2. Orientation of the Saturnian Ring

We have performed a ring search following the method described in Section 3. The results for all 168 Kepler objects are summarized in Tables 1–3.

We identify 29 candidate objects with p-values less than the threshold value of 0.05. For most of these systems, the ring model yields ${\chi }_{\mathrm{ring}}^{2}/\mathrm{dof}$ ∼ 1 (Figure 4). However, after inspecting the individual light curves of these systems, we conclude that none of them is a viable candidate for a ringed planet. Of the 29 candidates, 11 do not exhibit any convincing ringlike signatures in the light curves and so are excluded. The other 18 systems do show anomalous features in the light curves, but they are most likely ascribed to other mechanisms: gravity darkening (two systems), spot crossing (nine systems), disintegration of a planet (one system), artifacts generated during the folding process (three systems), and stellar activity (three systems). See Appendix B for the details of this process, as well as for individual light curves.

Zoom In Zoom Out Reset image size

4.2.1. Result

Given the null detection, we derive upper limits on the outer radius of the possible ring following the method described in Section 3.4. The resulting upper limits, ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Aligned}}$, ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Saturn}}$, and ${({R}_{\mathrm{out}}/{R}_{\star })}_{\mathrm{upp}}$, are listed in Tables 1–3. If we cannot obtain upper limits due to poor S/Ns, we leave those values blank in the tables. The following discussions exclude 18 systems for ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Aligned}}$ and seven systems for ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Saturn}}$ that are identified as possible false positives on the Kepler Community Follow-up Program (CFOP) website.⁶

Figure 5 compares the upper limits ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp}}$ for the aligned and Saturn-like configurations against the physical planetary radii. The latter values are computed as ${({R}_{{\rm{p}}}/{R}_{\star })}_{\mathrm{ringless}}\times {R}_{\star }$, where the values of ${({R}_{{\rm{p}}}/{R}_{\star })}_{\mathrm{ringless}}$ are obtained from the ringless model and the stellar radii are taken from the Kepler catalog. Even assuming the ring aligned with the orbital plane, we find fairly tight limits on the ring size (several times R_p) for a few tens of systems.

Zoom In Zoom Out Reset image size

Figure 6 is similar to Figure 5 but against the equilibrium temperatures T_eq of the planets. The exhibited pattern does not reflect the physical dependence of ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp}}$ on T_eq but simply comes from the fact that the hotter planets have shorter orbital periods and hence larger S/Ns of the phase-folded light curve. With sufficient S/Ns for future data, however, such plots would provide interesting constraints on the physical properties of rings as a function of melting temperature of different compositions.

Zoom In Zoom Out Reset image size

As mentioned in Section 4.1, the light curves of some of the 29 systems with p < 0.05 include contributions from effects other than rings, such as gravity darkening and spot crossing. Nevertheless, we neglect them in deriving the upper limits on R_out/R_p. If we fit and remove those effects from the light curve, the upper limits may become more stringent. In this sense, the upper limits on R_out/R_p listed in Tables 1 and 2 would be a bit conservative.

Table 2.  Parameters and Statistics of 15 Systems with 0.003 < p < 0.05

KOI	Kepler	Porba (days)	tdamp (Gyr)	${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Aligned}}$	${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Saturn}}$	${({R}_{\mathrm{out}}/{R}_{\star })}_{\mathrm{upp}}$	${({R}_{{\rm{p}}}/{R}_{\star })}_{\mathrm{ringless}}$	p	(${\chi }_{\mathrm{ringless},\min }^{2}$, ${\chi }_{\mathrm{ring},\min }^{2}$, Nbin)	(S/N)	Commentb
4.01	⋯	3.85	1.76e–04	⋯	⋯	⋯	0.0394 ± 3e–04	3.44e–02	(538.39, 525.32, 500)	148.29	FP and small
5.01	⋯	4.78	6.05e–04	⋯	⋯	⋯	0.04 ± 1e–02	3.15e–02	(481.35, 469.45, 500)	455.35	FP and small
12.01	448 b	17.86	1.44e–02	1.86	1.23	0.107	0.09018 ± 5e–05	1.23e–02	(548.57, 532.47, 500)	792.22	Others
148.01	48 b	4.78	2.46e–03	9.51	1.44	0.026	0.0196 ± 1e–04	1.13e–02	(663.02, 643.29, 500)	102.75	Bad fold
212.01	⋯	5.70	1.04e–03	⋯	1.72	0.093	0.0649 ± 3e–04	2.14e–02	(525.10, 511.11, 500)	159.08	Small
214.01	424 b	3.31	1.46e–04	2.99	⋯	⋯	0.104 ± 3e–03	2.37e–02	(446.95, 435.27, 500)	448.51	Small
257.01	506 b	6.88	5.14e–03	9.25	2.22	0.033	0.0224 ± 2e–04	3.84e–02	(550.10, 537.05, 500)	161.18	Small
423.01	39 b	21.09	2.12e–02	3.89	1.85	0.129	0.0890 ± 6e–04	1.66e–02	(500.13, 486.18, 500)	225.79	Small
433.02	553 c	328.24	7.94e+01	6.10	⋯	⋯	0.120 ± 7e–03	2.07e–02	(250.09, 238.97, 303)	112.04	Small
531.01	⋯	3.69	2.97e–04	⋯	⋯	⋯	0.096 ± 4e–03	2.44e–02	(469.53, 457.33, 500)	183.64	Small
686.01	⋯	52.51	3.36e–01	4.68	⋯	⋯	0.118 ± 4e–03	4.38e–02	(46.09, 31.21, 40)	153.64	FP and small
872.01	46 b	33.60	1.36e–01	⋯	⋯	⋯	0.084 ± 2e–03	4.57e–02	(486.30, 475.20, 500)	173.00	Small
1131.01	⋯	0.70	2.08e–07	1.02	1.02	0.211	0.2 ± 7e–02	3.84e–02	(549.24, 536.21, 500)	727.17	FP and small
1353.01	289 c	125.87	4.46e+00	4.23	2.00	0.153	0.1048 ± 6e–04	9.04e–03	(544.28, 525.28, 442)	198.14	Spot
6016.01	⋯	4.55	7.76e–05	1.32	1.65	0.463	0.23 ± 2e–02	9.91e–03	(503.24, 487.04, 472)	1719.88	FP and spot

Notes. ^aValues from Kepler Object of Interest (KOI) Catalog Q1–Q17 DR 25 (https://exoplanetarchive.ipac.caltech.edu/). ^bFP = possible false positive (https://exofop.ipac.caltech.edu/cfop.php); GD = gravity darkening (Appendix B.1); Evap = evaporating planet (Appendix B.2); Spot = spot crossing (Appendix B.3); Bad fold = incorrect data folding (Appendix B.4); Small = nonsignificant signal (Appendix B.5); Others = Appendix B.6.

Download table as:  ASCIITypeset image

4.2.2. Comparison of Roche Radius and Upper Limits

To understand the implications of the upper limits physically, we compare the limits ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp}}$ with the Roche radii. If we consider ring formation by tidal destruction of incoming objects (e.g., satellites), the outer radius of the ring may be set by the Roche radius,

$$\begin{eqnarray}&&{R}_{\mathrm{out}}\sim 2.45\,{R}_{{\rm{p}}}{\left(\displaystyle \frac{{\rho }_{{\rm{p}}}}{{\rho }_{{\rm{s}}}}\right)}^{1/3}=1.6\,{R}_{{\rm{p}}}{\left(\displaystyle \frac{{\rho }_{{\rm{p}}}/1{\rm{g}}{\mathrm{cm}}^{-3}}{{\rho }_{{\rm{s}}}/3.5{\rm{g}}{\mathrm{cm}}^{-3}}\right)}^{1/3},\end{eqnarray} \tag{ 7 }$$

where ρ_p is the planetary density and ρ_s is that of the incoming object. Here we scale the result using ${\rho }_{{\rm{p}}}=1\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and ${\rho }_{{\rm{s}}}=3.5\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, which are the typical values for rocky components in the solar system.

This implies that, if the inferred upper limit on the ring size ${R}_{\mathrm{out}}/{R}_{{\rm{p}}}$ is much smaller than 1.6, the ring is unlikely to exist even inside that limit—unless ρ_s is unreasonably large. In our sample, six systems satisfy t_damp < 1 Gyr and ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Aligned}}\lt 1.6$, and one satisfies ${t}_{\mathrm{damp}}\gt 1$ Gyr and ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Saturn}}\lt 1.6$. We may exclude possible rings around these systems.

The above limits on R_out/R_p translate into the upper limit on the occurrence rate of rings q[>x] as a function of $x\equiv {R}_{\mathrm{out}}/{R}_{{\rm{p}}}$. Here q[>x] is the probability that a planet has a ring larger than x times the planetary radius. For example, q[>x = 1] is simply the occurrence rate of rings, and q[>x = 2] is that of rings larger than twice the planetary radii.

We attempt to estimate the upper limit on q[>x] as follows. For a given value of x, consider n samples extracted from systems with q[>x], for which the rings with R_out/R_p > x should have been readily detectable—so this may be chosen to be N[<x], the number of systems with ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp}}\lt x$. Then, the probability that we detect n_obs rings with R_out/R_p ≥ x out of the n samples is given simply by the binominal distribution:

$$\begin{eqnarray}&&\mathrm{Prob}{({n}_{\mathrm{obs}}| q{[\gt x],n){=}_{n}{C}_{{n}_{\mathrm{obs}}}q[\gt x]}^{{n}_{\mathrm{obs}}}(1-q[\gt x])}^{n-{n}_{\mathrm{obs}}}.\end{eqnarray} \tag{ 8 }$$

Without any prior knowledge of q[>x] or ${n}_{\mathrm{obs}}$, we assume the uniform distribution for $\mathrm{Prob}(q[\gt x])$ and $\mathrm{Prob}({n}_{\mathrm{obs}})$ with proper normalizations:

$$\begin{eqnarray}&&{\int }_{0}^{1}\mathrm{Prob}(q[\gt x]| n)\,{dq}[\gt x]=1\to \mathrm{Prob}(q[\gt x]| n)=1,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\displaystyle \sum _{{n}_{\mathrm{obs}}=0}^{n}\mathrm{Prob}({n}_{\mathrm{obs}}| n)=1\to \mathrm{Prob}({n}_{\mathrm{obs}}| n)=1/(n+1).\end{eqnarray} \tag{ 10 }$$

According to Bayes' theorem, we obtain

$$\begin{eqnarray}&&\mathrm{Prob}(q[\gt x]| {n}_{\mathrm{obs}}=0,n)\\ &&\quad =\,\displaystyle \frac{\mathrm{Prob}(q[\gt x]| n)\mathrm{Prob}({n}_{\mathrm{obs}}=0| q[\gt x],n)}{\mathrm{Prob}({n}_{\mathrm{obs}}=0| n)}\\ &&\quad =\,(n+1){(1-q[\gt x])}^{n}.\end{eqnarray} \tag{ 11 }$$

The corresponding cumulative distribution function for q[>x] is given by:

$$\begin{eqnarray}&&\mathrm{CDF}{(q[\gt x])=1-(1-q[\gt x])}^{n+1}.\end{eqnarray} \tag{ 12 }$$

Here we would like to obtain the 95% upper limits of q[>x]. Thus, the above equation gives

$$\begin{eqnarray}&&q{[\gt x]}_{\mathrm{upp}}=1-{(0.05)}^{\tfrac{1}{n+1}}.\end{eqnarray} \tag{ 13 }$$

Now we substitute the values of N[<x] plotted in Figure 6 into n and obtain the upper limits of q[>x] as a function of x.

Figure 7 shows $q{[\gt x]}_{\mathrm{upp}}$ using $N[\lt {({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Aligned}}]$ and $N[\lt {({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Saturn}}]$. Physically speaking, the limit ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Aligned}}$ is appropriate only for systems with small values of t_damp, which have likely achieved tidal alignment. On the other hand, ${({R}_{\mathrm{out}}/{R}_{{\rm{p}}})}_{\mathrm{upp},\mathrm{Saturn}}$ may be more relevant for those with large values of t_damp. Therefore, we distinguish the systems with ${t}_{\mathrm{damp}}\lt 1\,\mathrm{Gyr}$ and t_damp > 1 Gyr in the plot. The more relevant subset is shown with thick lines in each panel.

Zoom In Zoom Out Reset image size

Fast-rotating stars have higher (lower) effective surface temperatures in the polar (equatorial) regions because of the stronger centrifugal force along the equatorial plane. Thus, the transit light curve becomes asymmetric with respect to the central transit time, depending on the path of the planet. The anomaly due to gravity darkening is not confined preferentially around the ingress or egress phases, unlike the ring signature (e.g., see Figure 3), and can be distinguished easily by eye.

Figure 8 shows a light curve of our tentative candidate KOI-13.01 (Kepler-13 b), which cannot be well fitted anyway even by adding a ring. This system was analyzed first by Barnes et al. (2011), who found that the light curve is very well explained by gravity darkening. Masuda (2015) presented a further elaborated analysis of KOI-13 (Kepler-13), as well as another gravity-darkened system, KOI-2, in our targets.

Zoom In Zoom Out Reset image size

Another tentative candidate, KOI-3794.01 (KIC 12557548, Kepler-1520 b), is known as an evaporating planet (e.g., Rappaport et al. 2012), whose light curve is shown in Figure 9. Indeed, the transit depth of the light curves at different epochs (before phase-folded) exhibits significant time variation, which is inconsistent with the ring hypothesis.

Zoom In Zoom Out Reset image size

Stellar spots add nonnegligible anomalous features in the transit light curves. Among the 29 tentative candidates with p < 0.05, we find that nine systems are likely explained by spot-crossing events, not by a ring. As a significant example, we show the phase-folded light curve of KOI-1714.01 in Figure 10, where the entire flux is strongly affected by spot-crossing events.

Zoom In Zoom Out Reset image size

Spot-crossing features have already been reported for four out of nine systems; KOI-3.01 (Kepler-3 b) shows frequent spot-crossing anomalies at fairly similar phases, and its planetary orbit is estimated to be misaligned relative to the stellar spin (Sanchis-Ojeda & Winn 2011). Combining the spot anomalies and the Rossiter–McLaughlin effect of KOI-63.01 (Kepler-63 b), Sanchis-Ojeda et al. (2013) concluded that the system has a large spin–orbit misalignment of Ψ = 104°. Also, the variability of light curves due to spot-crossing events has been reported for KOI-676.01 (Kepler-210 c) by Sanchis-Ojeda et al. (2013) and KOI-1353.01 (Kepler-289 c) by Schmitt et al. (2014).

The other five systems—KOI-1416.01 (Kepler-850 b), 1539.01, 1714.01, 1729.01, and 6016.01—are classified as possible false positives on the Kepler CFOP website, and we confirmed that there are no ringlike signatures.

The phase-folded light curves of KOI-70.02, 102.01, and 148.01 show anomalous features around the egress and ingress phases. The transit depth of those three systems is very small, and we suspect that the anomalies are simply caused by inaccurate central transit epochs in phase folding.

Figure 11 shows an example for KOI-148.01. In the left panel, we show the light curve, which is folded as described in Section 3.1. As shown in the left panel, the anomalous features appear around the egress and ingress. Then, to find out the origin of the anomaly, we create a phase-folded light curve using a linear ephemeris. Specifically, when we fit the individual transit, we fix each transit center to ${t}_{\mathrm{cen},{\rm{i}}}\,={t}_{\mathrm{cen},0}+{{iP}}_{\mathrm{orb}}$, where ${t}_{\mathrm{cen},{\rm{i}}}$ is the transit center at the ith transit. Here we retrieve ${t}_{\mathrm{cen},0}$ and P_orb from the Kepler catalog. The refolded light curve is plotted in the right panel, which shows that the anomalous features disappear. We made sure that this is also the case for the other two systems, KOI-70.02 and 102.01.

Zoom In Zoom Out Reset image size

Out of the 168 targets, we find 11 systems that marginally favor the ring model at 2σ–3σ levels: KOI-4.01, 5.01, 212.01, 214.01, 257.01, 423.01, 433.02, 531.01, 686.01, 872.01, and 1131.01. Figure 12 shows their light curves, as well as the best-fit model with and without a ring (blue and red lines, respectively).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To examine their significance, we divide their individual transit light curves into two groups, as described in Section 3.3. If the anomaly is really caused by a ring, both p_even and p_odd should remain small.

We find that nine systems have both p_even and p_odd larger than 0.32 (i.e., 1σ), and two systems have both p_even and p_odd with merely 1σ significance: (0.11, 0.12) for KOI-4.01 and (0.029, 0.037) for KOI-257.01. Even though the two systems are likely to be statistical flukes, we examined the light curve visually in any case. The light curve of KOI-4.01 seems marginally consistent with the ring signature, but the amplitude is so tiny and can be easily produced by random noise. The features of KOI-257.01 are likely to be produced by the folding procedure we discussed in Appendix B.4 because the transit depth is so small.

We note that the rejection of the null hypothesis of a ring with a level of p = 0.05 implies that 168 × 0.05 =8.4 systems are expected to show 2σ signals even if there is no ring at all. Thus, 11 marginal systems, even if there is no ring system, are fairly consistent with our choice of the threshold.

Finally, we consider the remaining three systems that have not been discussed.

The light curve of KOI-12.01 (Kepler-448 b) shows anomalous features during the transit that are significant during −0.03 to 0.1 day with respect to the central transit epoch. We find that such large pulse-like signals also appear during out-of-transit. Thus, they are likely due to stellar activities.

The light curves of KOI-971.01 (KIC 11180361) show strong stellar activities that are typical for multiple-star systems (Niemczura et al. 2015), and the CFOP website also identifies this system as a false positive. Thus, the planetary rings are not the origin of the signals.

The light curve of KOI-868.01 shows an anomaly during the egress, which is shown in the left panel of Figure 13, along with the best-fit models. The fit yields ${\chi }_{\mathrm{ringless},\min }^{2}/\mathrm{dof}\,=202.0/190$, ${\chi }_{\mathrm{ring},\min }^{2}/\mathrm{dof}=171.6/195$, and $p=7.88\,\times {10}^{-6}$. The analysis based on the binned data supports a Neptune-sized ringed planet of an orbital period of 236 days. The best-fit ring model gives $\theta =25\buildrel{\circ}\over{.} 5\pm 10\buildrel{\circ}\over{.} 0$, ϕ = 124 ± 37, T = 0.46 ± 0.18, ${r}_{\mathrm{in}/{\rm{p}}}=1.88\pm 0.36$, and ${r}_{\mathrm{out}/\mathrm{in}}=1.63\pm 0.43$. The radius ratio ${R}_{{\rm{p}}}/{R}_{\star }=0.099\pm 0.012$ gives ${R}_{{\rm{p}}}/{R}_{{\rm{J}}}=0.63\pm 0.08$, assuming the stellar radius ${R}_{\star }={0.657}_{-0.032}^{+0.022}\,{R}_{\odot }$. The nonvanishing obliquity is consistent with the long alignment timescale t_damp = 2.95 Gyr.

Zoom In Zoom Out Reset image size

In order to check the consistency of the signals, we calculate the p-values for the two transits in the short-cadence data separately. As a result, we find p = 0.76 and 7.4e–07 for the first and second transits, respectively. Indeed, as indicated in the right panel of Figure 13, the light curves at the first and second transits are systematically different. Therefore, KOI-868.01 is unlikely to be a ringed planet. We do not understand the origin of the anomalies because there are only two transits, but we suspect that temporal stellar activities or spot-crossing events are responsible.