STABLE CONIC-HELICAL ORBITS OF PLANETS AROUND BINARY STARS: ANALYTICAL RESULTS

In this section we start by considering a model system consisting of two immobile stars of masses μ and μ', and a planet of a unit mass moving around a circle in the plane perpendicular to the interstellar axis on which the circle is centered. The mass μ is at the origin and the Oz axis is directed to the mass μ' located at z = R. We introduce notations

$$\begin{eqnarray}&&Z=G\mu ,Z^{\prime} =G\mu ^{\prime} .\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant. The Hamilton function can be written in the cylindrical coordinates $(z,\rho ,\phi )$ as follows.

$$\begin{eqnarray}&&H=(p_{z}^{2}+p_{\rho }^{2}+p_{\varphi }^{2}/{{\rho }^{2}})/2+U(z,\rho ).\end{eqnarray} \tag{ 2 }$$

The last term in Equation (2) is the potential energy:

$$\begin{eqnarray}&&U(z,\rho )=-Z/{{({{z}^{2}}+{{\rho }^{2}})}^{1/2}}-Z^{\prime} /{{[{{(R-z)}^{2}}+{{\rho }^{2}}]}^{1/2}}.\end{eqnarray} \tag{ 3 }$$

The Hamiltonian equations of the motion relate the momenta and the velocities:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} dz/dt & = & \partial H/\partial {{p}_{z}}={{p}_{z}},\quad d\rho /dt=\partial H/\partial {{p}_{\rho }}={{p}_{\rho }}, \\ d\varphi /dt & = & \partial H/\partial {{p}_{\varphi }}={{p}_{\varphi }}/{{\rho }^{2}}. \\ \end{array}\end{eqnarray} \tag{ 4 }$$

Since $\varphi$ is the cyclic coordinate, the corresponding momentum is conserved:

$$\begin{eqnarray}&&{{p}_{\varphi }}={{\rho }^{2}}d\varphi /dt\equiv \ M={\rm const}.\end{eqnarray} \tag{ 5 }$$

Physically, the quantity M is a projection of the planet angular momentum on the interstellar axis. Thus, we can separate the z- and ρ-motions from the $\varphi$-motion. So after solving the z- and ρ-motions, we can find the $\varphi$-motion using Equation (5) where M is the separation constant.

The separated z- and ρ-motions can be described by the following Hamilton function.

$$\begin{eqnarray}&&H=(p_{z}^{2}+p_{\rho }^{2})/2+{{U}_{{\rm eff}}}(z,\rho ).\end{eqnarray} \tag{ 6 }$$

The last term in Equation (6) is an effective potential energy:

$$\begin{eqnarray}&&{{U}_{{\rm eff}}}(z,\rho )={{M}^{2}}/(2{{\rho }^{2}})+U(z,\rho ).\end{eqnarray} \tag{ 7 }$$

Now we define scaled (dimensionless) variables w and v, a scaled projection of the angular momentum m, as well as a ratio of the star masses b:

$$\begin{eqnarray}&&w\equiv z/R,\quad v\equiv \;\;\rho /R,\quad m\equiv M/{{(ZR)}^{1/2}},\quad b\equiv \;\;\mu ^{\prime} /\mu .\end{eqnarray} \tag{ 8 }$$

In the scaled notations, the effective potential energy can be represented in the form

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{U}_{{\rm eff}}} & \equiv & (Z/R){{u}_{{\rm eff}}}(w,v,m,b),\ \;{{u}_{{\rm eff}}}(w,v,m,b)={{m}^{2}}/(2{{v}^{2}}) \\ {} & {} & -{{({{w}^{2}}+{{v}^{2}})}^{-1/2}}-b{{[{{(1-w)}^{2}}+{{v}^{2}}]}^{-1/2}}. \\ \end{array}\end{eqnarray} \tag{ 9 }$$

In equilibrium, the derivatives of U_eff by both scaled coordinates $(w,v)$ must vanish. By setting to zero the partial derivative of U_eff by w:

$$\begin{eqnarray}&&\partial {{U}_{{\rm eff}}}/\partial w=(Z/R)\partial {{u}_{{\rm eff}}}/\partial w=0,\end{eqnarray} \tag{ 10 }$$

we obtain the relation:

$$\begin{eqnarray}&&b/{{[{{(1-w)}^{2}}+{{v}^{2}}]}^{3/2}}=w/[(1-w){{({{w}^{2}}+{{v}^{2}})}^{3/2}}].\end{eqnarray} \tag{ 11 }$$

This relation determines a line v₀(w) in the plane (w, v), which is the locus of the equilibrium points:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{v}_{0}}(w,b) & = & \{[{{w}^{2/3}}{{(1-w)}^{4/3}}-{{b}^{2/3}}{{w}^{2}}]/ \\ {} & {} & \left[ {{b}^{2/3}} \right.-\left. {{w}^{2/3}}/{{(1-w)}^{2/3}} \right]{{\}}^{1/2}}. \\ \end{array}\end{eqnarray} \tag{ 12 }$$

The equilibrium points exist in the following regions of w:

• A.  
  for $b\gt 1$, the regions are 0 $\leqslant$ w $\leqslant$ 1/$(1+{{b}^{1/2}})$ and b/(1+b) < w $\leqslant$ 1;
