1. Introduction

Here we provide an overview of what is known and what is not known, about an interesting dynamical system known as the Kepler-Heisenberg problem. The main idea is to pose a version of the classical Kepler problem of planetary motion, but in a sub-Riemannian setting. The result is system which is surprisingly rich and beautiful, mysterious in some ways but tame in others, offering a substantial number of questions which seem non-trivial yet tractable.

1.1. History and motivation. The curved Kepler problem, typically set in a Riemannian manifold of constant curvature, is an active area with a fascinating history. In particular, while the co-discovery of hyperbolic space by Lobachevsky and Bolyai is famous mathematical lore, both men also independently attempted to pose the Kepler problem in this new geometry around 1835. Dirichlet and Schering worked on this problem shortly after. In 1873, Lipschitz posed and partially solved the Kepler problem on the 3-sphere, as did Killing (with a different notion of gravity) in 1885. In 1902, Liebmann proved that orbits are conics in both hyperbolic and and spherical 3-space. There is a parallel history of studying the Kepler problem on surfaces as well, and today there is active work on \(n\)-body problems in more general curved spaces. Details can be found in [1] and the references therein (while this paper is published, the arXiv version contains a much more elaborate history).

The original motivation for studying the Kepler-Heisenberg problem was to extend this idea to sub-Riemannian manifolds, which generalize Riemannian ones. In particular, the three dimensional Heisenberg group is the simplest non-trivial sub-Riemannian space and the natural place to begin such an investigation. While no knowledge of sub-Riemannian geometry is needed for this article, the interested reader can see [9] for an introduction to this subject in general, and details of the Heisenberg group in particular. As seen below, the resulting Kepler-Heisenberg system retains certain aspects of the classical Kepler problem in Euclidean space, while reflecting the unique geometry of the Heisenberg group. We are thus left with obvious questions about generalizing both the setting to other sub-Riemannian spaces and the force to other classical mechanical systems. See Section 3.

We can provide an alternative motivation as well. The classical Kepler problem enjoys the famous three Kepler's Laws of Planetary Motion. Do these hold in curved spaces? For spherical and hyperbolic 3-space, the first two laws hold but the third fails. Why? Kepler's third law is a consequence of homogeneity of the gravitational potential with respect to Euclidean dilations. But spherical and hyperbolic spaces do not admit such dilations. Moreover, by a result of Gromov (see [7]), there are no homogeneous Riemannian manifolds which admit dilations besides Euclidean spaces. Thus, if one wants to pose a curved Kepler problem in which Kepler's third law holds, one must leave the Riemannian realm. Certain important sub-Riemannian spaces known as Carnot groups do admit dilations. The Heisenberg group is the simplest, and we will see below that a version of Kepler's third law does indeed hold there.

1.2. The problem. We treat the Heisenberg group \(\mathcal H\) as the smooth manifold \(\mathbb R^3\) equipped with a certain non-Euclidean metric. Take the usual global coordinates \((x, y, z)\) and define two vector fields

\( X=\frac{\partial}{\partial x}-\frac{1}{2}y\frac{\partial}{\partial z}\quad \text{and}\ \quad Y=\frac{\partial}{\partial y}+\frac{1}{2}x\frac{\partial}{\partial z}. \)

These span a two-dimensional bracket-generating distribution \(\mathcal D\). The Kepler-Heisenberg problem is roughly the classical Kepler problem in \(\mathbb R^3\) with the constraint that all orbits are tangent to \(\mathcal D\); such curves are called horizontal. Declaring that \(X\) and \(Y\) are orthonormal defines an inner product on \(\mathcal D\). Then \(\mathcal H\) is now equipped with a sub-Riemannian structure whereby we can measure lengths of horizontal curves and (by the Chow-Rashevskii theorem) distances between points. The resulting geometry enjoys three interesting properties:

• the length of a horizontal curve is equal to the Euclidean length of its projection to the \(xy\)-plane (\(ds^2=(dx^2 +dy^2)|_{\mathcal D}\)),
• the \(z\)-coordinate of a horizontal curve grows like the area traced out by its projection (Stokes' theorem),
• geodesics are helices projecting to circles or lines (dual isoperimetric problem).

Our configuration space is \(\mathcal H\) equipped with this structure. Here we pose the Kepler-Heisenberg problem as a Hamiltonian system, taking \(T^*\mathcal H\) with coordinates \((x, y, z, p_x, p_y, p_z)\) as phase space. Consider the dual momenta to the vector fields \(X\) and \(Y\) given by

\( P_X=p_x-\tfrac{1}{2}yp_z \quad \text{and}\ \quad P_Y=p_y+\tfrac{1}{2}xp_z. \)

Then, as in the Riemannian case, we take our kinetic energy to be \(K=\tfrac{1}{2}(P_X^2+P_Y^2)\). This is induced by our metric and generates geodesic flow. Choosing the gravitational potential is more subtle. In the classical case, the potential \(1/r\) has many important properties. It is most obviously the reciprocal of distance from the origin, but for most (sub-)Riemannian manifolds there is no explicit expression for such distance, so this characterization is inadequate. Instead, we recognize \(1/r\) as (a multiple of) the fundamental solution to the Laplacian on \(\mathbb R^3\), and choose our potential accordingly. Thus we choose the Kepler-Heisenberg potential as \(U=-\frac{1}{8\pi}\Big((x^2+y^2)^2+16z^2\Big)^{-1/2}\), which is the fundamental solution to the Heisenberg sub-Laplacian \(\Delta_{\mathcal{H}}=X^2 +Y^2\) according to [6]. Note the singularity at the origin, which represents our planet colliding with the sun.

Thus we obtain our Hamiltonian

\( H=\underbrace{\tfrac{1}{2}((p_x-\tfrac{1}{2}yp_z)^2+(p_y+\tfrac{1}{2}xp_z)^2)}_K\underbrace{-\frac{1}{8\pi\sqrt{(x^2+y^2)^2+{16}z^2}}}_{U}. \)

The Kepler-Heisenberg problem is to analyze and eventually solve Hamilton's equations for \(H\).

2. Results

We first give some basic properties of the system, most of which can be found in [10]. An important tool will be the Heisenberg dilation \((x, y, z)\mapsto (\lambda x, \lambda y, \lambda^2 z)\) for \(\lambda>0\).

Proposition 1. The Kepler-Heisenberg system has two integrals of motion: the total energy \(H\) and the angular momentum \(p_{\theta}=xp_y-yp_x\). The dilational momentum \(J=xp_x+yp_y+2zp_z\) satisfies \(\dot J = 2H\). The dynamics are therefore integrable on the codimension one submanifold \(\{H=0\}\), where \(J\) is also conserved. Periodic orbits must satisfy \(H=J=0\).

These are the known symmetries -- there may be others. So we do not know whether this system is integrable or not. But we consider it `at least mostly integrable'.

Ostensibly, we should be able to integrate the zero-energy dynamics. However, the dilational action is non-compact, so the standard machinery of action-angle variables does not immediately apply. (While there is a generalized theory, as in [5], it is non-constructive.) Our invariant submanifolds can be visualized as cones of 2-tori, \(\mathbb T^2 \times \mathbb R^+\), with the cone point at \(0 \in \mathbb R^+\) representing the singularity at sun. If an orbit's constant dilational momentum is positive (resp. negative) then it winds away (resp. towards) the cone point; if \(J=0\) then the orbit stays on a compact \(\mathbb T^2\) and the motion is quasi-periodic or periodic. In [10], the integration of the \(H=0\) system is reduced to parametrizing a family of degree 6 algebraic curves, but there is no general method for obtaining such a parametrization.

As mentioned in the Introduction, we also have a version of Kepler's third law:

\( T^2 = ka^4 \)

where \(T\) represents an orbit's period, \(a\) is its size, and \(k\) is a universal constant. This is a consequence of our Hamiltonian being homogeneous of degree \(-2\) with respect to the Heisenberg dilation. Note that we do have Kepler's second law as well, since it is equivalent to the conservation of angular momentum. Kepler's first law, however, fails in the Heisenberg group (recall that it held in spherical and hyperbolic three-space, where the third law failed). We do have periodic orbits (see below), but they are not conic sections.

Finally, we mention a few remaining small but interesting results. All negative energy solutions are bounded. The only solutions contained in the \(xy\)-plane are lines through the origin. On the invariant submanifold \(\{z=p_z=p_{\theta}=0\}\) we have a classical central force problem in the plane, and the radial distance \(r(t)\) traces out a hyperbola, ellipse, or parabola if the energy is positive, negative, or zero, respectively. Lastly, there is one very strange family of orbits: each sits at a constant position on the \(z\)-axis while its momenta \(p_z\) grows linearly in time.

2.1. Periodic orbits. The existence of periodic orbits was established in [11]. However, the proof followed the direct method in the calculus of variations and provided very little information about such orbits beyond their existence. Numerical methods in [3] allowed the discovery of a very rich and beautiful symmetry structure.

Theorem 2. For any rational \(j/k \in (0,1]\), there is a periodic Kepler-Heisenberg orbit with rational rotation number \(j/k\).

Rational rotation numbers are treated in [8, 12]; here it simply means that after one \(k\)th of the period, the planet has rotated \(2\pi \cdot \frac{j}{k}\) radians counterclockwise about the \(z\)-axis. See Figures 1 and 2. All such orbits possess \(k\)-fold rotational symmetry about the \(z\)-axis, except the orbit with rotation number 1/1 (its projection does, but not \(z(t)\)).

Figure 1. Periodic orbits with rotation numbers 5/9, 6/41, 1/4, and 7/8. In all figures we show the projection of the orbit to the \(xy\)-plane.; the \(z\)-coordinate grows like the area traced out.

Figure 2. All periodic orbits with \(k=4\) or \(5\).

This symmetry structure was uncovered via numerical explorations. We strongly encourage other researchers to use our codebase [2]. All code is hosted publicly as free, open-source software; explicit instructions to reproduce each plot are provided there, with all relevant documentation and licensing. All results in [3] are thus 100% available and exactly reproducible.

Our Python program searches for a closed orbit nearby some given initial condition, which may be input directly or randomly generated. First, the program checks whether this starting initial condition \(X_0\) is close enough to a closed orbit. To do so we use a symplectic integrator and shooting method with objective function obj giving the smallest local minimum of distance from \(X_0\). If \(X_0\) is satisfactory, we use a Monte Carlo method to optimize obj as follows. Randomly pick a point in radially symmetric distribution on an annular neighborhood centered at \(X_0\), then symplectically integrate and evaluate obj. If obj does not decrease, try a new point; if obj decreases, use this point as new \(X_0\). Repeat this update step 1000 times. The output is an image, a binary data file, and a human readable summary containing all pertinent information about the orbit (initial conditions, angular momentum, period, etc.). The program also detects the rotation number using Fourier analysis.

Finally, our numerical investigations discovered another somewhat surprising phenomenon that we do not fully understand. Due to the symmetries and constraints of the system, we can reduce the search for periodic orbits to one dimension, which we can parametrize by angular momentum \(p_{\theta}\). Then as a function of \(p_{\theta}\), the rotation numbers are arranged according to the Farey sequence. See Table 1.

Table 1. As \(p_{\theta}\) grows, the rotation numbers are distributed like the Farey sequence.

2.2. Self-similarity. Zero-energy orbits are relatively well behaved and understood. Almost all exhibit an attractive self-similarity property (see [4]). To the best of our knowledge these solutions represent the first occurrence of this type of self-similarity within solutions to ODEs of Hamiltonian type.

Theorem 3. Suppose \(c\) is a curve in \(T^*\mathcal{H}\) such that:

• \(c\) solves Hamilton's equations
• \(H(c(t))=0\)
• \(z\) has consecutive zeros at \(t_0, t_1, t_2\)

Then \(c\) is entirely determined by its value on the fundamental domain \([t_0, t_2]\).

See Figure 4, which is explained in detail below. The self-similarity here involves a rotation, a dilation, and a time reparametrization. The first two assumptions in the theorem are clear, but the third requires some explanation. Numerical evidence suggests the following.

Conjecture 4. All zero-energy orbits avoiding the \(z\)-axis have \(z(t)\) oscillatory.

By oscillatory we mean that \(z(t)\) has infinitely many zeros without being identically zero on any interval, which of course implies the existence of three consecutive zeros needed for our theorem. We tried and failed to prove this conjecture using Sturm comparison methods. It is necessary that we avoid the \(z\)-axis in this conjecture; we have counterexamples otherwise. But the \(z\)-axis is special in the Heisenberg group -- it is the cut and conjugate locus for the origin -- and needs to be treated with care. However, numerical evidence suggests that even those zero-energy orbits which do intersect the \(z\)-axis are self-similar, although perhaps in a different sense than those which do not intersect the \(z\)-axis, with \(z(t)\) not necessarily oscillatory. We provide details at the end of this section. The upshot is the following.

Conjecture 5. All zero-energy orbits are self-similar.

The self-similarity in Theorem 3 is most clearly understood in new coordinates. Let

\[\begin{align*} s &= \frac{1}{4}\log ((x^2+y^2)^2+16z^2) \\ \theta &= \arg(x,y) \\ u &= \arg(x^2+y^2, 4z). \end{align*}\]

Then the corresponding momenta \(p_s\) and \(p_{\theta}\) correspond precisely to the dilational and angular momenta; that is, \(p_s=J\). This allows us to compute the dilational and rotational factors \(\lambda\) and \(\varphi\) for any given orbit satisfying our hypotheses as

\( \lambda = \exp (s(t_2) - s(t_0)) \qquad \text{and} \qquad \varphi = \theta(t_2) - \theta(t_0). \)

The rotational part is standard, but finding the correct coordinates for the dilational part requires solving a system of non-linear PDEs. Armed with the dilational factor \(\lambda\) we can now explicitly determine the time reparametrization \(\tau\) as follows (if \(\lambda=1\) then no reparametrization is necessary):

\[\begin{align*} \xi(t) &=\frac{1}{2}\log_{\lambda}\left(1-(t-t_0)\frac{1-\lambda^2}{t_2-t_0} \right) \\ \tau(t) &= t_0+(t_2-t_0)\frac{1-\lambda^{2\xi(t) + 2\lfloor \xi(t) \rfloor}}{1-\lambda^2}. \end{align*}\]

Suppose we are given a segment of a curve \(c\) on a fundamental time domain \([t_0, t_2]\), and we want to know \(c(t)\) at some future time \(t\). Then \(c\) can be extended from \([t_0, t_2]\) to next fundamental domain \([t_2, t_2 + \lambda^2(t_2-t_0)]\) by appropriately rotating, dilating, and reparametrizating time. We can then iterate to move forward to any future time (until possible collision), and can easily modify to move backward to past time. The function \(\xi\) plays an important discretizing role: its floor gives the number of iterations needed. That is, \(\lfloor \xi(t) \rfloor\) is the number of rotated and dilated copies of \(c\) which need to be pasted together in order to reach \(c(t)\). Note that consecutive "periods" scale by \(\lambda^2\), giving a geometric series, so we have the following corollary.

Conjecture 6. Orbits satisfying the hypotheses of Theorem 3 stratify according to Table 2. All collisions happen in finite time with

\( t_{\text{col}}=t_0+\frac{t_2-t_0}{1-\lambda^2}. \)

Table 2. Summary of behavior.

Figure 3 provides sample orbits from the stratification in Table 2.

Figure 3. Sample orbits with negative, positive, zero dilational momentum.

Most of the self-similarity discussed in this section is apparent in Figure 4, which displays an orbit with negative dilational momentum. On the left we see the projection of the orbit to the \(xy\)-plane, on the right we see \(z\) as a function of time. In both images the blue curve represents the orbit itself (numerically integrated), the purple segment represents the part over a chosen fundamental domain, and the orange segment is a rotated and dilated version of the purple. The idea is that we are given the purple segment, from which we can recover the orange segment, which in turn matches the blue curve. This can be iterated forward or backward in time. The vertical lines on the right show zeros of \(z(t)\); the given fundamental domain \([t_0, t_2]\) is shown with bold lines. We can see how these "periods" shrink geometrically and a collision occurs in finite future time. The codebase used to create Figure 4 will also be available at [2].

Figure 4. A self-similar orbit.

To conclude this section we revisit Conjecture 5 and the discussion preceding it. Recall that the \(z\)-axis is the conjugate and cut locus for the origin in the Heisenberg group. Zero-energy orbits meeting the \(z\)-axis need not have \(z(t)\) oscillatory, but we believe they are all still self-similar. However, numerical investigations suggest that such orbits bifurcate into two families with qualitatively different behavior, and this bifurcation occurs when dilational momentum \(|J|=\frac{1}{2\sqrt{\pi}}\).

Conjecture 7. Consider a zero-energy orbit meeting the \(z\)-axis. If \(|J|<\frac{1}{2\sqrt{\pi}}\) then we have self-similarity like a generic \(H=0\) orbit. In particular, \(z(t)\) oscillates and the orbit's projection looks like a dilating figure eight. If \(|J|>\frac{1}{2\sqrt{\pi}}\) then we have self-similarity like Heisenberg geodesic. In particular, \(z(t)\) is monotonic and the orbit looks like a dilating helix. See Figures 5 and 6.

Figure 5. Starting on \(z\)-axis with small dilational momentum.

Figure 6. Starting on \(z\)-axis with large dilational momentum.

3. Open Questions

Here we simply list some of the many questions we have about this system and its generalizations. Some of these we have investigated without success, others we have barely considered.

1. Is the Kepler-Heisenberg system integrable? Are there any other non-trivial symmetries?
2. How can we explicitly integrate the zero-energy system?
3. How do orbits with nonzero energy behave? Nearly all our efforts have focused on the \(H=0\) case, where the dynamics are integrable and where the periodic orbits live.
4. Are Conjectures 4, 5, 7 true?
5. Can we regularize collisions with the sun?
6. Can we quantize this system in a meaningful way?
7. The richness of this system gives us hope that other "sub-Riemannian mechanics" problems may be worth studying. In each case we would take the sub-Riemannian Hamiltonian as our kinetic energy and choose some appropriate potential energy. But this gives freedom to choose both the geometry and the potential:

• Can we study the Kepler problem on higher dimensional Heisenberg groups? Other Carnot groups? Other sub-Riemannian spaces? Is there a non-Euclidean setting where all three Kepler laws hold?
• Can we study the two body problem (not same as Kepler!) on the Heisenberg group? Three bodies? Other classical mechanics problems? Other potentials?
• Combine the previous two items. The \(n\)-body problem in jet space?

References

[1]   F. Diacu, E. Perez-Chavela, and M. Santoprete, The \(n\)-body problem in spaces of constant curvature, https://arxiv.org/abs/0807.1747, 2008.

[2]   V. Dods, Kepler-Heisenberg Problem Computational Tool Suite, available at: https://github.com/vdods/heisenberg, 2018.

[3]   V. Dods and C. Shanbrom, Numerical methods and closed orbits in the Kepler-Heisenberg problem, Exp. Math., 28 (2019), pp. 420-427.

[4]   V. Dods and C. Shanbrom, Self-similarity in the Kepler-Heisenberg problem, to appear in J. Nonlinear Sci., https://arxiv.org/abs/1912.12375.

[5]   E. Fiorani, G. Giachetta, and G. Sardanashvily, The Liouville-Arnold-Nekhoroshev theorem for non-compact invariant manifolds, J. Phys. A, 36 (2003), pp. L101-L107.

[6]   G. Folland, A fundamental solution to a subelliptic operator, Bull. Amer. Math. Soc., 79 (1973), pp. 373-376.

[7]   M. Gromov, Metric structures for Riemannian and Non-Riemannian Spaces, Birkhauser, 2007.

[8]   A. Katok and B. Hasselblatt, Introduction to the Modern Theory of Dynamical Systems, Cambridge University Press, 1995.

[9]   R. Montgomery, A Tour of Subriemannian Geometries, Their Geodesics and Applications, AMS Mathematical Surveys and Monographs, 91 (2002).

[10]   R. Montgomery and C. Shanbrom, Keplerian motion on the Heisenberg group and elsewhere, Geometry, mechanics, and dynamics: The Legacy of Jerry Marsden, Fields Inst. Comm., 73 (2015), pp. 319-342.

[11]   C. Shanbrom, Periodic orbits in the Kepler-Heisenberg problem, J. Geom. Mech., 6 (2014), pp. 261-278.

[12]   S. Tabachnikov, Geometry and Billiards, AMS Student Mathematical Library, 30 (2005).