# Stability interchanges in a curved Sitnikov problem

Luis Franco-Pérez<EMAIL_ADDRESS>,  Marian Gidea
<EMAIL_ADDRESS>,  Mark Levi<EMAIL_ADDRESS>and  Ernesto Pérez-
Chavela<EMAIL_ADDRESS>Departamento de Matemáticas Aplicadas y Sistemas,
UAM-Cuajimalpa, Av. Vasco de Quiroga 4871, México, D.F. 05348, México.
Department of Mathematical Sciences, Yeshiva University, 245 Lexington Ave,
New York, NY 10016, USA.  Mathematics Department, Penn State University,
University Park, PA 16802, USA.  Departamento de Matemáticas, ITAM México, Río
Hondo 1, Col. Progreso Tizapán, México D.F. 01080 .

###### Abstract.

We consider a curved Sitnikov problem, in which an infinitesimal particle
moves on a circle under the gravitational influence of two equal masses in
Keplerian motion within a plane perpendicular to that circle. There are two
equilibrium points, whose stability we are studying. We show that one of the
equilibrium points undergoes stability interchanges as the semi-major axis of
the Keplerian ellipses approaches the diameter of that circle. To derive this
result, we first formulate and prove a general theorem on stability
interchanges, and then we apply it to our model. The motivation for our model
resides with the $n$-body problem in spaces of constant curvature.

###### Key words and phrases:

Stability interchanges, qualitative theory, Sitnikov problem

##  1\. Introduction

###  1.1. A curved Sitnikov problem

We consider the following curved Sitnikov problem: Two bodies of equal masses
(primaries) move, under mutual gravity, on Keplerian ellipses about their
center of mass. A third, massless particle is confined to a circle passing
through the center of mass of the primaries, denoted by $P_{0}$, and
perpendicular to the plane of motion of the primaries; the second intersection
point of the circle with that plane is denoted by $P_{1}$. We assume that the
massless particle moves under the gravitational influence of the primaries
without affecting them. The dynamics of the massless particle has two
equilibrium points, at $P_{0}$ and $P_{1}$. We focus on the local dynamics
near these two points, more precisely, on the dependence of the linear
stability of these points on the parameters of the problem.

When the Keplerian ellipses are not too large or too small, $P_{0}$ is a local
center and $P_{1}$ is a hyperbolic fixed point. When we increase the size of
the Keplerian ellipses, as the distance between $P_{1}$ and the closest
ellipse approaches zero, then $P_{1}$ undergoes stability interchanges. That
is, there exists a sequence of open, mutually disjoint intervals of values of
the semi-major axis of the Keplerian ellipses, such that, on each of these
intervals the linearized stability of $P_{1}$ is strongly stable, and each
complementary interval contains values where the linearized stability is not
strongly stable, i.e., it is either hyperbolic or parabolic. The length of
these intervals approaches zero when the semi-major axis of the Keplerian
ellipses approaches the diameter of the circle on which the massless particle
moves. This phenomenon is the main focus in the paper.

It is stated in [16] and suggested by numerical evidence [11, 12] that the
linearized stability of the point $P_{0}$ also undergoes stability
interchanges when the size of the binary is kept fixed and the eccentricity of
the Keplerian ellipses approaches $1$.

Stability interchanges of the type described above are ubiquitous in systems
of varying parameters; they appear, for example, in the classical Hill’s
equation and in the Mathieu equation [24]. To prove the occurrence of this
phenomenon in our curved Sitnikov problem, we first formulate a general result
on stability interchanges for a general class of simple mechanical systems.
More precisely, we consider the motion of two bodies — one massive and one
massless — which are confined to a pair of curves and move under Newtonian
gravity. We let the distance between the two curves be controlled by some
parameter $\lambda$. We assume that the position of the infinitesimal particle
that achieves the minimum distance between the curves is an equilibrium point.
We show that, in the case when the minimum distance between the two curves
approaches zero, corresponding to $\lambda\rightarrow 0$, there existence a
sequence of mutually disjoint open intervals $(\lambda_{2n-1},\lambda_{2n})$,
whose lengths approach zero as $\lambda\to 0$, such that whenever
$\lambda\in(\lambda_{2n-1},\lambda_{2n})$ the linearized stability of the
equilibrium point is strongly stable, and each complementary interval contains
values of $\lambda$ where the linearized stability is not strongly stable.
From this result we derive the above mentioned stability interchange result
for the curved Sitnikov problem.

The curved Sitnikov problem considered in this paper is an extension of the
classical Sitnikov problem described in Section 1.2 (also, see e.g., [18]).
When the radius of the circle approaches infinity, in the limit we obtain the
classical Sitnikov problem — the infinitesimal mass moves along the line
perpendicular to the plane of the primaries and passing through the center of
mass. The equilibrium point $P_{1}$ becomes the point at infinity and is of a
degenerate hyperbolic type. Thus, stability interchanges of $P_{1}$ represent
a new phenomenon that we encounter in the curved Sitnikov problem but not in
the classical one. Also in the last case, it is well known that for
$\varepsilon=0$ the classical Sitnikov problem is integrable. In the case of
the curved one, numerical evidence suggests that it is not (see Figure 2).

The motivation for considering the curved Sitnikov problem resides in the
$n$-body problem in spaces with constant curvature, and with models of
planetary motions in binary star systems, as discussed in Section 1.3.

###  1.2. Classical Sitnikov problem

We recall here the classical Sitnikov problem. Two bodies (primaries) of equal
masses $m_{1}=m_{2}=1$ move in a plane on Keplerian ellipses of eccentricity
$\varepsilon$ about their center of mass, and a third, massless particle moves
on a line perpendicular to the plane of the primaries and passing through
their center of mass. By choosing the plane of the primaries the $xy$-plane
and the line on which the massless particle moves the $z$-axis, the equations
of motion of the massless particle can be written, in appropriate units, as

(1) $\ddot{z}=-\frac{2z}{(z^{2}+r^{2}(t))^{3/2}},$

where $r(t)$ is the distance from the primaries to their center of mass given
by

(2) $r(t)=1-\varepsilon\cos u(t),$

where $u(t)$ is the eccentric anomaly in the Kepler problem. By normalizing
the time we can assume that the period of the primaries is $2\pi$, and

(3) $r(t)=(1-\varepsilon\cos t)+O(\varepsilon^{2}),$

for small $\varepsilon$.

When $\varepsilon=0$, i.e., the primaries move on a circular orbit and the
dynamics of the massless particle is described by a 1-degree of freedom
Hamiltonian and so is integrable. Depending on the energy level, one has the
following types of solutions: an equilibrium solution, when the particle rests
at the center of mass of the primaries; periodic solutions around the center
of mass; escape orbits, either parabolic, that reach infinity with zero
velocity, or hyperbolic, that reach infinity with positive velocity.

When $\varepsilon\in(0,1)$, the differential equation (1) is non-autonomous
and the system is non-integrable. Consider the case $\varepsilon\ll 1$. The
system also has bounded and unbounded orbits, as well as unbounded oscillatory
orbits and capture orbits (oscillatory orbits are those for which
$\limsup_{t\to\pm\infty}|z(t)|=+\infty$ and
$\liminf_{t\to\pm\infty}|z(t)|<+\infty$, and capture orbits are those for
which $\limsup_{t\to-\infty}|z(t)|=+\infty$ and
$\limsup_{t\to+\infty}|z(t)|<+\infty$). In his famous paper about the final
evolutions in the three body problem, Chazy introduced the term oscillatory
motions [5], although he did not find examples of these, leaving the question
of their existence open. Sitnikov’s model yielded the first example of
oscillatory motions [23]. There are many relevant works on this problem,
including [1, 2, 3, 17, 18, 22, 6, 9, 10].

The curved Sitnikov problem introduced in Section 1.1 is a modification of the
classical problem when the massless particle moves on a circle rather than a
line. Here we regard the circle as a very simple restricted model of a space
with constant curvature. In Subsection 1.3 we introduce and summarize some
aspects of this problem.

###  1.3. The $n$-body problem in spaces with constant curvature

The $n$-body problem on spaces with constant curvature is a natural extension
of the $n$-body problem in the Euclidean space; in either case the
gravitational law considered is Newtonian. The extension was first proposed
independently by the founders of hyperbolic geometry, Nikolai Lobachevsky and
János Bolyai. It was subsequently studied in the late 19th, early 20th
century, by Serret, Killing, Lipschitz, Liebmann, Schering, etc. Schrödinger
developed a quantum mechanical analogue of the Kepler problem on the two-
sphere in 1940. The interest in the problem was revived by Kozlov, Harin,
Borisov, Mamaev, Kilin, Shchepetilov, Vozmischeva, and others, in the 1990’s.
A more recent surge of interest was stimulated by the works on relative
equilibria in spaces with constant curvature (both positive and negative) by
Diacu, Pérez-Chavela, Santoprete, and others, starting in the 2010’s. See [7]
for a history of the problem and a comprehensive list of references.

A distinctive aspect of the $n$-body problem on curved spaces is that the lack
of (Galilean) translational invariance results in the lack of center-of-mass
and linear-momentum integrals. Hence, the study of the motion cannot be
reduced to a barycentric coordinate system.

As a consequence, the two-body problem on a sphere can no longer be reduced to
the corresponding problem of motion in a central potential field, as is the
case for the Kepler problem in the Euclidean space. As it turns out, the two-
body problem on the sphere is not integrable [25].

Studying the three-body problem on spaces with curvature is also challenging.
Perhaps the simplest model is the restricted three-body problem on a circle.
This was studied in [8]. First, they consider the motion of the two primaries
on the circle, which is integrable, collisions can be regularized, and all
orbits can be classified into three different classes (elliptic, hyperbolic,
parabolic). Then they consider the motion of the massless particle under the
gravity of the primaries, when one or both primaries are at a fixed position.
They obtain once again a complete classification of all orbits of the massless
particle.

In this paper we take the ideas from above one step further, by considering
the curved Sitnikov problem, with the massless particle moving on a circle
under the gravitational influence of two primaries that move on Keplerian
ellipses in a plane perpendicular to that circle. In the limit case, when the
primaries are identified with one point, that is when the primaries coalesce
into a single body, the Keplerian ellipses degenerate to a point, and the
limit problem coincides with the two-body problem on a circle described above.

While the motivation of this work is theoretical, there are possible
connections with the dynamics of planets in binary star systems. About 20
planets outside of the Solar System have been confirmed to orbit about binary
stars systems; since more than half of the main sequence stars have at least
one stellar companion, it is expected that a substantial fraction of planets
form in binary stars systems. The orbital dynamics of such planets can vary
widely, with some planets orbiting one star and some others orbiting both
stars. Some chaotic-like planetary orbits have also been observed, e.g. planet
Kepler-413b orbiting Kepler-413 A and Kepler-413 B in the constellation
Cygnus, which displays erratic precession. This planet’s orbit is tilted away
from the plane of binaries and deviating from Kepler’s laws. It is
hypothesized that this tilt may be due to the gravitational influence of a
third star nearby [13]. Of a related interest is the relativistic version of
the Sitnikov problem [14]. Thus, mathematical models like the one considered
in this paper could be helpful to understand possible types of planetary
orbits in binary stars systems.

To complete this introduction, the paper is organized as follows: In Section 2
we go deeper in the description of the curved Sitnikov problem, studying the
limit cases and its general properties. In Section 3 we present a general
result on stability interchanges. In Section 4 we show that the equilibrium
points in the curved Sitnikov problem present stability interchanges. Finally,
in order to have a self contained paper, we add an Appendix with general
results (without proofs) from Floquet theory.

##  2\. The curved Sitnikov problem

###  2.1. Description of the model

We consider two bodies with equal masses (primaries) moving under mutual
Newtonian gravity on identical elliptical orbits of eccentricity
$\varepsilon$, about their center of mass. For small values of $\varepsilon$,
the distance $r(t)$ from either primary to the center of mass of the binary is
given by

(4) $\begin{split}r_{\varepsilon}(t;r)=&r\rho(t;\varepsilon),\qquad r>0,\\\
\rho(t;\varepsilon)=&\left(1-\varepsilon\cos(u(t))\right)=\left(1-\varepsilon\cos(t)\right)+\mathcal{O}(\varepsilon^{2}),\end{split}$

where $u(t)$ is the eccentric anomaly, which satisfies Kepler’s equation
$u-\varepsilon\sin u=(2\pi/\tau)t$, where $\tau/2\pi$ is the mean motion111The
mean motion is the time-average angular velocity over an orbit. of the
primaries. The expansion in (4) is convergent for
$\varepsilon<\varepsilon_{c}=0.6627...$; see [21].

A massless particle is confined on a circle of radius $R$ passing through the
center of mass of the binary and perpendicular to the plane of its motion. We
assume that the only force acting on the infinitesimal particle is the
component along the circle of the resultant of the gravitational forces
exerted by the primaries. The motion of the primaries take place in the
$xy$-coordinate plane and the circle with radius $R$ is in the $yz$-coordinate
plane. See Figure 1.

Figure 1. The curved Sitnikov problem.

We place the center of mass at the point $(0,R,0)$ in the $xyz$-coordinate
system. The position of the primaries are determined by the functions

$\displaystyle\mathbf{x}_{1}(t)$ $\displaystyle=$
$\displaystyle(r_{\varepsilon}(t;r)\sin t,R+r_{\varepsilon}(t;r)\cos t,0),$
$\displaystyle\mathbf{x}_{2}(t)$ $\displaystyle=$
$\displaystyle(-r_{\varepsilon}(t;r)\sin t,R-r_{\varepsilon}(t;r)\cos t,0).$

Note that $t=0$ corresponds to the passage of the primaries through the
pericenter at $y=R\pm r(1-\varepsilon)$ and $t=\pi$ to the passage of the
primaries through the apocenter at $y=R\pm r(1+\varepsilon)$; both peri- and
apo-centers lie on the plane of the circle of radius $R$ in the $y$-axis.

The position of the infinitesimal particle is $\mathbf{x}(t)=(0,y(t),z(t))$
(taking into account the restriction of motion for the infinitesimal particle
to the circle $y^{2}+z^{2}=R^{2}$). We will derive the equations of motion by
computing the gravitational forces exerted by the primaries:

$\displaystyle\mathbf{F}_{\varepsilon}(y,t;R,r)$ $\displaystyle=$
$\displaystyle-\frac{\mathbf{x}-\mathbf{x}_{1}}{||\mathbf{x}-\mathbf{x}_{1}||^{3}}-\frac{\mathbf{x}-\mathbf{x}_{2}}{||\mathbf{x}-\mathbf{x}_{2}||^{3}}$
$\displaystyle=$ $\displaystyle-\frac{(-r_{\varepsilon}(t;r)\sin t,y-R-
r_{\varepsilon}(t;r)\cos t,z)}{||\mathbf{x}-\mathbf{x}_{1}||^{3}}$
$\displaystyle-\frac{(r_{\varepsilon}(t;r)\sin t,y-R+r_{\varepsilon}(t;r)\cos
t,z)}{||\mathbf{x}-\mathbf{x}_{2}||^{3}}\,,$

where the distance from the particle to each primary is

$\displaystyle||\mathbf{x}-\mathbf{x}_{1}||=\left[r^{2}_{\varepsilon}(t;r)+2R^{2}+2Rr_{\varepsilon}(t;r)\cos
t-2y(R+r_{\varepsilon}(t;r)\cos t)\right]^{1/2},$
$\displaystyle||\mathbf{x}-\mathbf{x}_{2}||=\left[r^{2}_{\varepsilon}(t;r)+2R^{2}-2Rr_{\varepsilon}(t;r)\cos
t-2y(R-r_{\varepsilon}(t;r)\cos t)\right]^{1/2}.$

We note that when $r(1+\varepsilon)=2R$ the elliptical orbit of the primary
with the apo-center at $y<R$ crosses the circle of radius $R$, hence
collisions between the primary and the infinitesimal mass are possible.
Therefore we will restrict to $r<\frac{2R}{1+\varepsilon}$; when
$\varepsilon=0$, this means $r<2R$.

We write (2.1) in polar coordinates, that is $y=R\cos q$, $z=R\sin q$, and we
obtain

$\displaystyle\mathbf{F}_{\varepsilon}(q,t;R,r)$ $\displaystyle=$
$\displaystyle-\frac{(-r_{\varepsilon}(t;r)\sin t,R\cos q-R-
r_{\varepsilon}(t;r)\cos t,R\sin q)}{||\mathbf{x}-\mathbf{x}_{1}||^{3}}$
$\displaystyle-\frac{(r_{\varepsilon}(t;r)\sin t,R\cos
q-R+r_{\varepsilon}(t;r)\cos t,R\sin q)}{||\mathbf{x}-\mathbf{x}_{2}||^{3}}.$

The origin $q=0$ corresponds to the point $(0,R,0)$ in the $xyz$-coordinate
system. Thus, the primaries move on elliptical orbits around this point.

Next we will retain the component along the circle of the resulting force
(2.1). That is, we will ignore the constraint force that confines the motion
of the particle to the circle, as this force acts perpendicularly to the
tangential component of the gravitational attraction force. The unit tangent
vector to the circle of radius $R$ at $(0,R\cos(q),R\sin(q))$ pointing in the
positive direction is given by $\mathbf{u}(q)=(0,-\sin(q),\cos(q))$. The
component of the force $\mathbf{F}_{\varepsilon}(q,t;R,r)$ along the circle is
computed as

(7)
$\mathbf{F}_{\varepsilon}(q,t;R,r)\cdot\mathbf{u}(q)=-\frac{(R+r_{\varepsilon}(t;r)\cos(t))\sin(q)}{||\mathbf{x}-\mathbf{x}_{1}||^{3}}-\frac{(R-r_{\varepsilon}(t;r)\cos(t))\sin(q)}{||\mathbf{x}-\mathbf{x}_{2}||^{3}}.$

The motion of the particle, as a Hamiltonian system of one-and-a-half degrees
of freedom, corresponds to

(8) $\displaystyle\dot{q}$ $\displaystyle=$ $\displaystyle p\,,$
$\displaystyle\dot{p}$ $\displaystyle=$ $\displaystyle
f_{\varepsilon}(q,t;R,r),$

where

$\displaystyle f_{\varepsilon}(q,t;R,r)$ $\displaystyle:=$
$\displaystyle\mathbf{F}_{\varepsilon}(q,t;R,r)\cdot\mathbf{u}(q)\,,$
$\displaystyle||\mathbf{x}-\mathbf{x}_{1}||$ $\displaystyle=$
$\displaystyle\left[r^{2}_{\varepsilon}(t;r)+2R(1-\cos
q)(R+r_{\varepsilon}(t;r)\cos t)\right]^{1/2}\,,$
$\displaystyle||\mathbf{x}-\mathbf{x}_{2}||$ $\displaystyle=$
$\displaystyle\left[r^{2}_{\varepsilon}(t;r)+2R(1-\cos
q)(R-r_{\varepsilon}(t;r)\cos t)\right]^{1/2}\,,$

Hence

(9) $H_{\varepsilon}(q,p,t;R,r)=\frac{p^{2}}{2}+V_{\varepsilon}(q,t;R,r),$

where the potential is given by

(10) $\displaystyle V_{\varepsilon}(q,t;R,r)$
$\displaystyle=-\frac{1}{R}\left(\frac{1}{||\mathbf{x}-\mathbf{x}_{1}||}+\frac{1}{||\mathbf{x}-\mathbf{x}_{2}||}\right).$

###  2.2. Limit cases

The curved Sitnikov problem can be viewed as a link between the classical
Sitnikov problem and the Kepler problem on the circle, mentioned in Section 1.

####  2.2.1. The limit $R\rightarrow\infty$.

We express (7) in terms of the arc length $w=Rq$, obtaining

(11) ${f}_{\varepsilon}(w,t;R,r)=-\frac{(R+r_{\varepsilon}(t;r)\cos
t)\sin\left(w/R\right)}{\left[r^{2}_{\varepsilon}(t;r)+2R(1-\cos\left(w/R\right))(R+r_{\varepsilon}(t;r)\cos
t)\right]^{\frac{3}{2}}}\\\ -\frac{(R-r_{\varepsilon}(t;r)\cos
t)\sin\left(w/R\right)}{\left[r^{2}_{\varepsilon}(t;r)+2R(1-\cos\left(w/R\right))(R-r_{\varepsilon}(t;r)\cos
t)\right]^{\frac{3}{2}}}$

which we can write in a suitable form as

${f}_{\varepsilon}(w,t;R,r)=-\frac{w\frac{\sin\left(w/R\right)}{\left(w/R\right)}\left(1+\frac{r_{\varepsilon}(t;r)}{R}\cos
t\right)}{\left[r^{2}_{\varepsilon}(t;r)+2w^{2}\frac{(1-\cos\left(w/R\right))}{\left(w/R\right)^{2}}\left(1+\frac{r_{\varepsilon}(t;r)}{R}\cos
t\right)\right]^{\frac{3}{2}}}\\\
-\frac{w\frac{\sin\left(w/R\right)}{\left(w/R\right)}\left(1-\frac{r_{\varepsilon}(t;r)}{R}\cos
t\right)}{\left[r^{2}_{\varepsilon}(t;r)+2w^{2}\frac{(1-\cos\left(w/R\right))}{\left(w/R\right)^{2}}\left(1-\frac{r_{\varepsilon}(t;r)}{R}\cos
t\right)\right]^{\frac{3}{2}}}.$

Letting $R$ tend to infinity we obtain

$\lim_{R\rightarrow\infty}{f}_{\varepsilon}(w,t;R,r)=-\frac{2w}{\left(r^{2}_{\varepsilon}(t;r)+w^{2}\right)^{3/2}}\,,$

which is the classical Sitnikov Problem.

####  2.2.2. The limit $r\rightarrow 0$.

When we take the limit $r\rightarrow 0$ in (7) we are fusing the primaries
into a large mass at the center of mass and we obtain a two-body problem on
the circle. The component force along the circle corresponds to

(12) $\lim_{r\rightarrow
0}f_{\varepsilon}(q,t;R,r)=-\frac{\sin(q)}{\sqrt{2}R^{2}(1-\cos(q))^{3/2}}\,.$

This problem was studied in [8] with a different force given by

(13) $-\frac{1}{Rq^{2}}+\frac{1}{R(2\pi-q)^{2}},$

the distance between the large mass and the particle is measured by the arc
length (in that paper the authors assume that $R=1$). The potential of the
force (12) is

$V_{1}(q_{1})=-\frac{1}{R^{2}\sqrt{2}(1-\cos(q_{1}))^{1/2}},$

where $q_{1}$ denotes the angular coordinate, and the potential for (13) is

$V_{2}(q_{2})=-\frac{1}{Rq_{2}}-\frac{1}{R(2\pi-q_{2})},$

where $q_{2}$ denotes the angular coordinate. Each problem defines an
autonomous system with Hamiltonian

(14) $H_{i}(p_{i},q_{i})=\frac{1}{2}p_{i}^{2}+V_{i}(q_{i})\,,$

taking $p_{1}=dq_{1}/dt$, $p_{2}=dq_{2}/dt$ and $i=1,2$. Let $\phi^{i}_{t}$ be
the flow of the Hamiltonian $H_{i}$, and let $A_{i}$ denote the phase space,
$i=1,\ 2$

Using that all orbits are determined by the energy relations given by (14), it
is not difficult to define a homeomorphism $g:A_{1}\to A_{2}$ which maps
orbits of system (12) into orbits of system (13). In the same way we can
define a homeomorphism $h:A_{2}\to A_{1}$ which in fact is $g^{-1}$. This
shows the $C^{0}$ equivalence of the respective flows. One can show that the
two corresponding flows are $C^{0}$–equivalent.

We recall from [8] that the solutions of the two-body problem on the circle
(apart from the equilibrium antipodal to the fixed body) are classified in
three families (elliptic, parabolic and hyperbolic solutions) according to
their energy level. The elliptic solutions come out of a collision, stop
instantaneously, and reverse their path back to the collision with the fixed
body. The parabolic solutions come out of a collision and approach the
equilibrium as $t\rightarrow\infty$. Hyperbolic motions comes out of a
collision with the fixed body, traverse the whole circle and return to a
collision.

We remark that the two limit cases $R\to\infty$ and $r\to 0$ are not
equivalent. Indeed, in the case $r\to 0$ the resulting system is autonomous,
the point $q=0$ is a singularity for the system, and the point $q=\pi$ is a
hyperbolic fixed point, while in the case $R\to\infty$ the resulting system is
non-autonomous (for $\varepsilon\neq 0$), the point $q=0$ is a fixed point of
elliptic type, and the point $q=\infty$ is a degenerate hyperbolic periodic
orbit.

###  2.3. General properties

####  2.3.1. Extended phase space, symmetries, and equilibrium points

It is clear that, besides the limit cases $R\to\infty$ and $r\to 0$, the
dynamics of the system depends only on the ratio $r/R$, so we can fix $R=1$
and study the dependence of the global dynamics on $r$ where $0<r<2$. In this
case using (8) and (4) we get

(15)
${f}_{\varepsilon}(q,t;r)=-\frac{(1+r\rho(t;\varepsilon)\cos(t))\sin(q)}{\left[r^{2}\rho^{2}(t;\varepsilon)+2(1-\cos
q)(1+r\rho(t;\varepsilon)\cos t)\right]^{3/2}}\\\
-\frac{(1-r\rho(t;\varepsilon)\cos(t))\sin(q)}{\left[r^{2}\rho^{2}(t;\varepsilon)+2(1-\cos
q)(1-r\rho(t;\varepsilon)\cos t)\right]^{3/2}}.$

To study the non-autonomous system (8) we will make the system autonomous by
introducing the time as an extra dependent variable