• B.  
  for $b\lt 1$, the regions are $0\leqslant w\lt b/(1+b)$ and 1/$(1+{{b}^{1/2}})\;\leqslant$ w $\leqslant$ 1;
• C.  
  for b = 1, the region is the entire range of 0 $\leqslant$ w $\leqslant$ 1.

In the rest of the paper we call these intervals the "allowed ranges" of w. As an example, Figure 1 shows the dependence v₀(w) for b = 3.

Zoom In Zoom Out Reset image size

Now we set to zero the partial derivative of U_eff by v:

$$\begin{eqnarray}&&\partial {{U}_{{\rm eff}}}/\partial v=(Z/R)\partial {{u}_{{\rm eff}}}/\partial v=0.\end{eqnarray} \tag{ 13 }$$

On substituting v₀(w,b) from Equation (12) instead of v, we obtain

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} m=v_{0}^{2}(w,b)/\{{{(1-w)}^{1/2}}\left[ {{w}^{2}} \right.+{{\left. v_{0}^{2}(w,b) \right]}^{3/4}}\} \\ \equiv {{m}_{0}}(w,b). \\ \end{array}\end{eqnarray} \tag{ 14 }$$

In the course of the derivation of Equation (14), we employed Equation (11) for removing the explicit dependence on b. As a result, the stellar mass ratio b shows up in Equation (14) only as the argument of the function v₀(w, b). In the following we will also use Equation (11) for this purpose in similar derivations without further notice.

For each pair (w, b), where w belongs to the allowed ranges, Equation (12) specifies the equilibrium value v₀(w, b), while Equation (14) specifies the equilibrium value m₀(w, b). As an example, Figure 2 shows the dependence m₀(w) for b = 3.

Zoom In Zoom Out Reset image size

Equivalently, for a pair (b, m), solving Equation (14) produced one or more values of w_i(b, m). Subsequently, Equation (12) specifies the values v_i = v₀(w_i(b, m), b). The result is pairs (w_i, v$_{i})$ corresponding to equilibria. So, for each ratio of the stellar masses b, there are pairs of equilibrium values (w_i, v$_{0i}$, m$_{0i})$. Here v$_{0i}\;\;\equiv$ v₀(w$_{i})$, m$_{0i}\;\;\equiv$ m₀(w$_{i})$, where w_i is taken from the allowed ranges.

Now for some value of b, we consider small deviations from equilibrium values (w_i, v$_{0i}$, m$_{0i})$:

$$\begin{eqnarray}&&\delta w\equiv w-{{w}_{i}},\delta v\equiv v-{{v}_{0i}}.\end{eqnarray} \tag{ 15 }$$

By expanding the scaled effective potential energy u_eff in terms of δw and δv, we get

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{u}_{{\rm eff}}} & \approx & {{u}_{0}}+{{u}_{ww}}{{(\delta w)}^{2}}/2+{{u}_{vv}}{{(\delta v)}^{2}}/2 \\ {} & {} & +{{u}_{wv}}(\delta w)(\delta v),{{u}_{0}}\;\;\equiv {{u}_{{\rm eff}}}({{w}_{i}},{{v}_{0i}},{{m}_{0i}}), \\ \end{array}\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{u}_{ww}} & \equiv & {{({{\partial }^{2}}{{u}_{{\rm eff}}}/\partial {{w}^{2}})}_{0}}=[1/(1-w)-3wP/{{Q}^{2}}]/{{({{w}^{2}}+v_{0}^{2})}^{3/2}}, \\ {{u}_{vv}} & \equiv & {{({{\partial }^{2}}{{u}_{{\rm eff}}}/\partial {{v}^{2}})}_{0}}=[1/(1-w)+3wP/{{Q}^{2}}]/{{({{w}^{2}}+v_{0}^{2})}^{3/2}}, \\ {{u}_{wv}} & \equiv & {{({{\partial }^{2}}{{u}_{{\rm eff}}}/\partial w\partial v)}_{0}}=3{{v}_{0}}w(2w-1)/ \\ {} & {} & [{{({{w}^{2}}+v_{0}^{2})}^{3/2}}{{Q}^{2}}]. \\ \end{array}\end{eqnarray} \tag{ 17 }$$

Here the suffix 0 at the derivatives means that after the differentiation one should set v = v₀(w) and m = m₀(w). (For brevity we dropped the suffix i.) The following notations are used in Equation (17).

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} P & \equiv & w(1-w)+v_{0}^{2}, \\ Q & \equiv & {{({{w}^{2}}+v_{0}^{2})}^{1/2}}{{[{{(1-w)}^{2}}+v_{0}^{2}]}^{1/2}}. \\ \end{array}\end{eqnarray} \tag{ 18 }$$

For transforming the effective potential energy u_eff to so-called normal coordinates, which would diagonalize the matrix composed of the second derivatives of u_eff in the general case where u$_{wv}\;\;\ne$ 0, we rotate the reference frame:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \delta w^{\prime} & = & \delta w{\rm cos} \alpha +\delta v{\rm sin} \alpha , \\ \delta v^{\prime} & = & -\delta w{\rm sin} \alpha +\delta v{\rm cos} \alpha . \\ \end{array}\end{eqnarray} \tag{ 19 }$$

