Revisiting the Kepler problem with linear drag using the blowup method and normal form theory

TL;DR

This paper analyzes the Kepler problem with linear drag, revealing invariant manifolds and stable sets using blowup and normal form theory, and providing a geometric understanding of the dynamics with dissipation.

Contribution

It introduces a novel application of blowup and normal form theory to identify invariant manifolds and stable sets in the dissipative Kepler problem, including a zero eccentricity manifold.

Findings

Invariant manifolds are smooth and nonhyperbolic.

A zero limiting eccentricity manifold is characterized as a center manifold.

The blowup analysis offers a geometric view of the dissipative dynamics.

Abstract

In this paper, we revisit the Kepler problem with linear drag. With dissipation, the energy and the angular momentum are both decreasing, but in \cite{margheri2017a} it was shown that the eccentricity vector has a well-defined limit in the case of linear drag. This limiting eccentricity vector defines a conserved quantity, and in the present paper, we prove that the corresponding invariant sets are smooth manifolds. These results rely on normal form theory and a blowup transformation, which reveals that the invariant manifolds are (nonhyperbolic) stable sets of (limiting) periodic orbits. Moreover, we identify a separate invariant manifold which corresponds to a zero limiting eccentricity vector. This manifold is obtained as a generalized center manifold over the zero eigenspace of a zero-Hopf point. Finally, we present a detailed blowup analysis, which provides a geometric picture of the dynamics. %We will use this to shed light on a degenerate case of the limiting eccentricity vector. %The aforementioned invariant manifolds only make up a subset hereof. We believe that our approach and results will have general interest in problems with blowup dynamics.

Equations

$$\ddot{u}+\delta \dot{u}+c\frac{u}{|u{|}^3}=0,$$

$$E(u,\dot{u})=K(\dot{u})+P(u),K(\dot{u}):=\frac{1}{2}|\dot{u}{|}^2,P(u):=-\frac{1}{|u|},$$

$$L(u,\dot{u})=u\wedge \dot{u},$$

$$\mathcal{E}(u,\dot{u})=\dot{u}\wedge L(u,\dot{u})-\frac{u}{|u|},$$

$$\frac{d}{dt}L=-\delta L,$$

$$L(u(t),\dot{u}(t))=e^{-\delta t}L(u(0),\dot{u}(0)),$$

$${\mathcal{E}}_{\infty }(u,\dot{u})=\underset{t\to t_{\text{max}}^{-}}{\lim }\mathcal{E}({\varphi }_t(u,\dot{u})),$$

$$u:=r{e}^{i\theta },$$

$$l:=|L|=r^2\dot{\theta }\ge 0,$$

$$\ddot{r}$$

$$\dot{l}$$

$$r=l^2{r}_1,v=l\dot{r}$$

$$K(\dot{u})$$

$${\dot{r}}_1$$

$$\dot{v}$$

$$\dot{l}$$

$$r_1^{\prime }$$

$$v^{\prime }$$

$$l^{\prime }$$

$$r_1^{\prime }$$

$$v^{\prime }$$

$$H(r_1,v)=\frac{1}{2}{v}^2+\frac{(r_1-1{)}^2}{2r_1^2}.$$

$$H(l^{-2}r,l\dot{r})=\frac{1}{2}{l}^2{\dot{r}}^2+\frac{(l^2-r{)}^2}{2r^2},$$

$$\dot{u}=\frac{v}{l}{e}^{i\theta }+\frac{il}{r}{e}^{i\theta },$$

$$e^{-i\theta }\mathcal{E}(u,\dot{u})=\frac{1-r_1}{r_1}-iv,$$

$$H(r_1,v)=\frac{1}{2}|\mathcal{E}(u,\dot{u}){|}^2.$$

$$H_{\infty }(r_1,v,l):=\underset{t\to \infty }{\lim }H({\underline{r}}_1(t,r_1,v,l),\underline{v}(t,r_1,v,l))=\frac{1}{2}|{\mathcal{E}}_{\infty }(u,\dot{u}){|}^2,$$

$$\underline{z}(0,\dots )=z,$$

$$r_{1,-}(h)=\frac{1}{1+\sqrt{2h}},r_{1,+}(h)=\frac{1}{1-\sqrt{2h}}.$$

$$q_1:(r_1,v,l)=(1,0,0),$$

$$r_1$$

$$\frac{d^{2n}}{d{x}^{2n}}{F}_1(0)=0,\frac{d^{2n+1}}{d{x}^{2n+1}}{G}_1(0)=0,$$

$$r=l^2(1+l^3{F}_1(l^3)),\dot{r}=l^2{G}_1(l^3),$$

Full text

keywords. Invariant manifolds, nonhyperbolic sets, dynamical systems theory, blowup.

1. Introduction

In this paper, we consider the Kepler problem with linear drag [34, 35]

[TABLE]

with $u(t)\in\mathbb{R}^{3}\backslash\{0\}$, $c>0$ and for $\delta>0$. The singularity at $u=0$ corresponds to the collision limit and for $\delta=0$ (no drag/damping), we obtain the classical Kepler problem, whose orbits are conic sections. In fact, as in the classical case, we may scale $u$ and $t$ to achieve $c=1$ so we will assume this henceforth. (In this way, $\delta$ is replaced by $\delta c^{-1/2}$). For $\delta=0$, the energy:

[TABLE]

the angular momentum:

[TABLE]

and the eccentricity vector

[TABLE]

are all conserved quantities. For $\delta\neq 0$, we have $\frac{d}{dt}E=-\delta|\dot{u}|^{2}\leq 0$ and hence the energy is monotonically decreasing for $\delta>0$. Moreover, a simple calculation shows that

[TABLE]

so that

[TABLE]

and the angular momentum is therefore exponentially decreasing. Notice, however, that the direction of $L$ is constant and the motion $(u,\dot{u})$ is therefore contained in a plane (the orbital plane).

$\mathcal{E}$ is also not conserved for $\delta\neq 0$, but in [35] it was shown that there exists a limiting eccentricity vector: Let $\phi_{t}(u,\dot{u})$ denote flow associated with (1.1). Then

[TABLE]

exists for all $u,\dot{u}$. Here $t_{\text{max}}$ is the maximum time of existence (finite for $L=0$, infinite for $L\neq 0$, see [34]). In [35], it was shown that the components of $\mathcal{E}_{\infty}$ are functional independent and rotationally equivariant: $\mathcal{E}_{\infty}(Ru,R\dot{u})=R\mathcal{E}_{\infty}(u,\dot{u})$ for all $R\in SO(3)$, and that $|\mathcal{E}_{\infty}(u,\dot{u})|$ attains all values in $(0,1]$. If $|\mathcal{E}_{\infty}(u,\dot{u})|=1$ then $L=0$, see [35]. In [36], the same authors generalized their results to the case, where $\delta$ is a function $\delta(|u|)$ of $|u|$ satisfying $\delta(|u|)\geq c>0$. In [32], the drag $\delta(|u|)$ was singular at $u=0$, and they showed that in some cases, the limiting eccentricity vector can be discontinuous.

The question of smoothness of $\mathcal{E}_{\infty}$ was left open in [35]. In this paper, we will give a different characterization of $\mathcal{E}_{\infty}$, which will allow us to address the issue of smoothness for $|\mathcal{E}_{\infty}|\in(0,1)$. At the same time, using dynamical systems theory, we identify a new smooth invariant manifold corresponding to $\mathcal{E}_{\infty}=0$, which acts as a center of the oscillating orbits with $|\mathcal{E}_{\infty}|\in(0,1)$. It should be said that this invariant manifold, which will be one-dimensional in a reduced phase space, corresponding to $\mathcal{E}_{\infty}=0$, was actually derived as a formal series in [19], which sparked the interest of the present author. Essentially, our proof shows that this series is summable in the sense of Borel-Laplace. Separately, using blowup and compactification as our main tools, we provide a geometric description of the dynamics. This will shed further light on $|\mathcal{E}_{\infty}|\rightarrow 1$.

The study of dissipation in celestial mechanics has a long history. It even dates back to Jacobi [20], who introduced dissipative forces of the type $-\delta|\dot{u}|^{n-1}\dot{u}$; the case $n=1$ corresponds to the linear drag studied in the present paper. There are at least two important mechanisms for dissipation in celestial mechanics: particle collisions due to nebula (Stokes’ dissipation) and solar radiation (Poynting-Robertson dissipation), see [6, 34].

Corne and Rouche [8] and Diacu [9] were perhaps the first contributors towards the development of a qualitative theory of the Kepler problem with drag (1.1) for general families of non-constant drag forces $\delta=\delta(u,\dot{u})$, depending on $u$ and $\dot{u}$. [8] considered (1.1) with $\delta(u,\dot{u})=k(|\dot{u}|)/|\dot{u}|$ and showed, under some additional assumptions on the function $k$, that all solutions go to the singularity (potentially in finite time). On the other hand, [9] analyzed the qualitative dynamics of the dissipative Kepler problem within a generalized class of Stokes drag; this family includes the important Poynting-Robertson case. Many years later, Margheri, Ortega, and Rebelo in [34] studied the Kepler problem with linear drag (1.1) and provided a more thorough description of the dynamics. In particular, they proved that the system is complete, i.e. solutions exist globally in time, on the set of nonzero angular momentum. Later in [35], the same authors then provided a more geometric description, including the properties of the limiting eccentricity vector. Their results also showed that $\limsup_{t\rightarrow\infty}$ and $\liminf_{t\rightarrow\infty}$ of $|u(t)||L(t)|^{-2}$ both exist along orbits with $|\mathcal{E}_{\infty}|\in(0,1)$, see also [19].

In parallel, there has been some studies of dissipative versions of the restricted three body problem, see e.g. [6, 18, 33]. In [6] the authors used numerical methods (based upon Fast Lyapunov Indicators (FLI)) to provide information on the different regions of the phase space. They demonstrated both collision and non-collision trajectories. Interestingly, they also documented periodic orbit attractors, but only in the case of linear and Stokes drags. In contrast, in the case of the Poynting-Robertson dissipation, the authors found no other attractors beside the primaries (collisions). Following on from this research, [33] studied the existence (and nonexistence) of periodic orbits, including Hopf bifurcations around the libration points $L_{4}$ and $L_{5}$.

In this paper, we will use the blowup method and normal form theory to study the dynamics of (1.1). In the context of celestial mechanics and Hamiltonian systems, normal form theory has a long history, dating back to the work of Poincaré and Birkhoff. In fact, KAM theory [1] itself may be viewed as a normal form theory. On the other hand, the blowup method provides a general framework for studying degenerate equilibria in local dynamical systems theory, where the hyperbolic theory (e.g. Hartman-Grobman and center manifolds) does not apply. The rough idea of this approach is to apply a non-invertible transformation (like polar coordinates) that blows up the equilibrium to a sphere (or cylinders in case of lines of degenerate equilibria). By appropriately choosing weights associated to the transformation, it is possible (at least for analytic systems) to divide the resulting vector-field by a power of the radius, measuring the distance to the equilibrium. This gives rise to a new vector-field, only equivalent to the original one away from the equilibrium, for which hyperbolicity (or ellipticity) may be (partially) gained on the blowup of the singularity. Sometimes this approach of blowing up equilibria has to be used successively, see e.g. [14].

