SPIN–ORBIT ALIGNMENT OF EXOPLANET SYSTEMS: ENSEMBLE ANALYSIS USING ASTEROSEISMOLOGY

The angle ψ between the planetary orbital axis and the stellar spin axis (the true obliquity or spin–orbit angle) is a fundamental geometric property of planetary systems. Furthermore, it has been recognized as an important diagnostic of theories of planet formation, migration, and tidal evolution. For all these reasons, seeking empirical constraints on ψ is a matter of the utmost importance.

In the case of an exoplanet found through the radial-velocity (RV) method, no information is available about the degree of spin–orbit alignment. For transiting exoplanets, on the other hand, the Rossiter–McLaughlin (RM) effect has now been widely exploited (e.g., Queloz et al. 2000; Winn et al. 2005, 2009, 2010b, 2011; Hébrard et al. 2008; Triaud et al. 2010; Hirano et al. 2011; Albrecht et al. 2012, 2013). This technique is sensitive to the angle λ between the sky-projected orbital and spin axes (the projected spin–orbit angle). At the time of writing there are 87 planets with published measurements of the RM effect (see the online compilation¹⁵ by R. Heller), of which 36 (∼40%) show substantial misalignments according to at least one publication (see also Figure 7 of Xue et al. 2014). Other techniques that allow obliquity measurements of transiting systems include the analysis of planetary transits over starspots (e.g., Désert et al. 2011; Nutzman et al. 2011; Sanchis-Ojeda et al. 2011, 2012), Doppler tomography (e.g., Collier Cameron et al. 2010; Gandolfi et al. 2012; Albrecht et al. 2013; Johnson et al. 2014), the analysis of the effect of gravity darkening on the transit light curve (e.g., Barnes 2009; Barnes et al. 2011; Ahlers et al. 2014), and the analysis of the photometric amplitude distribution of stellar rotation (Mazeh et al. 2015).

Most obliquity measurements to date have been for systems harboring hot Jupiters, owing to the fact that the RM effect scales as the planet-to-star area ratio and to the increased opportunities for obtaining follow-up spectroscopic observations due to the frequent transit events. Empirical evidence has been found that the obliquities of hot-Jupiter systems are affected by tidal evolution (Schlaufman 2010; Winn et al. 2010a; Morton & Johnson 2011; Triaud 2011; Albrecht et al. 2012; although see Mazeh et al. 2015 for evidence against this theory): systems expected to undergo strong planet–star tidal interactions are preferentially found to have low obliquities, while systems with weaker tidal interactions display a wide range of obliquities that, besides well-aligned planets, also include planets in polar or even retrograde orbits. This suggests that the orbital plane has changed relative to the plane of the protoplanetary disk by the time hot Jupiters are formed and before tides have had any impact on shaping the system, which presumably happens due to the same mechanism responsible for their migration.

The above discussion assumes that the protoplanetary disk is coplanar with the stellar equator. The possibility remains, however, that primordial star–disk misalignments are ubiquitous, meaning that large obliquities could be a generic feature of planetary systems and not specific to hot-Jupiter formation. This hypothesis may in principle be tested by measuring the obliquities of systems with multiple transiting planets, since the planetary orbits in these systems are nearly coplanar and presumably trace the plane of the protoplanetary disk (Lissauer et al. 2011; Fabrycky et al. 2014). Accordingly, if multi-transiting systems tend to have low obliquities, then the high obliquities observed for hot-Jupiter systems are likely to be associated with planet migration. If, instead, the obliquity distribution of multi-transiting systems is similar to that of hot-Jupiter systems, then this would indicate more general processes of stellar and planet formation: processes such as chaotic star formation (e.g., Bate et al. 2010; Thies et al. 2011; Fielding et al. 2015), magnetic star–disk interactions (e.g., Lai et al. 2011), torques due to internal gravity waves (e.g., Rogers et al. 2012), or torques due to neighboring stars (e.g., Batygin 2012).

In order to study the dynamical histories of planetary systems across a comprehensive range of architectures and in a variety of environments, it is imperative to extend obliquity measurements to systems with smaller planets, longer-period planets, and multiple planets (note that, according to the current state of knowledge, hot Jupiters rarely have nearby planetary companions and may thus occur preferentially as single planets; Steffen et al. 2012). For long-period planets, however, the opportunities to observe transits occur less frequently, which limits the possibility of obtaining follow-up observations or identifying the transit geometry from starspot crossings. Furthermore, application of the RM technique becomes increasingly more challenging when dealing with multiple-planet systems, since these systems also tend to involve smaller planets (e.g., Latham et al. 2011).

An alternative technique for measuring the obliquities of planetary systems, one that does not depend on the signal-to-noise ratio (S/N) of the transit data and hence on planet size, makes use of a combination of the photometric stellar rotation period, P_rot, and the spectroscopically determined projected rotational velocity, $v\mathrm{sin}{i}_{{\rm{s}}}$, and stellar radius, R_s, to determine the sine of the angle i_s between the stellar spin axis and the line of sight (the stellar inclination angle). This technique evolved from the technique originally developed by Herbst et al. (1986) and Hendry et al. (1993), and has been revisited more recently by a number of authors (Jackson & Jeffries 2010; Schlaufman 2010), including a series of applications (e.g., Hirano et al. 2012, 2014; Walkowicz & Basri 2013; Morton & Winn 2014) to transiting systems observed with the NASA Kepler mission (Borucki et al. 2010; Koch et al. 2010).

Finally, asteroseismology provides an independent means of directly determining i_s. The asteroseismic estimation of i_s rests on our ability to resolve and extract signatures of rotation in the power spectra of non-radial modes of oscillation (e.g., Gizon & Solanki 2003; Ballot et al. 2006, 2008; Benomar et al. 2009; Campante et al. 2011). The asteroseismic technique requires bright targets and long-duration time series to attain the desired S/N and frequency resolution. The applicability of this technique depends entirely on the stellar properties and not on the planetary or orbital parameters, which is an advantage when measuring obliquities of systems with small and/or long-period planets. Following its application to host stars with single, non-transiting large planets discovered using the RV method (Wright et al. 2011; Gizon et al. 2013), the asteroseismic technique has been applied to several Kepler Sun-like hosts (Chaplin et al. 2013; Benomar et al. 2014; Lund et al. 2014a; Van Eylen et al. 2014). In addition, Huber et al. (2013a) used asteroseismology to measure a large obliquity for Kepler-56, a red giant hosting two transiting coplanar planets, thus showing that spin–orbit misalignments are not confined to hot-Jupiter systems. Another instance of an asteroseismic obliquity measurement of an evolved host is that of Kepler-432 (Quinn et al. 2015). Recently, the stellar inclination angles of the solar analogs 16 Cyg A and B (the B component hosts a Jovian planet) were determined using asteroseismology (Davies et al. 2015).

Here we present the first analysis of an ensemble of asteroseismic obliquity measurements obtained for solar-type stars with transiting planets. The asteroseismic sample consists of 25 Kepler planet-candidate host stars (also designated as Kepler Objects of Interest or KOIs), of which 20 are confirmed hosts. The host stars are distributed along the main sequence and subgiant branch, and all exhibit solar-like oscillations. The rest of the paper is organized as follows. In Section 2 we present a recap of the spin–orbit geometry. A detailed asteroseismic analysis of the individual planetary-system hosts follows in Section 3. In Section 4 we place statistical constraints on the spin–orbit alignment. Finally, a discussion of the results and conclusions make up Section 5.

Figure 1 shows an observer-oriented coordinate system (left panel) and an orbit-oriented coordinate system (right panel). In the former, the unit vector for orbital angular momentum, ${{\boldsymbol{n}}}_{{\rm{o}}}$, lies on the yz plane and is solely determined by the angle i_o between the planetary orbital axis and the line of sight (the orbital inclination angle). The angle i_o can be measured for a transiting planet via transit photometry (e.g., Charbonneau et al. 2000; Henry et al. 2000), in which case one necessarily has $\mathrm{sin}{i}_{{\rm{o}}}\approx 1$ (i.e., an edge-on orbit). Determination of the unit vector for stellar rotation angular momentum, ${{\boldsymbol{n}}}_{{\rm{s}}}$, requires a polar and an azimuthal angle, respectively i_s and λ. To avoid degeneracies, we restrict i_o and i_s to the interval $[0,\pi /2]$, while λ is allowed to vary in the interval $[-\pi ,\pi ]$. In the orbit-oriented coordinate system, ${{\boldsymbol{n}}}_{{\rm{s}}}$ is determined by the polar and azimuthal angles ψ and ϕ. The spin–orbit angle ψ is the angle between ${{\boldsymbol{n}}}_{{\rm{o}}}$ and ${{\boldsymbol{n}}}_{{\rm{s}}}$, taking values in the interval $[0,\pi ]$. Values of ψ lower (greater) than $\pi /2$ correspond to prograde (retrograde) orbits. The specific cases of a parallel, an antiparallel, and a polar orbit then correspond to ψ = 0, ψ = π, and ψ = π/2, respectively. The azimuthal angle ϕ is allowed to vary in the range $[-\pi ,\pi ]$ and takes the value π along the line of sight.