The angle of the rotation, achieving the diagonalization, satisfies the condition

$$\begin{eqnarray}&&{\rm tan} 2\alpha =2{{u}_{wv}}/({{u}_{ww}}-{{u}_{vv}})=(1-2w){{v}_{0}}/P,\end{eqnarray} \tag{ 20 }$$

so that

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {\rm cos} \alpha & = & {{[(1+P/Q)/2]}^{1/2}}, \\ {\rm sin} \alpha & = & {{[(1-P/Q)/2]}^{1/2}}{\rm sign}(1-2w). \\ \end{array}\end{eqnarray} \tag{ 21 }$$

The effective potential energy can be written in the normal coordinates as follows.

$$\begin{eqnarray}&&{{u}_{{\rm eff}}}\;\;\approx {{u}_{0}}+\delta w{{^{\prime} }^{2}}\;\;\omega _{-}^{2}/2+\delta v{{^{\prime} }^{2}}\;\;\omega _{+}^{2}/2,\end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray}&&{{\omega }_{\pm }}\equiv {{[1/(1-w)\pm 3w/Q]}^{1/2}}/{{({{w}^{2}}+v_{0}^{2})}^{3/4}}.\end{eqnarray} \tag{ 23 }$$

It should be emphasized that ${{\omega }_{+}}$, which is the scaled frequency of small oscillations in the direction of the normal coordinate δv', is always real. Hereafter, any frequency F and its scaled counterpart f are interconnected by the following relation.

$$\begin{eqnarray}&&f\equiv {{({{R}^{3}}/Z)}^{1/2}}F.\end{eqnarray} \tag{ 24 }$$

Concerning the other eigenfrequency ${{\omega }_{-}}$, it is real under the condition

$$\begin{eqnarray}&&Q\geqslant 3w(1-w).\end{eqnarray} \tag{ 25 }$$

On substituting Q from Equation (18), the inequality (25) can be reformulated as follows.

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{v}_{0}}(w,b) & \geqslant & \{w(1-w)-1/2+\left[ 9{{w}^{2}}{{(1-w)}^{2}} \right. \\ {} & {} & -w(1-w)+1/4{{]}^{1/2}}{{\}}^{1/2}}\;\;\equiv {{v}_{{\rm crit}}}(w). \\ \end{array}\end{eqnarray} \tag{ 26 }$$

When the inequality (26) is met, the quantity ${{\omega }_{-}}$ is the frequency of small oscillations in the direction of the normal coordinate δw'.

We arrive to the conclusion that, if v₀(w, b) > v_crit(w), the effective potential energy has a two-dimensional minimum at the equilibrium values of w and v = v₀(w, b). In other words, if v₀(w, b) > v_crit(w), the equilibrium is stable. After defining a scaled (dimensionless) time of

$$\begin{eqnarray}&&\tau \equiv {{(Z/{{R}^{3}})}^{1/2}}t,\end{eqnarray} \tag{ 27 }$$

the resulting expression for the small oscillations around the stable equilibrium can be represented by the following formulas.

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \delta w(\tau ) & = & {{a}_{w}}[{\rm cos} ({{\omega }_{-}}\;\;\tau +{{\psi }_{w}})]{\rm cos} \alpha \\ {} & {} & -{{a}_{v}}[{\rm cos} ({{\omega }_{+}}\;\;\tau +{{\psi }_{v}})]{\rm sin} \alpha , \\ \delta v(\tau ) & = & {{a}_{w}}[{\rm cos} ({{\omega }_{-}}\;\;\tau +{{\psi }_{w}})]{\rm sin} \alpha \\ {} & {} & +{{a}_{v}}[{\rm cos} ({{\omega }_{+}}\;\;\tau +{{\psi }_{v}})]{\rm cos} \alpha . \\ \end{array}\end{eqnarray} \tag{ 28 }$$

Here sin α and cos α are determined by Equation (21), while amplitudes a_w, a_v and phases ${{\psi }_{w}}$, ${{\psi }_{v}}$—by initial conditions.

Now we rewrite Equation (5) in the scaled notations as

$$\begin{eqnarray}&&d\varphi /d\tau =m/{{v}^{2}}\end{eqnarray} \tag{ 29 }$$

and substitute in $v(\tau )\;\;\approx$ v$_{0}+\delta v(\tau )$, where the quantity $\delta v(\tau )$ is defined by Equation (29). Then we integrate it over τ and find the following solution for the $\varphi$-motion

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \varphi (\tau ) & \approx & {{f}_{p}}\tau -2\{\omega _{-}^{-1}{{a}_{w}}[{\rm sin} ({{\omega }_{-}}\tau +{{\psi }_{w}})-{\rm sin} {{\psi }_{w}}] \\ {} & {} & \times {\rm sin} \alpha +\omega _{+}^{-1}{{a}_{v}}[{\rm sin} ({{\omega }_{+}}\tau +{{\psi }_{v}})-{\rm sin} {{\psi }_{v}})] \\ {} & {} & \times {\rm cos} \alpha \}/\{{{(1-w)}^{1/2}}{{[{{w}^{2}}+v_{0}^{2}]}^{3/4}}{{v}_{0}}\}. \\ \end{array}\end{eqnarray} \tag{ 30 }$$

