Solvable Many-Body Models of Goldfish Type with One-, Two- and Three-Body Forces

The class of solvable many-body problems"of goldfish type"is extended by including (the additional presence of) three-body forces. The solvable $N$-body problems thereby identified are characterized by Newtonian equations of motion featuring 19 arbitrary"coupling constants". Restrictions on these constants are identified which cause these systems - or appropriate variants of them - to be isochronous or asymptotically isochronous, i.e. all their solutions to be periodic with a fixed period (independent of the initial data) or to have this property up to contributions vanishing exponentially as $t\rightarrow\infty$.


Introduction
Over three decades ago [3] a class of solvable N -body problems featuring several free parameters ("coupling constants") was introduced by identifying the coordinates z n (t) of the moving particles with the N zeros of the time-dependent (monic) polynomial ψ(z, t) (of degree N in the independent variable z), itself evolving according to a linear partial differential equation (PDE) -suitably restricted to guarantee that it feature a polynomial solution of degree N in z. The simplest dynamical system belonging to this class displays remarkably neat properties and was therefore considered a "goldfish" (for a justification of this terminology, see [5]; subsequently this terminology has been often employed to identify N -body models belonging to this class, and this justifies its use also in the present paper, including its title). The behavior of the solutions of these N -body problems has been variously investigated and also used to arrive at related mathematical results, such as finite-dimensional representations of differential operators and Diophantine properties of the zeros of certain polynomials: see the two monographs [4] and [6] and the references quoted there (including the more recent ones added to the 2012 paperback version of [6]), and the more recent papers [1,2,7,8,9,10,11,12,13,14,15,16,17]. The class of linear PDEs satisfied by the polynomial ψ(z, t) which subtended these investigations was so far restricted to PDEs featuring derivatives up to second order. The restriction to time-differentiations of second order was motivated by the interest in N -body problems featuring equations of motion of Newtonian type (i.e., "acceleration equals force"); the more general restriction to differentiations of at most second-order implied that the N -body problems under consideration only involved one-body and two-body (generally velocity-dependent) forces. In the present paper we extend our consideration to third-order z-derivatives -while maintaining the restriction to second-order time-differentiations so as to only treat N -body problems of Newtonian type. This entails that the corresponding solvable N -body problems thereby identified involve, additionally, three-body forces.
In the following Section 2 we report our main results, namely we display the Newtonian equations of motion of the N -body problems identified in this paper and we indicate how they are solved by algebraic operations. In a terse Section 3 these findings are proven: these developments rely on identities reported and proven in Appendix A, constituting a substantial part of this paper. In Section 4 we highlight some special cases in which these N -body problemsor appropriate variants of them -are isochronous or asymptotically isochronous, i.e. all their solutions are periodic with a fixed period (independent of the initial data) or they feature this property up to contributions vanishing exponentially as t → ∞. A terse Section 5 entitled "Outlook" outlines possible future developments.

