ASTEROSEISMIC DETERMINATION OF OBLIQUITIES OF THE EXOPLANET SYSTEMS KEPLER-50 AND KEPLER-65

The obliquities of the host stars in exoplanetary systems display a surprising diversity, including low obliquities reminiscent of the solar system, strongly tilted stars, and retrograde systems in which the directions of stellar rotation and planetary orbital revolution are opposite. Most of these results have been obtained by detecting the Rossiter–McLaughlin (RM) effect, a spectroscopic anomaly that is observed during a planetary transit (Queloz et al. 2000; Winn et al. 2005). In addition, for some systems the obliquity has been determined through the detection and interpretation of transits of a planet over starspots (e.g., Deming et al. 2011; Désert et al. 2011; Nutzman et al. 2011; Sanchis-Ojeda et al. 2011), and for one system (Barnes et al. 2011; Szabó et al. 2011) it has been estimated using the signatures of gravity darkening from rapid stellar rotation (Barnes 2009).

Almost all of the previous results pertain to host stars with hot Jupiters. The diversity of obliquities seen in those systems has been taken as evidence that the process of converting a "normal" Jupiter into a hot Jupiter can tilt the inclination of the planetary orbit (see, e.g., Winn et al. 2010b; Triaud et al. 2010; Albrecht et al. 2012). It would be interesting to extend these measurements to systems with longer-period planets, and multiple-planet systems, to test whether the high obliquities are indeed confined to the hot-Jupiter systems. Unfortunately, the long-period and multiple-planet systems tend to involve smaller planets and intrinsically fainter host stars (Latham et al. 2011; Steffen et al. 2012), making it difficult to apply the RM and starspot techniques. This is why only two such systems have been examined to date (Sanchis-Ojeda et al. 2012; Hirano et al. 2012a). It would be advantageous to develop a technique that does not depend so critically on the signal-to-noise ratio (S/N) of the transit data. One possibility is to use a combination of the measured rotation period (P_rot), the projected rotation rate (vsin i_s), and the stellar radius (R) to determine sin i_s, the sine of the angle i_s between the stellar rotation axis and the line of sight (see, e.g., Hirano et al. 2012b). However, this method is usually limited by the relatively poor accuracy of vsin i_s measurements for cool stars.

Asteroseismology provides another potentially powerful method. The detection and interpretation of the solar-like oscillations shown by solar-type stars is well known to provide accurate fundamental properties of host stars (e.g., Bazot et al. 2005; Bouchy et al. 2005; Vauclair et al. 2008; Soriano & Vauclair 2010; Christensen-Dalsgaard et al. 2010; Ballot et al. 2011b; Batalha et al. 2011; Gilliland et al. 2011; Howell et al. 2012; Borucki et al. 2012; Carter et al. 2012; Escobar et al. 2012; Barclay et al. 2013). Less well known is that in some cases the rotationally induced splittings of oscillation modes can be used to determine i_s (Gizon & Solanki 2003). When the host star also has a transiting planet, the inclination i_p of the planetary orbit can be determined, and therefore the difference in inclination between the star and planetary orbit can be calculated.

In contrast to the RM and starspot techniques, the applicability of the asteroseismic method depends predominantly on the stellar parameters and hardly at all on the planetary parameters, giving the asteroseismic method a decisive advantage in measuring the stellar obliquities in systems with small planets or long-period planets. The asteroseismic analysis does, however, require bright targets and long-duration, high-cadence photometric time series to give the requisite S/N and frequency resolution for extracting clear signatures of rotation from the oscillation spectrum, and hence the stellar inclination angle.

Here we present the first application of asteroseismology to the problem of stellar obliquity determination for Sun-like exoplanet hosts with transiting planets. Both of the systems considered in this paper have solar-type stars hosting multiple, small (super-Earth sized) transiting planets. The identification of the two-planet Kepler-50 system (KOI-262, KIC 11807274), along with the validation of the transit signals as arising from planets, was previously reported by Steffen et al. (2013). Kepler-65 (KOI-85, KIC 5866724) is a three-planet system that is herein identified and validated for the first time. Both systems involve F-type stars at the brighter end of the Kepler target list, having apparent magnitudes of K_p = 10.42 and K_p = 11.02, respectively.