Here the quantity

$$\begin{eqnarray}&&{{f}_{p}}\;\;\equiv {{({{R}^{3}}/Z)}^{1/2}}\;\;{\Omega }=1/\{{{(1-w)}^{1/2}}{{[{{w}^{2}}+v_{0}^{2}]}^{3/4}}\}\end{eqnarray} \tag{ 31 }$$

is a scaled primary frequency of the $\varphi$-motion (the primary frequency of the $\varphi$-motion in the usual units of 1 s⁻¹ is denoted as ${\Omega }$).

Equations (30) and (31) demonstrate that the $\varphi$-motion is a rotation about the interstellar axis with the frequency f, slightly modulated by oscillations of the scaled radius of the orbit v at the frequencies ${{\omega }_{+}}$ and ω. Thus, for the stable motion, the planetary trajectory is a helix on the surface of a cone, with the axis coinciding with the interstellar axis. In this conic-helical state, the planet, while spiraling on the surface of the cone, oscillates between two end-circles that result from cutting the cone by two parallel planes perpendicular to its axis (Figure 3).

Zoom In Zoom Out Reset image size

The scaled eigenfrequencies ${{\omega }_{+}}$ and ${{\omega }_{-}}$ of the zρ-motion and the scaled frequency f_p of the $\varphi$-motion are related as follows.

$$\begin{eqnarray}&&(\omega _{+}^{2}+\omega _{-}^{2})/2=f_{p}^{2}.\end{eqnarray} \tag{ 32 }$$

Thus, the eigenfrequencies of the zρ-motion are generally of the same order of magnitude as the frequency of the $\varphi$-motion.

From this point on, the projection M of the planet orbital momentum on the interstellar axis will be fixed. We analyze the behavior of the energy at M = const > 0 (it is sufficient to consider only M > 0 because the results for M and -M are physically the same). In the analogous quantal problem of ORQ, this means studying "terms of the same symmetry."

Combining the definition of $m\equiv M/{{(ZR)}^{1/2}}$ from Equation (8) with (14), we get

$$\begin{eqnarray}&&R(w,b,M)={{M}^{2}}/[Zm_{0}^{2}(w,b)].\end{eqnarray} \tag{ 33 }$$

Now the coordinates z and $\rho$ can be represented as z(w, b, M) = wR(w, b, M) and $\rho$(w, b, M) = v$_{0}$(w, b)R(w, b, M). For any b > 0, for any w from the allowed ranges (controlled by the value of b), and for any M > 0, the functions R(w, b, M), z(w, b, M), and $\rho$(w, b, M) determine the equilibrium values of the interstellar distance R, the location of the orbital plane z, and the radius of the orbit $\rho$, respectively.

For simplicity, we neglect the small oscillations of the z- and $\rho$-coordinates and concentrate on the primary $\varphi$-motion, where the planet moves around the circle of the radius $\rho$(w, b, M) = v$_{0}$(w, b)R(w, b, M). Then the total energy E of the planet is the same as the effective potential energy presented by Equation (7):

$$\begin{eqnarray}&&E={{M}^{2}}/(2{{\rho }^{2}})+U(z,\rho ).\end{eqnarray} \tag{ 34 }$$

In the spirit of notations in Equation (9), we define a scaled total energy e in the following way:

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} E & \equiv & (Z/R)e(w,{{v}_{0}}(w,b),{{m}_{0}}(w,b),b), \\ {} & {} & e(w,{{v}_{0}}(w,b),{{m}_{0}}(w,b),b) \\ {} & {} & =m_{0}^{2}/(2v_{0}^{2})-{{({{w}^{2}}+v_{0}^{2})}^{-1/2}} \\ {} & {} & -b{{[{{(1-w)}^{2}}+v_{0}^{2}]}^{-1/2}}. \\ \end{array}\end{eqnarray} \tag{ 35 }$$

By employing Equations (11) and (14), the scaled total energy can be represented as

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} e(w,b) & = & -[w(1-w)+v_{0}^{2}(w,b)/2]/ \\ {} & {} & \{(1-w)\times \left. {{[{{w}^{2}}+v_{0}^{2}(w,b)]}^{3/2}} \right\}. \\ \end{array}\end{eqnarray} \tag{ 36 }$$

By substituting R(w, b, M) from Equation (33) in (35), we obtain:

$$\begin{eqnarray}&&E(w,b,M)={{(Z/M)}^{2}}m_{0}^{2}(w,b)e(w,b).\end{eqnarray} \tag{ 37 }$$

Here we take a fresh look at the combination of Equations (33) and (37). For any b > 0 and M $\geqslant$ 0, they determine in a parametric form (w being the parameter) the dependence of the energy E on the interstellar distance R. Borrowing the terminology from the analogous quantal problem, this dependence is called classical energy terms.

The plot of the scaled energy (M/Z)²E versus the scaled interstellar distance (Z/M²)R is presented in Figure 4 for the ratio of the stellar masses b = 3. There are three classical energy terms corresponding to the same value of M and two of them cross.

Zoom In Zoom Out Reset image size

