One-Dimensional Relativistic Self-Gravitating Systems

1. Introduction

One of the oldest problems in physics is the N-body problem: determination of the motion of a system of N particles mutually interacting through specified forces. This problem appears in a broad variety of subfields of physics, including cosmology, stellar dynamics, planetary motion, atomic physics, and nuclear physics. The N-body problem is a particular challenge if the interactions are purely gravitational. Although an exact solution is known for the 2-body problem in pure Newtonian gravity in three spatial dimensions, there is no closed form solution for large N, even for $N=3$ [1], though particular solutions exist in restricted cases [2]. No exact solution is known in the general-relativistic case even for $N=2$, since it experiences dissipation of energy in the form of gravitational radiation.

One-dimensional self-gravitating systems (OGS) have played an important role in advancing our understanding of the gravitational N-body problem [3]. Such systems have been of interest for over half a century, where they have played an important role in astrophysics and cosmology for more than 30 years [4]. Apart from being prototypes for studying the behaviour of gravity in higher dimensions, they also approximate the behaviour in 3 spatial dimensions of some physical systems. Examples include very long-lived core-halo configurations that model a dense massive core in near-equilibrium, surrounded by a halo of high kinetic energy stars that feebly interact with the core [4,5,6]. Other examples include cosmological models [7,8] the dynamics of stars in a direction orthogonal to the plane of a highly flattened galaxy [9], shells of matter interacting with a spherical globular cluster [10], and collisions of flat parallel domain walls moving in directions orthogonal to their surfaces.

Although the connection between the idealized OGS and natural astrophysical systems can be tenuous, the accuracy and ease with which their dynamical evolution may be simulated has remained the principal motivation for continued study of the OGS. Unlike 3-dimensional self-gravitating systems, in which the motion of the (point) masses must be numerically integrated, the OGS admits direct computation of the particle (or sheet, or shell) crossings. This provides accurate computation of the evolution of the system over many dynamical time scales. Furthermore, a number of interesting questions concerning the statistical properties of the OGS remain open, including whether or not it can attain a state of true equilibrium from arbitrary initial conditions [5], its ergodic behaviour, the circumstances (if any) under which equipartition of energy can be attained [11], and the appearance of fractal behaviour [3,12].

For three decades, studies of the OGS have been in a non-relativistic context, assuming Newtonian gravity with its standard causal structure [3,7,8,11,13,14,15,16,17,18,19,20,21]. Research into relativistic one-dimensional self-gravitating systems (ROGS) was generally ignored. In large part this was because relativistic effects do not play a dominant role in stellar dynamics, but it was also due to the lack of a theoretical framework for relativistic gravity in one spatial dimension. The Einstein tensor is identically zero in this $(1+1)$-dimensional spacetime context, and so Einstein’s equation at face value would simply imply vanishing stress-energy. However reduction of the number of spatial dimensions in a relativistic context can be expected to be quite useful since gravitational radiation (at least due to spin-2 gravitons) cannot exist in less than 3 spatial dimensions. However most (if not all) of the remaining conceptual features of relativistic gravity are retained in lower dimensions, and so one might hope to obtain insight into the nature of relativistic classical and quantum gravitation in a wide variety of physical situations by studying the ROGS.

It is straightforward to find a set of equations governing the motion of particles – these are furnished by the geodesic equations. In addition to this, what is needed to study ROGS is a set of equations governing the dynamics of the spacetime metric in a self-consistent way. Early versions of $(1+1)$–dimensional gravity [22,23] set the Ricci curvature scalar equal to a constant, yielding trivial dynamics for the spacetime metric (although containing sufficiently interesting features [24] such that this theory is still of interest today [25] ). Intensive investigation of a wide variety of gravitational theories ensued a few years later, primarily motivated by a quest to understand quantum gravity in a simplified context [26]. The overwhelming majority of such investigations were concerned with the (quantum) dynamics of the spacetime metric, and not with the dynamical motion of particles in such spacetimes.

At about the same time that interest $(1+1)$-dimensional quantum gravity began, investigations of the ROGS also began. The purpose of this article is to review the origins, results, and status of relativistic one-dimensional self-gravitating systems. After a brief review of the OGS, I begin by reviewing how the $D\to 2$ limit of D-dimensional general relativity [27] can be self-consistently coupled to point particles, thereby yielding the ROGS. The equations of motion for the particles are obtained using the canonical formalism, which I describe in some detail. I shall then consider in turn the 2-body, 3-body, 4-body, and N-body ROGS, discussing their distinctions from the OGS, their salient features, their chaotic behaviour, and their statistical properties, as relevant. I conclude by discussing a number of interesting open problems for relativistic one-dimensional self-gravitating systems. While other constants will retain their values throughout, the speed of light c will generally be set to unity, and only explicitly written where relevant for instructive purposes.

2. Non-Relativistic Self-Gravitating Systems

For a system of particles the Hamiltonian in Newtonian gravity in two dimensions is

$$H=\sum _a\frac{p_a^2}{2m_a}+\pi G\sum _{a,b=1}^N{m}_a{m}_b|z_a-z_b|$$

(1)

where $m_a$, $z_a$ and $p_a$ are the mass, the coordinate location, and the momentum of a-th particle, respectively, and G is the gravitational constant. The potential between any two particles is proportional to the product of their masses and the spatial separation between them, as expected from dimensionally continuing the well-known potential

$$V=-G_d\sum _{a\ne b}^N\frac{m_a{m}_b}{|z_a-z_b{|}^{d-2}}$$

(2)

of Newtonian gravity in d spatial dimensions, where $G_1=\pi G$. When $d=1$ the potential in (1) vanishes, and so the restriction $a\ne b$ in (1) is not required.

The equations of motion of the OGS (1) are given by Hamilton’s equations

$$\begin{array}{ccc}\hfill {\dot{z}}_a& =& \frac{\partial H}{\partial p_a}=\frac{p_a}{m_a}\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill {\dot{p}}_a& =& -\frac{\partial H}{\partial z_a}=-2\pi G\sum _b{m}_a{m}_b\frac{\partial |z_a-z_b|}{\partial z_a}\hfill \\ & =& -2\pi G\sum _b{m}_a{m}_b\mathrm{sgn}(z_a-z_b)\hfill \end{array}$$

(4)

yielding

$${\ddot{z}}_a=-2\pi G\sum _b{m}_b\mathrm{sgn}(z_a-z_b)$$

(5)

for the acceleration of the a-th particle.

We see that each particle experiences a constant force from each of the other particles, where the sign (- or +) of the force from any given particle depends on whether $z_a>z_b$ or vice-versa. The force is therefore always attractive: if $z_a$ is to the right of $z_b$ ($z_a>z_b$) then particle a will accelerate leftward toward $z_b$ (and b rightward toward a) until they meet, after which time $z_a$ is to the left of $z_b$ ($z_a<z_b$) and the force changes sign, accelerating particle a rightward. This scenario assumes that the particles can pass through each other without any influence, as would be appropriate for parallel sheets of particulate matter where collisions between the particles can be neglected. It is of course possible to include additional structure – for example, modelling the particles as impenetrable points would make them bounce off of each other – but this would detract from the study of pure gravitational effects. Such additional structure will not be considered in this article, apart from the inclusion of attractive and repulsive electromagnetic interactions.

Consider the case $N=2$. If the particles are initially separated by some distance d, they will move toward each other with constant acceleration until they cross, after which the acceleration of each flips sign. The particles fly apart increasingly slowly until they reach a maximal separation d, after which they move toward each other again with increasing speed. After crossing a second time, the particles separate, moving with decreasing speed until they return to their original positions. Assuming no other interactions, the motion then repeats perpetually. Prior to crossing, the entire scenario is equivalent to that of a body of mass m falling near the surface of the Earth.

For N particles, every particle initially undergoes constant acceleration until the first two particles cross. This causes a sign flip in the force between each particle in the crossing pair, thereby changing the magnitudes (and perhaps signs) of the accelerations of each due to the presence of the other particles. As more particles cross, this changes the accelerations of more and more particles, in general yielding chaotic motion.

The simplest example of this occurs for $N=3$. This 3-body OGS has been shown to exhibit a mild form of chaos [3]. This OGS can be mapped to a single-particle moving in 2 spatial dimensions in a hexagonal-well potential $V(x,y)$, whose sides are impenetrable flat sheets. If the masses of all particles are equal then the shape is that of a regular hexagon; unequal masses distort this symmetry so that the sides are of unequal length [28,29]. Alternatively, the 3-body OGS can be regarded as a single-particle under the influence of a constant gravitational field in 2 spatial dimensions that bounces off of a symmetric wedge of angle $2\theta$, where this angle parametrizes the relative inequality of the masses [3].

Constructing a relativistic of the Hamiltonian (1) is somewhat subtle. This is the subject of the next section.

3. Relativistic Gravity Coupled to Point Particles

In three spatial dimensions a self-gravitating system would consist of a set of N particles minimally coupled to Einstein gravity. Its action in n-dimensions is [30]

$$\begin{array}{c}\hfill I_n=\int d^{n}x\left[\frac{1}{2{\kappa }_n}\sqrt{-g}(R-2{\Lambda }_n)\right.-\left.\sum _{a=1}^N{m}_a\int d{\tau }_a{\left\{-g_{\mu \nu }\left(x\right)\frac{d{z}_a^{\mu }}{d{\tau }_a}\frac{d{z}_a^{\nu }}{d{\tau }_a}\right\}}^{1/2}{\delta }^{\left(n\right)}(x-z_a\left({\tau }_a\right))\right]\end{array}$$

(6)

where $|{\Lambda }_n|=\frac{(n-2)(n-1)}{2{\ell }^2}$ is the cosmological constant (whose sign $+/-$ respectively corresponds to asymptotically de Sitter/Anti de Sitter spacetime), ${\kappa }_n=8\pi G_n/c^4$ is the gravitational coupling and R is the Ricci scalar. Systems of astrophysical interest have $n=4$. Notwithstanding issues connected with collisions between the particles, the field equations that follow from this action embody what we expect from a self-gravitating system: the curvature of spacetime governs how the particles move along the trajectories $z_a\left({\tau }_a\right)$, and the masses and motion of these particles in turn govern how spacetime dynamically curves.

A ROGS that resembles a relativistic three-dimensional self-gravitating system as closely as possible should therefore have the following features.

• The stress-energy of the particles generates spacetime curvature in as simple a manner as possible.
• The curvature of spacetime guides the motion of each particle in accord with the equivalence principle, in the absence of any extraneous forces.
• The dynamics of the system is self-consistent.

One might expect to obtain the ROGS by simply setting $n=2$ in (6) (1 space dimension, 1 time dimension), but this will not work: the Ricci scalar is a total derivative in $(1+1)$ dimensions, in turn implying that the Einstein tensor is identically zero for any metric.

However there is a way of taking the $n\to 2$ limit of the action (6) [27]. Consider the action

$$\begin{array}{ccc}\hfill I_n^{EH}& =& \frac{1}{2{\kappa }_n}\left(\int d^{n}x\sqrt{-\tilde{g}}\tilde{R}-\int d^{n}x\sqrt{-g}R\right)\hfill \\ & =& \frac{1}{2{\kappa }_n}\int d^{n}x\left[e^{ϵ\Psi /2}\left[\left(R-(ϵ+1)\square \Psi \right)-\frac{1}{4}ϵ(ϵ+1){(\partial \Psi )}^2\right]-R\right]\hfill \end{array}$$

(7)

where $\tilde{g}=e^{\Psi }g$ is a conformally rescaled metric and $ϵ=n-2$. Expanding in powers of $ϵ$ yields

$$I_2^{EH}=\frac{n-2}{4{\kappa }_n}\int d^{2}x\sqrt{-g}\left[\Psi R-\Psi \square \Psi -\frac{1}{2}{(\nabla \Psi )}^2\right]$$

(8)

and so by setting ${\kappa }_2=\frac{2{\kappa }_n}{n-2}$ we obtain

$$I_2^{EH}\equiv \underset{n\to 2}{lim}{I}_n^{EH}=\frac{1}{2{\kappa }_2}\int d^{2}x\sqrt{-g}\left[\Psi R+\frac{1}{2}{(\nabla \Psi )}^2\right]$$

(9)

upon discarding total divergences. The limit $n\to 2$ in (6) is straightforward for the other two terms, yielding

$$\begin{array}{ccc}\hfill I& =& \int d^{2}x\sqrt{-g}\left[\frac{1}{2\kappa }\left(\Psi R+\frac{1}{2}{\left(\nabla \Psi \right)}^2-2\Lambda \right)\right.\hfill \\ & & \left.-\sum _{a=1}^N{m}_a\int d{\tau }_a{\left\{-g_{\mu \nu }\left(x\right)\frac{d{z}_a^{\mu }}{d{\tau }_a}\frac{d{z}_a^{\nu }}{d{\tau }_a}\right\}}^{1/2}{\delta }^{\left(2\right)}\left(x-z_a\left({\tau }_a\right)\right)\right]\hfill \end{array}$$

(10)

as the action for the ROGS, where $\kappa ={\kappa }_2$ and $\Lambda ={\Lambda }_2=\pm \frac{1}{{\ell }^2}$.

Note that this procedure incorporates an an additional field $\Psi$ in the gravitational action. This might seem to be in tension with the second desirable feature in the list above. However the field equations derived from the variations $\delta \Psi$ and $\delta g_{\mu \nu }$ are

$$R-g^{\mu \nu }{\nabla }_{\mu }{\nabla }_{\nu }\Psi =R-\frac{1}{\sqrt{-g}}{\partial }_{\mu }(\sqrt{-g}{g}^{\mu \nu }{\partial }_{\nu }\Psi )=0$$

(11)

$$\frac{1}{2}{\nabla }_{\mu }\Psi {\nabla }_{\nu }\Psi -\frac{1}{4}{g}_{\mu \nu }{\nabla }^{\lambda }\Psi {\nabla }_{\lambda }\Psi +g_{\mu \nu }{\nabla }^{\lambda }{\nabla }_{\lambda }\Psi -{\nabla }_{\mu }{\nabla }_{\nu }\Psi =\kappa T_{\mu \nu }-g_{\mu \nu }\Lambda$$

(12)

where

$$\begin{array}{ccc}\hfill T_{\mu \nu }& =& -\frac{2}{\sqrt{-g}}\frac{\delta {\mathcal{L}}_{\mathcal{M}}}{\delta g^{\mu \nu }}=\sum _a{m}_a\int d{\tau }_a\frac{1}{\sqrt{-g}}{g}_{\mu \sigma }{g}_{\nu \rho }\frac{d{z}_a^{\sigma }}{d{\tau }_a}\frac{d{z}_a^{\rho }}{d{\tau }_a}{\delta }^{\left(2\right)}(x-z_a\left({\tau }_a\right))\hfill \end{array}$$

(13)

where the last term in (10) is the matter Lagrangian ${\mathcal{L}}_{\mathcal{M}}$ . Taking the trace of (12) and inserting the result into (11) gives

$$R+2\Lambda =\kappa T_{\mu }^{\mu }.$$

(14)

which shows that the evolution of the metric depends only on the stress-energy, and decouples from the evolution of $\Psi$. The variation $\delta z_a^{\mu }$ yields

$$\frac{d}{d{\tau }_a}\left\{g_{\mu \nu }\left(z_a\right)\frac{d{z}_a^{\nu }}{d{\tau }_a}\right\}-\frac{1}{2}{g}_{\nu \lambda ,\mu }\left(z_a\right)\frac{d{z}_a^{\nu }}{d{\tau }_a}\frac{d{z}_a^{\lambda }}{d{\tau }_a}=0$$

(15)

which is the geodesic equation, or rather equations since there is one per particle.

The system (14,15) forms a closed relativistic self-gravitating system of N point particles. Spacetime curvature is determined from the stress-energy of the particulate matter from (14); in turn the evolution of the particles is determined by the spacetime curvature via (15).

The theory given by the action (10) is known as $R=T$ theory [31,32,33] . Its classical properties, including gravitational collapse, black holes, cosmological solutions, solitonic properties, and thermodynamics have been extensively studied [34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68]. Its chief interest lies in the fact that it captures the essence of classical general relativity in two spacetime dimensions, and has $(1+1)$–dimensional analogs of many of its properties [34,69,70]. Moreover, it has a well-defined Newtonian limit [32,34]. This is in contrast to generic scalar-tensor theories [71], where the dilaton does not decouple from the evolution of the gravitational field.

The quantum properties of $R=T$ theory have also received attention from a variety of perspectives [36,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88]. There is also a supersymmetric version of the theory [89], which has supersymmetric black hole solutions [90,91,92]. Recently experiments have been carried out [93,94] that test certain aspects of the theory, albeit in a simulated context.

The remainder of this article is concerned with the N-body problem and its solutions as determined from the equations that follow from the action (10). The canonical approach is the most useful way to obtain these equations. This is the subject of the next section.

4. Canonical Formalism for Particle Dynamics

To obtain the ROGS Hamiltonian, the canonical Arnowitt-Deser-Misner (ADM) formalism [95,96,97] can be employed, as with the ADM formalism in $(3+1)$–dimensional theory. The result of this procedure is that constraints are eliminated, coordinate conditions are imposed, and the reduced Hamiltonian is given as a spatial integral of the second derivative of $\Psi$, which is a function of the dynamical variables ($z_a$, $p_a$) of the particles [98]. Consequently the Hamiltonian is completely determined in terms of the coordinates and momenta of the particles.

4.1. Neutral Particles

Consider first the N-body system whose interactions are purely gravitational. Writing the metric as

$$d{s}^2=-N_0^{2}d{t}^2+\gamma {\left(dx+\frac{N_1}{\gamma }dt\right)}^2$$

(16)

its extrinsic curvature is

$$K={\left(2N_0\gamma \right)}^{-1}(2{\partial }_1{N}_1-{\gamma }^{-1}{N}_1{\partial }_1\gamma -{\partial }_0\gamma )=\sqrt{\gamma }\kappa (\pi -\Pi /\gamma )$$

(17)

where $\pi$ and $\Pi$ are the respective conjugate momenta to $\gamma$ and $\Psi$. Using this the action integral (10) can be written as

$$I=\int d^{2}x\left\{\sum _a{p}_a{\dot{z}}_a\delta (x-z_a\left(x^0\right))+\pi \dot{\gamma }+\Pi \dot{\Psi }+N_0{R}^0+N_1{R}^1\right\}$$

(18)

where

$$\begin{array}{ccc}\hfill R^0& =& -\kappa \sqrt{\gamma }\gamma {\pi }^2+2\kappa \sqrt{\gamma }\pi \Pi +\frac{1}{4\kappa \sqrt{\gamma }}{\left({\Psi }^{\prime }\right)}^2-\frac{1}{\kappa }{\left(\frac{{\Psi }^{\prime }}{\sqrt{\gamma }}\right)}^{\prime }\hfill \end{array}$$

$$\begin{array}{ccc}& & -\frac{\Lambda }{\kappa }\sqrt{\gamma }-\sum _a\sqrt{\frac{p_a^2}{\gamma }+m_a^2}\delta (x-z_a\left(x^0\right))\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill R^1& =& \frac{{\gamma }^{\prime }}{\gamma }\pi -\frac{1}{\gamma }\Pi {\Psi }^{\prime }+2{\pi }^{\prime }+\sum _a\frac{p_a}{\gamma }\delta (x-z_a\left(x^0\right))\hfill \end{array}$$

(20)

where ${\partial }_0$ is denoted by $\left(\dot{}\right)$ and ${\partial }_1$ by a symbol $\left({}^{\prime }\right)$.

Performing the canonical reduction is a somewhat involved procedure, but is fully analogous to the $(3+1)$-dimensional case [99,100]. Suppressing the details, the coordinate conditions

$$\gamma =1\hspace{1em}\mathrm{and}\hspace{1em}\Pi =0$$

(21)

can be consistently chosen, and subsequently the action integral reduces to [98]

$$I_R=\int d{x}^2\left\{\sum _a{p}_a{\dot{z}}_a\delta (x-z_a)-\mathcal{H}\right\}.$$

(22)

where, writing $▵\Psi \equiv {\Psi }^{″}$,

$$H=\int dx\mathcal{H}=-\frac{1}{\kappa }\int dx▵\Psi$$

(23)

which can be verified using a superpotential approach based on Noether’s theorem [101]. The auxiliary field $\Psi$ is determined by solving [98,102,103]

$$▵\Psi -\frac{1}{4}{\left({\Psi }^{\prime }\right)}^2+{\kappa }^2{\pi }^2+\Lambda +\kappa \sum _a\sqrt{p_a^2+m_a^2}\delta (x-z_a)=0,$$

(24)

$$2{\pi }^{\prime }+\sum _a{p}_a\delta (x-z_a)=0$$

(25)

which follow from reducing the constraint equations (19) and (20) via the canonical reduction procedure. Note that $\Psi$ is a function of the dynamical variables ($z_a$, $p_a$) of the particles, as stated above.

The expression (23) is analogous to the reduced Hamiltonian in $(3+1)$ dimensional general relativity [97,99,100]. Equation (24) can be shown to be an energy balance condition: the energy of the gravitational field (expressed in terms of $(\Psi ,\pi )$ and the cosmological constant $\Lambda$) plus the relativistic energy of motion of the particles must add to zero. Likewise, (25) expresses the fact that the total momentum of the gravitational field and the particles must add to zero.

Two interesting approximations can be obtained from (23) and the constraint equations. One is an expansion in powers of the coupling $\kappa$ – the weak-field expansion. The other is an expansion in inverse powers of c – the post-Newtonian expansion.

The $\kappa$-expansion can be carried out by solving (24) for $▵\Psi$ and inserting the result into (23). Setting $\Lambda =0$ for simplicity, this yields

$$\begin{array}{ccc}\hfill H& =& \sum _a\sqrt{p_a^2+m_a^2}\left\{1-\frac{1}{4}\Psi \left(z_a\right)+\frac{1}{32}\Psi \left(z_a\right)\Psi \left(z_a\right)\right\}\hfill \\ & & \hspace{1em}\hspace{1em}+\frac{\kappa }{2}\sum _a{p}_a\chi \left(z_a\right)\left\{1-\frac{1}{4}\Psi \left(z_a\right)\right\}+\frac{{\kappa }^2}{8}\sum _a\sqrt{p_a^2+m_a^2}\chi \left(z_a\right)\chi \left(z_a\right)\hfill \\ & =& \sum _a\sqrt{p_a^2+m_a^2}+\frac{\kappa }{8}\sum _a\sum _b\left(\sqrt{p_a^2+m_a^2}\sqrt{p_b^2+m_b^2}-p_a{p}_b\right)r_{ab}\hfill \\ & & +\frac{ϵ\kappa }{8}\sum _a\sum _b\left(\sqrt{p_a^2+m_a^2}{p}_b-p_a\sqrt{p_b^2+m_b^2}\right)(z_a-z_b)\hfill \\ & & +\frac{1}{4}{\left(\frac{\kappa }{4}\right)}^2\left\{\sum _a\sqrt{p_a^2+m_a^2}{\left[\sum _b{p}_b{r}_{ab}+ϵ\sum _b\sqrt{p_b^2+m_b^2}(z_a-z_b)\right]}^2\right.\hfill \\ & & -\sum _a{p}_a\left[\sum _b{p}_b{r}_{ab}+ϵ\sum _b\sqrt{p_b^2+m_b^2}(z_a-z_b)\right]\left[\sum _c\sqrt{p_c^2+m_c^2}{r}_{ac}+ϵ\sum _c{p}_c(z_a-z_c)\right]\hfill \\ & & +\sum _a\sum _b\left[\sqrt{p_a^2+m_a^2}\sqrt{p_b^2+m_b^2}{r}_{ab}-ϵp_a\sqrt{p_b^2+m_b^2}(z_a-z_b)\right]\hfill \\ & & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times \left[\sum _c\sqrt{p_c^2+m_c^2}{r}_{bc}+ϵ\sum _c{p}_c(z_b-z_c)\right]\hfill \\ & & \left.-\sum _{a,b}\left[\sqrt{p_a^2+m_a^2}{p}_b{r}_{ab}-ϵp_a{p}_b(z_a-z_b)\right]\left[\sum _c{p}_c{r}_{bc}+ϵ\sum _c\sqrt{p_c^2+m_c^2}(z_b-z_c)\right]\right\}\hfill \end{array}$$

(26)

to second order in $\kappa$, where $\chi$ is defined by ${\chi }^{\prime }\equiv \pi$. In order to obtain (26) the boundary term

$$S_b={\left[-\frac{1}{4\kappa }\Psi {\Psi }^{\prime }+\kappa \chi {\chi }^{\prime }-\frac{\kappa }{8}\left\{\Psi {\left({\chi }^2\right)}^{\prime }-{\Psi }^{\prime }{\chi }^2\right\}+\frac{1}{32\kappa }{\Psi }^2{\Psi }^{\prime }\right]}_{-\infty }^{\infty }$$

(27)

must vanish. Here there is a subtle problem in comparison to the $(3+1)$–dimensional setting because the dimensionless potential $Gm|x|/c^2$ becomes infinite at spatial infinity. It is straightforward to show that if the function $f\left(x\right)\equiv {\Psi }^2-4{\kappa }^2{\chi }^2$ and its first derivative vanish for $|x|>>|z_a|$ for all a then the surface term (27) vanishes, since it is a sum of terms proportional to $f^{\prime }$ and $(4-\Psi )f^{\prime }+f{\Psi }^{\prime }$.

The $\kappa$-expansion is the successive approximation of the ROGS in the background of Minkowskian space-time. At each order in $\kappa$ the relativistic form is preserved, and so this approximation is appropriate for describing relativistic motion of the particles in weak gravitational fields.

The post-Newtonian expansion is an approximation powers of $1/c$. Temporarily restoring $t\to ct$ we note that both $p_a^2/m_a^2$ and $\sqrt{\kappa }$ are of the order of $c^{-2}$. This yields

$$\begin{array}{ccc}\hfill H& =& \sum _a{m}_a{c}^2+\sum _a\frac{{\tilde{p}}_a^2}{2m_a}+\frac{\kappa }{8}\sum _a\sum _b{m}_a{m}_b{\tilde{r}}_{ab}\hfill \\ & & -\sum _a\frac{{\tilde{p}}_a^4}{8m_a^3{c}^2}+\frac{\kappa c^2}{8}\sum _a\sum _b{m}_a\frac{{\tilde{p}}_b^2}{m_b}{\tilde{r}}_{ab}-\frac{\kappa c^2}{8}\sum _a\sum _b{\tilde{p}}_a{\tilde{p}}_b{\tilde{r}}_{ab}\hfill \\ & & +\frac{1}{4}{\left(\frac{\kappa c^3}{4}\right)}^2\sum _a\sum _b\sum _c{m}_a{m}_b{m}_c\left\{{\tilde{r}}_{ab}{\tilde{r}}_{ac}-({\tilde{z}}_a-{\tilde{z}}_b)({\tilde{z}}_a-{\tilde{z}}_c)\right\}\hfill \end{array}$$

(28)

from (26), where the canonical variables $(z_a,p_a)$ have been redefined

$$\begin{array}{c}\hfill z_a⟶{\tilde{z}}_a=z_a{p}_a⟶{\tilde{p}}_a=p_a-\frac{ϵ\kappa }{4}\sum _b{m}_a{m}_b{r}_{ab}\end{array}$$