(16)
$\mathcal{X}_{\varepsilon}(q,p,s;r)=\left\\{\begin{array}[]{rcl}\dot{q}&=&p\\\
\dot{p}&=&{f}_{\varepsilon}(q,s;r)\\\ \dot{s}&=&1\end{array}\right..$

This vector field is defined on $[0,2\pi]\times\mathbb{R}\times[0,2\pi]$,
where we identify the boundary points of the closed intervals. The flow of
$\mathcal{X}_{\varepsilon}$ possesses symmetries defined by the functions

$\displaystyle\mathbb{S}_{1}(q,p,s)$ $\displaystyle=$
$\displaystyle(-q,-p,s)\,,$ $\displaystyle\mathbb{S}_{2}(q,p,s)$
$\displaystyle=$ $\displaystyle(q,p,s+2\pi)\,,$
$\displaystyle\mathbb{S}_{3}(q,p,s)$ $\displaystyle=$
$\displaystyle(q+2\pi,p,s)\,,$ $\displaystyle\mathbb{S}_{4}(q,p,s)$
$\displaystyle=$ $\displaystyle(q,-p,-s)\,$

in the sense that

  * •

$\mathbb{S}_{1}(\mathcal{X}_{\varepsilon}(q,p,s))=\mathcal{X}_{\varepsilon}(\mathbb{S}_{1}(q,p,s))$,

  * •

$\mathcal{X}_{\varepsilon}(q,p,s)=\mathcal{X}_{\varepsilon}(\mathbb{S}_{2}(q,p,s))$,

  * •

$\mathcal{X}_{\varepsilon}(q,p,s)=\mathcal{X}_{\varepsilon}(\mathbb{S}_{3}(q,p,s))$

  * •

$\mathbb{S}_{4}(\mathcal{X}_{\varepsilon}(q,p,s))=-\mathcal{X}_{\varepsilon}(\mathbb{S}_{4}(q,p,s))$,

as can be verified by a direct computation.

The function $\mathbb{S}_{1}$ describes the symmetry respect to the trajectory
$(0,0,s)$, $\mathbb{S}_{2}$ and $\mathbb{S}_{3}$ describe the bi-periodicity
of ${f}_{\varepsilon}(q,s;r)$ and $\mathbb{S}_{4}$ describes the reversibility
of the system.

System (8) has two equilibria $(0,0)$ and $(\pi,0)$, which correspond to
periodic orbits for $\mathcal{X}_{\varepsilon}$.

While the classical Sitnikov equation is autonomous for $\varepsilon=0$, our
equation (16) is not, and thus we expect it to be non–integrable, as is borne
out by numerical simulation. Figure 2 shows a Poincaré section corresponding
to $s=0$ (mod $2\pi$), with $\varepsilon=0$ and $r=1$. This simulation
suggests that the invariant KAM circles coexist with chaotic regions.

Figure 2. Poincaré section for the curved Sitnikov problem, for
$\varepsilon=0$ and $r=1$.

In the sequel, we will analyze the dynamics around the equilibrium points
$(\pi,0)$ and $(0,0)$. One important phenomenon that we will observe is that
both equilibrium points undergo stability interchanges as parameters are
varied. More precisely, when $\varepsilon$ sufficiently small is kept fixed
and $r\rightarrow\frac{2R}{1+\varepsilon}$, the point $(\pi,0)$ undergoes
infinitely many changes in stability, and when $r$ is kept fixed and
$\varepsilon\to 1$, the point $(0,0)$ undergoes infinitely many changes in
stability.

In the next section we will first prove a general result.

##  3\. A general result on stability interchanges

In this section we switch to a more general mechanical model which exhibits
stability interchanges. We consider the motion under mutual gravity of an
infinitesimal particle and a heavy mass each constrained to its own curve and
moving under gravitational attraction, and study the linear stability of the
equilibrium point corresponding to the closest position between the particles
along the curves they are moving on. In Section 4, we will apply this general
result to the equilibrium point $P_{1}$ of the curved Sitnikov problem
described in Section 1.1. The fact that in the curved Sitnikov problem there
are two heavy masses, rather than a single one as considered in this section,
does not change the validity of the stability interchanges result, since, as
we shall see, what it ultimately matters is the time-periodic gravitational
potential acting on the infinitesimal particle.

To describe the setting of this section, consider a particle constrained to a
curve ${\bf x}={\bf x}(s,\lambda)$ in ${\mathbb{R}}^{3}$, where $s$ is the arc
length along the curve and $\lambda$ is a parameter with values in some
interval $[0,\lambda_{0}]$, $\lambda_{0}>0$. Another (much larger)
gravitational mass undergoes a prescribed  periodic motion according to ${\bf
y}={\bf y}(t,\lambda)={\bf y}(t+1,\lambda)$; see Figure 3. We assume that the
mass of the particle at ${\bf x}(s)$ is negligible compared to the mass at
${\bf y}(t)$, treating the particle at ${\bf x}(s)$ as massless.

Figure 3.  An infinitesimal particle with coordinates ${\bf x}(s)$ is
constrained to a curve and moves under the influence of the larger mass with
coordinates ${\bf y}(t)$.

To write the equation of motion for the unknown coordinate $s$ of the massless
particle, let:

(17) ${\bf z}(s,t,\lambda)={\bf x}(s,\lambda)-{\bf y}(t,\lambda).$

Assume that $s=0,\ t=0$ minimize the distance between the two curves:

(18) $|{\bf z}(0,0,\lambda)|=\min_{s,t}|{\bf z}(s,t,\lambda)|\buildrel
def\over{=}\delta(\lambda),$

for all $\lambda\in[0,\lambda_{0}]$, and that this minimum point is non-
degenerate with respect to $t$, in the sense that

(19) $\frac{\partial^{2}}{\partial t^{2}}|z(0,t,\lambda)|_{\mid t=0}\neq 0.$

Moreover, we make the following orthogonality assumption:

(20) $\dot{\bf x}(0,\lambda)\cdot\dot{\bf y}(0,\lambda)=0,\,\textrm{ for all
}\lambda\in[0,\lambda_{0}].$

In the sequel we will study the case when the minimum distance
$\min_{s,t}|{\bf z}(s,t,\lambda)|\rightarrow 0$, that is, the shortest
distance from the orbit ${\bf y}(t,\lambda)$ of the massive body to the curve
${\bf x}(s,\lambda)$ drawn by the infinitesimal mass approaches $0$ as
$\lambda\to\lambda_{0}$. See Figure 3. In the curved Sitnikov problem, this
corresponds to the case when $r\to\frac{2R}{1+\varepsilon}$ ($r\to 2R$ when
$\varepsilon=0$).

To write the equation of motion for $s$, we note that the Newtonian
gravitational potential of the particle at ${\bf x}(s)$ is a function of $s$
and $t$ given by

(21) $U(s,t,\lambda)=-|{\bf z}(s,t,\lambda)|^{-1},$

and the evolution of $s$ is governed by the Euler–Lagrange equation
$\frac{d}{dt}L_{\dot{s}}-L_{s}=0$ with the Lagrangian

$L=\frac{1}{2}\dot{s}^{2}-U(s,t,\lambda),$

leading to222to explain this form of the Lagrangian, we note that the
equations for our massless particle are obtained by taking the limit of the
particle of small mass $m$; for such a particle, in the ambient potential $U$,
the Lagrangian is $\frac{1}{2}m\dot{s}^{2}-mU(s,t,\lambda)$ – the factor $m$
in front of $U$ is due to the fact that $U$ is the potential energy of the
unit mass. Dividing the Euler–Lagrange equation by $m$ gives (22).

(22) $\ddot{s}+U^{\prime}(s,t,\lambda)=0,\ \textrm{ where
}^{\prime}=\frac{\partial}{\partial s}.$

Note that $s=0$ is an equilibrium for any $\lambda$, since
$U^{\prime}(0,t,\lambda)=0$ for all $t$ and for all $\lambda$, according to
(18). Linearizing (22) around the equilibrium $s=0$ we obtain

(23) $\ddot{S}+a(t,\lambda)S=0,\ \ a(t+1,\lambda)=a(t,\lambda),$

where $a(t,\lambda)=U^{\prime\prime}(0,t,\lambda)$.

We have the following general result:

######  Theorem 1.

Assume that (18), (19), (20) hold, that ${\bf x}$ and ${\bf y}$ are both
bounded in the $C^{2}$–norm uniformly in $\lambda$, and

(24) $\min_{s,t}\|x(s,\lambda)-y(t,\lambda)\|\to 0\textrm{ as }\lambda\to 0.$

Then there exists an infinite sequence

(25)
$\lambda_{1}>\lambda_{2}\geq\lambda_{3}>\lambda_{4}\geq\cdots>\lambda_{2n-1}>\lambda_{2n}\geq\cdots\rightarrow
0$

such that the equilibrium solution $s=0$ of (22) is linearly strongly
stable333The equilibrium solution is linearly strongly stable if the Floquet
multipliers of (23) lie on the unit circle and are not real, or equivalently,
if the linearized system lies in the interior of the set of stable systems.
for all $\lambda\in(\lambda_{2n},\lambda_{2n-1})$. Furthermore, each
complementary $\lambda$–interval contains points where the linearized
equilibrium is not strongly stable, i.e. is either hyperbolic or parabolic.

The proof of Theorem 1 relies on Lemmas 2, 3 and 4 stated below.

######  Lemma 2.

Consider the linear system

(26) $\ddot{x}+a(t,\lambda)x=0,$

where $a(t,\lambda)$ is a continuous function of $t\in[0,1]$,
$\lambda\in(0,\bar{\lambda}]$, where $\bar{\lambda}>0$. Let
$z(t,\lambda)=x+i\dot{x}$, where $x=x(t,\lambda)$ is a nontrivial solution of
(26). Assume that there exists an interval
$[t_{0}(\lambda),t_{1}(\lambda)]\subset[0,1]$, possibly depending on
$\lambda$, such that

(27) $\lim_{\lambda\downarrow 0}\arg
z(t,\lambda)\biggl{|}_{t_{0}}^{t_{1}}\rightarrow-\infty.$

Here $\arg z(t,\lambda)$ is defined as a continuous function of $t$. Although
this choice of $\arg$ is unique only modulo $2\pi$, its increment as stated in
equation (27) is uniquely defined.

Then there exists a sequence $\\{\lambda_{k}\\}_{k=0}^{\infty}$ satisfying
(25) such that the Floquet matrix of (23) is strongly stable for all
$\lambda\in(\lambda_{2n},\lambda_{2n-1})$, for any $n>0$. Furthermore, every
complementary $\lambda$–interval contains values of $\lambda$ for which the
Floquet matrix is not strongly stable.

Proof of this lemma can be found in [15].

######  Lemma 3.

Consider the linear system

(28) $\ddot{x}+a(t)x=0,$

and assume that $a(t)>0$ on some interval $[t_{0},t_{1}]$. For any
(nontrivial) solution $x(t)$ the corresponding phase vector $z(t)=x+i\dot{x}$
rotates by

(29) $\theta[z]\buildrel def\over{=}\arg
z(t)\biggl{|}_{t_{0}}^{t_{1}}\leq-\min_{[t_{0},t_{1}]}\sqrt{a(t)}(t_{1}-t_{0})+\pi.$

###### Proof.

Writing the differential equation $\ddot{x}+a(t)x=0$ as a system $\dot{x}=y,\
\dot{y}=-a(t)x$, we obtain (using complex notation):

$\dot{\theta}=\frac{d}{dt}\arg(x+iy)=-\frac{d}{dt}\hbox{Im}(\ln
z)=-\hbox{Im}\frac{\dot{z}}{z}=-\hbox{Im}\frac{(y-iax)(x-iy)}{x^{2}+y^{2}}=-(a\cos^{2}\theta+\sin^{2}\theta).$

We conclude that for any solution of (28), the angle $\theta=\arg(x+iy)$
satisfies

$\dot{\theta}\leq-(a_{m}\cos^{2}\theta+\sin^{2}\theta),\ \ \hbox{for}\ \
t\in[t_{0},t_{1}],\ \ \hbox{where}\ \ a_{m}=\min_{[t_{0},t_{1}]}a(t).$

To invoke comparison estimates, consider $\bar{\theta}(t)$ which satisfies

(30)
$\frac{d}{dt}\bar{\theta}=-(a_{m}\cos^{2}\bar{\theta}+\sin^{2}\bar{\theta}),\
\ \bar{\theta}(t_{0})=\theta(t_{0}).$

By the comparison estimate, we conclude:

(31)
$\theta\biggl{|}_{t_{0}}^{t_{1}}\leq\bar{\theta}\biggl{|}_{t_{0}}^{t_{1}},$

and the proof of the lemma will be complete once we show that $\bar{\theta}$
satisfies the estimate (29). To that end we consider a solution of

$\ddot{\bar{x}}+a_{m}\bar{x}=0$

with the initial condition satisfying

(32) $\arg(\bar{x}(t_{0})+i\dot{\bar{x}}(t_{0}))=\bar{\theta}(t_{0}).$

This solution is of the form

${\bar{x}}(t)=A\cos(\sqrt{a_{m}}t-\varphi),\ \ A={\rm const.},$

where $\varphi$ is chosen so as to satisfy (32).

Since $\arg(\bar{x}+i\dot{\bar{x}})$ satisfies the same differential equation
as $\bar{\theta}$, and since the initial conditions match, we conclude that
$\bar{\theta}(t)=\arg(\bar{x}+i\dot{\bar{x}})$, so that

(33)
$\bar{\theta}(t)=\arg(\cos(\sqrt{a_{m}}t-\varphi)-i\sqrt{a_{m}}\sin(\sqrt{a_{m}}t-\varphi))=-(\sqrt{a_{m}}t-\varphi)+\widehat{\pi/2},$

where $\widehat{X}$ denotes a quantity whose absolute value does not exceed
$X$. In other words, $\bar{\theta}(t)$ is given by a linear function with
coefficient $-\sqrt{a_{m}}$, up to an error $<\pi/2$. The last inequality is
due to the fact that the complex numbers
$\cos(\sqrt{a_{m}}t-\phi)-i\sin(\sqrt{a_{m}}t-\phi)$ and
$\cos(\sqrt{a_{m}}t-\phi)-i\sqrt{a_{m}}t\sin(\sqrt{a_{m}}t-\phi)$ lie in the
same quadrant, so the difference between their arguments is no more than
$\pi/2$.

Therefore, over the interval $[t_{0},t_{1}]$ the function $\bar{\theta}$
changes by the amount $\sqrt{a_{m}}(t_{1}-t_{0})$ with the error of at most
$\frac{\pi}{2}+\frac{\pi}{2}=\pi$:

(34)
$\bar{\theta}\biggl{|}_{t_{0}}^{t_{1}}\leq-\sqrt{a_{m}}(t_{1}-t_{0})+\pi;$

restating this more formally, (33) implies

$\bar{\theta}(t_{1})<-(\sqrt{a_{m}}t_{1}-\varphi)+\frac{\pi}{2},\ \ \
\bar{\theta}(t_{0})>-(\sqrt{a_{m}}t_{0}-\varphi)-\frac{\pi}{2}.$

Subtracting the second inequality from the first gives (33).

Substituting (34) into (31) yields (29) and completes the proof of Lemma 3.
$\diamondsuit$ ∎

######  Lemma 4.

Consider the potential $U$ defined by equation (21). Assume that the functions
${\bf x}={\bf x}(s,\lambda)$, ${\bf y}={\bf y}(t,\lambda)$ satisfy

  * (i)

the minimum $\min_{s,t}|{\bf z}(\cdot,\cdot,\lambda)|=\delta(\lambda)$ is non-
degenerate and is achieved at $s=t=0$, where ${\bf z}(s,t,\lambda)={\bf
x}(s,\lambda)-{\bf y}(t,\lambda)$,

  * (ii)

$||{\bf z}(\cdot,\cdot,\lambda)||_{C^{2}}\leq M$ uniformly in
$0<\lambda<\lambda_{0}$, and

  * (iii)

$\delta(\lambda)\rightarrow 0$ as $\lambda\rightarrow 0$.

Then there exists time $\tau=\tau(\lambda)$ (approaching zero as
$\lambda\rightarrow 0$) such that

(35) $\lim_{\lambda\rightarrow
0}\left(\tau(\lambda)\cdot\min_{|t|\leq\tau(\lambda)}\sqrt{U^{\prime\prime}(0,t,\lambda)}\right)\rightarrow\infty.$

###### Proof.

Differentiating (21) with respect to $s$ twice, we get

(36) $U^{\prime\prime}(0,t,\lambda)=\left[\frac{({\bf z}^{\prime}\cdot{\bf
z}^{\prime}+{\bf z}\cdot{\bf z}^{\prime\prime})({\bf z}\cdot{\bf z})-3({\bf
z}\cdot{\bf z}^{\prime})^{2}}{({\bf z}\cdot{\bf z})^{5/2}}\right]_{\mid s=0}.$

We now estimate all the dot products in the above expression to obtain a lower
bound.

First,

(37) ${\bf z}^{\prime}\cdot{\bf z}^{\prime}=1,$

since $s$ is the arc length, and from here ${\bf z}^{\prime}\cdot{\bf
z}^{\prime\prime}=0$. Now, to estimate ${\bf z}\cdot{\bf z}$ and ${\bf
z}\cdot{\bf z}^{\prime}$ we observe that ${\bf z}\cdot{\bf
z}^{\prime}=\frac{1}{2}({\bf z}\cdot{\bf z})^{\prime}$ and we note that the
first expression, as a function of $t$ with $s=0$ fixed has a minimum at $t=0$
that we call $\delta^{2}$, and that the second function vanishes at $t=0$.
Applying Taylor’s formula with respect to $t$ we then have

(38) ${\bf z}\cdot{\bf z}=\delta^{2}+\widehat{kt^{2}},\hskip 19.91684pt{\bf
z}\cdot{\bf z}^{\prime}=\frac{1}{2}({\bf z}\cdot{\bf z})^{\prime}=\widehat{\
kt\ },$

where the constant $k$ is determined by the $C^{2}$–norm $M$ of ${\bf z}$. For
the remaining dot product we have (still keeping $s=0$ and $t$ arbitrary):

(39) $|{\bf z}\cdot{\bf z}^{\prime\prime}|\leq|{\bf z}||{\bf
z}^{\prime\prime}|\buildrel(\ref{eq:zz})\over{\leq}M\sqrt{\delta^{2}+kt^{2}}.$

Using the above estimates in (36), we obtain

(40)
$U^{\prime\prime}(0,t,\lambda)\geq\frac{(1-M\sqrt{\delta^{2}+kt^{2}})\delta^{2}-3k^{2}t^{2}}{(\delta^{2}+kt^{2})^{5/2}}.$

Now we restrict $t$ to have $\delta^{2}+kt^{2}\leq 2\delta^{2}$; this
guarantees that the denominator in (40) does not exceed $(2\delta)^{5/2}$; to
bound the numerator, we further restrict $t$ so that the dominant part
$\delta^{2}-3k^{2}t^{2}\geq\frac{1}{2}\delta^{2}$, thus bounding the numerator
from below by

$(\delta^{2}-3k^{2}t^{2})-M\sqrt{2\delta^{2}}\delta^{2}\geq\frac{1}{2}\delta^{2}-M\sqrt{2}\delta^{3}>\frac{1}{4}\delta^{2}$

if $\delta$ is sufficiently small. Summarizing, we restricted $t$ to

(41) $|t|\leq c\delta\buildrel def\over{=}\tau(\lambda),\ \ \hbox{where}\ \
c=\min(k^{-1/2},(k\sqrt{6})^{-1}),$

and showed that for all such $t$ and for $\delta$ small enough

(42)
$U^{\prime\prime}(0,t,\lambda)\geq\frac{\frac{1}{4}\delta^{2}}{(2\delta^{2})^{5/2}}=\frac{c_{1}}{\delta^{3}},$

where $c_{1}=2^{-9/2}$. With $\tau$ defined in (41) we obtain $\lim_{\tau\to
0}\left(\tau\cdot\min_{|t|\leq\tau}\sqrt{a(t,\lambda)}\right)=\infty$, thus
completing the proof of the lemma. $\diamondsuit$ ∎

Proof of Theorem 1 Consider the linearized equation

(43) $\ddot{x}+U^{\prime\prime}(0,t,\lambda)x=0,$

and consider the phase point $z=x+i\dot{x}$ of a nontrivial solution. Lemma 3
gives us the rotation estimate (29) for any time interval $[t_{0},t_{1}]$; let
this interval be $[-\tau(\lambda),\tau(\lambda)]$ where $\tau(\lambda)$ is
taken from the statement of Lemma 4. We then have from (29): (29)

(44) $\theta[z]\buildrel def\over{=}\arg
z(t)\biggl{|}_{t_{0}}^{t_{1}}\leq-2\tau(\lambda)\cdot\min_{|t|\leq\tau(\lambda)}\sqrt{U^{\prime\prime}(0,t,\lambda)}+\pi.$

According to the conclusion (35) of Lemma 4, $\theta\rightarrow-\infty$ as
$\lambda\rightarrow 0$. This satisfies the condition (27) Lemma 2, which now
applies and its conclusion comples the proof of Theorem 1. $\diamondsuit$

##  4\. Stability of the equilibrium points in the curved Sitnikov problem

The system (8) has two equilibria $(0,0)$ and $(\pi,0)$, which correspond to
periodic orbits for $\mathcal{X}_{\varepsilon}$. The associated linear system
around the fixed point $(q_{*},p_{*})$ can be written as

(45)
$\dot{\mathbf{v}}=A(t)\mathbf{v}\,,\qquad\mathbf{v}=\left(\begin{array}[]{c}x\\\
y\\\ \end{array}\right)\,,$

with

$A(t)=\left.\left(\begin{array}[]{cc}0&1\\\
\frac{\partial{f}_{\varepsilon}}{\partial q}&0\\\
\end{array}\right)\right|_{q=q_{*}}.$

Let $X(t)$ be a fundamental matrix solution of system (45) given by

$X(t)=\left(\begin{array}[]{cc}x_{1}(t)&x_{2}(t)\\\ y_{1}(t)&y_{2}(t)\\\
\end{array}\right),$

with the initial condition $X(0)=I$, the identity matrix; $x_{1}$ is an even
function and $x_{2}$ and odd one since $\partial{f}_{\varepsilon}/\partial q$
is an even function with respect to $t$. The monodromy matrix is given by
$X(2\pi)$ and we denote $\lambda_{1},\lambda_{2}$ its eigenvalues, the Floquet
multipliers associated to (45). These are given by

(46)
$\lambda_{1},\lambda_{2}=\frac{x_{1}(2\pi)+y_{2}(2\pi)\pm\sqrt{\left(x_{1}(2\pi)+y_{2}(2\pi)\right)^{2}-4}}{2}\,,$

the trace $\textrm{Tr}(X(2\pi))=x_{1}(2\pi)+y_{2}(2\pi)$ determines the
linearized dynamics around the fixed point. Moreover, since the function
$(\partial{f}_{\varepsilon}/\partial q)_{\mid q=q^{*}}$ is an even function,
we know $x_{1}(2\pi)=y_{2}(2\pi)$ (see [24]) and then

(47)
$\lambda_{1},\lambda_{2}=y_{2}(2\pi)\pm\sqrt{\left(y_{2}(2\pi)\right)^{2}-1}\,.$

To emphasize the dependence on the parameters $\varepsilon,r$ we write

(48) $y_{2}(2\pi;\varepsilon,r)=y_{2}(2\pi).$

Thus the linear stability of the equilibrium is

  1. (1)

Elliptic type: $|y_{2}(2\pi;\varepsilon,r)|<1$.

  2. (2)

Parabolic type: $|y_{2}(2\pi;\varepsilon,r)|=1$.

  3. (3)

Hyperbolic type: $|y_{2}(2\pi;\varepsilon,r)|>1$.

Since the Wronskian is equal to $1$ for all $t$, then

$W(\varepsilon,r)=(y_{2}(2\pi;\varepsilon,r))^{2}-x_{2}(2\pi;\varepsilon,r)y_{1}(2\pi;\varepsilon,r)=1,$

and we have:

(49)
$\left(y_{2}(2\pi;\varepsilon,r)\right)^{2}=1+x_{2}(2\pi;\varepsilon,r)y_{1}(2\pi;\varepsilon,r)\,.$

From this last expression it follows that, in the parabolic case, the periodic
orbit corresponding to the equilibrium point $(\pi,0)$ is associated to
$x_{1}(t)$ or to $x_{2}(t)$ (or to both).

###  4.1. Stability of the equilibrium point $(\pi,0)$. Case $\varepsilon=0$.

In this subsection we consider the system (8)-(15) for the case
$\varepsilon=0$, namely when the primaries are following circular
trajectories. We also fix $R=1$. The extended vector field under study is

(50) $\mathcal{X}_{0}(q,p,t;\gamma)=\left\\{\begin{array}[]{rcl}\dot{q}&=&p\\\
\dot{p}&=&{f}_{0}(q,t;r)\\\ \dot{t}&=&1\end{array}\right.$

where

(51)
${f}_{0}(q,t;r):=-\frac{(1+r\cos(t))\sin(q)}{||\mathbf{x}-\mathbf{x}_{1}||^{3}}-\frac{(1-r\cos(t))\sin(q)}{||\mathbf{x}-\mathbf{x}_{2}||^{3}}\,,\\\
$ $\displaystyle||\mathbf{x}-\mathbf{x}_{1}||$ $\displaystyle=$
$\displaystyle\left[r^{2}+2+2r\cos(t)-2\cos(q)-2r\cos(q)\cos(t))\right]^{1/2}\,,$
$\displaystyle||\mathbf{x}-\mathbf{x}_{2}||$ $\displaystyle=$
$\displaystyle\left[r^{2}+2-2r\cos(t)-2\cos(q)+2r\cos(q)\cos(t))\right]^{1/2}\,,$