This example, given for the ratio of the stellar masses μ'/μ = 3, actually represents a typical situation. We found that for any pair of μ and $\mu ^{\prime} \ne \mu$ there are three classical energy terms of the same symmetry and the upper term always crosses the middle term. (For μ' = μ there is only one term in the corresponding plot and no crossing—as should be expected.)

The origin of these three classical energy terms can be determined for any ratio of the stellar masses $\mu ^{\prime} /\mu \ne 1$. At R $\to \;\;\infty$, the lower term corresponds to the energy of the planet around a single star of the mass ${{\mu }_{{\rm max} }}\equiv {\rm max} (\mu ^{\prime} ,\mu )$, (E $\to -{{(G{{\mu }_{{\rm max} }}/M)}^{2}}$/2), slightly perturbed by the star of the mass ${{\mu }_{{\rm min} }}\equiv$ min $(\mu ^{\prime} ,\mu )$. The planet is in the state corresponding to the zero projection A$_{0}$ of the Runge–Lenz vector [13] on the interstellar axis. At R $\to$ 0, the lower term corresponds to the energy of the planet around an effective single star of the mass $\mu ^{\prime} +\mu$, (E $\to -{{[G(\mu +\mu ^{\prime} )/M]}^{2}}$/2).

At R $\to \;\;\infty$, the middle term corresponds to the energy of the planet around a single star of the mass ${{\mu }_{{\rm min} }}\equiv {\rm min} (\mu ^{\prime} ,\mu )$, (E $\to$ −${{(G{{\mu }_{{\rm min} }}/M)}^{2}}$/2), slightly perturbed by the star of the mass ${{\mu }_{{\rm max} }}\equiv$ max $(\mu ^{\prime} ,\mu )$. The planet is in the state corresponding to A$_{0}$ = 0.

At R $\to \;\;\infty$, the upper term corresponds to a near-zero-energy state where the planet is almost free. In terms of the parameter w, this term at $R\to \infty$ corresponds to $w\to 1/(1+{{b}^{1/2}})\equiv {{w}_{0}}$.

