Extremely Inclined Orbit of the S-type Planet γ Cep Ab Induced by the Eccentric Kozai–Lidov Mechanism

The γ Cep system is known as a close binary at a distance of 13.79 pc (Hatzes et al. 2003), and it is a candidate target of the CHES mission (Ji et al. 2022). The primary star γ Cep A is a planet-hosting bright star of spectral type K1III–IV, with a stellar mass of m₀ = 1.40 ± 0.12 M_⊙ (Neuhäuser et al. 2007). Neuhäuser et al. (2007) presented a direct detection of the companion γ Cep B, where the parameters of the secondary star are m₂ = 0.409 ± 0.018 M_⊙, a₂ = 20.18 ± 0.66 au, I₂ = 1193, Ω₂ = 18.04 ± 098, and the orbital period is P₂ = 67.5 ± 1.4 yr. Reffert & Quirrenbach (2011) conducted a fitting only for the inclination I₁ and the ascending node Ω₁, adopting P, a, e, and K from the best-fitting solution from Butler et al. (2006). They thereby obtained the best-fitting values of I₁ = 57, Ω₁ = 375 or I₁ = 1731, Ω₁ = 3561 (see Table 1), where I₂ and Ω₂ are constrained by Neuhäuser et al. (2007).

Table 1. Orbital Solutions of the γ Cep System

Object	Parameters	This Work		Hatzes et al. (2003)	Torres (2007)	Neuhäuser et al. (2007)
		χ2 = 1.48	χ2 = 1.44
	K1 (m s−1)	${28.08}_{-1.23}^{+1.45}$	${26.40}_{-1.18}^{+1.41}$	27.5 ± 1.5	27.1 ± 1.5	27.0 ± 1.5
	P1 (days)	905.64 ± 2.83	${901.46}_{-2.83}^{+2.85}$	905.574 ± 3.08	902.8 ± 3.5	902.9 ± 3.5
	a1 (AU)	2.1459 ± 0.0048	${2.1376}_{-0.0047}^{+0.0048}$	2.13 ± 0.05	1.94 ± 0.06	2.044 ± 0.057
	e1	${0.0724}_{-0.0575}^{+0.0879}$	${0.0856}_{-0.0624}^{+0.0866}$	0.12 ± 0.05	0.113 ± 0.058	0.115 ± 0.058
γ Cep Ab	ω1 (deg)	${48.47}_{-1.81}^{+4.56}$	${55.37}_{-4.57}^{+8.81}$	49.6 ± 25.6	63.0 ± 27.0	63.0 ± 27.0
	T0,1 (HJD-2400000)	${53140.16}_{-38.31}^{+34.15}$	${53107.63}_{-30.90}^{+25.47}$	53121.9 ± 66.9	53146.0 ± 72.0	53146.0 ± 71.0
	${m}_{1}\sin {I}_{1}\left({M}_{\mathrm{Jup}}\right)$	${1.7420}_{-0.0743}^{+0.0733}$	${1.6306}_{-0.0694}^{+0.0698}$	1.70 ± 0.40	1.43 ± 0.13	1.60 ± 0.13
	Ω1 (deg)	37.5 a	356.1 a	...	...	...
	I1 (deg)	5.7 a	173.1 a	...	>4.9	...
	K2 (m s−1)	${1821.70}_{-1.01}^{+0.95}$	${1699.94}_{-3.19}^{+3.45}$	1820.0 ± 49.0	1925.0 ± 14.0	1932.0 ± 14.0
	P2 (days)	${20754.66}_{-57.09}^{+59.01}$	${20731.68}_{-59.53}^{+57.18}$	20750.658 ± 1568.6	24392 ± 522	24654.375 ± 511
	a2 (AU)	${18.6421}_{-0.0381}^{+0.0392}$	${18.6217}_{-0.0389}^{+0.0381}$	18.5 ± 1.1	19.02 ± 0.64	20.18 ± 0.66
	e2	0.3603 ± 0.0012	0.3605 ± 0.0026	0.3610 ± 0.0230	0.4085 ± 0.0065	0.4112 ± 0.0063
γ Cep AB	ω2 (deg)	158.92 ± 0.20	158.90 ± 0.20	158.76 ± 1.20	160.96 ± 0.40	161.01 ± 0.40
	T0,2 (HJD-2400000)	${48435.22}_{-0.37}^{+0.80}$	${48435.04}_{-0.88}^{+0.46}$	48429.03 ± 27.0	48479.0 ± 12.0	48444.8 ± 11.16
	${m}_{2}\left({M}_{\odot }\right)$	0.3991 ± 0.0005	0.3986 ± 0.0005	...	0.362 ± 0.022	0.409 ± 0.018
	Ω2 (deg)	18.04 ± 0.98 b	18.04 ± 0.98 b	...	13.0 ± 2.4	18.04 ± 0.98
	I2 (deg)	119.3 ± 1.0 b	119.3 ± 1.0 b	...	118.1 ± 1.2	119.3 ± 1.0

Notes.

^a Astrometric data fitting results adopted from Reffert & Quirrenbach (2011). ^b Orbital solutions adopted from the imaging observations in Neuhäuser et al. (2007).

Download table as:  ASCIITypeset image

By using the most recent spectroscopic observations from the Canada–France–Hawaii Telescope (CFHT; Walker et al. 1992) and the McDonald Observatory Planetary Search (MOPS) program (Hatzes et al. 2003), two-Keplerian solutions of the γ Cep system have previously been given (Hatzes et al. 2003; Neuhäuser et al. 2007; Torres 2007). In these studies, an uncoupled two-Keplerian orbital fitting approach is utilized with the standard nonlinear least-squares technique. However, for a hierarchical system with a massive secondary companion, the significant mutual interaction between the planet and the star in the pair is supposed to be considered under the Jacobi frame, in order for the real dynamics to be represented (Lee & Peale 2003).