In recent years, this blowup approach has gained prominence within the area of singular perturbation theory, because here degenerate equilibria occur naturally, see [12, 13, 15, 30]. In combination, Fenichel’s geometric singular perturbation theory and the blowup method has been very successful in describing global phenomena in slow-fast models, see e.g. [23, 24, 25, 26, 27]. More recently, the blowup approach has been generalized with the purpose of “gaining smoothness”, rather than hyperbolicity, in the context of smooth systems approaching nonsmooth ones, see [21, 22, 29, 31].

Obviously, blowup can also take on a different meaning in mathematics, namely (finite time) blowup of solutions of (ordinary or partial) differential equations. In dynamical systems theory, blowup solutions can be studied by Poincaré or Poincaré-Lyapunov compactification (which in fact bear some resemblance to the blowup method), see [14]. Upon compactification one can study equilibria, again after proper desingularization of the vector-field, at infinity and these can be analyzed by local methods of dynamical systems theory. In particular, such points at infinity are sometimes completely degenerate, which can then be resolved by the blowup method.

Related to blowup of solutions is the existence of solutions approaching true singularities of ordinary differential equations; like the collision limit $u=0$ for (1.1) where the associated vector-field is ill-defined. The analysis of collision (as well as near-collision) solutions in the $n$-body problem has a long history, also dating back to Poincaré, see also [5, 10, 11]. For the two-body problem, which can be reduced to (1.1)δ=0, there is a change of coordinates (the Levi-Cevita transformation) and a nonlinear transformation of time (related to desingularization), that transforms the collision into a regular point of the equations. The Levi-Cevita transformation is – similar to a blowup transformation – not invertible; it is in fact a double cover. More generally in the $n$-body problem, the Levi-Cevita transformation can be applied to show that solutions of the $n$-body problem can be analytically continued through isolated binary collisions (this is also known as regularization in the context of celestial mechanics). Interestingly, Mcgehee [38] used the blowup method to study the more complicated triple collision in the context of the collinear three-body problem. Indeed, the author blew up the collision set to a sphere and upon applying desingularization, he obtained hyperbolic equilibria points on a collision manifold. This led the author (through hyperbolic invariant manifolds) to conclude that the triple collision cannot be regularized. Subsequently, this approach was used in [16] to show that the simultaneous binary collision scenario was $C^{0}$-regularizable. Interestingly, [37] proved that it is exactly $C^{8/3}$-regularizable in the collinear case. A separate geometric prove – based upon normal form theory and the blowup approach of [16, 38] – was given in [10]. This approach led the same authors in [11] to prove that the same results hold in the planar case, a result that was initially conjectured by [37].

1.1. Outline

The paper is organized as follows: In Section 2, we lay out our approach (based upon certain blowup transformations) for characterizing $\mathcal{E}_{\infty}$ and present two theorems on the existence of smooth invariant manifolds of (1.1): Theorem 2.6 for the existence of the orbit corresponding to $\mathcal{E}_{\infty}=0$ and Theorem 2.7 for the existence of a smooth invariant manifold corresponding to $|\mathcal{E}_{\infty}|\in(0,1)$, see also Theorem 2.8. We prove Theorem 2.6 in Section 3 using general results on Gevrey-1 invariant manifolds $y=Y(x)$ for analytic systems of the form $x^{2}\frac{dy}{dx}=F(x,y)$, with $F(0,0)=0$ and $D_{y}F(0,0)$ non-singular. It is author’s impression that these results are not so well-known. We follow [4], which proves the existence of such manifolds (and certain normal forms) in perhaps the most accessible way. Theorem 2.7 is proven in Section 4 using normal form theory (based upon averaging) to set up an appropriately regular equation for the invariant manifolds that can be solved upon application of the implicit function theorem. This approach may have general interest. In Section 5, we apply a sequence of blowup transformations along with an appropriate compactification in order to provide a geometric description of the dynamics. We summarize this in Fig. 1. From [35], it is known that $|\mathcal{E}_{\infty}|$ cannot exceed $1$. The results of our blowup analysis, will provide a different characterization of this fact that also allow us to address the subtleties of $|\mathcal{E}_{\infty}|\rightarrow 1$. We lay this out in further details in our final discussion section, Section 6.

2. Existence of invariant manifolds

Since the direction of $L$ is preserved, it is without loss of generality to consider $u(t)\in\mathbb{R}^{2}$. We will do so henceforth. Upon identifying $\mathbb{R}^{2}$ with $\mathbb{C}$ in the usual way, we put

[TABLE]

and let

[TABLE]

denote the magnitude of the angular momentum. Then (1.1) becomes

[TABLE]

For the purpose of this section, we suppose that $l>0$, and define the coordinates $r_{1}$ and $v$ by

[TABLE]

These coordinates appear in a systematic way through our blowup approach, see Section 5. Then the kinetic energy, $K$, takes the following form:

[TABLE]

in the $(r_{1},v,l)$-coordinates. Moreover, we have the following differential equations

[TABLE]

Now, since $l^{-3}$ is a common factor of the right hand side, we make a nonlinear transformation of time that corresponds to multiplication by $l^{3}$:

[TABLE]

(2.5) and (2.6) are equivalent for $l>0$. However, (2.6) is defined for $l=0$ (it defines an invariant set) and we can therefore use dynamical systems theory to infer properties from $l=0$ to $l>0$ by working on (2.6). In turn, this then carries over to the equivalent system (2.5).

Remark 2.1**.**

In the following we will reserve $\dot{()}$ to denote differentiation with respect to the original time. In comparison, we will use $()^{\prime}$ repeatedly to refer to differentiation with respect to different times. It should be clear form the context how different $()^{\prime}$ may be related.

Setting $l=0$ in (2.6) gives

[TABLE]

which is Hamiltonian with Hamiltonian function:

[TABLE]

Remark 2.2**.**

If $\delta=0$, then (1.1) reduces to the conservative Kepler problem and the Hamiltonian function (2.8), written in terms of the original variables $r,\dot{r}$ and $l$:

[TABLE]

*is, along with $l$, a conserved quantity of (2.3). In fact, $H=\frac{1}{2}|\mathcal{E}|^{2}$ in this case. *

Lemma 2.3**.**

Consider (2.1), (2.2) and (2.4). Then

[TABLE]

and

[TABLE]

so that

[TABLE]

Proof.

Direct calculation. ∎

As a corollary, we have that

[TABLE]

where we use $(\underline{r}_{1}(t,r_{1},v,l),\underline{v}(t,r_{1},v,l),\underline{l}(t,r_{1},v,l))$ with

[TABLE]

$z=r_{1},v,$ and $l$, as our notation for the flow of (2.6).

Lemma 2.4**.**

$(r_{1},v)=(1,0)$* is a center for (2.7), surrounded by periodic orbits $\Gamma_{1}(h)$, $h\in\left(0,\frac{1}{2}\right)$, given by the level sets $H(r_{1},v)=h$. $\Gamma_{1}(h)$ intersects the $r_{1}$-axis in two points $(r_{1,\pm}(h),0)$ with*

[TABLE]

The unbounded orbit $\Gamma_{1}\!\left(\frac{1}{2}\right)$ given by the level set $H(r_{1},v)=\frac{1}{2}$ (corresponding to $|\mathcal{E}_{\infty}|=1$ by (2.10)), intersects the $r_{1}$-axis once in $r_{1}=r_{1,-}(\frac{1}{2})=\frac{1}{2}$, and is a separatrix, separating the bounded orbits ($H(r_{1},v)<\frac{1}{2}$) from unbounded orbits ($H(r_{1},v)>\frac{1}{2}$).

Proof.

Direct calculation. Notice in particular that $(r_{1},v)=(1,0)$ is an extremum (minimum) of the Hamiltonian function and therefore also a center of (2.7). In fact, the linearization around $(r_{1},v)=(1,0)$ produces $\pm i$ as the eigenvalues. ∎

We illustrate the phase portrait of (2.7) in Fig. 2.

Linearization of (2.6) at the corresponding equilibrium point $q_{1}$ of the full system, given by

[TABLE]

clearly produces eigenvalues $\pm i,0$. It therefore corresponds to a zero-Hopf point [17]. Similarly, all periodic orbits $\Gamma_{1}(h)$, $h\in\left(0,\frac{1}{2}\right)$ are also degenerate when embedded within the full system (2.6). The description of the stable sets of these sets of points is therefore complicated by the lack of hyperbolicity.

Remark 2.5**.**

For $\delta=0$, recall Remark 2.2, $q_{1}$ corresponds to circular orbits of the (conservative) Kepler problem with zero eccentricity, whereas $\Gamma_{1}(h)$, $h\in(0,\frac{1}{2})$, correspond to elliptic ones with $|\mathcal{E}|\in(0,1)$. Finally, $\Gamma_{1}\!\left(\frac{1}{2}\right)$ corresponds to the parabolic orbit whereas $\Gamma_{1}(h)$, $h>\frac{1}{2}$ corresponds to hyperbolic ones.

2.1. Main results

We now state our main results.

Theorem 2.6**.**

(The local stable manifold $W_{loc}^{s}(q_{1})$ of $q_{1}$.)* The local stable set $W_{loc}^{s}(q_{1})$ is a smooth one-dimensional manifold, taking the following graph form*

[TABLE]

for $l_{0}>0$ small enough and where $F_{1},G_{1}:[0,l_{0}]\rightarrow\mathbb{R}$ are Gevrey-1 smooth functions. Moreover,

[TABLE]

for all $n\in\mathbb{N}_{0}$.

Upon using (2.4), we find that $W^{s}(q_{1})$ in (2.11) takes the following form

[TABLE]

with respect to $r$ and $\dot{r}$. These quantities are smooth functions of $l$ along this orbit and decay like $e^{-2\delta t}$. Since $F_{1}(0)=0$, see (2.12) with $n=0$, we obtain from (1.2), (2.2) and (2.4) that

[TABLE]

and consequently that the sum

[TABLE]

being equal to

[TABLE]

by (2.13), is a smooth function of $l^{3}$ along the orbit in (2.11). This result complements results of [35, Eq. (43)] which showed that $\limsup_{t\rightarrow\infty}$ and $\liminf_{t\rightarrow\infty}$ of $\dot{\theta}l^{-3}$ both exist whenever $|\mathcal{E}_{\infty}|\in(0,1)$.

We now turn our attention to the stable sets of $W^{s}(\Gamma_{1}(h))$ with $h\in(0,\frac{1}{2})$.

Theorem 2.7**.**