Previously, asteroseismic methods have been applied to host stars with single, non-transiting large planets discovered using the Doppler method—HD 52265, a solar-type host with asteroseismic data from CoRoT (Ballot et al. 2011b; and HR 8799, an A-type host showing γ Doradus pulsations in ground-based observations (Wright et al. 2011)—with only moderate constraints returned on the stellar inclinations.

The rest of the paper is organized as follows. We begin in Section 2 by estimating the fundamental stellar properties, using the solar-like oscillations detected in the Kepler light curves and complementary spectroscopic data. Section 3 presents the planet properties of both systems, including validation of the planets orbiting Kepler-65 and discussion of the mutual inclinations of the planetary orbits of both systems. The asteroseismic estimation of the stellar obliquities, which depends on extracting signatures of rotation from the oscillation spectra, is presented in Section 4. Section 5 compares the asteroseismic results on rotation with independent estimates of the surface rotation based on the quasi-periodic variations seen in the Kepler light curves, and measurements of the sky-projected surface rotational velocity based on spectroscopic line broadening. We finish in Section 6 with a discussion of the implications of our results for theories of planetary formation.

We determined the fundamental stellar properties of Kepler-50 and Kepler-65 by comparing a few key asteroseismic and spectroscopic observables to the outputs of stellar-evolutionary models.

The asteroseismic results are based on the Kepler short-cadence (SC) data (Gilliland et al. 2010), whose one-minute sampling is needed to detect the short-period oscillations observed in solar-type stars (see also Chaplin et al. 2011b). The light curve for Kepler-50 spans 18 months, from Kepler observing quarters 6 through 11 inclusive. The light curve for Kepler-65 spans 27 months, from quarters 3 through 11.

Before computing power spectra, the planetary transit signals were removed from the time series by applying a median high-pass filter of width appropriate for the transit durations (see, e.g., Christensen-Dalsgaard et al. 2010). The clear separation of the relevant timescales—i.e., periods of days associated with the transits versus periods of minutes associated with the dominant oscillations—means that this approach cleans the frequency-power spectrum in such a way as to allow the asteroseismic analysis to proceed unhindered. Figure 1 shows frequency-power spectra of the light curves of Kepler-50 (top panel) and Kepler-65 (bottom panel). The spectra were computed using a Lomb–Scargle periodogram (Scargle 1982), and calibrated to satisfy Parseval's theorem. Both stars present clear patterns of peaks due to solar-like oscillations, which are small-amplitude pulsations that are stochastically excited and intrinsically damped by the near-surface convection. Many acoustic (pressure, or p) modes of high radial order, n, are excited to observable amplitudes. Solar-type stars oscillate in both radial and non-radial modes. The modes may be decomposed onto spherical harmonic functions of degree l. Both stars show detectable overtones of modes with l ⩽ 2.

Zoom In Zoom Out Reset image size

2.1. Spectroscopic Data and Analysis

2.2. Asteroseismic Estimation of Stellar Properties

It is important to establish whether the transit-like photometric signals represent actual transits of a system of planets across the disk of the intended target star, as opposed to a "false positive" such as a system of eclipsing stars blended with the intended target star. For Kepler-50, transit timing variations (TTVs) have been observed for both of its planets and are anti-correlated, a clear sign that the planets are interacting with each other and hence orbit the same star (Steffen et al. 2013). Kepler-65 has not been confirmed in this manner; in the following section we validate the system by other means.

3.1. Validation of Kepler-65

To validate the Kepler-65 system, in this section we will demonstrate that: (1) background eclipsing binaries are unlikely to be responsible for any of the three candidate transit signals; (2) all three transiting objects are likely orbiting the same star, which must have a mean density very similar to that of the intended target star; and (3) planets c and d are near a 7:5 mean-motion commensurability, and the smaller planet in this pair (planet d) exhibits a significant TTV signal of the nature expected for such a configuration.

3.1.1. Excluding Background Binary Scenarios

Lissauer et al. (2012) considered the question of how many of Kepler's multiple-planet candidates actually represent true multiple-planet systems, as opposed to unresolved blends of systems each having only one eclipsing object. For example, a candidate two-planet system could actually be a single-planet system along nearly the same line of sight to a background eclipsing binary, or there could be two eclipsing binaries along the same line of sight whose eclipses are diluted to planet-like proportions by the constant light of a foreground star. Lissauer et al. (2012) recognized that false positives of this nature would be randomly distributed among the target stars, and that the number of multiple-planet candidates is much larger than would be expected if the candidates were assigned randomly to target stars. From this analysis they concluded that the vast majority of Kepler's multiple-planet candidates do not represent superpositions of singly eclipsing systems. For the population of three-transit candidates such as Kepler-65, Lissauer et al. (2012) estimated the chance that at least one of the candidates represents an unrelated eclipsing system is 0.07% (an expectation of 0.13 such false positives out of 178 candidates).

These general considerations show that Kepler-65 is very likely a true multiple-planet system, as opposed to unrelated singly eclipsing systems that are blended together in the Kepler photometric aperture. In the remainder of this section we examine the specific circumstances and follow-up observations of Kepler-65 that also support this conclusion.

The photometric aperture used for the star changes from quarter to quarter, but in all cases has a size of approximately 4 × 4 pixels. With a detector scale of 3.98 arcsec pixel⁻¹, stars within a radius of about 12 arcsec from Kepler-65 could contribute light to the aperture and could in principle be the source of some of the transiting signals. We checked for possible contaminating stars using two different data sets.

First, we consulted the catalog by Adams et al. (2012) of adaptive optics (AO) images of a large sample of Kepler Objects of Interest (KOIs). The range of star magnitudes that can be detected depends on the distance to the star, such that one loses the ability to detect faint stars very close to the main star. These images have a range of 6 arcsec, outside of which no information was provided¹⁷ Only one other star was detected on the AO image, at a separation of 2.9 arcsec from Kepler-65. Adams et al. (2012) estimated that the Kepler apparent magnitude of this star is K_p ≈ 21, i.e., about 10 mag fainter than Kepler-65.

Second, to seek companions outside the 6 arcsec radius, we consulted the Naval Observatory Merged Astrometric Dataset (Zacharias et al. 2004). In this catalog, 11 stars are detected within a box of 30 arcsec centered on the position of Kepler-65. Only three stars were found that could be candidates for a background blend; one at a separation of 7.2 arcsec with an R-band magnitude of 18.6 (compared to 10.5 for Kepler-65), a second star at a separation of 11.5 arcsec with an R-band magnitude of 14.1, and a third star at a separation of 11.8 arcsec with a B-band magnitude of 19.3 (11.6 for Kepler-65).

If one of these objects were a background binary star mimicking a transiting planet, then the spatial, first-moment centroid of the light gathered on the aperture would be displaced during eclipses by an amount approximately equal to the observed transit depth multiplied by the projected distance from the object to Kepler-65. Therefore, since the transit depth is known independently with high precision, an upper bound on the centroid displacement can be used to set a maximum distance at which a contaminating binary can be located (also known as the radius of confusion). This notion has been applied to detect background binaries among the KOIs (see, e.g., Batalha et al. 2010), to estimate false alarm probabilities (FAPs) for particular KOIs (Morton & Johnson 2011) and to validate individual candidates using the BLENDER technique (Torres et al. 2011).

In the case of Kepler-65, one pixel of the stellar image is saturated, and consequently the distribution of light does not follow the standard point-spread function. Rather than attempting to model the saturated point-spread function, we used the flux-weighted column and row centroids produced by the Kepler pipeline. With this method of computing centroids we are only sensitive to displacements larger than ≈1 pixel (≈4 arcsec), but this is sufficient for our purpose. For each transit observed at SC, we selected a window in time of width 4.8 hr centered on each transit. Transits that occurred within 6 hr of another were excluded. To eliminate the effects of outliers we omitted data points differing by more than 3σ from a median-smoothed version of the time series, where the smoothing was performed over 30 minute intervals. We found that the centroid motion was approximately a linear function of time, presumably because of the continuous pointing drift of the telescope. We corrected for this effect by fitting the out-of-transit portions of the data set with a linear function of time. All the centroid information for each candidate in a given quarter was then phase-folded using a linear transit ephemeris from the KOI input catalog (Batalha et al. 2013). The sequences of in-transit and out-of-transit centroids were approximated as Gaussian distributions, and the centroid displacement was computed as the difference between the means of the distributions, with an uncertainty based on standard error propagation.

The left-hand panel of Figure 2 shows the measured centroid displacements after dividing by the corresponding transit depths, so that implied physical distances are plotted. The signals cannot originate from a source outside the 6 arcsec radius (shown in blue) since no points lie outside that range. The right-hand panel of Figure 2 shows this constraint, in combination with constraints from other considerations, which together limit the radius of confusion to 0.7 arcsec. One constraint is the AO imaging described previously. Another is that binary stars just over 9 mag fainter than Kepler-65 cannot decrease the total amount of light by 100 ppm (the depth of the shallowest transit), even were they to have eclipses of 50% depth, i.e., as given by ΔK_p = −2.5log₁₀[100 × 10⁻⁶/0.5] ≃ 9.2. This limit is marked as the horizontal dotted line in right-hand panel of Figure 2. No stars brighter than that limit and within 0.7–6 arcsec from Kepler-65 were detected within the AO image. With such a small radius of confusion, the probability that any of the signals come from background binaries is small. We estimate this probability to be ⩽0.15% for each candidate of Kepler-65, based on the work by Morton & Johnson (2011), who calculated the local surface density of eclipsing binaries whose properties could mimic those of each Kepler candidate. Specifically we took their estimated false-positive probability of ⩽1%, which assumed a radius of confusion of 2 arcsec, and scaled it by (0.7/2.0)² (because we have demonstrated that the true radius of confusion for Kepler-65 is 0.7 arcsec).

Zoom In Zoom Out Reset image size

Finally there is the "multiplicity boost," as discussed at the beginning of this section. Since there are three-transit candidates for a single Kepler target, the false alarm probability of each individual transiting object is further reduced, according to the statistical argument of Lissauer et al. (2012). Here the boost factor is approximately 50, which would reduce the individual false-positive probabilities from ⩽0.15% to ⩽3 × 10⁻⁵.

3.1.2. Evidence That the Three Planets Orbit Kepler-65

We have demonstrated that it is unlikely that any of the transit candidates arises from a background eclipsing binary. We next ask whether the signals represent three planets all orbiting the intended target star Kepler-65, or whether any of them could actually be orbiting a companion star that is gravitationally bound to it. To address this question we searched for a pattern in the transit observables that would suggest that all the planets orbit the same star; namely, transit durations scaling as the cube root of the orbital period, which is a sign that they transit a star with similar density (see, e.g., Fabrycky et al. 2012; Lissauer et al. 2012). If some of the transit signals represented transits across a different star, then such a pattern would occur only by coincidence.

To measure transit durations, we constructed phase-folded transit light curves for each of the three candidates using the SC data and assuming a constant orbital period. Transits that occurred within 6 hr of another were excluded, removing all possible overlapping transits. The data were binned into 7.5 s intervals to increase the speed of subsequent computations. We used a standard description of the loss of light due to a transiting planet (Mandel & Agol 2002) to model the binned light curves simultaneously. We adopted a quadratic limb-darkening law, with the two coefficients left as free parameters (and shared by all three candidates). The free parameters describing each light curve were the squared planet-to-star radius ratio (R_p/R)², the impact parameter, b, and the stellar radius divided by the orbital distance, R/a. We then found the best-fitting model parameters that minimized the standard χ² function, with uncertainties on the measurements defined as the standard deviation of the points outside transit for each of the folded light curves. A Markov Chain Monte Carlo (MCMC) code was then used to explore the range of allowed parameters.

From the best-fitting model parameters we computed the transit durations, defined as the interval over which the center of the planet is projected in front of the stellar disk. This parameter is generally well constrained, and does not change much in the presence of a small TTV signal (whereas the ingress duration would experience larger fractional variations). The transit durations are plotted in Figure 3, as a function of orbital period. We compared these values with those expected for planets in circular orbits around a star with a mean density equal to 0.621 g cm⁻³, which is the mean density of Kepler-65 as estimated from the asteroseismic analysis (see Table 1 and Section 2). The measured durations agree well with a model in which all planets have the same orbital inclination, which is good evidence that the planets have nearly coplanar and circular orbits around a single star with a density similar to that of Kepler-65.

Zoom In Zoom Out Reset image size

Another way to perform this test is to use the transit observables to compute the implied mean density of the host star, and compare the result to the mean density obtained from asteroseismology. To this end we performed a second fit to the data in which the R/a value for Kepler-65c (the candidate with the highest S/N) was a free parameter, and the R/a values for the other two planets were fixed according to the assumption of circular orbits around the same star (i.e., by scaling according to orbital period and Kepler's third law). This effectively introduces a constraint that all three planetary signals agree on the mean stellar density (Seager & Mallén-Ornelas 2003). We found the photometrically derived mean density to be $0.57^{+0.06}_{-0.07} \,\rm g\,cm^{-3}$, in agreement with the asteroseismically derived mean density of 0.621  ±  0.011 g cm⁻³ (see Table 1 and Section 2). Therefore, the transit observables are consistent with a system of three planets on nearly circular orbits around a star with the same mean density as the star that is the source of the observed p-mode oscillations.

A devil's advocate would raise the possibility that this agreement is a coincidence, and that one or more of the planets actually orbit a secondary companion star. This seems unlikely indeed although we do not attempt to assign a quantitative false-positive probability to this scenario. To establish the probability of such a coincidence one would need to consider a realistic distribution of companions, along with their planets, transit probabilities (which may be correlated with the transit probabilities of the primary star), and transit durations. One would then need to exclude cases in which the companion would have been detectable in the optical spectrum, the spectrum of p-mode oscillations¹⁸ (i.e., by contributing signatures of its own oscillations), or through excessively diluted transit depths, and then compute the integrated probability of the allowed phase space. This is beyond the scope of this study.

3.1.3. Detection of a TTV Signal for Kepler-65d

The detection of TTVs have proven to be useful for validating planets as well as constraining the masses of the transiting planets. We performed a transit-timing analysis of Kepler-65 as follows. To measure individual transit times for each planet we employed a phase-folded light curve as a template function. Specifically we used a phase-folded light curve that was obtained by fitting all of the transits under the assumption of a circular orbit with a constant period. The template was then fitted to the data from each transit observed at SC, with three free parameters: the central time of the transit, the out-of-transit flux level, and a constant gradient in the out-of-transit flux level. An MCMC code was used to obtain the posterior distribution for the time of transit. Since no SC data were available in Q0, Q1, or Q2, for those quarters we used the transit times from the Kepler TTV catalog (J. F. Rowe et al. 2013, in preparation), measured as described by Ford et al. (2011).

The individual transit times were then fitted with a linear function of epoch, and this function was subtracted from the timing data to isolate any timing residuals. A visual inspection showed no obvious TTV signal in the residuals. We fitted these residuals with a sinusoidal model with three parameters: a TTV period, phase, and amplitude. To facilitate the exploration of the parameter space, we divided the range of periods into small intervals covering periods from 10 to 500 days, and for each period we optimized the other two parameters. For the three planets, the best-fitting sinusoids gave unacceptably high χ² values relative to the number of degrees of freedom (2200, with 356 points for planet b, 198 with 121 points for planet c, and 200 with 81 points for planet d). The uncertainties on the individual transit times are likely underestimated due to correlated noise in the photometric time series (due to some combination of stellar granulation, p-modes, and instrumental noise). We proceeded by enlarging the uncertainties by a scale factor (see Table 2) such that the minimum χ² was equal to the number of degrees of freedom.

Table 2. TTV Signals of Kepler-65

Planet	N	Error Scale	Period	Amplitude	Δχ2	kln N
b	356	2.5	32.6 days	4.6 minutes	13.9	17.6
c	121	1.3	49.7 days	0.9 minutes	8.8	14.4
d	81	1.6	44.5 days	9.2 minutes	31.9	13.2

Notes. Summary of the search for sinusoidal TTV signals in the transit times of Kepler-65b, c, and d. After correcting the timing error bars to account for correlated noise and low S/N, we find that only Kepler-65d has a significant detection.

Download table as:  ASCIITypeset image

To search for a sinusoidal TTV signal for each of the planets, we used the Bayesian information criterion. The criterion requires that when k new parameters are introduced in a model, one needs to achieve a decrease in χ² larger than kln N, where N is the number of data points, to justify the addition of the extra parameters. Table 2 shows that the additional three parameters are only justified for Kepler-65d, and not for the other two planets. That planet d is singled out in this test is consistent with the hypothesis that all three planets orbit Kepler-65, as this planet has the longest orbital period and is closest to the largest and presumably most massive planet (Kepler-65c; see Section 3.2 and Table 3), factors which enhance the amplitude of the TTV signal. Its orbital period is 1.39 times that of Kepler-65c, making it near a 7:5 ratio. This has been observed in many other multiplanet systems (Fabrycky et al. 2012).

Table 3. Transit Parameters

Parameter	Kepler-50b	Kepler-50c	Kepler-65b	Kepler-65c	Kepler-65d
(Rp/R)2 (ppm)	$99^{+5}_{-11}$	$159^{+10}_{-11}$	$85.0^{+1.6}_{-1.1}$	$282^{+5}_{-2}$	$98^{+2}_{-2}$
Impact parameter b	$0.74^{+0.07}_{-0.36}$	$0.94^{+0.02}_{-0.06}$	$0.16^{+0.20}_{-0.11}$	$0.42^{+0.10}_{-0.02}$	$0.53^{+0.07}_{-0.02}$
R/a	$0.095^{+0.015}_{-0.025}$	$0.084^{+0.013}_{-0.022}$	$0.188^{+0.011}_{-0.002}$	$0.097^{+0.006}_{-0.001}$	$0.078^{+0.005}_{-0.001}$
LD coefficient u1	$0.28^{+0.08}_{-0.08}$	⋅⋅⋅	$0.25^{+0.07}_{-0.08}$	⋅⋅⋅	⋅⋅⋅
LD coefficient u2	$0.29^{+0.08}_{-0.07}$	⋅⋅⋅	$0.37^{+0.11}_{-0.10}$	⋅⋅⋅	⋅⋅⋅
Transit duration T1.5–3.5 (hr)	$3.80^{+0.05}_{-0.03}$	$2.06^{+0.03}_{-0.03}$	$3.077^{+0.07}_{-0.007}$	$3.928^{+0.008}_{-0.007}$	$4.10^{+0.02}_{-0.03}$
Orbital period (days)	7.81254(10)	9.37647(4)	2.154910(5)	5.859944(3)	8.13123(2)
Time of transit (BJD−2,454,900)	74.376(7)	69.958(3)	66.4990(13)	65.0391(3)	70.9905(16)
Planet radius (RE)	$1.71^{+0.05}_{-0.10}$	$2.17^{+0.07}_{-0.08}$	$1.42^{+0.03}_{-0.03}$	$2.58^{+0.06}_{-0.06}$	$1.52^{+0.04}_{-0.04}$
Semi-major axis (AU)	$0.077^{+0.012}_{-0.020}$	$0.087^{+0.014}_{-0.023}$	$0.035^{+0.002}_{-0.001}$	$0.068^{+0.004}_{-0.002}$	$0.084^{+0.006}_{-0.002}$

Notes. Summary of planetary parameters. The first five parameters were estimated from the folded light curve analysis, and the durations were obtained from those parameters. Uncertainties come from an MCMC analysis. Orbital periods and times of transit come from a fit to the transit times, with a sinusoidal component in the case of Kepler-50 and only a linear term for Kepler-65. Estimation of planetary radii and semi-major axes also made use of the stellar radii from the asteroseismic analysis.

Download table as:  ASCIITypeset image

Even though the TTV signal of Kepler-65c did not satisfy the Bayesian information criterion for detection, the measured TTV period of the best-fitting sinusoid is close to the value expected from the formula of Agol et al. (2005), which for the c and d pair is 50 days (in agreement with the detected period). With more data one might be able to establish this signal more securely. Using the 0.9 minute amplitude of this hypothetical signal as a reference, one would estimate a mass (or upper bound; see Lithwick et al. 2012) for Kepler-65d of approximately 10 M_E.

We conclude that Kepler-65 is indeed transited by a system of three planets, based on the low FAP for each individual transit, the unlikely coincidence that would be required for a spurious system to produce the observed trend of transit durations versus orbital periods, the agreement between the photometric and asteroseismic estimates of the mean stellar density, and the detection of a physically reasonable TTV signal for at least one of the planets. However, as in many other cases, we acknowledge the fact that we cannot completely rule out the unlikely companion scenario, in which one or more of the planets orbits a fainter, bound companion having a higher density than Kepler-65.

3.2. Transit Parameters for Kepler-50 and Kepler-65

TTVs are quite significant for Kepler-50. For this reason, care was needed in producing a phase-folded transit light curve for subsequent analysis. In addition to the three-transit model parameters (R/a, b, (R_p/R)²) for each candidate and the two limb-darkening coefficients, we modeled the interval between transits as a constant plus a sinusoidal function of time. We fitted this model to SC data (Q6–Q11) for transits separated from each other by at least 6 hr and, using the best-fitting model, we folded the data and binned it to a cadence of 7.5 s. This template was then used to obtain the transit timings with uncertainties as described in the previous section. The new TTV signal, including the long-cadence (LC) timings from the Kepler catalog, was then fitted to improve the sinusoidal component of the ephemeris, which in turn was used to properly fold the data. This iterative process converged when the sinusoidal component did not change significantly from one step to the next.

The final phase-folded light curves and the best-fitting models are shown in Figure 4. The two planets have very similar orbital periods, with a period ratio close to 1.2. We see in the figure that the transits of the outer planet are much shorter in duration than those of the inner planet. This indicates that the transits of planet c have a high impact parameter. In this situation there is a risk of bias in the determination of the planet radius due to poorly constrained limb-darkening coefficients. To avoid this, we introduced Gaussian priors on each coefficient with values of 0.3  ±  0.1, based on the theoretical coefficients given by Claret & Bloemen (2011) for stars similar to Kepler-50. Since the orbits are so close to each other, we assumed circular orbits around the same star, essentially linking all of the R/a parameters for the two planets (see Section 3.1.2). The resulting planet parameters are given in Table 3. The stellar density derived from this transit model has a large uncertainty due to the low S/N of the transits, but the final value $0.40^{+0.6}_{-0.10}\,\rm g\,cm^{-3}$ is nevertheless compatible with the much more precisely determined asteroseismic density of 0.441 ± 0.004 g cm⁻³ (see Table 1 and Section 2).

Zoom In Zoom Out Reset image size

For Kepler-65 we used the analysis discussed in the previous section to construct the phase-folded light curves. Since no TTV signal was detected for planets b and c, a constant period was assumed in constructing the phase-folded light curves based on SC data. For planet d, the best-fitting sinusoidal TTV model was used to fold the transits. An individual analysis for each planet showed that the ingress duration of planet d was still larger than expected, by a factor of about two. This long ingress duration implies a large impact parameter, which again leads to large uncertainty and possible bias in the planetary radius. As for Kepler-50, we assumed the planets to be on circular orbits around the same star. The phase-folded light curves of all three transiting planets are shown in Figure 5, along with the best-fitting models.

Zoom In Zoom Out Reset image size

Final values of the planet parameters are presented in Table 3. With the improved folded light curve for planet d, the stellar mean density (assuming circular orbits) is $0.61^{+0.02}_{-0.10}\,\rm g\,cm^{-3}$, which also agrees with the value obtained from asteroseismology. One could use the asteroseismic density as a prior on our model, but this would not necessarily lead to greater accuracy because non-zero eccentricities cannot be ruled out for this system.

Finally, we note that in the transit analysis for both Kepler-50 and Kepler-65 we have assumed that the light from blended stars is negligible (i.e., a contamination factor of zero). This is well justified by the AO images (Adams et al. 2012). The contamination factors given in the Kepler Input Catalog are very low for both systems, and the uncertainties in the planetary radii are dominated by statistical uncertainties rather than the systematic effects of possible contamination.

3.3. Discussion of the Coplanarity of the Systems

4.1. Principles of the Method

Asteroseismic estimation of the stellar angle of inclination, i_s, rests on our ability to resolve and extract signatures of rotation in the non-radial modes from the oscillation spectrum. Detailed descriptions of the principles of the asteroseismic method may be found in Gizon & Solanki (2003) and Ballot et al. (2006, 2008). Here, we summarize the key points.

Rotation lifts the degeneracy in the oscillation frequencies ν_nl, so that the frequencies of non-radial modes (l > 0) depend on the azimuthal order, m. For the fairly modest rates of rotation typical of solar-like oscillators we may ignore, to first order, the effects of the centrifugal distortion (e.g., see Reese et al. 2006; Ballot 2010). The 2l + 1 rotationally split frequencies may then be written:

$$\begin{equation} \nu _{{nlm}} \equiv \nu _{{nl}} + \delta \nu _{{nlm}}, \end{equation} \tag{ 1 }$$

with

$$\begin{equation} \delta \nu _{{nlm}} \simeq \frac{m}{2\pi } \int _0^R \int _0^\pi K_{{nlm}}(r,\theta) \Omega (r,\theta)r\,{d}r\,{d}\theta. \end{equation} \tag{ 2 }$$

Here, Ω(r, θ) is the position-dependent internal angular velocity (in radius r, and co-latitude θ), and K_nlm is a weighting kernel that reflects the sensitivity of the mode to the internal rotation as a function of depth. For modest rates of differential rotation (in latitude and radius) and absolute rotation, the splittings δν_nlm of the observable high-n, low-l p-modes will take very similar values, hence tending to the approximation of solid-body rotation (Ledoux 1951). Here, we found no evidence for significant mode-to-mode variation of the frequency splittings in the oscillation spectrum of either star. In what follows we therefore modeled all splittings as being equal, i.e., δν_nlm = δν_s. The above also neglects any contributions to the splittings from near-surface magnetic fields, which give rise to frequency asymmetries of the observed splittings. The levels of activity in both stars—as revealed by signatures of rotational modulation of spots and active regions in the Kepler light curves—are notably lower than those displayed by the active Sun (see later, in Section 5). Since magnetic contributions to the solar low-l splittings are small in size and very hard to measure in Sun-as-a-star data of much higher S/N (e.g., see Gough & Thompson 1990; Chaplin 2011, and references therein), asymmetries here should not be a cause for concern for the analysis.

The determination of the inclination of the stellar rotation axis relies on the fact that the mode patterns of the non-radial modes are not spherically symmetric. The disk-integrated amplitudes of the m components in any given non-radial multiplet will therefore depend on the viewing angle. Figure 6 shows a snapshot of the intensity perturbations of the m = 1 (left-hand column) and m = 0 (right-hand column) components of an l = 1 mode viewed at different angles, i_s. The perturbations are shown at a phase corresponding to extreme displacement of each oscillation mode. The filled circles mark the pole of the rotation axis and the lines show the stellar equator. Note that m = −1 perturbations are π out of phase with the m = 1 (and have not been plotted here). When the rotation axis lies in the plane of the sky (i_s = 90°), the m = ±1 components presents their strongest observable amplitudes. In contrast, the m = 0 component cannot be detected because the intensity perturbations in the northern and southern hemispheres cancel at all phases of the pulsation cycle, giving no disk-averaged signal. The situation is reversed at i_s = 0°, when the rotation axis lies along the line of sight and perturbations due to the m = ±1 components are no longer visible owing to geometric cancellation.

Zoom In Zoom Out Reset image size

This dependence (measured in power) may be written explicitly as

$$\begin{equation} {\cal E}_{lm}(i_{\rm s}) = \frac{(l-|m|)!}{(l+|m|)!} \left[ P_l^{|m|} (\cos i_{\rm s}) \right]^2, \end{equation} \tag{ 3 }$$

where $P_l^{|m|}$ is the Legendre function, and the sum over ${\cal E}_{lm}(i_{\rm s})$ is normalized to unity. Measuring the relative power of the azimuthal components of different |m| in a non-radial multiplet therefore provides a direct estimate of the stellar angle of inclination, i_s, or more properly |i_s| since symmetries inherent in Equation (3) mean we cannot discriminate between i_s and −i_s, and π − i_s, and π + i_s.

The above discussion rests on two assumptions. First, that contributions to the observed stellar intensity across the visible stellar disk depend only on the angular distance from the disk center. This is valid for photometric observations, where limb darkening controls the weighting. Second, that there is equipartition of energy between the different m components.¹⁹ This should be valid except in very rapid rotators where rotation can affect convection, which excites and damps the modes. The predicted power asymmetries (Belkacem et al. 2009) of our stars are of the order of 1%, which are negligible for our analysis.

Figure 7 shows the appearance in the frequency-power spectrum of an idealized l = 1 multiplet as a function of the angle i_s (see also Gizon & Solanki 2003). Panels in the left-hand column correspond to the cases shown in Figure 6. The l = 1 modes are approximately three times more prominent in the frequency-power spectrum than the l = 2 modes (see, e.g., Ballot et al. 2011a). Hence, it is these modes that largely constrain our ability to infer i_s.

Zoom In Zoom Out Reset image size

The individual components in Figure 7, which are plotted in blue (m  =   ± 1) and red (m = 0), were modeled as Lorentzian functions, the underlying function used to describe the damped p-modes. The width Γ of each Lorentzian—which is proportional to the mode damping rate—is 3.0 μHz, which corresponds approximately to the linewidths observed in the most prominent l = 1 modes of Kepler-50 and Kepler-65. The splitting is δν_s = 1.5 μHz, which corresponds to a rotation period of 7.7 days, and so matches approximately what we observe for the two stars. The thick black lines show the combined multiplet profiles.

Given sufficient resolution in frequency and good S/N in the modes, it is the ratio δν_s/Γ of intrinsic stellar properties that determines whether it is possible to resolve the components, hence to infer the true underlying ${\cal E}_{lm}(i_{\rm s})$ and hence the value of i_s. As noted previously, at angles close to 90° the m = 0 component has insignificant visibility and the overall appearance is dominated by the |m| = 1 components (the converse being true at angles close to 0°). Uncertainties in the inferred angle will be largest when the i_s matches these extreme cases, all other factors being equal (Ballot et al. 2008). This is because there are then only modest variations in the overall appearance of the mode multiplet with changing i_s.

4.2. Estimation of Stellar Inclination Angles

Extracting the required information from the rotationally split components proceeds via a careful fitting of the modes in the observed frequency-power spectrum, sometimes referred to as peak-bagging (see Appourchaux et al. 2012, and references therein, for further results on Kepler targets).

Frequency splittings δν_s due to rotation are clearly visible in the oscillation spectra of both stars. Figure 8 shows two prominent l = 1 modes in each star. The light gray lines plot the observed spectra after applying a light amount of smoothing. The thick dark gray lines follow the spectra after they have been smoothed with a filter of width 1.5 μHz, which provides an approximate representation of the underlying (noise-free) profiles. Even without a detailed analysis it is apparent that the observed modes bear a striking resemblance to the high-inclination cases in Figure 7. The dark-blue lines follow the best-fitting Lorentzian models, which we describe below.

Zoom In Zoom Out Reset image size

The observed frequency-power spectrum $\mathcal {P}(\nu)$ of each star was modeled as

$$\begin{equation} \mathcal {P}(\nu) = \mathcal {O}(\nu) + B(\nu), \end{equation} \tag{ 4 }$$

where $\mathcal {O}(\nu)$ describes the oscillations and B(ν) contains background terms due to granulation, activity, and photon shot noise. The oscillations $\mathcal {O}(\nu)$ were modeled as a series of Lorentzian profiles describing the stochastically excited and intrinsically damped modes. We adopted a global description, in which we modeled simultaneously all the observable modes instead of modeling and analyzing the spectrum one order at a time. This approach improves the accuracy of the modeling because it takes proper account of the power from the slowly decaying Lorentzian peaks that bleeds in frequency between the neighboring modes. The modeled oscillations' spectrum was thus described by

$$\begin{equation} \mathcal {O}(\nu) = \sum _{n^{\prime },l} \sum _{m=-l}^{l} \frac{{\cal E}_{lm}(i_{\rm s}) H_{n^{\prime }l} }{1 + 4/\Gamma _{n^{\prime }}^2 (\nu - \nu _{n^{\prime }l} - m \delta \nu _{\rm {s}})^2}, \end{equation} \tag{ 5 }$$

The inner sum in the above runs over the m components of each rotationally split multiplet, while the outer sum runs over all observed modes, in radial order n, and degree l. Note that the dummy variable n', which tags the radial order, takes values n' = n for l = 0 and l = 1 modes, and n' = n − 1 for l = 2 modes (which lie adjacent in frequency to l = 0 modes of n' = n). The angle i_s and single splitting parameter δν_s are two of the parameters to be optimized, along with the frequencies ν_itnl used to estimate the fundamental stellar properties (Section 2).

The parameters H_nl and Γ_n describe the height (maximum power spectral density) and linewidth of each mode. We fit a single linewidth parameter to each order. The relative heights of the components in each non-radial multiplet are controlled by i_s through the function ${\cal E}_{lm}(i_{\rm s})$ (Equation (3)). The heights are constrained by the relation $H_{n^{\prime }l} = H_{n^{\prime }0} V^2_{l}$, where the parameter $V^2_{l}$ describes the visibilities of modes of different l, relative to l = 0. The visibilities are given by integrating the spherical harmonic functions over the visible disk with suitable allowance made for limb darkening and the spectral bandpass of the observations (Ballot et al. 2011a). We adopted fixed values of $V^2_0=1.0$, $V^2_1=1.5$, and $V^2_2=0.5$ in our analysis.

The background was modeled as the sum of three components: a flat photon shot-noise component, W, and two frequency-dependent components to describe contributions from granulation and activity. For the latter components, we used functions based on the Lorentzian-like forms proposed by Harvey (1985), which provide a good description of the observed backgrounds in solar-type stars (e.g., see Metcalfe et al. 2012). The composite background was then described by

$$\begin{equation} B(\nu) = W + \sum _{k=1}^2\frac{4 \sigma _k^2 \tau _k}{1 + (2\pi \tau _k\nu)^2 + (2\pi \tau _k\nu)^4}, \end{equation} \tag{ 6 }$$

with k = 1 associated with the granulation component, and k = 2 associated with the activity component. The granulation and activity components each have two free parameters to be optimized: σ_k is related to the rms amplitude of the signal in the time domain, while τ_k is the characteristic timescale of the decaying autocorrelation function. For granulation, σ and τ are smaller than the corresponding activity parameters.

We adopted two different approaches to the fitting, and hence to estimation of i_s. In the first approach the parameters of the model in Equation (4) were optimized using a MCMC approach, as described by Handberg & Campante (2011). We adopted a flat prior for i_s between 0° and 90°. In order to avoid the rejection of sample jumps close to the boundaries—i.e., those that would jump beyond the range set by the prior—in practice we selected from samples in the range −90° to 180° and modified accepted jumps that went beyond the allowed range by reflecting about the i_s = 0° and 90° boundaries. A flat prior was imposed on δν_s, running between 0 and 2.5 μHz for Kepler-50, and 0 and 5 μHz for Kepler-65. Besides making it possible to incorporate relevant prior information through Bayes' theorem, the MCMC approach also gave the marginal probability density function (PDF) of each of the model parameters (e.g., see discussion in Appourchaux 2011). In order to provide a cross-check we also used maximum likelihood estimation (MLE) to fit the spectrum (e.g., see Duvall & Harvey 1986; Toutain & Appourchaux 1994), using the so-called pseudo-global fitting described by Fletcher et al. (2009). Rather than fit the i_s directly with MLE, we instead performed a series of MLE fits with the angle fixed at different values, the aim being to sample the maximized likelihood of the best-fitting model as a function of the chosen i_s. Even though i_s is independent of the splitting parameter δν_s, the measured values are often highly correlated (e.g., see Ballot et al. 2006, 2008). In order to constrain the two individual parameters, or their combination the so-called reduced splitting (i.e., δν_s sin i_s), it is desirable to have access to the corresponding maximized likelihood in two-dimensional parameter space. We therefore performed fits with both i_s and δν_s taking values on a dense grid. This yielded a two-dimensional grid of maximized likelihoods, making possible inference on the inclination and splitting from construction of confidence intervals based on the likelihood surface. The MLE approach had the advantage of being more computationally efficient than the MCMC analysis. However, given that the input values for the inclination and splitting are fixed prior to the fitting, one cannot extract a bona fide posterior probability distribution. The MCMC and MLE approaches returned results in excellent agreement. Here, we present those given by the MCMC approach.

Table 4 lists the final MCMC estimates of the inclinations and splittings, together with their corresponding 1σ credible regions. The estimated i_s of both stars are consistent with 90°, to within the uncertainties. We note that the final values for the splittings were given by the median of each posterior distribution, while for the angles we opted to use the mode of each distribution. The rationale behind this decision was twofold. First, the PDF of the inclination is truncated at i_s = 90° and the median is thus not a representative statistic of the distribution. Second, in each case the model of the oscillations' spectrum built by using the mode for the inclination, together with the median for all the remaining parameters (including the splitting), has a higher posterior likelihood than the models made using exclusively either the median or the mode for all parameters.

Table 4. Estimated Stellar Inclinations and Rotational Splittings

Star	is	sin  i	δνs ≡ Ω/2π
	(deg)		(μHz)
Kepler-50	$82^{+8}_{-7}$	$0.99^{+0.01}_{-0.02}$	$1.51^{+0.09}_{-0.08}$
Kepler-65	$81^{+9}_{-16}$	$0.99^{+0.01}_{-0.08}$	$1.30^{+0.19}_{-0.16}$

Download table as:  ASCIITypeset image

As noted above, the thick black lines in Figure 8 show best-fitting models from the MCMC analysis across frequency ranges occupied by two l = 1 modes in each star. Figures 9 and 10 show the correlation maps in the angle and splitting, as well as the PDFs obtained after marginalization. Binwidths for the PDFs were fixed by following the procedure given in Scott (1979; see also Handberg & Campante 2011).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have compared the asteroseismic results from Section 4.2 with two independent estimates of the surface rotation: one extracted from signatures of rotational modulation in the Kepler light curves, and another extracted from spectroscopic data on both stars.

If a star has spots on its surface then rotation will carry the spots in and out of view, inducing quasi-periodic flux variations. Such variations have been detected for many stars, and it is not unusual to see activity in stars as hot as our host stars (Basri et al. 2011). The rotation period of Kepler-50 has already been detected in Kepler data (Hirano et al. 2012b) and Kepler-65 also shows clear signs of rotational modulation in its light curve.

The raw Kepler data are known to suffer from systematic trends that appear to be shared by most of the stars on a given CCD detector module. The absolute effect of these trends on the measured stellar fluxes is much larger than the activity levels for both stars, so we needed to choose an appropriate detrending algorithm that would suppress the unwanted systematics without removing the astrophysical signal. One such algorithm available to us is PDC-MAP (Stumpe et al. 2012; Smith et al. 2012), which was developed by the Kepler team. Principal component analysis is used to extract a basis of co-trending vectors that capture the systematic trends in each module. Each stellar flux datum may be decomposed into a linear combination of the astrophysical variability and the co-trending vectors. One could perform a least-squares fit to extract those coefficients but PDC-MAP goes one step further. During a first pass on the data it applies a least-squares approach to obtain the coefficients for all stars that fall on a given detector module. The distributions are then used as priors in a Bayesian sense, thereby mitigating possible overcorrection of the time series for any individual star. PDC-MAP appears to do a reasonable job of eliminating the systematic trends in almost every quarter of data, but does fail on a few occasions. From the available LC data (Jenkins et al. 2010) collected through Q11, we discarded the Q3 and Q11 data for Kepler-50, and the Q7 and Q8 data for Kepler-65. Outliers were also removed by performing 3σ clipping using a 10 hr long moving-median filter. The remaining data were quite adequate to estimate robust rotation periods. The left-hand panels of Figure 11 show three-month segments of the data, in which intrinsic stellar variability on timescales of days is evident. This variation is smaller than that displayed by the active Sun (e.g., see Basri et al. 2010, 2011), for which the semi-amplitude of the variability reaches levels close to 1 part in 10³.

Zoom In Zoom Out Reset image size

To extract estimates of the surface rotation periods we followed the analysis described by Hirano et al. (2012b). We first calculated the Lomb–Scargle periodogram of each set of processed data, which are plotted in the right-hand panels of Figure 11. Both stars show significant peaks around 8 days due to rotational modulation of spots. We checked that our results did not depend on the detrending, as follows. First, we used the simple least-squares fit to the co-trending vectors described above. Second, we applied a 20 day median smoothing filter to divide the raw data. Both approaches led to periodograms in good agreement with our main results.

Both stars show several significant, closely spaced peaks in their periodograms. The spread in period of the peaks will have contributions from the finite spot lifetimes and might also suggest the presence of surface differential rotation. Figure 11 also shows periodograms obtained from particular subsets of the data, which suggest that the rotation periods might be evolving with time. Again, this might be associated with stochastic variability due to the spot lifetimes, or it could have a contribution due to changes in the spot latitudes. These signatures are not unexpected for such hot stars (Collier Cameron 2007) but the low S/N of the stellar variability in the light curves makes it hard to test the results in greater detail. It is worth adding that we also checked that the spread of significant periods was not an artifact of the observational window function. We sampled commensurate sine waves of period 8 days at the same time stamps as the real, cleaned PDC-MAP light curves, and found that the resulting periodograms displayed a much smaller width of significant periods than the real data (of order 0.2–0.4 days).

As in Hirano et al. (2012b), for each star the continuous range of periods containing all peaks with more than half the power of the highest peak was taken as the uncertainty on the surface rotation period, with the center of the range adopted as the quoted period, P_rot. The rotation periods obtained in this way were 8.4 ± 1.0 days for Kepler-50 and 8.4 ± 0.3 days for Kepler-65. The rotational frequency splittings δν_s given by the asteroseismic analysis are equivalent to rotation periods P_rot = 1/δν_s of $7.7^{+0.5}_{-0.4}\,\rm days$ for Kepler-50 and $8.9^{+1.3}_{-1.1}\,\rm days$ for Kepler-65. The close agreement between the surface and internal rotation rates is consistent with the expectation that the frequency splittings of p-modes observed in main-sequence stars are largely determined by the rotation profile in the stellar envelope (as in the case of the Sun).

We can also compare our results with vsin  i_s estimates extracted from the spectroscopic data (see Section 2). For Kepler-50 the spectroscopic analyses gave 8.6 ± 0.8 km s⁻¹ (SME) and 10.3 ± 0.5 km s⁻¹ (SPC), while for Kepler-65 they gave 9.8 ± 0.8 km s⁻¹ (SME) and 11.9 ± 0.5 km s⁻¹ (SPC). To compare with the asteroseismic results, we converted the estimated rotational frequency splittings to projected rotational velocities using

$$\begin{equation} v \sin \,i_{\rm s} \equiv 2\pi R\,\delta \nu _{\rm s}\sin \,i_{\rm s}, \end{equation} \tag{ 7 }$$

with the stellar radii R given by the asteroseismic analysis discussed in Section 2. These conversions gave equivalent asteroseismically determined projected velocities of $8.0^{+1.2}_{-1.0}\,\rm km\,s^{-1}$ for Kepler-50 and 10.4 ± 0.6 km s⁻¹ for Kepler-65. Again, we find agreement between the asteroseismic and surface estimates. We note that the SPC vsin  i_s values are higher than the SME values by about 20%, in agreement with the findings by Torres et al. (2012).

Finally, we may combine the spot modulation and vsin  i_s results to provide independent estimates of i_s, via

$$\begin{equation} \sin \,i_{\rm s} = P_{\rm rot} (v\sin \,i_{\rm s})/(2\pi R). \end{equation} \tag{ 8 }$$

When we used the SME results for vsin  i_s we found that sines of the angles were constrained at the 1σ level to lie above 0.89 for Kepler-50 and 0.90 for Kepler-65, again implying that both stars have their rotation axes nearly perpendicular to the line of sight. When the SPC vsin  i_s were used we obtained sin  i_s > 1.0 for both systems suggesting that (at least for these systems) the SPC results are overestimated.

Our central result is that the host stars of the Kepler-50 and Kepler-65 planetary systems have their rotation axes nearly perpendicular to the line of sight, with sin i_s constrained at the 1σ level to lie above 0.97 and 0.91, respectively. Expressed in terms of angles, we have |90° − i_s| < 15° for Kepler-50 and <25° for Kepler-65. The orbital inclinations of the planets in each system are also near 90°, with a deviation of only ≈5° for the planets of Kepler-50 and ≈2° for the planets of Kepler-65. Therefore our observations are consistent with small differences in the stellar and orbital inclination angles, and low stellar obliquities.

A limitation of the results is that it is possible for the obliquity to be large but for the difference in inclination angles to be small. Spherical geometry dictates that the three-dimensional obliquity angle ψ between the stellar spin axis and the planetary orbital axis is given by

$$\begin{equation} \cos \psi = \cos i_{\rm s} \cos i_{\rm p} + \sin i_{\rm s} \sin i_{\rm p} \cos \lambda , \end{equation} \tag{ 9 }$$

where λ is the projected angle on the sky between the orbital and rotational angular momentum vectors, which the asteroseismic method cannot provide. This is the converse of the situation with, for example, the RM and spot-occultation methods, for which the data reveal λ but are generally insensitive to i_s. To overcome this limitation, one would need to combine the various measurement techniques to determine both λ and i_s for the same system, or else conduct observations of many systems and perform a statistical analysis of the ensemble (see, e.g., Fabrycky & Winn 2009). Measuring the RM effects given by the planets in the two systems considered in this paper would be very challenging, due to the small sizes of the planets. The highest-amplitude RM signal for Kepler-50 would be around ≃ 0.3 m s⁻¹ (given by Kepler-50c), while for Kepler-65 it would be ≃ 1.8 m s⁻¹ (given by Kepler-65c).

While the two stars we have analyzed are special—in the sense that they have transits—it is important to test the null hypothesis that they are members of a population of stars randomly oriented in space. Put another way: had the stars been randomly oriented, what is the chance that the inclination angles would have been as close to 90°, as observed? For an isotropic distribution, the probability of observing |90° − i_s| < x is simply sin x. Evaluating this for Kepler-50 and Kepler-65 we find the chance to be 26% and 42% individually, and the chance of observing both of them so close to 90° is 11%. In this sense there is only an 11% chance we would have obtained our results in the absence of any correlation between the orientations of the stellar rotation and planetary orbits.

Despite these limitations, let us consider the implications of low obliquities. The number of extrasolar multiplanet systems for which the stellar obliquity has been measured is now four, with Kepler-50 and Kepler-65 joining the previously studied systems Kepler-30 (Sanchis-Ojeda et al. 2012) and KOI-94 (Hirano et al. 2012a). In all four cases, the obliquity is consistent with zero. We have already summarized our results for Kepler-50 and Kepler-65, and for the other two systems λ was found to be consistent with 0 to within about 10°.

Several years ago these results would have seemed mundane. Low obliquities are expected in the standard picture in which the star and planets have the same direction of angular momentum originating from a common accretion disk. Furthermore, up until 2008, low obliquities had been observed in all the exoplanetary systems that had been examined (all involving single hot Jupiters). Since that time, systems with hot Jupiters have been found to have host stars with a wide range of obliquities (Albrecht et al. 2012). Those results have been taken as evidence that the process that produces hot Jupiters also tilts their orbits relative to the initial plane of their formation. Specifically, the results have been taken to support theories for the origin of hot Jupiters involving few-body dynamics and tidal circularization (alternatively referred to as high-eccentricity migration) as opposed to the formerly dominant paradigm of disk migration.

It has been pointed out, however, that this conclusion depends critically on the assumption that the current stellar equatorial plane is aligned with the original plane of the planetary orbits. Several papers have questioned this assumption. Bate et al. (2010) suggested that the chaotic environment of a star-forming region might lead to large mismatches between the direction of stellar rotation and the orientation of the late-stage protoplanetary disk. Thies et al. (2011) proposed that inclined planets arise from capture of gas from a neighboring star. Lai et al. (2011) presented a theory of magnetic interaction between a young star and the inner edge of its accretion disk that can tip the star by a significant angle. Rogers et al. (2012) considered stars with convective cores and radiative envelopes, and found that they might be susceptible to a directional wandering of the outermost layers of the star due to transport of angular momentum by internal gravity waves from the convective–radiative boundary.

In a few cases, it has been possible to compare the orientation of stellar rotation and the orientation of a surrounding disk (Le Bouquin et al. 2009; Watson et al. 2011), and the results have favored the hypothesis of close alignment and low obliquities. However, more stringent tests are warranted, and are provided by the results for the four multiplanet systems Kepler-30, 50, 65, and KOI-94. The low obliquities suggest that high obliquities are confined to the hot-Jupiter systems, or at least have not yet falsified that hypothesis. The results are therefore starting to provide support to the argument that most or all hot Jupiters are formed through inclination-lifting processes and not via disk migration.

All four of the stars are cool enough to have outer convective zones (indeed, Kepler-50 and Kepler-65 must have convection because they exhibit solar-like oscillations). It is unclear how effective the model of Rogers et al. (2012) might be for stars with thin convective envelopes (e.g., like Kepler-50, Kepler-65, and KOI-94, which are all hotter than the Sun). It is worth noting that this theory already has difficulty accommodating the hot-Jupiter results involving cool, convective stars with high obliquities, namely, HAT-P-11 (Winn et al. 2010a; Sanchis-Ojeda & Winn 2011), HD 80606 (Winn et al. 2009; Pont et al. 2009; Hébrard et al. 2010), and WASP-8 (Queloz et al. 2010).

Of course, one should not be satisfied with a sample of only four systems, especially given the limitations caused by projection effects as noted above. The asteroseismic technique that was deployed in this work has the advantage that the detectability of the signal is chiefly a function of the stellar properties, as opposed to the planetary properties, and therefore has no intrinsic difficulty with small planets or long-period planets. We expect it will be possible to apply this technique to a sample of at least 10 other Kepler systems in the near future. It will be possible to draw stronger conclusions with these results in hand. Definitive conclusions will also be possible on individual systems with transiting exoplanets when asteroseismology demonstrates that i_s is significantly different from 90°.

Funding for this Discovery mission is provided by NASA's Science Mission Directorate. The authors wish to thank the entire Kepler team, without whom these results would not be possible. W.J.C., T.L.C., G.R.D., Y.E., and A.M. acknowledge the support of the UK Science and Technology Facilities Council (STFC). J.N.W. was supported by the NASA Kepler Participating Scientist program through grant NNX12AC76G. S.B. acknowledges NSF grant AST 1105930. Funding for the Stellar Astrophysics Centre (SAC) is provided by The Danish National Research Foundation. The research is supported by the ASTERISK project (ASTERoseismic Investigations with SONG and Kepler) funded by the European Research Council (grant agreement no.: 267864). S.H. acknowledges financial support from the Netherlands Organisation for Scientific Research (NWO). Computational time on Kraken at the National Institute of Computational Sciences was provided through NSF TeraGrid allocation TG-AST090107. We acknowledge the Pale Blue Dot Project, hosted by White Dwarf Research Corporation (whitedwarf.org/palebluedot), and we are also grateful for support from the International Space Science Institute (ISSI).

In the second stage of the stellar properties estimation (Section 2) three members of the team (S.B., J.C.D., and T.M.) performed a detailed modeling of the stars using estimates of the individual oscillation frequencies and the revised spectroscopic data as inputs.

S.B. made use of the Yale stellar evolution code, YREC (Demarque et al. 2008) to model both stars. The input physics included the OPAL equation of state tables of Rogers & Nayfonov (2002), and OPAL high-temperature opacities (Iglesias & Rogers 1996) supplemented with low-temperature opacities from Ferguson et al. (2005). All nuclear reaction rates were from Adelberger et al. (1998), except for the rate of the ¹⁴N(p, γ)¹⁵O reaction, which was fixed at the value of Formicola et al. (2004). Models were constructed for two values of core overshoot, 0 and 0.2H_p. Two families of models were constructed, one that included the diffusion and settling of helium and heavy elements as per the formulation of Thoul et al. (1994), and one that did not include any diffusion and settling.

YREC was used in an iterative mode whereby the final T_eff and radius for a star of a given mass and metallicity was specified, and for a given mixing length parameter α the code iterated over the initial helium abundance Y₀ until a model with the specified T_eff and radius was found. This is similar to the construction of a standard solar model, although in the solar case iterations are made over both the mixing length parameter and Y₀ with solar age a fixed independent constraint. Since the ages of the Kepler stars are not known independently, iteration over Y₀ were performed for many different values of the mixing length parameter. All solutions for which the initial helium abundance was less than the primordial helium abundance, Y_p were rejected. Y_p was assumed to be 0.245.

Corrections for near-surface effects (the so-called surface term) were handled in the following manner. The first step was the construction of a standard solar model with exactly the same physics as that used to model the Kepler stars. This yielded the set ν_nl☉ of solar model frequencies. These were then used to estimate a set of "surface term" frequency offsets, δν_nl☉, for the Sun by computing differences between the solar model frequencies and the solar low-degree frequencies observed by the Birmingham Solar Oscillations Network (BiSON; as listed in Basu et al. 2009).

For each stellar model, ${\cal M^{\prime }}$, ν_nl☉, and δν_nl☉ were then scaled to the mass and radius of ${\cal M^{\prime }}$ using the homology scaling $r = \langle \Delta \nu ({\cal M^{\prime }}) \rangle / \langle \Delta \nu _{{nl}\odot } \rangle$. The resulting rν_nl☉−rδν_nl☉ relation was then used to correct the stellar model for the surface term. Using a least-squares minimization a factor β was selected so as to minimize $\sum (\nu ^{\rm obs}_{{nl}}-\nu ^{\rm corr}_{{nl}} )^2/ (\sigma ^{\rm obs}_{{nl}} )^2$ over all observed modes, where $\nu ^{\rm corr}_{{nl}}=\nu _{{nl}}^{\cal M^{\prime }}+ \beta \;r\delta \nu _{{nl}\odot }$, with rδν_nl☉ evaluated at $r\nu _{{nl}\odot }=\nu ^{\rm obs}_{{nl}}$.

J.C.D. followed a prescription that has previously been applied to the Hubble observations of HD 17156 (Gilliland et al. 2011), and several Kepler exoplanet host stars, i.e., HAT-P-7 (Christensen-Dalsgaard et al. 2010), Kepler-10 (Batalha et al. 2011), and Kepler-36 (Carter et al. 2012). Stellar-evolutionary models were computed with the ASTEC code (Christensen-Dalsgaard 2008a). The calculations used the OPAL equation of state tables (see Rogers et al. 1996) and OPAL opacities at temperatures above 10⁴ K (Iglesias & Rogers 1996); at lower temperature the Ferguson et al. (2005) opacities were used. Nuclear reactions were calculated using the NACRE parameters (Angulo et al. 1999). Convection was treated using the Böhm-Vitense (1958) mixing-length formulation. Frequencies were computed for the models using ADIPLS (Christensen-Dalsgaard 2008b) and then corrected for surface effects following the prescription of Kjeldsen et al. (2008).

For each evolutionary sequence in the grid of ASTEC models, the model ${\cal M}^{\prime }_{\rm min}$ whose surface-corrected frequencies provided the best χ² match to the observations was selected. The best match was obtained from application of homology scaling, under the assumption that in the vicinity of ${\cal M}^{\prime }_{\rm min}$ frequencies could be calculated using $r\nu _{{nl}}({\cal M}^{\prime }_{\rm min})$, where $r=[R/R({\cal M}^{\prime }_{\rm min})]^{-1.5}$, R being the surface radius of the model. A best-fitting model was then determined by minimizing the sum $\sum (\nu ^{\rm obs}_{{nl}} - r\nu _{{nl}}({\cal M}^{\prime }_{\rm min}) )^2/ (\sigma ^{\rm obs}_{{nl}} )^2$ over all observed modes, as a function of r. The resulting minimum value of r defined an estimate of the radius of the best-fitting model along the given sequence. The other properties of this best-fitting model were determined by linear interpolation in R, to the minimum of R. Statistical analysis of the ensemble of best-fitting properties from all evolutionary sequences then yielded the final stellar properties, and their uncertainties (see Christensen-Dalsgaard et al. 2010; Carter et al. 2012).

T.M. used the Asteroseismic Modeling Portal (AMP), a web-based tool tied to TeraGrid computing resources that uses a parallel genetic algorithm (Metcalfe & Charbonneau 2003) to optimize, in an automated manner, the match to observational data (see Metcalfe et al. 2009; Woitaszek et al. 2009 for more details). AMP employs the Aarhus stellar evolution code ASTEC (Christensen-Dalsgaard 2008a) and adiabatic pulsation code ADIPLS (Christensen-Dalsgaard 2008b). Models were made using the OPAL 2005 equation of state and the most recent OPAL opacities supplemented by Ferguson et al. (2005) opacities at low temperature, nuclear reaction rates from NACRE (Angulo et al. 1999), and helium diffusion and settling following Michaud & Proffitt (1993). Convection was treated with standard mixing-length theory without overshooting (Böhm-Vitense 1958).

Each AMP model evaluation involved the computation of a stellar evolution track from the zero-age main sequence through a mass-dependent number of internal time steps, terminating prior to the beginning of the red giant stage. The asteroseismic age was optimized along each evolutionary track using a binary decision tree under the assumption that <Δν > is a monotonically decreasing function of age (see Metcalfe et al. 2009, and references therein). The Kjeldsen et al. (2008) prescription was again applied in an attempt to deal with the surface term. The optimal model was then subjected to a local analysis that uses singular value decomposition to quantify the uncertainties of the final model parameters (see Creevey et al. 2007).