Here, we solve the orbits of the planet and the star companion with the MCMC ensemble sampler emcee (Foreman-Mackey et al. 2013) in the full N-body model (Ford 2006; Nelson et al. 2016). 14 parameters of $\{{m}_{1}\sin {I}_{1}$, P₁, $\sqrt{{e}_{1}}\sin {\omega }_{1}$, $\sqrt{{e}_{1}}\cos {\omega }_{1}$, M₁, ${m}_{2}\sin {I}_{2}$, P₂, $\sqrt{{e}_{2}}\sin {\omega }_{2}$, $\sqrt{{e}_{2}}\cos {\omega }_{2}$, M₂}, plus the RV offsets {RV_0,1, RV_0,2, RV_0,3, RV_0,4} of four time series (CFHT, MOPS I, MOPS II, and MOPS III), are adopted for fitting at the first observation epoch (HJD-2444754.129) in the Jacobi reference frame (Lee & Peale 2003). $\sqrt{e}\sin \omega$ and $\sqrt{e}\cos \omega$ are more efficient than e and ω for low-eccentricity planets, and can avoid a situation of multiple solutions (Ford 2006).

We then utilize the N-body integrator IAS15 (Rein & Spiegel 2015), which is an integrator with an adaptive step-size control, to compute the perturbed orbits of the planet and the secondary star in each step of the emcee sampler. The initial orbital inclinations and ascending nodes are assumed to be known (Reffert & Quirrenbach 2011): I₁ = 57, Ω₁ = 375 or I₁ = 1731, Ω₁ = 3561. Other initials include the planetary mass m₁ and the secondary mass m₂, the SMA a_1,2, the eccentricity e_1,2, and the argument of periastron ω_1,2, which can be derived from the resultant fitting parameters. The integration precision is given as 10⁻¹⁵ and the minimum time step is set to 0.001 days. The theoretical RV signals of the observation epochs are calculated from the Jacobian orbital elements, with the RV semi-amplitude K_1,2 being defined in Lee & Peale (2003):

$$\begin{eqnarray}&&\begin{array}{l}{K}_{1}={\left(\displaystyle \frac{2\pi G}{{P}_{1}}\right)}^{1/3}\displaystyle \frac{{m}_{1}\sin {I}_{1}}{{\left({m}_{0}+{m}_{1}\right)}^{2/3}}\displaystyle \frac{1}{\sqrt{1-{e}_{1}^{2}}},\\ {K}_{2}={\left(\displaystyle \frac{2\pi G}{{P}_{2}}\right)}^{1/3}\displaystyle \frac{{m}_{2}\sin {I}_{2}}{{\left({m}_{0}+{m}_{1}+{m}_{2}\right)}^{2/3}}\displaystyle \frac{1}{\sqrt{1-{e}_{2}^{2}}},\end{array}\end{eqnarray} \tag{ 1 }$$

where P_1,2, I_1,2, and e_1,2 are time-variable orbital elements at each observation epoch.

To derive the N-body best-fitting solutions for this system, we run 100,000 steps each for 28 walkers in the 14-dimensional parameter space. The mean acceptance fraction of the sampler walkers is 0.36, which is in the range from 0.2 to 0.5 that is suggested to produce representative samples (Foreman-Mackey et al. 2013), and the autocorrelation time for each fitting parameter ranges from 200 to 500 steps, which is smaller than our number of sampling steps to ensure the convergence of the MCMC chains. Here, the autocorrelation time represents the number of steps required for the chain to produce an independent sample (Foreman-Mackey et al. 2013).

Finally, we report two set of solutions of the minimum mass $m\sin I$, the RV semi-amplitude K, the SMA a, the eccentricity e, the argument of periastron ω, and the epoch of periastron passage T₀ in Table 1, with χ² = 1.48 for I₁ = 57 and χ² = 1.44 for I₁ = 1731, respectively. In Figures 2 and 3, we plot the best-fitting RV signals from the N-body solutions (the black solid lines), using the measurements versus the epoch with χ² = 1.44, with respect to ${{RV}}_{\mathrm{0,1}}={1290.0}_{-17.6}^{+18.5}$ ms⁻¹, ${{RV}}_{\mathrm{0,2}}={2035.2}_{-18.8}^{+20.3}$ ms⁻¹, ${{RV}}_{\mathrm{0,3}}={2231.1}_{-16.8}^{+19.2}$ ms⁻¹, and ${{RV}}_{\mathrm{0,4}}={865.0}_{-19.4}^{+22.4}$ ms⁻¹, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By subtracting the theoretical RV induced by the secondary, Figure 4 presents the Lomb–Scargle periodogram analysis of the planet b, where the strongest signal reveals an orbital period of roughly 900 days. Figure 5 shows the posterior distributions of $\{{m}_{1}\sin {I}_{1}$, P₁, $\sqrt{{e}_{1}}\sin {\omega }_{1}$, $\sqrt{{e}_{1}}\cos {\omega }_{1}$, M₁, ${m}_{2}\sin {I}_{2}$, P₂, $\sqrt{{e}_{2}}\sin {\omega }_{2}$, $\sqrt{{e}_{2}}\cos {\omega }_{2}$, M₂}. We report the median (50th percentile) of the posterior distribution as the best-fitting value, with the 1σ uncertainties being provided using the 16th and 84th percentiles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To derive more reliable orbital solutions, it is necessary to utilize high-precision RV and astrometric measurements in the fitting procedure at the same time. As is well known, Gaia Data Release 2 (Gaia Collaboration et al. 2018) was first accessible in 2018, and the science team recently announced Gaia Early Data Release 3 (Gaia Collaboration et al. 2021). For most of the sources in the Gaia catalogs, sequential astrometric data are not yet available, and will be released in Gaia Data Release 4. Thus, improvements of the planetary masses from Gaia will be expected in the future.