(The local stable manifold of $\Gamma_{1}(h)$.)* Fix $h\in\left(0,\frac{1}{2}\right)$. Then the stable set $W^{s}(\Gamma_{1}(h))$ of $\Gamma_{1}(h)$ is defined by $H_{\infty}(r_{1},v,l)=h$. Moreover, fix any $k\in\mathbb{N}$. Then there exists a neighboorhood $N(k,h)$ of $\Gamma_{1}(h)$ such that $W_{loc}^{s}(\Gamma_{1}(h))=W^{s}(\Gamma_{1}(h))\cap N(k,h)$ is a $C^{k}$-smooth two-dimensional submanifold. The dependency on $h$ is also $C^{k}$-smooth.*

Although the domain $N(k,h)$ depends upon $k$ and $h$, the flow is regular away from $l=0$ and consequently, we can globalize the local manifolds by application of the backward flow. Since the system is real analytic, it follows that these global stable manifolds are in fact $C^{k}$ also. Similarly, if we fix a compact interval $I\subset\left(0,\frac{1}{2}\right)$ then we have a uniform description of all local manifolds within $N(k):=\cap_{h\in I}N(k,h)$ and by working on this set, our approach also shows that the stable manifolds are also $C^{k}$. In turn, seeing that $k\in\mathbb{N}$ is arbitrary, we obtain the following:

Theorem 2.8**.**

The global stable manifolds $W^{s}(q_{1})$ and $W^{s}(\Gamma_{1}(h))$ of $q_{1}:\,(r_{1},v,l)=(1,0,0)$ and $\Gamma_{1}(h)$, $h\in\left(0,\frac{1}{2}\right)$, respectively, are each $C^{\infty}$. The dependency of $W^{s}(\Gamma_{1}(h))$ on $h$ is also $C^{\infty}$.

Due to (2.10), we have obtained a complete description of the smoothness of $|\mathcal{E}_{\infty}(u,\dot{u})|\in(0,1)$ which was left open in [35]. It seems likely that the global stable manifolds of $\Gamma_{1}(h)$ are also Gevrey-1 (as the stable manifold $W^{s}(q_{1})$ of $q_{1}:\,(r_{1},v,l)=(1,0,0)$), but this would require better normal forms. [3] considers such normal forms but the condition regarding the trace is violated in the present context.

Finally, we remark that the existence of the invariant set $W^{s}(q_{1})$ was proven more indirectly in [35] using degree theory. Although this approach does not address the smoothness, the result in [36] holds true for a more general class of dissipations defined by $\delta=\delta(|u|)$, with $\delta(|u|)\geq c>0$. However, our results can be extended to this general case also. We leave this to the interested reader.

3. Proof of Theorem 2.6

To prove Theorem 2.6, we consider (2.6) in terms of

[TABLE]

This gives

[TABLE]

or

[TABLE]

upon eliminating time. To zoom in on $q_{1}$ at $r_{1}=1,v=0$, $x=0$, we perform a blowup transformation:

[TABLE]

leaving $x$ fixed. This gives the final system:

[TABLE]

setting $y=(r_{11},v_{11})$ and where

[TABLE]

The eigenvalues of $A$ are $\pm\frac{i}{\sqrt{3\delta}}$ and $f$ is real analytic. It is standard, that there exists a formal series solution

[TABLE]

of (3.5), see e.g. [4] and also [19] for the details of the expansion in the present case. In particular,

[TABLE]

and by induction on $n$

[TABLE]

for all $n\in\mathbb{N}$. Here $*$ is an unspecified quantity that depends upon $n$ and $\delta$.

Let $S(\phi,r)\subset\mathbb{C}$ be the open sector region in $\mathbb{C}$ centered along the positive real axis, with radius $r$ and opening $\phi$:

[TABLE]

Let $\overline{S}(\phi,r)$ denote its closure. The following result shows that the series (3.6) is Gevrey-1 and 1-summable [2]:

Proposition 3.1**.**

There exists $x_{0}>0$, $\phi>0$ both sufficiently small and a Gevrey-1 function $Y:\overline{S}(\phi+\pi,r)\rightarrow\mathbb{C}$ with $Y(0)=0$, which is real analytic on $S(\phi+\pi,r)$, such that the graph

[TABLE]

solves (3.5) with $Y(0)=0$. Here $Y$ has $\widehat{Y}$, see (3.6), as a Gevrey-1 asymptotic series, such that $Y^{(n)}(0)=n!Y_{n}$

Proof.

The result follows from [4, Theorem 3]; although this result does not address the existence of an invariant manifold directly (instead [4] proves existence of a certain normal form), this can be obtained as a corollary, using the invariance of the set $y=0$ of [4, Eq. (8)]. For completeness, we include a version of the proof that only addresses the existence of $y=Y(x)$, see also [28, App. A], which includes a similar proof.

Following [4] we proceed by first (a) transforming (3.3) into an equation on the “Borel-plane” (through the Borel-transform $\mathcal{B}$), then (b) apply a fixed-point argument there and finally (c) obtain our desired solution by applying the Laplace transform.

For our purposes, the Borel transform is defined in the following way: If $h(x)=\sum_{n=1}^{\infty}h_{n}x^{n}$ is a Gevrey-1 formal series:

[TABLE]

then the Borel transform of $h$ is given by

[TABLE]

Clearly, $\mathcal{B}(h)$ is analytic on $|u|<b^{-1}$ if (3.8) holds true. For the Laplace transform, on the other hand, we need analytic functions $\alpha(u)$ that are at most exponentially growing $|\alpha(u)|\leq\mathcal{O}(1)e^{\zeta|u|}$, $\zeta>0$, in an infinite sector. With this in mind, let $S(\phi)\subset\mathbb{C}$ be the (infinite) sector centered along the positive real $x$-axis with opening $\phi\in(0,\pi)$, and let $B(R)$ be the open ball of radius $R$ centered at [math]. Finally, set

[TABLE]

Then for any $\zeta>0$, we define the norm

[TABLE]

see [4], on the space of analytic functions on $\Delta$:

[TABLE]

The normed space $(\mathcal{G},\|\cdot\|_{\zeta})$ is a complete space. We will need $\zeta$ sufficiently large in the following. The factor $1+\zeta^{2}|w|^{2}$ in the norm $\|\cdot\|_{\zeta}$ ensures that the convolution:

[TABLE]

is continuous as a bilinear operator on $\mathcal{G}$. In particular, we have

[TABLE]

see [4, Proposition 4]. The Laplace transform (along the positive real axis)

[TABLE]

is then well-defined for any $\alpha\in\mathcal{G}$. In fact, we have the following result.

Lemma 3.2**.**

[4, Proposition 3]** The Laplace transform defines a linear continuous mapping, with operator norm $\|\mathcal{L}\|\leq 1$, from $\mathcal{G}$ to the set of analytic functions on a local sector $S(\pi+\phi,R_{0})=S(\pi+\phi)\cap B(R_{0})$ for $\phi\in(0,\pi)$ and $R_{0}>0$ sufficiently small. Moreover,

[TABLE]

and

[TABLE]

with $u\alpha$ being the function $u\mapsto u\alpha(u)$, for every $\alpha,\beta\in\mathcal{G}$.

Following (3.10), we are now led to write the left hand side of (3.5) with $y=Y(x)$ as $[uI-A]\Phi(u)$ with $\Phi=\mathcal{B}(Y)$. To set up the associated right hand side, we need to deal with the nonlinearity $f(x,Y(x))$. This is described in [4, Proposition 5]:

Lemma 3.3**.**

Write $f$ as the convergent series $f(x,y)=\sum_{n=1}^{\infty}f_{n}(x)y^{n}$, $f_{n}(x):=\sum_{m=1}^{\infty}f_{mn}x^{m}$ and let $F_{n}$ be the Borel transform of $f_{n}$. Then $F_{n}\in\mathcal{G}$ for each $n$. Fix $C_{0}>0$ and for large values of $\zeta$, recall (3.9), consider $\alpha\in\mathcal{G}$ with $\|\alpha\|_{\zeta}\leq C_{0}$. Then $\alpha\mapsto f^{*}(\alpha)$ defined by

[TABLE]

which converges in $\mathcal{G}$, is differentiable and satisfies the following estimates

[TABLE]

for some constant $C_{1}>0$ depending only on $f$ and $C_{0}$. Moreover,

[TABLE]

Following (3.12), we are therefore finally led to consider

[TABLE]

where $\Phi$ is the Borel transform of $Y$. The equation (3.13) has the form of a fixed point equation. Since the eigenvalues of $A$ are imaginary, if we take $\phi\in(0,\pi)$ then $[uI-A]^{-1}$ is uniformly bounded on $S(\phi)$. Using (3.11), it therefore follows that there is some $M>0$, depending on $f$, $\phi$ and $R$, such that the right hand side of (3.13) defines a contraction on the subset of $\mathcal{G}$ with $\|\cdot\|_{\zeta}\leq M$ for $\zeta>0$ large enough. Consequently, by Banach’s fixed point theorem there is a unique solution $\Phi\in\mathcal{G},\|\Phi\|_{\zeta}\leq M$, solving (3.13). By applying the Laplace transform, we obtain the desired solution

[TABLE]

of (3.5), using (3.10) and (3.12). The function $Y$ is defined on the domain $S(\pi+\phi,R_{0})=S(\pi+\phi)\cap B(R_{0})$ and has the properties specified by Lemma 3.2. Finally, we emphasize that the Borel transform is real when the argument is. Consequently, (3.13) is real when $\Phi$ is real. Since the Laplace transformation is real upon integrating along the positive real axis, $Y=\mathcal{L}(\Phi)$ is real analytic on $S(\pi+\theta,R_{0})$ as claimed. ∎

Upon transforming the manifold in Proposition 3.1 back to the $(r_{1},v,l)$-coordinates, using (3.1) and

[TABLE]

cf. (3.4), we obtain (2.11). (2.12) then also follows from (3.7). $W_{loc}^{s}(q_{1})$ cannot contain points not in (2.11); this follows from Theorem 2.7, which we prove in the following section. This completes the proof of Theorem 2.6.

4. Proof of Theorem 2.7

To prove Theorem 2.7, we will proceed in three steps: First we introduce appropriate (action-angle) coordinates to parameterize $H(r_{1},v)=h$, see Section 4.1. Subsequently in Section 4.2, we bring our system into a normal form (based upon averaging). Finally in Section 4.3, we use the normal form to set up an equation for the invariant manifold $W^{s}(\Gamma_{1}(h))$, which we solve using the implicit function theorem. Essentially, our approach is reminiscint of a flow-box argument; we will show that there are smooth coordinates $(h,\phi,l)$, $\phi\in\mathbb{T}=\mathbb{R}/(2\pi\mathbb{Z})$, with $h^{\prime}=0$.

4.1. Action-angle coordinates

We first consider the planar Hamiltonian system (2.7), repeated here for convinience:

[TABLE]

with Hamiltonian function:

[TABLE]

