A 38 Million Year Old Neptune-sized Planet in the Kepler Field

Kepler 1627A is a G8V star previously known to host a 3.8 R ⊕ planet on a 7.2 day orbit. The star was observed by the Kepler space telescope because it is nearby (d = 329 pc) and it resembles the Sun. Here, we show using Gaia kinematics, TESS stellar rotation periods, and spectroscopic lithium abundances that Kepler 1627 is a member of the 38−5+6 Myr old δ Lyr cluster. To our knowledge, this makes Kepler 1627Ab the youngest planet with a precise age yet found by the prime Kepler mission. The Kepler photometry shows two peculiarities: the average transit profile is asymmetric, and the individual transit times might be correlated with the local light-curve slope. We discuss possible explanations for each anomaly. More importantly, the δ Lyr cluster is one of ∼103 coeval groups whose properties have been clarified by Gaia. Many other exoplanet hosts are candidate members of these clusters; their ages can be verified with the trifecta of Gaia, TESS, and ground-based spectroscopy.

1. INTRODUCTION While thousands of exoplanets have been discovered orbiting nearby stars, the vast majority of them are several billion years old.This makes it difficult to test origin theories for the different families of planets, since many evolutionary processes are expected to operate on timescales of less than 100 million years.
For instance, the "mini-Neptunes", thought to be made of metal cores, silicate mantles (Kite et al. 2020), and extended hydrogen-dominated atmospheres, are expected to shrink in size by factors of several over their first 10 8 years.Specifically, in the models of Owen & Wu (2016) and Owen (2020), the ≈5 M ⊕ planets start with sizes of 4 -12 R ⊕ shortly after the time of disk dispersal ( 10 7 years), and shrink to sizes of 2 -4 R ⊕ by 10 8 years.While the majority of this change is expected to occur within the first few million years after the disk disperses (Ikoma & Hori 2012), stellar irradiation and internal heat can also power gradual outflows which, if strong enough, can deplete or entirely strip the atmosphere (Lopez et al. 2012;Owen & Wu 2013;Ginzburg et al. 2018).The photoevaporative and core-powered outflows are thought to persist for ≈10 8 to ≈10 9 years, though the details depend on the planetary masses, the irradiation environments, and the initial atmospheric mass fractions (Owen & Wu 2017;Gupta & Schlichting 2020;Rogers & Owen 2021;King & Wheatley 2021).Discovering young planets, measuring their masses, and detecting their atmospheric outflows are key steps toward testing this paradigm, which is often invoked to explain the observed radius distribution of mature exoplanets (Fulton et al. 2017;Van Eylen et al. 2018).
The K2 and TESS missions have now enabled the detection of about ten close-in planets younger than 100 million years, all smaller than Jupiter (Mann et al. 2016;David et al. 2016David et al. , 2019;;Newton et al. 2019;Bouma et al. 2020;Plavchan et al. 2020;Rizzuto et al. 2020;Martioli et al. 2021).The Kepler mission however has not yielded any planets with precise ages below one gigayear (Meibom et al. 2013).The reason is that during the prime Kepler mission (2009)(2010)(2011)(2012)(2013), only four open clusters were known in the Kepler field, with ages spanning 0.7 Gyr to 9 Gyr (Meibom et al. 2011).Though isochronal, gyrochronal, and lithium-based analyses suggest that younger Kepler planets do exist (Walkowicz & Basri 2013;Berger et al. 2018;David et al. 2021), accurate and precise age measurements typically require an ensemble of stars.Fortunately, recent analyses of the Gaia data have greatly expanded our knowledge of cluster memberships (e.g., Cantat-Gaudin et al. 2018;Zari et al. 2018;Kounkel & Covey 2019;Meingast et al. 2021;Kerr et al. 2021).As part of our Cluster Difference Imaging Photometric Survey (CDIPS, Bouma et al. 2019), we concatenated the available analyses from the literature, which yielded a list of candidate young and agedated stars (see Appendix A).
Matching our young star list against stars observed by Kepler revealed that Kepler observed a portion of the δ Lyr cluster Theia 73).More specifically, a clustering analysis of the Gaia data by Kounkel & Covey (2019) reported that Kepler 1627 (KIC 6184894; KOI 5245) is a δ Lyr cluster member.Given the previous statistical validation of the close-in Neptune-sized planet Kepler 1627b (Tenenbaum et al. 2012;Morton et al. 2016;Thompson et al. 2018), we begin by scrutinizing the properties of the cluster (Section 2).We find that the δ Lyr cluster is 38 +6 −5 Myr old, and in Section 3 show that Kepler 1627 is both a binary and also a member of the cluster.Focusing on the planet (Section 4), we confirm that despite the existence of the previously unreported M3V companion, hereafter Kepler 1627B, the planet orbits the G-dwarf primary, Kepler 1627A.We also analyze an asymmetry in the average transit profile, and a possible correlation between the individual transit times and the local light curve slope.We conclude by discussing broader implications for our ability to age-date a larger sample of planets (Section 5).
2. THE CLUSTER To measure the age of the δ Lyr cluster, we first selected a set of candidate cluster members (Section 2.1), and then analyzed these stars using a combination of the isochronal and gyrochronal techniques (Section 2.2).

Selecting Cluster Members
Kounkel & Covey (2019) applied an unsupervised clustering algorithm to Gaia DR2 onsky positions, proper motions, and parallaxes for stars within the nearest kiloparsec.For the δ Lyr cluster (Theia 73), they reported 3,071 candidate members.We matched these stars against the latest Gaia EDR3 observations using the dr2_neighbourhood table from the ESA archive, taking the stars closest in proper motion and epoch-corrected angular distance as the presumed match (Gaia Collaboration et al. 2021a).For plotting purposes, we focused only on the stars with parallax signal-to-noise exceeding 20.We calculated the tangential velocities for each of these stars relative to Kepler 1627 (∆µ * ) by subtracting the observed proper motion from what the proper motion at each star's position would be if it were co-moving with Kepler 1627.
Figure 1 shows that the reported cluster members (gray and black points) extend over a much larger volume in both physical and kinematic space than the cluster previously identified by Stephenson (1959) and later corroborated by Eggen (1968).While the non-uniform "clumps" of stars might comprise a bona fide cluster of identically-aged stars, they could also be heavily contaminated by field stars.One reason to suspect this is that the spread in tangential velocities exceeds typical limits for kinematic coherence (Meingast et al. 2021).We therefore considered stars only in the immediate kinematic and spatial vicinity of Kepler 1627 as candidate cluster members.We performed this selection cut manually, by drawing lassos with the interactive glue visualization tool (Beaumont et al. 2014) in the four projections shown in Figure 1.The overlap between the Kepler field and the resulting candidate cluster members is shown in Fig- ure 2. While this method will include some field interlopers in the "cluster star" sample, and viceversa, it should suffice for our aim of verifying the existence of the cluster in the vicinity of Kepler 1627.