Table 1 summarizes our derived orbital parameters of the γ Cep system. The minimum planetary mass is fitted as ${1.7420}_{-0.0743}^{+0.0733}$ M_Jup or ${1.6306}_{-0.0694}^{+0.0698}$ M_Jup. Note that the real mass of planet b is significantly dependent upon the accuracy of the RV observations and the dynamical integrations; here, we simply estimate the mass range of γ Cep Ab. When the observed inclination of the planet is I₁ ∈ [38, 208] (Reffert & Quirrenbach 2011), we obtain an estimated planetary mass m₁ ∈ [5.0, 26.6] M_Jup, while for I₁ ∈ [1666, 1748] (Reffert & Quirrenbach 2011), m₁ ∈ [7.1, 26.2] M_Jup can be derived. In addition, based on the planetary mass limit of 16.9 M_Jup (Torres 2007), we here assume m₁ ∈ [5, 16.9] M_Jup.

To determine the current observed mutual inclination i_mut in the γ Cep system, we use the law of cosines for angles of a spherical triangle (Gellert et al. 1977):

$$\begin{eqnarray}&&\cos {i}_{\mathrm{mut}}=\cos {I}_{1}\cos {I}_{2}+\sin {I}_{1}\sin {I}_{2}\cos ({{\rm{\Omega }}}_{1}-{{\rm{\Omega }}}_{2}),\end{eqnarray} \tag{ 2 }$$

deriving i_mut = 1139 for I₁ = 57. Thus, a question naturally arises—what kind of dynamical process could trigger such high mutual inclinations of the two bodies, and could the EKL mechanism play a role in the secular evolution of γ Cep Ab? In the following section, we describe the EKL mechanism and explore the secular evolution of γ Cep Ab under EKL.

In the secular evolution of triple-body systems, the perturbation arising from the third body acts on a timescale much longer than its orbital period. One typical origin of secular resonance is from the perturbation potential from adjacent orbits. When the mass of the perturbed body is negligibly small in comparison to those bodies in a hierarchical system, the test particle approximation comes into play. At this time, if the perturbation of a circular outer orbit works, the Hamiltonian of this system can be expanded to the quadrupole level, known as the KL mechanism, as aforementioned.

The KL mechanism addresses the inclination and eccentricity of the inner test particle oscillating over secular evolution in a hierarchical system, where the test body and the primary star are surrounded by a distant companion (von Zeipel 1910; Kozai 1962; Lidov 1962). With the secular approximation, the inner and outer orbits only exchange angular momentum, so the SMAs of the orbits do not change. When the SMA ratio α_r = a₁/a₂ is a small parameter, the perturbation term of the complete Hamiltonian can be expanded as a power series in α_r (Naoz 2016):

$$\begin{eqnarray}\displaystyle \begin{array}{rcl}{ \mathcal H } & = & \displaystyle \frac{{k}^{2}{m}_{0}{m}_{1}}{2{a}_{1}}+\displaystyle \frac{{k}^{2}{m}_{2}\left({m}_{0}+{m}_{1}\right)}{2{a}_{2}}\\ & & +\,\displaystyle \frac{{k}^{2}}{{a}_{2}}\sum _{j=2}^{\infty }{\alpha }_{r}^{j}{M}_{j}{\left(\displaystyle \frac{{r}_{1}}{{a}_{1}}\right)}^{j}{\left(\displaystyle \frac{{a}_{2}}{{r}_{2}}\right)}^{j+1}{P}_{j}(\cos {\rm{\Phi }}),\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{M}_{j}={m}_{0}{m}_{1}{m}_{2}\displaystyle \frac{{m}_{0}^{j-1}-{\left(-{m}_{1}\right)}^{j-1}}{{\left({m}_{0}+{m}_{1}\right)}^{j}},\end{eqnarray} \tag{ 4 }$$

where m₀ is the mass of the primary, m₁ and m₂ are the masses of the inner and the outer body, k² is the gravitational constant (with the mass unit of M_⊙ and the length unit of au), r₁ is the distance between m₀ and m₁, r₂ is the distance between the center of mass of the inner binary and m₂, P_j is the Legendre polynomial, and Φ is the angle between the vectors r₁ and r₂ (the subscripts j = 1, 2 represent the inner and outer orbits, respectively).

The invariable plane reference frame is imported here to define three important Delaunay elements (Valtonen & Karttunen 2006): l, g, h, and their conjugate momenta L, G, H, where l, g, and h are the notations of the mean anomaly M, the argument of periastron ω, and the longitude of ascending node Ω, respectively. As shown in Figure 6, G _tot is the total angular momentum vector of the system, which is conserved over the secular evolution.

Zoom In Zoom Out Reset image size

The three conjugate momenta L, G, H are expressed as (Naoz et al. 2013):

$$\begin{eqnarray}\begin{array}{rcl}{L}_{1} & = & \displaystyle \frac{{m}_{0}{m}_{1}}{{m}_{0}+{m}_{1}}\sqrt{{k}^{2}\left({m}_{0}+{m}_{1}\right){a}_{1}},\\ {L}_{2} & = & \displaystyle \frac{{m}_{2}\left({m}_{0}+{m}_{1}\right)}{{m}_{0}+{m}_{1}+{m}_{2}}\sqrt{{k}^{2}\left({m}_{0}+{m}_{1}+{m}_{2}\right){a}_{2}},\\ {G}_{1} & = & {L}_{1}\sqrt{1-{e}_{1}^{2}},\quad {G}_{2}={L}_{2}\sqrt{1-{e}_{2}^{2}},\\ {H}_{1} & = & {G}_{1}\cos {i}_{1},\quad {H}_{2}={G}_{2}\cos {i}_{2},\\ {G}_{\mathrm{tot}} & = & {H}_{1}+{H}_{2},\end{array}\end{eqnarray} \tag{ 5 }$$