At the crossing of the middle and upper energy terms, the eigenfrequency ${{\omega }_{-}}$, defined by Equation (23), vanishes. The middle term corresponds to stable equilibria of the three-dimensional planetary motion, while the upper term corresponds to unstable equilibria of the three-dimensional planetary motion. The value of w_c corresponding to the crossing can be found in any of the following three: E'$({{w}_{c}})$ = 0, M'$({{w}_{c}})$ = 0 (where the symbol ' stands for the derivative), or

$$\begin{eqnarray}&&{{v}_{0}}({{w}_{c}},b)={{v}_{{\rm crit}}}({{w}_{c}}),\end{eqnarray} \tag{ 38 }$$

where the critical value of the scaled radius of the planetary orbit v_crit(w) was defined by Equation (26). While solving any of those three equations, only the root w_c within the interval (0, 1) should be accepted; in other words, the values 0 and 1 (the endpoints) should be excluded.

The Kepler frequency ω of the two stars orbiting their barycenter is

$$\begin{eqnarray}&&\omega \;\;={{[G(\mu +\mu ^{\prime} )/{{R}^{3}}]}^{1/2}}.\end{eqnarray} \tag{ 39 }$$

This is also the frequency of oscillations of the interstellar distance in the case of eccentric stellar orbits. According to Equation (24), the scaled, dimensionless counterpart of the Kepler frequency is

$$\begin{eqnarray}&&{{f}_{s}}={{({{R}^{3}}/Z)}^{1/2}}\omega ={{(1+b)}^{1/2}},\end{eqnarray} \tag{ 40 }$$

where we used the previously introduced notations Z = G$\mu$ and b = $\mu$'/$\mu$. The ratio of the scaled primary frequency f$_{p}$ of the planetary motion (given by Equation (31)) to the scaled Kepler frequency f$_{s}$ of the stars is

$$\begin{eqnarray}&&\begin{array}{ccccccccccccccc} r(w,b) \\ =1/\{{{(1+b)}^{1/2}}{{(1-w)}^{1/2}}{{[{{w}^{2}}+{{v}_{0}}{{(w,b)}^{2}}]}^{3/4}}\}. \\ \end{array}\end{eqnarray} \tag{ 41 }$$

As an example, Figure 5 shows the dependence r(w) for the ratio of the stellar masses b = 3. It is seen that this ratio becomes sufficiently large if the projection of the planetary orbit on the interstellar axis is either close to the star of smaller mass ($w\ll 1$) or close to the star of larger mass ((1 − w) $\ll \;1$).

Zoom In Zoom Out Reset image size

We denote by ${{D}_{1}}$ the scaled critical distance between the projection of the planetary orbit on the interstellar axis and the star of the smaller mass, such that at $w\lt {{D}_{1}}$ we have the ratio of the two frequencies r > 10. Similarly we denote by ${{D}_{2}}$ the scaled critical distance between the projection of the planetary orbit on the interstellar axis and the star of the larger mass, such that at (1 − w) $\lt \;{{D}_{2}}$ we have the ratio of the two frequencies r > 10. Figures 6 and 7 show the dependence of ${{D}_{1}}$ and ${{D}_{2}}$, respectively, on the ratio of the stellar masses b.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is seen that as the ratio b of the stellar masses increases, the scaled critical distance ${{D}_{1}}$ also increases, whereas the scaled critical distance ${{D}_{2}}$ decreases. The most important thing to note here is that there are ranges of parameters where the separation into the rapid and slow subsystems is justified, and those ranges are clearly shown in Figures 6 and 7 for any ratio of the stellar masses.

In the reference frame rotating together with the stars with the Kepler frequency $\omega$, the Hamilton function of the planet acquires an additional term –${\boldsymbol{r}} \;\cdot \;[{\boldsymbol{p}} \;\times \;{\boldsymbol{ \omega }} ]$ and becomes (see, e.g., Landau & Lifshitz 1960)

$$\begin{eqnarray}&&H={{p}^{2}}/2+U+{{H}_{1}}.{{H}_{1}}=-{\boldsymbol{r}} \;\cdot \;({\boldsymbol{p}} \;\times \;{\boldsymbol{ \omega }} ),\end{eqnarray} \tag{ 42 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{p}} ={\bf v}+{\boldsymbol{ \omega }} \;\times \;{\boldsymbol{r}} .\end{eqnarray} \tag{ 43 }$$

This leads to an additional (Coriolis) force

$$\begin{eqnarray}&&{{{\boldsymbol{F}} }_{1}}=-\partial {{H}_{1}}/\partial {\boldsymbol{r}} ={\boldsymbol{p}} \;\times \;{\boldsymbol{ \omega }} ={\boldsymbol{v}} \;\times \;{\boldsymbol{ \omega }} +{\boldsymbol{ \omega }} \;\times \;({\boldsymbol{ \omega }} \;\times \;{\boldsymbol{r}} ).\end{eqnarray} \tag{ 44 }$$

Because v $\sim \;\;{\Omega }\rho$—where ${\Omega }$ is the primary frequency of the planetary motion and $\rho$ is the average radius of the planetary orbit—and ${\Omega }\gg \omega$, then the additional force is approximately

$$\begin{eqnarray}&&{{{\boldsymbol{F}} }_{1}}\;\;\approx {\boldsymbol{v}} \;\times \;{\boldsymbol{ \omega }} .\end{eqnarray} \tag{ 45 }$$

Expression (45) has a clear physical meaning for the quantal counterpart-problem of ORQ. Namely, it is a Lorentz electric field ${\boldsymbol{v}} \times {{{\boldsymbol{B}} }_{{\rm eff}}}$/c in the effective magnetic field ${{{\boldsymbol{B}} }_{{\rm eff}}}=c{\boldsymbol{ \omega }}$ in atomic units (or ${{{\boldsymbol{B}} }_{{\rm eff}}}={{m}_{e}}c{\boldsymbol{ \omega }}$ in the CGS units, m$_{e}$ being the electron mass).

We choose the Ox axis along vector ${\boldsymbol{ \omega }}$, which is obviously perpendicular to the interstellar axis chosen as the Oz axis. Because the planetary velocity ${\boldsymbol{v}}$ is primarily in the xy-plane perpendicular to the interstellar axis, then the additional force ${{{\boldsymbol{F}} }_{1}}$ is primarily along the interstellar axis.

A more accurate calculation is as follows. Representing ${\boldsymbol{r}} ={{{\boldsymbol{e}} }_{x}}{\rm cos} {\Omega }t+{{{\boldsymbol{e}} }_{y}}{\rm sin} {\Omega }t$, we get ${\boldsymbol{v}} =(-{\Omega }{\rm sin} {\Omega }t){{{\boldsymbol{e}} }_{x}}$ + $({\Omega }{\rm cos} {\Omega }t){{{\boldsymbol{e}} }_{y}}$. Because ${\boldsymbol{ \omega }} =\omega {{{\boldsymbol{e}} }_{x}}$, we obtain from Equation (44):

$$\begin{eqnarray}&&{{{\boldsymbol{F}} }_{1}}=(-\omega {\Omega }{\rm cos} {\Omega }t){{{\boldsymbol{e}} }_{z}}+({{\omega }^{2}}{\rm sin} {\Omega }t){{{\boldsymbol{e}} }_{y}}.\end{eqnarray} \tag{ 46 }$$

In the ranges of parameters where ${\Omega }\gg \omega$, Equation (46) becomes

$$\begin{eqnarray}&&{{{\boldsymbol{F}} }_{1}}\;\;\approx (-\omega {\Omega }{\rm cos} {\Omega }t){{{\boldsymbol{e}} }_{z}}.\end{eqnarray} \tag{ 47 }$$

Thus, for the z$\rho$-motion (which in the scaled coordinates is wv-motion), the situation represents a two-dimensional oscillator, having the scaled eigenfrequencies ${{\omega }_{+}}$ and ${{\omega }_{-}}$ defined by Equation (23), which is driven by the force ${{{\boldsymbol{F}} }_{1}}$ oscillating at the frequency ${\Omega }$. Using the well-known formulas for driven oscillators (see, e.g., Jose & Saletan 1998), the solution in the coordinates w', v' rotated by the angle $\alpha$ (defined by Equation (20)) compared to the coordinates w, v (i.e., in the coordinates, where the two oscillators are decoupled), can be written as follows.

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \delta w^{\prime} (\tau ) & = & [({{\omega }_{s}}{{f}_{p}}{\rm cos} \alpha )/(f_{p}^{2}-\omega _{-}^{2})]{\rm cos} ({{f}_{p}}\tau ), \\ \delta v^{\prime} (\tau ) & = & [({{\omega }_{s}}{{f}_{p}}{\rm sin} \alpha )/(f_{p}^{2}-\omega _{+}^{2})]{\rm cos} ({{f}_{p}}\tau ). \\ \end{array}\end{eqnarray} \tag{ 48 }$$

Here f$_{p}$ is the scaled primary frequency of the planetary motion defined by Equation (31), $\tau$ is the scaled time defined by Equation (27), and ${{\omega }_{s}}$ is the scaled Kepler frequency of the stars' rotation:

$$\begin{eqnarray}&&{{\omega }_{s}}={{({{R}^{3}}/Z)}^{1/2}}\;\;\omega .\end{eqnarray} \tag{ 49 }$$

In the original scaled coordinates w, v, the forced oscillations are

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \delta w(\tau ) & = & {{\omega }_{s}}{{f}_{p}}\left[ ({{{\rm cos} }^{2}}\alpha )/(f_{p}^{2}-\omega _{-}^{2}) \right. \\ {} & {} & +\left. ({{{\rm sin} }^{2}}\alpha )/(f_{p}^{2}-\omega _{+}^{2}) \right]{\rm cos} ({{f}_{p}}\tau ), \\ \delta v(\tau ) & = & {{\omega }_{s}}{{f}_{p}}({\rm sin} \alpha {\rm cos} \alpha )\left[ 1/(f_{p}^{2}-\omega _{-}^{2}) \right. \\ {} & {} & -\left. 1/(f_{p}^{2}-\omega _{+}^{2}) \right]{\rm cos} ({{f}_{p}}\tau ). \\ \end{array}\end{eqnarray} \tag{ 50 }$$

Using the relation (32), the latter formulas can be rewritten as follows.

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} \delta w(\tau ) & = & [2{{\omega }_{s}}{{f}_{p}}/(\omega _{+}^{2}-\omega _{-}^{2})]{\rm cos} 2\alpha {\rm cos} ({{f}_{p}}\tau ), \\ \delta v(\tau ) & = & [2{{\omega }_{s}}{{f}_{p}}/(\omega _{+}^{2}-\omega _{-}^{2})]{\rm sin} 2\alpha {\rm cos} ({{f}_{p}}\tau ). \\ \end{array}\end{eqnarray} \tag{ 51 }$$