Main results
In this section we report the main results of this paper, which are then proven in the following section. But firstly let us specify our notation.
Notation 2.1. The N coordinates ("dependent variables") of the point, unit-mass, moving particles which are the protagonists of the N -body problem treated in this paper are generally denoted as z n ≡ z n (t), with t ("time") the "independent variable". As usual differentiation with respect to time is denoted by a superimposed dot, henceż n ≡ dz n /dt,z n ≡ d 2 z n /dt 2 . As we just did here, often the indication of the t-dependence is not explicitly displayed. We generally assume that these coordinates z n are complex numbers, so that the points with coordinates z n (t) move in the complex z-plane; but special subcases in which the coordinates z n are real are of course possible; and it is of course also possible to reinterpret motions taking place in the complex zplane as instead taking place in the real horizontal plane by identifying the real and imaginary parts of the complex numbers z n = x n + iy n as the Cartesian components of the real 2-vectors r n = (x n , y n ) (as explained in Chapter 4 of [4], entitled "Solvable and/or integrable many-body problems in the plane, obtained by complexification"). Here and hereafter i is the imaginary unit, i 2 = −1, N is an arbitrary positive integer (generally N ≥ 2), and it is understood that subscripts such as n (and also m, k, , but not j; see below) run over the positive integers from 1 to N (unless otherwise specified); the reader is often (but not always) explicitly reminded of this fact. We also use occasionally the notation z to denote the (generally complex ) N -vector of components z n , z ≡ (z 1 , . . . , z N ); and likewise for other underlined letters (see below). A key role in our treatment is played by the time-dependent (monic) polynomial ψ(z, t) of degree N in the (scalar, generally complex ) variable z which features the N coordinates z n (t) as its N zeros, see (1.1).
The solvable N -body problem treated in this paper is characterized by the following Newtonian equations of motion: Here and hereafter the 19 upper-case letters A j , B j , D j , E, F j , G j denote time-independent parameters ("coupling constants"). They are a priori arbitrary (possibly complex ) numbers; special assignments are occasionally considered below. This notation is chosen consistently with equation (2.3.3-2) of [4] (with C = 1), to which this system of N equations of motion reduces when the 9 "new" parameters vanish, i.e.
Note that the terms associated with the 4 new parameters F 0 , F 1 , F 2 , F 3 , represent velocityindependent three-body forces, the terms associated with the 3 new parameters G 0 , G 1 , G 2 represent velocity-dependent three-body forces, and the terms associated with the new parameters F 4 respectively G 3 represent velocity-independent respectively velocity-dependent one-body and three-body forces. (Of course terms representing one-body forces can always be absorbed in those representing many-body forces: for instance the one-body term −2(N − 1)A 3 z 2 n in the right-hand side of (2.1) could be eliminated by replacing z 3 n with z 2 n z in the two-body term multiplying A 3 in the first sum in the right-hand side of (2.1)).
This Newtonian N -body system is solvable by algebraic operations because -as shown in the following Section 3 -the coordinates z n (t) evolving according to this system of N Ordinary Differential Equations (ODEs) coincide with the N zeros of the time-dependent polynomial (1.1) of degree N in z, itself evolving according to the following linear PDE: Here the 19 upper-case letters are of course the same time-independent parameters featured by the Newtonian equations of motion (2.1). Again, this notation is consistent with that used in [4]: indeed this PDE is a natural generalization of equation (2.3.3-1) of [4] (with C = 1) to which it clearly reduces when the new parameters vanish, see (2.2). This PDE implies that the N coefficients c m ≡ c m (t), see (1.1), evolve according to the following system of N linear ODEs: of course with c 0 = 1 and c m = 0 for m < 0 and for m > N .
Remark 2.1. This system of N ODEs satisfied by the N coefficients c m (t) justifies the assertion made above that the linear PDE (2.3) admits a polynomial solution: specifically, the polynomial (1.1) of degree N in z. Of course this system of N ODEs reduces to equation (2.3.3-8) of [4] (with C = 1) when the new parameters vanish, see (2.2). Likewise, it reduces to equation (4.54) of [6] (up to the correction of a trivial misprint in that equation, and to an obvious notational change). Note however that the analogous equation has been wrongly reported (as equation (3)) in [13]: due to a trivial misprint (a multiplicative factor c m−2 missing in the first term in the second line) and the mistake of inserting two terms in the right-hand side (which should instead just be zero); fortunately this mistake has no consequence on the remaining part of that paper -except for the mistaken Remark 1.1 which should of course be ignored.
This system, (2.4), of N autonomous linear ODEs is of course solvable by algebraic operations. Indeed its general solution reads m are a priori arbitrary -to be fixed a posteriori in order to satisfy the 2N initial condi- m are the 2N eigenvectors, respectively the 2N eigenvalues, of the following (time-independent) generalized matrix-vector eigenvalue problem: Here of course I is the N × N unit matrix (I mn = δ m,n where, here and below, δ m,n is the Kronecker symbol) and the two N × N matrices U and V are defined componentwise as follows: This implies of course that the 2N eigenvalues λ (±) m are the 2N roots of the following polynomial equation (of degree 2N in λ): These findings show that the solution of the system (2.4) is achieved by the algebraic operation of determining the eigenvalues and eigenvectors of the matrix-vector generalized eigenvalue problem (2.6a).
The algebraic equation (2.6d) can be explicitly solved (for arbitrary N ) in the two special cases in which the two N × N matrices U and V are either both upper triangular or both lower triangular.
The first of these two special cases obtains if, of the 19 parameters in (2.4), the following 4 vanish: The second of these two special cases obtains if instead, of the 19 parameters in (2.4), the following 9 vanish: It is then easily seen that -in both these two cases -the 2N eigenvalues λ These findings imply that -at least in these cases -specific predictions on the actual behavior of the solutions of the Newtonian N -body problem (2.1) can be easily made. In Section 4 we identify in particular the cases in which this Newtonian N -body problem is isochronous or asymptotically isochronous, or an appropriate variant of it is isochronous; while the findings described in this Section 2 are proven in the following Section 3.