The set $H^{-1}(I)$ with $I:=\left(0,\frac{1}{2}\right)$ is filled with periodic orbits centered around the point $(r_{1},v)=(1,0)$. It is standard, see e.g. [39], that there exists a real-analytic symplectic diffeomorphism on this set

[TABLE]

with action-angle coordinates $(A,\phi_{1})$, such that (4.1) becomes

[TABLE]

Here $\Omega(A)=H\circ\Psi^{-1}(A,\phi_{1})$ is the Hamiltonian function expressed in the new coordinates; the important point is obviously that $\Omega$ is independent of $\phi_{1}$. Consider $\Gamma_{1}(h)=\{(r_{1},v)\,:\,H(r_{1},v)=h\}$, $h\in I$, and write $\Psi=(\Psi_{A},\Psi_{\phi})$. Then $A=\Psi_{A}(r_{1},v)$ is just the area of the oval $\{(r_{1},v)\,:\,H(r_{1},v)\leq h\}$, whereas $\phi_{1}=\Psi_{\phi}(r_{1},v)\in\mathbb{T}$ is so that $dr_{1}\wedge dv=d\phi_{1}\wedge dA$. By (4.2), we see that the angle $\phi_{1}$ is a scaling of time ($\operatorname{mod}2\pi$) such that each $\Gamma_{1}(h)$ is $2\pi$-periodic. Consequently, if $p(h)$ is the period of $\Gamma_{1}(h)$ then

[TABLE]

For our purposes, it will be more convinient to use $H=\Omega(A)$ rather than $A$ as an action-variable. Using (4.3), we have that

[TABLE]

is a real-analytic diffeomorphism with the inverse function $H\mapsto\Omega^{-1}(H)$ defined on $H\in\left(0,\frac{1}{2}\right)$. We compose $\Psi_{A}$ with (4.4) and obtain the following real-analytic diffeomorphism:

[TABLE]

on $H^{-1}(I)$, transforming (4.1) into

[TABLE]

with

[TABLE]

We now turn our attention to the $l>0$-system. We will again prefer to work with $x=l^{3}$ rather than $l$ and therefore consider (3.2), repeated here for convinience

[TABLE]

Lemma 4.1**.**

*Let $X_{1}$ denote the vector-field of (4.6) and define the following function *

[TABLE]

Then

[TABLE]

where $\mathcal{L}_{X_{1}}H_{1}=\nabla H_{1}\cdot X_{1}$ is the Lie-derivative.

Proof.

Follows from a direct calculation:

[TABLE]

so that

[TABLE]

∎

Lemma 4.2**.**

Fix an open interval $J$ such that $\overline{J}\subset I$. Then there exists an $x_{0}>0$ such that

[TABLE]

defined on $H^{-1}(J)\times(-x_{0},x_{0})$, is a real analytic diffeomorphism.

Proof.

$\widetilde{\Psi}(r_{1},v,x)=(\widetilde{\Psi}_{1}(r_{1},v,x),\widetilde{\Psi}_{2}(r_{1},v,x),x)$ is an $x$-fibered (polynomial) perturbation of the diffeomorphism defined by $(r_{1},v,x)\mapsto(H(r_{1},v_{1}),\phi_{1}(r_{1},v_{1}),x)$. The result therefore follows by the inverse function theorem. ∎

Lemma 4.3**.**

Let $X_{1}$ denote the vector-field of (4.6) and consider $\widetilde{\Psi}$ from Lemma 4.2. Then $\widetilde{\Psi}_{*}X_{1}$ takes the following form:

[TABLE]

with $\Omega_{0}$ defined in (4.5), for all $H_{1}\in J$. The functions $R_{1}$ and $P_{1}$ are both real-analytic functions defined on $J\times\mathbb{T}\times(-x_{1},x_{1})$ for $x_{1}>0$ sufficiently small.

Proof.

The result follows directly from (4.2), Lemma 4.1 and Lemma 4.2. In particular, by (4.7) we have that

[TABLE]

with $r_{1}=r_{1}(H_{1},\phi_{1},x)$, cf. Lemma 4.2. ∎

4.2. A normal form

We will now use an averaging approach to normalize (4.8). This will consist of pushing the angle-dependency to higher order with respect to $x$.

Lemma 4.4**.**

(The Iterative Lemma)* Consider the real-analytic system*

[TABLE]

defined on $J_{n}\times\mathbb{T}\times(x_{n},x_{n})$ and with $n\in\mathbb{N}$. Here $\Omega_{n}(H_{n},0)=\Omega_{0}(H_{n})\geq c>0$. Then for any $J_{n+1}\subset J_{n}$ there exists a constant $x_{n+1}>0$ and an $x$-fibered, real-analytic diffeomorphism $(H_{n},\phi_{n},x)\mapsto(H_{n+1},\phi_{n+},x)$ of the near-identify form

[TABLE]

such that

[TABLE]

with the right hand side being real-analytic on $J_{n+1}\times\mathbb{T}\times(-x_{n+1},x_{n+1})$.

Proof.

The result follows from a modification of the classical averaging theorem, see e.g. [17]. We write

[TABLE]

where

[TABLE]

is the mean of $\phi\mapsto R_{n}(H_{n},\phi,x)$, and where $\widetilde{R}_{n}$ has zero mean:

[TABLE]

We then take

[TABLE]

in (4.10). Notice that $T_{n}$ is well-defined for $\phi_{n}\in\mathbb{T}$ by (4.13). Moreover,

[TABLE]

We therefore have by (4.10):

[TABLE]

with

[TABLE]

using (4.14) to conclude that $\{\cdots\}=0$.

Subsequently, we define

[TABLE]

such that

[TABLE]

It follows from (4.10) that $P_{n+\frac{1}{2}}$ extends smoothly to $x=0$. We write

[TABLE]

where

[TABLE]

is the mean of $\phi\mapsto P_{n+\frac{1}{2}}(H_{n+1},\phi,x)$, and where $\widetilde{P}_{n+\frac{1}{2}}$ has zero mean:

[TABLE]

Then we apply the same procedure on the $\phi_{n}$-equation. In particular, we consider a transformation defined by

[TABLE]

fixing $H_{n+1}$ and $x$, with

[TABLE]

This leads to the following equation for $\phi_{n+1}$:

[TABLE]

with

[TABLE]

using (4.15) to conclude that $\{\cdots\}=0$. Finally, we put

[TABLE]

This completes the proof. ∎

Since (4.8) satisfies the conditions of The Iterative Lemma with $n=1$ and $\Lambda_{1}(H_{1},x)\equiv 0$, $\Omega_{1}(H_{1},x)\equiv\Omega_{0}(H_{1})$, we conclude that for every $J\subset I$ and every $N\in\mathbb{N}_{0}$, we can transform (4.8) into “the normal form”:

[TABLE]

on $J\times\mathbb{T}\times(-x_{N+1},x_{N+1})$ for $x_{N+1}>0$ sufficiently small, by a near-identify transformation of $H$ and $\phi$. We drop the subscripts henceforth and write it in the equivalent form

[TABLE]

for $x>0$.

Lemma 4.5**.**

The flow of (4.17), $\underline{H}(\tau,H_{0},\phi_{0},x_{0}),\underline{\phi}(\tau,H_{0},\phi_{0},x_{0})$, $\underline{x}(\tau,H_{0},\phi_{0},x_{0})=x_{0}-\tau$ with

[TABLE]

for $z=H,\phi$, is $C^{\infty}$ on the set

[TABLE]

*for $\xi>0$ sufficiently small. *

Proof.

The right hand side is smooth, even real-analytic, and the existence of a smooth local flow $\underline{H}(\tau,H_{0},\phi_{0},x_{0}),\underline{\phi}(\tau,H_{0},\phi_{0},x_{0}),\underline{x}(\tau,H_{0},\phi_{0},x_{0})=x_{0}-\tau$, with $\tau\in I(H_{0},\phi_{0},x_{0}):=(0,\tau_{\text{max}})$, therefore follows. We then integrate both sides of (4.17):

[TABLE]

with $x=x_{0}-\tau$, and $[\cdots]$, $\{\cdots\}$ being the brackets in (4.17) evaluated at $H=\underline{H}(x_{0}-x,H_{0},\phi_{0},x_{0}),\phi=\underline{\phi}(x_{0}-x,H_{0},\phi_{0},x_{0})$. It is then standard to arrive at the following estimate for $\tau\in I(H_{0},\phi_{0},x_{0})$:

[TABLE]

for some $C>0$ large enough and all $(H_{0},\phi_{0},x_{0})\in J\times\mathbb{T}\times(0,\xi)$, provided that $\xi>0$ is small enough. This shows that $\tau_{\text{max}}=x_{0}$ and completes the proof.

∎

4.3. The existence of smooth invariant manifolds

The inequalities in (4.20) provide $C^{0}$-estimates of $\underline{H}$ and $\underline{\phi}$, respectively.

Lemma 4.6**.**

$\underline{H}(\tau,H_{0},\phi_{0},x_{0})$* extends continuously to the closure $\overline{V(\xi)}$, recall (4.18), with $\underline{H}(0,H_{0},\phi_{0},0)=H_{0}$ for all $(H_{0},\phi_{0})\in J\times\mathbb{T}$.*

Remark 4.7**.**

This is true even though $\underline{\phi}$ itself does not extend to $\tau=x_{0}$.

Proof.

Consider $\underline{\phi}(\tau,H_{0},\phi_{0},x_{0})$ given for $\tau\in[0,x_{0})$. It is continuous on $V(\xi)$. Moreover, $(\tau,x_{0})\mapsto\underline{H}(\tau,H_{0},\phi_{0},x_{0})$ is absolutely continuous, uniformly in $(H_{0},\phi_{0})\in J\times\mathbb{T}$, see (4.19). Consequently, $\underline{H}$ extends continuously and uniquely to the closure $\overline{V(\xi)}$. Moreover, $\underline{H}(0,H_{0},\phi_{0},x_{0})=H_{0}$ by definition for all $x_{0}\in(0,\xi)$ and therefore also $H(0,H_{0},\phi_{0},0)=H_{0}$.∎

Let $h\in J_{0}\subset J$. We then consider the resulting equation

[TABLE]

with $(H,\phi,x)\mapsto\underline{H}(x,H,\phi,x)$ being defined and continuous by Lemma 4.6 on $J\times\mathbb{T}\times[0,\xi]$.

Lemma 4.8**.**

The equation (4.21) defines an invariant set $W^{s}(\Gamma_{1}(h))$ in the $(H,\phi,x)$-space.

Proof.

$(x_{0},H_{0},\phi_{0})\in W^{s}(\Gamma_{1}(h))\Longrightarrow$

[TABLE]

using the group properties of the flow. Consequently, $(\underline{H}(s,H_{0},\phi_{0},x_{0}),\underline{\phi}(s,H_{0},\phi_{0},x_{0}),x_{0}-s)\in W^{s}(\Gamma_{1}(h))$ for all $s\in(0,x_{0})$. ∎