Zoom In Zoom Out Reset image size

The various angles are related according to the following equations (for a derivation see, e.g., Fabrycky & Winn 2009):

$$\begin{eqnarray}&&\mathrm{sin}{i}_{{\rm{s}}}\mathrm{sin}\lambda \;=\;\mathrm{sin}\psi \mathrm{sin}\phi ,\end{eqnarray} \tag{ 1a }$$

$$\begin{eqnarray}&&\mathrm{cos}\psi \;=\;\mathrm{sin}{i}_{{\rm{s}}}\mathrm{cos}\lambda \mathrm{sin}{i}_{{\rm{o}}}+\mathrm{cos}{i}_{{\rm{s}}}\mathrm{cos}{i}_{{\rm{o}}},\end{eqnarray} \tag{ 1b }$$

$$\begin{eqnarray}&&\mathrm{sin}\psi \mathrm{cos}\phi \;=\;\mathrm{sin}{i}_{{\rm{s}}}\mathrm{cos}\lambda \mathrm{cos}{i}_{{\rm{o}}}-\mathrm{cos}{i}_{{\rm{s}}}\mathrm{sin}{i}_{{\rm{o}}}.\end{eqnarray} \tag{ 1c }$$

Equation (1b) is of particular interest, as it allows determination of the spin–orbit angle ψ provided that measurements of i_o, i_s, and λ are available. A joint analysis of asteroseismology, the transit light curve, and the RM effect made it possible to determine ψ for the hot-Jupiter system HAT-P-7 (Kepler-2; Benomar et al. 2014; Lund et al. 2014a) and the multi-transiting system Kepler-25 (Benomar et al. 2014). Both these systems are revisited in this work.

For an individual transiting system, we would still expect to place mild constraints on ψ even when lacking a measurement of λ. In Figure 2 we show the posterior probability distribution (PPD; after normalization) for ψ conditioned on i_s and i_o, $p(\psi | {i}_{{\rm{s}}},{i}_{{\rm{o}}})$ (see Appendix A for a derivation of the analytical expression). We have assumed an edge-on orbit (i.e., ${i}_{{\rm{o}}}\;=\;90^\circ$), having selected three error-free values for the stellar inclination angle (${i}_{{\rm{s}}}\;=\;30^\circ$, ${i}_{{\rm{s}}}\;=\;60^\circ$, and ${i}_{{\rm{s}}}\;=\;85^\circ$). The main conclusions to be drawn from this simple exercise follow immediately from an inspection of Figure 2: for a transiting system, a small value of i_s implies a spin–orbit misalignment. The converse is not true, since a large value of i_s is consistent with, but does not necessarily imply, a spin–orbit alignment. The interpretation of the spin–orbit alignment in terms of the measured i_s can thus be ambiguous. This provides one of the main motivations for this work, namely that, in order to draw general inferences about spin–orbit alignment, a statistical analysis of an ensemble of such measurements is needed.

Zoom In Zoom Out Reset image size

3.1. Sample Characterization

Our asteroseismic sample consists of 25 solar-type KOIs, of which 20 are confirmed hosts and thus have been assigned a KeplerID (see Table 1 for a complete list). At the time of writing, all planets in the systems awaiting confirmation have been designated as "CANDIDATE" in the cumulative table of the NASA Exoplanet Archive¹⁶ (Akeson et al. 2013). Fundamental properties (e.g., surface gravity, radius, mass, and age) from a detailed asteroseismic analysis are available for all the KOIs in the sample (Silva Aguirre et al. 2015).

Table 1.  Asteroseismic Sample of Kepler Objects of Interest

KIC	KOI	Kepler ID	mKep	Quarters	No. of	Rp	Po	$v\mathrm{sin}{i}_{{\rm{s}}}$	Prot	References
					Planets	(R⊕)	(days)	(km s−1)	(days)
3425851	268	...	10.56	Q6.1–Q8.3	1	3.03	110.4	9.5 ± 0.5	7.873 ± 0.001	1, 9
3544595	69	Kepler-93	9.93	Q2.3–Q17.2	2	1.48	>3650	2.0 ± 0.5	...	1
3632418	975	Kepler-21	8.22	Q5.1–Q17.2	1	1.64	2.8	7.75 ± 1.0	12.591 ± 0.036	2, 9
4141376	280	...	11.07	Q6.1–Q17.2	1	1.94	11.9	3.5 ± 0.5	15.78 ± 2.12	1, 9
4349452	244	Kepler-25	10.73	Q5.1–Q17.2	3	2.60	123	8.2 ± 0.2	23.147 ± 0.039	3, 9
4914423	108	Kepler-103	12.29	Q3.1–Q12.3	2	3.37	179.6	3.8 ± 0.7	...	1
5866724	85	Kepler-65	11.02	Q3.1–Q17.2	3	1.42	8.1	10.4 ± 0.5	7.911 ± 0.155	1, 10
6278762a	3158	Kepler-444	8.72	Q15.1–Q17.2	5	0.40	9.7	∼2.2	...	4
6521045	41	Kepler-100	11.20	Q3.1–Q17.2	3	1.32	35.3	2.9 ± 0.5	24.988 ± 2.192	1, 9
7670943	269	...	10.93	Q6.1–Q17.2	1	1.75	18.0	13.5 ± 0.6	5.274 ± 0.033	1, 9
8077137	274	Kepler-128	11.39	Q6.1–Q17.2	2	1.13	22.8	7.7 ± 0.5	...	1
8292840	260	Kepler-126	10.50	Q5.1–Q17.2	3	1.52	100.3	10.4 ± 0.5	...	1
8478994	245	Kepler-37	9.71	Q5.1–Q17.2	4	0.30	51.2	1.1 ± 1.1	28.79 ± 3.29	5, 11
8494142	370	Kepler-145	11.93	Q7.1–Q17.2	2	2.65	42.9	8.9 ± 0.5	...	1
8866102b	42	Kepler-410 A	9.36	Q3.1–Q17.2	>1	2.84	17.8	15.0 ± 0.5	20.850 ± 0.007	1, 9
9414417	974	...	9.58	Q6.1–Q17.2	1	2.49	53.5	9.3 ± 0.5	10.847 ± 0.002	1, 9
9592705	288	...	11.02	Q6.1–Q17.2	1	3.17	10.3	9.4 ± 0.5	13.380 ± 0.099	1, 9
9955598	1925	Kepler-409	9.44	Q5.1–Q17.2	1	1.19	69.0	1.6 ± 0.5	...	1
10586004	275	Kepler-129	11.70	Q6.1–Q7.3	2	2.37	82.2	3.4 ± 0.5	...	1
10666592	2	Kepler-2	10.46	Q0–Q17.2	1	15.04	2.2	3.8 ± 0.5	...	6
10963065	1612	Kepler-408	8.77	Q2.3–Q17.2	1	0.82	2.5	3.4 ± 0.5	12.444 ± 0.172	1, 9
11295426	246	Kepler-68	10.00	Q5.1–Q17.2	3	0.95	580	0.5 ± 0.5	...	7
11401755	277	Kepler-36	11.87	Q6.1–Q17.2	2	1.49	16.2	4.9 ± 1.0	...	8
11807274	262	Kepler-50	10.42	Q6.1–Q17.2	2	1.71	9.4	11.7 ± 0.6	7.553 ± 0.755	1, 9
11904151	72	Kepler-10	10.96	Q2.2–Q17.2	2	1.47	45.3	1.9 ± 0.5	...	1

Notes. The number of planets refers to the total number of confirmed or else candidate transiting planets. R_p is the planetary radius (source: NASA Exoplanet Archive); for multiple-planet systems, the smallest planet in the system is considered. P_o is the orbital period (source: NASA Exoplanet Archive); for multiple-planet systems, the longest-period planet in the system is considered. See references below for sources of $v\mathrm{sin}{i}_{{\rm{s}}}$ and P_rot.

^aA recent spectroscopic follow-up of this star with the Subaru High Dispersion Spectrograph (HDS) returned an upper limit for $v\mathrm{sin}{i}_{{\rm{s}}}$ (at the 2σ level) of 0.56 km s⁻¹ (T. Hirano 2016, private communication). ^bThe P_o value is for the transiting planet (the TTV signal has a period of 957 days). The P_rot value is associated with a blended star (Van Eylen et al. 2014).

References. (1) Huber et al. (2013b), (2) Howell et al. (2012), (3) Albrecht et al. (2013), (4) Campante et al. (2015), (5) Barclay et al. (2013), (6) Pál et al. (2008), (7) Gilliland et al. (2013), (8) Carter et al. (2012), (9) McQuillan et al. (2013), (10) Hirano et al. (2012, 2014), (11) Walkowicz & Basri (2013).

Download table as:  ASCIITypeset image