Color-Absolute Magnitude Diagram
We measured the isochrone age using an empirical approach.The upper left panel of We also compared against µ-Tau (62 ± 7 Myr; Gagné et al. 2020) and the Upper-Centaurus-Lupus (UCL) component of the Sco OB2 association (≈16 Myr; Pecaut & Mamajek 2016).We adopted the UCL members from Damiani et al. (2019).For visual clarity, the latter two clusters are not shown in Figure 3.We cleaned the membership lists following the data filtering criteria from Gaia Collaboration et al. (2018a, Appendix B), except that we weakened the parallax precision requirement to ϖ/σ ϖ > 5.This also involved cuts on the photometric signal to noise ratio, the number of visibility periods used, the astrometric χ 2 of the singlesource solution, and the G BP − G RP color excess factor.These filters were designed to include genuine binaries while omitting instrumental artifacts.
To correct for extinction, we queried the 3dimensional maps of Capitanio et al. (2017) and Lallement et al. (2018) 1 , and applied the extinction coefficients k X ≡ A X /A 0 computed by Gaia Collaboration et al. (2018a) assuming that A 0 = 3.1E(B − V ).For UCL, IC 2602, the Pleiades, and the δ Lyr cluster, this procedure yielded a respective mean and standard deviation for the reddening of E(B −V ) = {0.084± 0.041, 0.020 ± 0.003, 0.045 ± 0.008, 0.032 ± 0.006}.These values are within a factor of two of previously reported values in the literature (Pecaut & Mamajek 2016;Gaia Collaboration et al. 2018a;Kounkel & Covey 2019;Bossini et al. 2019), and are all small enough that the choice of whether to use them vs. other extinction estimates does not affect our primary conclusions.
Figure 3 shows that the δ Lyr cluster and IC 2602 overlap, and therefore are approximately the same age.The pre-main-sequence M dwarfs of a fixed color in the δ Lyr cluster are more luminous than those of the Pleiades, and were also seen to be less luminous than those in UCL.To turn this heuristic interpolation into a quantitative age measurement, we used the empirical method developed by Gagné et al. (2020).In brief, we fitted the pre-main-sequence loci of a set of reference clusters, and the locus of the target δ Lyr cluster was then modeled as a piecewise linear combination of these reference clusters.For our reference clusters, we used UCL, IC 2602, and the Pleiades.We removed binaries by requiring RUWE < 1.3, radial_velocity_error below the 80 th percentile of each cluster's distribution, and excluded stars that were obvious photometric binaries in the CAMD. 2 We then passed a moving box average and standard deviation across the CAMD in 0.10 mag bins, fitted a univariate spline to the binned values, and assembled a piecewise grid of isochrones spanning the ages between UCL and the Pleiades using Equation 6 from Gagné et al. (2020).To derive a probability distribution function for the age of δ Lyr cluster, we then assumed a Gaussian likelihood that treated the interpolated isochrones as the "model" and the δ Lyr cluster's isochrone as the "data" (Equation 7 from Gagné et al. 2020).The cluster's age and its statistical uncertainty are then quoted as the the mean and standard deviation of this age posterior.
The ages returned by this procedure depend on the ages assumed for each reference cluster.We adopted a 115 Myr age for the Pleiades (Dahm 2015), and a 16 Myr age for UCL (Pecaut & Mamajek 2016).The age of IC 2602 however is the most important ingredient, since it receives the most weight in the interpolation.Plausible ages for IC 2602 span 30 Myr to 46 Myr, with older ages being preferred by the lithium-depletion-boundary (LDB) measurements (Dobbie et al. 2010;Randich et al. 2018) and younger ages by the mainsequence turn-off (Stauffer et al. 1997;David & Hillenbrand 2015;Bossini et al. 2019).If we were to adopt the 30 Myr age for IC 2602, then the δ Lyr cluster would be 31 Of the 3,071 candidate δ Lyr cluster members reported by Kounkel & Covey (2019), 924 stars were amenable to rotation period measurements (G < 17 and (G BP − G RP ) 0 > 0.5) using the TESS full frame image data.As a matter of scope, we restricted our attention to the 391 stars discussed in Section 2.1 in the spatial and kinematic proximity of Kepler 1627.We extracted light curves from the TESS images using the nearest pixel to each star, and regressed them against systematics with the causal pixel model implemented in the unpopular package (Hattori et al. 2021).We then measured candidate rotation periods using a Lomb-Scargle periodogram (Lomb 1976;Scargle 1982;Astropy Collaboration et al. 2018).To enable cuts on crowding, we queried the Gaia source catalog for stars within a 21.0 radius of the target star (a radius of 1 TESS pixel).Within this radius, we recorded the number of stars with greater brightness than the target star, and with brightness within 1.25 TESS magnitudes of the target star.We then cleaned the candidate TESS rotation period measurements through a combination of automated and manual steps.First, to validate the TESS rotation periods, we compared against 28 stars from McQuillan et al. (2014) that were also observed by Kepler.Of the 23 stars with Kepler periods below 10 days, 21 of the TESS periods agreed with the Kepler rotation periods; the other 2 were measured at the double-period harmonic.Of the remaining 5 stars with Kepler rotation periods above 10 days, none were correctly recovered by TESS, and 3 were near the half-period harmonic.We were therefore wary of any TESS rotation periods exceeding 10 days, and used the Kepler rotation periods whenever possible.For the remaining stars with only TESS data, we focused only on the stars for which no companions were known with a brightness exceeding one-tenth of the target star in a 21.0 radius.There were 192 stars that met these crowding requirements, and that had TESS data available.For plotting purposes we then im-  (Rebull et al. 2016;Douglas et al. 2017).Most candidate δ Lyr cluster members are gyrochronally younger than the Pleiades; outliers are probably field interlopers.posed a selection based on the strength of the signal itself: we required the Lomb Scargle power to exceed 0.2, and the period to be below 15 days.
The lower panel of Figure 3 shows the resulting 145 stars.The majority of these stars fall below the "slow sequence" of the Pleiades, consistent with a gyrochronal age for the δ Lyr cluster below 100 Myr.In fact, the rotation-color distributions of other 30 Myr to 50 Myr clusters (e.g., IC 2602 and IC 2391) are indistinguishable (Douglas et al. 2021).Approximately 10 of the δ Lyr cluster stars appear as outliers above the "slow sequence".Assuming that they are all false positives (i.e., field interlopers), our rotation period detection fraction would be 135/192 ≈ 70%.Although some of these outlier stars might be unresolved F+K binaries that are in the cluster (Stauffer et al. 2016), assuming that they are field contaminants provides a more secure lower bound of the rotation period detection fraction.A final possible confounding factor -binarity -is known to affect the "fast sequence" of stars beneath the slow sequence (Meibom et al. 2007;Gillen et al. 2020;Bouma et al. 2021).We do not expect it to change the central conclusion regarding the cluster's age.Based on the spatial and kinematic association of Kepler 1627 with the δ Lyr cluster, and the assumption that the planet formed shortly after the star, it seems likely that Kepler 1627 is the same age as the cluster.There are two consistency checks on whether this is true: rotation and lithium.Based on the Kepler light curve, the rotation period is 2.642 ± 0.042 days, where the quoted uncertainty is based on the scatter in rotation periods measured from each individual Kepler quarter.This is consistent with comparable cluster members (Figure 3).
To infer the amount of Li I from the 6708 Å doublet (e.g., Soderblom et al. 2014), we acquired an iodine-free spectrum from Keck/HIRES on the night of 2021 March 26 using the standard setup and reduction techniques of the California Planet Survey (Howard et al. 2010).Following the equivalent width measurement procedure described by Bouma et al. (2021), we find EW Li = 233 +5 −7 mÅ.This value does not correct for the Fe I blend at 6707.44Å.Nonetheless, given the stellar effective temperature (Table 1), this measurement is in agreement with expectations for a ≈ 40 Myr star (e.g., as measured in IC 2602 by Randich et al. 2018).It is also larger than any lithium equivalent widths measured by Berger et al. (2018) in their analysis of 1,301 Kepler-star spectra.

Stellar Properties
The adopted stellar parameters are listed in Table 1.The stellar mass, radius, and effective temperature are found by interpolating against a 38 Myr MIST isochrone (Choi et al. 2016).The statistical uncertainties are propagated from the absolute magnitude (mostly originating from the parallax uncertainty) and the color; the systematic uncertainties are taken to be the difference between the PARSEC (Bressan et al. 2012) and MIST isochrones.Reported uncertainties are a quadrature sum of the statistical and systematic components.As a consistency check, we analyzed the aforementioned Keck/HIRES spectrum from the night of 2021 March 26 using a combination of SpecMatch-Emp for stellar properties, and SpecMatch-Synth for v sin i (Yee et al. 2017).This procedure yielded T eff = 5498 ± 100 K, log g = 4.6 ± 0.1, [Fe/H] = 0.15 ± 0.10 from SpecMatch-Emp, and v sin i = 18.9 ± 1.0 from SpecMatch-Synth.These values are within the 1-σ uncertainties of our adopted values from the isochrone interpolation.

