Robust Transiting Exoplanet Radii in the Presence of Starspots from Ingress and Egress Durations

We typically measure the radii of transiting exoplanets from the transit depth, which is given by the ratio of cross-sectional areas of the planet and star. However, if a star has dark starspots (or bright regions) distributed throughout the transit chord, the transit depth will be biased towards smaller (larger) values, and thus the inferred planet radius will be smaller (larger) if unaccounted for. We reparameterize the transit light curve to account for"self-contamination"by photospheric inhomogeneities by splitting the parameter $R_p/R_\star$ into two parameters: one for the radius ratio -- which controls the duration of ingress and egress -- and another which measures the possibly contaminated transit depth. We show that this is equivalent to the formulation for contamination by a second star (with positive or negative flux), and that it is sensitive to time-steady inhomogeneity of the stellar photosphere. We use synthetic light curves of spotted stars at high signal-to-noise to show that the radius recovered from measurement of the ingress/egress duration can recover the true radii of planets transiting spotted stars with axisymmetric spot distributions if the limb-darkening parameters are precisely known. We fit time-averaged high signal-to-noise transit light curves from Kepler and Spitzer of ten planets to measure the planet radii and search for evidence of spot distributions. We find that this sample has a range of measured depths and ingress durations which are self-consistent, providing no strong evidence for contamination by spots. However, there is suggestive evidence for occultation of starspots on Kepler-17, and that relatively bright regions are occulted by the planets of Kepler-412 and HD 80606. Future observations with the James Webb Space Telescope may enable this technique to yield accurate planetary radii in the presence of stellar inhomogeneities.


