New many-body problems in the plane with periodic solutions

In this paper we discuss a family of toy models for many-body interactions including velocity-dependent forces. By generalizing a construction due to Calogero, we obtain a class of N-body problems in the plane which have periodic orbits for a large class of initial conditions. The two- and three-body cases (N=2, 3) are exactly solvable, with all solutions being periodic, and we present their explicit solutions. For N⩾4 Painlevé analysis indicates that the system should not be integrable, and some periodic and non-periodic trajectories are calculated numerically. The construction can be generalized to a broad class of systems, and the mechanism which describes the transition to orbits with higher periods, and eventually to aperiodic or even chaotic orbits, could be present in more realistic models with a mixed phase space. This scenario is different from the onset of chaos by a sequence of Hopf bifurcations.


Introduction
The subject of this paper is a family of toy models for many-body interactions including velocity-dependent forces. We restrict ourselves to the consideration of motions lying in a plane, and by choosing a particular type of interaction we find the surprising feature that for a large class of initial data, all the motions are periodic, while for other initial conditions the trajectories show aperiodic, perhaps even chaotic, behaviour. While the construction of our model N-body systems is based on mathematical considerations rather than physical motivations, we believe that they shed light on the behaviour of more physically realistic models where both regular and irregular motions coexist in the phase space. In particular, the mathematical construction is valid for a very broad class of dynamical systems. This particular family of N-body problems should give a better understanding of systems where there is a transition in phase space from periodic orbits, via orbits of increasing periods, until aperiodic and irregular motion is reached. The mechanism described here is different from other well-known scenarios for the onset of chaos [1]. Although periodic motion is a characteristic of so-called integrable systems, we must emphasize that the models we study are not integrable (although the two-and three-body problems are exactly solvable), and thus are representative of a much wider class of physical models.
The fact that a particular dynamical system possesses many periodic orbits is regarded as being very special, and it is normally related to the existence of a high degree of symmetry. For example, for the motion of a single particle in three-dimensional space under the action of a central potential V(r), the only potentials for which all bounded orbits are closed are the harmonic oscillator V(r) = r 2 and the Kepler potential V(r) = 1/r. It is well known in celestial mechanics that small deviations from the Kepler potential cause the precession of the perihelion of the planets, the orbits being no longer closed ellipses, but rosetta-like figures.
The Kepler potential and the harmonic oscillator are two examples of superintegrable systems. A system with N degrees of freedom is said to be (maximally) superintegrable if 2N − 1 functionally independent first integrals I i (x, p) exist. Since the existence of k independent such quantities causes the flow in phase space to be contained in the submanifold of codimension k determined by the level sets of I i (x, p), i = 1, . . . , k, it is clear that for maximally superintegrable systems (k = 2N − 1) almost all bounded orbits must be periodic [2]. There has been considerable interest in knowing all superintegrable potentials for systems with a low number of degrees of freedom [3,4]. The systematic program consists of classifying all systems that would admit the required number of first integrals, normally assuming a polynomial dependence of the conserved quantities in the momenta. Conversely, a dynamical system, all of whose bounded trajectories are periodic, should have a maximal number of conserved quantities, yet these quantities could be very difficult to construct, or may not even be expressible in terms of elementary functions.
In this direction, a systematic method has been recently proposed by Calogero to modify evolution equations in such a way that they will have periodic solutions for a wide class of initial data, and sometimes for the whole space of initial conditions [5]. This 'trick'(essentially a change of variables) can be applied to a wide variety of equations [6], and interesting examples thereof are a deformation of the classical equations of gravitation [7], systems of coupled non-linear harmonic oscillators in arbitrary dimensions [8] and a many-body rotator problem in the plane [9]. Some of the models constructed in this way are proved to have periodic orbits for arbitrary initial conditions, yet the method does not provide any way of finding the constants of motion. Most often the 'trick' produces dynamical systems that will have periodic orbits only for a certain subset of initial conditions having finite measure in phase space, while apparently non-periodic or possibly chaotic orbits exist for initial conditions outside that set. This work follows closely the results of Cologero et al [9], where a new many-body problem in the plane having many periodic orbits was introduced. The equations of motion of that model describe the motion of N particles subject to pairwise interactions with velocity-dependent forces that vanish asymptotically when the particles are far apart. They are also subject to an external force that could be interpreted as a Lorentz force produced by a constant magnetic field perpendicular to the plane of motion if all the particles had (positive) unitary charge. The model analysed in this paper has the same kind of interactions, plus the addition of an anharmonic external potential.
Another important step in the construction of the model is the relation between the evolution of the time-dependent zeros and coefficients of a polynomial, a technique that has been in use for some time [10] for constructing solvable many-body problems. If the coefficients evolve in time in a controlled way (typically they will satisfy a system of linear ODEs), then the evolution equations for the zeros of the polynomial can be understood as the equations of motion of a many-body problem. This dynamical system is solvable by construction, since the position of the particles at any instant of time t can be calculated by simply finding the roots of a polynomial whose coefficients at time t are known. An extension of this idea has recently been applied to quantum many-body problems [11].
The difference in this work with respect to previous results in this direction [9,12] is that the evolution of the coefficients is non-linear, which gives rise to richer behaviours but makes the analysis considerably more difficult. The fact that these models with complicated non-linear interactions have many periodic orbits is a very remarkable property from the dynamical point of view. However, once the mathematical subtleties in the construction of these models are understood, the periodicity property emerges in a completely natural way. This paper is structured as follows: in section 2 the equations of motion of an N-body problem in the plane are derived by using the relation between the zeros and coefficients of a polynomial, and the complex transformation of time that gives rise to periodic orbits is explained, with the main result summarized in theorem 2.1. In section 3, the non-linear evolution equations for the coefficients of the polynomial are solved for N = 2 and 3, and some special solutions are given for higher N. In those cases where no explicit general solution is available, we performed a Painlevé analysis of the equations in section 4, since it is essentially the analyticity properties of the solutions that determine the periodicity (or otherwise) of the orbits. Finally, we have checked all these predictions by performing a numerical integration of the equations of motion, and we display the results in section 5. The explicit estimates on the initial conditions which ensure periodic orbits are presented in theorem A.1 in the appendix.