Kepler 1627B
We first noted the presence of a close neighbor in the Kepler 1627 system on 2015 July 22 when we acquired adaptive optics imaging using the NIRC2 imager on Keck-II.We used the narrow camera (FOV = 10.2 ) to obtain 8 images in the K filter (λ = 2.12 µm) with a total exposure time of 160 s.We analyzed these data following Kraus et al. (2016), which entailed using PSF-fitting to measure the separation, position angle, and contrast of the candidate companion.The best-fitting empirical PSF template was identified from among the near-contemporaneous observations of single stars in the same filter.The mean values inferred from the 8 images are reported in Table 1.To estimate the detection limits, we analyzed the residuals after subtracting the empirical PSF template.Within each residual image, the flux was measured through 40 mas apertures centered on every pixel, and then the noise as a function of radius was estimated from the RMS within concentric rings.Finally, the detection limits were estimated from the strehl-weighted sum of the detection significances in the image stack, and we adopted the 6-σ threshold as the detection limit for ruling out additional companions.We also observed Kepler 1627 on Gemini-North using the 'Alopeke speckle imager on 2021 June 24.'Alopeke is a dual-channel speckle interferometer that uses narrow-band filters centered at 0.83 µm and 0.56 µm.We acquired three sets of 1000 × 60 msec exposures during good seeing (0.45 ), and used the autocorrelation function of these images to reconstruct a single image and 5σ detection limits (see Howell et al. 2011).This procedure yielded a detection of the companion in the 0.83 µm notch filter, but not the 0.56 µm filter.
The measured projected separation and magnitude difference are given in Table 1. Figure 4 summarizes the results of the highresolution imaging.The Gaia EDR3 parallax for the primary implies a projected separation of 53 ± 4 AU, assuming the companion is bound.Although the companion is unresolved in the Gaia source catalog (there are no comoving, codistant candidate companions brighter than G < 20.5 mag within ρ < 120 ), its existence was also suggested by the primary star's large RUWE relative to other members of the δ Lyr cluster (RUWE≈2.9;roughly the 98 th percentile of the cluster's distribution).Based on the apparent separation, the binary orbital period is of order hundreds of years.The large RUWE is therefore more likely to be caused by a PSF-mismatch skewing the Gaia centroiding during successive scans, rather than true astrometric motion.Regardless, given the low geometric probability that a companion imaged at ρ ≈ 0. 16 is a chance line-of-sight companion, we proceed under the assumption that the companion is bound, and that Kepler 1627 is a binary.Given the distance and age, the models of Baraffe et al. (2015) imply a companion mass of M B ≈ 0.30 M and companion temperature of T eff,B ≈ 3408 K.The corresponding spectral type is roughly M3V (Pecaut & Mamajek 2013).These models combined with the NIRC2 contrast limits imply physical limits on tertiary companions of M ter < 50M Jup at ρ = 50 AU, M ter < 20M Jup at ρ = 100 AU, and M ter < 10M Jup at ρ = 330 AU.

THE PLANET 4.1. Kepler Light Curve
The Kepler space telescope observed Kepler 1627 at a 30-minute cadence from 2009 May 2 until 2013 April 8. Data gaps during quarters 4, 9, and 13 led to an average duty cycle over the 3.9 year interval of 67%.Kepler 1627 was also observed at 1-minute cadence from 2012 Oct 5 until 2013 Jan 11.The top panel of Figure 5 shows a portion of the 30-minute cadence PDCSAP light curve.Nonastrophysical variability has been removed using the methods discussed by Smith et al. (2017); the default optimal aperture was assumed (Smith et al. 2016).Cadences with non-zero quality flags (9% of the data) have been omitted.The resulting photometry is dominated by a quasi-periodic starspot signal with a peak-to-peak amplitude that varies between 2% and 8%.Given that the secondary companion's brightness in the Kepler band is ≈1.5% that of the primary, source confusion for the rotation signal is not expected to be an issue.Previous analyses have identified and characterized the smaller transit signal (Tenenbaum et al. 2012;Thompson et al. 2018), validated its planetary nature (Morton et al. 2016), and even searched the system for transit timing variations (Holczer et al. 2016).Nonetheless, since the cluster membership provides us with more precise stellar parameters than those previously available, we opted to reanalyze the light curve.

Transit and Stellar Variability Model
We fitted the Kepler long cadence time series with a model that simultaneously included the planetary transit and the stellar variability.The stellar variability was modeled with the RotationTerm Gaussian Process kernel in exoplanet (Foreman-Mackey et al. 2020).This kernel assumes that the variability is generated by a mixture of two damped simple harmonic oscillators with characteristic frequencies set by 1/P rot and its first harmonic.We additionally included a jitter term to inflate the flux uncertainties in a manner that accounted for otherwise unmodeled excess white noise, and let the eccentricity float.For the limb-darkening, we assumed a quadratic law, and sampled using the uninformative prior suggested by Kipping (2013).
Our model therefore included 10 free parameters for the transit ({P,t 0 , δ, b, u 1 , u 2 , R , log g, e, ω}), 2 free parameters for the light curve normalization and a white noise jitter ({ f , σ f }), and 5 hyperparameters for the GP ({σ rot , P rot , Q 0 , dQ, f }).We also considered including an additive SHOTerm kernel to account for stochastic noise, but found that this did not affect the results, and so opted for the simpler GP kernel.We fitted the models using PyMC3 (Salvatier et al. 2016; Theano Development Team 2016), and accounted for the finite integration time of each exposure in the numerical integration when evaluating the model light curve (see Kipping 2010).We assumed a Gaussian likelihood, and after initializing each model with the parameters of the maximum a posteriori model, we sampled using PyMC3's gradient-based No-U-Turn Sampler (Hoffman & Gelman 2014) in the bases indicated in Table 2.We used R as our convergence diagnostic (Gelman & Rubin 1992).
Figure 5 shows the resulting best-fit model in orange (top) and purple (bottom).The model parameters and their uncertainties, given in Table 2, are broadly consistent with a Neptune-sized planet (3.82 ± 0.16 R ⊕ ) on a close-in circular4 orbit around a G8V host star (0.88 ± 0.02R ).This best-fit planet size is consistent with those previously reported by Morton et al. (2016) and Berger et al. (2018), and corrects for the small amount of flux dilution from Kepler 1627B.

Transit Asymmetry
The transit fit however is not perfect: the lower panels of Figure 5  Best fit 1 2 Figure 6.Weak evidence for a prograde orbit of Kepler 1627 Ab.The time of each Kepler transit was measured, along with the local slope of the light curve.The two quantities might be anti-correlated (≈2-σ), which could be caused by starspot crossings during the first (second) half of transit inducing a positive (negative) TTV, provided that the orbit is prograde (Mazeh et al. 2015).The units along the abscissa can be understood by considering that the stellar flux changes by ∼60 ppt per half rotation period (∼1.3 days).
the second half.The semi-amplitude of this deviation is ≈ 50 ppm, which represents a ≈ 3% distortion of the transit depth (δ = 1759 ± 62 ppm).Note that although this asymmetry is within the 2σ model uncertainties, the model has a jitter term that grows to account for excess white noise in the flux.The significance of the asymmetry is therefore best assessed in comparison against the intrinsic out-of-transit scatter in the data (≈ 16 ppm), not the model uncertainties.The lower right panel of Figure 5 demonstrates that the scatter during transit is higher than during all other phases of the planet's orbit.
To determine whether the asymmetry could be a systematic caused by our stellar variability model, we explored an alternative approach in which we isolated each transit window, locally fitted out polynomial trends, and then binned all the observed transits; the asymmetry was still present at a comparable amplitude.Appendix B describes a more detailed analysis, which finds that the asymmetry also seems to be robust to different methods of data binning in time and by local light curve slope.Possible astrophysical explanations are discussed in Section 5.