where L is only determined by constant parameters, including the masses m₀, m₁, and m₂, the SMAs a₁, a₂, and the gravitational constant k². Thus, L is a constant for a specific system in the evolution, while G and H are time-varying. G represents the magnitude of each orbit's angular momentum and H is the component of G along the z-axis.

According to the geometric relations and the assumption of h₁ − h₂ = π, the mutual inclination between the inner and outer orbit i_mut could be derived as (Naoz et al. 2013):

$$\begin{eqnarray}&&\cos {i}_{\mathrm{mut}}=\cos ({i}_{1}+{i}_{2})=\displaystyle \frac{{G}_{\mathrm{tot}}^{2}-{G}_{1}^{2}-{G}_{2}^{2}}{2{G}_{1}{G}_{2}}.\end{eqnarray} \tag{ 6 }$$

Generally, the equations of motion can be expressed by canonical relations of the three conjugate momenta, L, G, H, and the three Delaunay elements, l, g, h, as the mean anomaly can be eliminated under the double-averaged secular approximation, and h₁ and h₂ in equations of motion have been removed by the relation h₁ − h₂ = π. The time evolution for ω, e, and i can easily be derived from the reduced canonical relations (Naoz 2016):

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{G}_{j}}{\,{\rm{d}}t}=\displaystyle \frac{\partial { \mathcal H }}{\partial {g}_{j}},\quad \displaystyle \frac{{\rm{d}}{g}_{j}}{\,{\rm{d}}t}=-\displaystyle \frac{\partial { \mathcal H }}{\partial {G}_{j}},\end{eqnarray} \tag{ 7 }$$

where j = 1, 2. The original nonplanar three-body system can thus be reduced to a dynamical system of 2 degrees of freedom (dof).

In Naoz et al. (2013), i₂ is set to be 0, thus i_mut is equal to i₁. In this work, we treat the orbital plane of the secondary as the invariable plane, where i_mut = i₁. In the following sections, we redefine the orbital flip of the planet b as the variation of i_mut around 90°, instead of the real observed orbital inclination.

In the γ Cep Ab B system, e₂ is close to 0.36 and the estimated maximum mass of the planet is near the deuterium-burning limit. Thus, γ Cep Ab B is a nonrestricted triple system with an eccentric outer orbit, where the octupole-level terms in the Hamiltonian can become important, and the eccentricities of the two orbits are coupled and oscillate simultaneously over secular evolution. Naoz et al. (2013) called it the EKL mechanism in nonrestricted triple systems and derived the complete Hamiltonian of the system, including the octupole term in addition to the quadrupole term:

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal H }}_{3}({\rm{\Delta }}h\to \pi ) & = & {{ \mathcal H }}_{\mathrm{quad}}+{{ \mathcal H }}_{\mathrm{oct}}\\ & = & {C}_{2}[(2+3{e}_{1}^{2})(3{\cos }^{2}{i}_{\mathrm{mut}}-1)\\ & & +\,15{e}_{1}^{2}{\sin }^{2}{i}_{\mathrm{mut}}\cos (2{g}_{1})]\\ & & +\,{C}_{3}{e}_{1}{e}_{2}[A\cos \phi +10\cos {i}_{\mathrm{mut}}{\sin }^{2}{i}_{\mathrm{mut}}\\ & & \times \,(1-{e}_{1}^{2})\sin {g}_{1}\sin {g}_{2}],\end{array}\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&{C}_{2}=-\displaystyle \frac{{k}^{4}}{16}\displaystyle \frac{{\left({m}_{0}+{m}_{1}\right)}^{7}}{{\left({m}_{0}+{m}_{1}+{m}_{2}\right)}^{3}}\displaystyle \frac{{m}_{2}^{7}}{{\left({m}_{0}{m}_{1}\right)}^{3}}\displaystyle \frac{{L}_{1}^{4}}{{L}_{2}^{3}{G}_{2}^{3}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\epsilon }_{M}=\left(\displaystyle \frac{{m}_{0}-{m}_{1}}{{m}_{0}+{m}_{1}}\right)\left(\displaystyle \frac{{a}_{1}}{{a}_{2}}\right)\displaystyle \frac{{e}_{2}}{1-{e}_{2}^{2}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{C}_{3}=\displaystyle \frac{15}{4}{\epsilon }_{M}{C}_{2},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}\begin{array}{rcl}A & = & 4+3{e}_{1}^{2}-\displaystyle \frac{5}{2}B\sin {i}_{\mathrm{mut}}^{2},\\ B & = & 2+5{e}_{1}^{2}-7{e}_{1}^{2}\cos \left(2{g}_{1}\right),\end{array}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\cos \phi =-\cos {g}_{1}\cos {g}_{2}-\cos {i}_{\mathrm{mut}}\sin {g}_{1}\sin {g}_{2}.\end{eqnarray} \tag{ 13 }$$

The secular perturbation theory in this specific case is called the EKL mechanism, as previously mentioned. The time evolution for ω, e, and i_mut can be derived through Equation (8) as well. Here, _M can be used as an indicator to characterize the strength of the octupole-level effect on the system's perturbation potential.

The investigation of the amplitude of the inclination oscillation reveals that γ Cep Ab could evolve into an extremely inclined orbit. However, one may encounter difficulties in accurately predicting the theoretical correlation between the extent of the orbital inclination excitation and the initial conditions. Katz et al. (2011) defined a complex function of the analytical critical flip condition by averaging over the quadrupole-level effect.