with $0<r<2$. We observe that in this case, function (51) is $\pi$–periodic
(remember that ${f}_{\varepsilon}(q,t;r)$ is $2\pi$–periodic if
$\varepsilon>0$).

The linear system associated to (50) around the fixed point $(\pi,0)$ is
defined by the function

(52) $\frac{\partial{f}_{0}}{\partial
q}(\pi,t;r)=\frac{(1+r\cos(t))}{\left[r^{2}+4+4r\cos(t)\right]^{3/2}}+\frac{(1-r\cos(t))}{\left[r^{2}+4-4r\cos(t)\right]^{3/2}}\,,$

which is $C^{1}$ respect to $t$ and $r$.

######  Lemma 5.

The function $\frac{\partial{f}_{0}}{\partial q}(\pi,t;r)$ is monotone
decreasing with respect to $r$, that is, for all $t\in[0,\pi)$

$\frac{\partial{f}_{0}}{\partial
q}(\pi,t;r_{1})>\frac{\partial{f}_{0}}{\partial q}(\pi,t;r_{2})$

if $r_{1}<r_{2}$ (see Figure 4).

###### Proof.

Let be $F(t,r)=\frac{\partial{f}_{0}}{\partial q}(\pi,t;r)$, then by
straightforward computation we get

$\displaystyle\frac{\partial F}{\partial r}(t,r)$ $\displaystyle=$
$\displaystyle\frac{(\cos(t))[r^{2}+4r\cos(t)+4]-3(1+r\cos(t))[r+2\cos(t)]}{\left[r^{2}+4+4r\cos(t)\right]^{5/2}}$
$\displaystyle+$
$\displaystyle\frac{(-\cos(t))[r^{2}-4r\cos(t)+4]-3(1-r\cos(t))[r-2\cos(t)]}{\left[r^{2}+4-4r\cos(t)\right]^{5/2}}$
$\displaystyle\leq$
$\displaystyle\frac{-4r\cos(t)^{2}-6r}{\min\\{\left[r^{2}+4+4r\cos(t)\right]^{5/2},\left[r^{2}+4-4r\cos(t)\right]^{5/2}\\}}<0.$

$\diamondsuit$ ∎

Figure 4. Plots for $\partial\bar{f}_{0}/\partial q$ varying $r$ from $0$ to
$2$.

######  Proposition 6.

The equilibrium point $(\pi,0)$ is of hyperbolic type if
$r\leq\left(\sqrt{17}-3\right)^{1/2}=1.059\dots$.

###### Proof.

We first show that $F(t,r)\geq 0$ for all $t$ if and only if
$r\leq\left(\sqrt{17}-3\right)^{1/2}$. We compute

(53) $\begin{split}\frac{\partial F}{\partial t}(t,r)=r\sin
t\left[-\frac{1}{(r^{2}+4r\cos t+4)^{1/2}}+\frac{1}{(r^{2}-4r\cos
t+4)^{1/2}}\right.\\\ \left.+\frac{6(1+r\cos t)}{(r^{2}+4r\cos
t+4)^{5/2}}-\frac{6(1-r\cos t)}{(r^{2}-4r\cos t+4)^{5/2}}\right].\end{split}$

We have $F(t,r)=F\left(\pi-t,r\right)$ so it is enough to restrict
$t\in[0,\frac{\pi}{2}]$. Note that $\frac{\partial F}{\partial t}(t,r)=0$ for
$t=0,\pi/2$.

For $0<r<2$, we have

$F\left(\frac{\pi}{2},r\right)=\frac{2}{(r^{2}+4)^{3/2}}>0,$

and

$F(0,r)=\frac{1+r}{(2+r)^{3}}+\frac{1-r}{(2-r)^{3}}=\frac{-2(r^{4}+6r^{2}-8)}{(4-r^{2})^{3}}.$

It follows immediately that $F(0,r)<0$ if $r>(\sqrt{17}-3)^{1/2}$. If $r\cos
t\leq 1$ then (52) implies $F(t,r)\geq 0$. Hence $F(t,r)\geq 0$ for all $t$
provided $r\leq 1$. Let $1<r\leq(\sqrt{17}-3)^{1/2}$. If
$t\in[0,\cos^{-1}(1/r)]$, which is equivalent to $r\cos t>1$, then

$\begin{split}\frac{1}{(r^{2}-4r\cos t+4)^{1/2}}-\frac{1}{(r^{2}+4r\cos
t+4)^{1/2}}\geq 0,\\\ \frac{6(1+r\cos t)}{(r^{2}+4r\cos
t+4)^{5/2}}-\frac{6(1-r\cos t)}{(r^{2}-4r\cos t+4)^{5/2}}\geq 0.\end{split}$