Transit Timing and the Local Slope
The previous analysis by Holczer et al. (2016) did not find any significant long-term transit timing or duration variations (TTVs or TDVs) for Kepler 1627.Quantitatively, the mean and standard deviation of the TTVs and TDVs they mea-sured were −1.1±13.8min and −3.3±22.1 min.In an earlier analysis however, Holczer et al. (2015) studied correlations between TTVs and local light curve slopes, and for Kepler 1627 found a weak correlation of −29 ± 13 min day −1 between the two quantities.Given the possible connection between such correlations and the unresolved starspot crossings that we expect to be present in the Kepler 1627 light curve (Mazeh et al. 2015), we opted to reexamine the individual transit times.We therefore isolated each of the 144 observed transits to within ±4.5 hr of each transit, and fitted each window with both i) a local second or fourth-order polynomial baseline plus the transit, and ii) a local linear trend plus the transit.We let the mid-time of each transit float, and then calculated the residual between the measured midtime and that of a periodic orbit.This residual, the transit timing variation, is plotted in Figure 6 against the local linear slope assuming the fourthorder polynomial baseline.The slope of −21 ± 10 min day −1 is similar to that found by Holczer et al. (2015).The best-fit line yields χ 2 = 306.1, with n = 140 data points.An alternative model of just a flat line yields χ 2 = 315.6.The difference in the Bayesian information criterion between the two models is BIC flat − BIC line = 4.5, which corresponds to a Bayes factor of ≈9.4.According to the usual Kass & Raftery (1995) criteria, this is "positive" evidence for the model with a finite slope.We view it as suggestive at best, particularly given the poor reduced χ 2 .A separate concern we had in this analysis was whether our transit fitting procedure might induce spurious correlations between the slope and transit time.In particular, assuming the secondorder polynomial baseline yielded a larger anticorrelation between the TTVs and local slopes, of −79 ± 14 min day −1 .We therefore performed an injection-recovery procedure in which we injected transits at different phases in the Kepler 1627 light curve and repeated the TTV analysis.This was done at ≈50 phases, each separated by 0.02 P orb while omitting the phases in transit.For the second-order polynomial baseline, this procedure yielded a similar anti-correlation in the injected transits as that present in the real transit; assuming this baseline would therefore bias the result.However for the fourth-order baseline, the correlation present in the data was stronger than in all but one of the injected transits.Possible interpretations are discussed below.Given the lack of statistical significance, this analysis should be interpreted as suggestive at best.

Planet Confirmation
If the Kepler 1627Ab transit signal is created by a genuine planet, then to our knowledge it would be the youngest planet yet found by the prime Kepler mission. 5Could the transit be produced by anything other than a planet orbiting this nearsolar analog?Morton et al. (2016) validated the planet based on the transit shape, arguing that the most probable false positive scenario was that of a background eclipsing binary, which had a modeldependent probability of ≈10 −5 .However, this calculation was performed without knowledge of the low-mass stellar companion (M B ≈ 0.30 M ).Validated planets have also previously been refuted (e.g., Shporer et al. 2017).We therefore reassessed false positive scenarios in some detail.
As an initial plausibility check, Kepler 1627B contributes 1% to 2% of the total flux observed in the Kepler aperture.For the sake of argument, assume the former value.The observed transit has a depth of ≈0.18%.An ≈18% deep eclipse of Kepler 1627B would therefore be needed to produce a signal with the appropriate depth.The shape of the transit signal however requires the impact parameter to be below 0.77 (2-σ); the tertiary transiting the secondary would therefore need to be non-grazing with R 3 /R 2 ≈ 0.4.This yields a contradiction: this scenario requires an ingress and egress phase that each span ≈40% of the transit duration (≈ 68 min).The actual measured ingress and egress duration is ≈17 min, a factor of four times too short.The combination of Kepler 1627B's brightness, the transit depth, and the ingress duration therefore disfavor the scenario that Kepler 1627B might host the transit signal.
Beyond this simple test, a line of evidence that effectively confirms the planetary interpretation is that the stellar density implied by the transit duration and orbital period is inconsistent with an eclipsing body around the M-dwarf companion.We find ρ = 2.00 ± 0.24 g cm −3 , while the theoretically expected density for Kepler 1627B given its nominal age of 38 Myr and mass of 0.30 M is ≈4.6 g cm −3 (Baraffe et al. 2015).The transit duration is therefore too long to be explained by a star eclipsing the M dwarf secondary at 10σ.Adversarially assuming a younger age (32 Myr) and a larger companion mass (0.35 M ), the companion's density could be as low as ≈ 3.5 g cm −3 .This is still in tension with the fitted stellar density.While the planet might hypothetically still orbit a hidden close and bright companion, this possibility is implausible given i) the lack of secondary lines in the HIRES spectra, ii) the lack of secondary rotation signals in the Kepler photometry, and iii) the proximity of Kepler 1627 to the δ Lyr cluster locus on the Gaia CAMD (Figure 3).
The correlation noted in Section 4.1.3between the TTVs and the local light curve slope might be an additional line of evidence in support of the planetary interpretation.Unless it is a statistical fluke (a ≈5% possibility), then the most likely cause of the correlation is unresolved starspot crossings (Mazeh et al. 2015).These would only be possible if the planet transits the primary star, which excludes a background eclipsing binary scenario.The correlation would also suggest that the planet's orbit is prograde.The latter point assumes that the dominant photometric variability is induced by dark spots, and not bright faculae.Given the observed transition of Sun-like stellar variability from spot to faculae-dominated regimes between young and old ages, we expect this latter assumption to be reasonably secure (Shapiro et al. 2016;Montet et al. 2017;Reinhold & Hekker 2020).A third supporting line of evidence for the planetary interpretation also exists.We observed a transit of Kepler 1627Ab on the night of 2021 Aug 7 spectroscopically with HIRES at the Keck-I telescope and photometrically in griz bands with MuS-CAT3 at Haleakalā Observatory.Details of the observation sequence are discussed in Appendix C; Figure 7 shows the results.Although we did not detect the Rossiter-McLaughlin (RM) anomaly, the multi-band MuSCAT3 light curves show that the transit is achromatic.Quantitatively, when we fitted the MuSCAT3 photometry with a model that lets the transit depths vary across each bandpass, we found griz depths consistent with the Kepler depth at 0.6, 0.3, 0.3, and 1.1-σ respectively.The achromatic transits strongly favor Kepler 1627A as the transit host, since Kepler 1627B is a much redder star.Conditioned on the ephemeris and transit depth from the Kepler data, the MuSCAT3 observations also suggested a transit duration 17.3 ± 4.3 min shorter than the Kepler transits.However, given both the lack of TDVs in the Kepler data and the relatively low signal-to-noise of the MuS-CAT3 transit, further photometric follow-up would be necessary to determine whether the transit duration is actually changing.
For our RM analysis, the details are discussed in Appendix C. While the velocities are marginally more consistent with a prograde or polar orbit than a retrograde orbit, the spot-corrected exposure-toexposure scatter (σ RV ≈ 30 m s −1 ) is comparable to the expected RM anomaly assuming an aligned orbit (∆v RM ≈ 20 m s −1 ).We are therefore not in a position to claim a spectroscopic detection of the RM effect, nor to quantify the stellar obliquity.lation.Separately, transit spectroscopy aimed at detecting atmospheric outflows could yield insight into the evolutionary state of the atmosphere (e.g., Ehrenreich et al. 2015;Spake et al. 2018;Vissapragada et al. 2020).Observations that quantify the amount of high-energy irradiation incident on the planet would complement these efforts, by helping to clarify the expected outflow rate (e.g., Poppenhaeger et al. 2021).Finally, a challenging but informative quantity to measure would be the planet's mass.Measured at sufficient precision, for instance through a multi-wavelength radial velocity campaign, the combination of the size, mass, and age would yield constraints on both the planet's composition and its initial entropy (Owen 2020).
More immediately, the Kepler data may yet contain additional information.For instance, one possible explanation for the transit asymmetry shown in Figure 5 is that of a dusty asymmetric outflow.Dusty outflows are theoretically expected for young mini-Neptunes, and the amplitude of the observed asymmetry is consistent with predictions (Wang & Dai 2019).A second possibility is that the planetary orbit is slightly misaligned from the stellar spin axis, and tends to transit starspot groups at favored stellar latitudes.This geometry would be necessary in order to explain how the starspot crossings could add up coherently, given that the planetary orbital period (7.203 days) and the stellar rotation period (2.642 days) are not a rational combination.Other possibilities including gravity darkening or TTVs causing the asymmetry are disfavored (see Appendix B).
Beyond the asymmetric transits, Appendix D highlights an additional abnormality in the shortcadence Kepler data, in the arrival time distribution of stellar flares.We encourage its exploration by investigators more versed in the topic than ourselves.
It seems though that smaller planets could have been detected in the Kepler 1627 system: based on the per-target detection contours, the Kepler pipeline's median completeness extended to 1.6 R ⊕ at 10 day orbital periods, and 3.3 R ⊕ at 100 days (Burke & Catanzarite 2021).These limits account for the spot-induced variability in the system through a correction based on the Combined Differential Photometric Precision in the light curve over the relevant transit timescales (Burke & Catanzarite 2017).The large size of Kepler 1627Ab relative to most Kepler mini-Neptunes might therefore support a picture in which the typical 5 M ⊕ mini-Neptune (Wu 2019) loses a significant fraction of its primordial atmosphere over its first gigayear (Owen & Wu 2013;Ginzburg et al. 2018).It could also be consistent with a scenario in which an earlier "boil-off" of the planet's atmosphere during disk dispersal decreases the entropy of the planetary interior, leading to a rather long ∼10 8 year Kelvin-Helmholtz contraction timescale (Owen 2020).Confirming either of these scenarios would require a measurement of the planetary mass; otherwise, alternative explanations for its large size could include that it is just abnormally massive, or that it has an abnormally large envelope to core mass ratio.
Ultimately, the main advance of this work is a precise measurement of the age of Kepler 1627Ab.This measurement was enabled by identifying the connection of the star to the δ Lyr cluster using Gaia kinematics, and by then using the Gaia colorabsolute magnitude diagram and TESS stellar rotation periods to verify the cluster's existence.Table 3 enables similar cross-matches for both known and forthcoming exoplanet systems (e.g., Guerrero et al. 2021).Confirming these candidate associations using independent age indicators is essential because their false positive rates are not known.A related path is to identify new kinematic associations around known exoplanet host stars using positions and tangential velocities from Gaia, and to verify these associations with stellar rotation periods and spectroscopy (e.g., Tofflemire et al. 2021).Each approach seems likely to expand the census of planets with precisely measured ages over the coming years, which will help in deciphering the early stages of exoplanet evolution.also grateful to K. Collins for helping resolve the scheduling conflict that would have otherwise prevented the MuSCAT3 observations.L.G.B. acknowledges support from a Charlotte Elizabeth Procter Fellowship from Princeton University, as well as from the TESS GI Program (NASA grants 80NSSC19K0386 and 80NSSC19K1728) and the Heising-Simons Foundation (51 Pegasi b Fellowship).Keck/NIRC2 imaging was acquired by program 2015A/N301N2L (PI: A. Kraus).In addition, this paper is based in part on observations made with the MuSCAT3 instrument, developed by the Astrobiology Center and under financial support by JSPS KAKENHI (JP18H05439) and JST PRESTO (JPMJPR1775), at Faulkes Telescope North on Maui, HI, operated by the Las Cumbres Observatory.This work is partly supported by JSPS KAK-ENHI Grant Numbers 22000005, JP15H02063, JP17H04574, JP18H05439, JP18H05442, JST PRESTO Grant Number JPMJPR1775, the Astrobiology Center of National Institutes of Natural Sciences (NINS) (Grant Number AB031010).This paper also includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST).Funding for the TESS mission is provided by NASA's Science Mission directorate.We thank the TESS Architects (G.Ricker, R. Vanderspek, D. Latham, S. Seager, J. Jenkins) and the many TESS team members for their efforts to make the mission a continued success.Finally, we also thank the Keck Observatory staff for their support of HIRES and remote observing.We recognize the importance that the summit of Maunakea has within the indigenous Hawaiian community, and are deeply grateful to have the opportunity to conduct observations from this mountain.
(3) The eccentricity vectors are sampled in the (e cos ω, e sin ω) plane.(4) The true planet size is a factor of ((F1 + F2)/F1) 1/2 larger than that from the fit because of dilution from Kepler 1627B, where F1 is the flux from the primary, and F2 is that from the secondary; the mean and standard deviation of Rp = 3.817 ± 0.158 R⊕ quoted in the text includes this correction, assuming (F1 + F2)/F1 ≈ 1.015.3 is published in its entirety in a machine-readable format.This table is a concatenation of the studies listed in Table 4.One entry is shown for guidance regarding form and content.In this particular example, the star has a cold Jupiter on a 16 year orbit, HD 150706b (Boisse et al. 2012).An infrared excess has been reported (Cotten & Song 2016), and the star was identified by Ujjwal et al. (2020) as a candidate UMa moving group member (≈ 400 Myr; Mann et al. 2020).The star's RV activity and TESS rotation period corroborate its youth.4 describes the provenances for the young and age-dateable stars in Table 3. NGaia: number of Gaia stars we parsed from the literature source.NAge: number of stars in the literature source with ages reported.NG RP <16 : number of Gaia stars we parsed from the literature source with either GRP < 16, or a parallax S/N exceeding 5 and a distance closer than 100 pc.The latter criterion included a few hundred white dwarfs that would have otherwise been neglected.Some studies are listed multiple times if they contain multiple tables.Wenger et al. (2000) refers to the SIMBAD database.