We provide the planetary mass m₁ ∈ [5, 16.9] M_Jup in Section 2.2. Here, we assume the mass of γ Cep Ab to be m₁ = {5, 9, 11, 15} M_Jup, then calculate the ranges of the perturbation Hamiltonian ${ \mathcal H }$ and the total angular momentum G_tot for each value of m₁ with known parameters {m₀, m₂, a₁, a₂} and e₁ ∈ [0, 1], e₂ ∈ [0.35, 0.45], and i_mut ∈ [0°, 180°]. Furthermore, we set g₁ = 0° and g₂ = 0° in the initial conditions, since g₁ and g₂ can always go through either 0° or 180° over secular evolution.

The contour maps of the total angular momentum G_tot and the perturbation Hamiltonian ${ \mathcal H }$ are simultaneously plotted in the (e₁, i_mut) plane. In Figure 7, we show two examples of m₁ = {5, 15} M_Jup with e₂ ∈ [0.35, 0.45], which is around the current observation and will not vary significantly over the evolution. For other planetary masses, the structures of the ${ \mathcal H }$–G_tot contour maps have similar features with different ranges of G_tot and ${ \mathcal H }$. When m₁ = 5 M_Jup, G_tot ∈ [1.7050, 1.8025], ${ \mathcal H }\in [-1.05\times {10}^{-6}$, 3.0 × 10⁻⁷], and when m₁ = 9 M_Jup, G_tot ∈ [1.6994, 1.8137], ${ \mathcal H }\in [-1.75\times {10}^{-6}$, 5.0 × 10⁻⁷]. When m₁ = 11 M_Jup, G_tot ∈ [1.6975, 1.8187], ${ \mathcal H }\in [-2.10\times {10}^{-6}$, 6.0 × 10⁻⁷], and when m₁ = 15 M_Jup, G_tot ∈ [1.696, 1.824] and ${ \mathcal H }\in [-2.8\times {10}^{-6}$, 8.0 × 10⁻⁷]. The cross points of the G_tot and ${ \mathcal H }$ contours indicate all the possible initial conditions of i_mut and e₁ over the secular evolution. The range of values of G_tot and ${ \mathcal H }$ in Figure 7 will change with variational planetary masses, while the structure of the contours will remain similar.

Zoom In Zoom Out Reset image size

To further derive the evolution results of these general initial conditions, we uniformly choose specific values of ${ \mathcal H }$ and G_tot between the upper and lower limits in Figure 7. For a given ${ \mathcal H }$, based on the conserved Hamiltonian and the total angular momentum shown in Equations (5) and (8) of the octupole perturbation theory, we choose at least 40 cases of initial e₁ and i_mut. Considering the conservation of energy and the total angular momentum, we further explore the extreme value of the orbital inclination and orbital stability under these selected initial conditions in detail.

In the non–test particle approximation under the classical KL mechanism, the eccentricity and inclination oscillate regularly over a well-defined timescale t_quad (Antognini 2015):

$$\begin{eqnarray}&&{t}_{\mathrm{quad}}\sim \displaystyle \frac{16}{15}\displaystyle \frac{{a}_{2}^{3}{\left(1-{e}_{2}^{2}\right)}^{3/2}\sqrt{{m}_{0}+{m}_{1}}}{{a}_{1}^{3/2}{m}_{2}k}.\end{eqnarray} \tag{ 14 }$$

This relationship was derived under consideration of the equation of motion of ω, by integrating between the maximum and minimum eccentricities. Here, t_quad can be applied to estimate the timescale in the EKL scenario.

According to the values of m₀, m₁, m₂, e₂, a₁, and a₂, the quadrupole period t_quad of the γ Cep system is estimated to be ∼1000 yr, which is consistent with our numerical simulation results. Here, we investigate the secular evolution of γ Cep Ab by considering a diverse planetary mass and performing the simulation for 100 Myr (∼10⁵ t_quad) using the RKF7(8) integrator. The observed orbital inclinations I₁ and I₂ are required to calculate the constant total angular momentum and Hamiltonian. Hereafter, i₁ and i₂ denote the angles of the orbital plane relative to the invariable plane. However, the true orbital inclinations of the inner and outer orbits over the evolution are not well known.

Panels (a)–(d) of Figure 8 each plots three sets of typical outcomes, with ${{ \mathcal H }}_{1}$, ${{ \mathcal H }}_{2}$, and ${{ \mathcal H }}_{3}$. If m₁ = 5 M_Jup, both the prograde and retrograde orbits will flip when e₁ < 0.5 for ${{ \mathcal H }}_{1}=-8.06\times {10}^{-8}$, e₁ > 0.65 for ${{ \mathcal H }}_{2}=-4.95\times {10}^{-7}$, and e₁ > 0.8 for ${{ \mathcal H }}_{3}=-6.08\times {10}^{-7}$. We note that most intense orbital flips take place when the initial i_mut is close to 90°.

Zoom In Zoom Out Reset image size

From the estimated ${m}_{1}\sin {I}_{1}$ in Table 1, and the law of cosines for the angles of a spherical triangle in Equation (6), if m₁ = 5 M_Jup and the observed inclination of planet b is roughly 20°, the mutual inclination i_mut ∼ 100° is I₂ = 1193. As a consequence, it is clear that under the selected values of ${{ \mathcal H }}_{1}$, ${{ \mathcal H }}_{2}$, and ${{ \mathcal H }}_{3}$, such a high mutual inclination could always be achieved for original prograde orbits. In fact, as shown in Figure 8, the evolution of the mutual inclination of prograde and retrograde orbits is approximately symmetric about 90°, which can also be obtained from Equation (8).