A systematic study of Kepler planet-candidate hosts using asteroseismology was presented by Huber et al. (2013b), in which fundamental properties were determined for 66 host stars based on their average asteroseismic parameters. This raised the number of Kepler hosts with asteroseismic solutions to nearly 80 stars. Whether or not a given host star is included in the present sample was determined by our ability to resolve and extract signatures of rotation in the oscillation spectrum, which required relatively bright targets (Kepler-band magnitude ${m}_{{\rm{Kep}}}\lesssim 12$) and multi-quarter observations. The intrinsic stellar properties have also played a crucial role in this regard, since it is well known that the signatures of rotation tend to be substantially blended in the power spectra of main-sequence solar-like oscillators hotter than about $6400\ {\rm{K}}$ (e.g., Appourchaux et al. 2012). Figure 3 displays the sample on a $\mathrm{log}g$ versus ${T}_{{\rm{eff}}}$ diagram. They are predominantly positioned along the main sequence and range in spectral type from early K to late F (i.e., $5000\ {\rm{K}}\lesssim {T}_{{\rm{eff}}}\lesssim 6500\ {\rm{K}}$). A number of stars in the sample seem to have evolved slightly beyond the main-sequence phase, one example being the bright subgiant Kepler-21 (Howell et al. 2012). There is also a varying level of evidence of mixed (e.g., Osaki 1975; Aizenman et al. 1977) quadrupole modes in the power spectra of Kepler-36, Kepler-100, Kepler-128, and Kepler-129, an indication that these stars may have already left the main sequence. The fact that central hydrogen has been depleted in models of these stars (Silva Aguirre et al. 2015) supports this scenario.

Zoom In Zoom Out Reset image size

The sample contains 16 multiple-planet systems, of which all except Kepler-93 (Ballard et al. 2014) and Kepler-410 A (Van Eylen et al. 2014) are also multi-transiting systems.¹⁷ Moreover, a non-transiting planet was revealed by RV measurements orbiting beyond the two transiting planets in the Kepler-25 (Marcy et al. 2014) and Kepler-68 (Gilliland et al. 2013) systems. Most of the multi-transiting systems in our sample have had the eccentricities of their planets measured from transit photometry (Van Eylen & Albrecht 2015), which revealed a tendency toward low eccentricities that are consistent with nearly circular orbits. Of the remaining nine single-planet systems, only one (HAT-P-7; Pál et al. 2008) is a hot-Jupiter system. From Table 1, we also see a clear prevalence of systems that contain small planets (i.e., ${R}_{{\rm{p}}}\leqslant 4\;{R}_{\oplus }$) and long-period planets (i.e., ${P}_{{\rm{o}}}\gt 10\ {\rm{days}}$).

3.2. Data Preparation

3.3. Estimation of the Stellar Inclination Angle

3.3.1. Principle

Solar-like oscillations are predominantly acoustic global standing waves. Commonly called p modes, owing to the fact that pressure plays the role of the restoring force, these modes are intrinsically damped and stochastically excited by near-surface convection (for a review see, e.g., Christensen-Dalsgaard 2004; Cunha et al. 2007; Chaplin & Miglio 2013). The oscillation modes are characterized by the radial order n, the spherical degree l, and the azimuthal order m. Radial modes have l = 0, whereas non-radial modes have l > 0. Values of m range from $-l$ to l, meaning that there are $2\;l+1$ azimuthal components for a given multiplet of degree l. Observed oscillation modes are typically high-order modes of low spherical degree, with the associated power spectrum showing a pattern of peaks with near-regular frequency separations (Vandakurov 1967; Tassoul 1980).

The asteroseismic estimation of i_s rests on our ability to resolve and extract signatures of rotation in the power spectra of non-radial modes of oscillation. Rotation lifts the degeneracy in the frequencies, ν_nl, of non-radial modes and introduces a dependence of the oscillation frequencies on m, with prograde (retrograde) modes (with m > 0 and m < 0, respectively) having frequencies slightly higher (lower) than the axisymmetric mode (m = 0) in the observer's frame of reference. For the fairly modest values of the stellar angular velocity Ω that are typical of solar-like oscillators, the effect of rotation can be treated following a perturbative analysis (e.g., Reese et al. 2006) and the star is generally assumed to rotate as a solid body (i.e., Ω = const.). In the limit of solid-body rotation, the frequency ν_nlm of a mode, as observed in an inertial frame, can be expressed to first order as (Ledoux 1951)

$$\begin{eqnarray}&&{\nu }_{{nlm}}\;=\;{\nu }_{{nl}0}+m\displaystyle \frac{{\rm{\Omega }}}{2\pi }(1-{C}_{{nl}})\approx {\nu }_{{nl}0}+m{\nu }_{{\rm{s}}}.\end{eqnarray} \tag{ 2 }$$

The kinematic splitting $m{\rm{\Omega }}/(2\pi )$ is corrected for the effect of the Coriolis force through the dimensionless Ledoux constant, C_nl. For high-order p modes, ${C}_{{nl}}\ll 1$, with the rotational splitting being dominated by advection and given approximately by the angular velocity, i.e., ${\nu }_{{\rm{s}}}\approx {\rm{\Omega }}/(2\pi )$ (e.g., Lund et al. 2014b; Davies et al. 2015).

To a second order of approximation, centrifugal effects that disrupt the equilibrium structure of the star can be taken into account through an additional frequency perturbation (e.g., Ballot 2010). This perturbation scales as the ratio of the centrifugal to the gravitational forces at the stellar surface, i.e., ${{\rm{\Omega }}}_{{\rm{surf}}}^{2}{R}_{{\rm{s}}}^{3}/({{GM}}_{{\rm{s}}})$, where G is the gravitational constant. We made use of the available values of P_rot in Table 1 to compute the ratios of the surface angular velocity to the Keplerian break-up rotation rate, i.e., ${{\rm{\Omega }}}_{{\rm{surf}}}/{{\rm{\Omega }}}_{{\rm{K}}}\equiv 2\pi /({P}_{{\rm{rot}}}\sqrt{{{GM}}_{{\rm{s}}}/{R}_{{\rm{s}}}^{3}})$. We obtained ${({{\rm{\Omega }}}_{{\rm{surf}}}/{{\rm{\Omega }}}_{{\rm{K}}})}^{2}\lesssim 9\times {10}^{-4}$ for the stars in the asteroseismic sample and have thus decided to neglect second-order effects in the present analysis.

Assuming energy equipartition between multiplet components with different azimuthal order,¹⁹ the dependence of mode power on m is given by (Dziembowski 1977; Dziembowski & Goode 1985; Gizon & Solanki 2003)

$$\begin{eqnarray}&&{{\mathcal{E}}}_{{lm}}({i}_{{\rm{s}}})=\displaystyle \frac{(l-| m| )!}{(l+| m| )!}{[{P}_{l}^{| m| }(\mathrm{cos}{i}_{{\rm{s}}})]}^{2},\end{eqnarray} \tag{ 3 }$$

where ${P}_{l}^{m}(x)$ are the associated Legendre functions and the sum over m of ${{\mathcal{E}}}_{{lm}}({i}_{{\rm{s}}})$ has been normalized to unity. Measurement of the relative power of the azimuthal components in a non-radial multiplet thus provides a direct estimate of the stellar inclination angle. The above formalism further relies on the assumption that contributions to the observed intensity across the visible stellar disk depend only on the angular distance from the disk center, which is valid for photometric observations. According to Equation (3), when the stellar spin axis points toward the observer (pole-on configuration), only the axisymmetric mode is visible and no inference can thus be made about rotation. When the spin axis lies on the plane of the sky (edge-on configuration), as is approximately the case for the Sun, observations are essentially sensitive only to modes with even $| l-m|$. Given sufficient frequency resolution and S/N, it will be the intrinsic ratio ${\nu }_{{\rm{s}}}/{\rm{\Gamma }}$ (where Γ is the full width at half maximum, or linewidth, of the mode profile of each azimuthal component; see Equation (12)) that determines whether it is possible to resolve the azimuthal components (Ballot et al. 2006; Bedding et al. 2010). Moreover, dipole (l = 1) modes are approximately three times more prominent in the power spectra of intensity observations than quadrupole (l = 2) modes of similar frequency, and consequently it is the former modes that ultimately determine our ability to constrain i_s.

3.3.2. Results