APPENDIX
A. YOUNG, AGE-DATED, AND AGE-DATEABLE STAR COMPILATION The v0.6 CDIPS target catalog (Table 3) includes stars that are young, age-dated, and age-dateable.By "age-dateable", we mean that the stellar age should be measurable at greater precision than that of a typical FGK field star, through either isochronal, gyrochronal, or spectroscopic techniques.As in Bouma et al. (2019), we collected stars that met these criteria from across the literature.Table 4 gives a list of the studies included, and brief summary statistics.The age measurement methodologies adopted by each study differ: in many, spatial and kinematic clustering has been performed on the Gaia data, and ensemble isochrone fitting of the resulting clusters has been performed (typically focusing on the turn-off).In other studies however, the claim of youth is based on the location of a single star in the color-absolute magnitude diagram, or on spectroscopic information.
One major change in Table 3 relative to the earlier iteration from Bouma et al. ( 2019) is that the extent of Gaia-based analyses has now matured to the point that we can neglect pre-Gaia cluster memberships, except for a few cases with spectroscopically confirmed samples of age-dated stars.The membership lists for instance of Kharchenko et al. (2013) and Dias et al. (2014) (MWSC and DAML) are no longer required.This is helpful for various post-processing projects, since the field star contamination rates were typically much higher in these catalogs than in the newer Gaia-based catalogs.
The most crucial parameters of a given star for our purposes are the Gaia DR2 source_id, the cluster or group name (cluster), and the age.Given the hierarchical nature of many stellar associations, we do not attempt to resolve the cluster names to a single unique string.The Orion complex for instance can be divided into almost one hundred kinematic subgroups (Kounkel et al. 2018).Based on Figure 1, the δ Lyr cluster may also be part of a similar hierarchical association.Similar complexity applies to the problem of determining homogeneous ages, which we do not attempt to resolve.Instead, we simply merged the cluster names and ages reported by various authors into a comma-separated string.
This means that the age column can be null, for cases in which the original authors did not report an age, or for which a reference literature age was not readily available.Nonetheless, since we do prefer stars with known ages, we made a few additional efforts to populate this column.When available, the age provenance is from the original analysis of the cluster.In a few cases however we adopted other ages when stringbased cross-matching on the cluster name was straightforward.In particular, we used the ages determined by Cantat-Gaudin et al. (2020)  The catalogs we included for which ages were not immediately available were those of Cotten & Song (2016), Oh et al. (2017), Zari et al. (2018), Gagné et al. (2018a), Gagné et al. (2018b), Gagné & Faherty (2018), and Ujjwal et al. (2020).While in principle the moving group members discussed by Gagné et al. (2018a,b); Gagné & Faherty (2018) and Ujjwal et al. (2020) have easily associated ages, our SIMBAD crossmatch did not retain the moving group identifiers given by those studies, which should therefore be recovered using tools such as BANYAN Σ (Gagné et al. 2018b).We also included the SIMBAD object identifiers TT * , Y * O,Y * ?, TT?, and pMS * .Finally, we included every star in the NASA Exoplanet Archive planetary system (ps) table that had a Gaia identifier available (Akeson et al. 2013).If the age had finite uncertainties, we also included it, since stellar ages determined through the combination of isochrone-fitting and transit-derived stellar densities typically have higher precision than from isochrones alone.
For any of the catalogs for which Gaia DR2 identifiers were not available, we either followed the spatial (plus proper-motion) cross-matching procedures described in Bouma et al. (2019), or else we pulled the Gaia DR2 source identifiers associated with the catalog from SIMBAD.We consequently opted to drop the ext_catalog_name and dist columns maintained in Bouma et al. (2019), as these were only populated for a small number of stars.The technical manipulations for the merging, cleaning, and joining were performed using pandas (McKinney 2010).The eventual cross-match (using the Gaia DR2 source_id) against the Gaia DR2 archive was performed asychronously on the Gaia archive website.7 shows an anomaly at roughly the same transit phase.Year 2 correspondingly shows the strongest anomaly out of any year in Figure 10; the asymmetry is visually apparent however in each of the four years.
To bin by quartiles in local slope, we used our measurements of the local linear slopes in each of the observed transit windows (144 transits total).Four outlier transits were removed, leaving 140 transits.These were then divided into quartiles, so that each panel shows 35 transits binned together.The exact light curve slope intervals are listed in the lower left panels of Figure 11.Binned by local slope quartiles (Figure 11), the asymmetry is visually present in three of the four quartiles: the only bin in which it does not appear is d f /dt ∈ [3.0, 32.5] ppt day −1 .
Within the theory presented by Mazeh et al. (2015), unresolved starspot crossings cause the weak correlation between TTVs and the local light curve slope (Figure 6).In this model, we would expect the light curves with the most negative local slopes to have the largest positive TTVs, due to spot crossing events during the latter half of transit.The upper-left panel of Figure 11 agrees with this expectation.However, we would also expect the sign of the effect to reverse when considering the most positive local slopes (most negative TTVs).The lower-right panel of Figure 11 contradicts this expectation: the residual in both cases maintains the same parity!On the one hand, this shows that the residual is not dependent on the local light curve slope, which lowers the likelihood that it might be an artifact of our detrending methods.On the other, it raises the question of whether unresolved starspot crossings are indeed the root cause of the correlation shown in Figure 6.While we do not have a solution to this contradiction, the injection-recovery tests discussed in Section 4.1.3provide some assurance that the TTV-slope correlation is not simply a systematic artifact.