Equations of motion
Consider the following monic polynomial in one complex variable ζ, with τ-dependent coefficients where we also take τ to be a complex variable. The previous polynomial can also be written in terms of its zeros as If we impose that the polynomial P(ζ, τ) satisfies the following linear equation: then by inserting (1) and (2) in (3), it can be shown that the coefficients c m (τ) and the zeros ζ j (τ) must satisfy where the primes denote derivatives with respect to τ, and in (4) the convention c m ≡ 0 for m > N holds. The second set of equations (5), after some modifications (see section 2.1), will be considered as the Newtonian equations of motion for an N-body problem in the plane with pairwise interactions and velocity-dependent forces in the presence of an external field. It is clear that this system will be in some sense solvable if one is able to solve the first set of ODEs (4) for the coefficients c m (τ), since the position of the particles ζ j (τ) at any time can then be calculated by a purely algebraic procedure, namely finding the roots of an Nth degree polynomial.
In previous studies [13], the equivalent of (4) for the time evolution of the coefficients c m (τ) was always a system of linear ODEs and therefore always solvable by standard techniques for linear systems. In this case, although the presence of the last term in the left-hand side of (4) produces quadratic non-linearities, we are still able to solve the system explicitly for the two-and three-body cases (see section 3), and to provide some results for the structure of the singularities of the solutions of (4) for higher N (see section 4).