A detailed fitting of the modes of oscillation was conducted to extract signatures of rotation from the power spectra of stars in the asteroseismic sample (see Appendix B for a description of the method; in Appendix C we test the robustness of the returned uncertainties on i_s). This allowed us to map the joint PPD of i_s and ν_s. An example is shown in Figure 4 (bottom left panel) for Kepler-410 A, together with the corresponding marginalized PPDs (top left and bottom right panels). Also shown in the bottom left panel is an estimate of the projected splitting (black solid line), given by²⁰ $\;{\nu }_{{\rm{s}}}\mathrm{sin}{i}_{{\rm{s}}}\approx v\mathrm{sin}{i}_{{\rm{s}}}/(2\pi {R}_{{\rm{s}}})$. This assumes that the internal rotation rates probed by the asteroseismic splitting are similar to the surface rate of rotation. For a typical 1.2 M_⊙ star in our sample (such as Kepler-410 A), however, the asteroseismic splitting is nearly equally sensitive to rotation in the radiative and convective zones (Lund et al. 2014b; Benomar et al. 2015), and the assumption that the surface and asteroseismic rotation periods are similar should thus be regarded as approximate. We also note that different authors very often disagree on the value of the measured $v\mathrm{sin}{i}_{{\rm{s}}}$. The bottom left panel of Figure 4 displays such an example, where besides the projected splitting based on the $v\mathrm{sin}{i}_{{\rm{s}}}$ value from Huber et al. (2013b) given in Table 1 (black solid line), we also use the other available $v\mathrm{sin}{i}_{{\rm{s}}}$ value for this star from Molenda-Żakowicz et al. (2013; gray solid line). Accurate and precise $v\mathrm{sin}{i}_{{\rm{s}}}$ measurements from spectroscopy are difficult to obtain for solar-type stars, since the rotational line broadening is often comparable to the effects of instrumental broadening and macroturbulence (${v}_{{\rm{mac}}}$). In a recent work, Doyle et al. (2014) used accurate $v\mathrm{sin}{i}_{{\rm{s}}}$ values from asteroseismology to break the degeneracy between the spectroscopic $v\mathrm{sin}{i}_{{\rm{s}}}$ and ${v}_{{\rm{mac}}}$ in spectral line profiles, thus obtaining a calibration for stellar macroturbulence. The marginalized PPDs are presented as histograms, with dotted and dotted–dashed lines respectively enclosing the 68.3% and 95.4% highest posterior density (HPD) credible regions. Our results for the stellar inclination angle of Kepler-410 A are in very good agreement with those obtained by Van Eylen et al. (2014). The stellar inclination is tightly constrained and its posterior distribution indicates an edge-on configuration. This was to be expected from inspection of the l = 1 mode profiles in Figure 16.

Zoom In Zoom Out Reset image size

Next, we present the results for the two targets in the asteroseismic sample that have an available λ measurement in the literature, namely, HAT-P-7 and Kepler-25 (see Figures 5–8; see also Table 2). Of particular interest are Figures 6 and 8, showing the derived PPDs of the spin–orbit angle ψ (black histograms). These were sampled by means of Monte Carlo simulations (via Equation (1b)) using the PPD of i_s from our analysis, and assuming normal and uniform distributions, respectively, for λ and i_o around their adopted literature values (see Table 2). Also shown are the PPDs obtained by instead assuming an isotropic distribution in λ (gray histograms). We find that the orbit of the hot Jupiter HAT-P-7b is likely to be retrograde, while that of Kepler-25c is well aligned with the stellar spin axis. Our findings for Kepler-25 are not in agreement with the statement made by Benomar et al. (2014) that this system is mildly misaligned and this point will be discussed further in Section 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Results on the Spin–orbit Angle ψ for HAT-P-7 and Kepler-25

			ψ
Kepler ID	λ	io	Median	68.3% HPD Lower Bound	68.3% HPD Upper Bound	References
	(deg)	(deg)	(deg)	(deg)	(deg)
HAT-P-7b	182.5 ± 9.4a/$-{132.6}_{-16.3}^{+10.5}$/155 ± 37	${80.8}_{-1.2}^{+2.8}$/83.111 ± 0.030	116.4	93.5	138.4	1, 2, 3, 4
Kepler-25c	−0.5 ± 5.7	87.27 ± 0.05	12.6	1.6	19.3	5

Notes. When several measurements of λ or i_o are available (see references below), the adopted value is given in boldface. Values for λ all come from measurements of the RM effect.

^aEquivalent to $\lambda \;=\;-177\buildrel{\circ}\over{.} 5\pm 9\buildrel{\circ}\over{.} 4$ if brought into the interval $[-\pi ,\pi ]$.

References. (1) Winn et al. (2009), (2) Narita et al. (2009), (3) Albrecht et al. (2012), (4) Morris et al. (2013), (5) Albrecht et al. (2013).

Download table as:  ASCIITypeset image

Asteroseismic results for the remainder of the stars in the sample are shown in Figures 18–39. Even though i_s is independent of the rotational splitting ν_s, their measured values are highly correlated, as can be seen from the joint PPDs. Consequently, even when a constrained solution for i_s and ν_s cannot be obtained, as is often the case in the present analysis, it is still possible to estimate the projected splitting ${\nu }_{{\rm{s}}}\mathrm{sin}{i}_{{\rm{s}}}$. Table 3 reports the 68.3% and 95.4% HPD credible regions for the stellar inclination angle i_s for all the stars in the asteroseismic sample. Note that the PPD of i_s is often too asymmetric to be adequately summarized by a single estimate and we have thus refrained from tabulating any measures of central tendency, such as the median or the mode of the distributions. We stress that Bayesian inference does not provide point estimates of parameters. Instead the Bayesian solution to the problem of parameter estimation is the full PPD.²¹ The several statistical analyses conducted in the next section are based on the full posteriors. Our results for i_s suggest that there are five other systems in the asteroseismic sample besides HAT-P-7 that could potentially be misaligned, namely, Kepler-50, Kepler-93, Kepler-145, Kepler-409, and KOI-280. And while this interpretation is certainly valid when considering the 68.3% HPD credible regions, we note that our results are consistent with alignment at the 2σ level for all systems in the sample.

Table 3.  Results on the Stellar Inclination Angle i_s for the Stars in the Asteroseismic Sample

						is
KOI	Kepler ID	No. of	Frequency Range	S/N	νs/Γ	68.3% HPD Credible Region	95.4% HPD Credible Region	MW14 Upper Bound
		Orders	(μHz)			(deg)	(deg)	(deg)
268	...	6	1615–2170	0.30	0.21	[38.5, 85.0]	[18.6, 90.0]	...
69	Kepler-93	6	2735–3610	2.7	0.49	[48.6, 69.8]	[45.1, 87.4]	...
975	Kepler-21	8	810–1300	12.4	0.31	[70.0, 90.0]	[53.4, 90.0]	89.4
280	...	5	2520–3160	1.3	0.73	[41.9, 74.8]	[38.4, 90.0]	88.8
244	Kepler-25	6	1715–2300	1.0	0.46	[68.7, 90.0]	[54.3, 90.0]	...
108	Kepler-103	5	1355–1765	1.0	0.28	[41.7, 90.0]	[17.6,90.0]	...
85	Kepler-65	6	1475–2010	1.4	0.56	[74.3, 90.0]	[60.8, 90.0]	90.0
3158	Kepler-444	6	3730–4800	2.2	0.31	[31.3, 90.0]	[22.7, 90.0]	...
41	Kepler-100	7	1120–1660	4.3	0.32	[65.0, 90.0]	[46.2, 90.0]	89.4
269	...	6	1540–2075	1.0	0.54	[58.5, 90.0]	[39.1, 90.0]	63.6
274	Kepler-128	6	1065–1480	1.4	0.27	[56.9, 90.0]	[36.5, 90.0]	...
260	Kepler-126	6	1615–2175	1.9	0.44	[71.8, 90.0]	[59.5, 90.0]	...
245	Kepler-37	6	3710–4780	1.0	0.35	[52.0, 90.0]	[27.1, 90.0]	...
370	Kepler-145	5	960–1270	1.3	0.36	[17.6, 65.9]	[14.1, 90.0]	...
42	Kepler-410 A	8	1550–2305	3.3	0.55	[82.8, 90.0]	[76.8, 90.0]	...
974	...	8	800–1280	5.2	0.27	[55.4, 90.0]	[34.1, 90.0]	68.1
288	...	7	770–1150	1.8	0.27	[51.2, 90.0]	[25.9, 90.0]	90.0
1925	Kepler-409	7	2880–3950	5.1	0.42	[39.7, 78.4]	[33.2, 90.0]	...
275	Kepler-129	5	1150–1505	2.1	0.40	[32.4, 86.6]	[15.5, 90.0]	...
2	Kepler-2	7	840–1255	2.2	0.17	[19.1, 69.8]	[15.8, 90.0]	...
1612	Kepler-408	8	1705–2530	10.1	0.28	[53.5, 90.0]	[36.5, 90.0]	...
246	Kepler-68	8	1690–2500	10.0	0.21	[36.9, 90.0]	[18.7, 90.0]	...
277	Kepler-36	6	980–1390	1.8	0.38	[53.6, 90.0]	[24.8, 90.0]	...
262	Kepler-50	7	1170–1700	2.7	0.53	[50.7, 73.9]	[48.9, 90.0]	89.4
72	Kepler-10	5	2200–2790	1.9	0.37	[39.7, 90.0]	[12.7, 90.0]	...

Note. The number of orders refers to the number of observed orders entering the fit (see Appendix B). Frequency ranges are approximate. The height-to-background ratio of the dipole mode in the vicinity of the frequency of maximum oscillation amplitude (${\nu }_{{\rm{max}}}$) is taken as a measure of the S/N in the p modes. The ratio ${\nu }_{{\rm{s}}}/{\rm{\Gamma }}$ refers to the order around ${\nu }_{{\rm{max}}}$ (for reference, ${\rm{\Gamma }}\approx 1.0\ \mu \mathrm{Hz}$ and ${\nu }_{{\rm{s}}}/{\rm{\Gamma }}\approx 0.4$ for the Sun). The rightmost column gives the 95% confidence upper bound on i_s obtained by Morton & Winn (2014).

Download table as:  ASCIITypeset image

