Torus knot choreographies in the $n$-body problem

TL;DR

This paper introduces a systematic, computer-assisted method to prove the existence of spatial choreographies in the gravitational n-body problem, focusing on torus knot solutions with specific symmetries.

Contribution

It develops a constructive approach using delay differential equations and Fourier analysis to establish choreographies with prescribed topological and symmetry properties.

Findings

Proved existence of spatial choreographies for n=4,5,7,9 bodies.

Utilized computer-assisted proofs with analytic estimates.

Reduced the problem to solving a delay differential equation in a Fourier coefficient space.

Abstract

We develop a systematic approach for proving the existence of choreographic solutions in the gravitational $n$ body problem. Our main focus is on spatial torus knots: that is, periodic motions where the positions of all $n$ bodies follow a single closed which winds around a $2$-torus in $\mathbb{R}^3$. After changing to rotating coordinates and exploiting symmetries, the equation of a choreographic configuration is reduced to a delay differential equation (DDE) describing the position and velocity of a single body. We study periodic solutions of this DDE in a Banach space of rapidly decaying Fourier coefficients. Imposing appropriate constraint equations lets us isolate choreographies having prescribed symmetries and topological properties. Our argument is constructive and makes extensive use of the digital computer. We provide all the necessary analytic estimates as well as a working implementation for any number of bodies. We illustrate the utility of the approach by proving the existence of some spatial choreographies for $n=4,5,7$, and $9$ bodies.

Figures

Equations

$$s_1=\frac{1}{4}\sum _{j=1}^{n-1}\frac{1}{\sin (j\zeta /2)},\zeta =\frac{2\pi }{n}.$$

$$w_j(t)=e^{ij\zeta }{w}_n(t+jk\zeta ),$$

$$z_j(t)=z_n(t+jk\zeta ).$$

$$\omega =\sqrt{s_1}\frac{p}{q},$$

$$\begin{array}[c]{rl}\ddot{w}(t)+2\sqrt{s_{1}}i\dot{w}(t)&=\displaystyle s_{1}w(t)-\sum_{j=1}^{n-1}\frac{w(t)-e^{ij\zeta}w(t+jk\zeta)}{\left(\left|w(t)-e^{ij\zeta}w(t+jk\zeta)\right|^{2}+\left|z(t)-z(t+jk\zeta)\right|^{2}\right)^{3/2}}\\ \ddot{z}(t)&=\displaystyle-\sum_{j=1}^{n-1}\frac{z(t)-z(t+jk\zeta)}{\left(\left|w(t)-e^{ij\zeta}w(t+jk\zeta)\right|^{2}+\left|z(t)-z(t+jk\zeta)\right|^{2}\right)^{3/2}}.\end{array}$$

$$\bar{I}=diag(1,1,0)\text{ and }\bar{J}=diag(J,0)\text{,}$$

$$q_j(t)=e^{\sqrt{s_1}t\bar{J}}{u}_j(\omega t),$$

$${\omega }^2{\ddot{u}}_j+2\omega \sqrt{s_1}\bar{J}{\dot{u}}_j-s_1\bar{I}{u}_j=-\sum _{i=1(i\ne j)}^n\frac{u_j-u_i}{{\Vert u_j-u_i\Vert }^3},$$

$$u_j(t)=e^{j\bar{J}\zeta }{u}_n(t+jk\zeta )\text{.}$$