The trick
Let us rewrite the system of equations (5) in terms of a new independent real variable t (the physical time) defined by τ = [e iωt − 1]/iω (6) and rescaled dependent variables System (5) then becomes where the dot denotes a derivative with respect to t. It is now clear that if a solution ζ(τ) ≡ (ζ 1 (τ), . . . , ζ N (τ)) of (5) is known for the set of initial data ζ j (0), ζ j (0), then one can determine the solution to (8) corresponding to the initial data If the complex plane is identified with the real physical plane via the relation then the system of ODEs (8) can be interpreted as the equations of motion of N points in the plane subject to mutual two-body forces under the influence of an external field. The symmetric nature of the two-body interactions can be seen more clearly by rewriting system (8) in the form which includes the centre of mass N k=1 z k on the left-hand side. However, the physical interpretation is limited due to the fact that the equations (8) are not invariant under rotations. In fact, rotational invariance corresponds to the complex system (8) being invariant under the transformation z j → e iθ z j , which is clearly not the case. For a complete account of the rotationally invariant models in the plane obtained by complexification see chapter 4 in [13]. For a simpler model, which is invariant under both rotations and translations, a similar treatment to the one performed in this article can be found in [9]. In the case of the model in [9] it is apparent that the centre of mass performs a circular trajectory, which is clearly not the case in (11).
Despite the fact that the interaction (8) is rather complicated, involving non-linear and velocity-dependent forces, we shall shortly see that the behaviour of the system features nice and surprising properties. For instance, in the N = 3 case we will see that for generic initial conditions, all the trajectories are periodic. Here 'generic' initial conditions are arbitrary except for the exclusion of a set of measure zero where singularities can appear. For higher N, periodic motion occurs for a subset of initial conditions having the full dimension of the phase space.
Let us now explain the consequences that the change of variables (6) has for the behaviour of the solutions of (8). The transformation (6) introduces the constant ω to which we associate the fundamental period We now observe that as the real variable t (the physical time) varies from 0 to T , the (complex) variable τ goes from τ = 0 back to τ = 0 describing a circular contourC on the upper half plane with its centre at i/ω and radius 1/ω. Another way to see this is in terms of the inverse of (6), namely which defines t as a multi-valued function of τ, with the branch point of the logarithm at τ = i/ω, the centre of the circleC. Thus for each circuit made aroundC in the τ-plane, a branch cut is crossed and t increases by 2π/ω. We can now state the following important theorem relating the analyticity properties of the solutions of (5) to the periodicity of the solutions of (8): (5) is a holomorphic function of τ both inside and on the circular contourC, then the corresponding solution z(t) ≡ (z 1 (t), . . . , z N (t)) of the system (8) is non-singular and completely periodic in the real time t, with period T . Similarly, if a solution of the system (5) is meromorphic on this disk, with no poles on the boundaryC, then the same conclusion holds for the corresponding solution z(t). Moreover, if the only singularities of ζ(τ) inside the disk enclosed byC are a finite number of algebraic branch points, then the corresponding solution z(t) will again be completely periodic with period an integer multiple of T .
The proof of theorem 2.1 should be fairly obvious, in the light of the preceding discussion, but more details can be found in earlier works of Calogero et al [5]- [7], [9]. Theorem 2.1 enables one to find many periodic solutions to the dynamical system (8) provided it is known that the solutions ζ(τ) of (5), for a certain set of initial conditions, are holomorphic or meromorphic functions of τ in some region of the complex τ-plane enclosing the contourC. It is useful to understand this relation directly from the analyticity of the coefficients c m (τ). If the solutions c m (τ) of (4) are analytic or meromorphic (and hence single-valued) inside and onC, then (up to a factor of e −iωt ) the physical co-ordinates z j (t) are the zeros of a polynomial whose coefficients are periodic functions of t with period T , namelyc m (t) = c m (τ(t)). This entails that the complete set of zeros is also periodic with period T , but after one period some of the zeros might have interchanged their positions. The system will therefore be periodic with a period given by pT where is the least common multiple of the elements of some partition {p 1 , . . . , p k } N, corresponding to subsets of p j particles interchanging their positions after each fundamental period T . Since the number of possible partitions of any integer N is finite, the upper bound on the largest possible period pT is easily obtained by maximizing (14) over all partitions.
In fact, all solutions ζ(τ) of (5) corresponding to arbitrary non-singular initial conditions ζ(0) and ζ (0) (i.e. excluding a set of zero measure in phase space, namely the set of diagonals {ζ j (0) = ζ k (0) | j = k; j, k = 1, . . . , N}) are holomorphic in the neighbourhood of τ = 0. This is a consequence of the standard theorem [14] which guarantees the existence, uniqueness and analyticity of the solutions of analytic ODEs, in a sufficiently small disk D around the origin in the complex τ-plane (see theorem A.1). The size of the disk D is determined by the location of the singularity of ζ(τ) closest to the origin. If for some initial conditions ζ(0), ζ (0) (which determine z(0) andż(0) via (9)) the disk D encloses the circleC, then we know from theorem 2.1 that the trajectories z j (t) will be periodic with period pT . We shall come back to these considerations in section 5, while estimates on the initial conditions that ensure periodicity are given in the appendix.

Solving the system for the coefficients c m (τ)
In this section we solve system (4) in all known cases where an explicit solution can be obtained. In the rest of the cases where no explicit solution is available, we perform Painlevé analysis to determine the singularity structure of the solutions (see section 4).