In the ranges of parameters where ${\Omega }\gg \omega$, and consequently f$_{p}\gg {{\omega }_{s}}$, and given that ${{\omega }_{+}},{{\omega }_{-}}$  and f$_{p}$ are generally of the same order of magnitude (according to Equation (32)), it is seen from Equation (50) that the forced oscillation of the planetary orbit occurs with a relatively small amplitude. Thus, the rotation of the stars does not affect the stability of the planetary orbit.

In the reference frame rotating together with the stars with the Kepler frequency $\omega$, a non-zero eccentricity $\varepsilon$ of the stars' orbits results in the oscillation of the interstellar distance R with the frequency $\omega$. In the ranges of parameters where ${\Omega }\gg \omega$, the oscillation of the interstellar distance is an adiabatic perturbation of the planetary motion. According to the principle of adiabatic invariance, the planetary motion will adjust to the slowly varying R while keeping as the constant the projection of the planetary angular momentum M on the interstellar axis (M is the adiabatic invariant).

The classical energy terms, such as those presented in Figure 4, correspond to the situation where M is a constant. Consequently, as a result of the eccentricity of the stars' orbits, the state of the planetary motion will adiabatically oscillate along a particular energy term E(R).

The entire lower energy term, which relates to the situation where the center of symmetry of the conic-helical planetary orbit is relatively close to the star of the larger mass, corresponds to a stable planetary motion. So, the eccentricity of the stars' orbits cannot affect the stability of the planetary motion in this case.

The middle energy term, which relates to the situation where the center of symmetry of the conic-helical planetary orbit is relatively close to the star of the smaller mass, is also the zone of stability as long as R > R_c, where R_c is the interstellar distance where the middle energy term crosses the upper energy term. For the state of the planetary motion corresponding to the middle energy term, the planetary motion would remain stable if, in the course of the oscillation along the middle energy term, the system could not reach R = R_c. However, if in the course of the oscillation along the middle energy term the system did reach R = R_c, the planetary motion would switch to the upper energy term corresponding to the unstable mode and the planet would be able to escape the system.