$$\mathcal{G}(u_n,\omega )\stackrel{\text{mbox}[}{=}{\omega }^2{\ddot{u}}_n+2\omega \sqrt{s_1}\bar{J}{\dot{u}}_n-s_1\bar{I}{u}_n+\sum _{j=1}^{n-1}\frac{u_n-e^{j\bar{J}\zeta }{u}_n(t+jk\zeta )}{{\Vert u_n-e^{j\bar{J}\zeta }{u}_n(t+jk\zeta )\Vert }^3}:X\times \mathbb{R}\to Y$$

$$Q_j(t)\stackrel{\text{mbox}[}{=}{q}_j(t/\omega )=e^{t\bar{J}\sqrt{s_1}/\omega }{u}_j(t)$$

$$Q_j(t)=Q_n(t+j\tilde{k}\zeta )\text{,}$$

$$Q_n(t-2\pi )=e_n^{-\bar{J}2\pi q/p}{Q}_n(t),$$

$$r=(kq-p)/n\in \mathbb{Z},$$

$$Q_j(t)=e^{-\bar{J}2\pi j(r/p)}{Q}_n(t+jk\zeta )\text{.}$$

$$Q_n(t-2\pi jr\tilde{q})=e^{-\bar{J}2\pi j(r/p)}{Q}_n(t).$$

$$Q_j(t)=e^{-\bar{J}2\pi j(r/p)}{Q}_n(t+jk\zeta )=Q_n(t+j(k-rn\tilde{q})\zeta ).\text{qed}$$

$$u=(u_1,u_2,u_3).$$

$$(\theta ,\phi ,\tau )u(t)=e^{\bar{J}\theta }u(t+\phi )+(0,0,\tau )\text{,}$$

$$\mathcal{G}((\theta ,\phi ,\tau )u_0)=(\theta ,\phi ,\tau )\mathcal{G}(u_0)=0.$$

$$A_1(u)$$

$$A_2(u)$$

$$A_3(u)$$

$$0={\partial }_j\mathcal{A}(u)={\partial }_j\mathcal{A}((\theta ,\phi ,\tau )u)=⟨\nabla \mathcal{A}(u),A_j(u)⟩=⟨\mathcal{G}(u),A_j(u)⟩\text{,}$$

$$\mathcal{G}(u,\omega )+\sum _{j=1}^3{\lambda }_j{A}_j(u)=0.$$

$$\sum _{j=1}^3{\lambda }_j⟨A_j(u),A_i(u)⟩=0.$$

$$I_j(u)=⟨u-\tilde{u},A_j(\tilde{u})⟩:X\to {\mathbb{R}}^3\text{,}$$

$$I_3(u)={\int }_0^{2\pi }u(t)\cdot (0,0,1)dt={\int }_0^{2\pi }{u}_3(t)dt.$$

$$I_1(u)={\int }_0^{2\pi }\left(u-\tilde{u}\right)\cdot \bar{J}\tilde{u}dt={\int }_0^{2\pi }u\cdot \bar{J}\tilde{u}dt$$

$$I_2(u)={\int }_0^{2\pi }\left(u-\tilde{u}\right)\cdot {\tilde{u}}^{\prime }(t)={\int }_0^{2\pi }u(t)\cdot {\tilde{u}}^{\prime }(t)dt\text{.}$$

$${\omega }^2\dot{v}+2\omega \sqrt{s_1}\bar{J}v-s_1\bar{I}u+\sum _{j=1}^{n-1}\frac{u-e^{j\bar{J}\zeta }u(t+jk\zeta )}{{\Vert u-e^{j\bar{J}\zeta }u(t+jk\zeta )\Vert }^3}=0.$$

Full text

1 Introduction

A choreography is a periodic solution of the gravitational $n$-body problem, where $n$ equal masses follow the same path. Circular choreographies with masses located at the vertices of a regular $n$-gon were already studied by Lagrange in the Eighteenth Century. The first choreography differing from a polygon was discovered by Moore in [1] and has three bodies moving around the now famous figure-eight. Chenciner and Montgomery in [2] gave a rigorous mathematical proof of the existence of this figure eight orbit by minimizing the action for Newton’s equation. The name choreographies was adopted after the work of Simó [3] on numerical computation of choreographic solutions.

The variational approach to the existence of choreographies consists of finding critical points of the classical Newtonian action subject to appropriate symmetry constraints. The main obstacle to this approach is the existence of paths with collisions. Terracini and Ferrario in [4] gave conditions on the symmetries which imply that a minimizer is free of collisions (this is called the rotating circle property). Although a lot of simple choreographies have been found numerically since Simó [3], rigorous proofs using only analytical methods are difficult. Notable exceptions include works on: the figure-eight of three bodies [2], the rotating $n$-gon [5], the figure-eight type for odd bodies [4] and the super-eight of four bodies [6]. Other variational approaches related to existence of planar choreographies can be found in [7, 8, 9, 10, 11, 12] and the references therein.

The difficulties just mentioned have led some authors to develop mathematically rigorous computer assisted proofs (CAPs) for choreographies. This is a natural alternative to pen-and-paper analysis since both the discovery and many subsequent studies of choreographies employ numerical methods. The interested reader will want to consult for example the works of Kapela, Simó, and Zgliczyński [13, 14, 15] for both CAPs of existence for planar choreographies and mathematically rigorous stability analysis. See also Remark 2 below.

Recall now that a $(p,q)$-torus knot is an embedding of $\mathbb{S}^{1}$ into a two torus $\mathbb{T}^{2}\subset\mathbb{R}^{3}$, winding $p$ times around one generating circle of the torus and $q$ times around the other, with $p$ and $q$ coprime and neither equal to zero. The embedding of the two torus is required to be unknotted in $\mathbb{R}^{3}$. A torus knot may or may not be a trivial when viewed as a knot in $\mathbb{R}^{3}$. Indeed, it is trivial if and only if either $p$ or $q$ is equal to $\pm 1$. The idea is illustrated in Figure 1.

A difficult problem in this area is to prove the existence of spatial torus knot choreographies. Indeed when both topological and symmetric constraints are involved, it is difficult to prove the coercitivity of the action. For this reason few results with topological constraints are available. A notable exception is a torus knot choreography for $3$-bodies obtained by Arioli, Barutello, and Terracini in [16], where the authors localize a mountain pass solution of the Newtonian action in a rotating frame. Again the result is obtained by means of CAP, not variational methods. In general it is hard to determine whether a critical point of the action is a spatial torus-knot choreography. We provide a systematic procedure to obtain countable families of torus knots for any number of bodies.

Contribution: The main result of the present work is to give mathematically rigorous existence proofs for $(p,q)$-torus knot choreographies in the $n$-body problem for several different values of $n$.

Our approach is functional analytic (a choreography is a zero of a nonlinear operator posed on a Banach space) and computer-assisted. When it succeeds it produces countably many verified results. For example we establish the existence of the $5$-body trefoil knot choreography illustrated in Figure 2, and the existence of countable many choreographies close to it. We describe the pen and paper estimates for any number of bodies and, while we illustrate the method for only few explicit examples, our setup and resulting implementation apply (in principle) to any spatial choreography.

Before describing our approach in detail we recall several related developments. In [17] it is observed that choreographies appear in dense sets along the vertical Lyapunov families attached to the relative equilibrium solutions given by the planar $n$-gon . Existence of vertical Lyapunov families follows from the Weinstein-Moser theory and, when the frequency varies continuously, the authors obtain the existence of an infinite number of choreographies along these vertical families. This hypothesis however has been verified only for some families with $n=3,4,5,6$ and even though similar computations can be carried out for other values of $n$, it is an open problem to establish the hypothesis for all $n$.

The existence of global Lyapunov families arising from the polygonal relative equilibrium of the rotating problem was established in [18, 19] for all $n$. By saying that these families are global what we mean that, in the space of normalized $2\pi$ periodic solutions, the families form a continuum set with at least one of the following properties: either the Sobolev norm of the orbits in the family goes to infinity, the period of the orbits goes to infinity, the family ends in an orbit with collision, or the family returns to another equilibrium solution. This fact is proved using $G$-equivariant degree theory [20] where $G=\mathbb{Z}_{n}\times\mathbb{Z}_{2}\times SO(2)\times S^{1}$ acts as permutations, $z$-reflection and $(x,y)$-rotations of bodies, and time shift respectively. In addition the analysis of [18, 19] concludes that the Lyapunov families have the symmetries of a twisted subgroup of $G$.

Specifically, let $(w_{j},z_{j})\in\mathbb{C}\times\mathbb{R}$ represents the planar and spatial coordinates of the $j$-th body in a rotating coordinate frame with frequency $\sqrt{s_{1}}$, where

[TABLE]

The $n$-polygon consisting of $n$ bodies on the unit circle $w_{j}=e^{ij\zeta}$ is an equilibrium solution of Newton’s equations in a co-rotating frame. After normalizing the period to $2\pi$, the planar Lyapunov families arising from this equilibrium polygon have the planar symmetries,

[TABLE]

and the spatial symmetries

[TABLE]

For $1\leq j\leq n-1$ the $j$-th body follows an identical path as the $n$-th body, after a rotation in space and a shift in time. It is proved in [19] that taking $k=2,..,n-2$ in the planar case gives the $n-3$ planar Lyapunov families, and that taking $k=1,...,n-1$ in the spatial case gives the $n-1$ vertical Lyapunov families.

We stress that the $G$-equivariant degree theory provides only an alternative concerning the global behavior of the Lyapunov families. Without additional information we do not know what actually happens along a given branch. This question is considered in [21], where the authors conduct a numerical exploration of the global behavior of the Lyapunov families using the software package AUTO (e.g. see [22]).

Let $p,q\in\mathbb{Z}$ be relatively prime such that $kq-p\in n\mathbb{Z}$. It is proved in [21] that an orbit with the symmetries defined in equations (2) and (3), and frequency

[TABLE]

is a simple choreography when converted back to the inertial reference frame. In the case that $p$ and $q$ do not satisfy this diophantine equation, the solution in the inertial frame corresponds to a multiple choreographic solution [8], while the case that $\omega/\sqrt{s_{1}}$ is irrational implies that the solution is quasiperiodic. Since the set of rational numbers $p/q$ satisfying the diophantine relation (4) is dense, one has the following: when the frequency $\omega$ varies continuously along the Lyapunov family, there are infinitely many orbits in the rotating frame that correspond to simple choreographies in the inertial frame.

The authors of [21] give compelling numerical evidence which suggests that an axial family of solutions appears after a symmetry-breaking bifurcation from the vertical Lyapunov family in the rotating $n$-body problem. The numerics suggest that this axial family has the symmetries of equations (2) and (3). It is shown further in the same reference that, if the hypothesized axial family exists, then orbits in this family correspond to choreographies in the inertial frame which wind $p$ and $q$ times around the generators of a $2$-torus. That is, the periodic orbits in this alleged axial family give rise to $(p,q)$-torus knot choreographies for the $n$-body problem.

A more refined description of our contribution is that we prove the existence of this axial family. Using the symmetries (2, 3) in Newton’s laws we reduced the equations of motion to a single equation describing the motion of the $n$-th body $u_{n}=(w,z)\in\mathbb{C\times R}$. The equation is a delay differential equation (DDE) with multiple constant delays. More explicitly, we have

[TABLE]

For any number of bodies, these reduced equations (5) represents a system of six scalar equations with multiple constant delays.

Our computer assisted arguments are in the functional analytic tradition of Lanford, Eckamnn, Koch, and Wittwer [23, 24, 25, 26], and build heavily on the earlier work of [27, 28, 29] on DDEs. More precisely, we formulate the existence proofs on a Banach space of Fourier coefficient sequences. The delay operator acts as a multiplicative (diagonal) operator in Fourier coefficient space, and the regularity of periodic solutions translates into rapid decay of the Fourier coefficients. Indeed, as was shown in [30], a periodic solution of a delay differential equation with analytic nonlinearity is analytic when the delays are constant. Then we know a-priori that the Fourier coefficients of a periodic solution of Equations (5) decay exponentially fast.

An important feature of Equations (5) is the conservation of energy, which allows us to fix a desired frequency for the periodic solution a-priori. This reduction greatly simplifies the analysis of the delay differential equation in Fourier space, but requires adding an unfolding parameter to balance the system. In addition we utilize automatic differentiation as in [31, 32, 33], and reformulate (5) as a problem with polynomial nonlinearities. The polynomial problem is amenable to straight forward analysis exploiting the Banach algebra properties of the solution space and we use the FFT algorithm as in [34]. The cost of this simplification is that each additional body augments the system with a single additional scalar equation and a single additional unfolding parameter. Finally we validate the existence of solutions by means of a Newton-Kantorovich argument exploiting the radii polynomial approach as in [35].

We conclude this introduction by mentioning some interesting problems for future study. The zero finding problem studied in the present work is amenable to validated continuation techniques as discussed in [29, 36, 37, 38]. A follow up study will investigate global properties of continuous families of spatial choreographies in the $n$ body problem, and study bifurcations encountered along the branches. In this way we hope to prove for example the conjecture of Marchal/Chenciner [17] that the Lagrange triangle is connected with the figure-eigth choreography trough Marchal’s P-12 family [39]. We also remark that all the choreographies shown to exist in the present work are unstable. Actually, the only known stable choreographies are close to the figure eight for $n=3$. Stability of torus knots in the $n=3$ is being investigated in a forthcoming paper.

Let us also mention that the procedure developed in this paper could be adapted to prove existence of asymmetric planar or spatial choreographies. These choreographies do appear in dense sets of symmetry-breaking families from planar and spatial Lyapunov families. Furthermore, this procedure could be adapted to study choreography solutions in problems with other potentials, such as $r^{-\alpha}$ (with $\alpha=1$ being the gravitational case, $\alpha<1$ the weak force case, and $\alpha>2$ the strong force case). It could also be adapted to Hamiltonian systems with different radial potentials, as long as the polynomial embedding (see Section 2.3) can be done. An interesting problem would be to adapt the method to validate choreographies in families that bifurcate from the polygonal equilibrium in DNLS equations [40] or the $n$-vortex problem on the plane, disk, or sphere [41].

Remark 1** (CAPs in celestial mechanics and dynamics of DDEs).**

Numerical calculations have been central to the development of celestial mechanics since the late Nineteenth and early Twentieth centuries. The reader interested in historical developments before the age of the digital computer can consult the works of George Darwin, Francis Ray Moulton, and the group in Copenhagen led by Elis Strömgren [42, 43, 44]. Problems in celestial navigation and orbit design helped drive the explosion of scientific computing during the space race of the mid Twentieth century. A fascinating account and a much more complete bibliography are found in the book [45].**

As researchers developed computer assisted methods of proof for computational dynamics it was natural to look for challenging open problems in celestial mechanics. The relevant literature is rich and we direct the interested reader to the works of [46, 47, 16, 48, 49] for a much more complete view of the literature. Other authors have studied center manifolds [50], transverse intersections of stable/unstable manifolds [31, 51], Melnikov theory [52], Arnold diffusion and transport [53, 54, 55], and existence/continuation/bifurcation of Halo orbits [32, 56] – all in gravitational $n$-body problems and all using computer assisted arguments. Especially relevant to the present work are the computer assisted existence and KAM stability proofs for $n$-body choreographies in [13, 14, 15, 16]. (See also Remark 2 below). Again, the references given in the preceding paragraph are meant only to point the reader in the direction of the relevant literature. A more complete view of the literature is found in the references of the cited works.**

*The present work grows out of the existing literature on CAPs for dynamics of DDEs, the foundations of which were laid in [27]. The work just cited studied periodic solutions – as well as branches of periodic solutions – for scalar DDEs with a single delay and polynomial nonlinearities. Extensions to multiple delays appear in [28], and more recent work considers systems of DDEs with non-polynomial nonlinearities [33]. The interested reader can consult the works of [29, 57, 58, 59] for more complete discussion of this area. We mention also the recent Ph.D. thesis of Jonathan Jaquette, who settled the decades old conjectures of Wright and Jones about the global dynamics of Wright’s equation [60, 61] using ideas from this field. Another approach to computer assisted proof for periodic orbits of DDEs – based on rigorous integration of the induced flow in function space – is found in [62]. ***

In spite of the picture painted above, computer assisted methods of proof are regularly applied outside the boundaries of celestial mechanics and delay differential equations. For a broader perspective on the area, still focusing on nonlinear dynamics, we refer to the review articles [63, 64] and to the book of Tucker [65].**

Remark 2** (Phase space and functional analytic approaches).**

The existence proofs for planar choreographies in [13] and [15], the proof of the spatial mountain pass solution in [16], and the proof of KAM stability of the figure eight choreography in [14] use a different setup from that developed in the present work. More precisely, the works just mentioned study directly the Newtonian equations of motion in phase space. The works of [13, 15, 14] exploit the powerful CAPD library for rigorous integration of ODEs to construct mathematically rigorous arguments in appropriate Poincaré sections. See [66, 67] for more complete discussion and references to the CAPD library. The work of [16] utilizes a functional analytic method akin to that of the present work, but applied directly to periodic orbits for the Hamiltonian vector field rather than reducing to the delay differential equation as in the present work.

---

*In the case of the planar choreography problem the phase space is of dimension $4n$, while the spatial choreography problem scales like $6n$. These figures are in some sense conservative, as applying the topological arguments of [13, 15] require integration of the equations of first variation (and equations of higher variation in the case of the KAM stability argument).

---

The setup of the present work considers six scalar equations, independent of the number of bodies considered. This is a dramatic reduction of the dimension of the problem. This dimension reduction facilitates consideration of – in principle – choreographies involving any number of bodies. A technical remark is that our implementation uses automatic differentiation to reduce to a polynomial nonlinearity, adding one additional scalar equation for each body being considered. This brings our count to $6+(n-1)$ scalar equations. While this quantity scales with $n$ much better than the $6n$ mentioned above, we stress that our implementation could be improved using techniques similar to those discussed in [16, 68, 69] for evaluation of non-polynomial nonlinearities on Fourier data. With such an improvement our approach would consider only 6 scalar equations no matter the number of bodies.**

*For the sake of simplicity we do not pursue this option at the present time, as we believe that the reduction to a polynomial nonlinearity makes both the presentation and implementation of the method more transparent. We also believe that the polynomial version of the problem is more amenable to high order branch following methods and bifurcation analysis to be pursued in a future work. We remark that, since we work in a space of analytic functions, our argument produces useful by-products such as bounds on coefficient decay rates, and lower bounds on the domain of analyticity/bounds on the distances to poles in the complex plane. This information can be used to obtain a-posteriori bounds on derivatives via the usual Cauchy bounds of complex analysis. ***

The paper is organized as follows. In Section 2, we introduce the Fourier map $F:X\to Y$ defined on a Banach space $X$ of geometrically decaying Fourier coefficients, whose zeros are choreographies having prescribed symmetries and topological properties. In Section 3, we introduce the ideas of the a-posteriori validation for the Fourier map, that is on how to demonstrate the existence of true solutions of $F(x)=0$ close to numerical approximations. In Section 4, we present explicit formulas for the bounds necessary to apply the a-posteriori validation of Section 3. We conclude the paper by presenting the results in Section 5, where we present proofs of existence of some spatial torus knot choreographies for $n=4,5,7$, and $9$ bodies. The computer programs used in the paper are available at [70].

2 Formulation of the problem

Let $q_{j}(t)\in\mathbb{R}^{3}$ be the position in space of the body $j\in\{1,...,n\}$ with mass $1$ at this $t$. Define the matrices

[TABLE]

where $J$ is the symplectic matrix in $\mathbb{R}^{2}$. In rotating coordinates and with the period rescaled to $2\pi$,

[TABLE]

the Newton equations for the $n$ bodies are

[TABLE]

where $\omega$ is the frequency and $s_{1}$ is defined by (1).

Using that $u_{j}=(w_{j},z_{j})$, the symmetries (2) and (3) correspond to the symmetry

[TABLE]

Therefore, the solutions of the equation (6) with symmetries (7) are zeros of the map

[TABLE]

defined in spaces $X$ and $Y$ of analytic $2\pi$-periodic functions, which we will specify later in Fourier components. The equation $\mathcal{G}(u_{n},\omega)=0$, with $\mathcal{G}$ defined in (8) is a delay differential equation (DDE).

2.1 Choreographies

We say that a solution of $\mathcal{G}(u_{n},\omega)=0$, i.e. a solution of the $n$-body problem with symmetry (7), is $p:q$* resonant *when it has frequency $\omega=\sqrt{s_{1}}p/q$ and (a) $kq-p=0$ or (b) $p$ and $q$ are relatively prime and $kq-p\in n\mathbb{Z}$. In [21] is proven that $p:q$ resonant orbits are choreographies in the inertial frame, see also [17]. For sake of completeness, here we reproduce a short version of this result.

Proposition 3**.**

Let

[TABLE]

be a reparameterization of a periodic solution in the inertial frame. An $p:q$ resonant solution u_{n}\of $\mathcal{G}(u_{n},\omega)=0$ is a choreography in inertial frame, satisfying that $Q_{n}(t)$ is $2\pi p$-periodic and

[TABLE]

where $\tilde{k}=k-(kq-p)\tilde{q}$ with $\tilde{q}$ the $p$-modular inverse of $q$. The orbit of the choreography is symmetric with respect to rotations by an angle $2\pi/p$ and the $n$ bodies form groups of $h$-polygons, where $h$ is the biggest common divisor of $n$ and $k$.

Proof.

Since $u_{n}(t)$ is $2\pi$-periodic and $e^{t\bar{J}\sqrt{s_{1}}/\omega}=e^{t\bar{J}q/p}$ is $2\pi p$-periodic, then the function $Q_{n}(t)=e^{t\bar{J}\sqrt{s_{1}}/\omega}u_{n}(t)$ is $2\pi p$-periodic. Furthermore, since

[TABLE]

the orbit of $Q_{n}(t)$ is invariant under rotations of $2\pi/p$. The fact that the $n$ bodies form $h$-polygons follows from symmetry (7) and the definition of $Q_{j}(t)$.

By assumption

[TABLE]

then symmetry (7) implies that the solution in inertial frame satisfies

[TABLE]

In the case (a) that $kq-p=0$, the symmetry (9) gives straightforward that $Q_{j}(t)=Q_{n}(t+kj\zeta)$. In the case (b) that $p$ and $q$ are relatively prime, we can find $\tilde{q}$ such that $q\tilde{q}=1$ mod $p$. It follows from the symmetry (9) that

[TABLE]

Therefore,

[TABLE]

Corollary 4** ($(p,q)$-torus knots).**

In the case that $u_{n}(t)$ is a $p:q$ resonant orbit in the axial family that does not cross the $z$-axis, then $Q_{n}(t)$ winds (after the period $2\pi p$) around a toroidal manifold with winding numbers $p$ and $q$, i.e., the choreography path is a $(p,q)$-torus knot. In the case that $u_{n}(t)$ is a $p:q$ resonant orbit in the vertical Lyapunov family that does not cross the $z$-axis, then the choreography $Q_{n}(t)$ winds $p$ times in a cylindrical surface.

We conclude that the solution $q_{j}(t)=Q_{j}(t\omega)$ is a $2\pi q/\sqrt{s_{1}}$-periodic choreography satisfying the properties discussed above for $Q_{j}(t)$. Therefore, by validating solutions of $\mathcal{G}(u_{n},\omega)=0$ in the axial family we prove rigorously the existence of choreography paths that are $(p,q)$-torus knots.

2.2 Symmetries, integrals of movement and Poincaré conditions

Here after we omit the index $n$ that represents the $n$th body in the map $\mathcal{G}(u)$ and denote the components of $u$ by

[TABLE]

The map $\mathcal{G}(u)$ that gives the existence of choreographies is the gradient of the action $\mathcal{A}(u):X\rightarrow\mathbb{R}$ of the $n$-body problem reduced to paths with symmetries (7). The action $\mathcal{A}(u)$ is invariant under the action of the group $(\theta,\varphi,\tau)\in G\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,T^{2}\times\mathbb{R}$ in $u\in X$ given by

[TABLE]

which corresponds to $z$-translations and $(x,y)$-rotations of bodies, and time shift.

Given that the gradient $\mathcal{G}=\nabla\mathcal{A}$ is $G$-equivariant, $\mathcal{G}((\theta,\varphi,\tau)u)=(\theta,\varphi,\tau)\mathcal{G}(u)$, if $u_{0}$ is a critical point of $\mathcal{A}$, then $(\theta,\varphi,\tau)u_{0}$ is a critical point for all $(\theta,\varphi,\tau)\in G$, because

[TABLE]

Therefore, if $u_{0}$ is not fixed by the elements of $G$, then its orbit under the action of the group forms a $3$-dimensional manifold of zeros of $\mathcal{G}$. Taking derivatives respect the parameters $\theta,$ $\varphi$ and $\tau$ of equation (11) and evaluating the parameter at [math], we obtain by the chain rule that $D\mathcal{G}(u_{0})A_{j}(u_{0})=0$, where $A_{j}$ are the generator fields of the group $G$,

[TABLE]

Therefore $D\mathcal{G}(u_{0})$ has the zero eigenvalues $A_{j}(u_{0})$ for $j=1,2,3$ corresponding to tangent vectors to the $3$-dimensional manifold generated by the action of $G$. This property holds for any equivariant field even if it is not gradient.

In addition, for gradient maps $\mathcal{G}=\nabla\mathcal{A}$, we have also conserved quantities generated by the action of the group $G$ (Noether theorem). That is, since the action is invariant, $\mathcal{A}((\theta,\varphi,\tau)u)=\mathcal{A}(u)$, deriving respect $\theta,$ $\varphi$ and $\tau$ and evaluating the parameters at [math], we have by chain rule that

[TABLE]

i.e. the field $\mathcal{G}$ is orthogonal to the infinitesimal generators $A_{j}(u)$ for $j=1,2,3$.

In summary, we have that the map $\mathcal{G}$ has $3$-dimensional families of zeros and also $3$-restrictions given by (12). To prove the existence of solutions, we could take $3$-restrictions in the domain and range of $\mathcal{G}$. But given that the range is a non-flat manifold, it is simpler to augment the delay differential equation $\mathcal{G}=0$ with the three Lagrangian multipliers $\lambda_{j}$ for $j=1,2,3$,

[TABLE]

An important observation is that the solutions of equation (13) are equivalent to the solutions of the original equations of motion.

Proposition 5**.**

If $A_{j}(u)$ are linearly independent for $j=1,2,3$, then a solution $u$ to $\mathcal{G}(u,\omega)=0$ is a solution to the equation (13) if and only if $\lambda_{j}=0$ for $j=1,2,3$.

Proof.

Taking the product of (13) with respect to a generator $A_{j}(u)$ and using the orthogonality we obtain

[TABLE]

The result follows from the linear independence of $A_{j}(u)$, see [21] for details. ∎

Also the restriction in the domain forms a non-flat manifold, and it is simpler to augment the equation (13) with three equations that represent the respective Poincaré sections $I_{j}(u)=0$. Each geometric condition $I_{j}(u)=0$ with

[TABLE]

implies that $u$ is in the orthogonal plane to the orbit of $\tilde{u}$ under the action of $G$, where $\tilde{u}$ is a reference solution, which typically is the solution in the previous step of the continuation.

Taking as reference $\tilde{u}=(1,0,0)$ for the generators $A_{3}(\tilde{u})=(0,0,1)$, then

[TABLE]

Given a reference solution $\tilde{u}$, the other geometric conditions are given explicitly by

[TABLE]

and

[TABLE]

The generators $A_{j}(u)$ are linearly independent in the solutions that we are looking. In other cases the solutions are relative equilibria, which represents a simpler problem than the map $\mathcal{G}$.

2.3 Automatic differentiation: obtaining a polynomial problem

Setting $\dot{u}=v$, equation $\mathcal{G}(u,\omega)=0$ becomes

[TABLE]

In this section, we turn the non-polynomial DDE (13) into a higher dimensional DDE with polynomial nonlinearities, using the automatic differentiation technique as in [31, 32, 33]. For this, we define for $j=1,...,n-1$ the variables

[TABLE]

Then $w_{j}$ satisfy

[TABLE]

Therefore, the augmented system of equations (13) is

[TABLE]

for $j=1,\dots,n-1$, where $e_{3}=(0,0,1)$. We supplement these equations with the conditions

[TABLE]

which are balanced by the unfolding parameters $\alpha_{1},\dots,\alpha_{n-1}$ (e.g. see [32]), similarly to the manner in which the phase conditions $I_{1}(u)=I_{2}(u)=I_{3}(u)=0$ (given respectively by (15), (16) and (14)) are balanced by the unfolding parameters $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$. Indeed, we can prove that a solution of this system is necessarily a solution of the $n$-body problem similarly to Proposition 5.

Proposition 6**.**

A $2\pi$-periodic solution $(u,v,w)$ of the system (17)-(19) with the conditions (20) satisfies that $\alpha_{j}=0$ for $j=1,...,n$, i.e. $u$ is a $2\pi$-periodic solution of $\mathcal{G}(u,\omega)=0$.

Proof.

Dividing the equation for $w_{j}$ by $w_{j}^{3}$ and using that $v=\dot{u}$, we obtain that

[TABLE]

Since $(u,v,w)$ is $2\pi$-periodic, integrating over the period $2\pi$, we obtain that $2\pi\alpha_{j}=0$, see [32] for details. Given that $\alpha_{j}=0$, the initial condition (20) implies that $w_{j}(t)=\left\|u(t)-e^{\ j\bar{J}\zeta}u(t+jk\zeta)\right\|^{-1}$. Therefore, $u$ is a solution to the augmented system (13) and, by Proposition 5, to the equation $\mathcal{G}(u,\omega)=0$. ∎

In the next section, equations (17), (18), (19) and (20) are combined with Fourier expansions to set up the Fourier map whose zeros corresponds to choreographies having the prescribed symmetry (7) and the topological property of a torus knot.

2.4 Fourier map for automatic differentiation

The goal of this section is to look for periodic solutions of the delay differential equations (17), (18) and (19) satisfying the extra conditions (20) using the Fourier series expansions

[TABLE]

Based on the fact that periodic solutions of analytic DDEs are analytic [30], we consider the following Banach space of geometrically decaying Fourier coefficients

[TABLE]

where $\nu\geq 1$. If $\nu>1$ and $a=(a_{\ell})_{\ell\in\mathbb{Z}}\in\ell_{\nu}^{1}$, then the function $t\mapsto\sum_{\ell\in\mathbb{Z}}e^{i\ell t}a_{\ell}$ defines a $2\pi$-periodic analytic function on the complex strip of width $\ln(\nu)>0$. Another useful property of the space $\ell_{\nu}^{1}$ is that it is a Banach algebra under discrete convolution $*:\ell_{\nu}^{1}\times\ell_{\nu}^{1}\to\ell_{\nu}^{1}$ defined as

[TABLE]

where $a,b\in\ell_{\nu}^{1}$. More explicitly, $\|a*b\|_{\nu}\leq\|a\|_{\nu}\|b\|_{\nu}$, for all $a,b\in\ell_{\nu}^{1}$ and $\nu\geq 1$.

The unknowns of the DDEs (17), (18) and (19) are given by the unfolding parameters $\lambda\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,(\lambda_{j})_{j=1}^{3}\in\mathbb{C}^{3}$ and $\alpha\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,(\alpha_{j})_{j=1}^{n-1}\in\mathbb{C}^{n-1}$, and the Fourier coefficients $u=(u_{j})_{j=1}^{3}\in(\ell_{\nu}^{1})^{3}$, $v=(v_{j})_{j=1}^{3}\in(\ell_{\nu}^{1})^{3}$ and $w=(w_{j})_{j=1}^{n-1}\in(\ell_{\nu}^{1})^{n-1}$. The total vector of unknown $x$ and the Banach space $X$ are then given by

[TABLE]

The Banach space $X$ is endowed with the norm

[TABLE]

where

[TABLE]

In order to define the Fourier map problem $F(x)=0$, we plug the Fourier expansions (21) in (17), (18), (19) and (20), and solve for the corresponding nonlinear map. First note that

[TABLE]

where $M_{j\ell}$ is defined as

[TABLE]

since $\bar{J}=diag(J,0)$ with $J=\begin{pmatrix}0&-1\\ 1&0\end{pmatrix}$.

In Fourier space, the phase conditions $I_{1}(u)=I_{2}(u)=I_{3}(u)=0$ (see (15), (16) and (14), respectively) are given by

[TABLE]

where $\tilde{u}_{1}$, $\tilde{u}_{2}$ and $\tilde{u}_{3}$ have only finitely many non zero terms.

Hence, setting $\eta:(\ell_{\nu}^{1})^{3}\to\mathbb{C}^{3}$ as

[TABLE]

we get that $\eta(u)=0$ implies that $I_{1}(u)=I_{2}(u)=I_{3}(u)=0$. Given $j=1,\dots,n-1$ and $u\in(\ell_{\nu}^{1})^{3}$, denote $M_{j}u\in(\ell_{\nu}^{1})^{3}$ component-wise by

[TABLE]

In Fourier space, the extra initial condition (20) (given $j=1,\dots,n-1$) is simplified as

[TABLE]

Set $\gamma:(\ell_{\nu}^{1})^{3}\times(\ell_{\nu}^{1})^{n-1}\to\mathbb{C}^{n-1}$ as

[TABLE]

Hence, $\gamma(u,w)=0$ implies that (20) holds.

For sake of simplicity of the presentation, given any $N\in\mathbb{N}$, denote the differentiation operator $D$ acting on $u\in(\ell_{\nu}^{1})^{N}$ as

[TABLE]

Remark 7**.**

The linear operator $D$ is not bounded on $(\ell_{\nu}^{1})^{N}$. However, it is bounded when considering the image to be slightly less regular. More explicitly, letting

[TABLE]

we can easily verify that $D:(\ell_{\nu}^{1})^{N}\to(\tilde{\ell}_{\nu}^{1})^{N}$ is a bounded linear operator.

Let $f:(\ell_{\nu}^{1})^{3}\times(\ell_{\nu}^{1})^{3}\to(\tilde{\ell}_{\nu}^{1})^{3}$ be defined by

[TABLE]

Note that $f(u,v)=0$ ensures that (17) holds. Let $g:\mathbb{C}^{3}\times(\ell_{\nu}^{1})^{3}\times(\ell_{\nu}^{1})^{3}\times(\ell_{\nu}^{1})^{n-1}\times\mathbb{C}\to(\tilde{\ell}_{\nu}^{1})^{3}$ be defined by

[TABLE]

where $(M_{j}u)*w_{j}^{3}\in(\ell_{\nu}^{1})^{3}$ is given component-wise by

[TABLE]

and where $\hat{e}_{3}\in(\ell_{\nu}^{1})^{3}$ is given component-wise by

[TABLE]

with $\delta_{i,j}$ being the Kronecker delta. Note that $g(\lambda,u,v,w,\omega)=0$ ensures that (18) holds.

Let $h_{j}:\mathbb{C}\times(\ell_{\nu}^{1})^{3}\times(\ell_{\nu}^{1})^{3}\times\ell_{\nu}^{1}\to\tilde{\ell}_{\nu}^{1}$ be defined by

[TABLE]

and let $h:\mathbb{C}^{n-1}\times(\ell_{\nu}^{1})^{3}\times(\ell_{\nu}^{1})^{3}\times(\ell_{\nu}^{1})^{n-1}\to(\tilde{\ell}_{\nu}^{1})^{n-1}$ be defined by

[TABLE]

Hence, $h(\alpha,u,v,w)=0$ implies that (19) hold.

Defining

[TABLE]

the Fourier map $F:X\times\mathbb{R}\to Y$ is defined by

[TABLE]

For a fixed $\omega>0$, we introduce in Section 3 an a-posteriori validation method for the Fourier map, that is we develop a systematic and constructive approach to prove existence of $x\in X$ such that $F(x,\omega)=0$. By construction, the solution $x$ yields a choreography having the prescribed symmetry (7) and the topological property of a torus knot.

3 A-posteriori validation for the Fourier map

The idea of the computer-assisted proof of existence of a spatial torus-knot choreography is to demonstrate that a certain Newton-like operator is a contraction on a closed ball centered at a numerical approximation ${\bar{x}}$. To compute ${\bar{x}}$, we consider a finite dimensional projection of the Fourier map $F:X\times\mathbb{R}\to Y$. Given a number $m\in\mathbb{N}$, and given a vector $a=(a_{\ell})_{\ell\in\mathbb{Z}}\in\ell_{\nu}^{1}$, consider the projection

[TABLE]

We generalize that projection to get $\pi_{N}^{m}:(\ell_{\nu}^{1})^{N}\to\mathbb{C}^{N(2m-1)}$ defined by

[TABLE]

and $\Pi^{(m)}:X\to\mathbb{C}^{2m(n+5)-3}$ defined by

[TABLE]

Often, given $x\in X$, we denote

[TABLE]

Moreover, we define the natural inclusion $\iota^{m}:\mathbb{C}^{2m-1}\xhookrightarrow{}\ell_{\nu}^{1}$ as follows. For $a=(a_{\ell})_{|\ell|<m}\in\mathbb{C}^{2m-1}$ let $\iota^{m}a\in\ell_{\nu}^{1}$ be defined component-wise by

[TABLE]

Similarly, let $\iota_{N}^{m}:\mathbb{C}^{N(2m-1)}\xhookrightarrow{}(\ell_{\nu}^{1})^{N}$ be the natural inclusion defined as follows. Given $a=(a^{(1)},\dots,a^{(N)})\in(\mathbb{C}^{2m-1})^{N}\cong\mathbb{C}^{N(2m-1)}$,

[TABLE]

Finally, let the natural inclusion $\bm{\iota}^{(m)}:\mathbb{C}^{2m(n+5)-3}\xhookrightarrow{}X$ be defined, for $x\in\mathbb{C}^{2m(n+5)-3}$ as

[TABLE]

Finally, let the finite dimensional projection $F^{(m)}:\mathbb{C}^{2m(n+5)-3}\to\mathbb{C}^{2m(n+5)-3}$ of the Fourier map be defined, for $x\in\mathbb{C}^{2m(n+5)-3}$, as

[TABLE]

Also denote $F^{(m)}=\left(\eta^{(m)},\gamma^{(m)},f^{(m)},g^{(m)},h^{(m)}\right)$.

Assume that, using Newton’s method, a numerical approximation ${\bar{x}}\in\mathbb{C}^{2m(n+5)-3}$ of (35) has been obtained at a parameter (frequency) value $\omega$, that is $F^{(m)}({\bar{x}},\omega)\approx 0$. We slightly abuse the notation and denote ${\bar{x}}\in\mathbb{C}^{2m(n+5)-3}$ and $\bm{\iota}^{(m)}{\bar{x}}\in X$ both using ${\bar{x}}$.

We now fix an $\omega_{0}\in\mathbb{R}$ and consider the mapping $F\colon X\to Y$ defined by $F(x)=F(x,\omega_{0})$. The following result is a Newton-Kantorovich theorem with a smoothing approximate inverse. It provides an a-posteriori validation method for proving rigorously the existence of a point ${\tilde{x}}$ such that $F({\tilde{x}})=0$ and $\|{\tilde{x}}-{\bar{x}}\|_{X}\leq r$ for a small radius $r$. Recalling the norm on $X$ given in (24), denote by

[TABLE]

the ball of radius $r$ centered at $y\in X$.

Theorem 8** **(Radii Polynomial Approach).

For $\bar{x}\in X$ and $r>0$ assume that $F:X\to Y$ is Fréchet differentiable on the ball $B_{r}(\bar{x})$. Consider bounded linear operators $A^{\dagger}\in B(X,Y)$ and $A\in B(Y,X)$, where $A^{\dagger}$ is an approximation of $DF({\bar{x}})$ and $A$ is an approximate inverse of $DF({\bar{x}})$. Observe that

[TABLE]

Assume that $A$ is injective. Let $Y_{0},Z_{0},Z_{1},Z_{2}\geq 0$ be bounds satisfying

[TABLE]

Define the radii polynomial

[TABLE]

If there exists $0<r_{0}\leq r$ such that

[TABLE]

then there exists a unique ${\tilde{x}}\in B_{r_{0}}({\bar{x}})$ such that $F({\tilde{x}})=0$.

Proof.

Details of the elementary proof are found in Appendix A of [71]. The idea is to first show that $T(x)\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,x-AF(x)$ satisfies $T(B_{r_{0}}({\bar{x}}))\subset B_{r_{0}}({\bar{x}})$, and then to show the existence of $\kappa<1$ such that $\|T(x)-T(y)\|_{X}\leq\kappa\|x-y\|_{X}$ for all $x,y\in B_{r_{0}}({\bar{x}})$. These facts follow from the inequalities of Equations (37), (38), (39), (40), and from the hypothesis that $p(r_{0})<0$. The proof then follows from the contraction mapping theorem and the injectivity of $A$. ∎

The following corollary provides an additional useful byproduct.

Corollary 9** **(Non-degeneracy at the true solution).

Given the hypotheses of Theorem 8, the linear operator $ADF(\tilde{x})$ is boundedly invertible with

[TABLE]

Proof.

From

[TABLE]

we obtain

[TABLE]

or

[TABLE]

Since $Y_{0}$ and $r_{0}$ are both positive it follows that

[TABLE]

Since $\tilde{x}\in B_{r_{0}}(\bar{x})$ we have that $\tilde{x}=\bar{x}+b$ for some $b\in B_{r_{0}}(0)$, and by applying the inequalities of Equations (38), (39), and (40) we have that

[TABLE]

Then

[TABLE]

is invertible by the Neuman theorem and

[TABLE]

as desired. ∎

Returning to the parameter dependent problem, suppose that $\tilde{x}$ is a zero of $F(x)=F(x,\omega_{0})$ and that $ADF(\tilde{x})=AD_{x}F(\tilde{x},\omega_{0})$ is boundedly invertible as above. Notice that $F(x,\omega)$ is differentiable with respect to $\omega$ near $\omega_{0}$. Define the mapping $G(x,\omega)=AF(x,\omega)$ and observe that $G$ and $F$ have the same zero set as $A$ is injective. Observe also that $D_{x}G(x,\omega)=AD_{x}F(x,\omega)$. So $(\tilde{x},\omega_{0})$ is a zero of $G$ with $D_{x}G(\tilde{x},\omega_{0})$ an isomorphism, it follows from the implicit function theorem that $G$ has a smooth branch of zeros through $\tilde{x}$. More precisely there exists an $\epsilon>0$ and a smooth function $x\colon(\omega_{0}-\epsilon,\omega_{0}+\epsilon)\to X$ with $x(\omega_{0})=\tilde{x}$ and

[TABLE]

for all $\omega\in(\omega_{0}-\epsilon,\omega_{0}+\epsilon)$. It follows again from the injectivity of $A$ that $F(x(\omega),\omega)=0$ for all $\omega\in(\omega_{0}-\epsilon,\omega_{0}+\epsilon)$. Finally, as discussed in the introduction, we obtain that for any rational number $\sqrt{s_{1}}p/q\in(\omega_{0}-\epsilon,\omega_{0}+\epsilon)$ , the solution $x(\sqrt{s_{1}}p/q)$ produces spatial torus knot choreography orbit near $\tilde{x}$ by proposition 4. Taken together the results of this section show that our method produces the existence of countably many spatial torus knot choreographies as soon as Theorem 8 succeeds at a given $\omega_{0}$.

3.1 Isolated solutions yield real periodic solutions

In this short section, we show how the output ${\tilde{x}}\in B_{r_{0}}({\bar{x}})$ of Theorem 8 (if any) yields a real periodic solution, provided the numerical approximation is chosen to represent a real periodic solution.

Define the operator $\sigma:\ell_{\nu}^{1}\to\ell_{\nu}^{1}$ by $(\sigma(a))_{\ell}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,a_{-\ell}^{*},\quad\text{for }\ell\in\mathbb{Z}$, where $z^{*}$ denotes the complex conjugate of $z\in\mathbb{C}$. Define the symmetry subspace $\ell_{\nu}^{1,\real}\subset\ell_{\nu}^{1}$ by

[TABLE]

Note that if $(u_{\ell})_{\ell\in\mathbb{Z}}\in\ell_{\nu}^{1,\real}$, then the function $u(t)\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\sum_{\ell\in\mathbb{Z}}u_{\ell}e^{i\ell t}$ is a real $2\pi$-periodic function. Define the operator $\Sigma:X\to X$ acting on $x=(\lambda,\alpha,u,v,w)\in X$ as

[TABLE]

where $\lambda^{*}\in\mathbb{C}^{3}$ and $\alpha^{*}\in\mathbb{C}^{n-1}$ denote the component-wise complex conjugate of $\lambda\in\mathbb{C}^{3}$ and $\alpha\in\mathbb{C}^{n-1}$, respectively. Define the subspace $X\subset X$ as

[TABLE]

It follows by definition that $X=\mathbb{R}^{n+2}\times(\ell_{\nu}^{1,\real})^{n+5}$.

Proposition 10**.**

Fix a frequency $\omega>0$ and assume that the numerical approximation denoted ${\bar{x}}=(\bar{\lambda},\bar{\alpha},{\bar{u}},{\bar{v}},{\bar{w}})$ satisfies ${\bar{x}}\in X$ and that the reference solution ${\tilde{u}}=({\tilde{u}}_{1},{\tilde{u}}_{2},{\tilde{u}}_{3})$ satisfies ${\tilde{u}}\in(\ell_{\nu}^{1,\real})^{3}$. Assume that there exists a unique $x\in B_{r}({\bar{x}})$ such that $F(x,\omega)=0$. Then $x\in X$.

Proof.

Denote the solution $x=(\lambda,\alpha,u,v,w)\in B_{r}({\bar{x}})$. The proof is twofold: (1) show that $F(\Sigma(x),\omega)=0$; and (2) show that $\Sigma(x)\in B_{r}({\bar{x}})$. The conclusion $\Sigma(x)=x$ (that is $x\in X$) then follows by unicity of the solution. First, we have that $F(\Sigma(x),\omega)=\Sigma\left(F(x,\omega)\right)$, since the operator $F$ corresponds to the complex extension of a real equation. Since $F(x,\omega)=0$, then $F(\Sigma(x),\omega)=\Sigma\left(F(x,\omega)\right)=\Sigma\left(0\right)=0$. Second, to prove that $\Sigma(x)\in B_{r}({\bar{x}})$, it is sufficient to realize that $|z^{*}|=|z|$ and that given any $c\in\ell_{\nu}^{1}$,

[TABLE]

which shows that for any $\xi\in X$, $\|\Sigma(\xi)\|_{X}=\|\xi\|_{X}$. Hence, since $\Sigma({\bar{x}})={\bar{x}}$, we conclude that

[TABLE]

3.2 Definition of the operators $A^{\dagger}$ and $A$

To apply the radii polynomial approach of Theorem 8, we need to define the approximate derivative $A^{\dagger}$ and the smoothing approximate inverse $A$. Consider the finite dimensional projection $F^{(m)}:\mathbb{C}^{2m(n+5)-3}\to\mathbb{C}^{2m(n+5)-3}$ and assume that at a fixed frequency $\omega>0$ we computed ${\bar{x}}\in\mathbb{C}^{2m(n+5)-3}$ such that $F^{(m)}({\bar{x}},\omega)\approx 0$. Denote by $DF^{(m)}({\bar{x}},\omega)\in M_{2m(n+5)-3}(\mathbb{C})$ the Jacobian matrix of $F^{(m)}$ at $({\bar{x}},\omega)$. Given $x\in X$, define

[TABLE]

where $\Pi^{(m)}A^{\dagger}x=DF^{(m)}({\bar{x}},\omega)x^{(m)}$ and

[TABLE]

Recalling the definition of the Banach space $Y$ in (33), we can verify that the operator $A^{\dagger}:X\to Y$ is a bounded linear operator, that is $A^{\dagger}\in B(X,Y)$. For $m$ large enough, it acts as an approximation of the true Fréchet derivative $D_{x}F({\bar{x}},\omega)$. Its action on the finite dimensional projection is the Jacobian matrix (the derivative) of $F^{(m)}$ at $({\bar{x}},\omega)$ while its action on the tail keeps only keep the unbounded terms involving the differentiation $D$ as defined in (27).

Consider now a matrix $A^{(m)}\in M_{2m(n+5)-3}(\mathbb{C})$ computed so that $A^{(m)}\approx{DF^{(m)}({\bar{x}},\omega)}^{-1}$. In other words, this means that $\|I-A^{(m)}DF^{(m)}({\bar{x}},\omega)\|\ll 1$. This step is performed using a numerical software (MATLAB in our case). We decompose the matrix $A^{(m)}$ block-wise as

[TABLE]

so that it acts on $x^{(m)}=(\lambda,\alpha,u^{(m)},v^{(m)},w^{(m)})\in\mathbb{C}^{2m(n+5)-3}$. Thus we define $A$ as

[TABLE]

where the action of each block of $A$ is finite (that is they act on $x^{(m)}=\Pi^{(m)}x$ only) except for the three diagonal blocks $A_{u,u}$, $A_{v,v}$ and $A_{w,w}$ which have infinite tails. More explicitly, for each $p=1,2,3$,

[TABLE]

and for each $j=1,\dots,n-1$,

[TABLE]

Having defined the operators $A$ and $A^{\dagger}$, we are ready to define the bounds $Y_{0}$, $Z_{0}$, $Z_{1}$ and $Z_{2}$ (satisfying (37), (38), (39) and (40), respectively), required to build the radii polynomial defined on (41).

4 The technical estimates for the Fourier map

In this section, we introduce explicit formulas for the theoretical bounds (37), (38), (39) and (40). While most of the work is analytical, the actual definition of the bounds still requires computing and verifying inequalities. In particular, there are many occasions in which the most practical means of obtaining necessary explicit inequalities is by using the computer. However, as floating point arithmetic is only capable of representing a finite set of rational numbers, round off errors in the computation of the bounds can be dealt with by using interval arithmetic [72] where real numbers are represented by intervals bounded by rational numbers that have floating point representation. Furthermore, there is software that performs interval arithmetic (e.g. INTLAB [73]) which we use for completing our computer-assisted proofs. With this in mind, in this section, when using phrases of the form we can compute the following bounds, this should be interpreted as shorthand for the statement using the interval arithmetic software INTLAB we can compute the following bounds.

4.1 $Y_{0}$ bound

Denote the numerical approximation ${\bar{x}}=(\bar{\lambda},\bar{\alpha},{\bar{u}},{\bar{v}},{\bar{w}})\in X$ with ${\bar{u}}=({\bar{u}}_{1},{\bar{u}}_{2},{\bar{u}}_{3})\in(\ell_{\nu}^{1})^{3}$, ${\bar{v}}=({\bar{v}}_{1},{\bar{v}}_{2},{\bar{v}}_{3})\in(\ell_{\nu}^{1})^{3}$ and ${\bar{w}}=({\bar{w}}_{1},\dots,{\bar{w}}_{n-1})\in(\ell_{\nu}^{1})^{n-1}$. Recalling (29), (30) and (31), one has that

[TABLE]

since the product of $p$ trigonometric functions of degree $m-1$ is a trigonometric function of degree $p(m-1)$. For instance, recalling (30), the highest degree terms in $g(\bar{\lambda},{\bar{u}},{\bar{v}},{\bar{w}},\omega)$ are of the form $(M_{j}{\bar{u}})*{\bar{w}}_{j}^{3}$ which are convolutions of degree four, and therefore have zero Fourier coefficients for all frequencies $\ell$ such that $|\ell|>4m-4$. This implies that $F({\bar{x}},\omega)$ has only a finite number of nonzero terms. Hence, we can compute $Y_{0}$ satisfying (37).

4.2 $Z_{0}$ bound

Let $B\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,I-AA^{\dagger}$, which we denote block-wise by

[TABLE]

Note that by definition of the diagonal tails of $A$ and $A^{\dagger}$, the tails of $B$ vanish, that is all $B_{\delta,\tilde{\delta}}$ ($\delta,\tilde{\in}\{u,v,w\}$) are represented by $2m-1\times 2m-1$ matrices. We can compute the bound

[TABLE]

By construction, letting

[TABLE]

we get that

[TABLE]

4.3 $Z_{1}$ bound

Recall from (39) that the $Z_{1}$ bound satisfy

[TABLE]

For the computation of this bound, it is convenient to define, given any $h\in B_{1}(0)\in X$

[TABLE]

Denote

[TABLE]

The construction of $Z_{1}$ hence requires computing an upper bound for $\|Az\|_{X}$ for all $h\in B_{1}(0)\in X$. This is done by splitting $Az$ as

[TABLE]

and by handling each term separately.

Remark 11**.**

We choose the Galerkin projection number $m$ greater than the number $m_{1}$ of nonzero Fourier coefficients of the previous orbit $({\tilde{u}}_{1},{\tilde{u}}_{2},{\tilde{u}}_{3})$. Then $z_{\lambda}=0\in\mathbb{C}^{3}$. This is because the phase conditions $\eta(u)$ defined in (25) only depend on the modes of the finite dimensional approximation and therefore $A^{\dagger}$ contains all contribution from $D\eta({\bar{u}})h$.

As $\Pi^{(m)}Az=A^{(m)}z^{(m)}$, we compute a uniform component-wise upper bound

[TABLE]

for the complex modulus of each component of

[TABLE]

The computation of the bounds $\hat{z}_{\alpha}$, $\hat{z}^{(m)}_{u}$, $\hat{z}^{(m)}_{v}$ and $\hat{z}^{(m)}_{w}$ is done in Sections 4.3.1, 4.3.2, 4.3.3 and 4.3.4, respectively. Using these uniform bounds (i.e. for all $h\in B_{1}(0)$), let

[TABLE]

where the entries of the matrix $|A^{(m)}|$ are the component-wise complex magnitudes of the entries of $A^{(m)}$. By construction, the bound $\xi^{(m)}$ of (50) provides a uniform component-wise upper bound for the first term $\bm{\iota}^{(m)}\Pi^{(m)}Az$ of the splitting (49) of $Az$. To handle the second term $(I-\bm{\iota}^{(m)}\Pi^{(m)})Az$ of (49), we compute the uniform (i.e. for all $h\in B_{1}(0)$) tail bounds $(\delta_{u})_{p}$, $(\delta_{v})_{p}$ (for $p=1,2,3$) and $(\delta_{w})_{j}$ (for $j=1,\dots,n-1$) satisfying

[TABLE]

The computation of the bounds $\delta_{u}$, $\delta_{v}$ and $\delta_{w}$ is presented in Sections 4.3.2, 4.3.3 and 4.3.4, respectively. Combining the above bounds, we get that

[TABLE]

4.3.1 Computation of the bound $\hat{z}_{\alpha}$

Recalling (48) and the definition of $A^{\dagger}$ in (45), one can verify that for any $j=1,\dots,n-1$,

[TABLE]

Straightforward calculations (e.g. using Lemma 2.1 in [35]) involving bounding linear functionals on $\ell_{\nu}^{1}$ and using that $(h_{u})_{p}\in B_{1}(0)\subset\ell_{\nu}^{1}$ for $p=1,2,3$ yield that

[TABLE]

We therefore get the component-wise bound (given $j=1,\dots,n-1$)

[TABLE]

4.3.2 Computation of the bounds $\hat{z}^{(m)}_{u}$ and $\delta_{u}$

From (48) and (45), on can verify that for each $p=1,2,3$,

[TABLE]

Hence, since $z_{u}$ only has a tail and since the blocks $A_{\lambda,u}$, $A_{\alpha,u}$, $A_{v,u}$ and $A_{w,u}$ only acts on the finite part, then $A_{\delta,u}z_{u}=0$ for $\delta=\lambda,\alpha,v,w$ and for $p=1,2,3$

[TABLE]

Now,

[TABLE]

We can then set

[TABLE]

4.3.3 Computation of the bound $\hat{z}^{(m)}_{v}$ and $\delta_{v}$

The following technical lemma (which is a slight modification of Corollary 3 in [35]) is the key to the truncation error analysis of $\hat{z}_{v}^{(m)}$ and $\hat{z}_{w}^{(m)}$.

Lemma 12**.**

Fix a truncation Fourier mode to be $m$. Given $h\in\ell_{\nu}^{1}$, set

[TABLE]

Let $N\in\mathbb{N}$ and let $\bar{\alpha}=(\ldots,0,0,\bar{\alpha}_{-N},\ldots,\bar{\alpha}_{N},0,0,\ldots)\in\ell_{\nu}^{1}$. Then, for all $h\in\ell_{\nu}^{1}$ such that $\|h\|_{\nu}\leq 1$, and for $|\ell|<m$,

[TABLE]

Now, from (48) and (45), on can verify that for each $p=1,2,3$,

[TABLE]

Using Lemma 12, we obtain that for $|\ell|<m$ and $p=1,2,3$,

[TABLE]

which provides a component-wise definition of the vector $\hat{z}_{v}^{(m)}\in\mathbb{R}_{+}^{3(2m-1)}$. Finally, one can verify using the fact that $\ell_{\nu}^{1}$ is a Banach algebra, that

[TABLE]

4.3.4 Computation of the bound $\hat{z}^{(m)}_{w}$ and $\delta_{w}$

From (48) and (45), on can verify that for each $j=1,\dots,n-1$,

[TABLE]

Using Lemma 12, we obtain that for $|\ell|<m$ and $j=1,\dots,n-1$,

[TABLE]

Moreover, for $j=1,\dots,n-1$,

[TABLE]

Combining (52), (53), (56) and (60), we define the uniform bound $\hat{z}^{(m)}$ which is then used to compute $\xi^{(m)}$ in (50). Moreover, combining (54), (57), (58), (59) and (61) provides the explicit bounds $\delta_{u}$, $\delta_{v}$ and $\delta_{w}$. All of these uniform bounds combined are finally used to compute the bound $Z_{1}$ in (51) which by construction satisfy (39).

4.4 $Z_{2}$ bound

Recall that we look for a bound $Z_{2}$ satisfying (40). Consider $Z_{2}$ satisfying

[TABLE]

Then, for any $b\in B_{r}(0)$, applying the Mean Value Inequality yields

[TABLE]

Given $\xi\in B_{r}({\bar{x}})$ and $h^{(1)},h^{(2)}\in B_{1}(0)$, we aim at bounding $\|D^{2}_{x}F(\xi,\omega)(h^{(1)},h^{(2)})\|_{X}$. Let

[TABLE]

which we denote by $z=(z_{\lambda},z_{\alpha},z_{u},z_{v},z_{w})=(0,z_{\alpha},0,z_{v},z_{w})$, where $z_{\lambda}$ and $z_{u}$ are both zero since $\eta$ and $f$ are linear. Denote

[TABLE]

Then, for $j=1,\dots,n-1$,

[TABLE]

Consider $r_{*}>0$ such that $r\leq r_{*}$. For $j=1,\dots,n-1$ and $i=1,2$,

[TABLE]

Then, for $j=1,\dots,n-1$,

[TABLE]

One verifies that

[TABLE]

and hence using the Banach algebra structure of $\ell_{\nu}^{1}$, we get that (for $p=1,2,3$)

[TABLE]

For $j=1,\dots,n-1$,

[TABLE]

and hence,

[TABLE]

Combining (62), (63) and (64), set

[TABLE]

and therefore, for all $b\in B_{r}(0)$,

[TABLE]

5 Results

In this section, we present several computer-assisted proofs of existence of spatial torus-knot choreographies. First fix the number of bodies $n$, a prescribed symmetry (7) (determined by the integer $k$), a resonance $(p,q)$, the frequency $\omega$ given in (4), and a Galerkin projection number $m$. Then compute a real numerical approximation ${\bar{x}}\in X$ of the finite dimensional projection $F^{(m)}$ defined in (35), where $X$ is defined in (43). Define the operators $A^{\dagger}$ and $A$ as in Section 3.2. Since the tail of the diagonal blocks of the approximate inverse $A$ (which is defined in (46)) involves the operator $D^{-1}$, we can easily show (using that $\ell_{\nu}^{1}$ is a Banach algebra under discrete convolutions) that the hypothesis (36) of Theorem 8 holds, that is $AF\colon X\times\mathbb{R}\to X$. Having described how to compute the bounds $Y_{0}$ in Section 4.1, $Z_{0}$ in (47), $Z_{1}$ in (51) and $Z_{2}$ in (65), we have all the ingredients to compute the radii polynomial defined in (41). The proof of existence then reduces to verify rigorously the hypothesis (42) of Theorem 8. This is done with a computer program in MATLAB implemented with the interval arithmetic package INTLAB, and available at [70]. All computations are performed with 16 decimal digits’ precision.

Let us present in details the computer-assisted proof resulting in the constructive existence of the torus-knot choreography of Figure 2.

Theorem 13**.**

Fix $n=5$ and consider the symmetry (7) with $k=3$. Let $(p,q)=(3,1)$ be the resonance. Let $s_{1}=\frac{1}{4}\sum_{j=1}^{4}\frac{1}{\sin(j\pi/5)}$ be given by (1) and the frequency $\omega=3\sqrt{s_{1}}$ be as in (4). Fix the Galerkin projection number $m=25$ and the decay rate parameter $\nu=1.03$. Consider the numerical approximation

[TABLE]

where the real and the imaginary part of the Fourier coefficients $({\bar{u}}_{j})_{\ell}$ can be found in the Appendix in Table 1. Then there exist sequences ${\tilde{u}}_{1},{\tilde{u}}_{2},{\tilde{u}}_{3}\in\ell_{\nu}^{1}$ such that

[TABLE]

with

[TABLE]

and such that $\mathcal{G}({\tilde{u}},\omega)=0$, with $\mathcal{G}$ defined in (8). Then $\{Q_{j}\}_{j=1}^{5}$ defined in the inertial frame by

[TABLE]

is a (renormalized) $6\pi$-periodic choreography that is symmetric by $2\pi/3$-rotations. Moreover, there exist countably many choreographies with frequencies near $\omega=3\sqrt{s_{1}}$.

Proof.

First denote by ${\bar{x}}=(\bar{\lambda},\bar{\alpha},{\bar{u}},{\bar{v}},{\bar{w}})\in\mathbb{C}^{2m(n+5)-3}=\mathbb{C}^{497}$ a numerical approximation of the finite dimensional reduction $F^{(497)}:\mathbb{C}^{497}\to\mathbb{C}^{497}$ defined in (35). The approximation satisfies ${\bar{x}}\in X$ and can be found in the file pt_five_bodies.mat available at [70]. Note that ${\bar{u}}\in\mathbb{C}^{3(2m-1)}=\mathbb{C}^{147}$ is recovered from the coefficients in Table 1 of the Appendix. Fix $\nu=1.03$. The MATLAB computer program proof_five_bodies.m available at [70] computes $Y_{0}$ as in Section 4.1, $Z_{0}$ in (47), $Z_{1}$ in (51) and $Z_{2}$ in (65), and verifies rigorously (using INTLAB) the hypothesis (42) of Theorem 8 with $r_{0}=4.7\times 10^{-10}$. Combining Theorem 8 and Proposition 10, there exists ${\tilde{x}}=(\tilde{\lambda},\tilde{\alpha},{\tilde{u}},\tilde{v},\tilde{w})\in X$ such that $F({\tilde{x}},\omega)=0$ and $\|{\tilde{x}}-{\bar{x}}\|_{X}\leq r_{0}=4.7\times 10^{-10}$. Hence, for a given $j\in\{1,2,3\}$,

[TABLE]

By construction of the Fourier map $F$ introduced in Section 2.4, the solution ${\tilde{x}}$ yields a $2\pi$-periodic solution $({\tilde{u}},\tilde{v},\tilde{w})$ of the delay equations (17), (18) and (19), which also satisfies the extra condition (20). By Proposition 6, ${\tilde{u}}$ satisfies $\mathcal{G}({\tilde{u}},\omega)=0$. The result follows from Proposition 3. The existence of countably many choreographies with frequencies near $\omega=3\sqrt{s_{1}}$ follows from Corollary 9 and the discussion thereafter. ∎

In the two left subfigures of Figure 4, we can visualize (in red) the $2\pi$-periodic solution ${\tilde{u}}$ satisfying the reduced delay equations (5). The initial condition ${\tilde{u}}(0)=(x_{0},y_{0},z_{0})$ of that red orbit can be found in Table 2 of the Appendix. This orbit is in the rotating frame. Still in the rotating frame, the position of the other bodies (in blue) can be recovered via the symmetry (7). In the two right subfigures of Figure 4, we can visualize the position of the bodies $Q_{1},\dots,Q_{5}$, which are now in the inertial frame. Since 3 and 5 are relative prime, the factor 3 in the equality $Q_{j}(t)=Q_{n}(t+3j\zeta)$ is just a re-ordering of the numbering of the bodies $j=1,2,3,4$.

Remark 14** (Resonance numbers versus torus winding).**

When $u_{n}(t)$ is a $p:q$ resonant orbit in the axial family with zero winding with respect to the $z$-axis, the choreography $Q_{n}(t)$ is a $(p,q)$-torus knot. See Corollary 4. In this case the resonance order $p:q$ in our functional analytic set-up corresponds exactly to the windings in the definition of a $(p,q)$-torus knot.

If on the other hand $u_{n}(t)$ has winding number one with respect to the $z$-axis – as in the case of the orbit $\tilde{u}(t)$ in Figure 4 (see the far left frame of that figure)– then the choreography $Q_{n}(t)$ (whose normalized period is $6\pi$) has toroidal winding in the second component one more than the $q$ value of the resonance. So even though the choreography illustrated in Figure 4 is resonant with order $p=3$ and $q=1$, the corresponding choreography is a $(3,2)$-torus knot after taking into account the non-trivial winding about the $z$-axis. We conclude that the choreography –illustrated in the center right and far right frames of Figure 4 –is a $(3,2)$-torus knot: that is, a trefoil knot. **

Following exactly the same approach as in Theorem 13, we prove the existence of several choreographies for $n=4$, $n=7$ and $n=9$ bodies. Results from several of our proofs are illustrated in Figures 3, 5, and 6 for four, seven, and nine bodies respectively. The computer-assisted proofs are obtained by running the codes proofs_four_bodies.m, proofs_seven_bodies.m and proofs_nine_bodies.m. The approximations can be found in the data files pts_four_bodies.mat, pts_seven_bodies.mat and pts_nine_bodies.mat. All files are available at [70]. In Table 2 of the Appendix, the initial conditions ${\tilde{u}}(0)=(x_{0},y_{0},z_{0})$ for the $n$-th body of each proven choreography is available. In Table 3 of the Appendix, some data for the proofs are given. For each of these proofs, the existence of countably many choreographies with near frequencies follows from Corollary 9 and the discussion thereafter. The reader interested in reproducing the choreographies via numerical integration will find at [70] the data files containing initial conditions – in inertial coordinates– for each of the $4,7$, and $9$ body choreographies illustrated in the Figures.

6 Acknowledgments

The authors wish to sincerely thank Sebius Doedel, who provided the numerical data for the $n$-body spatial choreographies computed in [21]. This data is massaged (via Fourier space Newton/continuation schemes we implemented in MatLab) to get appropriate approximate periodic solutions of the DDE, and is the starting point for all the analysis in the present manuscript. This material is based upon work supported by the National Science Foundation under Grant No. 1440140, while R.C. was in residence at the Mathematical Sciences Research Institute in Berkeley, California, during the Fall of 2018. R.C., J.-P.L., and J.D.M.J. were partially supported by a UNAM-PAPIIT project IA102818 and IN101020. C.G.A was partially supported by a UNAM-PAPIIT project IN115019. J.D.M.J was partially supported by NSF grant DMS-1813501. J.-P. L. was partially supported by NSERC.

Appendix

Tables 1, 2, and 3 in this appendix contain numerical data needed in the proofs discussed in the main body of the present work.