B.2. Interpretation
The transit asymmetry seems robust against most methods of binning the data, though with some caveats (e.g., the "middle quartile" in local flux, d f /dt ∈ [3.0, 32.5] ppt day −1 , where the asymmetry does not appear).Nonetheless, if the asymmetric were systematic we might expect its parity to reverse as a function of the sign of the local slope, and it does not.We therefore entertained four possible astrophysical explanations: gravity darkening, transit timing variations, spot-crossing events, and a persistent asymmetric dusty outflow.
Gravity darkening is based on the premise that the rapidly rotating star is oblate, and brighter near the poles than the equator (e.g., Masuda 2015).The fractional transit shape change due to gravity darkening is on the order of (P break /P rot ) 2 , for P break the break-up rotation period, and P rot the rotation period.Using the parameters from Table 2, this yields an expected 0.14% distortion of the ≈1.8 ppt transit depth: i.e., an absolute deviation of ≈2.5 ppm.The observed residual has a semi-amplitude of ≈ 50 ppm.Since the expected signal is smaller than the observed anomaly by over an order of magnitude, gravity darkening seems to be an unlikely explanation.
The scenario of transit timing variations (TTVs) producing the asymmetry seems unlikely because the transit timing variations we do observe are correlated with the local light curve slope, which increases roughly as much as it decreases.From our analysis, the mean TTV and its standard deviation are 0.66 ± 5.53 min; similarly the mean local slope and its standard deviation are 0.59 ± 45.50 ppt day −1 .There is therefore little expectation for TTVs to produce the asymmetry.A separate line of argument comes from Figure 11.If the local slope were essential to producing the transit asymmetry, we would expect that in the largest d f /dt bin, d f /dt ∈ [3.0, 32.5] ppt day −1 , the sign of the asymmetry would reverse.We do not see evidence for this being the case.
The third and related possibility is that of starspot crossings.Young stars have higher spot-covering fractions than old stars (e.g., Morris 2020).Young solar-type stars may also host dark starspots at high stellar latitudes (e.g., EK Dra; Strassmeier 2009), though interferometric imaging of spotted giant stars has shown different starspot latitude distributions than those inferred from Doppler imaging (Roettenbacher et al. 2017).Regardless, for any spot-crossing anomalies to add coherently over the 140 Kepler transits, it seems likely that we would need either for spots to be persistent at a particular latitude (and for the planetary orbit to be somewhat misaligned), or for a "stroboscopic" longitudinal phasing (e.g., Dai et al. 2018).For our system, P orb /P rot ≈ 2.76, which means that every 4 transits and 11 stellar rotations, the planet crosses over roughly the same stellar longitude, which might enable the necessary phasing if the spot-groups are large and long-lived.Unfortunately, the S/N per Kepler transit is ≈ 8, which renders individual spot-crossing events unresolved.This explanation seems marginally plausible, mainly because the expected spot-crossing anomaly amplitudes (≈ 100 ppm) resemble the observed amplitude of the asymmetry (≈ 50 ppm).One issue with this explanation however is that there is no reason to expect starspot crossing events to last exactly half the transit duration.A persistent feature of the planet itself might therefore be needed to explain the transit asymmetry.An asymmetric outflow from the planet's atmosphere could at least geometrically meet the requirements (e.g., McCann et al. 2019).To explain the asymmetric transit, a small, dense component would lead the planet, and a long, more rarefied (and variable) component would trail it.This might also explain the slight flux decrement visible for ∼1 hour pre-ingress (Figure 5).The amplitude of the asymmetry is roughly in line with theoretical expectations for dusty outflows (Wang & Dai 2019), and based on the planet's size, its mass is likely in a regime where such outflows are possible.Out of the four explanations discussed, this one at least theoretically seems the most plausible.By composition, the expectation would be that the envelope is mostly hydrogen and helium gas, with a dust or haze component providing the broadband opacity in the Kepler bandpass.A natural path for testing this idea would be to observe additional transits of the planet in hydrogen absorption, metastable helium absorption, or across a broad wavelength range in the near-infrared.
C. SPECTROSCOPY AND PHOTOMETRY DURING THE TRANSIT OF 2021 AUG 7 We used the ephemeris of Holczer et al. (2016) to observe a transit of Kepler 1627Ab on the night of 2021 Aug 7 both spectroscopically and photometrically.We used the HIRES echelle spectrograph at the Keck-I telescope and the MuSCAT3 photometer at Haleakalā Observatory on Maui, HI (Narita et al. 2020).For the HIRES wavelength calibration, we used the iodine cell, and extracted the 1-D spectra using the standard California Planet Survey pipeline (Howard et al. 2010).Given the faintness of the target (V = 13.1),we observed using the C2 decker, which yielded an instrument resolution of ≈50,000.The airmass ranged between 1.1 and 2.2 from the start through the end of observations; the seeing ranged from 1. 1 at the beginning to 1. 5 at the end.The HIRES exposure time was set at ≈15 minutes in order to resolve the 2.8 hour transit event, which yielded a S/N per resolution element of ≈ 75 (low for precision radial velocity standards).For the MuSCAT3 observations, we observed simultaneously across the griz bands.The exposure times in each bandpass ranged between 23 and 46 seconds, and were chosen in order to yield a S/N in the peak pixel that exceeded 130 while also preventing saturation.Performing aperture photometry on the latter image stack yielded the data given in Table 5.We considered two approaches to measuring the velocities: in the first, hereafter "Method 1", we crosscorrelated against a template found via spectral classification with SpecMatch-Emp (Yee et al. 2017).In "Method 2", we used a high S/N template of V1298 Tau.Although V1298 Tau is cooler than Kepler 1627A by ≈500 K, it has a comparable amount of line-broadening (v sin i = 23 km s −1 ), and a comparable level of stellar activity.The mean and standard deviation of the internal RV uncertainties averaged over all epochs were 16.2 ± 1.1 m s −1 from Method 1, and 12.6 ± 0.6 m s −1 from Method 2. The corresponding time-averaged reduced χ 2 from the template match was 1.57 ± 0.04 (Method 1) and 1.30 ± 0.02 (Method 2).Given these diagnostics, we adopted the velocities from the second approach, which are reported in Table 6.
Figure 7 shows the results.The MuSCAT3 photometry shows the expected starspot trend, along with the transit and what is likely a chromatic starspot crossing event at JD − 2459433 = 0.955.The radial velocities decrease by ≈250 m s −1 over the six hour window.This decrease in RV is correlated with a decrease in the S-indices derived from the Ca HK lines.One outlying RV point is apparent shortly before egress; it is temporally coincident with an outlying value in the S-index time series.
Overall, we expect the dominant trends in both the photometry and radial velocities to be caused by starspots on the stellar photosphere rotating into and out of view.The plasma in the leading and receding limbs of the stellar disk has an apparent line-of-sight velocity of ±20 km s −1 .Over 10% of a rotation cycle (P rot = 2.6 days), spots near these limbs come into and out of view, modulate the stellar velocity profile, and can thereby produce the overall ≈250 m s −1 trend.The Ca HK and Hα emission profiles support this interpretation; Figure 12 shows that each line gradually decreases in intensity over the course of the six hour sequence.
The expectation however is for the starspot-induced signals to be smooth, at worst with contributions at 0.5 P rot or 0.25 P rot (Klein & Donati 2020).We therefore fitted the RVs using the Hirano et al. (2010Hirano et al. ( , 2011) ) models for the Rossiter-McLaughlin (RM) effect, and allowed for an optional linear and quadratic trend in time to fit the ≈250 m s −1 spot-induced trend.We followed the methodology developed by Stefansson et al. (2020).We allowed the sky-projected obliquity, the projected stellar equatorial velocity, and the Gaussian dispersion of the spectral lines to vary, and fixed the limb-darkening using the V -band tabulation from Claret & Bloemen (2011).We assumed a Gaussian prior on v sin i and a/R from Table 1, and also allowed for a white-noise jitter term to be added in quadrature to the measurement uncertainties.We used a 15 minute exposure time to numerically evaluate the model.
The quadratic model with the RM effect is shown in Figure 7; the jitter term is incorporated in the model uncertainties, but not the plotted measurement uncertainties.The plotted measurement uncertainties are the internal uncertainties on the RVs (≈ 13 m s −1 ), and are dominated by the v sin i broadening.However, between exposures, the RVs show significant additional scatter that is not captured by the slow quadratic trend.The white-noise jitter for this particular model is σ RV = 27 +6 −5 m s −1 , which is comparable to the expected RM anomaly of ∆v , assuming a perfectly aligned orbit.The presence of this additional scatter prevents a convincing detection of the RM effect.The reason can be understood via model comparison.If we compare the model with a quadratic trend and the RM effect against a model with a linear trend and the RM effect, or even a model with no RM effect at all, then the respective Bayesian Information Criterion (BIC) values are as follows.(C1) There is therefore no evidence to prefer the model with the RM effect against a model that only accounts for the stellar variability.The "only quadratic" model does particularly well because it can inflate the jitter term to account for scatter during the transit (even if the scatter contains astrophysics!), and it has fewer free parameters.However, we cannot justify a physical prior on the jitter term, because we do not understand the origin of the exposure-to-exposure scatter.As noted above, the velocity deviations from starspots are expected to have contributions at the stellar rotation frequency, or harmonics thereof.This jitter is present on the exposure timescale (15 minutes), which is only 0.4% of the stellar rotation period; it is not obvious that starspots would be the culprit.
The amplitude of both the spot-induced trend and the jitter are somewhat larger than recent comparable measurements in systems such as AU Mic (Palle et al. 2020), DS Tuc (Montet et al. 2020;Zhou et al. 2020) and TOI 942 (Wirth et al. 2021).One possible explanation for the jitter is that it is astrophysical in origin, and that it is caused by some novel process operating on the surface of Kepler 1627A.Another possibility is that our RV analysis underestimates our measurement uncertainties; in order to achieve the requisite timesampling the S/N per resolution element in our spectra was 70 to 80, which is lower than desired for deriving high-precision velocities.In addition, the rapid rotation of the star could affect accuracy of the uncertainties from the velocity extraction.Observations at higher S/N are necessary to distinguish these two possibilities, and remain worthwhile in order to clarify the orbital geometry of Kepler 1627Ab.Useful next steps would include transit observations with a stabilized spectrograph in the optical (e.g., Gibson et al. 2016;Seifahrt et al. 2018), or in the near-infrared (e.g., Feinstein et al. 2021).D. FLARE ANALYSIS In addition to the 3.9 years of long cadence data, short cadence (1-minute) Kepler observations were acquired over 97.7 days during Quarter 15.The short cadence light curve shows a higher rate of flaring than visible in the long cadence data (Figure 13).We analyzed the short cadence light curve and its flares according to the following procedure.
1. Fit the starspot-induced variability using a Gaussian Process with a SHOTerm kernel, a white-noise jitter term, and the mean flux.
2. Select points more than twice the median absolute deviation from the residual, and exclude them from the light curve (these points include the flares).Repeat Step 1.