In terms of the time $s$ (say) used in (4.16), $\tau\rightarrow x_{0}^{-}$ corresponds to $s\rightarrow\infty$, and since

[TABLE]

recall (2.8), see also Lemma 4.2 and Lemma 4.4, we conclude that $W^{s}(\Gamma_{1}(h))$ is the stable set of $\Gamma_{1}(h)$, $h\in J_{0}\subset(0,\frac{1}{2})$, as desired.

Proposition 4.9**.**

Consider (4.17) with $N\in\mathbb{N}$ fixed and suppose that there is an $M\in\mathbb{N}$ so that $\underline{H}$ extends as a $C^{M}$-smooth function to the closure $\overline{V(\xi)}$, recall (4.18), satisfying

[TABLE]

Then for any $J_{0}\subset J$ there exists a $\xi>0$ such that the following holds. (4.21) with $h\in J_{0}$ has a unique solution for $(H,\phi,x)\in J\times\mathbb{T}\times[0,\xi]$ of the following graph form

[TABLE]

with $F\in C^{M}$ on $J_{0}\times\mathbb{T}\times[0,\xi]$.

Proof.

Follows directly from the implicit function theorem. Indeed, we have

[TABLE]

for all $\phi\in\mathbb{T}$ by (4.22), and the right hand side of (4.21) is well-defined and $C^{M}$ with respect to $(H,\phi,x)\in J\times\mathbb{T}\times[0,x_{0}]$ by assumption. Consequently, we can solve (4.21) for $H$ as a function of $h,\phi,x$. This gives (4.23) and $F$ is $C^{M}$-smooth since $\underline{H}$ is so. ∎

Under the assumptions of Proposition 4.9, we then have the following by returning to $l=x^{3}$: For any $h\in J_{0}$, (4.23) parametrizes $W^{s}(\Gamma_{1}(h))$ locally in the $(H,\phi,l)$-space as a $C^{M}$-smooth graph $H=F(h,\phi,l^{3})$, $\phi\in\mathbb{T}$, $l\in(0,\xi^{1/3})$. Now, by Lemma 4.2 and Lemma 4.4, it follows that $(r_{1},v,l)\mapsto(H,\phi,l)$, with $l\in(0,\xi^{1/3})$ for $\xi>0$ small enough, is a smooth diffeomorphism (on the relevant set). Therefore we obtain a $C^{M}$ invariant manifold $W^{s}(\Gamma_{1}(h))$ as the stable set of $\Gamma_{1}(h)$ in the original $(r_{1},v,l)$-space, as desired.

Consequently, in order to finish the proof of Theorem 2.7, it suffices to verify the conditions of Proposition 4.9 and to note that $M$ can be taken to be arbitrary (upon increasing $N$). We will show that we can take $M=\lfloor\frac{N}{2}\rfloor$. Here $\lfloor x\rfloor$ for $x\in\mathbb{R}$ is the floor function, i.e. $n=\lfloor x\rfloor$ is the largest integer such that $n\leq x$. $\underline{H}$ clearly satisfies (4.22) once we have shown that it extends $C^{M}$-smoothly to $\overline{V(\xi)}$.

Define

[TABLE]

for $z=H,\phi,x$ and where $\mathbf{\nu}=(\nu_{1},\nu_{2},\nu_{3},\nu_{3})\in\mathbb{N}_{0}^{4}$, $|\nu|=\nu_{1}+\cdots+\nu_{4}\geq 0$.

Lemma 4.10**.**

Fix $N\in\mathbb{N}$ with $N\geq 3$ and let $M=\lfloor\frac{N}{2}\rfloor$. Then there exists a constant $C_{N}>0$ large enough, and a constant $\xi>0$ small enough such that

[TABLE]

for $(\tau,H_{0},\phi_{0},x_{0})\in\overline{V(\xi)}$, and all $|\nu|\leq M$.

Proof.

The proof is delayed to the Appendix A. It rests upon careful estimation of the higher order variational equations of (4.17). The main difficulty lies in estimating $\underline{\phi}_{\mathbf{\nu}}$, due to the singular nature of (4.17) at $x=0$. We find that

[TABLE]

for $(\tau,H_{0},\phi_{0},x_{0})\in V(\xi)$, and all $|\nu|\leq M$. Upon using that the $\phi$-dependent term in the $H$-equation has a $x^{N+1}$-factor, this allow us to control and extend (as in the proof of Lemma 4.6) $\underline{H}_{\mathbf{\nu}}$ to the closure $\overline{V(\xi)}$ provided that $N>\nu_{1}+\nu_{4}$.

∎

5. Blowing up the linearly damped Kepler problem

In Section 2, we characterized constant values of $|\mathcal{E}_{\infty}|\in(0,1)$ as a smooth cylinder in the $(r_{1},v,l)$-space, see Theorem 2.7. $\mathcal{E}_{\infty}=0$, on the other hand, became a one-dimensional manifold, see Theorem 2.6. At the same time, since $|\mathcal{E}_{\infty}|=1\Rightarrow L=0$, see [35], $|\mathcal{E}_{\infty}|=1$ is the plane defined by $(r,v,0)$ in the $(r,v,l)$-space. $|\mathcal{E}_{\infty}|=0$ and $|\mathcal{E}_{\infty}|=1$ are therefore special cases where the associated invariant sets bifurcate. From [35, 36], it is also known that $|\mathcal{E}_{\infty}|\in[0,1]$ and that all values are attained in this set. In our characterization of $|\mathcal{E}_{\infty}|$ through the Hamiltonian function $H$, see (2.9) and (2.10), this means the following:

Lemma 5.1**.**

The set $\Gamma_{1}(h)$ defined by $H(r_{1},v)=h,\rho_{1}=0$ with $h\geq\frac{1}{2}$, is not an $\omega$-limit set.

Proof.

$H(r_{1},v)>\frac{1}{2}$ is obvious since then $|\mathcal{E}_{\infty}|>1$, which would contradict [35]. Moreover, although $H(r_{1},v)=\frac{1}{2}$ corresponds to $|\mathcal{E}_{\infty}|=1$ it cannot be an $\omega$-limit set either, because of $|\mathcal{E}_{\infty}|=1\Rightarrow L=0$, cf. [35], which is an invariant set (that does not contain $\Gamma_{1}(h)$). ∎

We obtained this result as a corollary of [35]. It cannot be understood directly from the perspective in Section 2. However, in this section, we will perform a thorough geometric description of the dynamics of (1.1) within the orbital plane, i.e. $u=re^{i\theta}$, by using the blowup method (as well as desingularization and compactification). In this way, we obtain a system where all singularities have associated eigenvalues with nonzero real part (with the exception of $q$ which only has imaginary eigenvalues), and this allows us to interpret Lemma 5.1 in a separate geometric way, which will also shed light on $|\mathcal{E}_{\infty}|\rightarrow 1$ (corresponding to $h\rightarrow\frac{1}{2}$ cf. (2.10)), see also Section 6.

Our starting point for our blowup approach is to use the coordinates $(r,v,l)$, recall (2.4) and (2.2). This produces the following system

[TABLE]

We now define a reparametrization of time, corresponding to multiplication of the right hand side of (5.1) by $lr^{3}$:

[TABLE]

(5.1) and (5.2) are equivalent for $r>0$, $l>0$, but (5.2) has the advantage of being well-defined on $r=0,l=0$. This enables the use of dynamical systems theory to infer properties from $r=0,l=0$, to $r>0,l>0$ for (5.2) and therefore also (5.1).

The system (5.2) has two invariant planes: $\{r=0\}$ where

[TABLE]

and $\{l=0\}$ where

[TABLE]

see Fig. 3. Their intersection $V=\{r=l=0\}$ and $P=\{v=l=0\}$ consist entirely of completely degenerate equilibria $v^{\prime}=0$, insofar that the linearization about any point in $V$ or $P$ has only zero eigenvalues. We therefore apply the blowup method, see [14, 15]. We proceed as follows: First, we blow up $V$ by application of the following cylindrical blowup transformation

[TABLE]

which fixes $v$. See Fig. 4. Let $X$ denote the vector-field in (5.2). Then the blowup weights are chosen so that $\overline{X}=\Phi^{V}_{*}X$ has $\rho^{4}$ as a common factor. It is therefore the desingularized vector-field $\rho^{-4}\overline{X}$, that we study in the following. To perform calculations, we work in two separate charts, that we will denote by $(\bar{l}=1)_{1}$ and $(\bar{r}=1)_{2}$, respectively, with chart-specific coordinates $(r_{1},v,\rho_{1})$ and $(\rho_{2},v,l_{2})$ defined by

[TABLE]

See also Fig. 4 for an illustration of these coordinates. The following expressions

[TABLE]

define the smooth change of coordinates for $\rho_{1}\geq 0,r_{1}>0$. We will achieve the desingularization through division of the local vector-fields by $\rho_{1}^{4}$ and $\rho_{2}^{4}$, respectively.

Remark 5.2**.**

We used the $r_{1}$-coordinate of (5.4) in Section 2, see (2.4), to prove Theorem 2.6 and Theorem 2.7. The coordinates in the $(\bar{r}=1)_{2}$-chart allow us to follow and analyze the unbounded orbits $\Gamma(h)$ for $h\geq\frac{1}{2}$.

Subsequently, in the $(\bar{r}=1)_{2}$-chart, we will find that the set $P_{2}$ defined by

[TABLE]

and $\rho_{2}\geq 0$, is a set of completely degenerate equilibria; this set clearly corresponds to $P$. We therefore blowup $P_{2}$ by application of the following cylindrical blowup transformation:

[TABLE]

leaving $\rho_{2}$ fixed. See Fig. 5. Let $X_{2}$ denote the local vector-field in the $(\bar{r}=1)_{2}$-chart. Then $\overline{X}_{2}=\Phi^{P_{2}}_{*}X_{2}$ has $\mu$ as a common factor. It is therefore the desingularized vector-field $\mu^{-1}\overline{X}_{2}$, that we study in the following. To perform calculations, we work in three separate charts $(\bar{r}=1,\bar{l}_{2}=1)_{21}$, $(\bar{r}=1,\bar{v}=1)_{22}$ and $(\bar{r}=1,\bar{v}=-1)_{23}$ with chart-specific coordinates $(\rho_{2},v_{1},\mu_{1})$, $(\rho_{2},\mu_{2},l_{22})$ and $(\rho_{2},\mu_{3},l_{23})$ defined by

[TABLE]

See also Fig. 5 for an illustration of these coordinates. The change of coordinates between $(\bar{r}=1,\bar{l}_{2}=1)_{21}$ and $(\bar{r}=1,\bar{v}=1)_{22}$ is given by

[TABLE]

for $\mu_{1}\geq 0$, $v_{1}>0$. Similarly, between $(\bar{r}=1,\bar{l}_{2}=1)_{21}$ and $(\bar{r}=1,\bar{v}=-1)_{22}$ we have

[TABLE]