INTRODUCTION
bmmorris@uw.edu * Guggenheim Fellow Precise and accurate radii of exoplanets are critical to understanding their habitability. Radii are often used as a proxy for bulk composition when masses aren't available, informed by population studies which suggest that planets smaller than R 1.5R ⊕ are more likely rocky than gaseous (Weiss & Marcy 2014;Rogers 2015;Fulton et al. 2017). Since bulk densities are sensitive to the radius to the third power, accurate radii are required in order to infer accurate bulk densities. The radii of transiting exoplanets are typically constrained by the transit depth, which is interpreted as the ratio of the cross-sectional areas of the planet and star ∆F/F = (R p /R ) 2 (ignoring limb-darkening). The transit depth accuracy as a measure of the planet's radius, however, is sensitive to stellar surface inhomogeneities. When planets pass over dark star spots, for example, the transit light curve is biased towards shallower transit depths during the transit (see e.g. Csizmadia et al. 2013;Oshagh et al. 2013). As we discover more planetary systems of late type stars like TRAPPIST-1 (Gillon et al. 2016(Gillon et al. , 2017Luger et al. 2017), the time-and wavelength-dependent variability of these active stars and non-uniformity is likely to limit the accuracy and precision of our measurements of the planet radii, especially as a function of wavelength (Roettenbacher & Kane 2017;Rackham et al. 2018;Zhang et al. 2018; Morris et al. 2018a).
Particularly insidious are nearly axisymmetric stellar inhomogeneities and/or inhomogeneities which are comprised of odd harmonics. These do not lead to variability in the total flux of the star, and do not necessarily cause spot-crossing features in the transit light curve; yet they will still affect the depth of transit, and hence the inferred planet-star radius ratio and transmission spectrum. Consequently, even without indications of star-spot activity in the light curve of a star, spots (or other inhomogeneities) can still affect the depths of transits; this phenomenon is referred to as the "transit light source effect" in the recent work by Rackham et al. (2018). These forms of time-steady "activity" might be due to a multitude of star spots (such as along a particularly active latitude), due to gravity darkening, or due to a strong global magnetic field. The message of this paper is that if one measures the duration of ingress as compared to the total transit duration, this time-dependent quantity will be less affected by the presence of spots on the rest of the photosphere, and thus may give a more precise, albeit less accurate, measurement of the planet-star radius ratio compared with measuring the radius ratio from the depth of transit. However, in most cases it does require an independent constraint upon the impact parameter of the transit, as well as a careful handling of limb-darkening.
In Section 2 we reparameterize the transit light curve formalism of Mandel & Agol (2002) for more robust planet radii in the presence of significant starspot coverage, and elaborate on the degeneracies in the model in Section 3. In Section 4, we validate the new transit model with fits to synthetic light curves of simulated spotted stars. We fit real transit light curves of a variety of planets including TRAPPIST-1 and GJ 1214 in Section 5, and discuss the results and conclusions in Sections 6 and 7.

SELF-CONTAMINATION AND THE
TRANSIT LIGHT CURVE Following Mandel & Agol (2002), the ratio of unocculted flux to occulted flux F is given by F e (p 0 , z, u 1 , u 2 ) = 1 − λ e (p 0 , z, u 1 , u 2 ), (1) where p 0 = R p /R , z is the projected sky separation between the centers of the star and planet (in units of the stellar radius), and u 1 and u 2 are the quadratic limb-darkening parameters. One way to interpret the transit light curve -for a star known to have little contamination from other nearby stars -is to create a flux contamination term which will account for bright or dark regions of the surface of the star known to be transited. We introduce the new parameter p 1 and renormalize the λ term as: F e (p 1 , p 0 , z, u 1 , u 2 ) = 1− p 1 p 0 2 λ e (p 0 , z, u 1 , u 2 ).
(2) It follows that p 1 ≈ √ δ where δ is the transit depth (the approximation becomes an equivalence for the uniform case where u 1 = u 2 = 0). We emphasize here that p 1 is not a radius ratio parameter like p 0 ; p 2 1 is related to the observed transit depth in the presence of bright or dark spot contamination. A measured difference between p 0 and p 1 is in principle sensitive to starspots, even with an axisymmetric distribution, unlike rotational modulation.
This renormalization is equivalent to the analysis often performed on Kepler planet candidate host stars in order to search for contaminating light from a second, unresolved star (Torres et al. 2011;Morton 2012;Teske et al. 2018). Typically in the external contamination case, the light from the background star dilutes the flux of the exoplanet host star, decreasing the transit depth, modifying the inferred impact parameter, and leading to underestimated the planet radii, if the contamination is unaccounted for.
In equation 2, we are essentially allowing for a contaminating light source on the host star itself, which can be positive or negative due to bright or dark regions on the stellar surface, thus increasing or decreasing the transit depth p 2 1 with respect to the expectation given the planet's radius p 0 . Typically for contaminating nearby sources, the ratio of the unocculted flux to occulted flux is given by where is the flux of the contaminating source. In this parameterization, We focus our attention in this work on the "selfcontamination" effect of bright or dark regions on the host star itself. We have implemented the algorithm in CPython called robin, based on a fork of the batman code by Kreidberg (2015), which is publicly available 1 .

CONSTRAINING THE RADIUS RATIO FROM DURATIONS OF THE TRANSIT AND OF INGRESS/EGRESS
Given the untoward effect of contamination on measurement of the planet-star radius ratio, we require some additional constraint to derive this ratio in the presence of a heterogenous stellar photosphere or blending with other sources of light. This can be provided by a purely geometrical constraint: the ratio of ingress duration to transit duration. However, to implement this requires a knowledge of (or constraint upon) the impact parameter of the transit, b, which we describe in this section.
The duration of ingress for an eccentric orbit is given approximately by: where θ is the angle between the tangent of the limb of the star and the path of the planet (projected onto the sky plane) and v is the sky velocity of the planet relative to the star during transit, which we take to be constant in this section. The duration of the transit, T , from mid-ingress to mid-egress is given by 1 Open source, available online: https://github.com/ bmorris3/robin These formulae assume that a R * and R p R * , neglecting the curvature of the orbit and the stellar limb (see Winn et al. (2010) or the appendix for more complete relations).
The transit and ingress/egress durations are simply a function of time, and do not depend upon the transit depth, and hence will not be affected by dilution by star spots or flux from a blend, such as a companion star. Thus, a precise measurement of τ and T can in principle give another means of constraining the radius ratio of the planet to the star which will be less affected by blending or a heterogeneous stellar photosphere.
There are two ways to solve for this radius ratio, the first of which is derived from the ratio of ingress to duration of transit, This requires a constraint upon the impact parameter, b. The shape of ingress and egress depends weakly on the impact parameter; however, this is generally too subtle to measure, even with the highest signal-to-noise transits. Note that a given measured value of τ and T imply a maximum value of R p /R * , which is the value this ratio would have if b = 0. This could be useful for placing upper limits on the radii of planets given measured properties of the transit; this constraint is independent of the eccentricity and period of an orbit and independent of the density of the star. If the depth of transit (as parameterized by p 2 1 ) implies a radius ratio that is greater than this maximum value, then this likely implies that both a small impact parameter and that non-transited star spots must be present (or that the planet transits a bright chord) to cause a larger transit depth, or that the planet is significantly oblate (see next section).
The second means of solving for the radius ratio from time-dependent quantities is given by which requires a measurement of the normalized sky velocity, given by: where ρ * is the stellar density, e is the orbital eccentricity, and ω is the longitude of periastron of the star measured along the orbital path from when the star crosses the sky plane away from the observer. 2 Thus, if a/R * is constrained, for example, from asteroseismology or asterodensity profiling (which requires additional transiting planets), and if the eccentricity vector, e = {e cos ω, e sin ω} is constrained, for example, from radial velocity measurements or transit-timing variations, then the transit chord, d/R * = 2 √ 1 − b 2 (in units of R * ), and hence the impact parameter, can be constrained from the transit duration.
In practice, then, the radius ratio may only be determined in cases for which there is a prior on the impact parameter, and generally this is derived from the dependence of the transit duration upon the stellar density and the eccentricity of the planet's orbit, as well as other wellmeasured quantities such as the planet's orbital period. The formulae above are only approximate as they neglect the curvature of the orbit and the curvature of the limb of the star (see the Appendix A for more general expressions), but these relations are generally good approxima-2 For nearly edge-on orbits, transits occur when the true longitude of the star is θ = π/2. The relative velocity of the star and planet is maximum at periastron, and if the time of periastron coincides with the time of transit, then ω = π/2, at which time the velocity is (1 + e)/(1 − e) of the mean orbital velocity.
tions, while the full relation can be accounted for with a full orbital transit model.

Planetary Oblateness
The preceding analysis assumes that a planet is spherical. However, the transit depth and the duration of ingress and egress are also affected by the shape of a planet. If a planet is distorted, for example, due to oblateness, then the area of a planet and the distance a planet moves during ingress and egress are not simply related to the mean radius of the planet. Oblateness can be induced by rotation of a planet, causing a bulging of the equator due to centripetal acceleration Barnes & Fortney (2003), or by thermal structure of a planet, causing a larger scale height at the hotter equator relative to the poles (Dobbs-Dixon et al. 2012). There are several references in the literature which fully consider the transits of oblate planets, and we refer the interested reader to Hui & Seager (2002); Barnes & Fortney (2003); Carter & Winn (2010); Zhu et al. (2014); Biersteker & Schlichting (2017).
In addition to bright or dark active regions and contaminating light sources, planetary oblateness can cause a mismatch between the observed values of p 0 and p 1 . Here we derive a relationship between p 0 , p 1 , and the oblateness of the planet, The projected ellipsoid of an oblate planet with semimajor axis α (in units of R ) and semiminor axis β makes first contact with the star at point P , see Figure 1, which is an angle θ from the semimajor axis of the planet (or the +X axis in the planet-centered coordinate system). The angle from the +X-axis to the normal of the tangent of the ellipse φ = sin −1 (b), where b is the impact parameter.
The definition of the elliptical outline of the planet r is and the slope of a tangent line which meets the ellipse at point P is We assume for simplicity the major axis of the oblate planet is aligned along its direction of motion relative to the star, so α = R eq /R * and β = R pole /R * Now we can solve for p 0 given where ψ is the angle measured up from the +X axis to p 0 . Solving for p 0 in terms of α, β, Note that in the oblate planet case f > 0, p 0 = α for non-zero b. The depth of the transit will be the ratio of the cross-sectional area of the planet to the star, so the general form of the relation between p 0 , p 1 , b and f is: where the last expression is valid for f 1. In the special case where b = 0, this simplifies to implying that p 2 0 = α 2 = R 2 eq /R 2 * and p 2 1 = αβ = R pole R eq /R 2 * .

VALIDATION ON SIMULATED LIGHT CURVES
To estimate the ability to measure the planetstar radius ratio in the presence of a heterogeneous stellar photosphere, we used the reparameterized transit light curve to fit simulated light curves generated with STSP (Hebb et al. 2018 3 ). STSP synthesizes transit light curves of spotted stars by analytically computing the overlap between the planet, starspots, and the star, with limb darkening approximated by concentric circles of constant surface brightness.
We consider two axisymmetric starspot distributions: the first in which the planet's orbital angular momentum is aligned with that of the star in the sky plane, and the planet crosses an active latitude with dark spots, and the second in which active latitudes are present with dark spots, but not transited by the planet. This first case causes shallower transits relative to a uniform star, which is equivalent to contamination by an additional source with > 0, while the second case causes deeper transits, and is equivalent to dilution by a source with negative flux, < 0, which is the deficit of stellar flux caused by the darker spots.

Occultation of an active latitude
We synthesize a transit light curve of a star with a dense band of spots that blankets the stellar equator using STSP with c = 0.7 -see Figure 2. We observe the transit at one second cadence with uncertainties of 48 ppm for each flux measurement, which is similar in scale to the oscillations in flux at the "bottom" of the transit light curve, such that fits to the transit light curve with the standard (Mandel & Agol 2002) model yield a reduced-χ 2 = 0.95. This is a very high signal-to-noise example, even exceeding reasonable S/N estimates for bright targets observed with JWST. The corresponding transit light curve for a planet with impact parameter b = 0 and projected spin-orbit alignment λ = 0 • (i.e. the planet transits the stellar equator) is shallower than the transit of the same system without starspots. The duration of ingress and egress independently encode the radius of the planet. Figure 3 shows the posterior distributions of fits for the parameters p 0 , p 1 for the spotless (left) and spotted (right) light curves, allowing the impact parameter and transit duration T 14 to vary, and holding the limb darkening parameters fixed at their true values. In the spotless light curve case, as expected, we find that p 0 = p 1 . In the case of transiting an active latitude of dark spots, p 0 > p 1 , indicating that the transit depth is being decreased by occultations of dark starspots. The solution for p 0 is consistent with the injected transit radius p 0 = R p /R = 0.05971.
The degeneracy between p 0 , the limb-darkening parameters, and the impact parameter becomes increasingly important as one decreases the S/N of the observations. For example, we inflated the uncertainties on each observation by a factor of three, and we found that the uncertainties on the p 0 , impact parameter, and limb-darkening parameters increased by factors of 5, 2, and 2, respectively. This exercise demonstrates the importance of having independent measurements of the impact parameter or a/R , as discussed in Section 3.

Occultation of a bright latitude
We approximate a bright latitude centered on the stellar equator with STSP by creating active latitudes of spots centered on ±30 • latitude, see Figure 4, while keeping the same planet properties and uncertainties as in Section 4.1. The planet transits a relatively bright region of the stellar surface, so in this example the transit   Figure 2, showing agreement between the radius measurement from the ingress/egress duration, p 0 = R p /R , and the radius measurement from the transit depth p 1 ≈ √ δ where δ is the transit depth. Right: posterior distributions from a fit to the spotted transit light curve in Figure 2, demonstrating p 0 > p 1 , and that the durationdependent radius measurement p 0 recovers the true radius (vertical blue line).
light curve for the spotted star is deeper than expected for the spotless case. Figure 5 shows the posterior distribution of fits for the parameters p 0 , p 1 for the spotless (left) and spotted (right) light curves, allowing the impact parameter and transit duration T 14 to vary, and holding the limb darkening pa-rameters fixed at their true values, as in Section 4.1. In the spotless light curve case, we again find that p 0 = p 1 . In the case of transiting a bright latitude, p 0 < p 1 , indicating that the transit depth is being increased by the spot distribution. The solution for p 0 is   Figure 4, showing agreement between the radius measurement from the ingress/egress duration, p 0 = R p /R , and the radius measurement from the transit depth p 1 ≈ √ δ where δ is the transit depth. Right: posterior distributions from a fit to the spotted transit light curve in Figure 4, demonstrating p 0 < p 1 , and that the durationdependent radius measurement p 0 recovers the true radius (vertical blue line).
again consistent with the injected transit radius We consider further limiting cases with extreme spot distributions in Appendix B, which suggests that plausible limits for selfcontamination are −0.5 ≡ 1 − (p 1 /p 0 ) 2 0.5.

Required Precision
It is important to note that extremely high S/N light curves are required in order to produce meaningful constraints on the difference between p 0 and p 1 , similar to the S/N assumed in the previous two subsections. Descriptions of the uncertainties on p 0 and p 1 in terms of planetary system and observing parameters are detailed in the Appendix. Here we will take the examples of the previous two subsections and ask what photometric precision is necessary to measure significant deviations between p 0 and p 1 . Figure 6 shows the confidence interval for detecting p 0 < p 1 when crossing a dark latitude with the spot map from Section 4.1, as a function of photometric precision on onesecond cadence observations of TRAPPIST-1 g (left); and the confidence interval for detecting p 0 > p 1 when crossing a bright latitude with the spot map from Section 4.2 (right). For the TRAPPIST-1 g system in particular, extremely high S/N precision -10s of ppm uncertainty per flux measurement -would be required to significantly detect the occultation of a dark latitude. The situation is a bit less severe for the occultation of a bright latitude -100s of ppm uncertainty per flux measurement are sufficient to significantly detect p 0 < p 1 in this case. In each case, we note that the required precision is similar to or better than the precision achievable with phase-folded Kepler observations, so we anticipate that any detections of spot distributions with the self-contamination technique from Kepler observations (Section 5.1) will yield only marginal significance, and higher S/N observations, like those of JWST, may be able to achieve the precision necessary to unambiguously separate the posterior distributions of p 0 and p 1 .

Kepler Light Curves
We apply the new parameterization to the highest S/N transit light curves availablethose of the transiting planets Kepler-1, Kepler-2 and Kepler-3 (a.k.a. TrES-2, HAT-P-7 and HAT-P-11). Each of these light curves have short cadence observations throughout the Kepler mission, short period planets, and very bright host stars (K p = 9.84, 9.33, 7.00, respec-tively). As a result, these light curves should yield some of the best constraints on ingress and egress durations available. We also explore some exceptional cases: the highly spotted star Kepler-17, the potentially oblate planet Kepler-39 b, and a planet with a vary high impact parameter, and thus long ingress/egress durations, Kepler-412.
For Kepler-1 through -3, the Kepler Input Catalog estimated contamination in the Kepler aperture is 1%, so a third light source can likely be ruled out as the cause for discrepancies between p 0 and p 1 for these targets (Brown et al. 2011). For Kepler-17, -39 and -412, the contamination is 6%, so discrepancies larger than 6% are required to invoke stellar activity as the culprit for disagreement between p 0 and p 1 .
For each Kepler transit light curve, we normalize the SAP flux by a quadratic fit to the out-of-transit flux. We assume the periods, midtransit epochs, eccentricities and arguments of periastron for each target listed in the NASA Exoplanet Archive.
We fit the transit light curve in Figure 7 for p 0 , p 1 , i o , a/R , holding the quadratic limb darkening parameters fixed at u 1 , u 2 = 0.39256, 0.29064 (Magic et al. 2015). The posterior distributions of p 0 and p 1 are overlapping, indicating that the planet radius inferred from the transit . Posterior probabilities -for p 0 < p 1 given the spot distribution in Section 4.1 (crossing a dark latitude, left), and for p 0 > p 1 given the spot distribution in Section 4.2 (crossing a bright latitude, right) -as a function of the photometric uncertainty on one-second cadence observations of TRAPPIST-1 g. Photometric uncertainties must be smaller than 50 ppm for a significant detection (3σ) of the spot distribution in Section 4.1, or smaller than 300 ppm for a significant detection of the spot distribution in depth is consistent with the radius inferred from the ingress and egress durations (see Figure 7). This measurement aligns with our expectation for a middle-aged, solar-like star such as TrES-2, since sunspots typically cover ∼ 0.03% of the solar surface (Howard et al. 1984). We note that the quadratic limb-darkening parameters of (Magic et al. 2015) appear to be sufficient for describing the transit of this G0 star even at extremely high S/N. Table 1 reports the results from this analysis.

HAT-P-7 (Kepler-2)
HAT-P-7 is a slightly evolved F6 star, orbited by a highly misaligned hot-Jupiter (Pál et al. 2008;Winn et al. 2009a). Asteroseismology indicates that the star is likely in a pole-on configuration (Lund et al. 2014;Benomar et al. 2014). Several authors have noted that there is a ∼ 20 ppm positive residual bump in the transit light curve between ingress and mid-transit, which may be attributed to stellar gravity darkening (Van Eylen et al. 2012;Morris et al. 2013;Masuda 2015).
We fit the phase-folded Kepler light curve of HAT-P-7 for p 0 , p 1 , a/R , i o , t 0 and the four nonlinear limb-darkening parameters. We put a Gaussian prior on a/R with the value inferred from asteroseismology by (Lund et al. 2014), ρ = 0.266 ± 0.008 g cm −3 or a/R = 4.090 ± 0.044. The results are shown in Figure 8. We find that p 0 ≈ p 1 , indicating that the transit chord of the planet is representative of the typical intensity of the star, despite the gravity-darkening signal evident in the transit residuals. This would appear to be a result of the stellar orientation, with its bright pole being occulted between mid-transit and egress, and its dimmer low latitudes being occulted between ingress and mid-transit. The net effect of the gravity darkening on the light curve residuals appears to cancel out since p 0 is consistent with p 1

HAT-P-11 (Kepler-3)
HAT-P-11 is a benchmark system for studying stellar activity with planetary transits. It is a 0.8M star with a rotation period similar to the sun (P rot = 29 d) (Bakos et al. 2010;Sanchis-Ojeda & Winn 2011;Southworth 2011;Deming et al. 2011). In Morris et al. (2017aMorris et al. ( ,b, 2018b we establish qualitative similarities between the Sun and HAT-P-11's starspot distributions and activity cycles. We fit the phase-folded Kepler light curve of HAT-P-11 for p 0 , p 1 , a/R , i o , t 0 and quadratic limb-darkening parameters. We place a Gaussian prior on a/R with the value inferred from asteroseismology by Christensen-Dalsgaard et al. (2010), ρ = 2.5127 ± 0.0009 g cm −3 ⇒ a/R = 14.6950 ± 0.0017, and fix the eccentricity and argument of periastron to the measurements of Winn et al. (2010). The results are shown in Figure 9. We find p 0 ≈ p 1 , indicating that the spot covering fraction is relatively small.
It is interesting to compare the spot covering fractions measured in-and out-of-transit for HAT-P-11 from Morris et al. (2017a). The Kepler light curve shows ∼ 3% rotational modulation each quarter. If we assume the starspots have contrast similar to sunspots, we can use Equations 9-14 of Morris et al. (2017a) to estimate the spot covering fraction within the transit chord from p 0 and p 1 by noticing that the depth of transit will be diluted by starspots with contrast c (where c = 1 means spots with the same intensity as the photosphere and c = 0 means spots that are perfectly dark), and covering fraction f S ,  Figure 9. Left: maximum-likelihood transit model for HAT-P-11 b (red) compared with Kepler short cadence observations (black), and binned residuals of the mode within each bin. Right: posterior distributions for p 0 and p 1 . Note that p 0 ≈ p 1 despite the spots known to be present on HAT-P-11. It would appear that the small spot covering fraction (f S ∼ 3% from Morris et al. 2017a) is insufficient to be measured from the Kepler photometry with this technique.
given p 0 and p 1 , For HAT-P-11, this yields the spot covering fraction in the transit chord f S,transit ≈ 0.01 ± 0.03. We can compare the in-transit spot coverage to the spot covering fraction that one would estimate from the flux deficit method whole star via or f S,asym 0.03 +0.06 −0.01 for HAT-P-11. Detailed spot occultation modeling from Morris et al. (2017a) yields spot covering fractions of 0-10% within the chord of any given transit, broadly consistent with the small covering fraction measured with the self-contamination technique.
We note that due to the significant obliquity of the orbit of the planet relative to the spin of the star, the time-averaged surface brightness contrast of the transit chord may match the stellar disk as the spots traverse about the same fraction of the stellar disk as they do of the transit chord. Consequently, the lack of detection in this case may be hampered by this geometry.

Estimating stellar densities in the absence of asteroseismology
When a star lacks an asteroseismic density measurement, the three-way degeneracy between p 0 , the limb-darkening parameters, and the impact parameter (or equivalently a/R or ρ ) must be broken in order to produce meaningful constraints on p 0 . For the following three planetary systems, Kepler-17, -39 and -412, no asteroseismic density is available in the literature.
We can estimate the stellar densities for stars without asteroseismic measurements by noting the tight correlation between stellar radius and density found in the asteroseismic sample of Huber et al. (2013) -see Figure 10. For bright Kepler targets, Huber et al. (2013) used asteroseismology to calculate stellar densities, and isochrone models to predict stellar radii and masses. We fit a power-law to the measured stellar densities as a function of radius, and find log 10 ρ [g cm −3 ] = (−2.540 ± 0.087) log 10 (R /R ) + (0.12 ± 0.001) ,

Kepler-17
Kepler-17 is an active G2V star with T eff = 5630 ± 100 K and age 3 ± 1.6 Gyr, with a hot Jupiter companion that produces starspot occultations in nearly every transit studied by Désert et al. (2011) and Davenport (2015). Thus Kepler-17's geometry is similar to the toy model in Section 4.1, where a planet occults an active latitude with many starspots, diluting the transit depth, and in which the stellar spin axis is aligned (in this case, to within 15 • ).
We fit the phase-folded Kepler light curve of Kepler-17 for p 0 , p 1 , b, t 0 , T 14 , M , R , and quadratic limb-darkening parameters. We place priors on M and R as described in Section 5.1.4. The results shown in Figure 11 indicate p 0 > p 1 at 93% confidence (∼ 2σ), in agreement with our expectation that starspot occultations dilute the apparent transit depth. We measure p 0 = R p /R = 0.1354 +0.0010 −0.0014 , which is significantly larger than the Désert et al. (2011) estimate (R p /R = 0.13031 +0.00022 −0.00018 ). As we did for HAT-P-11 in Section 5.1.3, we can compare the spot covering fractions measured in-and out-of-transit for Kepler-17. The Kepler light curve shows ∼ 4% rotational modulation each quarter. If we assume the starspots have contrast similar to sunspots, as in Davenport (2015), we can use Equation 23 to estimate the spot covering fraction within the transit chord and find f S,transit ≈ 0.27 +0.12 −0.18 . We can compare the in-transit spot coverage to the spot covering fraction that one would estimate from the flux deficit method, f S,asym 0.13. These two measurements are consistent, and perhaps suggest that the transit chord may be more spotted than the mean photosphere, consistent with the apparent equatorial active latitude deduced from detailed spot modeling by Davenport (2015).

Kepler-39
Kepler-39 is a young F7V host star with T eff = 6350 ± 100 K, and age 0.7 +0.9  Zhu et al. (2014) measure R p /R = 0.0889 ± 0.0006 with their oblate planet model, consistent with our maximum likelihood p 0 = R p /R = 0.090 ± 0.010. We find that p 0 > p 1 with 90% confidence, suggesting weak evidence for either starspot occultations in the transit light curve, or the degenerate signal of planetary oblateness. We find f = 0.25 +0.13 −0.06 , consistent with the oblateness reported by Zhu et al. (2014).
However, we caution that the Kepler light curve shows rotational modulation of the star consistent with starspot coverage (see Figure 13), albeit several times smaller than for Kepler-17. We can estimate the minimum spot covering fraction from the amplitude of the rotational modulation, which gives us a constraint on the asymmetry of the spot distribution as the star rotates, f S,min ≥ (1 − min flux)/(1 − c) where c = 0.7 is the spot contrast for sunspots (Solanki 2003). Thus assuming spots with the contrast of sunspots, f S,min ≥ 2% -a factor of 100 greater than the typical sunspot area coverage, comparable to the spot coverage of HAT-P-11 (Morris et al. 2017a). It is possible that the marginal detection of p 0 > p 1 is attributable to either or both planetary oblateness and occulted starspots.

Kepler-412
Kepler-412 is a G3V host to a hot Jupiter with M p = 0.939 ± 0.085M J , and an apparently inflated radius of 1.325 ± 0.043R J .
We fit the phase-folded Kepler light curve of Kepler-412 for p 0 , p 1 , b, t 0 , T 14 , M , R , and quadratic limb-darkening parameters. We place priors on M and R as described in Section 5.1.4. The results are shown in Figure 14; we find p 0 = R p /R = 0.0926 +0.0095 −0.0089 corresponding to R p = 1.16 ± 0.11R J -consistent with the literature value (R p /R = 0.1058 ± 0.0023, Deleuil et al. 2014). We measure p 0 ≈ p 1 .  Figure 11. Left: maximum-likelihood transit model for Kepler-17 b (red) compared with Kepler short cadence observations (black). Right: Posterior histograms for p 0 , p 1 . For this system we find weak evidence for p 0 > p 1 , consistent with an occultation of a dark region of the stellar surface. This is consistent with the detailed analyses of Désert et al. (2011) and Davenport (2015), which suggest that Kepler-17 occults an active latitude near the stellar equator.  Figure 14. Left: maximum-likelihood transit model for Kepler-412 b (red) compared with Kepler short cadence observations (black). Right: Posterior distributions for p 0 , p 1 and R p /R (equivalent to fixing p 0 = p 1 ). There is insignificant evidence for p 0 < p 1 . See Figure 13 for the out-of-transit rotational modulation of Kepler-412.
As with Kepler-39 in Section 5.1.6, p 0 p 1 could be interpreted as insignificant evidence (∼ 1σ) for dark starspots on the stellar photosphere outside of the transit chord; however, without the ability to further refine the uncertainty on p 0 , we cannot strengthen this claim. The Kepler light curve shows small variations in the stellar intensity that may be rotational modulation due to unocculted starspots (see Figure 13).

K2 Observations of TRAPPIST-1
TRAPPIST-1 is an M8V host to seven Earthsized planets (Gillon et al. 2016(Gillon et al. , 2017Luger et al. 2017;Delrez et al. 2018). Its spots seem to evolve on a timescale similar to the apparent rotation period (Roettenbacher & Kane 2017), the spot and faculae covering fractions might be quite high Zhang et al. 2018), and bright spots might be required to explain the apparent rotational modulation of TRAPPIST-1 (Morris et al. 2018a). Here we analyze the K2/EVEREST short-cadence TRAPPIST-1 light curves of the two innermost TRAPPIST-1 planets for evidence of photospheric inhomogeneities (Luger et al. 2016;Luger 2017).
We fit the transit light curves of TRAPPIST-1 b and c for p 0 , p 1 , T 14 and the quadratic limb darkening coefficients u 1 and u 2 , see Figures 15 and 16. We place Gaussian priors on the impact parameter based on the joint light curve analysis by , with impact parameters for planets b and c: b b = 0.093 +0.011 −0.015 and b c = 0.128 +0.015 −0.020 . We find p 0 ≈ p 1 for both planets, indicating insufficient evidence for a difference in spot coverage inside compared to outside of the overlapping transit chords of planets b and c. Combined with the lack of observed spot occultations, the indication that the star may not be highly spotted is somewhat at odds with the spot covering fraction estimates from Rackham et al. (2018); Zhang et al. (2018), who find a potentially large spot covering fraction, f S = 8 +18 −7 %. As best we can tell from Doppler imaging, fully-convective stars are likely highly spotted, and their spots may be randomly distributed, arranged into active latitudes, or concentrated at the poles (Barnes & Collier Cameron 2001;Donati et al. 2003;Morin et al. 2008). If we continue to measure p 0 ≈ p 1 with further transit observations of TRAPPIST-1, that may suggest that active regions on the star are small, low contrast, and/or uniformly distributed. It will be especially interesting to search for differences between p 0 and p 1 for the outer planets, which span a range of impact parameters, and thus a range of stellar latitudes. Upcoming Spitzer and JWST observations may improve this measurement, as the degeneracy between p 0 and the limb darkening is substantially diminished in the infrared, although the spot contrast is also diminished at longer wavelengths.

Spitzer Light Curves
We analyze two Spitzer light curves of interest. One is the transit of HD 80606 b, which is remarkable -among other reasons -for having a twelve hour transit observed at a very short cadence, providing us with exceptionally wellsampled transit ingress and egress. The other is GJ 1214 -an M dwarf host to a transiting super-Earth.
We fit the Spitzer transit light curve of HD 80606 b for p 0 , p 1 , b, t 0 and T 14 , fixing the quadratic limb-darkening parameters u 1 , u 2 = 0.0866, 0.1071 (Claret et al. 2013). We find p 0 < p 1 with 96% confidence, suggesting that either the planet may be occulting a relatively bright region of the stellar surface, or that we have to place a constraint upon the impact parameter. The transit light curve of HD 80606 shows evidence of at least one starspot occultation; perhaps unocculted starpsots are responsible for the tension between p 0 and p 1 (see Figure 17).
The density of the star inferred by Hébrard et al. (2010) is ρ * = 1.39 ± 0.07 g/cc, which is consistent with the value we infer from the Huber et al. (2013) relation that we derived above of ρ * ,Huber = 1.29 ± 0.11 g/cc. We conclude that if we had placed this density constraint as a prior on the stellar density, we would have obtained a smaller impact parameter, and hence a larger value of p 0 (see Figure 3) consistent with Hébrard et al. (2010).

GJ 1214
GJ 1214 is an M4.5 host to a transiting super-Earth (Charbonneau et al. 2009). The stellar activity of GJ 1214 has been studied extensively for its effect on the transmission spectrum of this potentially water-rich world (Fraine et al. 2013). Rotational modulation of the stellar flux is very limited, suggesting that either the star has either few or small spots, or more likely that the spots are arranged in a nearly-axisymmetric distribution (Berta et al. 2011;Narita et al. 2013). We re-analyze warm Spitzer 4.5 µm photometry of 13 transits of GJ 1214 b by Gillon et al. (2014) to verify that the spot distribution inside the transit chord is similar to the rest of the stellar surface.
We fit the Spitzer transit light curve of GJ 1214 for p 0 , p 1 , b, T 14 and fit for the quadratic limb darkening parameters q 1 and q 2 (Kipping 2013), see Figure 18. The resulting constraint on p 0 is very weak and is consistent with p 1 , suggesting an isotropic spot distribution, consistent with the conclusions of Berta et al. (2011) and Narita et al. (2013).

DISCUSSION
We have investigated the prospects for robust measurement the radius ratios of planets using the ratio of the ingress duration to total transit duration, which should be less affected by contamination by starspots than the transit depth measurements. We have then applied this to a sample of ten stars to look for differences be-tween the radius ratio measured from ingress (p 0 ) versus that measured from transit depth (p 1 ).
Given the probability measurements of p 0 < p 1 for the ten stars in our study, we can ask the question of whether there is evidence that the probability sample is drawn from a uniform probability distribution. We sort the probabilities of p 0 < p 1 for the ten stars in our study, and plot these versus a uniform probability distribution. Figure 19 shows the cumulative probability distribution of the sample versus the measured probabilities. If the distributions of p 0 and p 1 are statistically consistent for the entire sample, then we expect that ten measured probabilities of p 0 < p 1 will be consistent with being drawn from a uniform probability distribution. We apply the K-S test to the sample, finding a maximum distance of −0.2 for a sample size of 10, giving a probability that this is drawn from a uniform distribution of 28%. We view this as indicating that the distribution of p 0 < p 1 probabilities is likely drawn from a uniform distribution; i.e. that the sample as a whole is consistent with p 0 = p 1 for all planets.
Similar analyses have been conducted for the single systems CoRoT-2 (Bruno et al. 2016) and CoRoT-7 (Barros et al. 2014), and more general discussions appear in Oshagh et al. (2013) and Csizmadia et al. (2013). The main differences between this work and others is that we explicitly use the ingress and egress durations to measure the true planet radii in the presence of bright or dark regions on the stellar surface.
This work echoes Csizmadia et al. (2013) in particular, we found that fixing the limbdarkening parameters to their theoretical values for several systems introduced artificial constraints on the impact parameter and p 1 , occasionally producing apparent discrepancies where p 0 = p 1 , which become consistent once the limb-darkening parameters are allowed to float. We encourage users of this light curve Target  Table 1. Maximum-likelihood parameters for p 0 , p 1 , the self-contamination parameter = 1 − (p 1 /p 0 ) 2 , and the confidence interval (CI) for each detection of p 0 relative to p 1 .  Figure 19. Cumulative probability distribution for the probability that p 0 < p 1 for the ten stars in our sample (blue points). The light grey lines show 3000 Monte Carlo simulations of probabilities drawn from a uniform distribution between zero and one. The green curve shows a uniform probability distribution.
parameterization to fit for the limb-darkening parameters whenever possible, though we note the exceptionally good fit to TrES-2 using the quadratic limb-darkening parameters of (Magic et al. 2015).
We also note the importance of asteroseismic constraints on the orbital parameters for several systems. The asteroseismic stellar density provides us with an independent measurement of a/R , which is otherwise degenerate with i o , p 0 and the limb-darkening parameters. Bright TESS targets may also be amenable to asteroseismic density measurements, and therefore good candidates for this self-contamination analysis. Alternatively the stellar density can be measured to high precision in multi-planet systems from the planet periods and transit durations, which will be valuable for dim stars like TRAPPIST-1.
Diminished limb-darkening in the infrared promises that Spitzer and JWST observations of transiting exoplanets will be interesting subjects for analysis with the self-contamination analysis. Ground-based follow-up may also be a fruitful means of detecting the affects of stellar activity on transiting exoplanet light curves, especially with the high precision and time resolution enabled by holographic diffusers (see e.g. Stefansson et al. 2017;Morris et al. 2018b).

CONCLUSIONS
We have presented a reparameterization of the transit light curve of Mandel & Agol (2002) which splits the original "p" parameter into two parameters: p 0 = R p /R which defines the planet radius, and p 1 ≈ √ δ which defines the transit depth. This parameterization allows the transit model to account for significant contamination by bright or dark active regions on the stellar surface, or significant planetary oblateness in the planet's direction of motion. The resulting constraint on p 0 = R p /R is more accurate in the presence of stellar activity (and usually less precise).
We fit the light curves of several transiting planets to study the photospheric inhomogeneities of their host stars, see Table 1. We report p 0 ≈ p 1 (no significant detection of selfcontamination) for Kepler light curves of TrES-2, HAT-P-7, TRAPPIST-1 and Spitzer observations of GJ 1214. We find that the uncertainties on p 0 are typically an order of magnitude larger than p 1 , consistent with estimates of the noise (Appendix C). We find little evidence for transit depth dilution due to occulted starspots (p 0 > p 1 ) for the well-studied spotted pair of host stars HAT-P-11 and Kepler-17, and we find weak evidence that the transiting planets of Kepler-412 and HD 80606 likely occult relatively bright regions of the stellar photosphere (p 0 < p 1 ). We can recover the reported oblateness of Kepler-39 b using this parameterization, though its detection is not statistically significant. We note that in general the oblateness is likely degenerate with self-contamination -so the unusually large inferred oblateness may have an alternate explanation as being due to an inhomogeneous photosphere.
In the best cases we have studied in this paper, HAT-P-11b, Kepler-17b, and TrES-2b, we obtain sensitivity to self-contamination, at the 1-3% level (Table 1), which is close enough to the expected levels of contamination that we were unable to achieve a definitive detection. This ratio of precisions on p 0 range from 5 − 25 times the uncertainties on p 1 (Table 1) The main advantage, then, in measuring p 0 is to obtain the radius ratio of the planet to the star when the unknown contamination of the stellar flux is significant. The other circumstance which may favor measuring p 0 is when there is significant red noise which has a larger amplitude on the timescale of the transit than on the timescale of ingress/egress. Granulation noise tends to be strongly red, and thus will affect the measurement of the transit depth more strongly than the measurement of the ingress/egress duration.
With the James Webb Space Telescope (JWST) we expect to obtain much higher precision on p 0 due to the larger collecting area. As an example, the Trappist-1 system has ingress/egress durations of ≈ 2 − 5 minutes, and when observed with NIRSPEC, we expect precisions of 2 − 5 seconds for the precision of these measurements. This will enable a measurement of p 0 comparable to the current precisions obtained with Spitzer for which multiple transits have been observed for each planet. Thus, we expect that JWST will yield constraints upon p 0 which will give accurate radius ratios for these planets to diagnose the contamination by starspots. We find our simulations with extremely high signal-to-noise ( §4) demonstrate that the radius-ratio of the planet measured from p 0 can be recovered more accurately than the radius ratio measured from p 1 in the presence of star spot contamination, since p 1 is biased by the presence of star spots (although p 1 is measured more precisely).
We have shown that p 0 is degenerate with the limb-darkening parameters. As a result, this technique will work better at red wavelengths, where limb-darkening is less severe, than in the blue. However, the contrast of starspots against the stellar photosphere diminishes as one observes at longer wavelengths as well, so this technique might be best-suited to observations in the red optical, where the spot contrast may be significant but the limb-darkening is weaker.
This technique is complementary to other starspot measurement techniques -it is sensitive to time-independent spots, though it provides a weaker signal than techniques like the flux deficit. The flux deficit technique -which measures variations in flux as a star rotates -has been applied to vast numbers of Kepler stars (Walkowicz & Basri 2013;Notsu et al. 2013;Mathur et al. 2014). Or alternatively, spot distributions can be revealed by detailed modeling of individual spots throughout rotational modulation   For a transiting exoplanet in a circular orbit, the transit duration from first through fourth contact is and the duration from second through third contact is where p 0 = R p /R , P is the orbital period, b is the impact parameter, i is the planet's orbital inclination, and a/R is the scaled orbital semi-major axis (Winn 2010). We can invert these relations to find the planet-to-stellar radius ratio as a function of the transit durations, We can get a sense for the scaling of the terms in this equation by simplifying to the case where i = 90 • , b = 0 and T 14 , T 23 P using the small angle approximation, The constraint on the planetary radius (p 0 = R p /R in our reparameterization) primarily comes from the difference between the durations of first-through-fourth contact and second-through-third contact.
B. LIMITING CASES: EXTREME SPOT DISTRIBUTIONS In Sections 4.1 and 4.2 we assumed that the starspots had the same contrast as sunspots, c = 0.7. In this section, we revisit those toy-model analyses with more extreme spot contrasts, by setting the contrast to that of sunspot umbrae, c = 0.2.
First we simulate a dense band of spots with contrast c = 0.2 centered on the stellar equator -see Figure 20. The self-contamination parameter = 1 − (p 1 /p 0 ) 2 = −0.598 +0.014 −0.032 in this limiting case. Next we simulate a pair of polar spots that extend almost all the way from the pole to the edge of the transit chord, with contrast c = 0.2 -see Figure 21. The self-contamination parameter = 0.419 +0.011 −0.015 in this case.
Thus one might expect to vary roughly on −0.5 0.5 for extreme self-contamination due to occulted or unocculted starspots.

C. UNCERTAINTIES
Using the estimates by Carter et al. (2008) for the uncertainties on the ingress and total duration given Fisher information analysis of a piecewise-linear transit model, we can estimate the uncertainties with the transit of the same system without spots (gray dashed curve). At lower S/N, or for different spot geometries, the bottom of the spotted transit might simply appear flat and relatively shallow, but the ingress and egress durations are the same for both light curves, allowing us to recover the true planet radius from timing, independent of the transit depth which is affected by starspots.
on p 0 and p 1 . For flux observations with independent Gaussian uncertainties σ, for a transit with depth δ, ingress duration τ , mid-ingress to mid-egress duration T , defined as where τ 0 is with uncertainties on τ and T given by where Γ is the sampling rate (if N photometric measurements are taken uniformly over a total observing time T 0 , then Γ = N/T 0 ), so Γ −1 is the time between the start of successive exposures. Rearranging, we find and therefore the uncertainty, σ p 0 , is given by in the limit of τ T , where η = T /(T 0 − T − τ ) is the ratio of the duration of photometric measurements in-transit to out-of-transit (η → 0 for perfectly known out-of-transit flux). From this equation we see conclusions we discussed elsewhere in this work, for example, if τ 0 is known imprecisely by lack of prior constraints on the stellar density or transit impact parameter, then p 0 will have large uncertainties.