Using the residual from
Step 2, identify all flares, requiring them to be at least 20 cadences apart, at least 7 median absolute deviations above the median baseline, and lasting at least 2 cadences in duration.Build the mask spanning these times, from 5 minutes before each flare begins to 2.5 minutes after the final flare cadence.Repeat Step 1 a final time.
The final step of flare identification and fitting was performed using altaipony (Davenport 2016;Ilin et al. 2021).The analytic flare model is from Davenport et al. (2014) and it parametrizes the flare with a start time, an exponential lag time, and an amplitude.There were N f = 24 flares that exceeded 0.5% in relative flux during the short cadence observations.These 24 flares spanned a total of 6.5 hours (∼15 minutes per flare).Inspecting the data, we noticed a coincidence in the flare arrival times.The coincidence is that despite the low flare duty cycle, one orbital period after the brightest flare, a second flare followed.This and a similar event are shown in Figure 13.The timing error is good to a ≈ 0.2% difference from the orbital period, which given the duty cycle seems a priori unlikely.If we consider flares falling within 2% of the planet's orbital period after a previous flare, then 4 of the 24 flare events have candidate "successors".
As with any coincidence, if one does not have a firm prediction, it is difficult to assess whether a surprise is statistically significant.Since our surprise was specifically at the inter-arrival time of certain flares coinciding with special time intervals, we performed the following analysis.First, we considered all unordered pairs of flares.For N flares there are n 2 such pairs (for our case, 276 pairs).We then compared the distribution of the pair separations against that of a Poisson distribution.Specifically, we drew N f = 24 samples from a Poisson distribution with λ = ∆t/N f , for ∆t = 97.7 days the total duration of the observations, and repeated the draw 10 3 times with unique random seeds.
Figure 14 shows the results.The vertical lines in the figure show the planetary orbital period, the synodic period P syn = (P −1 rot − P −1 orb ) −1 , and linear combinations thereof.The tidal period (half the synodic period) is not shown.The bins are logarithmically spaced to give 100 bins between the minimum and maximum ordinate values.The gray bands express the range of values observed from the Poissonian draws.While it does seem like an odd coincidence for peaks in the observed flare arrival time distribution to coincide with the locations of these "special intervals", the statistical evidence for a non-Poissonian process driving the flares does not seem especially overwhelming.More quantitatively, the peaks observed at the orbital and synodic periods are within the ±2-σ range of a Poissonian process, and those at P orb + P syn and P orb + 2P syn are only slightly above this range.With that said, future analyses of these data by investigators with more knowledge of this topic could very well yield more quantitative insights.Such analyses should keep in mind an important caveat: the amplitude distribution of M-dwarf flares extends up to many times the quiescent flux (see Figure 7 of Günther et al. 2020).A flare on Kepler 1627B producing double its quiescent white-light flux would yield a ≈1% apparent amplitude.Such flares could represent a significant fraction of those in the Kepler observations.