Similarly, if m₁ = 9 M_Jup, the observed inclination of planet b is nearly 11°, thus the mutual inclination i_mut is estimated to be 108°. When ${{ \mathcal H }}_{1}=-1.45\times {10}^{-7}$, e₁ ∼ 0.5, ${{ \mathcal H }}_{2}=-6.35\times {10}^{-7}$, e₁ ∼ 0.65, and ${{ \mathcal H }}_{3}=-1.22\times {10}^{-6}$, e₁ ∼ 0.8, the amplitude of inclination for the orbital flip could be raised up to about 180°.

When m₁ is equal to 11 M_Jup, the calculated inclination of γ Cep Ab is about 9° and the targeted mutual inclination i_mut ∼ 110° can be achieved when ${{ \mathcal H }}_{1}=-1.98\times {10}^{-7}$, with e₁ ∼ 0.5, ${{ \mathcal H }}_{2}=-8.56\times {10}^{-7}$, with e₁ ∼ 0.7, and ${{ \mathcal H }}_{3}\,=-1.7\times {10}^{-6}$, with e₁ ∼ 0.85. Thus, the critical value of the initial eccentricity for the orbital flips of ${{ \mathcal H }}_{2}$ in Figure 8(c) is larger than the values for Figures 8(a), (b), and (d).

When m₁ is 15 M_Jup, the observed inclination of γ Cep Ab is about 6° and the critical conditions for the targeted mutual inclination i_mut ∼ 113° are ${{ \mathcal H }}_{1}=-2.4\times {10}^{-7}$, with e₁ ∼ 0.5, ${{ \mathcal H }}_{2}=-1.44\times {10}^{-6}$, with e₁ ∼ 0.65, and ${{ \mathcal H }}_{3}=-2.03\times {10}^{-6}$, with e₁ ∼ 0.8. In Figure 8(d), for ${{ \mathcal H }}_{2}=-1.44\times {10}^{-6}$, the orbital flip occurs when the initial i_mut < 60°. Thus, it is easier for γ Cep Ab to reach the targeted mutual inclination with m₁ = 15 M_Jup, e₁ < 0.7, and the critical initial i_mut being lower than 60°.

Above all, we conclude that the initial conditions for the orbital flips under investigation for the γ Cep system are e₁ > 0.5 and i_mut ∈ [60°, 120°], and that various planetary masses simply affect the critical eccentricity when e₁ > 0.6, with little influence on the maximum mutual inclination. The distribution tendency of the flipping conditions with various initial inclinations and eccentricities is similar to that in Lei (2022), while we perform the variation of the flipping conditions under different m₁, which affects the octupole-level factor _M .

To obtain an in-depth understanding of the rolling-over orbits reported in Section 4.1, we will further explore the orbital flip timescale over secular evolution, which may rely on the initial conditions, e.g., m₁, e₁, and i_mut. The duration of the flip is critical for the observational possibility of potentially extremely inclined S-type planets that are transforming between prograde and retrograde orbits. An approximate analytical timescale for the first flip for a nonchaotic orbit when m₁ → 0 has been given by Katz et al. (2011) and Antognini (2015):

$$\begin{eqnarray}&&t\sim \displaystyle \frac{128}{15\pi }\displaystyle \frac{{a}_{2}^{3}}{{a}_{1}^{3/2}}\displaystyle \frac{\sqrt{{m}_{0}}}{{{km}}_{2}}\sqrt{\displaystyle \frac{10}{\epsilon }}{\left(1-{e}_{2}\right)}^{3/2},\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&\epsilon \equiv \displaystyle \frac{{e}_{2}}{1-{e}_{2}^{2}}\displaystyle \frac{{a}_{1}}{{a}_{2}}.\end{eqnarray} \tag{ 16 }$$

This theoretical flip timescale is effective when the initial conditions meet e₁ → 0, ω₁ → 0, and i_mut → 90°. Note that the flip timescale for the circular test particle approximation is mainly determined by e₂ when the object masses and the orbital SMAs are fixed. Nevertheless, in the case of high-inclination oscillation, the timescale for the first flip is difficult to quantitatively assess, because this evolution is likely to be chaotic. In order to investigate the dependence of the orbital flip possibility on the integration timescale, we adopt the definition of the flip ratio f = t_flip/t_total (Teyssandier et al. 2013) in order to describe the timescale for the first flip and the observational possibility of the orbital flip process over secular evolution, where t_flip is the duration of the orbital flip process and t_total represents the total evolution timescale. We extract the cases from Section 4.1 in which orbital flips occur in Figure 9. For these rolling-over orbits, the integrations stop when e₁ becomes larger than 0.9999 and the planet falls into the Roche limit of the primary star.

Zoom In Zoom Out Reset image size

In Figure 9(a), we observe that when ${{ \mathcal H }}_{1}=-8.06\times {10}^{-8}$, the orbits turn over after 30 Myr and take more than 80 Myr to arrive at the equilibrium f. However, the flips under ${{ \mathcal H }}_{2}$ and ${{ \mathcal H }}_{3}$ occur from the very beginning of the evolution, and reach the equilibrium value before 1 Myr. This phenomenon also shows up in Figures 9(b) and (c), indicating that the flip timescale decreases with the increase of ${ \mathcal H }$, under most circumstances. The oscillation timescale of f over the flip procedure decreases with the rise of the perturbation Hamiltonian, as can be seen from the blue and green curves in the four subfigures.

In Figure 9(d), there is a special flipping case for ${{ \mathcal H }}_{2}$ when m₁ = 15 M_Jup. This rolling-over orbit has initial conditions of e₁ = 0.64 and i_mut = 589 in Figure 8(d). In comparison to the other cases under ${{ \mathcal H }}_{2}$ for different masses of planet, it can be concluded that a relatively low initial i_mut below 60° leads to a relatively large timescale for the first orbital flip, as well as more time to reach the equilibrium value of the flip ratio.