The stability of the conic-helical planetary orbit close to the star of the smaller mass can be also understood and analyzed quantitatively using Figure 2. This Figure depicts the dependence of the scaled projection of the planet angular momentum m ≡ M/(ZR)$^{1/2}$ on the interstellar axis versus the scaled averaged coordinate w ≡ z/R of the planet orbital plane along the interstellar axis. The left branch of this dependence reaches its maximum value m$_{{\rm max} }$ at w = w_c, corresponding to the crossing of the middle and upper energy terms in Figure 4. Values of m at w < w_c correspond to the stable motion. Under the condition that the orbital frequency of the planet is much greater than the Kepler frequency of the stars, at the average interstellar distance $\langle R\rangle$, the values of w and m are such that w < ${{D}_{1}}\ll 1$ and m < $m({{D}_{1}})\ll 1$ (${{D}_{1}}$ was defined before Figure 6 and in the caption to Figure 6). The point with coordinates w = ${{D}_{1}}$ and m = m(${{D}_{1}})$ in Figure 2 is very close to the origin. As R decreases to its minimum value R$_{{\rm min} }$, while M and Z remain constant, the value of m increases ${{[\langle R\rangle /{{R}_{{\rm min} }}]}^{1/2}}$ times: the corresponding point in Figure 2 moves up along the ascending part of the left branch. If m$_{{\rm max} }$/$m({{D}_{1}})\gt {{[\langle R\rangle /{{R}_{{\rm min} }}]}^{1/2}}$, then even at R = R$_{{\rm min} ,}$ the planetary orbit would remain stable.

Figure 8 presents the dependence of the ratio m$_{{\rm max} }$/$m({{D}_{1}})$ on the decimal logarithm of the ratio of the stellar masses b. It is seen that in the range 1 < b < 10, corresponding to the overwhelming majority of binary stars, this ratio is between 1.24 and 2.1.

Zoom In Zoom Out Reset image size

The average (over one revolution) interstellar distance is $\langle R\rangle$ = p/${{(1-{{\varepsilon }^{2}})}^{1/2}}$, where $\varepsilon$ is the eccentricity and p is the semi-latus rectum (which is the so-called harmonic average p = (1/R$_{{\rm min} }-1/{{R}_{{\rm max} }}$)/2). The ratio $\langle R\rangle$/R$_{{\rm min} }$, where R$_{{\rm min} }$ = p/(1+ $\varepsilon$), is

$$\begin{eqnarray}&&\langle R\rangle /{{R}_{{\rm min} }}={{(1+\varepsilon )}^{1/2}}/{{(1-\varepsilon )}^{1/2}}.\end{eqnarray} \tag{ 52 }$$

So, the above condition of the stability of the conic-helical planetary orbit close to the star of the smaller mass, m$_{{\rm max} }$/$m({{D}_{1}})\gt \ {{[\langle R\rangle /{{R}_{{\rm min} }}]}^{1/2}}$, can be reformulated as $\varepsilon \lt {{\varepsilon }_{c}}$(b) where the critical value of the eccentricity ${{\varepsilon }_{c}}$(b) is

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\varepsilon }_{c}} & = & \{{{[{{m}_{{\rm max} }}/m({{D}_{1}})]}^{4}}-1\}/ \\ {} & {} & \{{{[{{m}_{{\rm max} }}/m({{D}_{1}})]}^{4}}+1\}. \\ \end{array}\end{eqnarray} \tag{ 53 }$$

Figure 9 shows the dependence of the critical value of the eccentricity ${{\varepsilon }_{c}}$ on the decimal logarithm of the ratio of the stellar masses b. It is seen that in the range 1 < b < 10, corresponding to the overwhelming majority of binary stars, the critical value of the eccentricity is between 0.4 and 0.9. For most of those binary starts, the eccentricity is smaller than the critical value, so that the conic-helical planetary orbit would be stable even if it is close to the star of the smaller mass. (Remember that the conic-helical planetary orbit close to the star of the larger mass is stable for any eccentricity of the stars' orbit.)

Zoom In Zoom Out Reset image size

Let us consider, as an example, the Kepler-16 binary system. According to Doyle et al. (2011), the masses of the two stars are 0.6897 and 0.20255, respectively, in units of the mass of the Sun. So, the ratio of the star masses is b = 3.4. The semimajor axis of the stars' orbit is a = 0.22431 AU, while the eccentricity is $\varepsilon$ = 0.15944. So, the average interstellar separation is $\langle R\rangle$ = p/${{(1-\varepsilon )}^{1/2}}=a{{(1-\varepsilon )}^{1/2}}$ = 0.22144 AU.

If a planet in the conic-helical orbit was close to the star of the smaller mass, then the distance from the star of the smaller mass to the average orbital plane of the planet would be equal to or smaller than ${{D}_{1}}$(3.4)$\langle R\rangle$ = 0.0077 $\langle R\rangle$ = 0.0017 AU. The average radius of the planetary orbit would be equal to or smaller than v(0.0077)$\langle R\rangle$ = 0.029 AU. At b = 3.4, the critical value of the eccentricity ${{\varepsilon }_{c}}$(3.4) = 0.606. Because the actual eccentricity of the stars' orbit is 0.15944, then the conic-helical planetary orbit would be stable.

If a planet in the conic-helical orbit was close to the star of the larger mass, then the distance from the star of the larger mass to the average orbital plane of the planet would be equal to or smaller than ${{D}_{2}}$(3.4)$\langle R\rangle$ = 0.00216 $\langle R\rangle$ = 0.00048 AU. The average radius of the planetary orbit would be equal to or smaller than v(1 − 0.00216)$\langle R\rangle$ = 0.044 AU. The conic-helical planetary orbit would be stable because its stability close to the star of the larger mass is not affected by any value of the eccentricity of the stars' orbit.