Therefore (53) implies that $F(t,r)$ is increasing in $t$ for
$t\in[0,\cos^{-1}(1/r)]$, and, since $F(0,r)\geq 0$, it follows that
$F(t,r)\geq 0$ for all $t$.

Now let $x(t)=x_{1}(t)+x_{2}(t)$ be a particular solution of system (45); by
hypothesis it satisfies $x(0)=1,\dot{x}(0)=1$. From (45) with $\varepsilon=0$
we obtain

$\ddot{x}=\frac{\partial{f}_{0}}{\partial q}(\pi,t;r)x\quad{\rm
or}\quad\ddot{x}-F(t,r)x=0.$

Then, since $F(t,r)\leq 0$ for all $t$, we can apply directly Lyapunov’s
instability criterion (see for instance page 60 in [4]) to show that $(\pi,0)$
is of hyperbolic type. $\diamondsuit$∎

We can in fact estimate the first value $r_{1}$ of $r$ at which the
equilibrium point $(\pi,0)$ becomes of parabolic type for the first time. The
idea is to extend slightly the result beyond
$r=\left(\sqrt{17}-3\right)^{1/2}$, and then find the maximum of the $r$ for
which Proposition 6 holds. In the way we have to do straightforward analytic
but tedious computations that we decide to avoid in this paper. Finally we get
that the value of $r$ to have the first parabolic solution is $r_{1}\approx
1.2472\cdots$.

Figure 5. Graphs of $y_{2}(\pi;\varepsilon,r)$. The plot on the right is a
magnification of the one on the left for r close to 2

In Figure 5 we show a couple of numerical simulations which illustrate the
stability interchanges of the equilibrium point $(\pi,0)$ for $r\in(0,2)$ and
$\varepsilon=0$. The figure on the right hand side is a plot for $r$ close to
2.

###  4.2. Stability of the equilibrium point $(\pi,0)$. Case $\varepsilon\neq
0$.

In this case, for every $0<r<(\sqrt{17}-3)^{1/2}$, $\partial f_{0}/\partial
q(\pi,t;r)>0$, hence $\partial f_{\varepsilon}/\partial q(\pi,t;r)>0$ for all
$\varepsilon>0$ sufficiently small (depending on $r$). Therefore $(\pi,0)$
remains an equilibrium point of hyperbolic type in the case when the primaries
move on Keplerian ellipses of sufficiently small eccentricity $\varepsilon>0$.

We remark that Theorem 1 does not depend of the shape of the curves where
${\bf x}(x)$ or ${\bf y}(t)$ are moving, in other words we can apply Theorem 1
independently if the primaries are moving on a circle or on ellipses of
eccentricity $\varepsilon>0$. Thus we obtain the following result on stability
interchanges:

######  Theorem 7.

In the curved Sitnikov problem, let us fix any $\varepsilon\in[0,1)$, let
$R=1$, and consider $r$ (the semi-major axes of the Keplerian ellipses traced
out by the primaries) as the parameter. As $r$ approaches
$\frac{2}{1+\varepsilon}$, the distance between a Keplerian ellipse and the
circle of the massless particle approaches zero. There exists a sequence
$r_{n}\uparrow\frac{2}{1+\varepsilon}$ satisfying

(54) $r_{0}\leq r_{1}<r_{2}\leq r_{3}\cdots<r_{2n}\leq r_{2n+1}<\cdots,$

such that the equilibrium point $(\pi,0)$ of equation (45) is strongly stable
for $r\in(r_{2n-1},r_{2n})$ and not strongly stable for some $r$ in the
complementary intervals $(r_{2n},r_{2n+1})$. In other words, the equilibrium
point $(\pi,0)$ loses and then regains its strong stability infinitely many
times as $r$ increases towards $\frac{2}{1+\varepsilon}$.

###### Proof of Theorem 7.

We verify that Theorem 1 applies. The curve ${\bf x}(s,\lambda)$ of Theorem 1
is represented by the circle of radius $R=1$, the curve ${\bf y}(t,\lambda)$
is represented by the orbit of the primary that gets closer to $(\pi,0)$ (when
$\varepsilon=0$ the primaries are co-orbital), and the parameter $\lambda$
corresponds to $r$. The planes of the two curves are perpendicular, as in
(20), and the minimum distance between the curves, given by (18), is non-
degenerate as in (19), and it corresponds to the infinitesimal mass being at
$y=-1$ and the closest primary to this point being at $y=1-r(1+\varepsilon)$.
Thus, the minimum distance is $\delta(r)=2-r(1+\varepsilon)$, and it
approaches $0$ when $r\to\frac{2}{1+\varepsilon}$. To apply Theorem 1 we only
need to verify that ${\bf x}$ and ${\bf y}$ are bounded in the $C^{2}$-norm
uniformly in $r$. This is obviously true for ${\bf y}(t,r)$ since the motion
of the primary is not affected by the motion of the infinitesimal mass; it is
also true for ${\bf x}(s,r)$ since $s=$arclength and the motion lies on the
circle of radius $R=1$, hence $\|{\bf x}(s,r)\|=1$ (the radius of the circle),
$\|{\bf x}^{\prime}(s,r)\|=1$ (the unit speed of a curve parametrized by arc-
length), and $\|{\bf x}^{\prime\prime}(s,r)\|=1$ (the curvature of the
circle). Hence the conclusion of Theorem 1 follows immediately.
$\diamondsuit$∎

###  4.3. Stability of the equilibrium point $q=(0,0)$.

The linear stability of $(0,0)$ is the same as of the barycenter in the
classical Sitnikov problem, since

$\dfrac{\partial{f_{\varepsilon}}}{\partial
q}(0)=-\frac{2R}{r^{3}\rho(t;\varepsilon)^{3}}.$

The stability of this point can be treated very similarly to that of the
origin for Hill’s equation, so in this analysis we use results from that
theory.

As in the study of the other equilibrium point we start with the case
$\varepsilon=0.$ Here the function $f_{0}(q,t;r)$ defined in equation (8)
around $q=(0,0)$ is given by

(55) $\frac{\partial f_{0}}{\partial
q}(q)=-\frac{2}{r^{3}}q+\frac{9+r^{2}+9r^{2}(\cos(t))^{2}}{3r^{5}}q^{3}+\mathcal{O}(q^{5})\,.$

The local dynamics is determined by the linear part. The eigenvalues are on
the unit circle, given by $\sigma_{1,2}=\pm i\sqrt{2/r^{3}}$, and the Floquet
multipliers, which come from the monodromy matrix

$X(\pi)=\left(\begin{array}[]{cc}\cos\left(\sqrt{\frac{2}{r^{3}}}\pi\right)&\left(\sqrt{\frac{r^{3}}{2}}\right)\sin\left(\sqrt{\frac{2}{r^{3}}}\pi\right)\\\
-\left(\sqrt{\frac{2}{r^{3}}}\pi\right)\sin\left(\sqrt{\frac{2}{r^{3}}}\pi\right)&\cos\left(\sqrt{\frac{2}{{}^{3}}}\pi\right)\\\
\end{array}\right),\qquad X(0)=I\,,$

are $\lambda_{1,2}=e^{\pm i\sqrt{2/r^{3}}\pi}$.

Hence $q=(0,0)$ is of elliptic type if $\sqrt{2/r^{3}}\neq k$ for any
$k\in\mathbb{Z}$, and it is of parabolic type if $\sqrt{2/r^{3}}=k$ for some
$k\in\mathbb{Z}$. In fact, in the parabolic case, if $\sqrt{2/r^{3}}=2m$,
$m\in\mathbb{Z}$, then there exists a $\pi$-periodic solution, and if
$\sqrt{2/r^{3}}=2m+1$, $m\in\mathbb{Z}$, then there exists a $2\pi$-periodic
solution.

######  Proposition 8.

The equilibium point $q=(0,0)$ of the system defined by (55) is stable for all
$r\in(0,2)$.

###### Proof.

The linear part possess Floquet multipliers which place the system in either
the elliptic or the parabolic case. In the last case, when $\sqrt{2/r^{3}}=k$
for some $k\in\mathbb{Z}$, there exists two independent eigenvectors
associated to the eigenvalues $\lambda_{1,2}$. Then the conjugacy of $A$ and
$\pm I$ implies that the monodromy matrix is $A=\pm I$ and $q=(0,0)$ is stable
(see [20] for more details).

When the origin is of elliptic type for the linear system, we observe that the
coefficient of the term of order $3$ in equation (55) is greater than zero for
all $r\in(0,2)$, then we can apply directly Ortega’s theorem (see Appendix and
[19]), and therefore we obtain that the equilibrium point $(0,0)$ is stable
for the whole system, that is, including the nonlinear part. $\diamondsuit$ ∎

We note that, in the case when $\varepsilon=0$, stability interchanges in a
weak sense appear as $r\to 0$, since $(0,0)$ switches between elliptic type
when $r\neq(2/k^{2})^{1/3}$ and parabolic type when $r=(2/k^{2})^{1/3}$,
$k\in\mathbb{Z}^{+}$, as noted before.

In the case when $\varepsilon\neq 0$, as we mentioned earlier, the linear
stability of $(0,0)$ is the same as in the classical Sitnikov problem, so it
only depends on the eccentricity parameter $\varepsilon$. The papers [16, 11,
12] state that there are stability interchanges when the size $r$ of the
binary is kept fixed and the eccentricity $\varepsilon$ of the Keplerian
ellipses approaches $1$. We should point out that Theorem 1 does not apply to
this case since the function $r_{\varepsilon}(t,r)$ describing the motion of
the primaries — corresponding to $y(t,\lambda)$ in Theorem 1 — does not remain
bounded in the $C^{2}$ norm uniformly in $\varepsilon$.

##  Appendix A Floquet Theory

In order to have a self contained paper we add this appendix with the main
results on Floquet theory, must of them are very well known for people in the
field.

Consider the linear system

(56) $\dot{\mathbf{x}}=A(t)\mathbf{x}\,,\qquad\mathbf{x}\in\mathbb{R}^{2}\,,$

where $A(t)$ is a $T$-periodic matrix-valued function. Let $X(t)$ be the
fundamental matrix solution

(57) $X(t)=\left(\begin{array}[]{cc}x_{1}(t)&x_{2}(t)\\\ y_{1}(t)&y_{2}(t)\\\
\end{array}\right)$

with the initial condition $X(0)=I$, the identity matrix. Let $\lambda_{1}$
and $\lambda_{2}$ be the eigenvalues (Floquet multipliers) of the monodromy
matrix matrix $X(T)$ and let $\mu_{1},\mu_{2}$ (Floquet exponents) be such
that $\lambda_{1}=e^{\mu_{1}T},\lambda_{2}=e^{\mu_{2}T}$.

######  Theorem 9 (Floquet’s theorem).

Suppose $X(t)$ is a fundamental matrix solution for (56), then

$X(t+T)=X(t)X(T)$

for all $t\in\mathbb{R}$. Also there exists a constant matrix $B$ such that
$e^{TB}=X(T)$ and a $T$-periodic matrix $P(t)$, so that, for all $t$,

$X(t)=P(t)e^{Bt}.$

######  Lemma 10.

Let be $\lambda$ a Floquet multiplier for (56) and $\lambda=e^{\mu T}$, then
there exists a nontrivial solution $x(t)=e^{\mu t}p(t)$, with $p(t)$ a
$T$-periodic function. Moreover, $x(t+T)=\lambda x(t)$.

Thus, the Floquet multipliers lead to the following characterization:

  * •

If $|\lambda|<1\Leftrightarrow\textrm{Re}(\mu)<0$ then $x(t)\rightarrow 0$ as
$t\rightarrow\infty$.

  * •

If $|\lambda|=1\Leftrightarrow\textrm{Re}(\mu)=0$ then $x(t)$ is a pseudo-
periodic, bounded solution. In particular when $\lambda=1$ then $x(t)$ is
$T$-periodic and when $\lambda=-1$ then $x(t)$ is $2T$-periodic.

  * •

If $|\lambda|>1$ then $\textrm{Re}(\mu)>0$ and therefore
$x(t)\rightarrow\infty$ as $t\rightarrow\infty$, an unbounded solution.

######  Lemma 11.

If the Floquet multipliers satisfy $\lambda_{1}\neq\lambda_{2}$, then the
equation (56) has two linearly independent solutions

$x_{1}(t)=p_{1}(t)e^{\mu_{1}t}\,,\qquad x_{2}(t)=p_{2}(t)e^{\mu_{2}t}\,,$

where $p_{1}(t)$ and $p_{2}(t)$ are $T$-periodic functions and $\mu_{1}$ and
$\mu_{2}$ are the respective Floquet exponents.

In this way, the stability of the solution to (56) is

  * •

Asymptotically stable if $|\lambda_{i}|<1$ for $i=1,2$.

  * •

Lyapunov stable if $|\lambda_{i}|\leq 1$ for $i=1,2$, or if $|\lambda_{i}|=1$
and the algebraic multiplicity equals the geometric multiplicity.

  * •

Unstable if $|\lambda_{i}|>1$ for at least one $i$, or if $|\lambda_{i}|=1$
and the algebraic multiplicity is greater than the geometric multiplicity.

A particular case for the equation (56) is the so-called Hill’s equation,
namely the periodic linear second order differential equation:

(58) $\ddot{z}+f(t)z=0\,,$

where $f(t)$ is a $\pi$-periodic function. Equation (58), as a first order
system, is

(59)
$\dot{\mathbf{v}}=A(t)\mathbf{v}\,,\qquad\mathbf{v}=\left(\begin{array}[]{c}x\\\
y\\\ \end{array}\right)\,,\qquad A(t)=\left(\begin{array}[]{cc}0&1\\\
-f(t)&0\\\ \end{array}\right)\,.$

Let (57) be a fundamental matrix solution of (59). The monodromy matrix
corresponds to $X(\pi)$ and, as before, let $\lambda_{1}$ and $\lambda_{2}$ be
the Floquet multipliers associated to (59) and $\mu_{1}$ and $\mu_{2}$ be the
corresponding Floquet exponents.

In this paper we assume that $f(t)$ is an even function; from Floquet theory
we know that $x_{1}$ is an even and $x_{2}$ is an odd function. Also, the
trace of $A(t)$ vanishes and