4.1. Asteroseismic Sample

Figure 9 shows the concatenated PPDs of $\mathrm{cos}{i}_{{\rm{s}}}$, made by concatenating the individual posteriors and normalizing. This is shown for all stars in the asteroseismic sample (25 systems), as well as for the subsamples of single- (11) and multi-transiting (14) systems. The isotropic distribution is represented by a horizontal dashed line. We may have naively expected the concatenated posterior(s) not to exhibit the amount of structure seen in Figure 9. This is a consequence of the small size of the asteroseismic sample, with a few individual posteriors eventually dominating in specific i_s ranges.

Zoom In Zoom Out Reset image size

The stars in the asteroseismic sample appear to preferentially display large values of i_s (or, equivalently, small values of $\mathrm{cos}{i}_{{\rm{s}}}$), this being more accentuated than if the i_s had been drawn from an isotropic distribution. This departure from isotropy suggests that the directions of the stellar spin axis (${{\boldsymbol{n}}}_{{\rm{s}}}$) and planetary orbital axis (${{\boldsymbol{n}}}_{{\rm{o}}}$) are correlated. Given samplings of the PPDs of $\mathrm{cos}{i}_{{\rm{s}}}$ for the N = 25 stars in the asteroseismic sample, ${\{\mathrm{cos}{i}_{{\rm{s}}}\}}_{n=1}^{N}$, we would like to infer the distribution of the spin–orbit angle ψ. As a first step, we need to devise a parameterized model for the distribution function of ψ, ${f}_{{\boldsymbol{\alpha }}}(\psi )$, as well as a prior on the model parameters ${\boldsymbol{\alpha }}$. We could then, in principle, use the observed data to constrain these parameters by means of a hierarchical Bayesian analysis (see Appendix E for details). We note that the concatenated posteriors in Figure 9 are shown here only for comparison to Figure 9 of Hirano et al. (2014) and Figure 5 of Morton & Winn (2014). We do not use these concatenated posteriors directly for inference. Instead, we conduct a forward modeling of the observed data by defining a model distribution ${f}_{{\boldsymbol{\alpha }}}(\psi )$ and finding the model parameters α that best explain those data. Since we are modeling the distribution prior to observation, we are effectively performing a deconvolution that works on heteroscedastic data (i.e., with changing variance).

A Fisher distribution. An isotropic distribution for ψ is known to be inadequate (e.g., Winn et al. 2006), which also follows from our realization above that ${{\boldsymbol{n}}}_{{\rm{s}}}$ and ${{\boldsymbol{n}}}_{{\rm{o}}}$ are apparently correlated. We choose to model the distribution function of ψ as a Fisher distribution (Fisher 1953), the analog of a zero-mean Gaussian distribution on a sphere. The Fisher distribution has been previously proposed by Fabrycky & Winn (2009) to model ${f}_{{\boldsymbol{\alpha }}}(\psi )$. Its probability density function is given by

$$\begin{eqnarray}&&{f}_{{\boldsymbol{\alpha }}}(\psi )\equiv {p}_{{\rm{F}}}(\psi | \kappa )=\displaystyle \frac{\kappa }{2\mathrm{sinh}\kappa }\mathrm{exp}(\kappa \mathrm{cos}\psi )\mathrm{sin}\psi ,\end{eqnarray} \tag{ 4 }$$

where κ measures the probability concentration around ψ = 0. As $\kappa \to 0$, the distribution becomes the isotropic distribution, whereas for $\kappa \to \infty$ it becomes a Rayleigh distribution of width $\sigma \;=\;{\kappa }^{-1/2}$. However, our observable is $\mathrm{cos}{i}_{{\rm{s}}}$ and not ψ, and so we instead use the probability distribution for $\mathrm{cos}{i}_{{\rm{s}}}$ given κ as derived by Morton & Winn (2014, hereafter MW14):

$$\begin{eqnarray}{f}_{{\boldsymbol{\alpha }}}(\mathrm{cos}{i}_{{\rm{s}}}) & \equiv & {p}_{{\rm{F}}}(\mathrm{cos}{i}_{{\rm{s}}}| \kappa )=\displaystyle \frac{2\kappa }{\pi \mathrm{sinh}\kappa }{\displaystyle \int }_{\mathrm{cos}{i}_{{\rm{s}}}}^{1}\\ & & \times \displaystyle \frac{\mathrm{cosh}(\kappa \sqrt{1-{y}^{2}})}{\sqrt{1-{y}^{2}}}\displaystyle \frac{1}{\sqrt{1-{(\mathrm{cos}{i}_{{\rm{s}}}/y)}^{2}}}\;{dy}.\end{eqnarray} \tag{ 5 }$$

We refer to this parameterized model as the "single-Fisher model" and adopt a uniform prior for κ in the range [0, 200]. Figure 10 shows ${p}_{{\rm{F}}}(\psi | \kappa )$ and ${p}_{{\rm{F}}}(\mathrm{cos}{i}_{{\rm{s}}}| \kappa )$ for a set of concentration parameters κ. Finally, we sample the joint PPD of the model parameters ${\boldsymbol{\alpha }}$ (or simply κ in the present case) following a hierarchical Bayesian scheme.

Zoom In Zoom Out Reset image size

The PPD of the concentration parameter κ is shown in Figure 11. Distributions for the different cohorts of systems are color-coded. Table 4 presents a summary of the results. The posterior of κ is dominated by the (uniform) prior over most of the range [0, 200], the exception being at small κ, where the posterior falls sharply at $\kappa \approx 10$ (the logarithmic scale along the horizontal axis in Figure 11 emphasizes this point). Based on the lower bounds of the 95.4% HPD credible regions (see Table 4), it can be stated that the concentration parameter κ is significantly different from zero, thus ruling out isotropy.

Zoom In Zoom Out Reset image size

Table 4.  Results from the Hierarchical Bayesian Analysis

Model	Single-transiting			Multi-transiting			All
	κ	f	E	κ	f	E	κ	f	E
Asteroseismic Sample
Single-Fisher	${106}_{-94}^{+92}$	...	...	${101}_{-88}^{+96}$	...	...	${107}_{-89}^{+93}$	...	...
Combined Sample
Single-Fisher	${8.7}_{-5.4}^{+6.7}$	...	1.10 × 105	${118}_{-83}^{+82}$	...	5.11 × 107	${11.5}_{-5.7}^{+7.5}$	...	1.61 × 109
Mixture	${126}_{-86}^{+74}$	${0.21}_{-0.18}^{+0.24}$	1.24 × 107	${116}_{-82}^{+84}$	${0.08}_{-0.08}^{+0.24}$	6.16 × 106	${129}_{-80}^{+71}$	${0.13}_{-0.12}^{+0.16}$	7.97 × 1011

Note. We quote the median and 95.4% HPD credible region for both κ and f. The Bayesian model evidence, E, is reported for each cohort of systems in the combined sample.

Download table as:  ASCIITypeset image

The constraints on spin–orbit alignment made possible by an analysis of the asteroseismic sample are in qualitative agreement with the outcome of the analysis presented in Hirano et al. (2014), who have estimated i_s for a sample of 25 (coincidentally) KOIs based on measurements of their rotation period, rotational line broadening, and stellar radius. Furthermore, our results cannot be used to support the recent finding by MW14 that the obliquities of systems with a single transiting planet are systematically larger than those with multiple transiting planets. MW14 have, however, considered a substantially larger sample, having performed an ensemble analysis of a compilation of 70 KOIs whose obliquity measurements had previously been published (these included most of the KOIs analyzed by Hirano et al. 2014).

The small size of the asteroseismic sample prevents us from placing more stringent constraints on spin–orbit alignment and more systems would be needed to draw firmer conclusions. In the next section we explore the possibility of extending the MW14 sample by combining it with our asteroseismic sample, thus forming the largest ensemble to date of measured i_s for exoplanet-host stars. We will be particularly interested in assessing the impact of adding our asteroseismic sample (representing just over 1/4 of the combined sample) to the ensemble analysis.

4.2. Combined Sample

We combined the asteroseismic sample with the compilation by MW14 in order to place statistical constraints on the spin–orbit alignment of exoplanet systems. There are 62 additional (non-overlapping) systems in the sample of MW14, of which 23 are multi-transiting and 39 are single-transiting (we note that KOI-2636 has in the meantime been reclassified as multi-transiting). We thus obtain a combined sample of 87 systems, of which 37 are multi-transiting and 50 are single-transiting. There are at least three reasons for considering an analysis that uses the combined sample: (i) as already mentioned, the combined sample is, to date, the largest ensemble of measured i_s for exoplanet-host stars; (ii) the two samples being combined are complementary in terms of T_eff coverage (see Figure 12); (iii) the hierarchical Bayesian approach adopted in this work (compare Section 4.1) is able to handle the inhomogeneous nature of the combined sample in a straightforward manner.

Zoom In Zoom Out Reset image size