The two-body case
In the two-body case, system (4) reads whose general solution depending on four arbitrary constants τ 0 , g 3 , α and β is where ℘ ≡ ℘(τ − τ 0 ; 0, g 3 ) denotes the Weierstrass ℘-function with invariants g 2 = 0 and g 3 , and ζ denotes the Weierstrass ζ-function with the same arguments as ℘; this ζ should not be confused with the zeros of the polynomial P satisfying the dynamical system (5) above. The free constant τ 0 corresponds to the fact that system (4) is autonomous. From now on, unless otherwise stated, we shall write ℘ as the shorthand for ℘(τ − τ 0 ; 0, g 3 ), and similarly for the arguments of every other elliptic function below. We note that the only singularities of c 2 (τ) are a double pole at τ 0 and congruent points, and c 1 (τ) is a linear combination of one solution that goes as α(τ − τ 0 ) −2 and another that behaves as β(τ − τ 0 ) 3 when τ is near τ 0 . We know therefore from the considerations in the previous section that system (8) for N = 2 is periodic for a generic set of initial conditions. Generic initial conditions exclude the cases when a pole of the functions c 1 , c 2 lies on the circleC; this can happen only if τ 0 is congruent to some point oñ C modulo the periodic lattice of the elliptic function. The only possible periods are T and 2T . Some simulations corresponding to this type of motion can be found in section 5.1.

The three-body case
When N = 3, the system of equations for the coefficients (4) becomes whose general solution depending on six arbitrary constants τ 0 , g 3 , α 1 , β 1 , α 2 and β 2 can be written as Once again, in this case the only singularities c 1 (τ), c 2 (τ) and c 3 (τ) have are double poles at τ = τ 0 . Since the general solution (18) is meromorphic, we know from theorem 2.1 that all solutions of system (8) for N = 3 are periodic for any set of initial conditions, excluding the subset of measure zero when a pole of the elliptic function lies on the circleC. The only possible periods in this case are T , 2T and 3T . Some simulations corresponding to this type of motion can be found in section 5.2.

A special solution: the linear case
As we see from (4) each odd (even) coefficient is coupled to the subsequent odd (even) coefficient, and all of them are coupled to c 2 . A special solution can be easily found for the case c 2 = 0, but in this case (4) implies that all even coefficients must be zero too. The system then becomes linear in the odd coefficients. If we denote by x the smallest integer greater than or equal to x, and we let M = N /2, then the system for the M odd coefficients becomes whose solutions are the following polynomials: where k = 1, . . . , M. This solution depends on the 2M arbitrary constants . . , B 2M−1 , and it is obviously entire. This special solution occurs for the set of initial conditions z j (0),ż j (0) determined by the conditions

Another special solution
Suppose that all even coefficients except c 2 are zero; the system (4) now reads In this case we can still express the solution to the system by quadratures. As we already know, the solution to the last equation in (21) is The next equation for c 2M−1 is of Lamé type, and the rest of the equations are also a Schrödinger equation of Lamé type with an inhomogeneous term given by the solution of the previous equation of the sequence. We can use the Wronskian method to write the following solution by quadratures: where the functions n (τ, α, β) are defined by and and by convention F 0 [ϕ(τ)] = ϕ(τ). The first few functions are It is not difficult to see that all the functions k (τ, α, β) are analytic, except for 0 which has a double pole at τ = 0. Therefore, in this special case, the coefficients c m (τ) are all meromorphic functions of τ. A family of special meromorphic solutions to (4), which includes this particular solution, is also observed in the Painlevé analysis performed in the following section where it corresponds to a non-principal balance.

The N-body case: Painlevé analysis
In this section we carry out Painlevé analysis on the set of coupled equations for the evenindexed quantities c 2k (τ) that appear as coefficients in the polynomial (1) with zeros ζ j (τ). As noted previously, the equations for the odd coefficients c 2k+1 constitute a sequence of (generally inhomogeneous) linear Schrödinger equations all with the same potential 2c 2 , which are solved recursively by taking a multiple of the solution to the previous equation as the inhomogeneous term at each step. The potential is given in terms of c 2 which must first be obtained from the solution of the coupled non-linear system for the even coefficients, namely the k coupled equations where N = 2K or N = 2K + 1 depending on whether there is an even or an odd number of particles. Since singular points in linear differential equations can only arise from singularities in the coefficient functions, the singularity structure of the odd coefficients is entirely determined by the function c 2 , and hence by the solution of the non-linear system (26) of order 2K. For N = 2 and 3, it is clear from the explicit solutions (16) and (18) that the coefficients c m are all meromorphic functions of τ. For N 4 we are unable to find the general solution explicitly, so we resort to applying the Painlevé test for ordinary differential equations [15,16] to see whether the general solution has movable branch points.