$\lambda_{1}\lambda_{2}=e^{\int_{0}^{\pi}\textrm{tr}(A(t))dt}=1\,.$

Assuming $\mu_{1}=a+ib$ is the Floquet exponent corresponding to the Floquet
multiplier $\lambda_{1}$, the general solution to (59) is characterized as
follows:

  * •

Elliptic type: $\lambda_{1}\in\mathbb{C}\setminus\mathbb{R}$, with
$|\lambda_{1}|=1$ (and $\lambda_{2}=\bar{\lambda}_{1}$). The general solution
is pseudo-periodic and can be written as

$\mathbf{v}(t)=c_{1}\textrm{Re}\left(\mathbf{p}(t)e^{ibt/\pi}\right)+c_{2}\textrm{Im}\left(\mathbf{p}(t)e^{ibt/\pi}\right)\,.$

The origin is Lyapunov stable.

  * •

Parabolic type: $\lambda_{1}=\lambda_{2}=\pm 1$.

    * –

If $\lambda_{1}=1$ and there are two linearly independent eigenvectors of the
monodromy matrix, the general solution is

$\mathbf{v}(t)=c_{1}\mathbf{p}_{1}(t)+c_{2}\mathbf{p}_{2}(t)\,,$

and is $\pi$-periodic and Lyapunov stable.

    * –

If $\lambda_{1}=1$ and there is just one eigenvector associated to this
eigenvalue, thus the general solution is

$\mathbf{v}(t)=(c_{1}+c_{2}t)\mathbf{p}_{1}(t)+c_{2}\mathbf{p}_{2}(t)\,,$

and is unstable.

    * –

If $\lambda_{2}=-1$ and there are two linearly independent eigenvectors of the
monodromy matrix, the general solution is

$\mathbf{v}(t)=c_{1}\mathbf{p}_{1}(t)e^{it}+c_{2}\mathbf{p}_{2}(t)e^{it}\,,$

and is $2\pi$-periodic and Lyapunov stable.

    * –

If $\lambda_{1}=-1$ and there is just one eigenvector associated to this
eigenvalue, thus the general solution is

$\mathbf{v}(t)=(c_{1}+c_{2}t)\mathbf{p}_{1}(t)e^{it}+c_{2}\mathbf{p}_{2}(t)e^{it}\,,$

and is unstable.

  * •

Hyperbolic type: $\lambda_{1}\in\mathbb{R}$, but $|\lambda_{1}|\neq 1$ (and
$\lambda_{2}=1/\lambda_{1}$).

    * –

If $\lambda_{1}>1$, then the solution is

$\mathbf{v}(t)=c_{1}\mathbf{p}_{1}(t)e^{\mu_{1}t}+c_{2}\mathbf{p}_{2}(t)e^{-\mu_{1}t}\,.$

    * –

If $\lambda_{1}<-1$, then the solution is

$\mathbf{v}(t)=c_{1}\mathbf{p}_{1}(t)e^{\mu_{1}t}e^{it}+c_{2}\mathbf{p}_{2}(t)e^{-\mu_{1}t}e^{it}\,.$

Thus, the origin is unstable.

A useful result from Rafael Ortega [19] considers the nonlinear Hill equation

(60) $y^{\prime\prime}+a(t)y+c(t)y^{2n-1}+d(t,y)=0,$

with $n\geq 2$, where the functions $a,c:\mathbb{R}\rightarrow\mathbb{R}$ are
continuous, $T$-periodic, and $\int_{0}^{T}|c(t)|\neq 0$, and the function
$d:\mathbb{R}\times(-\epsilon,\epsilon)\rightarrow\mathbb{R}$, for
$\epsilon>0$, is continuous, has continuous derivatives of all orders respect
to $y$, is $T$-periodic function respect to $t$, and
$d(t,y)=\mathcal{O}\left(|y|^{2n}\right)$ as $y\rightarrow 0$ uniformly with
respect to $t\in\mathbb{R}$.

The linear part around the solution $y=0$ of (60) is

(61) $y^{\prime\prime}+a(t)y=0\,.$

######  Theorem 12 (Ortega’s Theorem).

Assume the following:

  1. (1)

The equation (61) is stable.

  2. (2)

$c\geq 0$ or $c\leq 0$.

Then $y=0$ is a stable solution of (60).

### Acknowledgements

We thank to the anonymous referees, their remarks and suggestions help us to
improve this paper. Research of M.G. was partially supported by NSF grant
DMS-1515851. M. L. gratefully acknowledges support by the NSF grant
DMS-1412542. The fourth author (EPC) has received partial support by the
Asociación Mexicana de Cultura A.C. Parts of this work have been done while
the authors visited CIMAT, Guanajuato, LFP and EPC visit Yeshiva University
and M.G. visited UAM-I in Mexico City. All authors are grateful for the
hospitality of these institutions.

## References

  * [1] V. Alekseev, Quasirandom dynamical systems I, Math. USSR Sbornik 5, 73-128 (1968).
  * [2] V. Alekseev, Quasirandom dynamical systems II, Math. USSR Sbornik, 6, 505-560 (1968).
  * [3] V. Alekseev, Quasirandom dynamical systems III, Math. USSR Sbornik, 7, 1-43 (1969).
  * [4] L. Cesari, _Asymptotic behavior and instability problems in ordinary differential equations_ , Springer-Verlag, New York, 1971.
  * [5] J. Chazy, _Sur l’allure finale du mouvement dans le problème des trois corps I, II, III_. Ann. Sci. Ecole Norm. Sup. $3^{e}$ Sér. 39 (1922); J. Math. Pures Appl. 8 (1929); Bull. Astron. 8 (1932).
  * [6] H. Dankowicz, P. Holmes, _The existence of transversal homoclinic points in the Sitnikov problem_. J. Differential Equations 116, 468-483 (1995).
  * [7] F. Diacu, _Relative Equilibria in the 3-Dimensional Curved-N-Body Problem_ , Memoirs of the American Mathematical Society, Vol. 228, 2014.
  * [8] L. Franco-Pérez, E. Pérez Chavela, _Global symplectic regularization for some restricted 3–body problems on $S^{1}$_. Nonlinear Analysis, 71, 5131-5143 (2009).
  * [9] A. Garcia, E. Pérez-Chavela, _Heteroclinic phenomena in the Sitnikov problem_ , Hamiltonian systems and celestial mechanics (Patzcuaro, 1998), 174-185, World Sci. Monogr. Ser. Math., 6, World Sci. Publ., River Edge, NJ, 2000.
  * [10] A. Gorodetski and V. Kaloshin, _Hausdorff dimension of oscillatory motions in the restricted planar circular three body problem and in Sitnikov problem_ , preprint.
  * [11] J. Hagel, C. Lhotka, _A High Order Perturbation Analysis of the Sitnikov Problem_ , Celestial Mechanics and Dynamical Astronomy 93, 201-228 (2005).
  * [12] V. O. Kalas, P. S. Krasil’nikov, _On equilibrium stability in the Sitnikov problem_ , Cosmic Research 49, 534-537 (2011).
  * [13] V.B. Kostov, P.R. McCullough, J.A. Carter, M. Deleuil, R.F. Díaz, D. C. Fabrycky, G. Hébrard, T.C. Hinse, T. Mazeh, J.A. Orosz, Z.I. Tsvetanov, W.F. Welsh, _Kepler-413b: a slightly misaligned, Neptune-size transiting circumbinary planet_ , Arxive: 1401.7275 (2014).
  * [14] T. Kovács, Gy. Bene and T. Tél, _Relativistic effects in the chaotic Sitnikov problem_ , Monthly Notices of the Royal Astronomical Society, 414 (3), 2275-2281 (2011).
  * [15] M. Levi, _Stability of the inverted pendulum - a topological explanation_. SIAM Review, 30 (4), 639-644 (1988).
  * [16] J. Martínez Alfaro, C. Chiralt, _Invariant rotational curves in Sitnikov’s Problem_ , Celestial Mechanics and Dynamical Astronomy, 55 (4), 351–367 (1993).
  * [17] R. McGehee, _A stable manifold theorem for degenerate fixed points with applications to Celestial Mechanics_ , J. Differential Equations 14, 70-88 (1973).
  * [18] J. K. Moser, _Stable and Random Motion in Dynamical Systems_. Annals Math. Studies 77, Princeton University Press, 1973.
  * [19] R. Ortega, _The stability of the equilibrium of a nonlinear Hill’s equation_. SIAM J. Math. Anal. 25 (5), 1393-1401 (1994).
  * [20] R. Ortega, _The stability of the equilibrium: a search for the right approximation_. Ten Mathematical Essays on Aproximation in Analysis and Topology, (J. Ferrera, J. López-Gómez and F.R. Ruiz del Portal, Eds.), Elsevier, 215-234, (2005).
  * [21] H.C. Plummer, _An Introductory Treatise on Dynamical Astronomy_ , New York, Dover, 1960.
  * [22] C. Robinson, _Homoclinic Orbits and Oscillation for the Planar Three-Body Problem_ , Journal of Differential Equations 52, 356-377 (1984).
  * [23] K. Sitnikov, _The Existence of Oscillatory Motions in the Three-Body Problem_. Translation from Doklady Akademii Nauk SSSR, 133 (2), 647-650 (1961).
  * [24] M. Wilhelm, W. Stanley, _Hill’s Equation_. First edition, Dover, 1979.
  * [25] S.L. Ziglin, Non-integrability of the restricted two-body problem on a sphere. Doklady RAN. 379  (4), 477-478 (2001). Engl. transl.: Physics-Doklady. 46  (8), 570-571 (2001).


[ type=editor, role=Researcher, orcid=0000-0002-0939-6932 ]

[1]

[ type=editor, role=Assistant, orcid=0009-0003-8681-3335 ]

1]organization=Leakey Foundation, addressline=1003B O’Reilly Avenue, city=San
Francisco, postcode=94129, state=California, country=United States

2]organization=Graduate Student, University of San Francisco, addressline=2130
Fulton St, city=San Francisco, postcode=94117, state=California,
country=United States

[cor1]Corresponding author

# A separate way to measure rate of return