Derivation of the equations of motion of the Newtonian N -body problem
Our task in this section is to obtain the equations of motion (2.1) characterizing the new N -body problem of goldfish type. The starting point is the linear third-order PDE (2.3) satisfied by the polynomial (1.1). We already saw in Section 2 the implication of this evolution PDE for the coefficients c m (t) of its polynomial solution (1.1), leading to the identification -via (1.1) and (2.5) -of its solution by algebraic operations. In this section we show that the fact that the polynomial (1.1) satisfies the PDE (2.3) implies that its zeros z n (t) indeed evolve according to the Newtonian equations of motion (2.1). This is in fact an immediate consequence -via trivial, if somewhat cumbersome, algebraof (some of) the identities reported in Appendix A of [6] and of the additional identities (A.2) and (A.3) reported (and proven) in Appendix A of this paper, see below. Indeed these identitiesvalid for an arbitrary time-dependent polynomial ψ(z, t) of degree N in z, see (1.1)) -allow to transform the linear PDE (2.3) satisfied by the polynomial ψ(z, t) into the system of nonlinear ODEs (2.1) satisfied by its zeros z n (t). Note that the fact that this outcome obtains is both a consequence and a confirmation of the fact that the PDE (2.3) admits as its solution the polynomial ψ(z, t), see (1.1), of degree N in z, hence featuring N zeros z n (t).

Isochronous and asymptotically isochronous cases
In this section we identify cases in which the Newtonian N -body problem (2.1) is isochronous or asymptotically isochronous, or variants of it are isochronous.
First of all we note, see (2.5), that if the 2N eigenvalues λ with the period independent of the initial data. Here of course q is the minimum common multiple of the 2N denominators of the 2N rational numbers r m . And it is plain that the same property of isochrony is then shared by the N coordinates z n (t), namely the Newtonian N -body problem (2.1) is then isochronous as well; with the possibility that in some open regions of its phase space the periodicity only holds for a period which is a (generally small ) integer multiple of T due to the fact that some of the zeros of the polynomial ψ(z, t) -itself evolving isochronously with period T , see (1.1) -might "exchange their roles" over the time evolution (for an analysis of this phenomenology, also explaining the meaning of the assertion made above that the integer multiple in question is generally small, see [18]).
Let us then focus on the two cases -as identified at the end of Section 2, see (2.7) and (2.8)in which the 2N eigenvalues λ (±) m can be explicitly obtained, see (2.9). It is then easy to identify the additional restrictions on the parameters which are necessary -and also sufficient, up to some minor additional restrictions to exclude the coincidence of eigenvalues -to guarantee that the Newtonian N -body problem (2.1) be isochronous with period T . They read: 3e) Here r 1 , r 2 , r 3 , r 4 are 4 arbitrary rational numbers and s 1 and s 2 are two arbitrary signs (+ or −). Indeed with these assignments λ (±) m clearly satisfies the condition (4.1) with It is moreover plain that, if the 2N eigenvalues λ (±) m satisfy, instead of the condition (4.1), the less restrictive condition -with ω again a positive real number, the 2N numbers r   Here T is given again by (4.2b), but now with q being the minimum common multiple of the denominators q (±) m 's of the rationals r (±) m 's associated with vanishing ρ (±) m 's, ρ (±) m = 0, see (4.5) and (4.1). And it is again plain that the same property is then shared by the N coordinates z n (t), namely that the Newtonian N -body problem (2.1) is then asymptotically isochronous as well (again, with a period which might be a, generally small, integer multiple of T [18]).
Next, less us investigate the cases in which -via a well-known trick, see for instance Section 2.1 (entitled "The trick") of [6] -isochronous variants can be manufactured of the Newtonian Nbody problem (2.1). To this end, it is convenient to re-write the equations of motion (2.1) via the formal replacement of dependent and independent variables z n (t) ⇒ ζ n (τ ), so that the equations of motion read as follows: Here and hereafter (in this section) appended primes denote differentiation with respect to the variable τ . We now perform the following change of dependent and independent variables ("the trick"): Here and below ω is a positive constant, and α a nonvanishing number that we reserve to assign, see below.
It is then plain that there hold the following formulas: And via these formulas one can easily obtain the equations of motion implied for the dependent variables z n (t) by the equations of motion (4.7) satisfied by the variables ζ n (τ ), and thereby ascertain for which assignments of the parameter α, and for which corresponding restrictions on the 19 coupling constants featured by these equations of motion, the resulting equations of motion satisfied by the N coordinates z n (t) are autonomous, i.e. they feature no explicit timedependence. We list below all these cases, on the understanding that all the coupling constants which do not appear in these new equations of motion have been set to zero (while those that do appear are arbitrary). For α = −2, these Newtonian equations of motion read For α = −1, these Newtonian equations of motion read For α = −2/3, these Newtonian equations of motion read . (4.11) For α = −1/2, these Newtonian equations of motion read For α = 1, these Newtonian equations of motion read For α = 2, these Newtonian equations of motion read n (4.14) And it is plain that all these Newtonian models are isochronous: see (4.8a), and note that the solutions ζ n (τ ) of the solvable N -body model (4.7) have at most algebraic singularities as functions of τ .