for $v_{1}<0$. We will achieve the desingularization in each of the charts through division of the local vector-fields by $\mu_{i}$, $i=1,2,3$, respectively.

Remark 5.3**.**

Using (5.5) and (2.4), we have that $v_{1}$ in (5.7) can be written in terms of $r$ and $\dot{r}$ as follows:

[TABLE]

The equations, we obtain in the $(\bar{r}=1,\bar{l}_{2}=1)_{21}$-chart, see Section 5.3 and (5.19), are therefore equivalent to (2.3) on $r=\rho_{2}^{2}>0$. However, the factor of $\rho_{2}$ in (5.12) induces a compactification of $\dot{r}$. See also Remark 5.6.

As we will see, orbits in $(\bar{r}=1,\bar{v}=1)_{22}$ go unbounded (with $\rho_{2}\rightarrow\infty$). We will therefore need to compactify the space $(\rho_{1},\mu_{2},l_{22})$. It turns out that the most convinient way to do this is as follows:

[TABLE]

with $\nu\geq 0,l_{222}\geq 0$, leaving $\mu_{2}$ fixed. In this way, $\nu=0$ corresponds to $\rho_{2}=\infty$ (or $r=\infty$ by (5.4)) and $l_{22}=0$. The latter property may seem unnatural, but upon using (5.10), we may realize that it leads to the following compactification of the $(\rho_{2},v_{1},\mu_{1})$-space associated with the $(\bar{r}=1,\bar{l}_{2}=1)_{21}$-chart:

[TABLE]

where $v_{11}=l_{222}^{-1}$, $\mu_{11}=\mu_{2}l_{222}$. Using (5.11), we then also obtain the following compactification in the $(\bar{r}=1,\bar{v}=-1)_{23}$-chart:

[TABLE]

with $\nu\geq 0,l_{233}\geq 0$, leaving $\mu_{3}$ fixed.

We present the final geometric picture of our compactified phase space in a schematic way in Fig. 6. This figure also illustrates the coordinates used at $r=\infty$.

In the following section, we study the dynamics in each of the charts.

Remark 5.4**.**

In fairness, orbits within $\{r=0\}$ are also unbounded in the $v$-direction for $l>0$, see Fig. 3, and from this perspective, it is also desirable to compactify the $v$-direction. However, for simplicity we have chosen not to include this.

5.1. Chart $(\bar{l}=1)_{1}$

By inserting (5.4) into (5.2), we obtain

[TABLE]

after desingularization, corresponding to division of the right hand side by $\rho_{1}^{4}$. Setting $\rho_{1}=l$ and dividing the right hand side by $r_{1}^{3}$, we obtain (2.6). Consequently, within $\rho_{1}=0$, we have the Hamiltonian system with Hamiltonian function $H(r_{1},v)=\frac{1}{2}v^{2}+\frac{(r_{1}-1)^{2}}{2r_{1}^{2}}$, recall (2.8), having periodic orbits $\Gamma_{1}(h)$, $h\in\left(0,\frac{1}{2}\right)$ within $\rho_{1}=0$, surrounding the center $(r_{1},v)=(1,0)$. Theorem 2.6 and Theorem 2.7 gave the existence of stable manifolds of $(r_{1},v,\rho_{1})=(1,0,0)$ and $\Gamma_{1}(h)$, $h\in\left(0,\frac{1}{2}\right)$. The orbit $\Gamma_{1}\!\left(\frac{1}{2}\right)$, defined by $H(r_{1},v)=\frac{1}{2}$ within $\rho_{1}=0$, is a separatrix, separating bounded (periodic) orbits from the unbounded ones ($H(r_{1},v)\geq\frac{1}{2}$), see Lemma 2.4 and Fig. 2.

5.2. Chart $(\bar{r}=1)_{2}$

By inserting (5.5) into (5.2), we obtain the equations:

[TABLE]

after desingularization (corresponding to division of the right hand side by $\rho_{2}^{4}$). Here we have two invariant planes defined by $\rho_{2}=0$ and $l_{2}=0$. Within the former, we rediscover the center at $(v,l_{2})=(0,1)$ and the periodic orbits $\Gamma_{2}(h)$, $h\in\left(0,\frac{1}{2}\right)$, given by

[TABLE]

surrounding the center. However, $\Gamma_{1}\!\left(\frac{1}{2}\right)$ now becomes a bounded orbit $\Gamma_{2}\!\left(\frac{1}{2}\right)$, which is homoclinic to the degenerate point $\gamma_{2}$ defined by $(\rho_{2},v,l_{2})=(0,0,0)$, see Fig. 7. All other points $(0,v,0)$ on the $v$-axis, are partially hyperbolic, the linearization having eigenvalues $\pm\frac{1}{2}$, [math]. All orbits $H_{2}(l_{2},v)=h$ with $h>\frac{1}{2}$ are heteroclinic connections within $\rho_{2}=0$ between points on the $v$-axis. Next, within $l_{2}=0$ we have

[TABLE]

The $\rho_{2}$-axis is the set of degenerate equilibria $P_{2}$, which is blown up by (5.6). Notice that $H_{2}(0,0)=\frac{1}{2}$, and since $H_{2}$ is independent of $\rho_{2}$, we conclude using (2.10) that

Lemma 5.5**.**

$|\mathcal{E}_{\infty}|=1$* on $P_{2}$.*

5.3. Chart $(\bar{r}=1,\bar{l}=1)_{21}$

By inserting (5.7) into (5.17), we obtain the following equations

[TABLE]

after desingularization (corresponding to division of the right hand side by $\mu_{1}$). $\mu_{1}=0$ corresponds to the blowup of $P_{2}$, and within this invariant subspace, we have that

[TABLE]

We have two hyperbolic equilibria $\gamma_{21}^{\pm}$ along $\rho_{2}=0$ given by $v_{1}=\pm\sqrt{2}$, the former is an unstable node while the latter is a stable node. As equilibria points $(0,\pm\sqrt{2},0)$ of the full system, they are hyperbolic saddles and the separatrix $\Gamma_{1}\!\left(\frac{1}{2}\right)$ from $(\bar{r}=1)_{1}$, now denoted by $\Gamma_{21}\!\left(\frac{1}{2}\right)$ and given by $H_{2}(\mu_{1},\mu_{1}v_{1})=\frac{1}{2}$ within $\rho_{2}=0$, is now a heteroclinic orbit connecting the two hyperbolic saddles. We summarize the findings in Fig. 8.

Remark 5.6**.**

Following Remark 5.3, (5.20) is equivalent to (2.3) with $l=0$:

[TABLE]

In particular, setting $v_{1}=\rho_{2}u$ with $u=\dot{r}$ and $\rho_{2}=\sqrt{r}$ transforms (5.20) into

[TABLE]

*which is studied in [35, Proposition 3.1]. The advantage of working with (5.20) is that the $u$-axis becomes compactified. In particular, the heteroclinic orbits, connecting $\gamma_{21}^{\pm}$ within $\mu_{1}=0$, see Fig. 8, (called ejection-collision orbits in [34]) become unbounded in the $(r,u=\dot{r})$-coordinates. *

5.4. Chart $(\bar{r}=1,\bar{l}=1)_{22}$

By inserting (5.8) into (5.17), we obtain the following equations

[TABLE]

after desingularization (corresponding to division of the right hand side by $\mu_{2}$). All of the three invariant planes defined by $\mu_{2}=0$, $\rho_{2}=0$ and $l_{22}=0$, respectively, are invariant. Within $\mu_{2}=0$, we find

[TABLE]

Here we rediscover $\gamma_{21}^{+}$ from the $(\bar{r}=1,\bar{l}=1)_{21}$-chart as a hyperbolic unstable node within $\mu_{2}=0$ given by

[TABLE]

see also (5.10). At the same time, $(\rho_{2},l_{22})=(0,0)$ is a hyperbolic saddle for (5.22), the linearization having eigenvalues $\pm\frac{1}{2}$. The two axes, $\rho_{2}$ and $l_{22}$, are the associated unstable and stable manifolds, respectively.

Next, within $l_{22}=0$, we have

[TABLE]

and the $\mu_{2}$-axis is therefore a line of saddle points of (5.21), having $l_{22}=0$ ($\rho_{2}=0$) as its unstable manifold (stable manifold, respectively). We summarize the findings in Fig. 9.

5.5. Chart $(\bar{r}=1,\bar{l}=1)_{23}$

By inserting (5.9) into (5.17), we obtain the following equations

[TABLE]

after desingularization (corresponding to division of the right hand side by $\mu_{3}$). The analysis in this chart is almost identical to the analysis in $(\bar{r}=1,\bar{l}=1)_{22}$. In particular, we rediscover the equilibrium point $\gamma^{-}_{21}$ from the $(\bar{r}=1,\bar{l}=1)_{21}$-chart, now given by

[TABLE]

At the same time, we have that the $\mu_{3}$-axis is line of saddle points of (5.23), having $l_{22}=0$ ($\rho_{2}=0$) as its stable manifold (unstable manifold, respectively). An illustration of the dynamics in this chart can be obtained by taking the time-reversal of the diagram Fig. 9 (replacing (a) $l_{22}$, $\mu_{2}$ by $l_{23}$ and $\mu_{3}$, respectively, (b) $\gamma_{22}^{+}$ by $\gamma_{23}^{-}$ and finally (c) replacing (in line with Fig. 6) the color green of the invariant plane $l_{22}=0$ to red).

5.6. Compactification in the $(\bar{r}=1,\bar{l}=1)_{21}$-chart

Upon applying the transformation of $(\rho_{2},v_{1})$ defined by (5.14) to the system (5.19), we obtain the following equations

[TABLE]

after desingularization (corresponding to multiplication of the right hand side by $\nu^{3}$). Here $\nu=0$, which corresponds to $r=\infty$, is an invariant set, upon which we find the following:

[TABLE]

We find two equilibria: $p_{21}^{+}:\,(\nu,v_{11},\mu_{12})=(0,0,0)$ and $p_{21}^{-}:\,(\nu,v_{11},\mu_{12})=(0,-\delta,0)$ of (5.24), both of which are hyperbolic for the reduced system (5.25) within $\nu=0$. Indeed, the linearization of (5.25) around $p_{21}^{+}$ produces the eigenvalues $-\delta,-\delta$ (semi-simple), whereas the linearization of (5.25) around $p_{21}^{-}$ produces $\delta,-2\delta$ as eigenvalues. Consequently, $p_{21}^{+}$ is a stable node for (5.25), whereas $p_{21}^{-}$ is a saddle.

$\mu_{11}=0$ is also an invariant set for (5.24), upon which we find the following:

[TABLE]

While $p_{21}^{-}$ is clearly an unstable node for these equations, $p_{21}^{+}$ is only semi-hyperbolic, the linearization having eigenvalues [math] and $-\delta$. It is a simple calculation, to show that the associated center manifold $W^{c}(p_{21}^{+})$ takes the following graph form:

[TABLE]