Figure 12 displays the stars in the combined sample in a diagram of $\mathrm{log}g$ versus T_eff. Multi- and single-transiting systems are seen to span similar ranges of effective temperature ($4500\;{\rm{K}}\lesssim {T}_{{\rm{eff}}}\lesssim 6500\ {\rm{K}}$) and surface gravity ($3.9\lesssim \mathrm{log}g\lesssim 4.6\ {\rm{dex}}$). Stars belonging to the asteroseismic sample are responsible for populating the hot (and hence massive) end of the parameter space, where flux modulations due to starspots are hardly detectable and a measurement of the photometric P_rot is thus made difficult (e.g., Radick et al. 1982). Conversely, the cool end of the parameter space is mainly populated by stars belonging to the MW14 sample, since their oscillation amplitudes are too small, and the stars in general too faint, to allow detection of the oscillations (e.g., Campante et al. 2014). Furthermore, planets in multi-transiting systems tend to be systematically smaller than planets in single-transiting systems (see Latham et al. 2011), a tendency made clearer in Figure 13. We also note that the fraction of multi-transiting systems is higher than in the full Kepler sample of planet-candidate host stars, where multi-transiting systems account for only 23% of the total (Burke et al. 2014). The fact that stars in the combined sample are relatively bright may partly explain the high observed fraction of multi-transiting systems, since the smaller planets in these systems would be preferentially detected around bright stars.

Zoom In Zoom Out Reset image size

Given samplings of the PPDs of $\mathrm{cos}{i}_{{\rm{s}}}$ for the N = 87 stars in the combined sample, we once more try to infer the distribution of the spin–orbit angle ψ following a hierarchical Bayesian scheme (compare Section 4.1). Having now access to a larger sample, we decided to also employ a slightly more complex model to represent the distribution function of ψ in addition to the single-Fisher model.

A sum of an isotropic and a Fisher distribution. An alternative to modeling the distribution function of ψ is to assume that all spin–orbit angles are drawn either from an isotropic distribution (with probability f) or from a Fisher distribution (with probability $1-f$). This so-called "mixture model" can be used to describe a scenario where two different channels exist by which planets migrate, ultimately giving rise to two different obliquity populations. One of these channels produces an isotropic distribution of spin–orbit angles and f can then be interpreted as the fraction of systems that follow such a migration route. We thus have

$$\begin{eqnarray}&&{f}_{{\boldsymbol{\alpha }}}(\psi )\equiv \displaystyle \frac{f}{2}\mathrm{sin}\psi +(1-f){p}_{{\rm{F}}}(\psi | \kappa ),\end{eqnarray} \tag{ 6 }$$

or, in terms of the observable $\mathrm{cos}{i}_{{\rm{s}}}$,

$$\begin{eqnarray}&&{f}_{{\boldsymbol{\alpha }}}(\mathrm{cos}{i}_{{\rm{s}}})\equiv f+(1-f){p}_{{\rm{F}}}(\mathrm{cos}{i}_{{\rm{s}}}| \kappa ).\end{eqnarray} \tag{ 7 }$$

The mixture model has two free parameters, namely, f and κ. We assign uniform priors to both.

The PPD of the concentration parameter κ in the single-Fisher model is shown in Figure 14. Figure 15 shows the PPDs of the fraction f of isotropic systems (top panel) and the concentration parameter κ (bottom panel) in the mixture model. Distributions for the different cohorts of systems are color-coded in each panel. Table 4 presents a summary of the results, where we also report the Bayesian model evidence, E, for each cohort of systems. The Bayesian evidence is computed by taking the integral of Equation (E1) over the parameter space spanned by the model parameters ${\boldsymbol{\alpha }}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

MW14 have also considered the single-Fisher model, and the addition of the asteroseismic sample in the present analysis leads to consistent results within the quoted credible regions. For multi-transiting systems, however, the posteriors of κ in both works differ in shape and appear to be dominated at large κ by their respective priors. Our choice of a uniform prior on κ, as opposed to $p(\kappa )\propto {(1+{\kappa }^{2})}^{-3/4}$ in MW14, was made to prevent the underestimation of κ when κ is large. We confirmed our suspicion that the prior may be dominating at large κ by repeating the analysis having used the latter prior on κ instead. We take this as an indication that the available data may not be able to fully constrain the posterior of κ for multi-transiting systems, in particular at large κ.

Figures 14 and 15 possess a number of other interesting features: (i) the posteriors of κ (in both models) and f are not compatible with an isotropic distribution for ψ; (ii) multi-transiting systems tend to be characterized by a large κ or, equivalently, by low obliquities; (iii) the posterior of κ is similar for multi-transiting systems irrespective of the model being considered, but a difference exists in the case of single-transiting systems. These points will be used in support of the discussion presented in the next section.

The main outputs of this work are twofold, namely, the presentation of individual asteroseismic results on the stellar inclination angle for a sample of Kepler host stars (Section 3) and the placement of statistical constraints on the spin–orbit alignment of exoplanet systems (Section 4).

5.1. Asteroseismic Analysis

5.2. Statistical Analysis

5.3. Outlook

This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission Directorate. Some of the data presented in this paper were obtained from the Kepler Asteroseismic Science Operations Center (KASOC) database. T.L.C., J.S.K., G.R.D., W.J.C., D.B., Y.P.E., and T.S.H.N. acknowledge the support of the UK Science and Technology Facilities Council (STFC). Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (Grant DNRF106). The research is supported by the ASTERISK project (ASTERoseismic Investigations with SONG and Kepler) funded by the European Research Council (Grant agreement no.: 267864). The research leading to the presented results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 338251 (StellarAges). Work by J.N.W. was supported by the NASA Origins program (NNX11AG85G). S.B. acknowledges support from NASA grant NNX13AE70G. C.K. and V.S.A. acknowledge the support of the Villum Foundation. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This study has made use of René Heller's Holt–Rossiter–McLaughlin Encyclopaedia. This research made use of the Python package obliquity developed by Timothy Morton and made available at https://github.com/timothydmorton/obliquity.

Facility: Kepler - The Kepler Mission.

Suppose an observer has measured i_s and i_o for a particular transiting system. What can then be inferred about ψ? Here we derive an analytical expression for $p(\psi | {i}_{{\rm{s}}},{i}_{{\rm{o}}})$, the PPD for the spin–orbit angle ψ conditioned on the stellar inclination angle i_s and on the orbital inclination angle i_o. As far as the authors are aware, such a derivation has not been presented in the literature.

We start by deriving an expression for $p({i}_{{\rm{s}}}| \psi ,{i}_{{\rm{o}}})$, i.e., the PPD for ${i}_{{\rm{s}}}$ conditioned on ψ and i_o. We assume that for a given ψ, the probability distribution of the azimuthal angle ϕ is uniformly distributed between $-\pi$ and $+\pi$ (see Fabrycky & Winn 2009). We then express i_s in terms of ψ, i_o, and ϕ by eliminating λ from Equations (1a)–(1c):

$$\begin{eqnarray}&&{i}_{{\rm{s}}}(\psi ,{i}_{{\rm{o}}},\phi )=\mathrm{arccos}(\mathrm{cos}\psi \mathrm{cos}{i}_{{\rm{o}}}-\mathrm{sin}\psi \mathrm{sin}{i}_{{\rm{o}}}\mathrm{cos}\phi ).\end{eqnarray} \tag{ 8 }$$

Since ${i}_{{\rm{s}}}(\psi ,{i}_{{\rm{o}}},-\phi )={i}_{{\rm{s}}}(\psi ,{i}_{{\rm{o}}},\phi )$, we need only consider ϕ in the interval $[0,\pi ]$. Next, we transform variables (e.g., Meyer 1970, chap. 5) from ϕ to i_s in order to arrive at an expression for $p({i}_{{\rm{s}}}| \psi ,{i}_{{\rm{o}}})$:

$$\begin{eqnarray}p({i}_{{\rm{s}}}| \psi ,{i}_{{\rm{o}}}) & = & p(\phi | \psi ,{i}_{{\rm{o}}})\left|\displaystyle \frac{\partial \phi }{\partial {i}_{{\rm{s}}}}\right|=\displaystyle \frac{1}{\pi }\left|\displaystyle \frac{\partial \phi }{\partial {i}_{{\rm{s}}}}\right|\\ & = & \displaystyle \frac{1}{\pi }\displaystyle \frac{\mathrm{sin}{i}_{{\rm{s}}}}{\sqrt{{\mathrm{sin}}^{2}\psi {\mathrm{sin}}^{2}{i}_{{\rm{o}}}-{(\mathrm{cos}\psi \mathrm{cos}{i}_{{\rm{o}}}-\mathrm{cos}{i}_{{\rm{s}}})}^{2}}},\\ & & \mathrm{for}\ \mathrm{all}\ | {i}_{{\rm{o}}}-\psi | \lt {i}_{{\rm{s}}}\lt \pi -| \pi -{i}_{{\rm{o}}}-\psi | .\end{eqnarray} \tag{ 9 }$$

Note that the angles i_s and ψ are allowed to vary in the interval $[0,\pi ]$ in the previous equation, whereas we restrict i_o to the interval $[0,\pi /2]$.

We now resort to Bayes' theorem to derive an analytical expression for $p(\psi | {i}_{{\rm{s}}},{i}_{{\rm{o}}})$:

$$\begin{eqnarray}&&p(\psi | {i}_{{\rm{s}}},{i}_{{\rm{o}}})\propto p({i}_{{\rm{s}}}| \psi ,{i}_{{\rm{o}}})p(\psi ),\end{eqnarray} \tag{ 10 }$$

where the prior probability $p(\psi )$ quantifies our assumptions on ψ. We will be adopting the uninformative prior probability $p(\psi )\propto \mathrm{sin}\psi$, which implies that ${{\boldsymbol{n}}}_{{\rm{s}}}$ and ${{\boldsymbol{n}}}_{{\rm{o}}}$ (see Figure 1) are uncorrelated. Finally, by multiplying Equation (A2) by $\mathrm{sin}\psi$ we arrive at an (unnormalized) expression for $p(\psi | {i}_{{\rm{s}}},{i}_{{\rm{o}}})$:

$$\begin{eqnarray}&&p(\psi | {i}_{{\rm{s}}},{i}_{{\rm{o}}})\propto \displaystyle \frac{\mathrm{sin}{i}_{{\rm{s}}}\mathrm{sin}\psi }{\sqrt{{\mathrm{sin}}^{2}\psi {\mathrm{sin}}^{2}{i}_{{\rm{o}}}-{(\mathrm{cos}\psi \mathrm{cos}{i}_{{\rm{o}}}-\mathrm{cos}{i}_{{\rm{s}}})}^{2}}},\\ &&\quad \mathrm{for}\ \mathrm{all}\ | {i}_{{\rm{o}}}-{i}_{{\rm{s}}}| \lt \psi \lt {i}_{{\rm{o}}}+{i}_{{\rm{s}}}.\end{eqnarray} \tag{ 11 }$$

Note that the angle ψ is allowed to vary in the interval $[0,\pi ]$ in the previous equation, whereas we restrict i_s and i_o to the interval $[0,\pi /2]$.

Extracting signatures of rotation from the power spectrum is accomplished by a detailed fitting of the observed modes of oscillation, a procedure often referred to as peak-bagging. We are primarily interested in performing a global fit to the power spectrum (e.g., Appourchaux et al. 2008, 2012; Campante et al. 2011 ), whereby a selection of the observed modes are fitted simultaneously over a broad frequency range. We modeled the limit (noise-free) oscillation spectrum as a series of standard Lorentzian profiles, ${\mathcal{O}}(\nu ;{{\boldsymbol{\Theta }}}_{{\rm{osc}}})$, which sit atop a background signal described by ${\mathcal{B}}(\nu ;{{\boldsymbol{\Theta }}}_{{\rm{bg}}}):$

$$\begin{eqnarray}{\mathcal{P}}(\nu ;{\boldsymbol{\Theta }}) & = & {\mathcal{O}}(\nu ;{{\boldsymbol{\Theta }}}_{{\rm{osc}}})+{\mathcal{B}}(\nu ;{{\boldsymbol{\Theta }}}_{{\rm{bg}}})\\ & = & \displaystyle \sum _{n^{\prime} ,l}\displaystyle \sum _{m=-l}^{l}\displaystyle \frac{{{\mathcal{E}}}_{{lm}}({i}_{{\rm{s}}}){H}_{n^{\prime} l}}{1+{\left[\displaystyle \frac{2(\nu -{\nu }_{n^{\prime} l0}-m{\nu }_{{\rm{s}}})}{{{\rm{\Gamma }}}_{n^{\prime} {lm}}}\right]}^{2}}+{\mathcal{B}}(\nu ;{{\boldsymbol{\Theta }}}_{{\rm{bg}}}).\end{eqnarray} \tag{ 12 }$$

The sets of mode and background model parameters are denoted by ${{\boldsymbol{\Theta }}}_{{\rm{osc}}}$ and ${{\boldsymbol{\Theta }}}_{{\rm{bg}}}$.

The inner sum in Equation (12) runs over the azimuthal components of each multiplet, while the outer sum runs over the selection of observed modes. The dummy variable n', which tags the observed order, takes values n' = n for l = 1, 2 modes, $n^{\prime} \;=\;n+1$ for l = 0 modes, and $n^{\prime} \;=\;n-1$ for l = 3 modes. To reduce the number of parameters entering our model, the heights and linewidths of non-radial modes are defined based on the heights and linewidths of the neighboring radial mode(s), the latter being allowed to vary with frequency. The set of mode parameters is ultimately given by ${{\boldsymbol{\Theta }}}_{{\rm{osc}}}\;=\;\{{\nu }_{n^{\prime} l0},{H}_{n^{\prime} 0},{{\rm{\Gamma }}}_{n^{\prime} 0},{i}_{{\rm{s}}},{\nu }_{{\rm{s}}}\}$.

The parameter ${H}_{n^{\prime} l}$ describes the height of a given mode of degree l, with the height of a multiplet component, ${H}_{n^{\prime} {lm}}$, given by ${{\mathcal{E}}}_{{lm}}({i}_{{\rm{s}}}){H}_{n^{\prime} l}$. ${H}_{n^{\prime} l}$ is defined relative to the height of the neighboring radial mode(s) according to ${H}_{n^{\prime} l}\;=\;{V}_{l}^{2}{H}_{n^{\prime} 0}$ for l = 2 modes or ${H}_{n^{\prime} l}\;=\;{V}_{l}^{2}[{H}_{n^{\prime} 0}+{H}_{(n^{\prime} -1)0}]/2$ for $l=1,3$ modes, where the V_l² describe the visibilities (in power) of modes of different l relative to those with l = 0 (we adopted fixed values of ${V}_{1}=1.22$, V₂ = 0.71, and V₃ = 0.14, obtained when taking into account the Kepler bandpass and the effect of limb darkening; Handberg & Campante 2011). The parameter ${{\rm{\Gamma }}}_{n^{\prime} {lm}}$ describes the mode linewidth, taken to be the same regardless of m for a given mode of degree l, i.e., ${{\rm{\Gamma }}}_{n^{\prime} {lm}}\equiv {{\rm{\Gamma }}}_{n^{\prime} l}$. ${{\rm{\Gamma }}}_{n^{\prime} l}$ is defined relative to the linewidth of the neighboring radial mode(s) according to ${{\rm{\Gamma }}}_{n^{\prime} l}\;=\;{{\rm{\Gamma }}}_{n^{\prime} 0}$ for l = 2 modes or ${{\rm{\Gamma }}}_{n^{\prime} l}\;=\;[{{\rm{\Gamma }}}_{n^{\prime} 0}+{{\rm{\Gamma }}}_{(n^{\prime} -1)0}]/2$ for $l=1,3$ modes.

The background signal was modeled as the superposition of three components:

$$\begin{eqnarray}&&{\mathcal{B}}(\nu ;{{\boldsymbol{\Theta }}}_{{\rm{bg}}})={B}_{0}+{\eta }^{2}(\nu )\left[\displaystyle \frac{{B}_{{\rm{gran}}}}{1+{(2\pi \nu {\tau }_{{\rm{gran}}})}^{a}}+\displaystyle \frac{{B}_{{\rm{act}}}}{{\nu }^{2}}\right].\end{eqnarray} \tag{ 13 }$$

A flat component B₀ is used to model photon shot noise. The contribution from granulation is described by a Harvey-like profile (Harvey & Plasma 1985; Harvey et al. 1993), where ${B}_{{\rm{gran}}}$ is the height at ν = 0 of the granulation component, τ_gran is the characteristic turnover timescale for granulation, and a calibrates the amount of memory in the process. Finally, a power-law component describes the contribution from activity (characterized by the scale factor B_act). This functional form results from considering a Harvey-like profile in the limit $2\pi \nu \;{\tau }_{{\rm{act}}}\gg 1$ with the exponent set to 2 as in the original work by Harvey & Plasma (1985). Such a power law is well suited to describing the decaying slope of the activity component. The attenuation factor η² is given by ${{\rm{sinc}}}^{2}\left[\frac{\pi }{2}\left(\frac{\nu }{{\nu }_{{\rm{Nyq}}}}\right)\right]$ for an integration duty cycle of 100% (e.g., Kallinger et al. 2014), where ${\nu }_{{\rm{Nyq}}}$ is the Nyquist frequency.

We started by fitting the background signal using a maximum likelihood estimator prior to the extraction of the mode parameters (e.g., Karoff et al. 2013, and references therein). The spectral range considered in this preliminary fit excluded the oscillation modes and the very low frequencies (Figure 16, top panel). Exclusion of the very low frequencies results from the power law being able to describe only the decaying slope of the activity component. The background model parameters, i.e., ${{\boldsymbol{\Theta }}}_{{\rm{bg}}}\;=\;\{{B}_{0},{B}_{{\rm{act}}},{B}_{{\rm{gran}}},{\tau }_{{\rm{gran}}},a\}$, determined in this way are subsequently held fixed in Equation (12) at their maximum likelihood estimates.

Mode parameters were then optimized in a Bayesian manner using an MCMC sampler that employs the Metropolis–Hastings algorithm (Figure 16, bottom panels; Handberg & Campante 2011; Campante 2012). In this way, we were able to map the joint PPD of ${{\boldsymbol{\Theta }}}_{{\rm{osc}}}$:

$$\begin{eqnarray}&&p({{\boldsymbol{\Theta }}}_{{\rm{osc}}}| D,I)\propto p({{\boldsymbol{\Theta }}}_{{\rm{osc}}}| I)p(D| {{\boldsymbol{\Theta }}}_{{\rm{osc}}},I),\end{eqnarray} \tag{ 14 }$$

where $p({{\boldsymbol{\Theta }}}_{{\rm{osc}}}| I)$ is the prior probability of ${{\boldsymbol{\Theta }}}_{{\rm{osc}}}$ based on any relevant prior information I and $p(D| {{\boldsymbol{\Theta }}}_{{\rm{osc}}},I)$ is the likelihood of the observed data D (a description of the likelihood function can be found in, e.g., Toutain & Appourchaux 1994). The PPD of a given mode parameter can then be simply arrived at through marginalization of the joint PPD.

A series of comprehensive fits to the power spectra of the stars in our sample have been conducted in a separate work (Davies et al. 2016). These fits considered the frequency range containing all observed modes and their output is used in the present analysis to define initial guesses for the mode parameters. Since we are mainly interested in estimating the stellar inclination angle, we focus here not on the full set of observed modes but instead on a selection of those modes, which necessarily span a somewhat narrower frequency range. The selection of the (continuous) frequency range used in the present analysis was done on a star-by-star basis, having taken into account the linewidths of the modes and their S/N as obtained from the above preliminary fits. This meant selecting those observed orders (a minimum of five) whose modes have the highest S/N and for which the intrinsic ${\nu }_{{\rm{s}}}/{\rm{\Gamma }}$ ratio is more favorable when it comes to resolving the azimuthal components (see Ballot et al. 2006). It should be noted that choosing the number of observed orders is a trade-off. The inclusion of more radial orders improves the precision in the final estimates. However, if orders are included where the modes show only modest or low S/N, improved precision may come at the cost of reduced accuracy (increased bias), rendering the outputs less robust. A measure of the S/N and ${\nu }_{{\rm{s}}}/{\rm{\Gamma }}$ at the frequency of maximum oscillation amplitude ${\nu }_{{\rm{max}}}$, as obtained from the present analysis, is given for each star in Table 3.

The main advantage of a Bayesian approach when compared to a frequentist approach is the ability to incorporate relevant prior information through Bayes' theorem (Equation (B3)). We imposed uniform priors on all mode parameters with the exception of ${i}_{{\rm{s}}}$, for which we assumed an isotropic orientation in the sphere (Abt 2001), i.e., $p({i}_{{\rm{s}}}| I)=\mathrm{sin}{i}_{{\rm{s}}}$ or, equivalently, a uniform prior on $\mathrm{cos}{i}_{{\rm{s}}}$. The stellar inclination angle was allowed to vary between $-\pi /2$ and π and all samples lying outside the interval $[0,\pi /2]$ were then reflected about the i_s = 0 and π/2 boundaries. This was done to avoid potential boundary effects in the MCMC sampling.

To validate our asteroseismic method, we have performed a series of tests with artificial data. The reduced visibility of the multiplet components in many of the stars in the asteroseismic sample (see Table 3) may raise concerns as to how robust the returned uncertainties on i_s are and how prone the MCMC fitting is to biases caused by the moderate-to-low S/N in the p modes. Fits to the power spectrum as described in Appendix B are computationally expensive. We have thus produced artificial power spectra for a representative set of stars in the asteroseismic sample, which include all six potentially misaligned systems (based on our results for i_s) plus Kepler-25. For each star, we have considered two input i_s with two noise realizations per input i_s, thus leading to a total of 28 fitted artificial spectra. For HAT-P-7 and Kepler-145 the input i_s was taken to be either 90° or 30°, whereas for the remainder of the stars we considered either 90° or 60°. Having fixed the input i_s, we then used the tabulated $v\mathrm{sin}{i}_{{\rm{s}}}$ and asteroseismic radius to determine the input ${\nu }_{{\rm{s}}}$. The power spectra were directly generated in the frequency domain, having preserved the same frequency resolution as in the real data and properly modeled the correlation between the background noise and the excitation function (see Chaplin et al. 2008). The adopted mode linewidths and S/N were based on the observed linewidths and S/N of radial modes (not split by rotation). We managed to retrieve the input i_s at the 1σ level in 21 of the fits (75%), at the 1.5σ level in 26 fits (∼93%), and at the 2σ level in 27 fits (∼96%). The returned uncertainties are based on the Bayesian HPD credible regions. Drawing a parallel between Bayesian credible regions and the more familiar frequentist confidence intervals, we are led to conclude that the returned uncertainties are robust.

We have also validated our asteroseismic method by performing a series of tests with degraded Sun-as-a-star data. Here, we are mostly interested in showing the effect of an increased temporal coverage in reducing the bias affecting the retrieved i_s. Using data provided by the green channel of the triple Sun PhotoMeter/Variability of solar IRradiance and Gravity Oscillations (SPM/VIRGO) instrument (Fröhlich et al. 1995) on board the Solar and Heliospheric Observatory (SOHO) spacecraft, with white noise added to levels comparable with Kepler-50, we have performed our asteroseismic analysis. This was done by treating the Sun as a Kepler target of magnitude m_Kep = 10 (and 10.5) observed during 270, 540, and 1080 days (the equivalent to 3, 6, and 12 Kepler quarters). The pristine time series were taken from around solar minimum and share the same starting time stamp. Figure 17 shows the output from the fits to the case with m_Kep = 10. The solar spin axis is inclined with respect to the normal of the ecliptic by $\sim 7^\circ$. The actual inclination for a particular block of SPM/VIRGO data will depend on the ephemerides of the Sun and is assumed to lie in the reference interval $[83^\circ ,90^\circ ]$. For reference, ${\nu }_{{\rm{s}}}/{\rm{\Gamma }}\approx 0.4$ at ν_max for the Sun. The bias affecting the retrieved i_s is reduced as we increase the temporal coverage (see Ballot et al. 2008). The solar inclination is retrieved at the 1σ level for the longest of the time series (or the equivalent to 12 Kepler quarters). We note that there are only three stars in the asteroseismic sample with a temporal coverage shorter than ten quarters (see Table 1): the posteriors of i_s for KOI-268 (Q6.1–Q8.3; Figure 18) and Kepler-129 (Q6.1–Q7.3; Figure 34) are dominated by the prior, whereas that of the bright (m_Kep = 8.72) star Kepler-444 (Q15.1–Q17.2; Figure 24) shows signs of bimodality.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The joint PPD of i_s and ν_s, as well as the corresponding marginalized PPDs, are shown in Figures 18–39 for the remainder of the stars in the asteroseismic sample.

We wish to sample the joint PPD of the model parameters ${\boldsymbol{\alpha }}$,

$$\begin{eqnarray}&&p({\boldsymbol{\alpha }}| {\{\mathrm{cos}{i}_{{\rm{s}}}\}}_{n=1}^{N})\propto p({\boldsymbol{\alpha }})\;{{ \mathcal L }}_{{\boldsymbol{\alpha }}},\end{eqnarray} \tag{ 15 }$$

where $p({\boldsymbol{\alpha }})$ is the prior on ${\boldsymbol{\alpha }}$ and ${{ \mathcal L }}_{{\boldsymbol{\alpha }}}$ is the marginalized likelihood for the model parameters. Here, we employ the hierarchical probabilistic method developed by Hogg et al. (2010), whereby the true distribution of a given quantity is inferred by taking as input a set of posterior samplings of that same quantity (obtained using an uninformative prior). According to this formalism, a sampling approximation to the marginalized likelihood for ${\boldsymbol{\alpha }}$ is given by

$$\begin{eqnarray}&&{{ \mathcal L }}_{{\boldsymbol{\alpha }}}\approx \displaystyle \prod _{n=1}^{N}\displaystyle \frac{1}{K}\displaystyle \sum _{k=1}^{K}\displaystyle \frac{{f}_{{\boldsymbol{\alpha }}}(\mathrm{cos}{i}_{{\rm{s}}({nk})})}{{p}_{0}(\mathrm{cos}{i}_{{\rm{s}}({nk})})},\end{eqnarray} \tag{ 16 }$$

where $\mathrm{cos}{i}_{{\rm{s}}({nk})}$ is the kth sample of the nth posterior of $\mathrm{cos}{i}_{{\rm{s}}}$. The sum is over K = 1000 posterior samples of each individual system, a number of samples large enough to guarantee that our analysis is not prone to sampling variance. The ratio inside the sum is between the model distribution and the uninformative prior on which the sampling was based (i.e., a uniform prior on $\mathrm{cos}{i}_{{\rm{s}}}$). Simply put, we are conducting a forward modeling of the observed data by defining a model distribution ${f}_{{\boldsymbol{\alpha }}}(\psi )$ for the true ψ and finding the model parameters ${\boldsymbol{\alpha }}$ that best explain those data. The use of input from several different observers/fitters, each of whom may have sampled their posteriors with different uninformative priors p₀, does not constitute a problem and can be straightforwardly accounted for in Equation (E2).