Gordon Getty<EMAIL_ADDRESS>Nikita Tkachenko
<EMAIL_ADDRESS>[  [

###### Abstract

Net profit is sometimes found from data for net operating surplus. We propose
a way to find it from data for consumption, pay and market-value capital, and
concomitantly to reveal the factor shares in consumption.

###### keywords:

Net profitMarket-value capital Consumption Pay

##  1 Introduction

Many economists are interested in net profit as the capital share of income,
and in rate of return as its ratio to capital invested. Rate of return to
equity and debt investments has been derived from market quotes and dividends,
from census methods as to real estate, and from real interest rates111See e.g.
Jordà et al. (2019), and Homer and Sylla (2005).. We suggest that these
methods can be supplemented by another which reasons from data for
consumption, pay222Compensation of employees plus a share of mixed income. and
market-value capital as reported in national accounts. We argue, that is, that
cash flow can be found as consumption less pay, and that this difference can
then be added to change in market-value capital to give net profit at market
value. Net profit can be divided by market-value capital to give rate of
return. The argument will meanwhile help clarify factor shares in consumption.

##  2 A separate measure of net profit

National accounts and most teaching are organized on the principles that net
output is realized in the forms of consumption and capital growth333Usually
stated as consumption and net investment, but with the understanding that net
investment is realized in capital growth. We measure capital at market value,
not at its cost in depreciated investment., and that this sum equals the
factor claims on it in pay and net profit. Thus net output and the equal
claims on it show as

$C+\Delta K=\Pi+P\quad\text{in national accounts,}$ (1)

where these notations respectively give consumption, capital growth, pay and
net profit. Appendix LABEL:appendix-a1 will question the teaching that these
two sums give net output and claims on net output. Eq. (1) will be confirmed
nonetheless in Section 5 below by another argument, and in Appendix
LABEL:appendix-a2 by another still.

Meanwhile Eq. (1) can be arranged as

$C-\Pi=P-\Delta K\ .$ (2)

$P-\Delta K$, or net profit less capital growth, gives cash flow, or the net
flow of value passed from assets to their owners for consumption or
reinvestment in other assets. We notate cash flow as $F(K)$, and write444We
find no standard notation for cash flow or cash flow rate, as might be
expected in that the cash flow concept is little used in economics outside the
specialty of finance. Our choice $F(K)$ risks misinterpretation in that the
form $F(X)$ tends to mean a unspecified function of $X$. We use $F(K)$, rather
than simply $F$, because we believe that the cash flow from human capital can
also be usefully defined, and can show as $F(H)$ (see Appendix
LABEL:appendix-a1).

$P-\Delta K=F(K)\ .\quad\text{By Eq. \eqref{eq-2}, then}$ (3) $C-\Pi=F(K)\ .$
(4)

By Eqs. (3) and (4), then, it is possible to measure net profit as

$P=\Delta K+F(K)=\Delta K+C-\Pi\ .$ (5)

We define

$r(K)=\frac{P}{K}\ ,\quad g(K)=\frac{\Delta K}{K}\quad\text{and}\quad
f(K)=\frac{F(K)}{K}=\frac{C-\Pi}{K}\ ,$

where $r(K)$ is rate of return, $g(K)$ is capital growth rate, and $f(K)$ is
our notation for cash flow rate555See footnote 4.. These definitions allow
division of Eq. (5) by time to show

$r(K)=g(K)+f(K)=g(K)+\frac{C-\Pi}{K}\ .$ (6)

Table 1 shows $f(K),g(K)$ and $r(K)$ as inferred from data for $C,\Pi$ and $K$
through Eqs. (5) and (6).

Table 1: Average $f(K)$, $g(K)$ and $r(K)$ in all countries, in descending
order of net income per capita(%).

Country Period $f(K)$ $g(K)$ $r(K)$ Luxembourg 1996 - 2018 1.93 4.50 6.43
Switzerland 1995 - 2019 0.88 3.70 4.58 Norway 1981 - 2020 3.55 5.64 9.19 USA
1970 - 2018 4.44 3.56 7.99 Denmark 1996 - 2020 2.16 5.42 7.59 Netherlands 1996
- 2019 2.42 3.35 5.77 Austria 1996 - 2019 2.67 2.59 5.27 Sweden 1950 - 2020
2.59 4.25 6.84 Saudi Arabia 2002 - 2009 10.88 5.73 16.61 Iceland 2000 - 2014
3.27 3.39 6.67 Belgium 1996 - 2019 1.54 2.80 4.33 Germany 1990 - 2020 3.53
3.11 6.64 New Zealand 1999 - 2017 2.46 5.65 8.12 France 1950 - 2019 2.39 4.28
6.67 Ireland 1996 - 2019 2.51 5.61 8.12 Finland 1996 - 2020 3.79 4.06 7.85
Australia 1960 - 2019 2.40 4.45 6.85 Canada 1972 - 2020 2.82 4.45 7.27 Cyprus
1996 - 2019 6.93 4.94 11.87 Israel 2000 - 2019 13.93 8.75 22.68 Malta 1996 -
2019 2.50 4.01 6.51 UK 1970 - 2019 5.11 3.00 8.11 Spain 1995 - 2019 2.72 3.57
6.28 Japan 1980 - 2017 2.28 2.81 5.09 Aruba 1996 - 2000 26.05 28.80 54.85
Venezuela 1997 - 2007 32.48 6.69 39.17 Slovenia 1996 - 2019 1.10 1.72 2.82
Lithuania 1996 - 2019 8.20 4.49 12.69 Czech Rep. 1994 - 2019 2.67 1.07 3.74
Italy 1980 - 2020 4.89 2.84 7.73 Estonia 1996 - 2019 2.98 4.29 7.26 Poland
1996 - 2019 10.72 5.86 16.58 Hungary 1996 - 2019 2.51 3.13 5.63 Portugal 1996
- 2020 2.77 2.54 5.32 Country Period $f(K)$ $g(K)$ $r(K)$ Slovakia 1996 - 2020
5.40 4.15 9.55 Romania 1996 - 2019 7.25 2.54 9.79 Greece 1996 - 2019 26.10
0.63 26.73 Turkey 2009 - 2017 28.14 7.20 35.34 Croatia 1996 - 2019 12.42 4.24
16.65 Latvia 1996 - 2019 4.39 4.30 8.69 Russian Fed. 2008 - 2018 2.82 -1.44
1.38 Uruguay 2012 - 2016 5.21 1.36 6.58 Kazakhstan 1996 - 2019 4.95 6.79 11.74
Malaysia 2006 - 2015 32.47 12.36 44.83 Chile 1996 - 2009 14.27 6.72 20.99
Costa Rica 2012 - 2017 8.84 5.21 14.04 Bulgaria 1996 - 2017 5.10 -0.09 5.00
Mexico 1996 - 2019 10.25 3.28 13.53 Domin. Rep. 1996 - 2016 33.68 8.44 42.12
Serbia 1997 - 2019 30.08 9.62 39.70 Brazil 1996 - 2018 8.97 4.44 13.42
Colombia 1996 - 2019 39.15 8.19 47.34 Azerbaijan 1996 - 2010 39.68 13.76 53.44
Iran 1996 - 2018 21.32 6.93 28.25 Peru 2007 - 2019 8.79 6.44 15.22 Egypt 1996
- 2015 88.06 14.11 102.17 Moldova 1998 - 2019 20.51 4.15 24.66 Mongolia 2005 -
2019 4.29 6.86 11.16 Botswana 1996 - 2000 20.19 6.70 26.89 Guatemala 2001 -
2019 24.90 6.98 31.89 Uzbekistan 2015 - 2017 4.53 6.99 11.52 Cabo Verde 2007 -
2017 22.72 3.82 26.54 Nicaragua 2006 - 2018 6.65 5.11 11.76 Honduras 2000 -
2015 20.88 8.54 29.42 Kyrgyzstan 1996 - 2019 20.79 4.09 24.87 Cote d’Ivoire
1996 - 2000 61.84 29.54 91.38 Senegal 2014 - 2015 28.95 8.68 37.63 Niger 1996
- 2019 67.46 11.87 79.32

##  3 High values of $f(K)$ and $r(K)$ in some countries

Table 1 shows average values of $f(K)$, calculated as $\frac{C-\Pi}{K}$, as
greater than 10% in 26 of the 68 countries reported. Evidence from stock
markets and business experience suggests that cash flow rates greater than
10%, with concomitantly high rates of return, are unlikely at a national
scale. Rather, we suspect that there are conditions under which national
accounts tend to overstate consumption, or understate pay and market value of
capital, or all at once. Our displays list countries in descending order of
net income per capita due to a first impression, which may be wrong, that
conditions for misgauging $C,\Pi$ and $K$ may have something to do with
technological advancement666Only 5 of those 26 countries appear in the top
half of all 68 countries ranked by net income per capita.. Results for those
26 countries remind us that measurements in accounts, including national
accounts, are generally imperfect, and that our derivation of $f(K),g(K)$ and
$r(K)$ from measurements of $C,\Pi$ and $K$ cannot give exact results in any
country at any time.

##  4 The possibility of negative values of $\Delta K$ and $P$

It is expected that $\Delta K$, even at national or global scales, will
occasionally be negative in times of recession, depression or external shocks.
When $-\Delta K$ exceeds $F(K)$ absolutely, then net profit $P$ will also be
negative. Table 1 and line charts in Section 1 of the web appendix777Web
Appendix: https://bridges-and-barriers.quarto.pub/web-appendix-for-a-separate-
way-to-measure-rate-of-return/ show these recurring realities of the
investment world.

##  5 Cash flow, consumption and pay

Net profit might be spent either on investment or consumption. Cash flow,
however, is defined as net of any concurrent investment, whether the source of
that investment is from pay, or from transfer payments received, or
reinvestment of cash flow received. Thus the only possible dispositions of
cash flow at the household scale, where that set includes all existing or
currently acquired assets of the household, are in consumption or in transfer
payments to other households. If net transfer $\Gamma$ means transfer paid out
less transfer received, then, this reasoning gives

$\Pi+F(K)=C+\Gamma\quad\text{at the household scale,}$ (7)

as the source and disposition of payments. At the collective scale, where
transfers offset to zero, then,

$\Pi+F(K)=C\ ,\quad\text{from which}\quad C-\Pi=F(K)\
,\quad\text{collectively.}$ (8)

This confirms Eq. (4), and also Eqs. (1) and (2) by the reasoning shown in Eq.
(3).

##  6 Where our method fits

We see our indirect derivation method for finding $f(K)$ and $r(K)$ as
complementary to methods of finding them by direct research as in Jordà et al.
(2019), which we choose as our current reference standard for that method.
Table 2 compares our findings with theirs. Their research approach gives finer
detail, while ours by inference from $C,\Pi$ and $K$, subject to
mismeasurement problems described in Section 3, gives broader coverage. Ours
describes for capital as a whole, including owner-occupied housing and
government property and personal effects, even if we are not sure what $f(K)$
and $r(K)$ mean for them. Meanwhile the research approach, which must focus on
areas where assets are commonly rented or traded, illuminates those areas in
ways which ours cannot, and with an accuracy that we cannot claim. Neither
consumption, nor pay, nor market-value capital can be measured as closely as
equity prices and dividend rates in finding $f(K)$, $g(K)$ and $r(K)$ for the
corporate sector.

Jordà et al. find $f(K)$ separately for investments in equities, rental
housing, and government bills and bonds as dividend income, rental income, and
real interest rates respectively. Table 2 shows these findings along with our
derivations of $f(K)$, $g(K)$, and $r(K)$ for each country for which Jordà et
al. report. In those 16 countries, our impression is that a weighted average
of the four measures of $f(K)$ in Jordà et al. would come close to our single
one derived from consumption, pay, and market-value capital, which in these
countries seems to hold to a range of some 3% to 6% per year.

Table 2: Comparison of findings from derivation and market research(%). | Derivation from $C,\Pi$ and $K$ | Market Research by Jordà et al.
---|---|---
|  $\frac{C-\Pi}{K}$ |  $\frac{\Delta K}{K}$ |  $\frac{C+\Delta K-\Pi}{K}$ | Equities | Housing | Gvt. Debt
Country | $f(K)$ | $g(K)$ | $r(K)$ |  Div Yld |  Cap Appr |  Total Rtn |  Rent Yld |  Cap Appr |  Total Rtn | Bills | Bonds
Australia | 2.40 | 4.45 | 6.85 | 4.90 | 3.06 | 7.79 | 3.99 | 2.53 | 6.37 | 1.29 | 2.24
Belgium | 1.54 | 2.80 | 4.33 | 3.83 | 2.53 | 6.23 | 6.15 | 1.95 | 7.89 | 1.21 | 3.01
Denmark | 2.16 | 5.42 | 7.59 | 4.95 | 2.71 | 7.49 | 7.13 | 1.26 | 8.22 | 3.08 | 3.58
Finland | 3.79 | 4.06 | 7.85 | 5.08 | 5.19 | 10.03 | 7.14 | 2.82 | 9.58 | 0.64 | 3.22
France | 2.39 | 4.28 | 6.67 | 3.73 | -0.37 | 3.21 | 5.09 | 1.55 | 6.39 | -0.47 | 1.54
Germany | 3.53 | 3.11 | 6.64 | 4.08 | 2.74 | 7.11 | 6.03 | 1.86 | 7.82 | 1.51 | 3.15
Italy | 4.89 | 2.84 | 7.73 | 3.61 | 3.78 | 7.25 | 3.49 | 1.45 | 4.77 | 1.20 | 2.53
Japan | 2.28 | 2.81 | 5.09 | 2.65 | 3.12 | 6.00 | 4.70 | 2.00 | 6.54 | 0.68 | 2.54
Netherlands | 2.42 | 3.35 | 5.77 | 4.87 | 3.38 | 8.10 | 5.96 | 1.75 | 7.51 | 1.37 | 2.71
Norway | 3.55 | 5.64 | 9.19 | 4.21 | 1.61 | 5.67 | 6.72 | 1.49 | 8.03 | 1.10 | 2.55
Portugal | 2.77 | 2.54 | 5.32 | 2.28 | 2.92 | 5.11 | 4.47 | 1.13 | 5.21 | 0.01 | 2.23
Spain | 2.72 | 3.57 | 6.28 | 4.53 | 1.80 | 5.83 | 4.16 | 1.26 | 5.21 | -0.04 | 1.41
Sweden | 2.59 | 4.25 | 6.84 | 4.12 | 4.08 | 8.02 | 7.12 | 1.39 | 8.30 | 1.77 | 3.25
Switzerland | 0.88 | 3.70 | 4.58 | 3.20 | 3.17 | 6.27 | 4.54 | 0.81 | 5.24 | 0.89 | 2.41
UK | 5.11 | 3.00 | 8.11 | 4.53 | 2.48 | 6.83 | 3.94 | 1.63 | 5.44 | 1.16 | 2.29
USA | 4.44 | 3.56 | 7.99 | 4.38 | 4.19 | 8.46 | 5.33 | 0.90 | 6.10 | 2.23 | 2.85
Comments:
1\. Data for equities, housing, bills, and bonds are reproduced from Jordà et
al. (2019), Tables IX and X.
2\. All data are shown from inception through 2015. Dates of inception are
generally earlier for the data from Jordà et al. than from ours. Ours differ
from those shown in Table 1, which do not end in 2015. Exact periods of
coverage of data from Jordà et al. (2019) are shown in Table 1 of that source.

##  7 Factor shares in consumption

At the scale of all assets together, where $F(K)$ can no longer be reinvested
in other capital, $F(K)$ simplifies to the cost of consumption afforded from
net profit, here notated $C(K)$, as distinct from the part of consumption
afforded from pay. Then

$F(K)=C(K)\quad\text{collectively.}$ (9)

It appears from Eq. (4), which can be arranged as $C=\Pi+F(K)$, that pay
equals the remaining part of consumption afforded from pay rather than from
cash flow $F(K)$. If this part is notated $C(\Pi)$, we have

$\Pi=C(\Pi)\quad\text{and}\quad C=C(\Pi)+C(K)\quad\text{collectively,}$ (10)

to give the factor shares in consumption. Division by $C$ gives

$\frac{C(\Pi)}{C}+\frac{C(K)}{C}=\frac{\Pi}{C}+\frac{F(K)}{C}=1\,\quad\text{collectively,}$
(11)

to give those factor shares as percents. Results are shown in Table 3 below.

Table 3: Labor and capital sharesa in consumption in 68 countries in
descending order of income per capita(%).

Country Period Labor Share Capital Share Luxembourg 1996 - 2018 81.53 18.47
Switzerland 1995 - 2019 94.07 5.93 Norway 1981 - 2020 79.97 20.03 United
States 1970 - 2018 81.00 19.00 Denmark 1996 - 2020 89.85 10.15 Netherlands
1996 - 2019 85.21 14.79 Austria 1996 - 2019 84.68 15.32 Sweden 1950 - 2020
90.37 9.63 Saudi Arabia 2002 - 2009 49.80 50.20 Iceland 2000 - 2014 83.21
16.79 Belgium 1996 - 2019 86.56 13.44 Germany 1990 - 2020 82.79 17.21 New
Zealand 1999 - 2017 81.38 18.62 France 1950 - 2019 89.54 10.46 Ireland 1996 -
2019 85.70 14.30 Finland 1996 - 2020 81.46 18.54 Australia 1960 - 2019 85.65
14.35 Canada 1972 - 2020 89.07 10.93 Cyprus 1996 - 2019 70.97 29.03 Israel
2000 - 2019 71.10 28.90 Malta 1996 - 2019 85.52 14.48 UK 1970 - 2019 73.17
26.83 Spain 1995 - 2019 81.21 18.79 Japan 1980 - 2017 83.49 16.51 Aruba 1996 -
2000 89.90 10.10 Venezuela 1997 - 2007 70.97 29.03 Slovenia 1996 - 2019 92.18
7.82 Lithuania 1996 - 2019 65.61 34.39 Czech Rep. 1994 - 2019 81.58 18.42
Italy 1980 - 2020 73.81 26.19 Estonia 1996 - 2019 83.33 16.67 Poland 1996 -
2019 79.39 20.61 Hungary 1996 - 2019 85.79 14.21 Portugal 1996 - 2020 77.11
22.89 Country Period Labor Share Capital Share Slovakia 1996 - 2020 73.14
26.86 Romania 1996 - 2019 65.39 34.61 Greece 1996 - 2019 67.18 32.82 Turkey
2009 - 2017 65.03 34.97 Croatia 1996 - 2019 88.42 11.58 Latvia 1996 - 2019
70.34 29.66 Russia 2008 - 2018 83.98 16.02 Uruguay 2012 - 2016 71.01 28.99
Kazakhstan 1996 - 2019 76.92 23.08 Malaysia 2006 - 2015 62.18 37.82 Chile 1996
- 2009 72.38 27.62 Costa Rica 2012 - 2017 72.33 27.67 Bulgaria 1996 - 2017
67.82 32.18 Mexico 1996 - 2019 56.91 43.09 Domin. Rep. 1996 - 2016 63.64 36.36
Serbia 1997 - 2019 71.66 28.34 Brazil 1996 - 2018 70.99 29.01 Colombia 1996 -
2019 64.40 35.60 Azerbaijan 1996 - 2010 63.15 36.85 Iran 1996 - 2018 74.00
26.00 Peru 2007 - 2019 68.15 31.85 Egypt 1996 - 2015 51.09 48.91 Moldova 1998
- 2019 67.91 32.09 Mongolia 2005 - 2019 73.33 26.67 Botswana 1996 - 2000 51.74
48.26 Guatemala 2001 - 2019 52.58 47.42 Uzbekistan 2015 - 2017 74.63 25.37
Cape Verde 2007 - 2017 78.99 21.01 Nicaragua 2006 - 2018 66.02 33.98 Honduras
2000 - 2015 67.63 32.37 Kyrgyzstan 1996 - 2019 64.66 35.34 Côte d’Ivoire 1996
- 2000 76.68 23.32 Senegal 2014 - 2015 68.10 31.90 Niger 1996 - 2019 38.32
61.68 a Derived at $\frac{\Pi}{C}$ and $\frac{C-\Pi}{C}$ respectively.

Section 2 of the web appendix888See footnote 7 for web address. to this paper
shows line charts plotting labor and capital shares in consumption for each
country over time. Note that only 10 of the 34 countries ranked in the top
half of these 68 countries by net income per capita, but 32 of the 34 ranked
in the lower half, show labor shares in consumption less than 80%. We suspect
that relative undermeasurement of pay, as suggested in explaining high
findings for $f(K)$ in less advanced economies in Table 1, may account for at
least part of this difference (See Appendix A.4.).

##  8 Rate of return in the stationary state

The ratio $C(K)/K$ can be shown as $c(K)$. Then Eq. (9) can be divided by $K$
to show

$f(K)=c(K)\quad\text{collectively.}$ (12)

Substitution in Eq. (6) now gives

$r(K)=g(K)+c(K)\quad\text{collectively.}$ (13)

The classical concept of the stationary state means a state where $g(K)$ holds
at zero. Eq. (13) shows that collective return in the stationary state holds
at $c(K)$, the ratio of consumption drawn from capital to capital. Thus

$r(K)=c(K)\quad\text{collectively, in the stationary state.}$ (14)

A qualifier is appropriate here. By rate of return, we and Jordà et al. mean
unleveraged return on capital owned free and clear of all debt. Some
economists have meant return in an opposite sense as return to an imaginary
entrepreneur who has invested no equity, but rather has borrowed the whole
cost of capital employed. Rate of return to such an entrepreneur would tend to
equal zero in the stationary state. Eq. (14) applies only in the deleveraged
sense used by us, by Jordà et al., and by most investors.

##  9 Data Sources

The data for this study are sourced from the World Inequality Database (WID),
accessible online at WID.world. This comprehensive database integrates data
from national accounts and taxation records across 105 countries,
standardizing these to align with the United Nations System of National
Accounts (SNA) guidelines. For our analysis, we specifically utilize data from
the subset of countries that provide comprehensive reports on net income,
labor income share, consumption and market-value capital.

The primary variables utilized from WID include net national income (mnninc),
the share of labor income (wlabsh), which is represented as $\Pi$, and
purchasing power parity conversions of local currencies to USD (xlcusp). WID
estimates labor share as compensation of employees plus 70% of mixed income.
Additionally, we use average national income per capita (anninc) for arranging
countries in our tables and graphs, and market-value national wealth (mnweal)
as a representation of wealth $K$. Consumption $C$, as referenced in our
analyses, aggregates Government Final Consumption Expenditure (mcongo) and
Private Expenditures of Households and NPISH (mconhn), encapsulating both
government expenditure (GCE) and personal consumption expenditure (PCE).

Data concerning the rates of returns are directly sourced from the tables
presented in Jordà et al. (2019).

##  10 Discussion and conclusions

Our method of deriving cash flow rate, capital growth rate, and rate of return
from data in national accounts serves as a complement to the traditional
method relying on direct research, as in the work of Jordà et al. (2019). Our
results tallied broadly with theirs in the 16 countries for which they report.

The labor and capital shares shown in Table 2 are broadly consistent with
expectations. The labor share of net product is usually estimated at 70% or
so. As consumption is smaller than net product whenever capital growth is
positive, the labor share of consumption applies the same numerator (pay) to a
smaller denominator (consumption versus net product). Also see Appendix A.4..

## Appendix A. Output in terms of human capital

### Appendix A.1. Necessary adjustments to the doctrine $Y=\Delta K+C$

We take value $V$ to mean the sum of human capital $H$ and physical capital
$K$. Then

$V=H+K\ .$ (A.1)

Value growth $\Delta V$ equals value creation $Y$, also called net output or
value added, plus investment from outside, less value yielded out. Value
yielded out less value invested in from outside is cash flow $F(V)$, which may
include cash flow from human capital (see Eq. A.3 below). Then

$\Delta V=Y-F(V)\ ,\quad\text{from which}\quad Y=\Delta V+F(V)\ .$ (A.2)

The sum $\Delta V+F(V)$, or value growth plus cash flow, is called total
return. Thus Eq. (A.2) shows the equivalence of value added and total return.
Although human capital is usually measured at cumulative cost, we follow Petty
(1664), Fisher (1907, 1930), Mincer (1958) and others in measuring it at
present value, or expected lifetime cash flow discounted at the time
preference rate. The flow discounted is pay, including imputed pay, as the
flow yielded out by human capital to its human owner, less invested
consumption999The terms invested consumption, pure consumption, self-invested
work and human depreciation were introduced in Schultz (1961). as the flow
invested in human capital from outside itself. If invested consumption shows
as $C_{s}$, then, we find

$F(H)=\Pi-C_{s}\ ,$ (A.3)

where $F(H)$ is cash flow of human capital. We also follow Ben-Porath (1967)
in describing the growth of human capital as the sum of invested consumption
plus self-invested work less human depreciation. If the last two flows show as
$W_{s}$ and $D(H)$, we have101010Our terms and notations need not follow those
of our sources. We are not aware, for example, that others have used the term
“human cash flow”.

$\Delta H=C_{s}+W_{s}-D(H)\ .$ (A.4)

Eq. (A.2) can be applied to human capital separately to find work $W$, the net
output produced by human capital, as

$\displaystyle W=\Delta H+F(H)=C_{s}+W_{s}-D(H)+\Pi-C_{s}=\Pi+W_{s}-D(H)\ ,$
(A.5)

showing that work exceeds pay by the amount $W_{s}-D(H)$.

Cash flow from assets generally can be spent on investment in other assets or
on consumption. Cash flow from all physical capital of all households together
is net of all investment in physical capital, and is spent wholly on
consumption (Eq. (9)). Cash flow from physical and human capital together is
net of all investment in human capital through invested consumption, and so is
spent on pure consumption alone111111The flow which maintains rather than
creates human capital, and so is exhausted from the economy in satisfying
tastes for vitality and survival.. If this is notated $C_{p}$, then preceding
equations allow

$\displaystyle Y$ $\displaystyle=\Delta V+C_{p}$ (A.6) $\displaystyle=\Delta
H+\Delta K+C_{p}$ $\displaystyle=C_{s}+W_{s}-D(H)+\Delta
K+C_{p}\quad\text{collectively.}$

As consumption $C$ is either invested in human capital or yielded from the
economy in satisfying tastes, we may write

$C_{s}+C_{p}=C\ .\quad\text{From Eq. \eqref{eq-a6}, then}$ (A.7)
$Y=C+W_{s}-D(H)+\Delta K\ ,\quad\text{collectively,}$ (A.8)

showing that $Y$ exceeds the sum $C+\Delta K$ by the amount $W_{s}-D(H)$ at
the collective scale.

### Appendix A.2. Confirmation of Eq. (1)

Net output can also be found at the sum of net outputs of the factors. That
is,

$Y=W+P\ .\quad\text{Substitution of Eq. \eqref{eq-a5} in this gives}$ (A.9)
$Y=\Pi+W_{s}-D(H)+P\ .\quad\text{This and Eq. \eqref{eq-a8}, with
rearrangement, allow}$ (A.10) $Y=C+\Delta K+W_{s}-D(H)=\Pi+P+W_{s}-D(H)\
,\quad\text{collectively,}$ (A.11)

from which

$C+\Delta K=\Pi+P\ ,\quad\text{collectively,}$ (A.12)

confirming Eq. (1) at the collective scale.

### Appendix A.3. Factor shares in net income

We have treated the larger factor as human capital measured as present value
of prospective lifetime pay less invested consumption, rather than as labor
measured, say, in headcount or hours or current pay, but we will follow custom
by terming its share in net income or in consumption as the labor share. Eq.
(A.11) shows that net output, and equivalently net income, exceed the sums
$C+\Delta K$, or equivalently $\Pi+P$, by the difference $W_{s}-D(H)$.

The labor share in consumption was reasoned in Eqs. (10) and (11) to equal
$\Pi$ absolutely, and $\Pi/C$ as a percentage. The capital share in
consumption was reasoned to equal $F(K)$ absolutely, and $F(K)/C$ as a
percentage, since this is the only remaining possible disposition of $C$. Thus
measurement of the factor shares in $C$ would require only measurement of
$\Pi$ and $C$, where measurement of $F(K)$ is inferred as their difference
rather than measured separately.

Factor shares in net income, however, would be defined as $W$ and $P$
absolutely and as $W/Y$ and $P/Y$ as percentages. Equations above allow

$\frac{W}{Y}=\frac{\Pi+W_{s}-D(H)}{C+\Delta
K+W_{s}-D(H)}\quad\text{and}\quad\frac{P}{Y}=\frac{\Delta K+C-\Pi}{C+\Delta
K+W_{s}-D(H)}\ ,$ (A.13)

showing that measurement of factor shares in net income require measurement of
$K,\ W_{s}$ and $D(H)$ as well as $C$ and $\Pi$. As human capital leaves
little market record apart from its rise in pay, and from an important
component of invested consumption in the cost of schooling, measurements of
$Y$ and $W$ are problematic. We suggest that distribution theory might
consider a refocus from shares in net income to shares in consumption on these
grounds of measurability and others. It can be argued that individuals vary in
appetite for wealth, and in talent for managing wealth, and in willingness to
accept the responsibility for that management, but that all agree in need for
consumption. Benevolence and wisdom of public policy, if so, might better be
measured by ample sufficiency of consumption than by equality of net income,
so that individual tastes and talents can account for the rest of differences
in net income.

### Appendix A.4. Summary and discussion

The bad news is that we cannot measure net output in the sense of value added
to both factors. The sums we had thought measured it, meaning $C+\Delta K$ or
equivalently $\Pi+P$, both overlook $W_{s}$ and $D(H)$ as positive and
negative components in value added to human capital. Neither $W_{s}$ (self-
invested work) nor $D(H)$ (human depreciation) is practical to measure.

The good news is that we can measure net profit at market value, including
imputed net profit on the housing and government sectors, subject to
measurement problems discussed earlier, from data for consumption, pay and
market value capital alone. Meanwhile we can measure the labor and capital
shares in consumption as $\Pi/C$ and $(C-\Pi)/C$, with the same qualifier,
without need for other data. This convenience gives some compensation for the
difficulty in finding net income, or factor shares in net income as adjusted
for $W_{s}$ and $D(H)$.

As arguments determining $f(K),g(K),r(K),c(\Pi)$ and $c(K)$ relied on no
simplifying assumptions, we could find those values exactly if we could
measure $C,\Pi$ and $K$ exactly. We cannot. It may be possible, however, to
focus more attention on measuring them, and in seeking consensus as to
appropriate imputations where measurement is impractical. It may be the case,
for example, that pay tends to be less fully recognized in agrarian economies,
where higher percentages of people are self-employed homemakers, farmers and
herders. Capital too might tend to be relatively underreported if such people
are more likely to make implements for use by themselves and their
communities, and to construct their own housing, rather than buy and sell
those things in formal markets. If any such differences in completeness of
measurement were verified, and quantified in standard imputations, it might be
possible to reach more realistic versions of Tables 1 and 2 for those
countries and others.

## References

  * Ben-Porath (1967) Ben-Porath, Y., 1967.  The production of human capital and the life cycle of earnings.  Journal of Political Economy 75, 352–365.
  * Fisher (1907) Fisher, I., 1907.  The Rate of Interest.  McMaster University Archive for the History of Economic Thought.
  * Fisher (1930) Fisher, I., 1930.  The Theory of Interest, as Determined by Impatience to Spend Income and Opportunity to Invest It.  Macmillan, New York.
  * Homer and Sylla (2005) Homer, S., Sylla, R., 2005.  A History of Interst Rates.  Wiley finance, Wiley.  Tex.lccn: 2005281578.
  * Jordà et al. (2019) Jordà, O., Knoll, K., Kuvshinov, D., Schularick, M., Taylor, A.M., 2019.  The rate of return on everything, 1870–2015.  The Quarterly Journal of Economics 134, 1225–1298.  doi:10.1093/qje/qjz012.
  * Mincer (1958) Mincer, J., 1958.  Investment in human capital and personal income distribution.  Journal of Political Economy 66, 281–302.
  * Petty (1664) Petty, W., 1664.  Verbum Sapienti.
  * Schultz (1961) Schultz, T.W., 1961.  Investment in human capital.  The American Economic Review 51, 1–17.