The case N = 4
In order to determine whether the solutions of system (26) can have branching, we must first look for the so-called dominant balances. For N = 4, the system for the even coefficients c 2 , c 4 reads simply To find the dominant balances, we look for leading-order singular behaviour of the form corresponding to a singularity in the solution at τ = 0 for at least one of µ, ν negative. The position of the singularity is movable, since we can always set τ → τ − τ 0 for an arbitrary point τ 0 in the complex τ plane. However, since the system is autonomous there is no loss of generality in carrying out all the analysis around τ = 0. There are three possible dominant balances for system (27), namely Another possible power-law behaviour around τ = 0 corresponds to µ, ν both being non-negative integers and leads to Taylor-series expansions, which are not relevant to our analysis of singular points.
The next step in applying the Painlevé test is to find the resonances, which correspond to the positions at which arbitrary constants appear in the Laurent expansions with leading terms given by (28). For system (27) to possess a strong Painlevé property we require that all resonances for all dominant balances must be integers, and at least one balance must have one resonance value of −1 with the rest being non-negative integers, in which case this is a principal balance for which the Laurent expansion should provide a local representation of the general solution. To find the resonance numbers r we substitute into the dominant terms of system (27) for each of the balances (i)-(iii), and take only the terms linear in δ and . This yields a pair of homogeneous linear equations for δ and (which correspond to the arbitrary coefficients appearing at the resonances). The determinant of this 2 × 2 system must vanish, which gives in each case a fourth-order polynomial in r.
Principal balance (i): It turns out that the balance (i) is the only principal balance, with resonances r = −1, 0, 5, 6.
The resonance −1 is always present, since it corresponds to the arbitrary position τ 0 of the pole, while r = 0 comes from the arbitrary constant b in the leading-order term of the expansion for c 4 ; the other two values arise from arbitrary coefficients higher up in the series for c 2 , c 4 , so that altogether there should be four arbitrary constants appearing in these Laurent series. However, for the Painlevé test to be satisfied we also require that all resonance conditions hold: so far we have only found the orders in the series where arbitary constants may appear, but it is necessary to check that all other terms vanish at this order when the series are substituted into the equations. Taking in each of the equations (27) we know already that the leading-order terms require At the next orders we further obtain k 2,2 = −3b 2 /5, k 4,2 = 7b 3 /15; k 2,3 = 0, k 4,3 arbitrary, so that the resonance condition at r = 5 corresponding to k 4,3 is satisfied. However, at the next order in the first equation of system (27), at the first appearance of the resonance coefficient k 2,4 , we find the additional relation which means that the resonance condition is not satisfied unless b = 0, contradicting the fact that b should be arbitrary. Thus the Painlevé test is failed by this principal balance. The only way to rectify the failure of the resonance condition and leave b as a free parameter is to modify (29) by adding logarithmic terms. More precisely, taking the resonance condition is resolved by taking However, the additional terms 2 , 4 in (30) must then consist of a doubly infinite series in powers of τ and log τ, with the leading order behaviours given by (31). Only in this way is it possible to represent the general solution of system (27) as an expansion in the neighbourhood of a singular point containing four arbitrary parameters. Such an infinite logarithmic branching is a strong indicator of non-integrability [15,16].
The presence of the negative integer value r = −5 means that this is a non-principal balance.
(For an extensive discussion of negative resonances see [17].) This gives Laurent expansions In this case all resonance conditions are satisfied and all higher coefficients in (32) are determined uniquely in terms of k and b. However, because it only contains three arbitrary constants (namely b, k and the position τ 0 of the pole), it cannot represent the general solution, but can correspond to a particular solution which is meromorphic.
Non-principal balance (iii): For the balance (iii) the resonances are given by r = −1 and the roots of the cubic equation r 3 − 15r 2 + 26r + 280 = 0, which turn out to be a real irrational number and a complex conjugate pair, approximately r = −3.2676, 9.1338 ± 1.5048i.
While non-integer rational resonances are allowed within the weak extension of the Painlevé test [18], whereby algebraic branching is permitted, irrational or complex resonances lead to infinite branching, and (as already evidenced by the principal balance (i)) the system (27) cannot possess the Painlevé property. This non-principal balance may be interpreted as a particular solution corresponding to a degenerate limit of the general solution, and perturbation of this particular solution (within the framework of the Conte-Fordy-Pickering perturbative Painlevé test [17]) will pick up the logarithmic branching present in the general solution.