(29)

to eliminate a spurious term of order $1/c$, with $r_{ab}=(z_a-z_b)$. Under this redefinition the Poisson brackets amongst the canonical variables remain unchanged. It is straightforward to show that the canonical equations of motion

$$\begin{array}{c}\hfill {\dot{\tilde{z}}}_a=\frac{\partial H}{\partial {\tilde{p}}_a}{\dot{\tilde{p}}}_a=-\frac{\partial H}{\partial {\tilde{z}}_a}\end{array}$$

(30)

are equivalent to the geodesic equation (15) [98].

For illustrative purposes, consider a single static source at the origin with $\Lambda =0$. The constraint equations (24) and (25) become

$$\begin{array}{ccc}\hfill -\frac{1}{\kappa }▵\Psi & =& \kappa {\pi }^2-\frac{1}{4\kappa }{\left({\Psi }^{\prime }\right)}^2+M\delta \left(x\right)\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill {\pi }^{\prime }& =& 0\hfill \end{array}$$

(32)

and the Hamiltonian equations that follow from (18) are

$$\begin{array}{ccc}\hfill \dot{\pi }+N_0\left[\frac{3\kappa }{2}{\pi }^2+\frac{1}{8\kappa }{\left({\Psi }^{\prime }\right)}^2\right]+N_1{\pi }^{\prime }+\frac{1}{2\kappa }{N}_0^{\prime }{\Psi }^{\prime }+N_1^{\prime }\pi & =& 0\hfill \end{array}$$

(33)

$$\begin{array}{ccc}\hfill \dot{\Psi }+2\kappa N_0\pi -N_1{\Psi }^{\prime }& =& 0\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill \kappa \pi N_0+N_1^{\prime }& =& 0\hfill \end{array}$$

(35)

$$\begin{array}{ccc}\hfill {\partial }_1\left(\frac{1}{2}{N}_0{\Psi }^{\prime }+N_0^{\prime }\right)& =& 0\hfill \end{array}$$

(36)

Solutions to these equations that ensure the boundary conditions $S_b=0$ hold are

$$\begin{array}{ccc}\hfill \Psi & =& -\frac{\kappa M}{2}|x|+\frac{ϵ\kappa M}{2}t\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill \pi & =& {\chi }^{\prime }=-\frac{ϵ}{4}M\hspace{1em}\hspace{1em}\hspace{1em}(\chi =-\frac{ϵM}{4}x)\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill N_0& =& e^{\frac{\kappa M}{4}|x|}\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill N_1& =& ϵ\frac{x}{|x|}\left(e^{\frac{\kappa M}{4}|x|}-1\right)\hfill \end{array}$$

(40)

where $ϵ$ is a constant of integration, with ${ϵ}^2=1$. Solutions with $ϵ=1$ are related by time-reversal to those with $ϵ=-1$. This factor guarantees the invariance of the whole theory under the time reversal. We also see from (37) that even for a static source the auxiliary field is not static.

The Hamiltonian is

$$\begin{array}{ccc}\hfill H& =& -\frac{1}{\kappa }\int dx▵\Psi =M\hfill \end{array}$$

(41)

showing that the total energy is the mass, and

$$\begin{array}{ccc}\hfill g_{00}& =& -N_0^2+N_1^2=1-2e^{\frac{\kappa M}{4}|x|}<-1\hfill \end{array}$$

indicating that this static point-particle solution has no event horizon.

4.2. Charged Particles

It is also possible to include non-gravitational interactions into the ROGS. Electromagnetism is an obvious choice to consider, and in fact the derviation of the canonically reduced Hamiltonian for charged particles is parallel to that of the previous subsection for the uncharged case [104]. Here I shall highlight aspects of the charged case that are distinct from the uncharged case.

The action integral for gravitational and electromagnetic fields coupled with N charged point masses is

$$\begin{array}{ccc}\hfill I& =& \int d^{2}x\left[\frac{1}{2\kappa }\sqrt{-g}{g}^{\mu \nu }\left\{\Psi R_{\mu \nu }+\frac{1}{2}{\nabla }_{\mu }\Psi {\nabla }_{\nu }\Psi -g_{\mu \nu }\Lambda \right\}\right.\hfill \\ & & \hspace{1em}\hspace{1em}-A_{\mu }{\mathcal{F}}_{,\nu }^{\mu \nu }+\frac{1}{4\sqrt{-g}}{\mathcal{F}}^{\mu \nu }{\mathcal{F}}^{\alpha \beta }{g}_{\mu \alpha }{g}_{\nu \beta }\hfill \\ & & \hspace{1em}\hspace{1em}\left.+\sum _a\int d{\tau }_a\left\{-m_a{\left(-g_{\mu \nu }\left(x\right)\frac{d{z}_a^{\mu }}{d{\tau }_a}\frac{d{z}_a^{\nu }}{d{\tau }_a}\right)}^{1/2}+e_a\frac{d{z}_a^{\mu }}{d{\tau }_a}{A}_{\mu }\left(x\right)\right\}{\delta }^2(x-z_a\left({\tau }_a\right))\right]\hfill \end{array}$$

(42)

where $A_{\mu }$ and ${\mathcal{F}}^{\mu \nu }$ are the vector potential and the field strength, and $e_a$ and ${\tau }_a$ are the charge and the proper time of a-th particle, respectively. The field equations that follow from the action (42) are

$$\begin{array}{ccc}& & R-g^{\mu \nu }{\nabla }_{\mu }{\nabla }_{\nu }\Psi =0\hfill \end{array}$$

(43)

$$\begin{array}{ccc}& & \frac{1}{2}{\nabla }_{\mu }\Psi {\nabla }_{\nu }\Psi -\frac{1}{4}{g}_{\mu \nu }{\nabla }^{\lambda }\Psi {\nabla }_{\lambda }\Psi +g_{\mu \nu }{\nabla }^{\lambda }{\nabla }_{\lambda }\Psi -{\nabla }_{\mu }{\nabla }_{\nu }\Psi =\kappa T_{\mu \nu }-g_{\mu \nu }\Lambda \hfill \end{array}$$

(44)

$$\begin{array}{ccc}& & {\mathcal{F}}_{,\nu }^{\mu \nu }=\sum _a{e}_a\int d{\tau }_a\frac{d{z}_a^{\mu }}{d{\tau }_a}{\delta }^2(x-z_a\left({\tau }_a\right))\hfill \end{array}$$

(45)

$$\begin{array}{ccc}& & \frac{1}{\sqrt{-g}}{\mathcal{F}}_{\mu \nu }={\partial }_{\mu }{A}_{\nu }-{\partial }_{\nu }{A}_{\mu },\hfill \end{array}$$

(46)

$$\begin{array}{ccc}& & m_a\left[\frac{d}{d{\tau }_a}\left\{g_{\mu \nu }\left(z_a\right)\frac{d{z}_a^{\nu }}{d{\tau }_a}\right\}-\frac{1}{2}{g}_{\nu \lambda ,\mu }\left(z_a\right)\frac{d{z}_a^{\nu }}{d{\tau }_a}\frac{d{z}_a^{\lambda }}{d{\tau }_a}\right]=e_a\frac{d{z}_a^{\nu }}{d{\tau }_a}\left\{A_{\nu ,\mu }\left(z_a\right)-A_{\mu ,\nu }\left(z_a\right)\right\}\hfill \end{array}$$

(47)

where the stress-energy is

$$T_{\mu \nu }=\sum _a{m}_a\int d{\tau }_a\frac{1}{\sqrt{-g}}{g}_{\mu \sigma }{g}_{\nu \rho }\frac{d{z}_a^{\sigma }}{d{\tau }_a}\frac{d{z}_a^{\rho }}{d{\tau }_a}{\delta }^2(x-z_a\left({\tau }_a\right))+\frac{1}{(-g)}\left\{{\mathcal{F}}_{\mu \alpha }{\mathcal{F}}_{\nu \beta }{g}^{\alpha \beta }-\frac{1}{4}{g}_{\mu \nu }{\mathcal{F}}_{\alpha \beta }{\mathcal{F}}^{\alpha \beta }\right\}$$

(48)

whose conservation is guaranteed from (44). Note that in (1+1) dimensions no magnetic component of the field exists. Inserting the trace of (44) into (43) yields (14) as before. This equation, along with (45), (46), and (47) form a closed system of equations for gravity, electromagnetism, and charged N-body matter.

Using the form (16) of the metric, the variational principle yields the set of equations

$$\begin{array}{ccc}\hfill \dot{\pi }& +& N_0\left\{\frac{3\kappa }{2}\sqrt{\gamma }{\pi }^2-\frac{\kappa }{\sqrt{\gamma }}\pi \Pi +\frac{1}{8\kappa \sqrt{\gamma }\gamma }{\left({\Psi }^{\prime }\right)}^2+\frac{1}{4\sqrt{\gamma }}(E^2+2\frac{\Lambda }{\kappa })\right.\hfill \\ & & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left.-\sum _a\frac{p_a^2}{2{\gamma }^2\sqrt{\frac{p_a^2}{\gamma }+m_a^2}}\delta (x-z_a\left(t\right))\right\}\hfill \\ & +& N_1\left\{-\frac{1}{{\gamma }^2}\Pi {\Psi }^{\prime }+\frac{{\pi }^{\prime }}{\gamma }+\sum _a\frac{p_a}{{\gamma }^2}\delta (x-z_a\left(t\right))\right\}+N_0^{\prime }\frac{1}{2\kappa \sqrt{\gamma }\gamma }{\Psi }^{\prime }+N_1^{\prime }\frac{\pi }{\gamma }=0,\hfill \end{array}$$

(49)

$$\begin{array}{ccc}& & \dot{\gamma }-N_0(2\kappa \sqrt{\gamma }\gamma \pi -2\kappa \sqrt{\gamma }\Pi )+N_1\frac{{\gamma }^{\prime }}{\gamma }-2N_1^{\prime }=0,\hfill \end{array}$$

(50)

$$\begin{array}{ccc}& & R^0=0,\hfill \end{array}$$

(51)

$$\begin{array}{ccc}& & R^1=0,\hfill \end{array}$$

(52)

$$\begin{array}{ccc}& & \dot{\Pi }+{\partial }_1(-\frac{1}{\gamma }{N}_1\Pi +\frac{1}{2\kappa \sqrt{\gamma }}{N}_0{\Psi }^{\prime }+\frac{1}{\kappa \sqrt{\gamma }}{N}_0^{\prime })=0,\hfill \end{array}$$

(53)

$$\begin{array}{ccc}& & \dot{\Psi }+N_0\left(2\kappa \sqrt{\gamma }\pi \right)-N_1(\frac{1}{\gamma }{\Psi }^{\prime })=0,\hfill \\ & & {\dot{p}}_a+\frac{\partial N_0}{\partial z_a}\sqrt{\frac{p_a^2}{\gamma }+m_a^2}-\frac{N_0}{2\sqrt{\frac{p_a^2}{\gamma }+m_a^2}}\frac{p_a^2}{{\gamma }^2}\frac{\partial \gamma }{\partial z_a}-\frac{\partial N_1}{\partial z_a}\frac{p_a}{\gamma }\hfill \end{array}$$

(54)

$$\begin{array}{ccc}& & +N_1\frac{p_a}{{\gamma }^2}\frac{\partial \gamma }{\partial z_a}+\int dx{N}_0\sqrt{\gamma }E\frac{\partial E}{\partial z_a}=0,\hfill \end{array}$$

(55)

$$\begin{array}{ccc}& & \dot{z_a}-N_0\frac{\frac{p_a}{\gamma }}{\sqrt{\frac{p_a^2}{\gamma }+m_a^2}}+\frac{N_1}{\gamma }=0.\hfill \end{array}$$

(56)

$$\begin{array}{ccc}& & \frac{\partial E}{\partial t}+\sum _a{e}_a\frac{d{z}_a}{dt}\delta (x-z_a\left(t\right))=0,\hfill \end{array}$$

(57)

$$\begin{array}{ccc}& & \sqrt{-g}E=-\frac{\partial \phi }{\partial x}-\frac{\partial A}{\partial t}\hfill \end{array}$$

(58)

where

$$A_{\mu }=(-\phi ,A)E={\mathcal{F}}^{01}$$

(59)

and the quantities

$$\begin{array}{ccc}\hfill R^0& =& -\kappa \sqrt{\gamma }\gamma {\pi }^2+2\kappa \sqrt{\gamma }\pi \Pi +\frac{1}{4\kappa \sqrt{\gamma }}{\left({\Psi }^{\prime }\right)}^2-\frac{1}{\kappa }{\left(\frac{{\Psi }^{\prime }}{\sqrt{\gamma }}\right)}^{\prime }-\frac{1}{2}\sqrt{\gamma }(E^2+2\frac{\Lambda }{\kappa })\hfill \end{array}$$

$$\begin{array}{ccc}& & -\sum _a\sqrt{\frac{p_a^2}{\gamma }+m_a^2}\delta (x-z_a\left(t\right))\hfill \end{array}$$

(60)

$$\begin{array}{ccc}\hfill R^1& =& \frac{{\gamma }^{\prime }}{\gamma }\pi -\frac{1}{\gamma }\Pi {\Psi }^{\prime }+2{\pi }^{\prime }+\sum _a\frac{p_a}{\gamma }\delta (x-z_a\left(t\right))\hfill \end{array}$$

(61)

$$\begin{array}{ccc}\hfill \frac{\partial E}{\partial x}& =& \sum _a{e}_a\delta (x-z_a\left(t\right))\hfill \end{array}$$

(62)

vanish due to the constraints. The set of equations (49) – (58) are equivalent to the set of Equations (43) – (47).

Note that in (55) and (56), all metric components ($N_0$, $N_1$, $\gamma$) are evaluated at the point $x=z_a$ with

$$\frac{\partial f}{\partial z_a}\equiv {\left.\frac{\partial f\left(x\right)}{\partial x}\right|}_{x=z_a}.$$

Furthermore, in (1+1) dimensions the electric field has no independent degrees of freedom. For charged particles moving within a finite region $(|x|<L)$, the electric field in an outer region $(|x|>L)$ is $E=\pm \frac{1}{2}{\sum }_a{e}_a$. $E=0$ in the outer region for a system of zero total charge.

The canonical reduction of this system proceeds in a manner similar to the uncharged case. The coordinate conditions (21) can again be chosen and the reduced action is still given by (22). The reduced Hamiltonian is again given by (23), but now the constraints become

$$▵\Psi -\frac{1}{4}{\left({\Psi }^{\prime }\right)}^2+{\kappa }^2{\pi }^2+\frac{1}{2}(\kappa E^2+2\Lambda )+\kappa \sum _a\sqrt{p_a^2+m_a^2}\delta (x-z_a)=0,$$

(63)

$$2{\pi }^{\prime }+\sum _a{p}_a\delta (x-z_a)=0.$$

(64)

The consistency of this canonical reduction is proven by showing that the canonical equations of motion derived from the reduced Hamiltonian (23) are identical with the Equations (55) and (56) [98,104].

5. The 2-Body Problem

As stated in the introduction, an exact solution to the two body problem is known in Newtonian theory, but not in general relativity due to dissipation of energy in the form of gravitational radiation. Hence analysis of a two body system in general relativity (e.g., binary pulsars) necessarily involves resorting to approximation methods such as a post-Newtonian expansion [99,105]. Quite remarkably, there are a number of exact solutions [103,104,104,106,107,108] to the 2-body problem in the ROGS theory described above. While this is in part understandable due to the lack of gravitational radiation (there are no gravitons in $(1+1)$ dimensions), the non-linearity of the system does not obviously admit exact solutions, and in fact a general dilaton theory of gravity will not have them.

5.1. Solution for 2 Charged Particles

Here we present the exact solution for the charged case [104], since from it the electrically neutral case is recovered by setting $e_1=e_2=0$. The standard approach for obtaining a solution is to find an explicit expression for the Hamiltonian (23) by solving (63) and (64). From this the equations of motion can be derived, which in turn can be solved to obtain the trajectories of the particles.

Noting that the electric field appears in the combination $(E^2+2\Lambda /\kappa )$ in all equations, it is convenient to set

$$V\left(x\right)\equiv E^2-C\mathrm{and}{\Lambda }_e\equiv -2\Lambda -\kappa C$$

(65)

with $C\equiv \frac{1}{4}{\left({\sum }_a{e}_a\right)}^2$. The quantity ${\Lambda }_e$ is an effective cosmological constant that includes the contribution from the electric field. This situation arises because, analogous to the way in which a 4-form behaves in $(3+1)$ dimensions [109], in $(1+1)$ dimensions the electromagnetic field strength is a 2-form, and so contributes to the stress-energy tensor in the same manner as a cosmological constant in compact spatial regions.

Defining $\varphi$ by

$$\Psi =-4ln|\varphi |,$$

(66)

the constraints (63) and (64) for a two-particle system become

$$\begin{array}{ccc}\hfill ▵\varphi -\frac{1}{4}\left\{{\kappa }^2{\left({\chi }^{\prime }\right)}^2+\frac{\kappa }{2}V-\frac{1}{2}{\Lambda }_e\right\}\varphi & =& \frac{\kappa }{4}\left\{\sqrt{p_1^2+m_1^2}\varphi \left(z_1\right)\delta (x-z_1)\right.\hfill \\ & & \left.+\sqrt{p_2^2+m_2^2}\varphi \left(z_2\right)\delta (x-z_2)\right\}\hfill \end{array}$$

(67)

$$▵\chi =-\frac{1}{2}\left\{p_1\delta (x-z_1)+p_2\delta (x-z_2)\right\}$$

(68)

along with (62) for the electric field, whose solution is

$$E=\frac{1}{2}{\partial }_x\sum _a{e}_a|x-z_a\left(t\right)|.$$

(69)

It is also straightforward to solve (68); the result is

$$\chi =-\frac{1}{4}\left\{p_1|x-z_1|+p_2|x-z_2|\right\}-ϵXx+ϵC_{\chi }$$

(70)

where X and $C_{\chi }$ are constants of integration and, as in (37), $ϵ$ (${ϵ}^2=1$) ensures that the time reversal properties of $\chi$ are explicitly manifest. By definition, $ϵ$ changes sign under time reversal and thus so does $\chi$.

Suppose that $z_2<z_1$. We can divide space into three regions: $z_1<x$ ((+) region), $z_2<x<z_1$ ((0) region) and $x<z_2$ ((-) region). From (69) and (70) one sees that in each region ${\chi }^{\prime }$ and V are constant:

$$V=\left\{\begin{array}{cc}0\hfill & (+)\mathrm{region},\hfill \\ -\frac{1}{2}{e}_1{e}_2\hfill & \left(0\right)\mathrm{region},\hfill \\ 0\hfill & (-)\mathrm{region},\hfill \end{array}\right.$$

(71)

$${\chi }^{\prime }=\left\{\begin{array}{cc}-ϵX-\frac{1}{4}(p_1+p_2)\hfill & \hspace{1em}\hspace{1em}\hspace{1em}(+)\mathrm{region},\hfill \\ -ϵX+\frac{1}{4}(p_1-p_2)\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\left(0\right)\mathrm{region},\hfill \\ -ϵX+\frac{1}{4}(p_1+p_2)\hfill & \hspace{1em}\hspace{1em}\hspace{1em}(-)\mathrm{region}\hfill \end{array}\right.$$

(72)

Obtaining the solution to (67) is somewhat more complicated [103], entailing a proper choice of boundary conditions that ensure finiteness of the Hamiltonian and the vanishing of surface terms which arise in transforming the action. The result is [104]

$$\begin{array}{ccc}\hfill {\varphi }_{+}& =& {\left(\frac{K_1}{{\mathcal{M}}_1}\right)}^{\frac{1}{2}}{e}^{-\frac{1}{4}(K_{01}{z}_1-K_{02}{z}_2)+\frac{1}{2}{K}_{+}(x-z_1)}\hfill \\ \hfill {\varphi }_0& =& \frac{e^{-\frac{1}{4}(K_{01}{z}_1-K_{02}{z}_2)}}{4K_0}\left\{{\left(K_1{\mathcal{M}}_1\right)}^{1/2}{e}^{-\frac{1}{2}{K}_0(x-z_1)}+{\left(K_2{\mathcal{M}}_2\right)}^{1/2}{e}^{\frac{1}{2}{K}_0(x-z_2)}\right\}\hfill \\ \hfill {\varphi }_{-}& =& {\left(\frac{K_2}{{\mathcal{M}}_2}\right)}^{\frac{1}{2}}{e}^{-\frac{1}{4}(K_{01}{z}_1-K_{02}{z}_2)-\frac{1}{2}{K}_{-}(x-z_2)}\hfill \end{array}$$

(73)

where

$$\left\{\begin{array}{ccc}{K}_{+}\hfill & =\sqrt{{\kappa }^2{\left(X+\frac{ϵ}{4}({\tilde{p}}_1+{\tilde{p}}_2)\right)}^2-\frac{1}{2}{\Lambda }_e}\hfill & (+)\mathrm{region}\hfill \\ K_0\hfill & =\sqrt{{\kappa }^2{\left(X-\frac{ϵ}{4}({\tilde{p}}_1-{\tilde{p}}_2)\right)}^2-\frac{\kappa }{2}{e}_1{e}_2-\frac{1}{2}{\Lambda }_e}\hfill & \left(0\right)\mathrm{region}\hfill \\ K_{-}\hfill & =\sqrt{{\kappa }^2{\left(X-\frac{ϵ}{4}({\tilde{p}}_1+{\tilde{p}}_2)\right)}^2-\frac{1}{2}{\Lambda }_e}\hfill & (-)\mathrm{region}\hfill \end{array}\right.$$

(74)

with ${\tilde{p}}_i=p_i\mathrm{sgn}(z_1-z_2)$, and

$$\begin{array}{ccc}\hfill K_1& \equiv & 2K_0+2K_{-}-\kappa \sqrt{p_2^2+m_2^2}\hfill \\ \hfill K_2& \equiv & 2K_0+2K_{+}-\kappa \sqrt{p_1^2+m_1^2}\hfill \\ \hfill K_{01}& \equiv & K_0-K_{+}+\frac{\kappa ϵ}{2}{\tilde{p}}_1,\hfill \\ \hfill K_{02}& \equiv & K_0-K_{-}-\frac{\kappa ϵ}{2}{\tilde{p}}_2,\hfill \\ \hfill {\mathcal{M}}_1& \equiv & \kappa \sqrt{p_1^2+m_1^2}+2K_0-2K_{+}\hfill \\ \hfill {\mathcal{M}}_2& \equiv & \kappa \sqrt{p_2^2+m_2^2}+2K_0-2K_{-}\hfill \end{array}$$

(75)

The Hamiltonian (23) becomes

$$\begin{array}{ccc}\hfill H& =& -\frac{1}{\kappa }\int dx▵\Psi =\frac{4}{\kappa }{\left[\frac{{\varphi }^{\prime }}{\varphi }\right]}_{-\infty }^{\infty }=\frac{2(K_{+}+K_{-})}{\kappa }\hfill \end{array}$$

(76)

and is explicitly determined in terms of the canonical variables once the solution X to

$$K_1{K}_2={\mathcal{M}}_1{\mathcal{M}}_2{e}^{K_0|z_1-z_2|}$$

(77)

is obtained. This relation is a (not obvious) consequence of solving (67), and can more explicitly be written as

$$\begin{array}{ccc}\hfill & & \left(4K_0^2+[\kappa \sqrt{p_1^2+m_1^2}-2K_{+}][\kappa \sqrt{p_2^2+m_2^2}-2K_{-}]\right)\tanh \left(\frac{1}{2}{K}_0|z_1-z_2|\right)\hfill \\ & & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=-2K_0\left([\kappa \sqrt{p_1^2+m_1^2}-2K_{+}]+[\kappa \sqrt{p_2^2+m_2^2}-2K_{-}]\right)\hfill \end{array}$$

(78)

which is the determining equation of the Hamiltonian.

Note that $K_{\pm }$ must both be real in order that the Hamiltonian (76) have a definite meaning. This imposes a restriction on X corresponding to a value of the effective cosmological constant ${\Lambda }_e$. But $K_0$ is not necessarily real; indeed for sufficiently strong electromagnetic repulsion (sufficiently large positive $e_1{e}_2$) $K_0$ will be imaginary. In this case in the (0) region the solution for $\varphi$ becomes

$${\varphi }_0\left(x\right)=A_s\sin \frac{1}{2}{\tilde{K}}_{0}x+A_c\cos \frac{1}{2}{\tilde{K}}_{0}x,$$

(79)

instead of (73), where

$$\begin{array}{ccc}\hfill {\tilde{K}}_0& =& -i{K}_0=\sqrt{\frac{\kappa }{2}{e}_1{e}_2+\frac{1}{2}{\Lambda }_e-{\kappa }^2{\left(X-\frac{ϵ}{4}({\tilde{p}}_1-{\tilde{p}}_2)\right)}^2}\hfill \end{array}$$

(80)

yielding a new determining equation of the Hamiltonian

$$\begin{array}{ccc}& & \left(4{\tilde{K}}_0^2-[\kappa \sqrt{p_1^2+m_1^2}-2K_{+}][\kappa \sqrt{p_2^2+m_2^2}-2K_{-}]\right)\tan \left(\frac{1}{2}{\tilde{K}}_0|z_1-z_2|\right)\hfill \\ & & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=2{\tilde{K}}_0\left([\kappa \sqrt{p_1^2+m_1^2}-2K_{+}]+[\kappa \sqrt{p_2^2+m_2^2}-2K_{-}]\right)\hfill \end{array}$$

(81)

instead of (78). Indeed (81) and (79) can respectively be obtained from (78) and (73) by formally replacing $K_0$ with $i{\tilde{K}}_0$.

Consequently (77) is applicable for all values of $K_0$. This transcendental equation determines H in terms of the momenta and positions of the particles. Although (77) cannot be explicitly written in terms of known functions, it can be used to obtain exactly the canonical equations of motion for $z_i$ and $p_i$, which are [104]

$$\begin{array}{ccc}\hfill {\dot{p}}_1& =& -\frac{\partial H}{\partial z_1}=-\frac{2}{\kappa }\left(\frac{Y_{+}}{K_{+}}+\frac{Y_{-}}{K_{-}}\right)\frac{K_0{K}_1{K}_2}{J}\hfill \end{array}$$

(82)

$$\begin{array}{ccc}\hfill {\dot{z}}_1& =& \frac{\partial H}{\partial p_1}=ϵ\frac{Y_{+}}{K_{+}}+\frac{8}{J}\left(\frac{Y_{+}}{K_{+}}+\frac{Y_{-}}{K_{-}}\right)\frac{K_0{K}_1}{{\mathcal{M}}_1}\left\{\frac{p_1}{\sqrt{p_1^2+m_1^2}}-ϵ\frac{Y_{+}}{K_{+}}\right\}\hfill \end{array}$$

(83)

$$\begin{array}{ccc}\hfill {\dot{p}}_2& =& \frac{2}{\kappa }\left(\frac{Y_{+}}{K_{+}}+\frac{Y_{-}}{K_{-}}\right)\frac{K_0{K}_1{K}_2}{J}\hfill \end{array}$$

(84)

$$\begin{array}{ccc}\hfill {\dot{z}}_2& =& -ϵ\frac{Y_{-}}{K_{-}}+\frac{8}{J}\left(\frac{Y_{+}}{K_{+}}+\frac{Y_{-}}{K_{-}}\right)\frac{K_0{K}_2}{{\mathcal{M}}_2}\left\{\frac{p_2}{\sqrt{p_2^2+m_2^2}}+ϵ\frac{Y_{-}}{K_{-}}\right\}\hfill \end{array}$$

(85)

where

$$\begin{array}{ccc}\hfill Y_{\pm }& \equiv & \kappa \left[X\pm \frac{ϵ}{4}(p_1+p_2)\right]\hfill \\ \hfill Y_0& \equiv & \kappa \left[X-\frac{ϵ}{4}(p_1-p_2)\right]\hfill \end{array}$$

(86)

and

$$\begin{array}{ccc}\hfill J& =& 2\left\{\left(\frac{Y_0}{K_0}+\frac{Y_{+}}{K_{+}}\right)K_1+\left(\frac{Y_0}{K_0}+\frac{Y_{-}}{K_{-}}\right)K_2\right\}\hfill \\ & & -2\left\{\left(\frac{Y_0}{K_0}-\frac{Y_{+}}{K_{+}}\right)\frac{1}{{\mathcal{M}}_1}+\left(\frac{Y_0}{K_0}-\frac{Y_{-}}{K_{-}}\right)\frac{1}{{\mathcal{M}}_2}\right\}{K}_1{K}_2-\frac{Y_0}{K_0}{K}_1{K}_2(z_1-z_2).\hfill \end{array}$$

(87)

These canonical equations ensure the conservation of the Hamiltonian and the total momentum $p_1+p_2$. They can also be shown [104] to be equivalent to the geodesic Equations (55) and (56), which become

$$\begin{array}{ccc}\hfill {\dot{p}}_a& =& -\frac{\partial N_0}{\partial z_a}\sqrt{p_a^2+m_a^2}+\frac{\partial N_1}{\partial z_a}{p}_a+\frac{1}{2}\sum _b{e}_a{e}_b{N}_0\frac{\partial }{\partial z_a}|z_a-z_b|,\hfill \end{array}$$

(88)

$$\begin{array}{ccc}\hfill \dot{z_a}& =& N_0\frac{p_a}{\sqrt{p_a^2+m_a^2}}-N_1.\hfill \end{array}$$

(89)

under the coordinate conditions (21), where

$$\frac{\partial N_{0,1}}{\partial z_i}\equiv \frac{1}{2}\left\{{\left.\frac{\partial N_{0,1}}{\partial x}\right|}_{x=z_i+0}+{\left.\frac{\partial N_{0,1}}{\partial x}\right|}_{x=z_i-0}\right\}.$$

(90)

defines the partial derivatives at $z_1,z_2$. Equations (49), (50), (53) and (54) yield the metric under the coordinate conditions (21) [104].

5.2. Test Particle Limit

Figure 1. A plot of the relativistic trajectory (blue) of a test particle in comparison to the non-relativistic case (red). The Hamiltonian for each case is set to $H=1.125m$ and $\mu =0.1m$.

Figure 1. A plot of the relativistic trajectory (blue) of a test particle in comparison to the non-relativistic case (red). The Hamiltonian for each case is set to $H=1.125m$ and $\mu =0.1m$.

To make sense of the solutions contained in (77) it is useful to first consider the test particle limit in which $m_1=\mu <<m=m_2$ and both particles have zero charge. Setting the latter to be a static source at the origin, with $z_2=0=p_2$, the defining Equation (77) becomes

$$(\sqrt{p^2+{\mu }^2}-H)(m+ϵ\tilde{p}-H)=(\sqrt{p^2+{\mu }^2}-ϵ\tilde{p})m{e}^{\frac{\kappa }{4}(H-ϵ\tilde{p})|z|}$$

(91)

where $z_1=z$ and $p_1=p$. The corresponding non-relativistic Hamiltonian is

$$H=m+\mu +\frac{\kappa m\mu }{4}|z|+\frac{p^2}{2\mu }$$

(92)

obtained by expanding (91) for $p<<m$, $\kappa m|z|<<1$. Typically the sum $m{c}^2+\mu c^2$ of the rest-energies (where we have set $c=1$ in the above) is subtracted from the non-relativistic Hamiltonian (92). However for proper comparison of the energy to the relativisitic case these terms need to be retained.

Figure 1 shows a comparison between the relativistic case (91) and non-relativistic case (92) in the phase space of the test particle. Even for relatively small values of the total energy, the non-relativistic motion (red) is notably distinct. For the given value of $H=1.125m$ and an initial condition of $z=0$, the initial momentum $p=0.075m$ in the relativistic case as compared to $p=0.071m$ in the non-relativisitic case. In both cases the static mass m attracts the test particle, but in the relativistic case its momentum initially decreases less rapidly. However this situation changes once the particle reaches its maximal separation from m. This occurs at $\kappa m|z|=1$ in the non-relativistic case, but at the shorter separation $\kappa m|z|=0.89$ in the relativistic case – relativistic gravitational attraction is stronger for the same total energy. One can observe the rather counter-intuitive feature that the particle is then attracted back to the source even though it has positive momentum in the relativistic case! The loss of momentum is more rapid than in the non-relativistic case, and continues to be so until the test particle returns to the origin, after which (in both cases) the motion repeats.

5.3. Exact Equal Mass Two-Body Motion for $\Lambda =0$

For neutral bodies of equal mass the determining Equation (77) can be solved explicitly for $\Lambda =0$, yielding [106,107]

$$H=\sqrt{p^2+m^2}+ϵp\mathrm{sgn}\left(r\right)-\frac{8W\left[-\frac{\kappa }{8}(|r|\sqrt{p^2+m^2}-ϵpr)\exp \left(\frac{\kappa }{8}(|r|\sqrt{p^2+m^2}-ϵpr\right)\right]}{\kappa |r|}$$

(93)

in the center of inertia system $p_1=-p_2=p$, where $r=z_1-z_2$. Here $W\left(x\right)$ is the Lambert W function [110]

$$y·e^y=x\Rightarrow y=W\left(x\right)$$

(94)

which has two real branches $W_0$ and $W_{-1}$ for real x.

Since the Hamiltonian (93) is exact for arbitrary values of m and p to infinite order in $\kappa$ the whole structure of the theory can be studied, from the weak field to the strong field limits. As in the test particle case, a phase space trajectory in $(r,p)$ space can be obtained by setting $H=H_0$. Alternatively, the equations of motion (82), (83), (84) and (85) yield

$$\begin{array}{ccc}\hfill \dot{p}& =& -\frac{\frac{\kappa }{4}(H-2ϵ\tilde{p})(H-ϵ\tilde{p}-\sqrt{p^2+m^2})}{2-\frac{\kappa r}{4}(H-ϵ\tilde{p}-\sqrt{p^2+m^2})}\mathrm{sgn}\left(r\right),\hfill \end{array}$$

(95)

$$\begin{array}{ccc}\hfill \dot{r}& =& 2ϵ\left\{1-\frac{H-2ϵ\tilde{p}}{2-\frac{\kappa r}{4}(H-ϵ\tilde{p}-\sqrt{p^2+m^2})}·\frac{1}{\sqrt{p^2+m^2}}\right\}\mathrm{sgn}\left(r\right),\hfill \end{array}$$

(96)

and can be solved numerically. However superficial coordinate singularities appear where the denominator $\left\{2-\kappa r(H-ϵ\tilde{p}-\sqrt{p^2+m^2})/4\right\}$ vanishes, corresponding to $W\left(x\right)=-1$. This is the transit point between two branches $W_0$ and $W_{-1}$. and are a consequence of t being a coordinate time.

This problem can be dealt with by describing the trajectories in terms of the proper time ${\tau }_a$ of each particle. Using (89) this is

$$\begin{array}{ccc}\hfill d{\tau }_a^2& =& d{t}^2\left\{N_0{\left(z_a\right)}^2-{(N_1\left(z_a\right)+{\dot{z}}_a)}^2\right\}=d{t}^2{N}_0{\left(z_a\right)}^2\frac{m_a^2}{p_a^2+m_a^2}(a=1,2)\hfill \end{array}$$

(97)

becoming

$$d\tau =d{\tau }_1=d{\tau }_2=\frac{(H-2ϵ\tilde{p})m}{\left\{2-\frac{\kappa r}{4}(H-2ϵ\tilde{p})\right\}(\sqrt{p^2+m^2}-ϵ\tilde{p})\sqrt{p^2+m^2}}dt$$

(98)

for both particles in the equal mass case. The canonical Equations (95) and (96) then become

$$\begin{array}{ccc}\hfill \frac{dp}{d\tau }& =& -\frac{\kappa }{4m}\sqrt{p^2+m^2}\left\{H(\sqrt{p^2+m^2}-ϵ\tilde{p})-m^2\right\}\mathrm{sgn}\left(r\right)\hfill \end{array}$$

(99)

$$\begin{array}{ccc}\hfill \frac{dr}{d\tau }& =& \frac{2ϵ}{m}(\sqrt{p^2+m^2}-ϵ\tilde{p})\left\{\frac{\left[2-\frac{\kappa r}{4}(H-ϵ\tilde{p}-\sqrt{p^2+m^2})\right]\sqrt{p^2+m^2}}{H-2ϵ\tilde{p}}-1\right\}\mathrm{sgn}\left(r\right)\hfill \end{array}$$

(100)

which remarkably have an exact solution.

This is procured as follows. Noting that H is constant, (99) is an ordinary differential equation that can be solved for $p\left(\tau \right)$. The trajectory $r\left(\tau \right)$ can then be obtained by solving (77) or by directly solving (100) after substituting the solution for p, yielding an exact expression for the proper separation of the two bodies as a function of their mutual proper time. The result is [107]

$$\begin{array}{ccc}\hfill p\left(\tau \right)& =& \frac{ϵm}{2}\left(f_0\left(\tau \right)-\frac{1}{f_0\left(\tau \right)}\right)\mathrm{sgn}\left(r\right)\hfill \end{array}$$

(101)

$$\begin{array}{ccc}\hfill r\left(\tau \right)& =& \frac{16\mathrm{sgn}\left(r\right)}{\kappa \left\{H-m\left(f_0-\frac{1}{f_0}\right)\right\}}{\tanh }^{-1}\left[\frac{H-m\left(f_0+\frac{1}{f_0}\right)}{H-m\left(f_0-\frac{1}{f_0}\right)}\right],\hfill \end{array}$$

(102)

with

$$f_0\left(\tau \right)=\frac{H}{m}\left[1-\frac{\sqrt{p_0^2+m^2}-ϵp_0\mathrm{sgn}\left(r\right)-\frac{m^2}{H}}{\sqrt{p_0^2+m^2}-ϵp_0\mathrm{sgn}\left(r\right)}{e}^{\frac{ϵ\kappa m}{4}(\tau -{\tau }_0)}\right],$$

(103)

where $p_0$ is the initial momentum at $\tau ={\tau }_0$.

In Figure 2 the trajectory r is plotted for various values of $H_0$ as a function of the mutual proper time of the two equal mass particles, with the corresponding phase space trajectories plotted underneath. For small $H_0$ we see the distorted oval-shaped trajectory similar to that in Figure 1. The corresponding motion (shown as the solid curve in the upper diagram) is oscillatory, corresponding to the two particles flying apart to some maximal separation and then merging together to repeat the motion with their positions interchanged. As $H_0$ increases we observe considerable distortion in the phase-space trajectory. The motion becomes increasingly asymmetric over a given half-period, with the maximal separation occurring at relatively smaller values of $\tau$. For most of the motion $p>0$; once $p<0$ the particles rapidly merge together to then repeat the motion with their positions switched.

As the energy $H_0$ of the 2-body system increases, the departure from non-relativistic motion is quite striking. This is illustrated in Figure 3 for $H_0=25m$. We see that the period of the relativistic motion is only slightly longer than a half-period of its non-relativistic counterpart. Furthermore, the relativistic system is much more tightly bound, with the maximal separation approximately half that of the non-relativistic case.

5.4. Exact Two-Body Motion with Equal Masses and Arbitrary Charges

Turning now to the charged system, in the center of inertia frame $p_1=-p_2=p$ the determining Equations (78) and (81) become

$$({\mathcal{J}}_{\Lambda }^2+B^2)\tanh \left(\frac{\kappa }{8}{\mathcal{J}}_{\Lambda }|r|\right)=2{\mathcal{J}}_{\Lambda }B,$$

(104)

and

$$({\tilde{\mathcal{J}}}_{\Lambda }^2-B^2)\tan \left(\frac{\kappa }{8}{\tilde{\mathcal{J}}}_{\Lambda }|r|\right)=-2{\tilde{\mathcal{J}}}_{\Lambda }B,$$

(105)

respectively, where

$$\begin{array}{ccc}& & {\mathcal{J}}_{\Lambda }=\sqrt{{\left(\sqrt{H^2+\frac{8{\Lambda }_e}{{\kappa }^2}}-2ϵ\tilde{p}\right)}^2-\frac{8e_1{e}_2}{\kappa }-\frac{8{\Lambda }_e}{{\kappa }^2}}\hfill \\ & & {\tilde{\mathcal{J}}}_{\Lambda }=\sqrt{\frac{8e_1{e}_2}{\kappa }+\frac{8{\Lambda }_e}{{\kappa }^2}-{\left(\sqrt{H^2+\frac{8{\Lambda }_e}{{\kappa }^2}}-2ϵ\tilde{p}\right)}^2}\hfill \\ & & B=H-2\sqrt{p^2+m^2}.\hfill \end{array}$$

(106)

Each case further subdivides into two parts:

$$\begin{array}{ccc}\hfill \tanh \left(\frac{\kappa }{16}{\mathcal{J}}_{\Lambda }|r|\right)& =& \frac{B}{{\mathcal{J}}_{\Lambda }}\left(\tanh \text{-}\mathrm{type}\mathrm{A}\right)\hfill \end{array}$$

(107)

$$\begin{array}{ccc}\hfill \tanh \left(\frac{\kappa }{16}{\mathcal{J}}_{\Lambda }|r|\right)& =& \frac{{\mathcal{J}}_{\Lambda }}{B}\left(\tanh \text{-}\mathrm{type}\mathrm{B}\right)\hfill \end{array}$$

(108)

for (104) and

$$\begin{array}{ccc}\hfill \tan \left(\frac{\kappa }{16}{\tilde{\mathcal{J}}}_{\Lambda }|r|\right)& =& -\frac{B}{{\tilde{\mathcal{J}}}_{\Lambda }}\left(\tan \text{-}\mathrm{type}\mathrm{A}\right)\hfill \end{array}$$

(109)

$$\begin{array}{ccc}\hfill \tan \left(\frac{\kappa }{16}{\tilde{\mathcal{J}}}_{\Lambda }|r|\right)& =& \frac{{\tilde{\mathcal{J}}}_{\Lambda }}{B}\left(\tan \text{-}\mathrm{type}\mathrm{B}\right)\hfill \end{array}$$

(110)

for (105). If both ${\Lambda }_e=0$ and $e_a=0$ then (108) has no solutions since ${\mathcal{J}}_{\Lambda }/B$ exceeds unity.

The canonical equations of motion

$$\begin{array}{ccc}\hfill \dot{p}& =& -\frac{\kappa {\mathcal{J}}_{\Lambda }^2({\mathcal{J}}_{\Lambda }^2-B^2)}{16\mathcal{C}}\mathrm{sgn}\left(r\right),\hfill \end{array}$$

(111)

$$\begin{array}{ccc}\hfill \dot{r}& =& 2ϵ\sqrt{1+\frac{8{\Lambda }_e}{{\kappa }^2{H}^2}}\left(1-\frac{{\mathcal{J}}_{\Lambda }^2}{\mathcal{C}}\right)\mathrm{sgn}\left(r\right)+\frac{2{\mathcal{J}}_{\Lambda }^2}{\mathcal{C}}\frac{p}{\sqrt{p^2+m^2}}.\hfill \end{array}$$

(112)

where

$$\mathcal{C}=\frac{1}{\sqrt{1+\frac{8{\Lambda }_e}{{\kappa }^2{H}^2}}}\left\{\sqrt{1+\frac{8{\Lambda }_e}{{\kappa }^2{H}^2}}{\mathcal{J}}_{\Lambda }^2-\left(\sqrt{H^2+\frac{8{\Lambda }_e}{{\kappa }^2}}-2ϵ\tilde{p}\right)\left(B+\frac{\kappa }{16}({\mathcal{J}}_{\Lambda }^2-B^2)r\right)\right\}$$

(113)

are the same for all four types. Coordinate singularities are present here as before, and these can again be dealt with by choosing the time coordinate to be the common proper time (97), which in this case is

$$d\tau =d{\tau }_1=d{\tau }_2=\frac{m}{\sqrt{p^2+m^2}}\frac{{\mathcal{J}}_{\Lambda }^2}{\mathcal{C}}dt$$

(114)

yielding

$$\begin{array}{ccc}& & \frac{dp}{d\tau }=-\frac{\kappa \sqrt{p^2+m^2}({\mathcal{J}}_{\Lambda }^2-B^2)}{16m}\mathrm{sgn}\left(r\right)\hfill \end{array}$$

(115)

$$\begin{array}{ccc}& & \frac{dr}{d\tau }=\frac{2ϵ}{m}\left\{\frac{\sqrt{p^2+m^2}\mathcal{C}}{{\mathcal{J}}_{\Lambda }^2}-(\sqrt{p^2+m^2}-ϵ\tilde{p})\right\}\mathrm{sgn}\left(r\right)\hfill \end{array}$$

(116)

for the canonical equations of motion (111) and (112).

Solving (115) gives a solution of the form (101), but where

$$f\left(\tau \right)=\left\{\begin{array}{cc}\frac{\frac{H}{m}\left(1+\sqrt{{\gamma }_H}\right)\left\{1-\eta e^{\frac{ϵ\kappa m}{4}\sqrt{{\gamma }_m}(\tau -{\tau }_0)}\right\}}{{\gamma }_e+\sqrt{{\gamma }_m}+\left(\sqrt{{\gamma }_m}-{\gamma }_e\right)\eta e^{\frac{ϵ\kappa m}{4}\sqrt{{\gamma }_m}(\tau -{\tau }_0)}}\hfill & {\gamma }_m>0\hfill \\ \hfill & \\ \frac{1+\sqrt{{\gamma }_H}}{\frac{m}{H}{\gamma }_e+\frac{\sigma }{m-\sigma \frac{ϵ\kappa H}{8}(\tau -{\tau }_0)}}\hfill & {\gamma }_m=0\hfill \\ \hfill & \hfill \\ \frac{\frac{H}{m}(1+\sqrt{{\gamma }_H})}{{\gamma }_e+\sqrt{-{\gamma }_m}\frac{\sigma +\frac{m^2}{H}\sqrt{-{\gamma }_m}tan\left[\frac{ϵ\kappa m}{8}\sqrt{-{\gamma }_m}(\tau -{\tau }_0)\right]}{\frac{m^2}{H}\sqrt{-{\gamma }_m}-\sigma tan\left[\frac{ϵ\kappa m}{8}\sqrt{-{\gamma }_m}(\tau -{\tau }_0)\right]}}\hfill & {\gamma }_m<0\hfill \end{array}\right.$$

(117)

with

$$\begin{array}{cc}\hfill & {\gamma }_H=1+\frac{8{\Lambda }_e}{{\kappa }^2{H}^2}{\gamma }_e=1+\frac{2e_1{e}_2}{\kappa m^2}{\gamma }_m={\gamma }_e^2+\frac{8{\Lambda }_e}{{\kappa }^2{m}^2}\hfill \end{array}$$

(118)

$$\begin{array}{cc}\hfill & \sigma =(1+\sqrt{{\gamma }_H})(\sqrt{p_0^2+m^2}-ϵp_0\mathrm{sgn}\left(r\right))-\frac{m^2}{H}{\gamma }_e\eta =\frac{\sigma -\frac{m^2}{H}\sqrt{{\gamma }_m}}{\sigma +\frac{m^2}{H}\sqrt{{\gamma }_m}}\hfill \end{array}$$

(119)

and $p_0$ is the initial momentum at $\tau ={\tau }_0$. These solutions are valid provided

$$1+\frac{8{\Lambda }_e}{{\kappa }^2{H}^2}\ge 0$$

(120)

which is satisfied for all ${\Lambda }_e>0$, and for ${\Lambda }_e<0$ imposes the constraint

$$H\ge \sqrt{-\frac{8{\Lambda }_e}{{\kappa }^2}}$$

(121)

on the Hamiltonian.

The solution for $r\left(\tau \right)$ is then

$$r\left(\tau \right)=\left\{\begin{array}{cc}\frac{16\mathrm{sgn}\left(r\right){\tanh }^{-1}\left[\frac{\kappa \left(H-m\left|f\left(\tau \right)+\frac{1}{f\left(\tau \right)}\right|\right)}{\sqrt{{\left(\sqrt{{\kappa }^2{H}^2+8{\Lambda }_e}-m\kappa (f\left(\tau \right)-\frac{1}{f\left(\tau \right)})\right)}^2-8\kappa e_1{e}_2-8{\Lambda }_e}}\right]}{\sqrt{{\left(\sqrt{{\kappa }^2{H}^2+8{\Lambda }_e}-m\kappa (f\left(\tau \right)-\frac{1}{f\left(\tau \right)})\right)}^2-8\kappa e_1{e}_2-8{\Lambda }_e}}\hfill & \left(\tanh \text{-}\mathrm{type}\mathrm{A}\right)\hfill \\ \hfill & \\ \frac{16\mathrm{sgn}\left(r\right){\tanh }^{-1}\left[\frac{\sqrt{{\left(\sqrt{{\kappa }^2{H}^2+8{\Lambda }_e}-m\kappa (f\left(\tau \right)-\frac{1}{f\left(\tau \right)})\right)}^2-8\kappa e_1{e}_2-8{\Lambda }_e}}{\kappa \left(H-m\left|f\left(\tau \right)+\frac{1}{f\left(\tau \right)}\right|\right)}\right]}{\sqrt{{\left(\sqrt{{\kappa }^2{H}^2+8{\Lambda }_e}-m\kappa (f\left(\tau \right)-\frac{1}{f\left(\tau \right)})\right)}^2-8\kappa e_1{e}_2-8{\Lambda }_e}}\hfill & \left(\tanh \text{-}\mathrm{type}\mathrm{B}\right)\hfill \end{array}\right.$$

(122)

for the tanh-type solutions, and

$$r\left(\tau \right)=\left\{\begin{array}{c}\frac{16\mathrm{sgn}\left(r\right)\left({\tan }^{-1}\left[\frac{\kappa \left(m\left|f\left(\tau \right)+\frac{1}{f\left(\tau \right)}\right|-H\right)}{\sqrt{8{\Lambda }_e+8\kappa e_1{e}_2-{\left(\sqrt{{\kappa }^2{H}^2+8{\Lambda }_e}-m\kappa (f\left(\tau \right)-\frac{1}{f\left(\tau \right)})\right)}^2}}\right]+n\pi \right)}{\sqrt{8{\Lambda }_e+8\kappa e_1{e}_2-{\left(\sqrt{{\kappa }^2{H}^2+8{\Lambda }_e}-m\kappa (f\left(\tau \right)-\frac{1}{f\left(\tau \right)})\right)}^2}}\left(\tan \text{-}\mathrm{type}\mathrm{A}\right)\hfill \\ \hfill & \\ \frac{16\mathrm{sgn}\left(r\right)\left({\tan }^{-1}\left[\frac{\sqrt{8{\Lambda }_e+8\kappa e_1{e}_2-{\left(\sqrt{{\kappa }^2{H}^2+8{\Lambda }_e}-m\kappa (f\left(\tau \right)-\frac{1}{f\left(\tau \right)})\right)}^2}}{\kappa \left(H-m\left|f\left(\tau \right)+\frac{1}{f\left(\tau \right)}\right|\right)}\right]+n\pi \right)}{\sqrt{8{\Lambda }_e+8\kappa e_1{e}_2-{\left(\sqrt{{\kappa }^2{H}^2+8{\Lambda }_e}-m\kappa (f\left(\tau \right)-\frac{1}{f\left(\tau \right)})\right)}^2}}\left(\tan \text{-}\mathrm{type}\mathrm{B}\right)\hfill \end{array}\right.$$

(123)

for the tan-type solutions.

5.4.1. Neutral Particle Motion

For electrically neutral particles, with $e_1=e_2=0$ the quantity ${\Lambda }_e=-2\Lambda$. In this case the gravitational attraction of the masses competes with cosmic expansion or attraction, depending on the sign of ${\Lambda }_e$.

It is instructive to examine the motion for a range of masses. To this end, an arbitrary fiducial mass $m_0$ can be chosen, which sets a calibration scale for all other quantities. Typical scenarios are illustrated in Figure 4 and Figure 5.

For ${\Lambda }_e<0$ the curvature $R<0$ in the absence of stress-energy (from (14)) corresponding to cosmic attraction, depicted in Figure 4. The particles always remain bound, and for large values of m undergo oscillation about their centre of inertia in a manner not too different from that shown in Figure 2 for small $H_0$. However for very small m a rather surprising situation arises: a sharp extremum develops early in each half-period followed by a second, broader extremum. The first extremum is visible in Figure 4 for the $m=0.05m_0$ curve, but also appears for smaller values of m; it is not visible for $m=0.001m_0$ due to plotting precision.

In all cases in Figure 4 the two particles start at $r=0$ with equal opposing momentum and and depart in opposite directions. For small enough m they quickly reach a maximum separation and then go back toward one another for a certain period of proper time. At some point they each reverse direction, moving apart more slowly until they attain a second maximal separation. They then reverse direction again, returning to their starting point where the motion then repeats itself with the particles interchanged.

This peculiar behavior, first observed in [102,103], is due to a subtle interplay between the gravitational attraction and relativistic motion of the particles in a spacetime with cosmic attraction. It is demonstrative of unexpected relativistic behaviour that self-gravitating systems can exhibit.

If ${\Lambda }_e>0$, $R>0$ in the absence of stress-energy, and a de Sitter-like cosmic expansion ensues. The situation for various values of m is illustrated in Figure 5. For large enough m the motion is bounded. As m decreases, the period of the motion becomes larger, until at a critical value of $m\lesssim 4.729$ the period becomes infinite. Here the cosmic expansion is balances the gravitational attraction of the particles, and for smaller values of m dominates. The bound behaviour is fully analogous to that of an overdense universe that expands to some maximal value of the scale factor and then recontracts; the unbound behaviour is analogous to that of an underdense universe whose deceleration parameter is insufficient to prevent full cosmic dilution.

All of the above behaviour is described by the tanh-type solutions. However there is a countably infinite set of unbounded motions described by the tan-type solutions. Some examples are shown in Figure 6 for the phase-space. For $0<{\Lambda }_e<{\Lambda }_c$ both bounded and unbounded motions exist for fixed value of $H=H_0$, shown in the left panel of Figure 6. The bounded motion consists of the distorted oscillations noted previously, whereas the unbounded motion consists of two particles simply approaching each other at some minimal value of $|r|$ after which they reverse direction and proceed toward infinity. The dotted curves come from the upper expression in (123) and the dashed curves from the lower one. As ${\Lambda }_e\to {\Lambda }_c$ the bulges of the solid tanh-type A curve and the dotted tan-type A curve come close, making contact at the critical value ${\Lambda }_e={\Lambda }_c$. For ${\Lambda }_e>{\Lambda }_c$ these two curves bifurcate into the solid curves shown in the right panel of Figure 6: the particles cross one another before receding toward infinity. The upper solid curve represents the motion in which

$$p\to p_{\pm }=\frac{\left(\pm \sqrt{{\kappa }^2{H}^2+8{\Lambda }_e}+\sqrt{8{\Lambda }_e}\right)}{2\kappa }$$

(124)

as $r\to \pm \infty$.

One peculiar feature of the unbounded motion is that the two particles diverge to infinite separation at finite proper time. The trajectory $x\left(t\right)$ of a null ray emitted from particle 2 at time T is governed by $d\tau =0$, yielding

$$\frac{dx}{dt}=\pm N_0(x\left(t\right),z_1\left(t\right),z_2\left(t\right),p\left(t\right))-N_1(x\left(t\right),z_1\left(t\right),z_2\left(t\right),p\left(t\right))$$

(125)

for the equation of the light signal directed toward (+) or away from (-) particle 1. Numerically solving (125) yields the trajectories shown in Figure 7, which plots the trajectories of null rays emitted from particle 2 at various times T. For small $\kappa m_{0}T\lesssim 3.4)$, the particles are in causal contact, as shown by the dotted curve of positive slope in the left panel of Figure 7. At $\kappa m_{0}T\approx 3.4$ the light ray just barely catches up to particle 1, which itself is almost in light-like motion; this is shown by the dashed curve of positive slope. Beyond this time the particles are no longer in causal contact. For larger T, the null ray of positive slope quickly asymptotes to a line parallel to the asymptotic worldline of particle 1, shown by the dashed curve. The light rays emitted away from particle 1 all asymptote to curves parallel to the trajectory of particle 2. As T becomes even larger, the null ray emitted in the direction of particle 1 experiences a repulsive effect, ultimately reversing direction as shown in the right panel of Figure 7. For $\kappa m_{0}T\approx 4.82$ the null ray emitted toward particle 1 asymptotes to $r=0$. For larger values of T null rays emitted in the positive direction eventually reverse their motion and asymptote to curves that are parallel to the asymptotic trajectory of particle 2, shown by the dashed and dot-dashed lines in the left panel of Figure 7. The trajectories of the light signals emitted away from particle 1 cannot be distinguished in the plot from the trajectory of particle 2.

This behaviour can be captured in a flat-space model as follows. Consider the 2-velocity

$$u^{\mu }=(h\left(\sigma \tau \right),\sqrt{h^2\left(\sigma \tau \right)-1})$$

(126)

where $h\left(\sigma \tau \right)$ is some function and

$$d{s}^2=-d{t}^2+d{x}^2$$

(127)

is the flat metric. The general expression for the 2- acceleration is

$$a^{\mu }=\frac{d{u}^{\mu }}{d\tau }=\sigma h^{\prime }(1,\frac{h}{\sqrt{f^2\left(\sigma \tau \right)-1}})$$

(128)

whose magnitude is

$$a·a=\frac{{\left(\sigma h^{\prime }\right)}^2}{h^2\left(\sigma \tau \right)-1}$$

(129)

where $h^{\prime }=dh\left(\tau \right)/d\tau$. Noting that $u·u=1$ and $a·u=0$ we have the following possible scenarios:

• As $\tau \to \infty$, $h\to h_0$ where $h_0$ is finite. In this case the particle never becomes lightlike.
• As $\tau \to \infty$, $h\to \infty$. In this case the particle becomes lightlike, but this happens in an infinite amount of proper time (and coordinate time). The standard example is $h=cosh\left(\sigma \tau \right)$, the constant acceleration example.
• The function $h\to \infty$ as $\tau \to {\tau }_0$, where ${\tau }_0$ is finite. In this case the particle asymptotes to a lightlike trajectory in a finite amount of proper time, but an infinite amount of coordinate time; an example is $h=sec\left(\sigma \tau \right)$. The acceleration increases as a function of proper time, diverging at $\tau ={\tau }_0$. This last situation is realized by the exact solutions (123) with $\Lambda >{\Lambda }_c$.

The proper time ${\tau }_{\infty }$ at which the particles attain infinite separation $r⟶\infty$ is

$${\tau }_{\infty }=\frac{4}{\kappa m\sqrt{{\gamma }_m}}\log \left(\frac{H_0(1+\sqrt{{\gamma }_H})-(p_{+}+\sqrt{p_{+}^2+m^2})(1+\sqrt{{\gamma }_m})}{\eta \left[H_0(1+\sqrt{{\gamma }_H})-(p_{+}+\sqrt{p_{+}^2+m^2})(\sqrt{{\gamma }_m}-1)\right]}\right)$$

(130)

for any given m, where $p_{+}$ is given by (124).

5.4.2. Charged Particle Motion with ${\Lambda }_e=0$

There are two qualtiatively distinct cases of charged particle motion: attractive, when the particles have opposite sign charges, and repulsive, when the charges have the same sign. The situation becomes further complicated for nonzero ${\Lambda }_e$, and so it is more instructive to first consider ${\Lambda }_e=0$. Since the charges appear in all solutions as the product $e_1{e}_2$, it is sufficient to set $e_1=-e_2=q$ for the attractive case and $e_1=e_2=q$ for the repulsive case.

For ${\Lambda }_e=0$ (117) simplifies to

$$\begin{array}{ccc}\hfill f_e\left(\tau \right)& =& \frac{H}{m{\gamma }_e}\left\{1-{\eta }_e{e}^{\frac{ϵ\kappa m}{4}{\gamma }_e(\tau -{\tau }_0)}\right\}\hfill \end{array}$$

(131)

where

$${\eta }_e=\frac{\sqrt{p_0^2+m^2}-ϵp_0\mathrm{sgn}\left(r\right)-\frac{m^2}{H}{\gamma }_e}{\sqrt{p_0^2+m^2}-ϵp_0\mathrm{sgn}\left(r\right)}$$

from (118) and (119) and $p\left(\tau \right)$ has the form (101) with $f_0\to f_e$; the relative distance $r\left(\tau \right)$ is likewise obtained from (107) - (110).

5.4.2.1 The attractive case: $e_1=-e_2=q$

The quantity ${\kappa }^2{\left(H-mf\left(\tau \right)+\frac{m}{f\left(\tau \right)}\right)}^2-8\kappa e_1{e}_2>0$ in the attractive case and the two particles always remain bound, with the motion described by the tanh-type solution (107). As expected, as $|q|$ increases the particles are more tightly bound. This is clear from Figure 8, where the maximal separation of the particles decreases as $|q|$ increases, due to the additional electromagnetic attraction. The frequency of the motion correspondingly increases.

The period is determined from the initial value of $p_0=\sqrt{{(H/2)}^2-m^2}$ at $r=0$:

$$T=\frac{16}{\kappa m{\gamma }_e}{\tanh }^{-1}\left(\frac{{\gamma }_e\sqrt{H^2-4m^2}}{(2-{\gamma }_e)H}\right).$$

(132)

Although the above expression diverges when ${\gamma }_e=\frac{2H}{2p_0+H}$, this situation is never realized in the attractive case, since ${\gamma }_e<1$ whereas $\frac{2H}{2p_0+H}>1$.

It is instructive to compare the exact relativistic motion to the motion in three approximations :

• the non-relativistic motion described by the Hamiltonian

  $$H=2m+\frac{p^2}{m}+\frac{q^2}{2}|r|+\frac{\kappa m^2}{4}|r|,$$

  (133)
• the linear approximation in $\kappa$, whose Hamiltonian is

  $$\begin{array}{ccc}\hfill H& =& 2\sqrt{p^2+m^2}+\frac{1}{2}{q}^2|r|+\frac{\kappa }{4}\left\{(2p^2+m^2-2ϵ\tilde{p}\sqrt{p^2+m^2})|r|\right.\hfill \\ & & \left.+\frac{1}{2}{q}^2(\sqrt{p^2+m^2}-ϵ\tilde{p})r^2+\frac{1}{24}{q}^4{|r|}^3\right\},\hfill \end{array}$$

  (134)
• the $\kappa \to 0$ limit, which is special-relativistic electrodynamics in (1+1)-dimensional flat space-time ; its Hamiltonian is

  $$H=2\sqrt{p^2+m^2}+\frac{q^2}{2}|r|.$$

  (135)

This comparison is illustrated in Figure 9 for $H_0=3m_0$ for $q/m_0=0.5,1,5$ and 10. For small $q=0.5$ both the exact solution (solid curve) and the linear approximation (dashed curve) follow S-shaped trajectories, whereas the non-relativistic (dotted curve) and flat electrodynamics (dot-dashed curve) trajectories have symmetric oval shapes. As $|q|$ the effect of gravity becomes relatively weak and all trajectories tend to coincide with the trajectory of flat electrodynamics.

5.4.2.2 The repulsive case: $e_1=e_2=q$, $|q|\le q_c$

In this case the electromagnetic force is repulsive and competes with the attractive gravitational force. There is a critical value

$$q_c=\sqrt{\frac{\kappa }{8}}\left(H-\sqrt{H^2-4m^2}\right),$$

(136)

of the charge, obtained from setting ${\mathcal{J}}^2={(H-2ϵ\tilde{p})}^2-8q^2/\kappa =0$. If $|q|<q_c$, bounded and unbounded motions can occur, described respectively by the tanh-type (${\mathcal{J}}^2>0$) and tan-type solutions (${\mathcal{J}}^2<0$) for $r\left(\tau \right)$. Alternatively, (136) gives the critical value of H for fixed $\kappa$ and q, or the critical value of $\kappa$ for fixed H and q, both corresponding to the transition from bounded to unbounded motion.

The two possibilities are illustrated in Figure 10, which plots $r\left(\tau \right)$ plots for fixed $H_0=3m_0$ (left) or fixed $|q|/m_0=0.1$ (right). We see for fixed $H_0$ that the period becomes larger as $|q|$ increases; once the critical value $q_c/m_0=0.2700907567$ is exceeded the motion becomes unbounded and the separation of the two particles diverges at finite $\tau$. Likewise, if $|q|$ is fixed, then a similar transition takes place at a threshold value of $H_0=H_c$; for the particular choice given here, $H_c/m_0=7.21249$.

The existence of two types of motion for fixed H and q has no non-relativistic analogue, and is a qualitatively new aspect of the relativistic gravitation. An illustration is provided in Figure 11, showing $r\left(\tau \right)$ for $H_0=3m_0$ and $|q|/m_0=0.25<q_c/m_0$. Both bounded motion and a sequence of possible unbounded motions can exist. The tan-type A and B motions have the specific feature that $r\left(\tau \right)$ diverges (the particles attain infinite separation) at finite proper time (but at infinite coordinate time).

5.4.2.3 The repulsive case: $e_1=e_2=q$, $|q|>q_c$

Only unbounded motion is permitted in this case since the electromagnetic repulsive force overwhelms the attractive gravitational force if $|q|>q_c$. The $r\left(\tau \right)$ trajectories are qualitatively similar to the unbounded motions depicted in Figure 11; to understand the distinctions between unbounded motion for this case and that for the $|q|\le q_c$ case it is instructive to consider the phase space plots.

Consider first the physical region of $(|q|,p)$ parameter space, shown in Figure 12 for $H_0=3m_0$. From (106), the shaded area is the region ${\mathcal{J}}^2>0$ and $B>0$, where tanh-type A and tanh-type B give the actual trajectories. The boundary of this area is fixed by $p=-\sqrt{2/\kappa }|q|+H_0/2$ and $p=\pm p_0=\pm \sqrt{{(H_0/2)}^2-m^2}$ ($=\pm \sqrt{5}{m}_0/2$ for the specific choice in Figure 12). The values of $|q|$ at the intersections of these boundary lines are $q_d$ (for the negative solution) and $q_c$ (the critical value (136)) for the positive one. The tanh-type B solution is realized in a very narrow region between the dashed ($\mathcal{J}-B=0$) and solid ($p=-\sqrt{2/\kappa }|q|+H_0/2$) curves in the shaded region. The tan-type A and B trajectories are realized in the unshaded ${\mathcal{J}}^2<0$ region.

The dotted line in Figure 12 denotes a constant $|q|/m_0$ line whose intercepts with the diagonal lines $p=\pm \sqrt{2/\kappa }|q|+H_0/2$ are $p_1/m_0$ and $p_2/m_0$. In the subcritical region $0<|q|<q_c$ this dotted line would be to the left of $q_c/m_0$. The tanh-type A solutions are realized in the shaded region with $-p_0<p<p_0$, and the associated $r\left(\tau \right)$ trajectory is given by the solid line in Figure 11. In the unshaded region $p_2<p<p_1$ between the diagonal lines and to the left of $q_c/m_0$, both tan-type A and B are realized; the associated $r\left(\tau \right)$ trajectories are respectively given by the dotted and dashed lines in Figure 11. The region $p_0<p<p_2$ is forbidden.

In the supercritical region $q_c<|q|<q_d$ all values of p in the range $-p_0<p<p_1$ are allowed, and the associated phase space trajectories are shown in the left panel of Figure 13. The single bounded motion N in the right panel of Figure 11 merges with the unbounded motion A0 at $|q|=q_c$ to yield two new unbounded trajectories $N1$ (corresponding to the unshaded region $p_2<p<p_1$), and $N2$ (corresponding to the shaded region $-p_0<p<p_2$) in Figure 13. The remaining unbounded motions $A_n$ and $B_n(n=1,2,\dots )$ are all described by respectively by the tan-type A and B solutions. In both Figure 11 and Figure 13 the analogous trajectories in flat space electrodynamics(the dot-dashed curves) are shown to illustrate the strongly deforming effects of gravity on the trajectories.

Finally for $|q|>q_d$ all solutions are tan-type A and B since $p_2<-p_0$. Sample phase space trajectories are shown in the right panel of Figure 13. A characteristic cusp appears at $r=0$ in the $N1$ and $N2$ trajectories. For the B trajectories the motion switches between B0 (for $-p_0\le p\le p_0$) and B1 (for $p_2\le p\le -p_0$ and $p_0\le p\le p_1$); the full motion is denoted B01. $B12$ is likewise composed of a combination of B1 and B2.

Another interesting aspect of the supercritical $|q|>q_c$ case is that $H<2m$ is possible – the total energy can be less than the rest energy of the particles! For $H<2m$ only ${\mathcal{J}}^2<0$ is possible and the shaded region in Figure 12 is absent. Consequently only unbounded motion is possible; some trajectories in phase space for this case are shown in Figure 14. All types of unbounded motions $A_n$ and $B_n$ are realized. Note that all trajectories curve more toward the r-axis (due to the additional effect of gravitational attraction) relative to that of flat space electrodynamics. They are also shifted in the positive p direction due to the p-dependence of the gravitational potential.

In 1+1 flat-space electrodynamics the total energy of the system is restricted to $H\ge 2m$ for attractive charges , but no such restriction on H exists for repulsive charges. We see that in a general relativistic theory the same situation prevails.

The time-reversed trajectories corresponding to these cases can all be obtained by setting $ϵ=-1$, and trajectories in the $r<0$ region are obtained from those in $r>0$ by reversing the signs of both r and p.

5.4.3. Charged Particle Motion with ${\Lambda }_e\ne 0$

This is the most general case for two-body motion in the electrovacuum. The dynamics is governed by a combination of gravitational attraction, electromagnetic attraction/repulsion, and the cosmological constant. The signs of ${\gamma }_m$ and ${\kappa }^2{\left(H_0-mf\left(\tau \right)+\frac{m}{f\left(\tau \right)}\right)}^2-8\kappa e_1{e}_2-8{\Lambda }_e$ characterize the solutions.

The possible range of motions is now rich and complicated [104], and a set of sample trajectories is presented in Figure 15. One of the noteworthy special cases occurs particular range of negative ${\Lambda }_e$ and small mass (or large $H/m_0$), where the trajectories have a double-peaked structure. This is shown in the top panel of Figure 15 and is particularly noticeable for ${\Lambda }_e/{\left(\kappa m_0\right)}^2=-0.5$, shown as the solid curve. The particles begin at $r=0$ with equal and opposite initial momenta, reach a maximum separation, and then reverse direction due to the combined attractive electromagentic and graviational forces. But they soon reverse direction again, moving apart until reaching a second maximum after which they return to the starting point. For small values of $|q|$ the double peak structure is present, but for sufficiently large $|q|$ it disappears.

This peculiar behavior takes place due to the induced momentum-dependent potential in conjunction with the dual force attraction in the cosmological vacuum. An expansion of the Hamiltonian in terms of $\kappa$ and ${\Lambda }_e/{\kappa }^2$

$$\begin{array}{ccc}\hfill H& =& 2\sqrt{p^2+m^2}+\frac{\kappa }{4}{(\sqrt{p^2+m^2}-ϵ\tilde{p})}^2|r|+\frac{{\kappa }^2}{4^2}{(\sqrt{p^2+m^2}-ϵ\tilde{p})}^3{r}^2\hfill \\ & & +\frac{7{\kappa }^3}{6\times 4^3}{(\sqrt{p^2+m^2}-ϵ\tilde{p})}^4{|r|}^3-\frac{{\Lambda }_e}{2\kappa }·\frac{ϵ\tilde{p}}{\sqrt{p^2+m^2}}|r|-\frac{{\Lambda }_e}{16}·\frac{ϵ\tilde{p}{m}^2}{p^2+m^2}{r}^2\hfill \\ & & +\frac{{\Lambda }_e^2}{4{\kappa }^3}·\frac{ϵ\tilde{p}}{{(p^2+m^2)}^{3/2}}|r|+···\hfill \end{array}$$

(137)

illustrates that negative ${\Lambda }_e$ acts effectively as an attractive force leading to bounded (periodic) motions, whereas positive ${\Lambda }_e$ acts effectively as a repulsive force.

For ${\Lambda }_e<0$ there can be both bounded and unbounded motion in the electromagnetically repulsive case depending on the size of $|q|$. For small $|q|$ cosmological attraction is smaller, whereas for large $|q|$ the opposite situation prevails. There is an intermediate range of $|q|$ for a given ${\Lambda }_e$ where both bounded and unbounded motion is possible [104]. Sample trajectories are shown in the middle panel of Figure 15.

In the electromagnetically attractive case, for ${\Lambda }_e>0$ there can again be both bounded and unbounded motion depending on the size of $|q|$. In this case attraction due to a sufficiently large value of $|q|$ will overwhelm cosmological repulsion. These effects depend on the relative size of H as compared to $\sqrt{{\kappa }^2{m}^2/2{\Lambda }_e+4m^2}$; if H is larger than this latter quantity then there is again an intermediate range of $|q|$ where both bounded and unbounded motions exist [104]. Sample trajectories are shown in the bottom panel of Figure 15.

Finally the case of joint electromagnetic and cosmological repulsion yield trajectories similar to those in Figure 11.

5.5. Exact Two-Body Motion with Unequal Masses

For unequal masses the proper time (97) of each particle

$$\begin{array}{cc}& d{\tau }_1=dt\frac{16Y{K}_0{K}_1{m}_1}{JK{M}_1\sqrt{p^2+m_1^2}}d{\tau }_2=dt\frac{16Y{K}_0{K}_2{m}_2}{JK{M}_2\sqrt{p^2+m_2^2}}\end{array}$$

(138)

differs, where $K\equiv K_{+}=K_{-}$ and $Y\equiv Y_{+}=Y_{-}$. While there are a variety of choices of time parameter one could use to describe the motion, the most natural one appears to be [103,104]

$$d\tilde{\tau }\equiv dt\frac{16Y{K}_0}{JK}{\left(\frac{K_1{K}_2{m}_1{m}_2}{M_1{M}_2\sqrt{p^2+m_1^2}\sqrt{p^2+m_1^2}}\right)}^{1/2},$$

(139)

which is symmetric with respect to $1\leftrightarrow 2$ and reduces to the proper time (114) when $m_1=m_2$.

The parameter space now has an additional variable and analysis of the various types of motions is correspondingly more detailed. However the general rubric of the four types of motion depending on cosmological attraction/repulsion and electromagnetic attraction/repulsion continues to hold.

A number of the more interesting sample cases are illustrated in Figure 16 and Figure 17. When $m_1$ is large, $m_2$ behaves like a test body and the two particles oscillate about their centre of inertia for all values of the parameters shown. The gravitational attraction is stronger and both the period and proper maximal separation between the particles becomes shorter.

For ${\Lambda }_e<0$, as the mass ratio decreases the maximal separation and period both get larger as shown in Figure 16. However eventually the attractive effect due to the cosmological constant begins to dominate, and for quite small mass ratios the period and proper maximal separation decrease and the double peak structures emerge. For larger $|q|$ (top panel) the electromagnetic attraction is stronger yielding correspondingly shorter periods and maximal separations as compared to the smaller $|q|$ case (bottom panel).

If there is either cosmological or electromagnetic repulsion (or both) then unbounded motion is possible, as shown in Figure 17. Large $m_1$ yields bounded motion as before, but for a small enough mass ratio the repulsive effect dominates and the particles fly apart to infinity. Maximal separations and periods are larger for cosmological and electromagnetic repulsion, as shown in the bottom panel of Figure 17.

5.6. Static Balance

In $(3+1)$ dimensions exact solutions to the 2-body problem in general relativity have been unobtainable, primarily because gravitational radiation carries energy away from the system. However a condition of static balance – in which gravitational attraction is exactly balanced by electromagnetic repulsion – was attained for two bodies [111,112] and later for N bodies on a line [113]. This condition is

$$e_i=\pm \sqrt{4\pi G}{m}_i,$$

(140)

which is more stringent than the corresponding condition

$$G{m}_1{m}_2-\frac{e_1{e}_2}{4\pi }=0$$

(141)

in Newtonian theory. Whether or not (140) is a is a necessary condition for static balance remains an open question.

The corresponding static balance condition in (1+1) dimensions can be straightforwardly obtained from the determining Equation (77): it is simply the extremum of H with respect to r. Setting $\partial H/\partial r=0$ yields, after some algebra,

$$\frac{\kappa }{2}(\sqrt{p^2+m_1^2}-ϵ\tilde{p})(\sqrt{p^2+m_2^2}-ϵ\tilde{p})-e_1{e}_2=0$$

(142)

which is the force-balance condition. If $p=0$ then

$$\frac{\kappa }{2}{m}_1{m}_2-e_1{e}_2=0$$

(143)

which is the static balance condition [108]. This condition is identical with that of (141) in Newtonian theory in (1+1) dimensions.

However there is a more general solution to (142), which is

$$p_c=\pm \frac{|{\left(\frac{\kappa }{2}\right)}^2{m}_1^2{m}_2^2-e_1^2{e}_2^2|}{\sqrt{2\kappa e_1{e}_2}\sqrt{(\frac{\kappa }{2}{m}_1^2+e_1{e}_2)(\frac{\kappa }{2}{m}_2^2+e_1{e}_2)}}.$$

(144)

provided $e_1{e}_2>0$, where $p=p_c$ is the (constant) value of the momentum for which the repulsive electromagnetic and attractive gravitational forces between both particles are the same.

Non-relativistically (143) is the force-balance condition that includes both the static case and uniform motion. However in the relativistic case these two are different. In general the force-balance condition (142) depends on the momentum and allows for uniform motion in the centre of inertia system in which both particles either approach or recede with constant momentum (144). The static balance condition (143) is the special case for which $p_c=0$.

6. The 3-Body Problem

The 3-body problem for a relativistic self-gravitating system is an interesting subject of study, in large part because its non-relativistic counterpart models several interesting physical systems. These include two elastically colliding billiard balls in a uniform gravitational field [114], perfectly elastic collisions of a particle with a wedge in a uniform gravitational field [3], and a “linear baryon” (a bound state of three quarks along a line) [115]. These systems can even be tested experimentally [116].

The first study of 3-body motion in a fully relativistic context was carried out over 20 years ago [28,29], with extensions to include unequal masses [117], a cosmological constant [118], and charge [119] subsequently undertaken. The most effective means to study the dynamics of the 3-body ROGS is to work in the canonical formalism [98]. This approach yields an exact expression for the Hamiltonian in terms of the four physical degrees of freedom of the system: the two proper separations of the point particles and their conjugate momenta.

The non-relativistic system can be shown to be equivalent to that of a single particle moving in a hexagonal-well potential in 2 spatial dimensions. The 3-body ROGS has an exact relativistic generalization of this hexagonal-well problem, providing insight into intrinsically non-perturbative relativistic effects. For equal mass particles, the cross-sectional shape of the well is that of a regular hexagon in the non-relativistic case. Unequal masses distort this symmetry to that of a ‘squashed’ hexagon, whose sides have differing length. The hexagonal symmetry is maintained relativistically, with the sides of the hexagon curved outward.

Poincaré sections can be used to examine the global structure of the phase space. The relativistic plots are qualitatively similar to their non-relativistic counterparts, but distorted toward the lower-right of the phase plane. This is due relativistic gravitational momentum that continuously distorts the basic structure of the Poincaré plot.

6.1. 3-Body Constraint Equations

The constraints (63) and (64) for a 3-body system become

$$\begin{array}{ccc}\hfill ▵\varphi -\frac{{\kappa }^2}{4}{\left({\chi }^{\prime }\right)}^2\varphi & =& \frac{\kappa }{4}\left\{\sqrt{p_1^2+m_1^2}\varphi \left(z_1\right)\delta (x-z_1)+\sqrt{p_2^2+m_2^2}\varphi \left(z_2\right)\delta (x-z_2)\right.\hfill \end{array}$$

$$\begin{array}{ccc}& & \left.+\sqrt{p_3^2+m_3^2}\varphi \left(z_3\right)\delta (x-z_3)\right\}\hfill \end{array}$$

(145)

$$\begin{array}{ccc}\hfill ▵\chi & =& -\frac{1}{2}\left\{p_1\delta (x-z_1)+p_2\delta (x-z_2)+p_3\delta (x-z_3)\right\}\hfill \end{array}$$

(146)

where the general solution to this latter equation is

$$\chi =-\frac{1}{4}\left\{p_1|x-z_1|+p_2|x-z_2|+p_3|x-z_3|\right\}-ϵXx+ϵC_{\chi }.$$

(147)

As before $ϵ$ (${ϵ}^2=1$) flips sign under time-reversal, and X and $C_{\chi }$ are constants of integration.

The strategy for solving (145) is similar to the 2-body case: spacetime is divided into four regions

$$z_1<x\hspace{1em}\hspace{1em}\hspace{1em}(+)\mathrm{region}$$

$$z_2<x<z_1\hspace{1em}\left(1\right)\mathrm{region}$$

$$z_3<x<z_2\hspace{1em}\left(2\right)\mathrm{region}$$

$$x<z_3\hspace{1em}\hspace{1em}\hspace{1em}(-)\mathrm{region}$$

assuming $z_3<z_2<z_1$. Within each of these regions ${\chi }^{\prime }$ is constant

$${\chi }^{\prime }=\left\{\begin{array}{cc}-ϵX-\frac{1}{4}(p_1+p_2+p_3)\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}(+)\mathrm{region}\hfill \\ -ϵX+\frac{1}{4}(p_1-p_2-p_3)\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left(1\right)\mathrm{region}\hfill \\ -ϵX+\frac{1}{4}(p_1+p_2-p_3)\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left(2\right)\mathrm{region}\hfill \\ -ϵX+\frac{1}{4}(p_1+p_2+p_3)\hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}(-)\mathrm{region}\hfill \end{array}\right.$$

(148)

and (145) has the solution

$$\left\{\begin{array}{c}{\varphi }_{+}\left(x\right)=A_{+}{e}^{\frac{\kappa }{2}{K}_{+}x}+B_{+}{e}^{-\frac{\kappa }{2}{K}_{+}x}\hfill \\ {\varphi }_1\left(x\right)=A_1{e}^{\frac{\kappa }{2}{K}_{1}x}+B_1{e}^{-\frac{\kappa }{2}{K}_{1}x}\hfill \\ {\varphi }_2\left(x\right)=A_2{e}^{\frac{\kappa }{2}{K}_{2}x}+B_2{e}^{-\frac{\kappa }{2}{K}_{2}x}\hfill \\ {\varphi }_{-}\left(x\right)=A_{-}{e}^{\frac{\kappa }{2}{K}_{-}x}+B_{-}{e}^{-\frac{\kappa }{2}{K}_{-}x}\hfill \end{array}\right.$$

(149)

in each region, where

$$\begin{array}{ccc}& & K_{+}\equiv X+\frac{ϵ}{4}(p_1+p_2+p_3)K_{-}\equiv X-\frac{ϵ}{4}(p_1+p_2+p_3)\hfill \\ & & K_1\equiv X-\frac{ϵ}{4}(p_1-p_2-p_3)K_2\equiv X-\frac{ϵ}{4}(p_1+p_2-p_3)\hfill \end{array}$$

(150)

In order that a full solution to (145) respect boundary conditions that ensure finiteness of the Hamiltonian, the function $\varphi$ must satisfy the following matching conditions

$$\begin{array}{ccc}& & {\varphi }_{+}\left(z_1\right)={\varphi }_1\left(z_1\right){\varphi }_1\left(z_2\right)={\varphi }_2\left(z_2\right){\varphi }_{-}\left(z_3\right)={\varphi }_2\left(z_3\right)\hfill \end{array}$$

(151a)

$$\begin{array}{ccc}& & {\varphi }_{+}^{\prime }\left(z_1\right)-{\varphi }_1^{\prime }\left(z_1\right)=\frac{\kappa }{4}\sqrt{p_1^2+m_1^2}\varphi \left(z_1\right)\hfill \end{array}$$

(151b)

$$\begin{array}{ccc}& & {\varphi }_1^{\prime }\left(z_2\right)-{\varphi }_2^{\prime }\left(z_2\right)=\frac{\kappa }{4}\sqrt{p_2^2+m_2^2}\varphi \left(z_2\right)\hfill \end{array}$$

(151c)

$$\begin{array}{ccc}& & {\varphi }_2^{\prime }\left(z_3\right)-{\varphi }_{-}^{\prime }\left(z_3\right)=\frac{\kappa }{4}\sqrt{p_3^2+m_3^2}\varphi \left(z_3\right).\hfill \end{array}$$

(151d)

at the locations of the particles $x=z_1,z_2,z_3$. Satisfying these conditions is a straightforward but tedious exercise, yielding the result

$$\begin{array}{ccc}\hfill L_1{L}_2{L}_3& =& {\mathfrak{M}}_{12}{\mathfrak{M}}_{21}{L}_3^{*}{e}^{\frac{\kappa }{4}{s}_{12}[(L_1+{\mathfrak{M}}_{12})z_{13}-(L_2+{\mathfrak{M}}_{21})z_{23}]}\hfill \\ & & +{\mathfrak{M}}_{23}{\mathfrak{M}}_{32}{L}_1^{*}{e}^{\frac{\kappa }{4}{s}_{23}[(L_2+{\mathfrak{M}}_{23})z_{21}-(L_3+{\mathfrak{M}}_{32})z_{31}]}\hfill \\ & & +{\mathfrak{M}}_{31}{\mathfrak{M}}_{13}{L}_2^{*}{e}^{\frac{\kappa }{4}{s}_{31}[(L_3+{\mathfrak{M}}_{31})z_{32}-(L_1+{\mathfrak{M}}_{13})z_{12}]}\hfill \end{array}$$

(152)

or more compactly

$$L_1{L}_2{L}_3=\sum _{ijk}\left|{ϵ}^{ijk}\right|{\mathfrak{M}}_{ij}{\mathfrak{M}}_{ji}{L}_k^{*}{e}^{\frac{\kappa }{4}{s}_{ij}[(L_i+{\mathfrak{M}}_{ij})z_{ik}-(L_j+{\mathfrak{M}}_{ji})z_{jk}]}$$

(153)

for the full determining equation of the Hamiltonian, where

$$\begin{array}{ccc}\hfill {\mathfrak{M}}_{ij}& =& M_i-ϵp_i{s}_{ij},M_i=\sqrt{p_i^2+m_i^2}\hfill \end{array}$$

(154)

$$\begin{array}{ccc}\hfill L_i& =& H-M_i-ϵ(\sum _j{p}_j{s}_{ji})L_i^{*}=(1-\prod _{j<k\ne i}{s}_{ij}{s}_{ik})M_i+L_i\hfill \end{array}$$

(155)

with $z_{ij}=(z_i-z_j)$ , $s_{ij}=$sgn$\left(z_{ij}\right)$, and ${ϵ}^{ijk}$ the 3-dimensional Levi-Civita tensor. Equation (153) reproduces the correct determining equation for any permutation of the particles.

The components of the metric can be determined from (49) – (54). Of greater relevance are the canonical equations of motion (30), which become

$$\begin{array}{ccc}& & {\dot{z}}_1\left\{L_2{L}_3+L_1{L}_3+L_1{L}_2\right.\hfill \\ & & -[M_2-ϵp_2{s}_{21}][M_1-ϵp_1{s}_{12}][1+\frac{\kappa }{4}{L}_3^{*}|z_{12}|]e^{\frac{\kappa }{4}{s}_{12}[(L_1+{\mathfrak{M}}_{12})z_{13}-(L_2+{\mathfrak{M}}_{21})z_{23}]}\hfill \\ & & -[M_3-ϵp_3{s}_{31}][M_1-ϵp_1{s}_{13}][1+\frac{\kappa }{4}{L}_2^{*}|z_{13}|]e^{\frac{\kappa }{4}{s}_{13}[(L_1+{\mathfrak{M}}_{23})z_{12}+(L_3+{\mathfrak{M}}_{32})z_{23}]}\hfill \\ & & \left.-[M_2-ϵp_2{s}_{23}][M_3-ϵp_3{s}_{32}][1+\frac{\kappa }{4}{L}_1^{*}|z_{23}|]e^{\frac{\kappa }{4}{s}_{23}[(L_3+{\mathfrak{M}}_{31})z_{13}-(L_2+{\mathfrak{M}}_{13})z_{12}]}\right\}\hfill \\ & & =[M_2-ϵp_2{s}_{21}][(\frac{\partial M_1}{\partial p_1}-ϵs_{12})L_3^{*}-(M_1-ϵp_1{s}_{12})(ϵs_{13}+\frac{\kappa }{4}{L}_3^{*}\left(ϵz_{12}\right))]\hfill \\ & & \times e^{\frac{\kappa }{4}{s}_{12}[(L_1+{\mathfrak{M}}_{12})z_{13}-(L_2+{\mathfrak{M}}_{21})z_{23}]}\hfill \\ & & +[M_3-ϵp_3{s}_{31}][(\frac{\partial M_1}{\partial p_1}-ϵs_{13})L_2^{*}-(M_1-ϵp_1{s}_{13})\{ϵs_{12}+\frac{\kappa }{4}{L}_2^{*}\left(ϵz_{13}\right)\}]\hfill \\ & & \times e^{\frac{\kappa }{4}{s}_{13}[(L_1+{\mathfrak{M}}_{23})z_{12}+(L_3+{\mathfrak{M}}_{32})z_{23}]}\hfill \\ & & +[M_2-ϵp_2{s}_{23}][M_3-ϵp_3{s}_{32}][-s_{12}{s}_{13}\frac{\partial M_1}{\partial p_1}+\frac{\kappa }{4}{s}_{23}{L}_1^{*}[ϵ|z_{12}|-ϵ|z_{13}|]]\hfill \\ & & \times e^{\frac{\kappa }{4}{s}_{23}[(L_3+{\mathfrak{M}}_{31})z_{13}-(L_2+{\mathfrak{M}}_{13})z_{12}]}\hfill \\ & & +\frac{\partial M_1}{\partial p_1}{L}_2{L}_3+ϵ(s_{12}{L}_1{L}_3+s_{13}{L}_2{L}_1)\hfill \end{array}$$

(156)

for ${\dot{z}}_1$, with similar expressions for ${\dot{z}}_2$ and ${\dot{z}}_3$. Differentiating (152) with respect to $z_1$ yields

$$\begin{array}{ccc}& & {\dot{p}}_1\left\{L_2{L}_3+L_1{L}_3+L_1{L}_2\right.\hfill \\ & & -[M_2-ϵp_2{s}_{21}][M_1-ϵp_1{s}_{12}][1+\frac{\kappa }{4}{L}_3^{*}|z_{12}|]e^{\frac{\kappa }{4}{s}_{12}[(L_1+{\mathfrak{M}}_{12})z_{13}-(L_2+{\mathfrak{M}}_{21})z_{23}]}\hfill \\ & & -[M_3-ϵp_3{s}_{31}][M_1-ϵp_1{s}_{13}][1+\frac{\kappa }{4}{L}_2^{*}|z_{13}|]e^{\frac{\kappa }{4}{s}_{13}[(L_1+{\mathfrak{M}}_{23})z_{12}+(L_3+{\mathfrak{M}}_{32})z_{23}]}\hfill \\ & & \left.-[M_2-ϵp_2{s}_{23}][M_3-ϵp_3{s}_{32}][1+\frac{\kappa }{4}{L}_1^{*}|z_{23}|]e^{\kappa /4s_{23}[(L_3+M_3-ϵp_3{s}_{32})z_{13}-(L_2+M_2-ϵp_2{s}_{23})z_{12}]}\right\}\hfill \\ & & =[M_2-ϵp_2{s}_{21}][M_1-ϵp_1{s}_{12}]\left[\frac{\kappa }{4}{s}_{12}{L}_3^{*}[H+ϵ(p_2-p_1)s_{12}+ϵp_3{s}_{13}]\right]\hfill \\ & & \times e^{\frac{\kappa }{4}{s}_{12}[(L_1+{\mathfrak{M}}_{12})z_{13}-(L_2+{\mathfrak{M}}_{21})z_{23}]}\hfill \\ & & +[M_3-ϵp_3{s}_{31}][M_1-ϵp_1{s}_{13}]\left[\frac{\kappa }{4}{s}_{13}{L}_2^{*}[H+ϵp_2{s}_{12}+ϵ(p_3-p_1)s_{13}]\right]\hfill \\ & & \times e^{\frac{\kappa }{4}{s}_{13}[(L_1+{\mathfrak{M}}_{23})z_{12}+(L_3+{\mathfrak{M}}_{32})z_{23}]}\hfill \\ & & +[M_2-ϵp_2{s}_{23}][M_3-ϵp_3{s}_{32}]\left[\frac{\kappa }{4}{s}_{23}{L}_1^{*}{p}_1(s_{12}-s_{13})\right]e^{\frac{\kappa }{4}{s}_{23}[(L_3+{\mathfrak{M}}_{31})z_{13}-(L_2+{\mathfrak{M}}_{13})z_{12}]}\hfill \end{array}$$

(157)

for ${\dot{p}}_1$, again with similar expressions for ${\dot{p}}_2$ and ${\dot{p}}_3$.

6.2. Effective Potential

The determining Equation (152) indicates that the Hamiltonian is a function of only four independent variables: the two separations between the particles and their conjugate momenta. Writing

$$\begin{array}{cc}\hfill z_1-z_2& =\sqrt{2}\rho \hfill \end{array}$$

(158)

$$\begin{array}{cc}\hfill z_1+z_2-2z_3& =\sqrt{6}\lambda \hfill \end{array}$$

(159)

$$\begin{array}{cc}\hfill z_1+z_2+z_3& =Z\hfill \end{array}$$

(160)

in turn implies the conjugate momenta

$$\begin{array}{ccc}\hfill p_{\rho }& =& \frac{1}{\sqrt{2}}(p_1-p_2)\hfill \end{array}$$

(161)

$$\begin{array}{ccc}\hfill p_{\lambda }& =& \frac{1}{\sqrt{6}}(p_1+p_2-2p_3)\hfill \end{array}$$

(162)

$$\begin{array}{ccc}\hfill p_Z& =& \frac{1}{3}(p_1+p_2+p_3)\hfill \end{array}$$

(163)

obtained by imposing the requirement

$$\{q_{\mathrm{A}},p_{\mathrm{b}}\}={\delta }_{\mathrm{AB}}$$

(164)

where A,B$=\rho ,\lambda ,Z$.

The coordinates $\rho$ and $\lambda$ describe the motions of the three particles about their center of mass. The remaining conjugate pair $(Z,p_Z)$ respectively correspond to the centre of mass and its conjugate momenta in the non-relativistic limit. The Hamiltonian is independent of these variables and so $p_Z$ can be set to zero without loss of generality; similarly the origin can be chosen to be $Z=0$.

The Hamiltonian can be therefore be regarded as a function $H=H\left(\rho ,\lambda ,p_{\rho },p_{\lambda }\right)$ of the four canonical degrees of freedom. The determining equation can be rewritten as

$$L_{+}{L}_{-}{L}_0={\mathfrak{M}}_{+-}{\mathfrak{M}}_{-+}{L}_0^{*}{e}^{\frac{\kappa }{4}{s}_0{H}_0}+{\mathfrak{M}}_{-0}{\mathfrak{M}}_{01}{L}_{+}^{*}{e}^{\frac{\kappa }{4}{s}_{-}{H}_{-}}+{\mathfrak{M}}_{0+}{\mathfrak{M}}_{+0}{L}_{-}^{*}{e}^{\frac{\kappa }{4}{s}_{+}{H}_{+}}$$

(165)

where

$$\begin{array}{ccc}\hfill H_0& \equiv & (L_1+{\mathfrak{M}}_{12})z_{13}-(L_2+{\mathfrak{M}}_{21})z_{23}\hfill \\ & =& \sqrt{2}H\rho -ϵ\left(2\left|\rho \right|p_{\rho }+\left[\left|\lambda +\frac{\rho }{\sqrt{3}}\right|-\left|\lambda -\frac{\rho }{\sqrt{3}}\right|\right]p_{\lambda }\right)\hfill \end{array}$$

(166)

$$\begin{array}{ccc}\hfill H_{-}& \equiv & (L_2+{\mathfrak{M}}_{23})z_{21}-(L_3+{\mathfrak{M}}_{32})z_{31}\hfill \\ & =& H\left(\lambda -\frac{\rho }{\sqrt{3}}\right)-ϵ\left(\left[\left|\rho \right|-\frac{\sqrt{3}}{2}\left|\lambda +\frac{\rho }{\sqrt{3}}\right|\right]\left(\frac{p_{\lambda }}{\sqrt{3}}+p_{\rho }\right)+\frac{3}{2}\left(p_{\lambda }-\frac{p_{\rho }}{\sqrt{3}}\right)\left|\lambda -\frac{\rho }{\sqrt{3}}\right|\right)\hfill \end{array}$$

(167)

$$\begin{array}{ccc}\hfill H_{+}& \equiv & (L_3+{\mathfrak{M}}_{31})z_{32}-(L_1+{\mathfrak{M}}_{13})z_{12}\hfill \\ & =& -H\left(\lambda +\frac{\rho }{\sqrt{3}}\right)-ϵ\left(\left[\left|\rho \right|-\frac{\sqrt{3}}{2}\left|\lambda -\frac{\rho }{\sqrt{3}}\right|\right]\left(\frac{p_{\lambda }}{\sqrt{3}}-p_{\rho }\right)-\frac{3}{2}\left(p_{\lambda }+\frac{p_{\rho }}{\sqrt{3}}\right)\left|\lambda +\frac{\rho }{\sqrt{3}}\right|\right)\hfill \end{array}$$

(168)

and

$$\begin{array}{ccc}\hfill M_0& =& \sqrt{\frac{2}{3}}\sqrt{p_{\lambda }^2+m_0^2}{M}_{\pm }=\frac{1}{\sqrt{2}}\sqrt{{\left(\frac{p_{\lambda }}{\sqrt{3}}\pm p_{\rho }\right)}^2+2m_{\pm }^2}\hfill \end{array}$$

(169)

$$\begin{array}{ccc}\hfill L_{\pm }& =& H-M_{\pm }\pm \frac{ϵ}{\sqrt{2}}\left[\left(\frac{p_{\lambda }}{\sqrt{3}}\mp p_{\rho }\right)s_0\pm \frac{2}{\sqrt{3}}{p}_{\lambda }{s}_{\pm }\right]\hfill \end{array}$$

(170)

$$\begin{array}{ccc}\hfill L_0& =& H-M_0-\frac{ϵ}{\sqrt{2}}\left[\left(\frac{p_{\lambda }}{\sqrt{3}}+p_{\rho }\right)s_{+}+\left(\frac{p_{\lambda }}{\sqrt{3}}-p_{\rho }\right)s_{-}\right]\hfill \end{array}$$

(171)

$$\begin{array}{ccc}\hfill {\mathfrak{M}}_{\pm \mp }& =& M_{\pm }\mp \frac{ϵ}{\sqrt{2}}\left(\frac{p_{\lambda }}{\sqrt{3}}\pm p_{\rho }\right)s_0,{\mathfrak{M}}_{\pm 0}=M_{\pm }-\frac{ϵ}{\sqrt{2}}\left(\frac{p_{\lambda }}{\sqrt{3}}\pm p_{\rho }\right)s_{\pm }\hfill \end{array}$$

(172)

$$\begin{array}{ccc}\hfill {\mathfrak{M}}_{0\pm }& =& M_o+\frac{ϵ\sqrt{2}}{\sqrt{3}}{p}_{\lambda }{s}_{\pm }{L}_0^{*}=(1-s_{+}{s}_{-})M_0+L_0\hfill \end{array}$$

(173)

$$\begin{array}{ccc}\hfill L_{\pm }^{*}& =& (1\mp s_0{s}_{\pm })M_{\pm }+L_{\pm }\hfill \end{array}$$

(174)

It is instructive to carry out an expansion of (165) in inverse powers of the speed of light c, which is the post-Newtonian (pN) expansion. This gives

$$\begin{array}{ccc}\hfill H_{pN}& =& 3m{c}^2+\frac{p_{\rho }^2+p_{\lambda }^2}{2m}+\frac{\kappa m^2{c}^4}{\sqrt{8}}\left[\left|\rho \right|+\frac{\sqrt{3}}{2}\left(\left|\lambda +\frac{\rho }{\sqrt{3}}\right|+\left|\lambda -\frac{\rho }{\sqrt{3}}\right|\right)\right]-\frac{{(p_{\rho }^2+p_{\lambda }^2)}^2}{16m^3{c}^2}+\frac{\kappa c^2}{\sqrt{8}}|\rho |p_{\rho }^2\hfill \\ & & +\frac{\kappa c^2}{16\sqrt{2}}\left[3\left(\left|\lambda +\frac{\rho }{\sqrt{3}}\right|+\left|\lambda -\frac{\rho }{\sqrt{3}}\right|\right)\left(\sqrt{3}{p}_{\lambda }^2+p_{\rho }^2\right)+6\left(\left|\lambda +\frac{\rho }{\sqrt{3}}\right|-\left|\lambda -\frac{\rho }{\sqrt{3}}\right|\right)p_{\rho }{p}_{\lambda }\right]\hfill \\ & & +\frac{{\kappa }^2{m}^3{c}^6}{16}\left[\frac{\left|\rho \right|\sqrt{3}}{2}\left(\left|\lambda +\frac{\rho }{\sqrt{3}}\right|+\left|\lambda -\frac{\rho }{\sqrt{3}}\right|\right)+\frac{3}{4}\left|\lambda +\frac{\rho }{\sqrt{3}}\right|\left|\lambda -\frac{\rho }{\sqrt{3}}\right|-\frac{3}{4}\left({\lambda }^2+{\rho }^2\right)\right]\hfill \end{array}$$

(175)

where factors of the speed of light c have been restored (recall $\kappa =\frac{8\pi G}{c^4}$). To leading order in $1/c$

$$H=H_N\equiv 3m{c}^2+\frac{p_{\rho }^2+p_{\lambda }^2}{2m}+\frac{\kappa m^2{c}^4}{\sqrt{8}}\left[\left|\rho \right|+\frac{\sqrt{3}}{2}\left(\left|\lambda +\frac{\rho }{\sqrt{3}}\right|+\left|\lambda -\frac{\rho }{\sqrt{3}}\right|\right)\right]+\cdots$$

(176)

which is the non-relativstic hexagonal-well Hamiltonian of a single particle [115] plus the total rest mass of the system. This latter quantity, while non-relativistically irrelevant, is useful to retain in order to straightforwardly compare energies and motions and energies of the relativistic (R) and non-relativistic (N) systems.

The Hamiltonian (176) describes the non-relativistic motion of a single particle of mass m (referred to as the hex-particle) in a linearly-increasing potential well in the $\left(\rho ,\lambda \right)$ plane, whose cross-sectional shape is that of a regular hexagon. To extend this to the pN and R cases, the potential can be defined by the relation $V\left(\rho ,\lambda \right)=H\left(p_{\rho }=0,p_{\lambda }=0\right)$. From (165) we then obtain

$$\begin{array}{ccc}& & \left(V_{\mathrm{R}}-m_1{c}^2\right)\left(V_{\mathrm{R}}-m_2{c}^2\right)\left(V_{\mathrm{R}}-m_3{c}^2\right)\hfill \\ & & =\left(V_{\mathrm{R}}-s_{+}{s}_{-}{m}_3{c}^2\right)m_1{m}_2{c}^{4}exp\left[\frac{\sqrt{2}\kappa R}{4}{V}_{\mathrm{R}}\left|sin\theta \right|\right]\hfill \\ & & +\left(V_{\mathrm{R}}-s_{\rho }{s}_{+}{m}_1{c}^2\right)m_2{m}_3{c}^{4}exp\left[\frac{\sqrt{2}\kappa R}{4}{V}_{\mathrm{R}}\left|sin\left(\theta -\frac{\pi }{3}\right)\right|\right]\hfill \\ & & +\left(V_{\mathrm{R}}+s_{\rho }{s}_{-}{m}_2{c}^2\right)m_3{m}_1{c}^{4}exp\left[\frac{\sqrt{2}\kappa R}{4}{V}_{\mathrm{R}}\left|sin\left(\theta +\frac{\pi }{3}\right)\right|\right]\hfill \end{array}$$

(177)

for the R potential, where

$$\rho =Rsin\theta \lambda =Rcos\theta$$

(178)

has been used to render the hexagonal symmetry manifest, and $s_{\pm }=\mathrm{sgn}\left(\sqrt{3}\lambda \pm \rho \right)$, $s_{\rho }=\mathrm{sgn}\left(\rho \right)$. The corresponding pN potential is

$$\begin{array}{cc}\hfill {\tilde{V}}_{pN}& =(m_1+m_2+m_3)c^2+\frac{\kappa c^4}{4\sqrt{2}}\left(2m_1{m}_2\left|\tilde{\rho }\right|+m_1{m}_3\left|\sqrt{3}\tilde{\lambda }+\tilde{\rho }\right|+m_2{m}_3\left|\sqrt{3}\tilde{\lambda }-\tilde{\rho }\right|\right)\hfill \\ \hfill & +\frac{1}{2c^2}{\left(\frac{\kappa c^4}{4}\right)}^2{m}_1{m}_2{m}_3\left(\left(1-{\tilde{s}}_{\rho }{\tilde{s}}_1\right)\left|\tilde{\rho }\right|\left|\sqrt{3}\tilde{\lambda }+\tilde{\rho }\right|+\left(1+{\tilde{s}}_{\rho }{\tilde{s}}_2\right)\left|\tilde{\rho }\right|\left|\sqrt{3}\tilde{\lambda }-\tilde{\rho }\right|\right.\hfill \\ \hfill & \left.+\frac{1}{2}\left(1-{\tilde{s}}_1{\tilde{s}}_2\right)\left|\sqrt{3}\tilde{\lambda }+\tilde{\rho }\right|\left|\sqrt{3}\tilde{\lambda }-\tilde{\rho }\right|\right)\hfill \end{array}$$

(179)

where $\tilde{\rho }$ and $\tilde{\lambda }$ are defined as in (158) and (159) using the ${\tilde{z}}_a$ coordinates of (29). As $c\to \infty$, $\kappa \to 0$, and we recover

$$\begin{array}{cc}\hfill V_N& =(m_1+m_2+m_3)c^2+\frac{2\pi G}{\sqrt{2}}\left(2m_1{m}_2\left|\rho \right|+m_1{m}_3\left|\sqrt{3}\lambda +\rho \right|+m_2{m}_3\left|\sqrt{3}\lambda -\rho \right|\right)\hfill \end{array}$$

(180)

which is hexagonal well potential of the N-system.

A comparison of $V_N$ and $V_R$ is given in Figure 18. At very low energies they are indistinguishable, but striking differences emerge with increasing energy. For all energies equipotential lines of $V_N$ form the shape of a regular hexagon in the $\left(\rho ,\lambda \right)$ plane, with the sides rising linearly in all directions, forming the hexagonal-well potential noted earlier. The post-Newtonian potential $V_{pN}$ retains the hexagonal symmetry, but distorts the sides to be parabolically concave. The relativistic potential $V_R$ also retains the hexagonal symmetry, but the sides of the hexagon become convex, even at energies only slightly larger than the rest mass. The overall size of the hexagon at a given value of $V_{\mathrm{R}}$ is considerably smaller since its growth is extremely rapid compared to the other two cases. Its cross-sectional size reaches a maximum at $V_{\mathrm{R}}=V_{\mathrm{R}c}=6.711968022m{c}^2$, after which it decreases in diameter like $ln\left(V_{\mathrm{R}}\right)/V_{\mathrm{R}}$ with increasing $V_{\mathrm{R}}$. The relativistic potential is therefore an hexagonal carafe, whose neck narrows as $V_{\mathrm{R}}$ increases. The part of the potential for which $V_{\mathrm{R}}>V_{\mathrm{R}c}$ is in an intrinsically non-perturbative relativistic regime: the motion for values of $V_{\mathrm{R}}$ larger than this cannot be understood as a perturbation from some classical limit of the motion. A comparison of the equipotential lines for each case in Figure 19 highlights these distinctions.

Furthermore, since there are couplings between the momentum and position of the hex-particle, the potential does not fully govern the motion in both the pN and R systems. In the pN system there is a momentum-dependent steepening of the walls of the hexagon to leading order in $1/c^2$.

6.3. Relativistic Equal Mass 3-Body Trajectories

An analysis of the 3-body system is best carried out by considering the motion of the hex-particle in the $\left(\rho ,\lambda \right)$ plane. In the Hamiltonian formalism the motion is given by two conjugate pairs of differential equations for $(\rho ,p_{\rho })$ and $(\lambda ,p_{\lambda })$ that are continuous everywhere except across the three hexagonal bisectors $\rho =0$, $\rho -\sqrt{3}\lambda =0$, and $\rho +\sqrt{3}\lambda =0$. These bisectors respectively correspond to the crossings of particles 1 and 2, 2 and 3, or 1 and 3, and divide the hexagon into sextants.

The non-relativistic analogue of this system is that of a ball that elastically collides with a wedge whilst experiencing a constant gravitational force [3,16]. A bounce at one of the edges of the wedge corresponds to a crossing one of the three hexagonal bisectors, and a discrete mapping can be constructed that describes the particle’s angular and radial velocities each time it collides with one of the edges. The systems are nearly identical in the equal mass case since a crossing between two equal mass particles cannot be distinguished from an elastic collision between them.

The interesting dynamics of the system arises due to these crossings, or alternatively due to the non-smoothness of the potential along the bisectors. Two types of motion can be distinguished [3]: A-motion, in which the hex-particle crosses a single bisector twice in succession (the same pair of particles cross twice in a row) and B-motion, corresponding to the hex-particle crossing two successive sextant boundaries (one particle crosses each of its compatriots in succession). Any given motion is characterized by a symbol sequence, a sequence of letters A and B, with a finite exponent n denoting n-repetitions and an overbar denoting an infinite repeated sequence. For example the expression $A^2{B}^3$ denotes three A-motions (two adjacent particles cross twice in a row twice in succession) followed by three B-motions (one particle crosses the other two in succession, which then cross each other). The expression ${\left(A^n{B}^m\right)}^p$ denotes the motion $A^n{B}^m$ repeated p times, and $\overline{\left(A^n{B}^m\right)}$ denotes infinitely many repeats of this motion. There is an ambiguity in classifying either the final or the initial crossing since the whether or not a motion is of A-type or B-type is contingent upon the previous crossing; this ambiguity can be resolved by taking the initial crossing of any sequence of motions as being unlabeled.

Since the same initial conditions for the N, pN, and R systems do not yield the same conserved energy H, comparison of these cases necessitates a choice: one can either compare at fixed values of the energy (FE conditions), modifying the initial conditions as appropriate (as required by the conservation laws for each system), or else fix the initial conditions, comparing trajectories at differing values of H (fixed-momenta (FM) conditions). Numerically it is useful to rescale

$$p_i=M_{tot}c\widehat{p_i}{z}_i=\frac{4}{\kappa M_{tot}{c}^2}\widehat{z_i}$$

(181)

in which case the equations of motion become

$$\begin{array}{ccc}\hfill \frac{\partial \eta }{\partial \widehat{p_i}}& =& \frac{1}{c}\frac{\partial H}{\partial p_i}=\frac{4}{\kappa M_{tot}{c}^3}\frac{d{\widehat{z}}_i}{dt}=\frac{d{\widehat{z}}_i}{d\widehat{t}}\hfill \end{array}$$

(182)

$$\begin{array}{ccc}\hfill \frac{\partial \eta }{\partial \widehat{z_i}}& =& \left(\frac{4}{\kappa M_{tot}^2{c}^4}\right)\frac{\partial H}{\partial z_i}=-\frac{4}{\kappa M_{tot}{c}^3}\frac{d\widehat{p_i}}{dt}=-\frac{d\widehat{p_i}}{d\widehat{t}}\hfill \end{array}$$

(183)

where $M_{tot}=3m$ is the total mass of the system, $t=\frac{4}{\kappa M_{tot}{c}^3}\widehat{t}$, and ${\widehat{z}}_i$ and ${\widehat{p}}_i$ are the respective dimensionless positions and momenta. The quantities $M_i$ and $L_i$ likewise rescale

$$\begin{array}{ccc}\hfill M_i& =& M_{tot}{c}^2\left(\sqrt{{\widehat{p_i}}^2+{\widehat{m_i}}^2}+\widehat{p_i}\right)=M_{tot}{c}^2{\widehat{M}}_i\hfill \end{array}$$

(184)

$$\begin{array}{ccc}\hfill L_i& =& M_{tot}{c}^2\left(\eta +1-\sqrt{{\widehat{p_i}}^2+{\widehat{m_i}}^2}-ϵ(\sum _j{\widehat{p}}_j{s}_{ji})\right)=M_{tot}{c}^2{\widehat{L}}_i\hfill \end{array}$$

(185)

where $\eta +1=H/M_{tot}{c}^2$, $\widehat{m_i}=\left(\frac{m_i}{M_{tot}}\right)$.

One consideration in describing the motion is that the proper time (97) differs for each particle, even in the equal mass case. The simplest choice (but not the only one) is to work with the coordinate time t.

The plots in this section were obtained numerically [29,117], with a time step in the numerical code having a value $\widehat{t}=1$. Absolute and relative error tolerances of ${ϵ}_{abs}={ϵ}_{rel}={10}^{-8}$ were imposed so that the estimated error in each of the dynamical variables $\rho \left(i\right),\lambda \left(i\right),p_{\rho }\left(i\right)$, and $p_{\lambda }\left(i\right)$ at each step i in the numerical integration is $ϵ\left(i\right)\le max\left({ϵ}_{rel}\left|y\left(i\right)\right|,{ϵ}_{abs}\right)$.

6.3.1. General Features of the Motion

For each of the N, pN, and R systems the motions fall into one of three principal classes: annuli, pretzel and chaotic. Within each class the orbits either (i) eventually densely cover the portion of $\left(\rho ,\lambda \right)$ space they occupy, or (ii) do not. A symbol sequence consisting of a finite sequence repeated infinitely many times would be in case (ii) whereas all chaotic orbits (by definition) are in case (i). Quasi-regular orbits are also in case (i); for these the symbol sequence consists of repeats of the same finite sequence, but with an A-motion added or removed at irregular intervals. In phase space the two types of orbits are separated by separatrixes (trajectories joining a pair of hyperbolic fixed points). Regular orbits lie inside the ‘elliptical’ region surrounding an elliptical fixed point whereas quasi-regular orbits lie outside such a region.

Quasi-periodic trajectories closely resemble their related periodic counterparts, except that the orbit does not exactly repeat itself. Consequently a quasi-periodic orbit eventually densely covers some region of phase space despite its high degree of regularity, manifest by its periodic symbol sequence. A particle moving on a torus $S^1\times S^1$ is a textbook example. Its motion is characterized by its angular velocity around each copy of $S^1$; if the ratio of these is rational, the motion will be periodic, whereas the motion will be quasi-periodic if the ratio is irrational. For the 3-body case non-periodic orbits with fixed symbol sequences are quasi-periodic, and appear as a collection of closed circles, ovals, or crescents in the Poincaré section (discussed in the next subsection). Quasi-regular orbits have symbol sequences that are not fixed.

6.3.2. Annulus Orbits

Annulus orbits have the symbol sequence $\overline{B}$ and consist of an annulus encircling the origin in the $\rho -\lambda$ plane. In these orbits the hex-particle never crosses the same bisector twice in succession.

Most annulus orbits are quasi-periodic and fill in the entire ring. However a few repeat themselves after some number of rotations about the origin, and a wide variety of patterns are possible contingent upon the initial conditions for the N, pN, and R systems. No qualitative distinctions between N and pN annuli were observed within numerically attainable values of $\eta$ [28,29].

An example of annulus orbits for the N and R systems (for FE conditions), along with the positions of each of the three bodies as a function of time, is shown in Figure 20 (periodic) and Figure 21 (near-chaotic). Periodic orbits are numerically difficult to find, so the orbits in Figure 20 are actually very close to periodic orbits; this allows the pattern of the periodic orbit to be visible. At similar energies the R hex-particle covers the $\left(\rho ,\lambda \right)$ plane more densely than its N counterpart, and has a higher frequency of oscillation. The higher frequency for the R system was seen in the previous section for two bodies and appears to hold generally for the 3-body system as well. The increased trajectory density for FE conditions in the R system consequently follows, since the same number of time steps were used for both.

These features of higher frequency and trajectory density are more apparent in Figure 22, which provides a comparison of orbits using FM initial conditions. For these conditions the R system (blue) has slightly higher energy than its N counterpart (red). The N system is less dense, covers a smaller region of the $(\rho ,\lambda )$ plane, and does not venture as close to the origin, characteristics that becoming increasingly pronounced for increasing $\eta$, provided that the R energy remains larger than its N counterpart. This is not guaranteed for FM conditions, as the bottom diagram in Figure 22 illustrates: here the N system has about 14% more energy than its R counterpart, and so covers a larger region of the $(\rho ,\lambda )$ plane.

6.3.3. Pretzel Orbits

Pretzel orbits have the generic symbol sequence ${\prod }_{ijk}{\left(A^{n_i}{B}^{3m_j}\right)}^{l_k}$ , where $n_i,m_j,l_k\in Z^{+}$, with some $l_k$ possibly infinite, and consist of orbits in which the hex-particle essentially oscillates back and forth about one of the three bisectors for some segments of its motion. A typical example is shown in Figure 23. For both the N and R systems, we see that two of the three bodies form a stable bound subsystem, which in turn orbits the third analogous to a 2-body system. The N system exhibits parabolic regularity for both the 2-body subsystem and the full system, whereas the R system has shoulder-like distortions observed previously in the 2-body case.

This formation of a stable (or quasi-stable) bound subsystem is characteristic of pretzel orbits, and the range of possible trajectories is extremely diverse. Many families of regular orbits exist. These generally have one base element in their symbol sequence (e.g., $A{B}^3$) and a sequence of elements formed by appending an A to each existing sequence of A’s (for example, $\{A{B}^3,A^2{B}^3,A^3{B}^3,...\}$. The $B^3$ sequence corresponds to a 180-degree rotation of the hex-particle about the origin, yielding a broad variety of twisted, pretzel-like figures. This is a notable distinction from the wedge system [3,16], for which B and $B^2$ sequences are also observed; only $B^3$ sequences are present in all pretzel orbits.

Distinctions between the R, pN and N systems are strongest for pretzel orbits. Both regular orbits (with an infinitely repeating symbol sequence) and irregular orbits that densely fill a cylindrical tube in the $\left(\rho ,\lambda \right)$ plane occur. Orbits in the R system generally have kinks about the $\lambda =0$ line that are absent in their N and pN counterparts; a cylindrical-shaped trajectory in the N system looks like an hourglass in the R system, for example. Furthermore, periodic and quasi-periodic orbits in the N system have counterparts with the same symbol sequence in the R system but not in the pN system, which exhibits chaotic behaviour not seen in the N and R systems.

A comparison of the time-evolution of trajectories in the N and R systems is shown in Figure 24 for FE conditions at small and large values of $\eta$. For small $\eta$ ($\eta =0.05$) there is very little distinction between the N and R motions, consistent with the smooth non-relativistic limit of (152). Significantly different trajectories occur for larger values $\eta$ ($\eta =0.85$). In the R system the oscillation frequency is higher and the pattern ‘weaves’ relative to the near-cylindrical shape in the N-system once enough time steps have occured. The R trajectory is more tightly confined, commensurate with the 2-body motion seen in the previous section.

In Figure 25 we display the sensitivity of trajectories to initial conditions. The fish-like diagram has an $A{B}^6$ symbol sequence: two of the particles oscillate quasi-regularly about each other (shown in the upper $z\left(t\right)$ plot), with this pair undergoing larger-amplitude and lower-frequency oscillations with the third. A slight change of the initial FE conditions yields the strudel-like figure on the left. Now one particle alternates its oscillations with the other two, maintaining a near-constant amplitude throughout.

By controlling the FM initial conditions interesting sequences of hex-particle orbits can be obtained. An example is given in Figure 26, which compares snake-like orbits in the N (red) and R (blue) systems. These quasi-regular orbits have symbol sequences $\left\{\overline{A^i{B}^3}\right\}$. In both systems the orbits have two sharp turning points separated by some number n of bumps. In the N system these have been shown to exist for arbitrary n [3], and it has been conjectured the same is true for the R and pN systems [29]. The figures in the R system develop an hourglass shape narrowing about $\lambda =0$, and cover a much narrower region in the $\rho$ direction (note the scale in the bottom sequence of plots). The N orbits, by contrast, are circumscribed by a cylinder.

The symbol sequence $A{B}^3$ results in boomerang-like figures, shown in Figure 27, where the qualitatively different physics due to relativistic effects is manifest. At low energies ($\eta =0.2$) the N (red) and R(blue) systems have similar boomerang shapes. However as $\eta$ increases, orbits in the R system develop two distinct turning points at different distances from the $\rho =0$ axis for $\lambda >0$, with symmetric counterparts for $\lambda <0$. This feature is particularly evident for $\eta =0.75$. A kink at the right-hand-side of the boomerang emerges, becoming increasingly pronounced with increasing $\eta$. The underlying reason behind the development of this structure is not clear.

6.3.4. Chaotic Orbits

Chaotic orbits are those for which the hex-particle wanders between A-motions and B-motions in an apparently irregular fashion. Unlike the annuli and pretzel orbits, chaotic orbits eventually wander into all areas of the $\rho -\lambda$ plane. Chaos emerges at the transition between annulus and pretzel orbits, where the hex-particle passes very close to the origin, for each system.

Figure 28 illustrates a typical case for the N and R systems. Two particles undergo a large-amplitude oscillation with the third one (the middle ‘m’ particle) mildly oscillating near the centre of momentum. At irregular intervals this third particle switches places with one of the other two, and the pattern repeats. The m-particle alternates in an irregular fashion, leading to chaos.

A comparison of the time development of chaotic trajectories in both the N and R systems is shown in Figure 29. The upper sequence shows how a chaotic trajectory can develop in the R system (blue) whilst the N system (red) forms a densely filled annulus for the same FE initial conditions with $\eta =0.5$. The lower sequence shows how a chaotic trajectory can develop in the N system (red) whilst the R system (blue) forms a densely filled pretzel, for the same FE initial conditions, again for $\eta =0.5$. In both cases the R trajectory attains its final pattern much more rapidly than its N counterpart, a manifestion of the difference in frequencies noted earlier.

The transition from an annulus to a pretzel orbit through a chaotic region for the R system is shown in Figure 30. Proceeding from from left-to-right and top-to-bottom with decreasing initial angular momentum for the hex-particle, the system begins as an annulus, passes through a set of chaotic orbits, and then becomes a pretzel. The chaotic trajectories through the origin (or very close to it), a characteristic feature for this region of chaos in all three systems. In the R system the transitional region shrinks as $\eta$ increases [29].

There is a striking distinction between the pN system and the (N,R) systems for the chaotic class of motions. Unlike the latter two, the pN exhibits an additional area of chaos in the pretzel region. This feature will become evident in the next subsection.

6.4. Poincaré Plots

Poincaré sections for each of the N, pN, and R systems can be constructed by plotting the square of the angular momentum $p_{\theta }^2$, against the radial momentum $p_R$ (both scaled to be dimensionless as per (181)) of the hex-particle each time it crosses a bisector. Since all bisectors are equivalent in the equal mass case, all crossings can be plotted on the same surface of section. The plots then indicate regions of chaos, as well as periodicity and quasi-periodicity.

Since the Hamiltonian is time-independent, the total energy is a constant of the motion, and the phase space is a 3-dimensional hypersurface in 4 dimensions for each system. The system is said to be integrable if an additional constant of the motion exists, in which case its trajectories are restricted to two-dimensional surfaces in the available phase space.

The types of motion that integrable systems can exhibit are either periodic or quasi-periodic. Periodic (1-dimensional) orbits have trajectories that always appear as lines or dots on the Poincaré section, since by definition they comprise the intersection of two 2-dimensional surfaces. By contrast, the extra degree of freedom for a non-integrable system permits orbits to visit all regions of phase space. In this case the system typically displays strongly chaotic behavior and the associated trajectories appear as filled in areas on the Poincaré plot.

Small perturbations of an integrable system admit small regions of chaos, though most of its orbits remain confined to two-dimensional surfaces. The chaotic regions grow as the perturbation increases in magnitude and eventually become connected areas on the Poincaré section, a phenomenon known as a Kolmogorov, Arnold and Moser (KAM) transition [120,121,122,123]. For sufficiently large perturbations, systems typically become almost fully ergodic [124], though islands of regularity may persist for quite some time prior to this and typically have an intricate fractal structure.

The Poincaré plots for the N, R and pN systems appear respectively in Figure 31, Figure 32, and Figure 33, using the same conventions as in [3] up to an overall normalization for each section. The energy conservation relation (176), which is

$$p_R^2+p_{\theta }^2\le \frac{2}{3}\eta$$

(186)

determines the outer boundary of each plot, where the presence of $\eta$ reflects the different normalizations for each system. When the potential energy is zero, equality in (186) holds, corresponding to the hex-particle being at the origin. Any departure from the origin reduces the values of $\left(p_R,p_{\theta }\right)$ and so yields the phase-space limit.

Figure 31 for the N system reproduces the results for the wedge problem [3]. Here $H-3m{c}^2$ is normalized to unity and the RHS of (186) is 1. In the equal mass case the 3-body problem corresponds to motion of a body falling toward a wedge whose sides are each at angles ${30}^o$ relative to the vertical axis. The energy constraint after an A-collision has taken place yields another boundary

$${\left(p_R-2\sqrt{3}|p_{\theta }|\right)}^2\le 1-p_{\theta }^2$$

(187)

whose satisfaction yields all points in phase space that have undergone an A-collision (the A-region). Points violating (187) are those for which a B-collision has taken place (the B-region). Since the interaction is gravitational, collisions with the third particle cannot ultimately be avoided. Consequently the A-region has no fixed points and any point in the A-region will inevitably venture into the B-region. However fixed points can occur in a subset of the B-region: here B-collisions are infinitely repeated, corresponding to the annulus orbits.

In the N system the centre of the plot in Figure 31 is a fixed point surrounded by a subregion of near-integrable curves. The annulus orbits are all contained within the large triangle surrounding this region. The closed circles in this annulus region correspond to quasi-periodic orbits about the periodic annuli with higher period, such as in Figure 20. The boundary of the annulus region is a thin region of chaos, most prominent at the corners as shown in the lower right inset. These chaotic regions are confined and not simply connected.

The region beyond this is the pretzel region, which has circles bounding the quasi-periodic near-integrable regions; these exhibit repeated self-similarity, shown in the upper right inset. The two large circles observed just below the annulus region correspond to the boomerang-shaped orbits ($\overline{A{B}^3}$) shown in Figure 27. The next set of circles will be $\overline{A^2{B}^3}$, and so on. Between these sets of circles are collections of crescents with sequences $\overline{A{B}^3{A}^2{B}^3}$, $\overline{A{B}^{3}A{B}^3{A}^2{B}^3}$, etc. Each circle is actually a continuum of possible circles, whose diameter depends on the initial conditions. At the center of this family of circles is a dot corresponding to the periodic orbit in question.

Somewhat remarkably, the highly nonlinear R system, shown for two different values of $\eta$ in Figure 32 retains all of the qualitative features of the N system, at least up to the range of $\eta$ numerically accessible. The annulus, pretzel and chaotic regions all retain their same basic structure, though asymmetrically deformed. This deformation increases as $\eta$ increases and occurs because the Hamiltonian given by Equation (152) is invariant under the discrete symmetry $\left(p_i,ϵ\right)\to \left(-p_i,-ϵ\right)$ instead of the $p_i\to -p_i$ symmetry in the N system. The discrete constant $ϵ=\pm 1$ is a measure of the flow of time of the gravitational field relative to the particle momenta. In Figure 32 $ϵ=+1$; the opposite choice would create the same distortion but towards the lower left. The situation is reminiscent of the 2-body case of the previous section, where the gravitational coupling to the kinetic-energy causes a distortion of the trajectory from an otherwise symmetric pattern. Whether or not KAM breakdown occurs for higher $\eta$ values remains an interesting open question.

As compared to the N and R plots, the pN system, shown in Figure 33 is notably different. It retains the $p_i\to -p_i$ symmetry of the N system but appears to undergo a KAM transition at $\eta \sim 0.3$, shown in Figure 34. For small $\eta$ the distinction with the N system is mild, but for $\eta =.21$, the lines across the bottom of the figure slightly widen. Larger regions of chaos become evident around the edges of the groups of ellipses in the lower regions of the figure for $\eta =0.26$. Further increasing $\eta \to .3$, the lower part of the Poincaré section becomes engulfed by a chaotic sea; only a few non-connected islands of regular motion remain. No such behavior is seen in R system for similar values of $\eta$. These distinctions are not artifacts due to differences in scalings between the systems. The underlying feature that enforces structure on the phase space in the R system but not in the pN system remains to be understood.

6.5. Unequal Masses

For unequal masses the hexagon becomes squashed, with two opposite corners moving inward, changing both the slopes of the straight edges and their relative lengths; relativistic corrections maintain this basic distortion, but with the straight edges becoming parabolic [117]. The exact R potential is given by (177), with its pN and N counterparts respectively given by (179) and (180).

The R potential is similar to the N potential except that the sides of the hexagon become concave, and have a much steeper slope as the radial variable in the $(\rho ,\lambda )$ space increases. For $V_R={\widehat{V}}_R$ such that

$$ln\left(\frac{\left({\widehat{V}}_R-m_j{c}^2\right)\left({\widehat{V}}_R-\left(M_{\mathrm{tot}}-m_j\right)c^2\right)}{\left(M_{\mathrm{tot}}-m_j\right)m_j{c}^4}\right)={\widehat{V}}_R\left[\frac{1}{\left({\widehat{V}}_R-m_j{c}^2\right)}+\frac{1}{\left({\widehat{V}}_R-\left(M_{\mathrm{tot}}-m_j\right)c^2\right)}\right]$$

(188)

(for $j=1,2$ or 3), the slope of the R potential becomes infinite. In the equal mass case this occurs at ${\widehat{V}}_R\approx 6.71197m{c}^2$, where $m=M_{\mathrm{tot}}/3$. The maximal possible critical value of the potential occurs when one of the masses $m_j=M_{\mathrm{tot}}/2$, for which ${\widehat{V}}_R\approx 6.886682m_j{c}^2$. For $m⟶0,M_{\mathrm{tot}}$ the potential ${\widehat{V}}_R⟶M_{\mathrm{tot}}{c}^2$ – no energy is available for motion. A plot of ${\widehat{V}}_R$ as a function of $m_j$ is shown in Figure 35. For values of $V_R>{\widehat{V}}_R$ the size of the distorted hexagon decreases like $\left(ln{V}_R\right)/V_R$.

The overall shape of the R potential is that of a distorted hexagonal carafe, whose distortion is analogous to that of the N potential, as shown in Figure 36. The relatively steeper growth of the R potential as a function of distance from the origin is manifest by the smaller scales for $\rho$ and $\lambda$ in the right-hand panel.

6.5.1. Trajectories

Expressing relative masses as a ratio $m_1:m_2:m_3$, annulus, pretezel and chaotic trajectories are obtained and the symbol sequence always takes the form

$$\prod _{i,j,k}{(A^{m_i},B^{3n_j})}^{l_k}$$

(189)

as in the equal mass case. Equation (189) implies that B motion always comes in multiples of three; no motions have been observed that depart from this situation [117].

The situation where two masses are equal ($m_1:m_2:m_3=1:1:\alpha$) is instructive, and plots of the relative motion of the particles for both the N and R systems is given in Figure 37 for decreasing values of $\alpha$. For equal mass ($\alpha =1$) there is annulus motion: no particle ever crosses another twice in a row. But as $m_3$ decreases ($\alpha$ decreases) its frequency of oscillation decreases while its amplitude increases with respect to the other two, which provide a gravitationally bound subsystem. The binding becomes tighter as $\alpha$ decreases, more so for the R system. Eventually the binding becomes so tight that the more massive particles execute an additional A motion before crossing the third particle, and the hex-particle transitions from annulus to pretzel motion. The greater the mass difference, the more difficult it is to set up initial conditions at a given energy so that particles 1 and 2 do not cross more than once during the long period oscillation of particle 3.

Large and small values of $\alpha$ are likewise instructive, and plots for $\alpha =100$ (particle 3 very massive) and $\alpha =0.01$ (particles 1 and 2 very massive) are respectively shown in Figure 38 and Figure 39. As expected, for large $\alpha$ the heavy particle 3 barely moves as the other two oscillate about it, depicted in Figure 38. However the passing of the other two particles causes small perturbations in the motion of the heavy body, shown in the insets. The perturbation is smooth and regular in the N system whereas in the R system the velocity of the large mass increases much more suddenly, leading to a more erratic and jerky trajectory.

For small $\alpha$ (Figure 39) a stable gravitationally bound subsystem is formed by the two heavy particles, with the third oscillating about their center of inertia. The oscillation amplitude is much larger and frequency much smaller in the N system than in the R system, commensurate with the two-body system in Section 5. The effect of the light particle 3 is to cause oscillatory motion of the center of mass of the two more massive particles, which is clear from the insets in the the N and R systems. This perturbation is almost imperceptible due to the two heavy particles being twice as massive as the single particle in the $\alpha =100$ case (Figure 38); consequently they are less susceptible to changes in motion from the smaller mass body.

6.5.2. Poincaré Plots

As with the equal mass case, the annulus and pretzel trajectories are in similar regions of the Poincaré plot, separated by a region of chaos, but the shapes and sizes of the different regions change. The R system is topologically similar to the N system, but with the various regions distorted in a manner similar to the equal mass case.

Figure 40 compares the Poincaré sections for the N and R systems for the mass ratio $1:1:0.1$ for $\eta =0.3$. The triangular annulus region moves towards the top the surface of section and becomes smaller in both systems. This latter effect is a manifestation of the difficulty in attaining annulus motion noted above, when one particle is much less massive than the other two.

If one particle is much more massive than the other two, the annulus region gets larger and extends towards the lower region of the plot, as shown in Figure 41 in both the R and N systems for the mass ratio $1:1:10$. In this situation the hex-particle needs less angular momentum to attain an annulus orbit in the $(\rho$-$\lambda )$ plane. Somewhat like a two-planet solar system in one-dimension, the two lighter bodies behave like two separate 2 body systems, with the heavy particle taking the role of the second body for each, as shown in Figure 38. Furthermore, additional regions of chaos appear that are absent in the equal mass case.

An example of a Poincaré plot when all three masses are unequal is shown in Figure 42 for the N and R systems, with mass ratio $1:5:10$. Since none of the bodies have the same mass, the symmetry about the $p_R=0$ axis in the N system is gone. Different regions are not as clearly segregated as in the $m_1=m_2$ case, and instead extend over a larger region of the plot. A region of chaos separating outer pretzel regions from the inner annulus region is marked by A in the left diagram, and above and below B are new regions of chaos amongst pretzel trajectories.

The R system further distorts the N diagram to the lower right. The chaotic region separating annulus and pretzel trajectories is now two loops that were created by a single trajectory, both marked by a 1. The annulus region is confined to the area inside both of these loops (marked by a 2), where a single annulus trajectory will visit both regions.

The key feature of the unequal mass case, for both the N and R systems, is the presence of additional chaotic regions that are absent in the equal mass case in the corresponding constant energy hyper-surface. These additional chaotic regions appear within the pretzel regions of the corresponding equal mass plot, and are characterized by broadened lines in the pretzel region, evident in each of Figure 40, Figure 41 and Figure 42. The origin of this additional chaos is not understood.

The unequal mass case is equivalent to the two-dimensional symmetric wedge billiard system in a uniform gravitational field [3], with the relative masses of the particles directly related to the wedge angle $\theta$ by

$$tan\theta =\frac{\sqrt{1+2{\alpha }^{-2}}}{1+2{\alpha }^{-1}}$$

(190)

where $\alpha$ is as defined as above, and $\theta =\pi /6$ corresponds to $\alpha =1$, the equal mass case. The angle of the wedge is related to the angle between the bisectors of the hexagonal well. The only distinction between the wedge system and the 3-body system is the absence of collisions in the the latter. If all masses are equal there is no distinction between a collision and a crossing of two bodies (apart from particle labelling) and so the Poincaré maps become identical in this case.

In both the N and R systems, there is an increase in the amount of chaos as the difference in the masses increases. Earlier studies explored a limited range of mass ratios, however, and it remains an open question as to whether or not one or both systems will undergo a transtion to global chaos, or if integrable and near integrable regions exist for all mass ratios.

6.6. Charge and Cosmological Constant

Including a cosmological constant [118] and endowing the particles with charges [125] significantly increases the parameter space. The solutions to (63) and (64) still yield (65), but now

$$V=\left\{\begin{array}{c}0\\ -e_1{e}_2-e_1{e}_3\\ -e_1{e}_3-e_2{e}_3\\ 0\end{array}\begin{array}{c}\mathrm{in}\mathrm{region}\left(1\right)\\ \mathrm{in}\mathrm{region}\left(2\right)\\ \mathrm{in}\mathrm{region}\left(3\right)\\ \mathrm{in}\mathrm{region}\left(4\right)\end{array}\right\}$$

(191)

where space is divided into four regions, chosen so that $z_1<z_2<z_3$ region (1) being to the left of particle 1 and region (4) to the right of particle 3. The determining equation becomes

$$\begin{array}{cc}\hfill & \left[\left(\left({\widehat{M}}_1+{\widehat{K}}_{1+}\right){\widehat{K}}_{3+}^2+\left({\widehat{M}}_3+{\widehat{K}}_4\right){\widehat{K}}_{2+}^2\right)tanh\left(\frac{{\widehat{K}}_{3+}{z}_{32}}{4}\right)tanh\left(\frac{{\widehat{K}}_{2+}{z}_{21}}{4}\right)\right.\hfill \\ \hfill & +\left({\widehat{M}}_1+{\widehat{K}}_{1+}\right)\left({\widehat{M}}_3+{\widehat{K}}_4\right){\widehat{M}}_{2}tanh\left(\frac{{\widehat{K}}_{3+}{z}_{32}}{4}\right)tanh\left(\frac{{\widehat{K}}_{2+}{z}_{21}}{4}\right)\hfill \\ \hfill & +\left(\left({\widehat{M}}_1+{\widehat{M}}_2+{\widehat{K}}_{1+}\right)\left({\widehat{M}}_3+{\widehat{K}}_4\right)+{\widehat{K}}_{3+}^2\right){\widehat{K}}_{2+}tanh\left(\frac{1}{4}{\widehat{K}}_{3+}{z}_{32}\right)\hfill \\ & +\left(\left({\widehat{M}}_1+{\widehat{K}}_{1+}\right)\left({\widehat{M}}_2+{\widehat{M}}_3+{\widehat{K}}_4\right)+{\widehat{K}}_{2+}^2\right){\widehat{K}}_{3+}tanh\left(\frac{1}{4}{\widehat{K}}_{2+}{z}_{21}\right)\hfill \\ & \left.+\left({\widehat{M}}_1+{\widehat{M}}_2+{\widehat{M}}_3+{\widehat{K}}_{1+}+{\widehat{K}}_4\right){\widehat{K}}_{2+}{\widehat{K}}_{3+}\right]\hfill \\ \hfill & =0\hfill \end{array}$$

(192)

where $z_{ij}\equiv (z_i-z_j)$, $s_{ij}\equiv \mathrm{sgn}\left(z_{ij}\right)$ and

$$\begin{array}{ccc}\hfill {\widehat{M}}_i& =& \kappa \sqrt{p_i^2+m_i^2}\hfill \end{array}$$

(193)

$$\begin{array}{ccc}\hfill {\widehat{K}}_{1\pm }& \equiv & \sqrt{{\kappa }^2{\left[X-\frac{ϵ}{4}\left(\sum _{a=1}^3{p}_a{s}_{a1}\pm p_1\right)\right]}^2-\frac{{\Lambda }_e}{2}}\hfill \end{array}$$

(194)

$$\begin{array}{ccc}\hfill {\widehat{K}}_{2\pm }& \equiv & \sqrt{{\kappa }^2{\left[X-\frac{ϵ}{4}\left(\sum _{a=1}^3{p}_a{s}_{a2}\pm p_2\right)\right]}^2-\frac{\kappa }{2}(e_1{e}_2+e_1{e}_3)-\frac{{\Lambda }_e}{2}}\hfill \end{array}$$

(195)

$$\begin{array}{ccc}\hfill {\widehat{K}}_{3\pm }& \equiv & \sqrt{{\kappa }^2{\left[X-\frac{ϵ}{4}\left(\sum _{a=1}^3{p}_a{s}_{a3}\pm p_3\right)\right]}^2-\frac{\kappa }{2}(e_1{e}_3+e_2{e}_3)-\frac{{\Lambda }_e}{2}}\hfill \end{array}$$

(196)

$$\begin{array}{ccc}\hfill {\widehat{K}}_4& =& -2\sqrt{{\kappa }^2{\left[X+\frac{ϵ}{4}\left(p_1+p_2+p_3\right)\right]}^2-\frac{{\Lambda }_e}{2}}\hfill \end{array}$$

(197)

with

$$H=\frac{4}{\kappa }\sqrt{{\kappa }^2{X}^2-\frac{{\Lambda }_e}{2}}\ge \frac{4}{\kappa }\sqrt{-\frac{{\Lambda }_e}{2}}$$

(198)

which defines the constant of integration X. This relation is the same as (76) when the charges are zero. The inequality follows since H and $X\ge 0$ must both be real, defining a negative critical value

$${\Lambda }_{negcrit}\equiv -\frac{{\kappa }^2{H}^2}{8}\le {\Lambda }_e$$

(199)

Permutation of the particles yields same determining equation with indicies appropriately switched.

There are six inequivalent configurations: (0 0 0), (+ 0 0), (+ + 0), (+ - 0), (+ + -) and (+ + +), where 0 denotes a neutral particle, since the potential is invariant for $e_i\to -e_i$ and the charge configurations interchange as the particles cross. The first of these has been considered in the previous subsections. The second of these, unless ${\Lambda }_e=0$, is equivalent to three neutral particles with a cosmological constant. If one particle is neutral the relative magnitudes of the charges are irrelevant since only products of their magnitudes matter. The effect of the electromagnetic energy is to yield a charge-dependent constant vacuum energy between the particles, with the size of the region and the magnitude of the vacuum energy changing as the particles move about. If one particle is neutral then the magnitude of the vacuum energy between the particles does not change.

The parameter space has been explored in considerable detail [118,125]. As it is so vast only a few key results shall be presented.

6.6.1. Neutral Configurations with ${\Lambda }_e\ne 0$

For a nonzero cosmological constant, the case of three neutral particles is equivalent to that of the (+ 0 0) configuration. In general ${\Lambda }_e$ significantly modifies the chaotic properties of the relativistic 3-body system in markedly different ways, depending on its sign. The following analysis sets all masses to be equal.

For ${\Lambda }_e<0$ there is a rapid decrease in the size of the chaotic regions. These become even even smaller than in the N system in the Poincaré plot. Given the high degree of non-linearity this is quite surprising. This is manifest even at fairly small energies, illustrated in Figure 43 for $H=1.2M_{\mathrm{tot}}{c}^2$.

The upper diagram in Figure 43 corresponds to the ${\Lambda }_e=0$ case, similar to the lower diagram in Figure 32 but with tracks (green curves) showing how different types of orbits move as ${\Lambda }_e$ changes from a negative value (lower left diagram) to a positive one (lower right diagram). These four different quasiperiodic (or stable) orbits are shown in Figure 44 and change as ${\Lambda }_e$ changes. Orbit (a) is located in the centre of the annulus region, whereas orbits (b), (c), and (d) produce a set of points on the Poincaré plot that follow the contours of the triangular shaped region. The stable and quasistable orbits remain so as ${\Lambda }_e$ becomes more negative, but ${\Lambda }_e>0$ stable orbits can become chaotic. The transition point depends on the initial conditions of the orbit or its specific location in phase space. These results support the intuitive higher-dimensional understanding of a negative cosmological constant as a parameter that provides stronger gravitational binding, leading to an increase of the integrability of the dynamics and thus an increase in the stability of trajectories.

Remarkably, as ${\Lambda }_e\to {\Lambda }_{negcrit}$ the chaotic regions nearly vanish. Since the area of the chaotic regions in the Poincaré section were found to be roughly proportional to $\left|\frac{{\Lambda }_e}{{\Lambda }_{negcrit}}\right|$ for the range of possible energies that could be numerically investigated, it has been conjectured that this holds for arbitrarily large values of H. Conversely, the area of the chaotic regions in the Poincaré section increases as ${\Lambda }_e$ becomes increasingly positive, shown in the lower right panel of Figure 43. This occurs within the regions corresponding to the pretzel orbits and in the regions between annulus and pretzel orbits. This phenomenon has likewise been conjectured to occur at all energies [118].

6.6.2. Charged Configurations

The charged 3-body case allows study of additional novel phenomena such as localized vacuum energy and the breaking the full hexagonal symmetry. Concerning the former, since the electromagnetic coupling between any pair of charges induces a vacuum energy between them, we can study how this localization of vacuum energy modifies effects with ${\Lambda }_e\ne 0$. Concerning the latter, the shape of the hexagonal potential becomes elongated unless all particles have the same charge, as in the unequal mass case [117].

Diagrams of representative cases for the potential (relative to the total rest mass) are shown in Figure 45 and Figure 46, taking (for simplicity) all charges to be equal in magnitude.

When all particles have identical charge the potential has hexagonal symmetry, but for different charges this symmetry becomes skewed. As shown in the first and third columns in Figure 45, the potential is stretched along the $\lambda =0$ axis in the (+ + ne) case, corresponding to the decrease in the electric potential as particles 1 and 2 separate. This results in an increase in the magnitude of $\rho$. As the magnitude of the charge increases the width of the potential at lower energies likewise increases and the sides of the hexagonal cross section become more concave. The value of $V_{Rc}$ at which the cross section of V is largest is also reduced. The same effects occur for the (+ 0 +) and (0 + +) configurations, but with the potential compressed along the $\rho +\sqrt{3}\lambda =0$ and $\rho -\sqrt{3}\lambda =0$ axis respectively.

By contrast in the (+ - 0) case, shown in the second and fourth columns in Figure 45, the potential is compressed along the $\lambda =0$ axis. This corresponds to a decrease in electric potential with increasing $\rho$.

When two of the particles have positive charges and the third negative (+ + -), shown in the first and third columns in Figure 46, the hexagon becomes elongated in the $\rho$ direction as the magnitude of the charge increases. If all three particles have equal positive charge (+ + +), hexagonal symmetry is preserved as shown in the second and fourth columns in Figure 46. The width of the potential increases at lower energies and its sides become more convex as the magnitude of the charge increases.

As before, annulus, pretzel, and chaotic trajectories are present depending on the initial conditions; samples are shown in Figure 47. We can gain more insight by considering Poincaré plots, shown in Figure 48 (+ + 0), Figure 49 (+ - 0), and Figure 50 (+ + -).

If one particle is neutral, only the product of the charges is relevant and not their individual magnitudes; their relative sign determines the sign of V from (191). The region between the two charged particles has constant vacuum energy throughout the motion, but the size of this region changes as the particles move. Consequently using equal magnitude charges has no loss of generality.

For the (+ + 0) case in Figure 48 there are two different Poincaré plots: one corresponding to the crossing of the two identical charged particles and the other corresponding to the neutral particle crossing either of the identical particles. Comparing to the lower right panel of Figure 48, we see that chaotic behaviour emerges more rapidly with increasing repulsive charge than increasing energy, with the chaotic regions much more widespread in Figure 48. The chaotic behaviour is notably enhanced, filling the pretzel regions and pushing the annulus region in the opposite direction as the charge increases, leaving almost no circular periodic motions.

For three bodies of distinct charge, the relative magnitudes of the charges do matter, since the magnitude of the vacuum energy between them now changes as they interchange positions. The (+ - 0) case depicted in Figure 49 is very different from that induced by a negative cosmological constant, shown in the lower left panel of Figure 43, where the amount of chaos is much less as compared to the non-relativistic case. In Figure 49 there is an increase in chaos throughout a band between the pretzel and annulus regions, even for comparatively small values of the total energy. As the energy increases (right panel) the expansion of the pretzel areas pushes out the annulus region. Some of the chaotic areas between the quasi annulus and pretzel regions transform into the substructures of the repeating circles.

The (+ + -) case in Figure 50 is yet again different. The recurring circles stay around the annulus region even for higher charges, and do not change into connected areas, though the annulus region is pushed to dissolve into the left pretzel region in the lower graphs and the lower pretzel in the upper graphs with a thin band of chaos between the two regions. The magnitude of the vacuum energy between the positive charges changes between 0 and $2e^2$, the former occuring whenever the negatively charged particle is between the other two. This leads to more widespread chaotic behaviour comparted to the (+ - ne) configuration, due to repeated changes in vacuum energy for a given total energy.

The behaviour of the charged case is quite rich and varied. A number of other scenarios have been studied at low energies [125], but exploration of high-energy behaviour has yet to be carried out.

7. The 4-Body Problem

The N particle OGS can be mapped to a single particle moving in $N-1$ dimensions in a linear potential whose equipotential surfaces are that of an $N-1$ simplex. Since the largest number of spatial dimensions that can be directly visualized is three, the $N=4$ system – the 4-body OGS – is of particular interest. Curiously, it has received almost no attention – only the non-relativistic system has been studied [126]. The relativistic 4-body system has yet to be investigated.

7.1. 4-Body Potential

The Hamiltonian for the non-relativistic 4-body problem is given by (1) with $N=4$. There are now six independent degrees of freedom: the three separations between the particles and their conjugate momenta. Writing

$$\begin{array}{cc}{z}_{12}=\sqrt{2}\rho & z_{34}=\sqrt{2}\alpha \\ z_{13}=\frac{1}{\sqrt{2}}(\rho +\sqrt{3}\beta -\alpha )& z_{23}=\frac{1}{\sqrt{2}}(-\rho +\sqrt{3}\beta -\alpha )\\ z_{24}=\frac{1}{\sqrt{2}}(-\rho +\sqrt{3}\beta +\alpha )& z_{14}=\frac{1}{\sqrt{2}}(\rho +\sqrt{3}\beta +\alpha )\end{array}$$

(200)

where $z_{ij}=z_i-z_j$, the conjugate momenta are

$$\begin{array}{cc}{p}_1=\frac{1}{\sqrt{2}}(p_{\rho }+\frac{\sqrt{3}}{2}{p}_{\beta })& p_2=\frac{1}{\sqrt{2}}(-p_{\rho }+\frac{\sqrt{3}}{2}{p}_{\beta })\\ p_3=\frac{1}{\sqrt{2}}(p_{\alpha }-\frac{\sqrt{3}}{2}{p}_{\beta })& p_4=\frac{1}{\sqrt{2}}(-p_{\alpha }-\frac{\sqrt{3}}{2}{p}_{\beta })\end{array}$$

(201)

with conservation of momentum allowing us to set $p_1+p_2+p_3+p_4=0$; the centre of mass can be fixed at the origin without loss of generality. When one of $z_{12}$, $z_{23}$, or $z_{13}$ vanish (two particles are placed directly on top of one another) this reduces to the 3-body case studied in the previous section.

In the equal mass case the Hamiltonian (1) becomes

$$\begin{array}{cc}\hfill H& =\frac{1}{2m}(p_{\rho }^2+p_{\alpha }^2+\frac{3}{2}{p}_{\beta }^2)+\frac{8\pi G{m}^2}{\sqrt{8}}\left[\left|\rho \right|+\left|\alpha \right|+\frac{1}{2}\left|\rho +\alpha +\sqrt{3}\beta \right|\right.\hfill \\ \hfill & \left.+\frac{1}{2}\left|\rho -\alpha +\sqrt{3}\beta \right|+\frac{1}{2}\left|\rho +\alpha -\sqrt{3}\beta \right|+\frac{1}{2}\left|\rho -\alpha -\sqrt{3}\beta \right|\right]\hfill \end{array}$$

(202)

which is the Hamiltonian of a single particle (the box-particle) moving in three spatial dimensions in a linear potential whose shape is that of a 3-simplex.

The potential

$$V(\rho ,\beta ,\alpha )=H(p_{\rho }=0,p_{\beta }=0,p_{\alpha }=0)$$

(203)

has equipotential surfaces which are that of a cube of pyramid-shaped sides, shown in Figure 51. As V increases the simplex likewise increases, as is clear from the comparing the two diagrams. A cross-section of this surface through any of the edges of one of these pyramids yields a hexagon with sides of unequal length. For such cross-sections the system reduces to that of 3-body case with unequal masses since two particles will occupy the same position.

7.2. Motion Classification

The 4-body system has an interesting structure that can be described in terms of braid operators. This generalizes the A-type and B-type motions of the 3-body case.

As in the 2-body and 3-body cases, each particle moves with a constant acceleration that is proportional to the difference between the total mass on its right and left sides prior to any collision. Assuming the particles move through each other, after a collision, the mass difference experienced by any given particle will in general change, and consequently the accelerations of the particles also changes. From the viewpoint of the box particle, any crossing of a pair of particles corresponds to the box particle crossing a plane bisecting the 3-simplex through its vertices and edges. There are a total of six such planes, obtained by setting any one of the six quantities in (200) to zero. The planes occur in pairs whose line of intersection is along each of the three principal axes.

Any sequence of crossings of N bodies can be described using braid Group notation [127] using the set $\{{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_{N-1}\}$, with ${\sigma }_j={\sigma }_j^{-1}$ since crossing direction is irrelevant. The positions of the particles (and not the particle themselves) are are ordered as $(1,2,3,...,N)$, where the left-most particle is at position 1, next 2, and so on with the right-most particle being at position N. Note that it is the particular sequence of collisions that is important; any permutation of the operators would result in loss of information about the motion in the system.

Applying this to the 3-body system, the braid operators are $\{{\sigma }_1,{\sigma }_2\}$ and the motion can be classified into

$$\begin{array}{cc}{\sigma }_1{\sigma }_1,{\sigma }_2{\sigma }_2& \mathrm{A}\mathrm{motion}\\ {\sigma }_1{\sigma }_2,{\sigma }_2{\sigma }_1& \mathrm{B}\mathrm{motion}\end{array}$$

(204)

or in other words the only interesting types of motion are when the same pair of particles crosses twice in a row (A-motion) or when one particle crosses each of its compatriots in succession (B-motion).

In the 4-body case, there are only 3 possible crossings – $\{z_{12},z_{23},z_{34}\}$ – at any given instant, and the braid operators $\{{\sigma }_1,{\sigma }_2,{\sigma }_3\}$, respectively correspond to an interchange between the right-most, middle, and left-most pair of bodies. In this case we have

$$\begin{array}{cc}{\sigma }_1{\sigma }_1,{\sigma }_2{\sigma }_2,{\sigma }_3{\sigma }_3& \mathrm{A}\mathrm{motion}\\ {\sigma }_1{\sigma }_2,{\sigma }_2{\sigma }_1,{\sigma }_2{\sigma }_3,{\sigma }_3{\sigma }_2& \mathrm{B}\mathrm{motion}\\ {\sigma }_1{\sigma }_3,{\sigma }_3{\sigma }_1& \mathrm{C}\mathrm{motion}\end{array}$$

(205)

which is depicted in Figure 52. The A and B motions represent the same physical situations as in the 3-body case, but the C motion is new: two particles cross one another and then the other two cross one another.

Multiple collisions can occur. In the 3-body system this happens when the hex particle crosses the origin, corresponding to all three bodies meeting at the same point at some instant of time. The analog of this in the 4-body case occurs when the box particle crosses the line of intersection of any two bisecting planes of the 3-simplexes. There are two kinds of these 3-body collisions described by $\{{\sigma }_{1^3},{\sigma }_{2^3}\}$. There are also two kinds of 4-body collisions described by $\{{\sigma }_{1^4},{\sigma }_{1^2{3}^2}\}$. The former corresponds to all 4 bodies meeting at a single point, equivalent to the box particle crossing the origin. The latter occurs when one pair of particles cross at one point and the other pair cross at a different point at the the same time, corresponding to the box particle crossing one of the three lines connecting opposite vertices of the pyramids in the simplex (see Figure 51).

7.3. Equal Mass Trajectories

It can be straightforwardly shown that if one of the box-particle’s position and momentum coordinates are initially zero, they will remain zero throughout the motion, and that all the phenomena seen in the 3-body case in the previous section are recovered [126]. The more interesting situation is when the box particle exhibits motion in all spatial directions.

A comparison is shown in Figure 53, where the $\alpha ,p_{\alpha }$ remains fixed in the upper plots, but in the lower plots either $\alpha$ (lower left) or $p_{\alpha }$ (lower right) deviates from zero. The motion in the $\alpha$ direction simply perturbs the patterns in the upper figures, effectively giving a “thickness” to the original hex-particle patterns.

By carefully choosing the initial conditions it is possible to obtain genuinely novel periodic orbits in three dimensions, as shown in Figure 54. This trajectory has a pretzel form when projected onto two of the planes $(\rho ,\beta )$ (upper left) and $(\rho ,\alpha )$ (not shown) and has annulus form when projected onto the third $(\beta ,\alpha )$ plane (upper right). The full three-dimensional orbit is shown in the lower left panel (with no perspective so that lines further away do not appear smaller). The motion here is $\overline{{\left(C{B}^{2}C{B}^{2}C{B}^6\right)}^{2}C{B}^6}$ and has no analogue in the 3-body system. Motions of each particle are shown in the lower right panel, and have two bodies undergoing small amplitude oscillations (solid and dotted lines) with the other two undergoing larger amplitude oscillation (dash and dot-dash lines).

The upper panels of Figure 55 show a situitation where the four particles begin at zero momentum and are equally spaced. As expected they are all attracted together and cross at the same point, repeating this pattern indefinitely. The symbol sequence is undefined since all four particles always ‘collide’ at the same time step. The box particle oscillates along a line in $(\rho ,\beta ,\alpha )$ space (upper left), with all particles crossing at the origin simultaneously. The outer two particles (solid, dot-dash) undergo large amplitude oscillations and the inner two (dash, dotted) undergo small amplitude oscillations. However the motion is unstable: a slight change in any of the initial conditions (via either a small perturbation in position or momentum) throws the system into chaos. This is shown in the bottom two panels of Figure 55 where the initial value of $\rho$ is slightly increased. All particle trajectories (lower right) continually vary their oscillation amplitudes. This is evident within 30 time steps, where we see the dashed line grow in amplitude whilst the dot-dash one shrinks.

7.4. Poincaré Plots

The most natural extension of a Poincaré plot to the 4-body case is use spherical coordinates, plotting the radial momentum $p_R$ against the squares of the two angular momenta $p_{\theta }^2$ and $p_{\varphi }^2$. Writing

$$\begin{array}{ccc}\hfill sin\varphi & =& \frac{\beta }{\sqrt{{\rho }^2+{\beta }^2}},cos\varphi =\frac{\rho }{\sqrt{{\rho }^2+{\beta }^2}}\hfill \\ \hfill sin\theta & =& \frac{\sqrt{{\rho }^2+{\beta }^2}}{\sqrt{{\rho }^2+{\beta }^2+{\alpha }^2}},cos\theta =\frac{\alpha }{\sqrt{{\rho }^2+{\beta }^2+{\alpha }^2}}\hfill \end{array}$$

(206)

R is the distance from the origin to the point of crossing in $(\rho ,\beta ,\alpha )$ space and $(\theta ,\varphi )$ are the polar and azimuthal angles of this point. The unit vectors for these spherical coordinates are

$$\widehat{R}=\left[\begin{array}{c}cos\varphi sin\theta \\ sin\varphi sin\theta \\ cos\theta \end{array}\right],\widehat{\varphi }=\left[\begin{array}{c}-sin\varphi \\ cos\varphi \\ 0\end{array}\right],\widehat{\theta }=\left[\begin{array}{c}cos\varphi cos\theta \\ sin\varphi cos\theta \\ -sin\theta \end{array}\right]$$

(207)

The associated momenta are

$$\begin{array}{ccc}\hfill p_R& =& \frac{p_{\rho }\rho +p_{\beta }\beta +p_{\alpha }\alpha }{\sqrt{{\rho }^2+{\beta }^2+{\alpha }^2}}=\widehat{R}·\overrightarrow{p}\hfill \\ \hfill p_{\varphi }& =& \frac{-p_{\rho }\beta +p_{\beta }\rho }{\sqrt{{\rho }^2+{\beta }^2}}=\widehat{\varphi }·\overrightarrow{p}\hfill \\ \hfill p_{\theta }& =& \frac{p_{\rho }\rho \alpha +p_{\beta }\beta \alpha -p_{\alpha }({\rho }^2+{\beta }^2)}{\sqrt{({\rho }^2+{\beta }^2+{\alpha }^2)({\rho }^2+{\beta }^2)}}=\widehat{\theta }·\overrightarrow{p}\hfill \end{array}$$

(208)

and (208) can be used to compute $p_R$ , $p_{\varphi }^2$ and $p_{\theta }^2$ whenever two of the four particles cross one another.

Constructing complete Poincaré plots is a challenge for two reasons. First, the standard approach of choosing a range of initial conditions that fill in the important regions is computationally much more formidable since with five independent variables the number of possible plots is very large. This can be dealt with by automating the generation of data over a specific range of initial conditions, but a significant reduction in the number of time steps must be employed for computational tractability. Visualizing the large number of discrete points in 3-dimensional space is the second challenge. The plots in Figure 56 were constructed by separating out the space into millions of minute three-dimensional boxes, assigning a value corresponding to the number of Poincare points that fall inside and a position corresponding to the location of the box. Although this limits the ability to zoom in to observe self-similar structures, highly saturated regions of chaos tend to show up well.

Figure 56 shows two special slices: $p_{\varphi }=0$ (the “bottom” slice) and $p_{\theta }^2<0.0005$ (the “side” slice). Out of 12 million points generated, 400,000 are in the bottom slice and 500,000 are in the side slice. In the latter case any trajectory with $p_{\theta }=0$ would remain in a cone rooted at the origin, and so the above bound on $p_{\theta }^2$ imposes the constraint of $p_{\theta }$ being “close” to zero.

The bottom slice in Figure 56 bears a resemblance to the non-relativistic 3-body case shown in Figure 31, exhibiting mixed regions of chaos and integrability. By constrast, the side slice does not display the same patterns and fractal-like properties as the bottom slice. Whether this is due to an insufficient number of time steps, a failure to cover a sufficient range of initial conditions or an intrinsic lack of any patterns is not yet known.

The non-relativistic 4-body problem has a number of other interesting aspects [126]. These include apparently chaotic motion in some projections with quasi-periodic motion in others, novel Poincaré plots for particular classes of orbits, and Lyapunov exponents that asymptote to constant values ranging between ${10}^{-2}$ for chaotic trajectories, ${10}^{-3}-{10}^{-4}$ for quasi-periodic trajectories and ${10}^{-5}$ (the limits of numerical precision) for periodic trajectories [126]. More complete studies remain to be carried out, not only for the non-relativistic system, but for its relativistic counterpart.

8. The N-Body Problem

Studies of the N-body OGS for $N>4$ generally have been concerned with its statistical properties. Many unanswered questions remain despite extensive study. Its ergodic and equipartition properties are still not well understood. Whether or not the OGS can attain a true equilibrium state from arbitrary initial conditions is also not clear. These issues remain outstanding in large part because the attractive interactions are cumulatively long range, unlike typical thermodynamic systems that have repulsive and short range interactions between their constituents.

However some statistical properties of the OGS are known. The single-particle distribution function in both the canonical and microcanonical ensemble has been derived [4]. These distribution functions reduce to the isothermal solution of the Vlasov equation in the large N limit.

Even less is known about the N-body ROGS, with only one study of its statistical properties having been carried out to date [128]. In this section a few of the basic results of this system will be summarized.

8.1. Motion Classification

Braid operators $\{{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_N\}$ can be used to classify motion in the N-body system. A sequence of m pair crossings will be described by

$${\sigma }_{f\left(1\right)}{\sigma }_{f\left(2\right)}{\sigma }_{f\left(3\right)}...{\sigma }_{f\left(m\right)}$$

(209)

where $1\le f\left(x\right)\le (N-1)$ for all $1\le x\le m$ is a discrete integer function and ${\sigma }_{f\left(x\right)}$ means that the bodies currently in positions $f\left(x\right)$ and $f\left(x\right)+1$ cross. Crossing directions are irrelevant, and for a given trajectory any given sequence of m braid operators forms a unique ordered list of crossings.

It is possible to define a a metric that describes the relative “distance” between any pair of crossings via

$$g\left(x\right)\equiv |\Delta f(x)|=|f(x+1)-f(x)|$$

(210)

which implies $0\le g\left(x\right)\le (N-2)$ for all $1\le x\le (m-1)$. The motion can be then classified as follows:

$$\begin{array}{cc}g\left(x\right)& \mathrm{Motion}\mathrm{Class}\\ 0& A\\ 1& B\\ 2& C\\ 3& D\\ ...\end{array}$$

(211)

denoting each type by increasing letters of the alphabet. A-motion corresponds to any 2 crossings in nearest proximity – two particles cross each other twice in succession. B-motion corresponds to any 2 crossings in next-nearest proximity – two particles cross each other, and then one of them crosses its other nearest neighbour. C-motion corresponds to any 2 crossings in next-to-next-nearest proximity: two particles cross each other and then a neighbouring pair cross each other. This continues until the left-most pair of particles cross one another followed by the crossing of the right-most pair (or vice-versa), which is the extreme case. In the 4-body case, for example, ${\sigma }_1{\sigma }_2{\sigma }_1{\sigma }_3{\sigma }_2$, yields from (205) the symbol sequence $BBCB$. Computing successive values of $g\left(x\right)$ ( $1,1,2,1$) gives same result.

A collision of m particles simultaneously corresponds to a single particle crossing through an $(N-m)$-dimensional surface in the interior of the $(N-1)$ simplex. Such collisions can be further classified by extending the braid group notation with the set $\{{\sigma }_{1^m},{\sigma }_{2^m},\dots ,{\sigma }_{{(N+1-m)}^m}\}$. The subscript denotes which set of particles is involved, beginning with the left-most; the superscripts denote the number of particles in the collision, with the superscript “2” dropped when pairwise collisions occur. For example ${\sigma }_{8^5}$ denotes a 5-particle collision that involves particles 8-12. All collisions yield crossings with the exception of the initial conditions causing m particles to occupy the same point throughout the motion (as in the upper left panel of Figure 55 in the 4-body case). In this situation the system reduces to that of an (unequal mass) $\left(N-m\right)$-body problem.

After a multiple particle collision it is always possible to predict the new order of particles given their preceding order. Defining rightward velocity as positive, for two adjacent particles a small time just before the collision, the left one must have a larger velocity than right one or else the latter would be moving away and not toward the left one and no collision would occur. Applying this reasoning to every adjacent pair implies that on moving from left to right the velocity of each particle decreases in the sequence, with the left-most particle having the largest velocity. Immediately after the collision the original order will be reversed, since the previously left-most particle will be travelling rightward faster than all other particles in the collision and emerge afterwards as the right-most particle, and so on for all particles in the collision. If any of the n particles does not satisfy the increasing velocity condition there will be fewer than n particles in the multiple collision.

8.2. Post-Newtonian Canonical Ensemble

The canonical one-particle distribution function can be shown to be

$$f_c^R(p,z)=\frac{1}{\mathcal{Z}N!}\int \int d\mathbf{p}d\mathbf{z}\delta \left(\overline{p}\right)\delta \left(\overline{z}\right)exp\left(-\beta H\right)N^{-1}\sum _a\delta \left(z-z_a\right)\delta \left(p-p_a\right)$$

(212)

by making use of momentum conservation and translation invariance [128], where H is the Hamiltonian of the system. This quantity is straightforward to compute for the N system [4], since its Hamiltonian (1) is at most quadratic in the canonical variables. However the Hamiltonian (76) for the R system is a highly non-linear function of these variables, and computation of (212) is not obvious. For practical reasons, the post-Newtonian Hamiltonian $H_{pN}$ in (28) has been used to gain insight into the statistical properties of the ROGS [128].

The quantity

$$\begin{array}{cc}\hfill \mathcal{Z}& =\frac{1}{N!}\int \int d\mathbf{p}d\mathbf{z}\delta \left(\overline{p}\right)\delta \left(\overline{z}\right)exp\left(-\beta H_{pN}\right)\hfill \end{array}$$

(213)

is the partition function. A somewhat tedious calculation yields the result [128]

$$\begin{array}{c}\hfill \mathcal{Z}=\frac{exp\left[-\beta M{c}^2-\frac{3\left(N-1\right)}{2}ln\left(\beta m{c}^2\right)-\left\{\frac{\left(5N+3\right)\left(N-1\right)+8N{\sum }_{k=1}^{N-1}{\sum }_{l=k+1}^{N-1}\frac{\left(l-k\right)}{l\left(N-k\right)}}{8N\beta m{c}^2}\right\}\right]}{\sqrt{N}{\left(\sqrt{2\pi }G/c^3\right)}^{\left(N-1\right)}{\left[\left(N-1\right)!\right]}^2}\end{array}$$

(214)

to lowest relativistic order, where $M={\sum }_{a=1}^N{m}_a$. The average energy is then

$$\begin{array}{ccc}\hfill 〈E〉& =& -\frac{\partial }{\partial \beta }ln\mathcal{Z}\hfill \\ & =& M{c}^2+\frac{3}{2\beta }\left(N-1\right)-\frac{\left(5N+3\right)\left(N-1\right)+8N{\sum }_{k=1}^{N-1}{\sum }_{l=k+1}^{N-1}\frac{\left(l-k\right)}{l\left(N-k\right)}}{8{\beta }^{2}M{c}^2}\hfill \end{array}$$

(215)

to the relevant order in $c^{-2}$. For fixed $M=Nm$ the relativistic correction grows quadratically with N and is negative. Consequently average energy of the ROGS is lower than its non-relativistic counterpart at the same temperature.

When the thermal energy $kT={\beta }^{-1}$ is sufficiently small relative to the rest energy $M{c}^2$ of the system, the average energy $\zeta \equiv \frac{\left(E-M{c}^2\right)}{M{c}^2}$ is not much different from its non-relativistic value of $\frac{3(N-1)}{2\beta M{c}^2}$. As $\beta$ decreases the value of $〈\zeta 〉$ increases more slowly than its non-relativistic counterpart, and reaches a maximum at

$$\beta ={\beta }_{\max }\equiv \frac{\left(5N+3\right)\left(N-1\right)+8N{\sum }_{k=1}^{N-1}{\sum }_{l=k+1}^{N-1}\frac{\left(l-k\right)}{l\left(N-k\right)}}{6\left(N-1\right)M{c}^2}\underset{N>>1}{\to }\left(\frac{7}{2}-\frac{2{\pi }^2}{9}\right)\frac{N}{M{c}^2}$$

(216)

which is half the value of its non-relativistic counterpart. This maximum value is plotted in Figure 57 as a function of N, and asymptotes to the constant value $〈\zeta 〉=.573940872$ as $N\to \infty$. For $\beta >{\beta }_{\max }$ the average energy $\zeta$ decreases with increasing $\beta$, vanishing at $\beta =\frac{1}{2}{\beta }_{\max }$. The post-Newtonian expansion (28) breaks down well before this value of $\beta$.

The canonical momentum distribution function can be obtained by integrating the single particle distribution function. The result is

$$\begin{array}{ccc}\hfill {\vartheta }_{cn}\left(p\right)& =& {\int }_{-\infty }^{\infty }dz{f}_{cn}(p,z)\hfill \\ & =& \sqrt{\frac{\left(N\beta \right)}{2\pi m\left(N-1\right)}}exp\left[-\frac{N\beta p^2}{2m\left(N-1\right)}\right]\hfill \\ & & \times \left(1+\frac{1}{\beta m{c}^2}\left(\frac{N{\beta }^2{p}^4\left(N^2-3N+3\right)}{8m^2{\left(N-1\right)}^3}-\frac{\beta p^2\left(4N^2-7N+6\right)}{4m{\left(N-1\right)}^2}+\frac{5N(N-1)+3}{8N\left(N-1\right)}\right)\right)\hfill \end{array}$$

(217)

and corrects the standard non-relativistic Gaussian expression by a polynomial in $p^2$.

The preceding expression can be rewritten in terms of dimensionless variables

$$\eta \equiv \frac{p}{mV}{V}^2\equiv \frac{4\left(E-M{c}^2\right)}{3M}=\frac{4\zeta c^2}{3}$$

(218)

and is plotted in Figure 58.

The central momentum density grows with increasing $\zeta$, but falls off more rapidly than its non-relativistic counterpart does. However for $\eta >2$, the momentum density grows exponentially relative to its non-relativistic counterpart, overtaking it for sufficiently large $\eta$. This is clearly seen in the right-hand panels of Figure 58. The differences become less pronounced as N increases, although the basic features remain the same for all N.

8.3. Other Statistical Features

Computation of other quantities, such as the canonical density distribution, the full single particle canonical distribution, and the microcanonical distribution functions involve a considerable amount of tedious algebra. Some general features emerge from this analysis [128]. One is that relativistic effects cool the system: at a given energy, the ROGS temperature is smaller than the OGS temperature. Another is that ROGS density and distribution functions become more sharply peaked than their OGS counterparts as $\zeta$ increases for any given N. For sufficiently large values of the position parameter z, ${\rho }_{\mathrm{OGS}}>{\rho }_{\mathrm{ROGS}}$ for both the canonical and microcanonical distribution functions. The momentum densities exhibit different behaviour, with ${\vartheta }_{\mathrm{OGS}}>{\vartheta }_{\mathrm{ROGS}}$ for small values to intermediate values $\eta$, but for large enough $\eta$, the inequality is reversed and ${\vartheta }_{\mathrm{OGS}}<{\vartheta }_{\mathrm{ROGS}}$.

Further exploration of the N-body ROGS will be challenging. A natural first step would be a consideration of the charged and cosmological systems at the post-Newtonian level. The unequal mass case is likely similarly tractable, though it will be considerably more difficult. However a full understanding almost certainly require a numerical approach, particularly to go beyond leading order corrections in $1/c^2$.

9. The Circular N-Body Problem

Thus far the discussion has been concerned with lineal topology. For non-relativistic gravity this is the only option. The OGS equations (4) for a general potential V are

$$\begin{array}{ccc}\hfill V^{\prime \prime }& =& 4\pi Gm\delta (x-z(t\left)\right)\hfill \end{array}$$

(219)

$$\begin{array}{ccc}\hfill \dot{p}& =& -V^{\prime }\left(z\right)\hfill \end{array}$$

(220)

$$\begin{array}{ccc}\hfill \dot{z}& =& \frac{p}{m}\hfill \end{array}$$

(221)

for a single particle. The first equation implies $V=2\pi Gm\left|x\right|$, which has a vanishing derivative at $x=0$ and yielding $p=z=0$ as a consistent single particle solution. However if the topology is circular then both $V\left(L\right)=V(-L)$ and $V^{\prime }\left(L\right)=V^{\prime }(-L)$ for some L, where $2L$ is the circumference of the circle. These matching conditions have no solution unless another point source of negative mass is introduced. For N bodies modelled as compact smeared sources the problem remains: the potential grows linearly with increasing distance from the source(s) and the matching conditions cannot be satisfied for physically reasonable (i.e., positive mass) sources.

However the ROGS, being in a dynamical spacetime, does not suffer from this problem, since the spacetime can expand or contract in response to the presence of sources. It is possible to solve the canonical equations of motion and obtain both single particle [129] and exact N-body static equilibrium solutions [130,131]. This latter solution corresponds to a spacetime that expands/contracts in response to N equal mass bodies at equidistant proper separations from one another. These are the first N-body dynamic equilibrium solutions in any relativistic theory of gravity.

The action is still given by (10) but the extrinsic curvature K in (17) is now taken to be a time coordinate $\tau \left(t\right)$. This allows the elimination of $\pi$ from all canonical field equations. The reduced canonical action (22) instead becomes

$$I_R=\int d^{2}x\left\{\sum _a{p}_a{\dot{z}}_a\delta (x-z_a)+\Pi \frac{\partial }{\partial t}\left(\Psi +ln\gamma \right)-\mathcal{H}\right\},$$

(222)

where now

$$H=\int dx\mathcal{H}=\frac{2\dot{\tau }}{\kappa }\int dx\sqrt{\gamma }$$

(223)

which is the circumference functional of the circle when $\dot{\tau }$ is constant. The N-body system is now time-dependent. The spatial metric can be chosen so that $\gamma =\gamma \left(t\right)$, as in general relativity on $(2+1)$-dimensional spatially compact manifolds [132]. The Hamiltonian will be time-independent if a time parametrization is chosen so that $\dot{\tau }\sqrt{\gamma }$ is constant.

Equilibrium solutions correspond to a situation in which all particles are motionless at various points around the circle. Consequently they are characterized by $\dot{z_a}=0=p_a$, and so the canonical field equations (49) – (54) become [130,131]

$$\begin{array}{ccc}\hfill 2{\Pi }^{\prime }-\Pi {\Psi }^{\prime }=0;& & \hfill \end{array}$$

(224)

$$\begin{array}{ccc}\hfill {\Psi }^{\prime \prime }-\frac{1}{4}{\left({\Psi }^{\prime }\right)}^2-{(\kappa \Pi )}^2+\gamma \left({\tau }^2-{\Lambda }_e\right)+\kappa \sum _a\sqrt{\gamma }m\delta (x-z_a\left(x^0\right))=0& & \hfill \end{array}$$

(225)

$$\begin{array}{ccc}\hfill N_0^{\prime \prime }-\dot{\tau }\gamma -N_0\left\{({\tau }^2-{\Lambda }_e)\gamma +\kappa /2\sum _a\sqrt{\gamma }m\delta (x-z_a\left(x^0\right))\right\}=0& & \hfill \end{array}$$

(226)

$$\begin{array}{ccc}\hfill N_1^{\prime }-\dot{\gamma }/2+\gamma \tau N_0=0& & \hfill \end{array}$$

(227)

$$\begin{array}{ccc}\hfill \dot{\Pi }+{\partial }_1(-\frac{1}{\gamma }{N}_1\Pi +\frac{1}{2\kappa \sqrt{\gamma }}{N}_0{\Psi }^{\prime }+\frac{1}{\kappa \sqrt{\gamma }}{N}_0^{\prime })=0& & \hfill \end{array}$$

(228)

$$\begin{array}{ccc}\hfill \dot{\Psi }+2N_0\left(\kappa \frac{\Pi }{\sqrt{\gamma }}+\tau \right)-N_1(\frac{1}{\gamma }{\Psi }^{\prime })=0& & \hfill \end{array}$$

(229)

with

$${\left.\frac{\partial N_0}{\partial x}\right|}_{x=z_a}=0N_1\left(z_a\right)=0$$

(230)

which follows from the geodesic equations.

Solutions to (224–229) must be appropriately matched at the locations $z_a$ of each particle and the identification point $\left|x\right|=L$. The solutions for $\Psi$ and $N_0$ must be continuous but are not differentiable at the particle locations; for example ${lim}_{ϵ\to 0}\left[{\Psi }^{\prime }(z_a+ϵ)-{\Psi }^{\prime }(z_a-ϵ)\right]=\kappa \sqrt{\gamma }m$. Setting $r=\left|z_a-z_{a+1}\right|=\frac{2L}{N}$ corresponding to N bodies of equal mass m at equal time-varying proper separations, the solution is

$$\sqrt{\gamma }=\frac{N}{L\sqrt{{\tau }^2-{\Lambda }_e}}\mathrm{arctanh}\left(\frac{\xi \sqrt{\frac{{\kappa }^2{M}^2}{16N^2}+({\xi }^2-1)({\tau }^2-{\Lambda }_e)}-\frac{\kappa M}{4N}}{\frac{{\kappa }^2{M}^2}{16N^2}+{\xi }^2({\tau }^2-{\Lambda }_e)}\sqrt{{\tau }^2-{\Lambda }_e}\right)$$

(231)

for the metric function $\gamma \left(t\right)$, where $\xi$ is an integration constant, ${\widehat{c}}^2\equiv \gamma \left({\tau }^2-{\Lambda }_e\right)=c_{+}^2>0$, and $M=mN$ is the total mass of the system.

The remaining functions are

$$\begin{array}{ccc}\hfill \kappa \Pi & =& \frac{\pm c_{+}\sqrt{{\beta }^2-1}}{\frac{sinh\left(\frac{c_{+}L}{N}\right)}{sinh\left(c_{+}L\right)}{\sum }_{a=1}^{n}cosh(c_{+}(|x-z_a|-L))-\beta }\hfill \end{array}$$

(232)

$$\begin{array}{ccc}\hfill \Psi & =& -2ln\left(\beta -\frac{sinh\left(\frac{c_{+}L}{N}\right)}{sinh\left(c_{+}L\right)}\sum _{a=1}^{n}cosh(c_{+}(|x-z_a|-L))\right)-2ln(\frac{\kappa {\Pi }_0\left(t\right)}{c_{+}\sqrt{{\beta }^2-1}})\hfill \end{array}$$

(233)

$$\begin{array}{ccc}\hfill N_0& =& \frac{\dot{\tau }\gamma }{c_{+}^2\beta }\left(\frac{sinh\left(\frac{c_{+}L}{N}\right)}{sinh\left(c_{+}L\right)}\sum _{a=1}^{n}cosh(c_{+}(|x-z_a|-L)\left)\right)-\beta \right)\hfill \end{array}$$

(234)

$$\begin{array}{ccc}\hfill N_1& =& \gamma \frac{\dot{c_{+}}}{c_{+}}x-\frac{{\gamma }^2\tau \dot{\tau }}{c_{+}^3\beta }\frac{sinh\left(\frac{c_{+}L}{N}\right)}{sinh\left(c_{+}L\right)}\sum _{a=1}^nϵ\left(x-z_a\right)\left[sinh(c_{+}(|x-z_a|-L))+sinh\left(c_{+}L\right)\right]\hfill \end{array}$$

(235)

where

$$\begin{array}{ccc}\hfill \beta & =& \frac{4N}{\kappa M}\sqrt{\frac{{\kappa }^2{M}^2}{16N^2}+({\xi }^2-1)({\tau }^2-{\Lambda }_e)}\hfill \end{array}$$

(236)

and $ϵ\left(x\right)=\frac{\left|x\right|}{x}$ is a step function with $ϵ\left(0\right)=0$.

The system is cyclically symmetric, with $N_1\left(L\right)=N_1(-L)$, and so we can choose the origin to be halfway between any two particles in the $N=$even case or on a particle in the $N=$odd case. The spatial periodicity of the solution can be better seen using the relation

$$\sum _{a=1}^{n}cosh(c_{+}(|x-z_a|-L))=cosh\left(c_{+}f\left(x\right)\right)\frac{sinh\left(c_{+}L\right)}{sinh\left(\frac{c_{+}L}{N}\right)}$$

(237)

where $f\left(x\right)$ is the saw-tooth function that peaks with a value of $L/N$ (in other words, $f\left(z_a\right)=L/N$) at the particle locations and vanishes half-way between the particles (i.e., $f\left(\frac{z_a-z_{a+1}}{2}\right)=0$ ). A simple shift of the origin in the cosh function and a subsequent manipulation of the sum yields the equivalence.

If ${\widehat{c}}^2=c_{-}^2<0$ the another class of solutions exists with

$$\sqrt{\gamma }=\frac{N}{L\sqrt{{\Lambda }_e-{\tau }^2}}\left[arctan\left(\frac{\xi \sqrt{\frac{{\kappa }^2{M}^2}{16N^2}+(1-{\xi }^2)({\Lambda }_e-{\tau }^2)}-\frac{\kappa M}{4N}}{\frac{{\kappa }^2{M}^2}{16N^2}-{\xi }^2({\Lambda }_e-{\tau }^2)}\sqrt{{\Lambda }_e-{\tau }^2}\right)+k\pi \right]$$

(238)

and $c_{+}\to i{c}_{-}$ in (232) – (235). The integration constant $\left|\xi \right|<\sqrt{1+\frac{{\left(\kappa M\right)}^2}{16{\Lambda }_e{N}^2}}$ as a due to the periodicity of $N_1(x,t)$ and k is an integer. No solutions exist for $c^2=0$.

The solution (231–235) (and its $c^2<0$ counterpart) corresponds to an expanding/contracting spacetime of a circle with N bodies at equal time-varying proper separations from one another around the circle. The solution for the spatial metric is equivalent to that of single-particle solution [129] upon rescaling $L\to L/N$ and $M\to M/N$. However this rescaling equivalence does not hold for the remaining functions. For ${\Lambda }_e=0$ the spacetime expands but perpetually decelerates due to the presence of the point masses. If ${\Lambda }_e<0$ the proper circumference of the circle expands from zero to some maximal size and then recontracts. The most interesting behaviour occurs if ${\Lambda }_e>0$. In this case the cosmological expansion opposes decelerating effects due to the point masses and the spacetime can expand from zero size to finite value, evolve from some minimal/maximal circumference to a maximal/minimal size, or undergo perpetual oscillation.

A generalization of this solution to one in which there are an even number of bodies with charge alternating in sign but equal in magnitude also has been obtained [131]. This solution, and its neutral counterpart in (231–235) almost certainly describe an unstable equilibrium, since the masses are all equal and the particles are evenly separated. Perturbations from equilibrium would be interesting to investigate, as they would form model inhomogeneous self-gravitating cosmological systems. This remains an interesting avenue for further study.

10. Conclusions

Relativistic one-dimensional self-gravitating systems provide an interesting (and in the view of this author) undervalued theoretical laboratory for studying a number of physical effects of which rather little is known. These include exact 2-body motion, static N-body equilibrium, relativistic thermodynamics and statistical mechanics, relativistic chaos, and the interplay between gravity, electromagnetism, and cosmological evolution. I shall close this review with a brief overview of seven different research avenues that warrant further study.

• Relativistic Chaos at High Energy
  All studies of relativistic 3-body chaos have been at energies below the cross-sectional maximum of the potential $V_R$ in (177). The regions of chaos in the R system are distortions of their N system counterparts, but do not increase. It would be of great interest to know what the chaotic behaviour is for energies larger than the cross-sectional maximum, where very strong relativistic effects are present. Will the chaotic regions in the Poincaré plots grow or shrink? Such studies would provide further insight into the effects of strong gravity on chaotic systems.
• 4-body Chaos
  The largest value of N in the N-body problem for which equipotential surfaces can be visualized is $N=4$. As noted in Section 7, only the N system has been investigated for its motions and chaotic behaviour. The R system has yet to be investigated along these lines. It is conceivable that qualitatively new features will be observed in this case.
• Fully relativistic Statistical Mechanics
  The scope for exploration here is very broad. Only the statistical properties of the neutral pN system have been studied. The effects of charge and cosmological expansion are not known, and a full study of the statistical properties of the R system remain to be carried out. This latter problem will be quite technically challenging, since the distributions functions cannot be analytically integrated. Some novel blend of analytic and numerical methods will need to be employed.
• Circular N-body Dynamics
  A novel feature of the N-body ROGS is that it admits two distinct spatial topologies: linear and circular. In the latter case there are no dynamical solutions for $N\ge 2$. These will likely need to be found numerically. It would be particularly interesting to investigate 3-body chaos in this setting to see what effects circular topology has on chaotic phenomena.
• The 2-dimensional N-body problem
  Since there is no gravitational radiation in 2 spatial dimensions the N-body problem in this setting is of considerable interest, all the more so since general relativity will provide the foundation for the field equations. This problem has been considered from a topological perspective [133], from which an implicit solution for the metric and the motion of N particles was obtained [134]. The solution becomes explicit for $N=2$. However the relationship between this approach and the canonical approach has only been explored to a limited extent [135,136]. A thorough analysis should be carried out, particularly since particle collisions can form black holes [137] and quite possibly lead to other interesting spacetime effects.
• Extensions to Dilaton Gravity
  The $R=T$ theory (10) has provided the context for exploring the relativistic N-body problem since it is the $D\to 2$ limit of general relativity [27]. However a broad class of 2-dimensional theories of gravity exist [71] and are of physical interest for a variety of reasons. Exploring the N-body problem in this broader context could lead to new physical insights into chaos, relativity, and quantum gravity.
• The Quantum N-body Problem
  The Hamiltonian (76) is the exact energy functional of all degrees of freedom in the relativistic charged 2-body system. Consequently its quantization will be tantamount to fully quantization of gravity coupled to charged matter in one spatial dimension. The N system can be fully quantized, with energy eigenfunctions given in terms of Airy functions, and the associated eigenvalues in terms of their zeroes. Perturbative solutions to the quantum pN system have been obtained [138], but a full analysis of the quantum R system has yet to be carried out. This problem is of considerable interest, since there is experimental evidence that the energy states of neutrons are given by the eigenstates of the N system [139], confirming the test-mass limit of the quantum N-system. A better understanding of the quantum R-system could conceivably lead to experimental tests of relativistic quantum gravity.