Moreover, we investigate the equilibrium value f for the flip cases for ${{ \mathcal H }}_{2}$ and ${{ \mathcal H }}_{3}$. We find that the final f of the dotted–dashed lines are all located above 0.5, while those of the solid lines are all lower than 0.5. Thus, the time duration of the flip process for those original retrograde orbits is larger. The maximum f occurs in Figure 9(b), with f ∼ 0.65 when m₁ = 9 M_Jup, ${{ \mathcal H }}_{2}=-6.35\times {10}^{-7}$, whereas the minimum f occurs in Figure 9(d), with f ∼ 0.4 when m₁ = 15 M_Jup, ${{ \mathcal H }}_{3}=-2.03\times {10}^{-6}$. Thus, the equilibrium value of the flip ratio under the EKL mechanism is related to the planetary mass, initial eccentricities, and mutual inclinations. When the equilibrium value of f is closer to 0.5, the observational possibility of the retrograde orbit transforming from the prograde is higher.

For a planetary system with given orbital parameters, the periodic orbits and stability are identified by the representative plane of (e₁, e₂), the level curves in the (e₁, g₁) plane, the Poincaré surface of section, and the long-term stability criterion. These methods are applied to our investigation of the stability of specific S-type planets in binary systems, while the perturbative treatment and the invariant manifolds characterize the stability with a fixed Hamiltonian (Lei 2022).

To derive a global and comprehensive view of the system dynamics, we first attempt to construct a parametrical analysis of the secular model. We first show the representative plane of (e₁, e₂) by Michtchenko & Malhotra (2004), which was then followed by studies in secular dynamics and resonances (Libert & Henrard 2006, 2007; Henrard & Libert 2008; Libert & Henrard 2008). This approach relies on the secular averaging of the full Hamiltonian and can be applied to construct global phase portraits of the systematic variables, thereby detecting resonances and discerning periodic orbits in the N-body dynamics.

According to Michtchenko & Malhotra (2004), the secular motion of planetary systems is mainly decided by the global quantities of total energy and the angular momentum deficit. The phase structures are determined by two constants: a₁/a₂ and m₁/m₂. The description of the secular behavior of this system is derived by the distribution of the initial values of e₁, e₂, g₁, g₂, and i_mut. To simplify the model, we fix m₁ = 15 M_Jup and i_mut = 60° in Figure 10, $\sqrt{1-{e}_{1}^{2}}$ cosg₁ > 0 with g₁ = 0°, while $\sqrt{1-{e}_{1}^{2}}$ cosg₁ < 0 with g₁ = 180°, and g₁ can always go through 0° and 180° over the evolution. As seen from Figure 10, we may come to the conclusion that periodic orbits occur when G_tot and ${ \mathcal H }$ have two cross points.

Zoom In Zoom Out Reset image size

For the secular evolution of the γ Cep Ab B system, one couple of free variables is e₁ and g₁ in the 2 dof averaged problem. Thus, we can first easily perform the qualitative analysis in the (e₁, g₁) plane (Tan et al. 2020). Given the limit of G_tot in Section 4.1, we select and fix the value of G_tot to present level curves of the Hamiltonian in the parameter space of (e₁, g₁), as shown in Figure 11. The green stream lines for circulating orbits and the blue circles for resonant orbits are easy to distinguish. With the changes of the perturbation Hamiltonian, the dynamical structure transitions from circulation to libration and the evolution interval of e₁ vary. The smallest circles have the largest magnitude of ${ \mathcal H }$, which indicates that the secular resonance effect is the most intense, as shown in Figure 12, with further details.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to the representative plane of ($\sqrt{1-{e}_{1}^{2}}$, $\sqrt{1-{e}_{2}^{2}}$) and the level curves in the (e₁, g₁) plane, we further investigate the stability of the rolling-over orbits under the selected Hamiltonian with the Poincaré surface of section. For a two-dimensional Hamiltonian system, the Poincaré surface of section is a commonly used method: under a given energy integral, the phase flow of the system is three-dimensional. By selecting a suitable section, the intersection point of the systematic phase flow and the section is projected in two dimensions. In other words, the Poincaré surface of section reflects the geometry of the system dynamics. Given the integration of energy (i.e., the octupole perturbation Hamiltonian) and the selection of the section (i.e., the g₁ = 0 section), the surface of section can be generated by plotting intersection points in the g₂ − e₂ frame for the plane of g₁ = 0.

The surface of section enables us to identify the order of orbital resonance and the dynamical stability of the system through a series of geometric structures. For a given fixed value of the Hamiltonian, the orbital mode can be derived from only two orbital parameters. Various orbital modes, including resonance, circulation, and chaos, are distinguished by the shape of the point distribution in the section.

We present typical structures in the surfaces of section of the γ Cep Ab B system in Figure 12. The surfaces of section are plotted in the g₁ = 0 plane for m₁ = {5, 9, 11, 15} M_Jup, with a corresponding perturbation Hamiltonian. To initially set e₂ and g₂ in the uniform grid, we solve other initial parameters, based on the conservation of the system's energy and the total angular momentum. There are two kinds of regions in these plots: the closed circles are resonant orbits (quasiperiodic orbits), where e₂ and g₂ oscillate in the bounded regions, while the stream lines denote circulative orbits, where at least one of the parameters circulates.