The $\mathcal{O}(\nu^{6})$-term is a smooth function of $\nu^{6}$ (and $\delta$). This follows from the fact that the first equation of (5.26) can be written as $(\nu^{6})^{\prime}=-3\nu^{6}v_{11}$. Inserting (5.27) into (5.26) gives

[TABLE]

and $\nu>0$ is therefore increasing on $W^{c}(p_{21}^{+})$. The center manifold is therefore unique as the unstable set of $p_{21}^{+}$ and $p_{21}^{+}$ is a nonhyperbolic saddle for (5.26). We illustrate our findings in Fig. 10.

Remark 5.7**.**

It is a simple calculation to show that the time used in (5.24) coincides with the original time in (5.1), which is why we use $\dot{()}$ instead of $()^{\prime}$, recall Remark 2.1.

Using (5.5), (5.12), and (5.14), we can write (5.27) in the $(r,\dot{r})$-plane as a graph

[TABLE]

over $r\gg 1$. We see that $\dot{r}\rightarrow 0$ on $W^{c}$ as $t\rightarrow-\infty$, in line with the results of [34, Proposition 3.1].

5.7. Compactification in the $(\bar{r}=1,\bar{l}=1)_{22}$-chart

Upon applying the transformation of $(\rho_{2},l_{22})$ defined by (5.13) to the system (5.21), we obtain the following equations

[TABLE]

Here $l_{222}=0$ defines an invariant set upon which we have $\nu^{\prime}=-\frac{1}{2}\nu$ and $\mu_{2}^{\prime}=0$. Since $l_{222}=0$ corresponds to $l_{22}=0$, see (5.13), these findings are obviously in agreement with the results in the $(\bar{r}=1,\bar{l}=1)_{22}$-chart whenever $\nu>0$, see Section 5.4 and Fig. 9 (green plane). On the other hand, $\nu=0$, corresponding to $r=\infty$, is now an invariant set of (5.28). In fact, $\nu=l_{222}=0$ is a line of saddle points; the linearization of (5.28) about any point in this set having eigenvalues $-\frac{1}{2},0,1$. The stable manifold of $\nu=l_{222}=0$ is the $(\nu,\mu_{2})$-plane, whereas the associated unstable set is the $(\mu_{2},l_{222})$-plane. Setting $\nu=0$ in (5.28) gives

[TABLE]

which has no equilibria within the first quadrant. In particular, $l_{222}>0$ ($\mu_{2}$) is monotonically increasing (decreasing, respectively). Finally, within $\mu_{2}=0$ we have

[TABLE]

having a hyperbolic saddle at the origin. We summarize the local findings in Fig. 11.

5.8. Compactification in the $(\bar{r}=1,\bar{l}=1)_{23}$-chart

Upon applying the transformation of $(\rho_{2},l_{22})$ defined by (5.13) to the system (5.21), we obtain the following equations

[TABLE]

The analysis of these equations is almost identical to the analysis performed in Section 5.6 and Section 5.7. We therefore only present a diagram, see Fig. 12.

Upon combining all of our findings in the local charts, we obtain the global perspective in Fig. 1 (using the geometric viewpoint in Fig. 6). Notice that from this perspective it follows directly that any $\Gamma(h)$ with $h>\frac{1}{2}$ cannot be the $\omega$-limit set. We elaborate further upon this in the following section.

6. Discussion

In this paper, we have revisited the linearly damped Kepler problem (1.1) with the main purpose of describing the smoothness of the invariant manifolds obtained in [35], see Theorem 2.7 and Theorem 2.8. In the process, we identified a separate invariant manifold $W^{s}(q)$, see Theorem 2.6. This one-dimensional manifold of (2.3), corresponding to orbits becoming more circular as $t\rightarrow\infty$, acts as the center of oscillations along the invariant manifolds of Theorem 2.7 and Theorem 2.8, see Fig. 1.

Finally, in Section 5 we performed a blowup analysis that led to the geometric description of the dynamics illustrated in Fig. 1. In future work, we aim to use a similar approach to study (1.1) with the general $(u,\dot{u})$-dependent nonlinear damping:

[TABLE]

for $\beta,\alpha\in\mathbb{N}_{0}$ and $k>0$. This family was considered by Poincaré, and he argued formally, see e.g. [32], that for $\alpha$ and $\beta$ sufficiently large, orbits tend to “circularize”, in the sense that $\mathcal{E}\rightarrow 0$ for $t\rightarrow\infty$ on an open set. To the best of our knowledge, this remains a conjecture to this today. We hope to solve this conjecture and determine all $\alpha$ and $\beta$ for which “circularization” occurs using a similar blowup approach. We emphasize that it is known that circularization does not occur for linear damping (indeed, it is exceptional occording to Theorem 2.6 and Theorem 2.7) nor does it occur for the Poynting-Plummer-Damby damping ($\alpha=0$, $\beta=2$). See [9] and [32] for a general analysis of the case $\alpha=0$, $\beta>0$.

In the following, we will use the geometric framework of the blowup approach to shed light on $|\mathcal{E}_{\infty}|\rightarrow 1$, recall Lemma 5.1.

Firstly, in Fig. 13 we have used Matlab’s ODE45 to simulate the system for initial conditions near the set $\Gamma_{2}(0.8)$ defined by $\rho_{2}=0$, $H_{2}(l_{2},v)=0.8>\frac{1}{2}$, recall (5.18). More specifically, in (a) we use the coordinates $(\rho_{2},v,l_{2})$ of the $(\bar{r}=1)_{2}$-chart, see (5.5), and select $50$ initial conditions (indicated by the cyan cylinder near $v=0$) with $(l_{2},v,\rho_{2})=(l_{20},0,\rho_{20})$, $H_{2}(l_{20},0)=0.8>\frac{1}{2}$, and $\rho_{20}$ in the interval $(0.15,0.25)$ (equispaced). The curve $\Gamma_{2}(0.8)$ is shown in red and the forward orbits (also in cyan) follow this red curve, but only up until a vicinity of the $v$-axis (due to the saddle-structure along this set, see Fig. 7 and the discussion of Fig. 14 below). From here $\rho_{2}$ increases and the orbits eventually contract towards the degenerate set $P_{2}$. In (b), we use the coordinates $(\rho_{2},v_{1},\mu_{1})$ of the $(\bar{r}=1,\bar{l}_{2}=1)_{21}$-chart, where $P_{2}$ has been blown up, to illustrate the same (computed) orbits. Here we see that the motion along $P_{2}$ in (a) is due to the contraction towards the center manifold $W^{c}(p_{21}^{+})$; this is also shown in (c) using a projection onto the $(\nu,v_{11})$-plane, recall (5.14). Beyond the motion along $W^{c}(p^{+})$ in Fig. 13 (b), the cyan curves come close together near $\gamma_{21}^{-}$ and follow the heteroclinic orbit $\Gamma_{21}\!\left(\frac{1}{2}\right)$, connecting $\gamma_{21}^{-}$ and $\gamma_{21}^{+}$. Once passing close to $\gamma_{21}^{+}$ the cyan curves move close to $\mu_{1}=0$ plane again before returning to $\gamma_{21}^{-}$ and $\Gamma_{21}\!\left(\frac{1}{2}\right)$. This process repeats itself. The computations were performed with low tolerances ($10^{-12}$) and ended when the value of $l$ reached $10^{-4}$.

The results in Fig. 13 are in agreement with the findings in Section 5 (see Fig. 1). We emphasize this point further in Fig. 14, where we illustrate a “singular” orbit of a point on the orange cylinder. Here singular refers to the fact that it lies on the blowup space (it is therefore not a true orbit of the system (1.1) since it occurs on $r=0$) and at the same time we use a concatenation of unstable and stable manifolds across saddle points. The forward orbit of initial conditions starting near the cyan curve on the orange cylinder (with $0<l\ll 1$) (as in Fig. 13 in the local charts $(\bar{r}=1)_{2}$ and $(\bar{r}=1,\bar{l}_{2}=1)_{21}$) will track this curve (meaning that the orbit remains close to local copies within compact subsets of the local charts) up until $\gamma^{+}$. This follows from Section 5, but we will leave out further details.

At $\gamma^{+}$ there is no unique forward (singular) orbit, since $\gamma^{+}$ is an unstable node on $\mu=0$ (the $l=0$ cylinder corresponding to the blowup of $P_{2}$, see (5.6)), see also Fig. 8. Nevertheless, since (a) $\gamma^{+}$ connects to $\gamma_{-}$, (b) $\gamma_{-}$ is contracting on $\rho=0$, and (c) $\Gamma\!\left(\frac{1}{2}\right)$ connects back to $\gamma^{+}$ this process repeats itself for an actual orbit, starting near the cyan disc with $0<l\ll 1$. However, since $l$ is monotonically decreasing, the excursions from $\gamma^{+}$ to $\gamma^{-}$ will get closer and closer to the purple connection (along the half-circle $(\bar{v},\bar{l}_{2})\in S^{1}$ $\bar{l}_{2}\geq 0$ within $\rho_{2}=0$). We observe this numerically (also indicated in Fig. 13 (b) by the blue segments close to $\mu_{1}=0$).

However, and this is the main point, since $\Gamma_{1}\!\left(\frac{1}{2}\right)$ cannot be the $\omega$-limit set, see Lemma 5.1, and since $l=0$ is invariant, the $\omega$-limit set of the initial conditions in Fig. 13 (or near the cyan disc in Fig. 14) has to be $\Gamma_{1}(h)$ with $h<\frac{1}{2}$ but $h\sim\frac{1}{2}$. In other words, these initial conditions belong to a $W^{s}(\Gamma(h))$ for some $h\in(0,\frac{1}{2})$, where $h\rightarrow\frac{1}{2}$ as $l_{0}\rightarrow 0$. The implication of this is that the cylinders $W^{s}(\Gamma(h))$ have a dramatic limit $h\rightarrow\frac{1}{2}$ (corresponding to $|\mathcal{E}_{\infty}|\rightarrow 1^{-}$ see (2.10)). Indeed, the two-dimensional cylindrical manifolds $W^{s}(\Gamma(h))$ with $h<\frac{1}{2}$ but $h\sim\frac{1}{2}$ open up like a flower, stretching (and collapsing for $h\rightarrow\frac{1}{2}$) both along the orange and the purple cylinder in Fig. 14 (corresponding to $r$ very small and very large, respectively).

At the same time, $l=0$ defines an invariant set (of collinear orbits) and the $\omega$-limit set of any point within this set is $\gamma^{-}$. Orbits having $\gamma^{+}$ as the $\alpha$-limit set (within $l=0$) are called ejection-collision orbits, whereas orbits having either $p^{-}$ or $p^{+}$ as the $\alpha$-limit set (going along $W^{s}(p_{21}^{+})$ in the latter case) are referred to as capture-collision orbits, see e.g. [9, 34]. Interestingly, the orbits in Fig. 13 with $0<l\ll 1$ exhibit both types of behaviour, first capture-collision by following close to $W^{c}(p_{21}^{+})$ and then subsequently (and repeatedly) ejection-collision type behaviour by passing close to $\gamma^{+}$ and $\gamma^{-}$. We illustrate this in our final figure Fig. 15 plotting $r$ as a function of $t$ along one of the cyan curves in Fig. 13(a).