Outlook
The search for more general solvable models of goldfish type should continue: they might still yield interesting results.
Another development likely to yield interesting findings is the investigation of the behavior, in the infinitesimal neighborhood of its equilibria, of the N -body model introduced above, especially in the isochronous cases. This investigation might yield new Diophantine findings for the zeros of interesting polynomials. We plan to pursue these results, which shall eventually be submitted to a journal devoted to special functions if they turn out to be sufficiently interesting to justify their publication.

A Appendix: identities involving the zeros of a polynomial
In this Appendix we report (and then prove) several identities for the time-dependent polynomial ψ ≡ ψ(z, t), see (1.1), of degree N in z. We of course use hereafter the notation introduced in Section 2, see Notation 2.1, and in addition the following convenient shorthand notation: as in [6] (see there equations (A.4) and (A.5)) -with D a differential operator acting on the independent variables z and t of the polynomial ψ(z, t), see (1.1) -stands for the identity Below we often, for notational simplicity, omit to indicate explicitly the dependence on their arguments of ψ ≡ ψ(z, t), z n ≡ z n (t) and c m (t) (see (1.1)). We now list the following identities, which complement those reported in Appendix A of [6] (see in particular the 2012 paperback version, where the formulas denoted in [6] as (A.8k) and (A.8l) are corrected): Of course in the last two, (A.3a) and (A.3b), superimposed dots indicate t-differentiations. Also note that in the third, (A.2c), of these identities, consistently with (1.1) But remarkably -due to a neat cancellation, see (A.2c) and equation (A.6b) of [6] -this quantity does not appear in the equations of motion (2.1). Our task in this Appendix is to prove these identities. To perform these proofs it is convenient to introduce the following shorthand notation denoting sums over (dummy) indices restricted to take different values (among themselves): (A.5) We moreover introduce the following shorthand notation for the "denominator" den(z n ,z k ,z ), which clearly has the property to be invariant under the cyclic exchange n → k → of the three indices n, k, , den(z n , z k , z ) = den(z k , z , z n ), (A.6b) and to be instead antisymmetric under the exchange of any two of its three arguments z n , z k , z , den(z n , z k , z ) = − den(z k , z n , z ) = − den(z , z k , z n ) = − den(z n , z , z k ). (A.6c) Likewise, we denote corresponding "numerators" as num(n, k, ) (i.e., num(z n , z k , z ) ≡ num(n, k, )) and, whenever one of them is the sum of an arbitrary number of terms num j (n, k, ), i.e. num(n, k, ) ≡ The validity of this assertion is an obvious consequence of the antisymmetry -see (A.6c) and (A.8) -under the exchange of two appropriately chosen dummy indices in the triple sums nk [num j (n, k, )/ den(z n , z k , z )], which therefore vanish (for all values of j). Next, let us report and prove the following identities, valid for any set of N arbitrary numbers z n (for convenience, we always assume them to be all different among themselves).
In the last three formulas of course c 1 is defined by (A.4). The proof of the (well-known) identity (A.10a) is trivial: The first step is justified by adding to the left-hand side of (A.10a) what is obtained by the exchange of the dummy indices k and (which does not change the result) and dividing by 2 (this operation is often repeated below without describing it in as much detail as done here); the second step is immediately implied by the first definition (A.5).
Then the proof of (A.10b) is no less trivial: it follows from (A.10a) and the second definition (A.5).
The proof of (A.10c) goes as follows. For p = 0 by adding to the right-hand side of this formula the sums obtained by performing sequentially two cyclic transformations of the three indices n, k, and by taking advantage of the invariance property (A.6a)) -and by dividing the sum of the three equal sums thereby obtained by 3 -we clearly get and it is then plain that this quantity vanishes. For p = 1 via (A.6a) nk z n (z n − z k )(z n − z ) = nk z n z k − z n z den(z n , z k , z ) (A. 13) and it is plain that this quantity vanishes thanks to Lemma A.1, since in the right-hand side the first term in the numerator is invariant under the exchange of dummy indices n ↔ k and in the second under the exchange n ↔ . The proof of (A.10d) goes through the following steps: Now we can finally proceed and prove the formulas (A.2) and (A.3). We start from reporting equations (A.1) (which coincides with (1.1)), (A.2), (A.3), (A.8a) and (A.9a) of [6]: The last two equations correspond of course to (A8.a) and (A.9a) via the convention defining the symbol ⇐⇒, see (A.1). Partial differentiation of (A.22) with respect to z gives Via (A.20) this becomes (after a convenient cancellation and change of dummy index from n to k) We then use the identity getting thereby We then exchange the dummy indices n and k in the second of the two sums over these indices, getting thereby which can also be written as follows: This implies (using the definition of the symbols ⇔ and nk , see above) and finally, taking advantage of the fact that the part of the summand antisymmetric under the exchange of the two dummy indices k and can be eliminated -so that (z n + z k − 2z )/(z k − z ) ≡ 3/2 + (2z n − z k − z )/[2(z k − z )] can be replaced by 3/2 -one gets (A.2a) with p = 0, which is thereby proven. This is the first of the new formulas analogous to those reported in Appendix A of [6].
To proceed it is convenient to write in longhand the formula we just proved: Multiplication by z p then yields By replacing the z p in the numerator with z p n + (z p − z p n ) we get z p ψ zzz = ψ (via the identity z 3 − z 3 n = (z − z n )(z 2 + zz n + z 2 n )). And it is then plain that this formula, via (A.10c) and (A.10d), yields (A.2b), which is thereby proven.
Likewise, to prove (A.2c) we set p = 4 in (A.30), which then reads z 4 ψ zzz = ψ Next, to prove (A.3a) (to begin with, with p = 0), we z-differentiate (A.23) getting thereby It is now plain that (A.3b) is proven if we show thatσ(3; z,ż) vanishes,σ(3; z,ż) = 0. To prove this we make the following steps: σ(3; z,ż) = nk [den(z n , z k , z )] −1 z 2 n [(ż n +ż )(z k − z ) + (ż k +ż )(z n − z )] −ż n (z n − z k )(z n − z )(z k − z ) = nk [den(z n , z k , z )] −1 ż n (z k − z ) z 2 n − (z n − z k )(z n − z ) +ż k z 2 n (z n − z ) +ż z 2 n (z n + z k − 2z ) = nk ż n den(z n , z k , z ) Here the first equality is justified by the definition (A.6a) and a bit of trivial algebra, the second equality obtains by trivial algebra, the third equality obtains by the exchange of the dummy indices n ↔ k in the term multiplyingż k and likewise the exchange n ↔ in the term multiplyingż (under these exchanges the denominator changes sign, see (A.6c)), the fourth equality obtains by trivial algebra, and the final equality to zero is yielded by Lemma A.1 since the numerator is invariant under the exchange k ↔ .