Balances for N > 4
For any N > 4 it is possible to repeat the above singularity analysis for the subsystem (26) for the even coefficients. This is a system of order 2K, for either an even number of particles, N = 2K, or an odd number, N = 2K + 1. In the previous subsection we have considered the case K = 2 in detail, and found one principal and two non-principal balances. However, we expect that the number of non-principal balances increases with K. Rather than attempting to enumerate all the possible balances for any K, we will concentrate on describing the analysis for the principal balance (the analogue of (i) above), which exists for any K 2 but fails the Painlevé test. We also find that for any K there is at least one non-principal balance which passes the test and should correspond to a special solution which is meromorphic.
Clearly, there are 2K − 1 non-negative integer resonances as well as the resonance −1 corresponding to the arbitrariness of the pole position τ 0 , so an expansion with leading order behaviour (33) should contain 2K arbitrary constants, corresponding to a local representation of the general solution.
Seeking local Laurent series expansions, denoted L 2k (τ), for the coefficient functions c 2k with leading orders given by (33), we find with the coefficients satisfying and 2k = 1 6 (2k (with the convention that any coefficient with index larger than 2K is zero). The relations (35) come from the first few terms obtained by substituting the expansions (34) into the first equation of the system, that is (26) with k = 1, while the relations (36) come from (26) for k 2. Along with τ 0 , the constants κ 2k , ξ 2k and f should constitute a further 2K − 1 arbitrary parameters, corresponding to the 2K − 1 non-negative resonances. However, although the resonance conditions for κ 2k and ξ 2k are identically satisfied, the last relation in (35) is the resonance condition for f , and this constitutes a constraint on the parameters. To see this, note that putting k = 2 into (36) yields while substituting for 4 in the final relation (35) implies Comparing (37) and (38) it is clear that the two different expressions for η 4 put a constraint between κ 4 , κ 6 and κ 8 , so that one of these coefficients is no longer arbitrary. Thus if the constraint is satisfied then one obtains a set of Laurent series (34) representing a particular solution depending on only 2K − 1 arbitrary parameters. As in the case K = 2 above, the only way to remove the constraint and get a local representation of the general solution is by adding infinitely many logarithmic terms to the Laurent series (34), so that for suitable coefficients a (j 1 ,j 2 ) 2k which can be determined recursively.
A non-principal balance: For any K 2 there is always a non-principal balance containing K + 1 arbitrary constants which passes the Painlevé test; this is the analogue of the balance (ii) described above. The leading-order behaviour is with resonances r = −5 (K − 1 times), −1, 0 (K − 1 times), 6. The resonance conditions for the leading-order coefficients ξ 2k as well as the resonance condition at r = 6 in the Laurent series for c 2 are all trivially satisfied. In fact it is simple to obtain this balance from the Laurent series (34) in the principal balance by setting Thus we expect that this non-principal balance corresponds to a particular solution which is meromorphic. The problem of finding all other non-principal balances in (26) for any K is more difficult, but for our purposes it is sufficient to realize that the principal balance contains logarithmic branching. We expect that this should be the origin of aperiodic solutions in the planar dynamics (8) of the functions z j (t), which are given in terms of the zeros ζ j (τ) of the polynomial by (7).

Analysis of various motions
In this section we analyse different orbits of the system (8) corresponding to different initial conditions. To test the predictions of the previous sections on the periodicity of the orbits, we have performed a numerical integration of the equations of motion. We have used computer software developed by one of the authors [19], which uses an embedded Runge-Kutta method of integration with variable time-step.
In order to explore the whole space of initial conditions in a systematic way, let us note the following proposition.
is a solution of the system which has the same form as (8) with ω replaced bỹ The transformation (39) corresponds to the following relation between the initial data: Let us fix, without loss of generality, the time-scale so that the fundamental period is unity by setting  As mentioned in section 2.1, the initial data z j (0),ż j (0) determine the singularity structure of the solution ζ(τ) of (5). From the previous proposition we also know that different initial data obtained from z j (0),ż j (0) by the one-parameter group of transformations (42) correspond via (6) and (7) to the same solution ζ(τ) of (5), by suitably choosing the parameter ω in (6) and (7) as ω = 2π/λ. In other words, we can obtain the solutions of (8) for different sets of initial data related by (42) by travelling on the Riemann surface associated to the same solution ζ(τ) of (5) on circular contoursC =C(λ) of radius λ/2π. The effect of scaling the initial data by (42) is to enlarge (or decrease) the radius of the disk enclosed byC by a factor of λ. We will consider an increasing sequence of values of λ in our numerical simulations, and we will see how the orbits transform. For sufficiently low values of λ, the disk becomes small and eventually ζ(τ) will be analytic inside and onC (and therefore we will observe periodic orbits with period T = 1, as stated in theorem 2.1). Then as λ increases, some of the singularities of ζ(τ) might enter the circleC(λ). The behaviour of the trajectories z j (t) will change drastically when this happens, and the periodic or aperiodic character will depend on the nature of the singularity (pole, algebraic branch point, logarithmic branch point, etc) that lies inside the circleC. In this manner, we can systematically explore the whole space of initial conditions.

Two-body problem
We already know from section 3.1 that for N = 2 the problem is solvable and all the orbits of the system corresponding to arbitrary initial conditions must be periodic, with period either T = 1 or T = 2. The numerical integration of the equations of motion confirms this result.
We will consider the sequence of motions with initial data given by for increasing values of λ. The corresponding motions have been collected in table 1. For λ 1.09 the motion is periodic with period 1 (this corresponds to the region around τ = 0 in which the solution ζ(τ) is ensured to be analytic). When λ = 1.09 the motion is still periodic with period 1, but the two particles almost collide at a certain instant of time (see figure 1(b)). Initial conditions for which there is a collision of particles correspond to a singularity of ζ(τ) lying exactly on the circular contourC, and in those cases the solution is only welldefined up to the finite collision time. When λ = 1.10, the radius ofC has increased and now it has encircled an algebraic branch point of order 2 of ζ(τ). As a consequence, the motion of the system is periodic with period T = 2 ( figure 1(c)). For higher values of λ it could happen that the circleC encloses two second-order branch points joined by a cut, so that the associated Riemann surface has only one sheet, and the period is again T = 1. This happens for instance when λ = 20, as can be seen in figure 1(d).

Three-body problem
The behaviour of the three-body system is not essentially different from the two-body case. We know from section 3.2 that this problem is solvable, all coefficients c m (τ), m = 1, 2, 3, are meromorphic functions of τ, and therefore the functions ζ j (τ) have only a finite number of algebraic branch points. The system is therefore periodic for arbitrary initial conditions, and the only possible periods are T = 1, 2 and 3.
We have explored the motion for initial conditions of the form for different values of the parameter λ, obtaining the results summarized in table 2. For values of λ near 1.65 and 2.32, the period passes from 1 to 2 and from 2 to 3 respectively, as the circleC in the complex τ-plane includes more rational branch points of ζ(τ). As for the two-body case, for higher values of λ the period decreases again to T = 2 and 1.

Four-body problem
The motions for N = 4 are much richer than in the previous cases due to the fact that the functions ζ(τ) in (5) in this case will have other singularities besides poles and algebraic branch points (figure 2). We have not been able to solve the system of coupled non-linear equations (28) for the coefficients c m (τ), but the Painlevé analysis in section 4 reveals that the set of equations does not pass the Painlevé test, and the general solution has logarithmic branching in the complex τ-plane.
Let us explore the motions that originate from the set of initial conditions given by for different values of the parameter λ. All orbits are completely periodic. For sufficiently low values of λ the system is periodic with period 1 as expected. Then for 1.3 λ 2.8 we observe periodic orbits with higher periods (the only possible periods are T = 1, 2, 3 and 4, see equation (14). But when λ 3.8 the trajectories seem to be no longer periodic (see table 3). This can hardly be appreciated in figure 3(g) where the whole trajectory from t = 0 to 4 has been plotted. It is much better seen in figure 3(h) which corresponds to a zoom on a region closer to the origin of the same motion for the same time interval. In this last figure the presence of attractor trajectories seem to appear. Of course, a numerical calculation does not constitute any proof that the trajectories are aperiodic, and the orbits have thus been named by the acronym HSL (Hic Sunt Leones) following the notation introduced in [8,9]. This emphasizes that it remains an open question to understand the nature of these types of orbits (periodic, quasi-periodic, chaotic, ergodic, etc).
However, Painlevé analysis of the non-linear ODEs for the coefficients c m (τ) suggests that the solutions ζ(τ) of (5) will have logarithmic branch points, which for real t should lead to aperiodic behaviour in the corresponding solution z(t) of (8). The transition from periodic orbits to HSL that takes place near λ 3.8 is probably due to the fact that the circleC for that value of λ encloses a logarithmic branch point of ζ(τ). This infinite branching is responsible for the aperiodic character of the orbits.
In the previous work [9], the aperiodic character of the orbits was thought to be caused by the fact that new rational branch points enter the circle in higher sheets of the Riemann surface, thus opening additional sheets in which other rational branching might occur, thus forming an infinite sequence of sheets. The situation in this problem can be different, since we know from  the Painlevé analysis of the equations that logarithmic branching should be present. Just one logarithmic branch point inside the circleC will cause the orbit to be non-periodic. The mechanism described in this section, by which higher-period orbits are originated as branch points of the solutions ζ(τ) of (5) enter the circleC is reminiscent of the Ruelle-Takens scenario of the onset of chaos by eigenvalue crossing [1]. However, in this scenario, the crossing of just one branch point can cause a transition from an orbit with period one to an orbit with period p > 2, or even to a chaotic orbit.

Conclusions and future work
We have constructed an N-body problem in the plane that models N particles subject to mutual velocity-dependent forces under the presence of an external anharmonic field and a Lorentz-like force (an external force similar to the Lorentz force caused by a magnetic field acting on the axis perpendicular to the plane if the particles had unitary (positive) charge). It has the following properties: • For N = 2 and 3 the trajectories are completely periodic, for arbitrary initial conditions apart from a set of measure zero (corresponding to singular trajectories). Moreover, the problem is actually solvable, and the positions can be calculated at any instant of time by computing the zeros of a polynomial whose coefficients are known periodic functions of the real time t. These theoretical results are confirmed by a numerical integration of the equations of motion. • For N 4 the problem is not solvable for arbitrary initial conditions. However, the trajectories that originate from a set of initial conditions that has non-vanishing measure in phase space are assured to be periodic. Painlevé analysis of the system of ODEs for the coefficients of the polynomial suggests that, in general, its zeros will have logarithmic branch points, and the motions will not be periodic. Both periodic and (seemingly) non-periodic orbits have been observed numerically.
A complete characterization of the chaotic nature of these orbits requires a better knowledge of the singularities of the solutions of system (5): these are more complicated than the singularities of the coefficients c m (τ) satisfying (4), as found in section 4. This more detailed analysis is currently being studied. Another extension of this work is to analyse the possible motions when different coupling constants are inserted for each pair force, as was done in [9,20]. In that case, the equations of motion are no longer related to the evolution of a polynomial, and Painlevé analysis should be performed directly on the equations analogous to (5). Numerical experiments have been conducted in the two-body case and suggest that there exist initial conditions for which the orbits are periodic (particularly, some islands of periodicity are observed for special values of the coupling constants). An investigation of the properties of these systems from the numerical analysis point of view seems interesting, since it is quite remarkable that the numerical integration of these complicated non-linear equations shows such good stability.

Appendix. Estimates for holomorphic solutions
Here we provide bounds on the initial conditions which ensure the existence of periodic solutions to the system (8). This is achieved in a straightforward way by finding sufficient conditions for the existence of holomorphic solutions of the associated system (5) inside a sufficiently large circle around the origin in the complex τ-plane. Provided that this circle is large enough to enclosẽ C, theorem 2.1 then guarantees periodicity of the corresponding solutions of (8), with the two sets of initial conditions for the different systems being related by (9). Thus for our purposes it is sufficient to derive bounds on the initial conditions for (5) that ensure holomorphicity in a suitably large region.
To apply standard results on the existence of holomorphic solutions to initial value problems, such as those in [14], it is necessary to rewrite the system (5) in the following first-order form: w j = c −1 w N+j + ζ j (0), j = 1, . . . , N, In the above, c is a positive constant which for the moment is arbitrary, and as the index j runs from 1 to N the 2N new independent variables in (A.1) are defined by w j (τ) = ζ j (τ) − ζ j (0), w N+j (τ) = c(ζ j (τ) − ζ j (0)), so that w j (0) = 0 for all j. Note also that compared with (5), we have absorbed the term 2(N − 1)ζ 3 j on the left-hand side into the sum appearing on the right of (A.1). We now introduce the bounds on the initial data for (5) and on the variables w j which ensure that the right-hand sides of (A.1) are bounded holomorphic functions of w j : 0 < Z 1 := min{|ζ j (0) − ζ k (0)| : j = k; j, k = 1, . . . , N}, Proof. This is a straightforward application of the result of section 12.21 in [14], setting m = 2N + 1 and a = ∞: because the system (5) is autonomous, a can be made arbitrarily large.
Remark. To ensure that the solutions of the transformed system (8) are periodic in real time t, it is necessary to take initial data such that ρ > 2/ω, so that the circleC lies inside the circle of radius ρ in the τ-plane, and the corresponding solutions of (5) are guaranteed to be holomorphic inside and onC. The bounds Z 2 , Z 3 on the initial conditions, as well as the two parameters r 1 , r 2 with the necessary constraints (A.8) and (A.9), can be varied suitably to make ρ > 2/ω.