Figure 12 simply exhibits the dynamical structure of the surfaces of section for m₁ = 11 and 15 M_Jup. With increasing ${ \mathcal H }$, the structures in the surfaces of section appear to be diverse. As shown in Figures 12(a), (c), and (e), when m₁ = 11 M_Jup, the orbits with respect to g₂ = 180° disappear, whereas both orbits with fixed points of g₂ = 0° and 180° remain. These libration regions perfectly match the orbits that would turn over. In Figures 12(b), (d), and (f), when m₁ = 15 M_Jup, the positions of the fixed angles of g₂ are the same as for m₁ = 11 M_Jup. Figure 12(a) shows that the orbital flips of γ Cep Ab are dominated by the octupole-level resonance, and the oscillation amplitude of e₂ in these flipping cases is close to 0.01. In Figures 12(e) and (f), when m₁ = 11 M_Jup, ${ \mathcal H }=-1.7\times {10}^{-6}$, and m₁ = 15 M_Jup, ${ \mathcal H }=-2.03\times {10}^{-6}$, respectively, most of the orbits around the initial values of e₂ could turn over in the secular evolution. Additional simulations for m₁ = 5 and 9 M_Jup reveal similar structures, but with various oscillation intervals of e₂, when compared to above results.

The orbital flip criterion of e₂ can be evaluated with the conserved total angular momentum in Equation (5), when i₂ = 0° and i_mut = 90°:

$$\begin{eqnarray}&&{e}_{2,{flip}}=\sqrt{1-{\left(\displaystyle \frac{{G}_{\mathrm{tot}}}{{L}_{2}}\right)}^{2}}.\end{eqnarray} \tag{ 17 }$$

We mark this flip criterion in Figure 12 with the horizontal red lines. The circles and stream lines crossing the horizontal red line represent those orbits that invert regularly, while the ellipses above the horizontal red line in Figures 12(a)–(b) stand for the orbits that do not invert. Here, we find that the structures in the phase portrait become much clearer with the decrease in the perturbation Hamiltonian ${ \mathcal H }$, where the chaotic orbits disappear, similar to those in Lei (2022).

Additionally, we compare the boundaries of the regions of orbital flips in the (e_1,0, i_mut,0) space for the nonrestricted model with those of Figure 3 (Lei 2022). By analyzing the results in the Poincaré surfaces of section, the inner orbit appears to be more stable in our hierarchical system when the planetary eccentricity 0.5 < e_1,0 < 0.6, which is simply occupied by circulating and librating orbits.

In order to further confirm whether the orbital flip in γ Cep Ab is stable, we adopt the long-term stability criterion given by Mardling & Aarseth (2001):

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{2}}{{a}_{1}}\gt 2.8{\left(1+\displaystyle \frac{{m}_{2}}{{m}_{0}+{m}_{1}}\right)}^{2/5}\displaystyle \frac{{\left(1+{e}_{2}\right)}^{2/5}}{{\left(1-{e}_{2}\right)}^{6/5}}\left(1-\displaystyle \frac{0.3{i}_{\mathrm{mut}}}{180^\circ }\right),\end{eqnarray} \tag{ 18 }$$

where e₂ and i_mut are time-varying, and the constants a₁, a₂, m₀, m₁, and m₂ can be moved to the left side of the expression. We then derive the new expression of this criterion by a new variable S:

$$\begin{eqnarray}&&S=\displaystyle \frac{{\left(1+{e}_{2}\right)}^{2/5}}{{\left(1-{e}_{2}\right)}^{6/5}}\left(1-\displaystyle \frac{0.3{i}_{\mathrm{mut}}}{180^\circ }\right).\end{eqnarray} \tag{ 19 }$$

For m₁ ∈ [5, 15] M_Jup, a₁ = 2.14 au, a₂ = 18.62 au, thus S should meet the requirement S < 2.807 to maintain the stability of the system. The calculated maximum S for the rolling-over orbits in Figures 12(a), (c), and (f) is 2.13, which is definitely within the stability criterion. In Figures 12(b), (d), and (f), the calculated maximum S is 2.64, providing stable cases using Equation (19). Hence, the orbital flips in Figure 12 are stable, without the chaotic excitation of the binary's eccentricity and inclination.

Comparing with previous work on the stability of the γ Cep system, Haghighipour (2004) conducted an extensive numerical study of the orbital stability of the γ Cep system and suggested that the system could remain steady for i_mut ∈ [0°, 60°] and e₂ < 0.5. This condition is also supported in this work, since we constrain the eccentricity e₂ ∈ [0.35, 0.45], which is well consistent with our simulation results. Taking Figure 12(b) as an example, the critical value of the initial i_mut that is necessary for orbital flips to occur is lower than 60°, indicating a stable initial status, according to Haghighipour (2004), with regular oscillations of e₁ and i_mut causing the system always to return to a stable situation. Satyal et al. (2013) examined the stability and quasiperiodicity of γ Cep, and explored the orbital stability for various inclinations and binary eccentricities e₂ through the reliability comparison of chaos indicator. They demonstrated that the planet γ Cep Ab can remain stable for e₂ being as high as 0.6 or for i₁ ≤ 25°. In this work, for rolling-over orbits, we demonstrate that γ Cep Ab can also remain stable when e₂ ∼ 0.4 and i_mut ∼ 90°. We further confirm in this work that there is a great possibility of the flipping cases of γ Cep Ab being locked into a Kozai resonance, based on the surfaces of section in Figure 12 and the long-term stability criterion S.

From the simulation results, we also see that the planetary eccentricity could be stirred up to 0.9999, due to secular perturbation from the binary, as close approaches may eject the planet out of the system, so that it would not be observed or its orbit could be shrunk, owing to tides over the evolution. While this process takes too much time to be observed on the timescale of the EKL mechanisms, we employ i_mut to explore the stability of the system, as the evolutions of e₁ and i_mut are coupled under this scenario. The extreme oscillations of inclination and eccentricity would enhance the rate at which the system would be brought into an unstable situation. Li et al. (2014a) showed that there could be chaotic behavior when the mutual inclination between the inner and outer orbit remained high.

The instability of S-type planets may be induced by large eccentricity excitations. The planet may fall into the region of Roche lobe and merge into the primary, or move away from the inner binary system (Eggleton et al. 1998; Kiseleva et al. 1998; Ford et al. 2000).