{'fraction_non_alphanumeric': 0.041739617620281534, 'fraction_numerical': 0.02881854471601653, 'mean_word_length': 3.9763024917232968, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Analytic Theory for the Orbits of Circumbinary Planets

Three transiting circumbinary planets (Kepler-16 b, Kepler-34 b, and Kepler-35 b) have recently been discovered from photometric data taken by the Kepler spacecraft. Their orbits are significantly non-Keplerian because of the large secondary-to-primary mass ratio and orbital eccentricity of the binaries, as well as the proximity of the planets to the binaries. We present an analytic theory, with the planet treated as a test particle, which shows that the planetary motion can be represented by the superposition of the circular motion of a guiding center, the forced oscillations due to the non-axisymmetric components of the binary's potential, the epicyclic motion, and the vertical motion. In this analytic theory, the periapse and ascending node of the planet precess at nearly equal rates in opposite directions. The largest forced oscillation term corresponds to a forced eccentricity (which is an explicit function of the parameters of the binary and of the guiding center radius of the planet), and the amplitude of the epicyclic motion (which is a free parameter of the theory) is the free eccentricity. Comparisons with direct numerical orbit integrations show that this analytic theory gives an accurate description of the planetary motion for all three Kepler systems. We find that all three Kepler circumbinary planets have nonzero free eccentricities.


INTRODUCTION
have recently discovered the first transiting circumbinary planet, Kepler-16 b, from photometric data taken by the Kepler spacecraft. The Saturn-mass planet orbits a pair of stars of 0.69M ⊙ and 0.20M ⊙ . Welsh et al. (2012) subsequently announced the discovery of two more circumbinary planets: Kepler-34 b and Kepler-35 b. Kepler-34 b is 0.22M J (where M J is the mass of Jupiter) and orbits two Sun-like stars, while Kepler-35 b is 0.13M J and orbits a pair of smaller stars (0.89M ⊙ and 0.81M ⊙ ). For all three systems, the orbits of the binary and planet are aligned to within 2 • . From the observed rate of circumbinary planets in their sample, Welsh et al. (2012) estimated that more than ∼ 1% of close binary stars have giant planets on nearly coplanar orbits.
Variations in eclipse times and transit durations, combined with radial velocity measurements, allow precise measurements of both physical and orbital parameters for all three systems. Table 1 Plots of the evolution of the osculating Keplerian orbital elements of the planets in the Supporting Online Material of Doyle et al. (2011) and Supplementary Information of Welsh et al. (2012) show significant variations on both orbital and secular timescales, with the eccentricity changing from nearly zero to 0.1 in the case of Kepler-16 b, and the precession period of the orbit is as short as ∼ 60 orbital periods in the case of Kepler-35 b (see figures below for more details). The nontrivial departures from unperturbed Keplerian orbits for these circumbinary planets are due to the large secondary-to-primary stellar mass ratio (m B /m A = 0.29-0.97), the large orbital eccentricity of the binary e AB = 0.14-0.52, and the proximity of the planet to the binary (orbital period ratio P b /P AB = 5.6-10.4).
In §2 we present an analytic theory for the orbit of a circumbinary planet that is valid in the limit that the planet has negligible mass and can be treated as a test particle. Lee & Peale (2006) have previously developed an analytic theory for the orbits of the small satellites of the Pluto-Charon system, assuming that the orbit of Charon relative to Pluto is circular. We generalize that theory to the case of an eccentric binary orbit. In §3 we present the results of direct numerical orbit integrations and compare them with the analytic theory. In §4 we discuss the limitations of the analytic theory and a simple modification that can significantly improve the analytic predictions when the orbital eccentricity of the binary is large. Our conclusions are summarized in §5.

ANALYTIC THEORY
In this section we follow a similar approach as Lee & Peale (2006) to develop an analytic theory for the orbit of a circumbinary planet. We extend their theory to first order in the eccentricity of the orbit of the binary. We assume that the planet has negligible mass and can be treated as a test particle. Then the orbit of the secondary (hereafter B) relative to the primary (hereafter A) is an elliptical Keplerian orbit with eccentricity e AB and semimajor axis a AB , and the distance between A and B is r AB = a AB (1 − e 2 AB )/(1 + e AB cos f B ), where f B is the true anomaly of B. We adopt a cylindrical coordinate system with the origin at the center of mass of the binary and the x-y plane being the orbital plane of the binary. The positions of B and A are r B = (R B , φ B , 0) and , m A is the mass of A, and m B is the mass of B.

Potential
The gravitational potential at r = (R, φ, z) due to the binary is Since the orbit of the planet is nearly coplanar with that of the binary, we expand 1/|r − r B | in powers of z: where The inverse powers of ρ can be expressed as cosine series using equation (6.62) of Murray & Dermott (1999) to give where δ k0 is the Kronecker delta and b k s (R B /R) are the Laplace coefficients. Similarly, Then the potential can be written as where To first order in e AB , where ̟ B and M B = n AB t + ϕ AB are, respectively, the longitude of periapse and mean anomaly of B relative A, n AB = [G(m A + m B )/a 3 AB ] 1/2 is the mean motion of the binary, and ϕ AB is a constant. Then To first order in e AB , and where a B = a AB m A /(m A + m B ), α B = a B /R, and D = d/dα. Equation (11) gives rise to terms involving which can be expressed as Terms involving R A /R can be found in a similar manner.
After re-grouping the terms, and -5 - The terms in equation (14) multiplied by e AB are new compared to the potential due to a binary on circular orbit derived by Lee & Peale (2006). As in the circular binary orbit case, the axisymmetric Φ j00 components of the potential are identical to those due to two rings: one of mass m A and radius a A and another of mass m B and radius a B .

Equations of Motion and Solution
The equations of motion in cylindrical coordinates arë As in Lee & Peale (2006), we approximate the orbit of the circumbinary planet as small deviations from the circular motion of a guiding center in the x-y plane: where the constant R 0 is the radius of the guiding center, |R 1 /R 0 | ≪ 1, |φ 1 | ≪ 1, and |z 1 /R 0 | ≪ 1.
Substituting equations (14) and (18) into equation (17), the only nontrivial equation at the zeroth order is which describes the circular motion of the guiding center. Its solution is where ϕ 0 is a constant and the mean motion n 0 is given by In the above equation n K = [G(m A + m B )/R 3 0 ] 1/2 is the Keplerian mean motion at R 0 , and α A and α B are evaluated at R = R 0 .
At the first order, the equations of motion arë where n = (R −1 dΦ 000 /dR) 1/2 is the mean motion at R, and the quantities in the square brackets with the subscript R 0 are evaluated at R = R 0 . Equation (23) can be integrated to giveφ 1 , which can then be substituted into equation (22) to yield where the epicyclic frequency κ 0 is given by Equation (25) is the equation of motion of a simple harmonic oscillator of natural frequency κ 0 that is driven at frequencies n AB , k|n 0 − n AB |, and |kn 0 − (k ± 1)n AB |, and its solution gives where e free and ψ are arbitrary constants and We can then solve equation (23) to give where The motion in R and φ is the superposition of the circular motion of the guiding center at R 0 at frequency n 0 , the epicyclic motion represented by the free eccentricity e free at frequency κ 0 , and the forced oscillations of fractional radial amplitudes C 0 , C 0 k and C ± k at frequencies n AB , k|n 0 − n AB | and |kn 0 − (k ± 1)n AB |, respectively. Note that The second and third combinations of frequencies correspond to the corotation and Lindblad resonances, respectively (e.g., Goldreich & Tremaine 1980). None of these resonances could be encountered if the planet is further away than the 3:1 mean-motion resonance with the binary.
For the motion in z, the solution to equation (24) is where i free and ζ are arbitrary constants and the vertical frequency ν 0 is given by Thus the motion in z decouples from that in R and φ and has only free oscillations represented by the free inclination i free at the vertical frequency ν 0 .
As we shall see, circumbinary planets have ν 0 > n 0 > n K > κ 0 . So the azimuthal period 2π/n 0 is shorter than the Keplerian orbital period 2π/n K , the periapse precesses in the prograde direction with the period 2π/|n 0 − κ 0 |, and the ascending node precesses in the retrograde direction with the period 2π/|n 0 − ν 0 |.
As in the circular binary orbit theory of Lee & Peale (2006), the motion is represented by the circular motion of the guiding center, the epicyclic motion, the forced oscillations and the vertical motion. The expressions for the mean motion n 0 (eq. [21]), the epicyclic frequency κ 0 (eq. [26]), and the vertical frequency ν 0 (eq. [36]) involve the axisymmetric Φ 000 and Φ 200 components of the potential and are identical to those in the circular binary orbit case (there are, however, corrections at the second order in e AB ; see §4). The motion in z is identical to the circular binary orbit case. The forced oscillations are composed of both terms identical to those in the circular binary orbit theory (C 0 k at frequencies k|n 0 − n AB |) and new terms (C 0 at frequency n AB and C ± k at frequencies |kn 0 − (k ± 1)n AB |). Note that the new terms C 0 and C ± k are proportional to e AB and that one of these terms, C − 1 , has frequency n 0 and can be identified as the forced eccentricity. The forced longitude of periapse is aligned with the binary's periapse because C − 1 ≥ 0.
If we expand the analytic expressions in powers of a AB /R 0 (note that a AB /R 0 0.32 for the three Kepler systems) and retain only the lowest order term, we find that the precession rates the forced eccentricity and the other forced oscillation terms (up to k = 3) We can see from equations (29) and (30) that C 0 k and C ± k involve Φ 0k0 , Φ 0k1 , and their derivatives with respect to R. According to equations (15) and (16)

COMPARISONS WITH NUMERICAL ORBIT INTEGRATIONS
For the numerical orbit integrations, we use Jacobi coordinates (where the position of the secondary B is relative to the primary A and the position of the planet is relative to the center of mass of the binary), with the invariable plane as the reference plane. This coordinate system reduces to that used in §2 when the mass of the planet is negligible. We perform direct numerical orbit integrations of the Kepler-16, 34, and 35 systems, using the Wisdom & Holman (1991) symplectic integrator with the modification described in Lee & Peale (2003). The modification allows the integration of circumbinary planets without an excessively small timestep, and we use a timestep of 0.1 day. We generate the initial positions and velocities of the binary and planet by using the best-fit osculating orbital parameters in Table 1.
For comparison with the analytic theory, we need to evaluate n K , n 0 , κ 0 , ν 0 , and the fractional amplitudes C 0 , C 0 k , and C ± k at a guiding center radius R 0 . For each system, we adopt the average of the maximum and minimum values of the cylindrical radius R of the planet's orbit in the numerical orbit integration over many precession cycles as R 0 . The adopted R 0 and the numerical values of n K , n 0 /n K , κ 0 /n K , ν 0 /n K , C 0 , C 0 k , and C ± k (k = 1, 2, and 3) are listed in Table 2.

Kepler-16
In Figure 1 we plot the evolution of the osculating Keplerian orbital elements of the planet Kepler-16 b over 100 years from the numerical orbit integration. The eccentricity e b shows variations on both orbital and apsidal precession timescales. The longitude of periapse ̟ b changes rapidly when e b is nearly zero, but the long-term trend is prograde precession with a period of 48.6 years. The precession of the longitude of ascending node Ω b is retrograde with a period of 41.0 years, and the inclination i b shows small oscillations around a constant value, with two oscillations per nodal precession period. Using n K , n 0 /n K , κ 0 /n K , and ν 0 /n K from Table 2, the analytic theory predicts that the apsidal and nodal precessions are at nearly equal rates in opposite directions, with the prograde apsidal precession having a period of 42.2 years and the retrograde nodal precession having a period of 42.8 years. These are in good agreement with the numerical results but slightly shorter for the apsidal precession and longer for the nodal precession.
In the bottom panels of Figure 2 we plot the variations in the orbital radius R b of the planet Kepler-16 b over 5.4 years and 100 years. Significant periodic variations in the amplitude of the oscillations in R b are observed from the 100-year plot. The variations have a period of 48.6 years, which is equal to the period of apsidal precession. The variations are the result of the superposition of the epicyclic motion at frequency κ 0 with amplitude e free and the forced oscillation at frequency n 0 with amplitude C − 1 . Without any loss of generality, we can assume that both e free and C − 1 ≥ 0. After some algebraic manipulation using sum and product formulae of trigonometric functions, these two terms in equation (27) can be written as where ̟ free = (n 0 − κ 0 )t + ϕ 0 − ψ is the free longitude of periapse. A maximum amplitude is reached when ̟ free − ̟ B = 2ℓπ (where ℓ is an integer), in which case the right-hand side of equation (43) becomes ±(C − 1 + e free ) cos{[(n 0 + κ 0 )t + ϕ 0 − ̟ B + ψ]/2}. Similarly, a minimum amplitude is reached when ̟ free − ̟ B = (2ℓ + 1)π, in which case the right-hand side of equation (43) becomes ±(C − 1 − e free ) sin{[(n 0 + κ 0 )t + ϕ 0 − ̟ B + ψ]/2}. Therefore, a maximum (or minimum) amplitude is reached every 2π/|n 0 − κ 0 |, which is equal to the apsidal precession period. The small minimum amplitude and large variations in the amplitude observed in the 100-year plot indicate that e free ∼ C − 1 . The 5.4-year plot in the lower left panel of Figure 2 clearly shows an increase in the amplitude of the oscillations in R b at the initial epoch due to the changing relative phase between the free and forced eccentricity terms, as well as higher-frequency forced oscillations. In order to study in more detail the forced oscillations and epicyclic motion, we plot in the upper panels of Figure 2 a transformed orbital radius defined by with R 0 , C 0 , C 0 k , and C ± k from Table 2 and φ 0 , M B , and ̟ B from the numerical integration itself (which eliminates phase errors due to small frequency errors and the very slow precession of the binary's periapse). It is clear from a comparison between the upper and lower panels of Figure 2 that the forced oscillations (including the forced eccentricity term) are sufficiently close to those predicted by the analytic theory that they are effectively eliminated in R ′ b . The largest forced oscillation term, other than the forced eccentricity term, is C − 2 with a fractional amplitude of 0.0024 and a period of 2π/(2n 0 − n AB ) = 64.5 days, and the other forced oscillation terms are at least a factor of 4 smaller in amplitudes. The free eccentricity can be easily obtained from R ′ b , which shows only periodic epicyclic variation at frequency κ 0 . The forced and free eccentricities of Kepler-16 b are C − 1 = 0.036 and e free = 0.030, respectively. Having obtained e free , which is a free parameter of the analytic theory, we can now directly plot the evolution of R according to equation (27) of the analytic theory (Fig. 3), as well as the evolution of the osculating Keplerian elements using R from equation (27), φ from equation (31), and their time derivatives (Fig. 4). They are in excellent agreement with the numerical orbit integration shown in Figures 1 and 2, except for the faster periapse precession (and hence faster periodic variations in the radial oscillation amplitude), the slower nodal precession, and the lack of periodic variations in the inclination. Our analytic theory only gives the free oscillations in the vertical direction and cannot explain the periodic variations in the inclination at twice the nodal precession rate. For the eccentricity variations, if we ignore the higher-frequency forced oscillations, equation (43) shows that the osculating eccentricity reaches a maximum of C − 1 + e free = 0.066 when the free longitude of periapse ̟ free is aligned with the forced longitude of periapse (which is equal to the longitude of periapse of the binary ̟ B ; see §2.2), and reaches a minimum of |C − 1 − e free | = 0.006 when the free longitude of periapse ̟ free is anti-aligned with the forced longitude of periapse. This behavior agrees with the usual definitions of the forced and free eccentricities and longitudes of periapse (see, e.g., Section 7.4 of Murray & Dermott 1999).

Kepler-34
Figures 5 and 6 show the evolution of the osculating Keplerian orbital elements, the orbital radius R b and the transformed orbital radius R ′ b of the planet Kepler-34 b. The periods of prograde apsidal precession and retrograde nodal precession are 62.9 and 67.9 years, respectively, from the numerical orbit integration. The analytic theory predicts 91.1 and 91.9 years, respectively, which are longer than the numerical results by more than 35%. The main reason for the large discrepancy is that the analytic theory is only accurate to first order in the binary eccentricity e AB and Kepler-34 has the largest e AB (= 0.52) among the three systems. We shall derive in §4 simple corrections at the second order in e AB that significantly improve the analytic predictions of the precession periods.
One might think that the large osculating Keplerian eccentricity (e b ∼ 0.2) is due to forcing by the eccentric binary. But the nearly identical plots of R b and R ′ b show that the variations in R b are primarily due to epicyclic motion and that the forced eccentricity C − 1 and other forced oscillation terms are small. Indeed, C − 1 = 0.0019, which is smaller than that for Kepler-16 b by more than an order of magnitude, and the next largest forced oscillation term is C − 2 = 0.00068 (see Table 2). The forced eccentricity C − 1 , as well as all C 0 k and C ± k terms with odd k, are small because of the nearly equal masses of the binary components of Kepler-34 (m B /m A = 0.97; see discussion in the last paragraph of §2.2). We find from the variations in R ′ b that the free eccentricity e free = 0.204.

Kepler-35
For Kepler-35 which has binary eccentricity (e AB = 0.14) similar to Kepler-16, the numerical and analytic apsidal and nodal precession periods of the planet are in good agreement. The numerically determined apsidal and nodal precession periods are 21.7 and 20.2 years, respectively (see Fig. 7), and the analytic ones are 20.4 and 20.8 years, respectively.
As in the case of Kepler-34, the binary components have nearly equal masses (m B /m A = 0.91) and the forced eccentricity of the planet C − 1 = 0.0025 is small. However, because the free eccentricity is much smaller than that for Kepler-34 b and comparable to that for Kepler-16 b, moderate variations in the amplitude of oscillations in R b with the same period as the apsidal precession are clearly observed in the 100-year plot in the lower right panel of Figure 8. The 4-year plot in the lower left panel of Figure 8 also shows the effects of higher-frequency forced oscillations terms. As for Kepler-16 b and 34 b, C − 2 is the largest forced oscillation term after the forced eccentricity term C − 1 (see Table 2). The forced oscillations are sufficiently close to those predicted by the analytic theory that the transformed orbital radius R ′ b shows only periodic epicyclic variation at frequency κ 0 (see upper panels of Fig. 8). The free eccentricity from the variations in R ′ b is e free = 0.038.

DISCUSSION
The analytic theory developed in §2 is accurate to first order in the binary eccentricity e AB and to first order in the deviations R 1 , φ 1 , and z 1 of the planetary motion from the circular motion of the guiding center. It also treats the planet as a test particle and ignores the gravitational effects of the planet on the motion of the binary. From the comparisons with direct numerical orbit integrations of the Kepler-16, 34, and 35 systems in §3, we have shown that the analytic theory gives an accurate description of the planetary motion in all three cases, except for the apsidal and nodal precession periods of Kepler-34 b with e AB = 0.52.
It was pointed out in §2 that the expressions for the mean motion n 0 (eq. [21]), the epicyclic frequency κ 0 (eq. [26]), and the vertical frequency ν 0 (eq. [36]) involve the axisymmetric Φ 000 and Φ 200 components of the potential and that the axisymmetric components of the potential are identical to those due to two rings: one of mass m A and radius a A and another of mass m B and radius a B . If we expand R B /a B to higher orders in e AB (see eq. [2.81] of Murray & Dermott 1999), which means that the time-averaged R B /a B = 1 + e 2 AB /2 and that it might be more appropriate for the axisymmetric components of the potential to be due to two rings of radii a A (1 + e 2 AB /2) and a B (1 + e 2 AB /2). This is achieved if we selectively include just the e 2 AB /2 term beyond first order in e AB and use R B /a B ≈ 1 − e AB cos M B + e 2 AB /2 (and similarly for R A /a A ). Then the only modifications to the analytic theory in §2 are that (1 + e 2 AB /2)], in equation (15) for Φ jk0 , and similarly for Db k (j+1)/2 (α A ) and Db k (j+1)/2 (α B ) in equation (16) for Φ jk1 . With this simple modification, the analytic predictions for the apsidal and nodal precession rates are faster than the unmodified values by only 2-3% for Kepler-16 b and 35 b, but by ∼ 29% for Kepler-34 b. The increase by approximately (1 + e 2 AB /2) 2 can be understood from the (a AB /R 0 ) 2 scaling of the lowest order expression for the precession rates in equation (37). The modified analytic precession periods for Kepler-34 b are 71.4 years for periapse and 72.1 years for ascending node, which are much closer to the numerical results (62.9 and 67.9 years, respectively) than the unmodified values (∼ 91 years).
The 1 + e 2 AB /2 modification also affects the amplitudes of the forced oscillation terms. The change in the largest of these, the forced eccentricity C − 1 , is small: 6% even for Kepler-34 b. The change in the second largest forced oscillation term, C − 2 , is approximately 1 + 5e 2 AB /6, which is only 2% for Kepler-16 b and 35 b but ∼ 24% for Kepler-34 b. However, as we saw in §3.2, both C − 1 and C − 2 are small compared to the free eccentricity and have no noticeable effect on the evolution of R b for Kepler-34 b.
The 1 + e 2 AB /2 modification that we have just described is not rigorously correct. We have attempted to construct an analytic theory that is accurate to O(e 2 AB ). Preliminary analysis indicates that there are no additional corrections at O(e 2 AB ) for C − 2 but that there are additional corrections to the precession periods and C − 1 . The full O(e 2 AB ) corrections to the precession periods and C − 1 of Kepler-16 b and Kepler-35 b, as well as C − 1 of Kepler-34 b, remain small (less than a few %). For Kepler-35 b, the apsidal (nodal) precession period may decrease (increase) by a few years beyond the 1 + e 2 AB /2 modification discussed above. The full O(e 2 AB ) theory also introduces new forced oscillation terms with frequencies |kn 0 − (k ± 2)n AB |, but they are likely small in amplitudes, as they are not observed in the direct numerical orbit integrations.
The most obvious effects of the gravitational force of the planet on the binary are the precession of the binary's periapse and ascending node. Due to our choice of the invariable plane as the reference plane, the longitude of ascending node of the binary must be 180 • from, and precesses at the same rate as, the longitude of ascending node of the planet. From the direct numerical orbit integrations, we find that the apsidal precession rates of the binary are 0.026, 0.0033 and 0.0086 degrees per year for Kepler-16, 34, and 35, respectively, which are much smaller than those of the planet. For the comparisons in §3, this very slow precession of the binary's periapse is taken into account by using ̟ B from the numerical integrations when we plot R ′ defined in equation (44).

CONCLUSIONS
We have developed an analytic theory to model the motion of the recently discovered circumbinary planets: Kepler-16 b, Kepler-34 b, and Kepler-35 b. Their orbits are significantly non-Keplerian due to the large secondary-to-primary mass ratio and orbital eccentricity of the binaries, as well as the proximity of the planets to the binaries. The analytic theory in §2 shows that the motion in R and φ can be represented by the superposition of the circular motion of a guiding center at mean motion n 0 , the epicyclic motion, and the forced oscillations, and that the motion in z decouples from that in R and φ and has only free oscillations. One of the forced oscillation terms has frequency n 0 and can be identified as the forced eccentricity, while the epicyclic motion at frequency κ 0 can be identified as the free eccentricity.
Comparisons with direct numerical orbit integrations in §3 show that the analytic theory (with the simple modification in §4) gives an accurate description of the planetary motion for all three Kepler systems, including the precession of the periapse and ascending node. The analytic theory explains the periodic variations in the amplitude of oscillations of the orbital radius (which is most obvious for Kepler-16 b and negligible for Kepler-34 b) by the superposition of the epicyclic motion and the forced eccentricity oscillation. The amplitude (and osculating eccentricity) variations have a period equal to that of the apsidal precession as predicted by the theory.
The amplitude, C − 1 , of the forced eccentricity term is an explicit function of the parameters of the binary and of the guiding center radius of the planet in the analytic theory. For Kepler-16 b, 34 b, and 35 b, C − 1 = 0.036, 0.0019, and 0.0025, respectively. The free eccentricity, e free , of the epicyclic motion is a free parameter in the analytic theory and can be obtained from, e.g., the variations in the orbital radius of the planet. For Kepler-16 b, 34 b, and 35 b, e free = 0.030, 0.204, and 0.038, respectively. Note that the Kepler-34 system with the largest binary eccentricity (e AB = 0.52) has the largest e free while the other two Kepler systems with comparable e AB have comparable e free . Since the free eccentricity is a free parameter that was set by dynamical processes during the formation and/or subsequent evolution of the circumbinary planet, the free eccentricity of the three Kepler circumbinary planets (and additional circumbinary planets in the future) should provide important clues to these processes.
While this paper was under review, three more circumbinary planetary systems were an-nounced: Kepler-38 and PH1 with one planet each and Kepler-47 with two planets (Orosz et al. 2012a,b;Schwamb et al. 2012). Direct numerical integrations and our analytic theory show that (i) Kepler-38 b is similar to Kepler-16 b in having e free ∼ C − 1 ∼ 0.024 and large variations in the amplitude of oscillations in R; (ii) Kepler-47 b is similar to Kepler-35 b in having e free larger than C − 1 ∼ 0.004 and moderate variations in the amplitude of oscillations in R; (iii) Kepler-47 c is similar to Kepler-34 in having e free much larger than C − 1 ∼ 0.001 and negligible variations in the amplitude of oscillations in R; and (iv) PH1 has e free ∼ 0.1 and C − 1 ∼ 0.04.
It is a pleasure to thank Josh Carter for furnishing the best-fit orbital parameters of the Kepler circumbinary planetary systems, Daniel Fabrycky for informative discussions, and the referee for helpful comments on the manuscript. C. K. L. also thanks K. H. Chan, X. Tan, and W. Y. Li for enlightening discussion. This work was supported in part by Hong Kong RGC grant HKU 7034/09P. Note. -The orbital parameters are the best-fit osculating Jacobian parameters relative to the invariable plane at the listed epoch.