Appendix A Proof of Lemma 4.10

We consider (4.17), repeated here for convinience

[TABLE]

for $x>0$. We then apply $\frac{\partial^{|\mathbb{\nu}|}}{\partial\tau^{\nu_{1}}\partial H_{0}^{\mathbb{\nu}_{2}}\partial\phi_{0}^{\mathbb{\nu}_{3}}\partial x_{0}^{\partial\nu_{4}}}$ on both sides of these equations; recall the notation defined in (4.24). We then obtain $\frac{d}{d\tau}\underline{H}_{\mathbb{\nu}}$, $\frac{d}{d\tau}\underline{\phi}_{\mathbb{\nu}}$ and $\frac{d}{d\tau}\underline{x}_{\mathbf{\nu}}$ on the left hand sides, respectively. Obviously,

[TABLE]

where $x_{(1,0,0,0)}=-1,\,x_{(0,0,0,1)}=1$.

To handle the associated right hand sides, we use the Faa di Bruno formula [7]: For $\mathbf{\nu}=(\nu_{1},\cdots,\nu_{d})\in\mathbb{N}_{0}$, $z=(z_{1},\cdots,z_{d})$, $d\in\mathbb{N}$, define

[TABLE]

Moreover, if we also consider $\mathbf{\mu}=(\mu_{1},\ldots,\nu_{d})$, then we write $\mathbf{\mu}\prec\mathbf{\nu}$ provided at least one of the conditions hold:

1. (1)

$|\mu|<|\nu|$; 2. (2)

$|\mu|=|\nu|$ and $\mu_{1}<\nu_{1}$; or 3. (3)

$|\mu|=|\nu|$, $\mu_{1}=\nu_{1},\ldots,\mu_{k}=\nu_{k}$ and $\mu_{k+1}<\nu_{k+1}$ for some $1\leq k<d$.

Finally, we write $\mathbf{0}=(0,\ldots,0)$ in $\mathbb{N}_{0}^{d}$.

Lemma A.1**.**

[7, Theorem 2.1]** Consider smooth functions $F:\mathbb{R}^{m}\rightarrow\mathbb{R}$ and $G=(G_{1},\ldots,G_{m}):\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}$, and define

[TABLE]

Then

[TABLE]

where $n=|\nu|$ and

[TABLE]

We now consider the following functions:

[TABLE]

of $(\tau,H_{0},\phi_{0},x_{0})\in V(\xi)$, which appear on the right hand side of (A.1). Here $V(\xi)$ is defined in (4.18), repeated here for convinience:

[TABLE]

Lemma A.2**.**

Let $M\in\mathbb{N}$ be so that $N-2M\geq 0$. Suppose that

[TABLE]

for all $(\tau,H_{0},\phi_{0},x_{0})\in V(\xi)$ and all $1\leq|\mathbf{\nu}|\leq M$. Then there exists a constant $K=K(M,|R|_{C^{M}},|P|_{C^{M}},C_{M})>0$, depending only on: (a) $M$, (b) uniform $C^{M}$ bounds $|R|_{C^{M}}$ and $|P|_{C^{M}}$ of $R$ and $P$, respectively, and (c) $C_{M}>0$, such that

[TABLE]

for all $(\tau,H_{0},\phi_{0},x_{0})\in V(\xi)$ and all $1\leq|\mathbf{\nu}|\leq M$.

Proof.

We focus on $W_{R}$. The estimate of $W_{P}$ can be obtained in a completely analogous way. First, we write

[TABLE]

using (A.2). Here we use $D^{\nu}\underline{R}$ to denote the partial differentiation of the composition function $R(\underline{H}(),\underline{\phi}(),\underline{x}())$. Then by the product rule, we have

[TABLE]

using (A.2) again. We then use Lemma A.1, in particular (A.3) with $\nu=(q_{1},\nu_{2},\nu_{3},q_{4})$, $d=4$, $m=3$, $D^{l_{j}}G_{1}=\underline{H}_{l_{j}}$, $D^{l_{j}}G_{2}=\underline{\phi}_{l_{j}}$, $D^{l_{j}}G_{3}=\underline{x}_{l_{j}}$, to estimate $D^{(q_{1},\nu_{2},\nu_{3},q_{4})}\underline{R}$:

[TABLE]

writing $k_{j}=(k_{j,1},k_{j,2},k_{j,3})$, $l_{j}=(l_{j,1},\ldots,l_{j,4})$, for some $K_{0}=K_{0}(M,|R|_{C^{M}})>0$ depending only on: (a) $M$ and (b) the $C^{M}$ bound $|R|_{C^{M}}$ on the smooth function $R$. In the final inequality, we have used the definition of $p_{s}$ (A.4). Specifically, from

[TABLE]

we have concluded that $\sum_{j=1}k_{j,2}(1+l_{j,1}+l_{j,2})\leq q_{1}+q_{2}$; notice that the first inequality in (A.6) is strict. Now using (A.5), we have

[TABLE]

for $K=K(M,|R|_{C^{M}},C_{M})>0$ large enough. ∎

Finally, we describe the the functions

[TABLE]

that also appear on the right hand side (A.1), but are independent of $\underline{\phi}$.

Lemma A.3**.**

Let $M\in\mathbb{N}$ and suppose that

[TABLE]

for all $(\tau,H_{0},\phi_{0},x_{0})\in V(\xi)$ and all $1\leq|\mathbf{\nu}|\leq M$. Then there exists a constant $K=K(M,|D|_{C^{M}},|\Omega|_{C^{M}},C_{M})>0$, depending only on: (a) $M$, (b) uniform $C^{M}$ bounds $|D|_{C^{M}}$ and $|\Omega|_{C^{M}}$ of $D$ and $\Omega$, respectively, and (c) on a constant $C_{M}>0$, such that

[TABLE]

for all $(\tau,H_{0},\phi_{0},x_{0})\in V(\xi)$ and all $1\leq|\mathbf{\nu}|\leq M$.

Proof.

Follows from a direct calculation. ∎

Lemma A.4**.**

Let $M=\lfloor\frac{N}{2}\rfloor$. Then for any $\mathbf{\nu}\in\mathbb{N}_{0}$ with $|\nu|\leq M$ there exists a constant $C_{\mathbf{\nu}}$ and a $\xi_{\mathbb{\nu}}>0$ such that

[TABLE]

for all $(\tau,H_{0},\phi_{0},x_{0})\in V(\xi_{\mathbf{\nu}})$.

Upon proceeding as in the proof of Lemma 4.6, (A.7) implies that $\underline{H}_{\mathbf{\nu}}$ extends continuously to $\overline{V(\xi)}$. Consequently, we complete the proof of Lemma 4.10 by proving Lemma A.4.

To prove Lemma A.4, we proceed by induction. It is true for $\mathbf{\nu}=\textbf{0}$, see (4.20) and the proof of Lemma 4.5. Next, suppose that it is true for $|\mathbf{\nu}|=n$. We then consider $\mathbf{\nu}^{\prime}:=\mathbf{\nu}+\mathbf{\delta}_{i}$, with $\mathbf{\delta}_{i}=(\delta_{i,1},\ldots,\delta_{i,4})$ for

[TABLE]

Lemma A.5**.**

For $\widetilde{C}_{\mathbf{\nu}^{\prime}}>0$ large enough and $\xi>0$ small small enough, we have that

[TABLE]

for all $(0,H_{0},\phi_{0},x_{0})\in V(\xi)$

Proof.

By definition $\underline{H}(0,H_{0},\phi_{0},x_{0})=H_{0}$, $\underline{\phi}(0,H_{0},\phi_{0},x_{0})=\phi_{0}$, $\underline{x}(0,H_{0},\phi_{0},x_{0})=x_{0}$. Consequently, if $\nu=(0,\nu_{2},\nu_{3},\nu_{4})$ then $\underline{z}_{\mathbf{\nu}}(0,\ldots)=0$ for $z=H,\phi,x$ (unless $\underline{H}_{(0,1,0,0)}=1$, $\underline{\phi}_{(0,0,1,0)}=1$, or $\underline{x}_{(0,0,0,1)}=1$). On the other hand, if $\nu_{1}\geq 1$ then by (A.1)

[TABLE]

The result then follows upon using the induction hypothesis, Lemma A.2 and Lemma A.3.

∎

We now consider

[TABLE]

Let $C_{H,\mathbf{\nu}^{\prime}}:=2\widetilde{C}_{\mathbf{\nu}^{\prime}}$ and $C_{\phi,\mathbf{\nu}^{\prime}}>\widetilde{C}_{\mathbf{\nu}^{\prime}}$. We may take $\widetilde{C}_{\mathbf{\nu}^{\prime}}$ large enough such that $\widetilde{C}_{\mathbf{\nu}^{\prime}}\geq C_{\nu}$ for all $|\nu|\leq n$. Then by Lemma A.5, we have that

[TABLE]

for all $\tau\in(0,\tau_{0})$, with $\tau_{0}(H_{0},\phi_{0},x_{0})>0$ small enough. In fact, due to Lemma A.2, and the induction hypothesis, we obtain that

[TABLE]

as $\tau,x_{0}\rightarrow 0$ for any $C_{H,\mathbf{\nu}^{\prime}},C_{\phi,\mathbf{\nu}^{\prime}}>0$. Consequently, by Lemma A.3, we find for $\xi>0$ small enough, that there is a constant $K=K(M,|\Lambda|_{C^{M}},|\Omega|_{C^{M}},\widetilde{C}_{\mathbf{\nu}^{\prime}})$ depending only on: (a) $M$, (b) uniform $C^{M}$ bounds $|\Lambda|_{C^{M}}$ and $|\Omega|_{C^{M}}$ of $\Lambda$ and $\Omega$, respectively, and on (c) $\widetilde{C}_{\mathbf{\nu}^{\prime}}$, such that

[TABLE]

for all $(\tau,H_{0},\phi_{0},x_{0})\in V(\xi)$ and $\tau\in(0,\tau_{0})$ . The main observation here is that by taking $\xi>0$ small enough, we ensure that $K$ is independent of $C_{H,\nu^{\prime}}$ and $C_{\phi,\mathbf{\nu}^{\prime}}$. We then integrate and use

[TABLE]

for any $q>1$. This produces the following

[TABLE]

where we have also used Lemma A.5. We now take

[TABLE]

and upon decreasing $\xi>0$ further, we ensure that $|x_{0}|\leq\frac{\widetilde{C}_{\mathbf{\nu}^{\prime}}}{K}$. Then it follows that (A.8) holds true for all $\tau\in(0,x_{0})$ as desired. This completes the proof of Lemma A.4.