Figure 1 .Figure 2 .
Figure 1.Galactic positions and tangential velocities of stars in the δ Lyr cluster.Points are reported cluster members fromKounkel & Covey (2019).The tangential velocities relative to Kepler 1627 (bottom right) are computed assuming that every star has the same three-dimensional spatial velocity as Kepler 1627.Our analysis considers stars (black points) in the spatial and kinematic vicinity of Kepler 1627 (yellow star).The question of whether the other candidate cluster members (gray points) are part of the cluster is outside our scope.The location of the Sun is ( ) is shown to clarify the direction along which parallax uncertainties are expected to produce erroneous clusters.
Figure 3 shows the color-absolute magnitude diagram (CAMD) of candidate δ Lyr cluster members, IC 2602 (≈ 38 Myr), the Pleiades (≈ 115 Myr), and the field.The stars from the Pleiades and IC 2602 were adopted from Cantat-Gaudin et al. (2018), and the field stars are from the Gaia EDR3 Catalog of Nearby Stars (Gaia Collaboration et al. 2021b).

Figure 3 .
Figure 3.The δ Lyr cluster is 38 +6 −5 Myr old.Top: Color-absolute magnitude diagram of candidate δ Lyr cluster members, in addition to stars in IC 2602 (≈ 38 Myr), the Pleiades (≈ 115 Myr), and the Gaia EDR3 Catalog of Nearby Stars (gray background).The zoomed right panel highlights the pre-main-sequence.The δ Lyr cluster and IC 2602 are approximately the same isochronal age.Bottom: TESS and Kepler stellar rotation period versus dereddened Gaia color, with the Pleiades and Praesepe (650 Myr) shown for reference(Rebull et al. 2016;Douglas et al. 2017).Most candidate δ Lyr cluster members are gyrochronally younger than the Pleiades; outliers are probably field interlopers.

Figure 4 .
Figure 4. Kepler 1627 is a binary.Left: High-resolution imaging from Gemini-North/'Alopeke and Keck/NIRC2shows an ≈M3V companion at ρ ≈ 0. 16, which corresponds to a projected separation of 53 ± 4 AU.The inset shows a cutout of the stacked NIRC2 image (North is up, East is left, scale is set by the separation of the binary).The lines show 5-σ contrast limits for the 'Alopeke filters, and 6-σ contrast limits for NIRC2 outside of 0. 15.Right: Gaia EDR3 renormalized unit weight error (RUWE) point estimates for candidate δ Lyr cluster members.Since other members of the cluster with similar brightnesses have comparable degrees of photometric variability, the high RUWE independently suggests that Kepler 1627 is a binary.

Figure 5 .
Figure 5.The light curve of Kepler 1627.Top: The Kepler data span 1,437 days (3.9 years), sampled at 30 minute cadence; a 100 day segment is shown.The top panel shows the PDCSAP median-subtracted flux in units of parts-perthousand (×10 −3 ).The dominant signal is induced by starspots.The stellar variability model (orange line) is subtracted below, revealing the transits of Kepler 1627Ab.The online Figure Set spans the entire 3.9 years of observations.Bottom: Phase-folded transit of Kepler 1627Ab with stellar variability removed.Windows over 20 hours (left) and the entire orbit (right) are shown, and the residual after subtracting the transit is in the bottom-most row.The 2-σ model uncertainties and the best-fit model are the light purple band and the dark purple line.Gray points are individual flux measurements; black points bin these to 15 minute intervals, and have a representative 1-σ error bar in the lower right of each panel.The asymmetric residual during transit is larger than the out-of-transit scatter.
show an asymmetric residual in the data relative to the model: the measured flux is high during the first half of transit, and low in

Figure 7 .
Figure 7. Keck/HIRES radial velocities and MuS-CAT3 photometry from the transit of 2021 Aug 7. Top: The radial velocity jitter across the 15 minute exposures (σRV ≈ 30 m s −1 ) prevented us from detecting the RM effect; a model including the RM anomaly and a quadratic trend in time to fit the spot-induced ≈ 250 m s −1 trend is shown (see Appendix C).Shaded bands show 2-σ model uncertainties.Middle: The RV variations are strongly correlated with varying emission in the Ca H and K lines.Bottom: The photometric transit depths are consistent across the griz bandpasses.The photometry is binned at 10 minute intervals.

Figure 8 .
Figure 8. Radii, orbital periods, and ages of transiting exoplanets.Planets younger than a gigayear with τ /στ > 3 are emphasized, where τ is the age and στ is its uncertainty.Kepler 1627Ab is shown with a star.The large sizes of the youngest transiting planets could be explained by their primordial atmospheres not yet having evaporated; direct measurements of the atmospheric outflows or planetary masses would help to confirm this expectation.Selection effects may also be important.Parameters are from the NASA Exoplanet Archive (2021 Sept 15).
B. THE TRANSIT ASYMMETRY B.1.How Robust is the Asymmetric Transit?As a means of exploring the robustness of the transit asymmetry, Figures 9, 10, and 11 show the Kepler data binned in three ways: over Kepler quarters, Julian years, and quartiles of local slope.Over Kepler quarters (Figure 9), Quarter 6 shows the strongest asymmetry: a deviation of about 3 ppt from expectation.Quarter2

Figure 9 .
Figure 9. Transit model residuals through time (binned by Kepler quarter).Left: Phase-folded transit of Kepler 1627b, with stellar variability removed.Black points are binned to 20 minute intervals.The 2-σ model uncertainties and the maximum a posteriori model are shown as the faint purple band, and the dark purple line.Right: As on the left, with the transit removed.
Figure 10.Transit model residuals through time (binned by year of observation).Left: Phase-folded transit of Kepler 1627b, with stellar variability removed.Points and models are as in Figure 9. Right: As on the left, with the transit removed.

Figure 11 .
Figure 11.Transit models and residuals, binned by quartiles in the local slope of the light curve.Representative uncertainties for the black points (binned at 20 minute intervals) are shown in the lower right of each panel.A similar transit asymmetry to that shown in Figure 5 seems to be present in three of the four bins.

Figure 12 .
Figure 12.Spectroscopic activity indicators during the transit of 2021 Aug 7. The top panels show the median line profiles Ca K, Ca H, and Hα line profiles from the HIRES spectra.The lower panels show the differences of each individual spectrum relative to the median spectrum.The bump in the red wing of Ca H is H .The spectra in the lower panels are smoothed for visualization purposes.
extreme of 46 Myr, the δ Lyr cluster would be 44 +8 −7 Myr old.We adopt an intermediate 38 Myr age for IC 2602, which yields an age for the δ Lyr cluster of 38 +6 −5 Myr. 3 Follow-up studies of the LDB or main-sequence turn-off in the δ Lyr cluster could help determine a more precise and accurate age for the cluster, and are left for future work.

Table 2 .
Priors and Posteriors for Model Fitted to the Long Cadence Kepler 1627Ab Light Curve.
NOTE-ESS refers to the number of effective samples.R is the Gelman-Rubin convergence diagnostic.Logarithms in this table are base-e.U denotes a uniform distribution, and N a normal distribution.(1) The ephemeris is in units of BJDTDB -2454833.

Table 3 .
Young, Age-dated, and Age-dateable Stars Within the Nearest Few Kiloparsecs.

Table 4 .
Provenances of Young and Age-dateable Stars.

Table 5 .
MuSCAT3 photometry of Kepler 1627.Table 5 is published in its entirety in a machine-readable format.Example entries are shown for guidance regarding form and content.

Table 6 .
Kepler 1627 radial velocities.Table 6 is published in its entirety in a machinereadable format.Example entries are shown for guidance regarding form and content.