A New Class of Solvable Many-Body Problems

A new class of solvable $N$-body problems is identified. They describe $N$ unit-mass point particles whose time-evolution, generally taking place in the complex plane, is characterized by Newtonian equations of motion"of goldfish type"(acceleration equal force, with specific velocity-dependent one-body and two-body forces) featuring several arbitrary coupling constants. The corresponding initial-value problems are solved by finding the eigenvalues of a time-dependent $N\times N$ matrix $U(t)$ explicitly defined in terms of the initial positions and velocities of the $N$ particles. Some of these models are asymptotically isochronous, i.e. in the remote future they become completely periodic with a period $T$ independent of the initial data (up to exponentially vanishing corrections). Alternative formulations of these models, obtained by changing the dependent variables from the $N$ zeros of a monic polynomial of degree $N$ to its $N$ coefficients, are also exhibited.


Introduction
In this paper a new class of solvable N -body problems is identified. They describe an arbitrary number N of unit-mass point particles whose time-evolution, generally taking place in the complex plane, is characterized by Newtonian equations of motion "of goldfish type" (acceleration equal force, with specific velocity-dependent one-body and two-body forces, see below) featuring several arbitrary coupling constants. The solvable character of these models is demonstrated by the possibility to ascertain their time evolution by purely algebraic operations. In particular it is shown below that their initial-value problems are solved by finding the eigenvalues of a timedependent N × N matrix U (t) explicitly defined in terms of the initial positions and velocities of the N particles. Some of these models are asymptotically isochronous, i.e. in the remote future they become completely periodic with a period T independent of the initial data (up to exponentially vanishing corrections). Alternative formulations of these models, obtained by changing the dependent variables from the N zeros of a monic polynomial of degree N to its N coefficients, are also exhibited.
The main idea to obtain these results is to identify the N coordinates z n (t) characterizing the positions of the particles of the N -body problem with the N eigenvalues of an N × N matrix U (t), itself evolving according to a solvable (or integrable) matrix ODE. This technique was invented long ago by Olshanetsky and Perelomov [1][2][3][4] and has been much exploited subsequently, identifying thereby many solvable (or integrable) many-body problems: see for instance [5] (in particular its Section 2.1.3.2 entitled "The technique of solution of Olshanetsky and Perelomov (OP)") and [6] (in particular its Section 4.2.2 entitled "Goldfishing"), and the more recent papers [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. The present paper provides new results of this kind, by taking as point of departure a solvable ODE characterizing the time-evolution of the N × N matrix U (t) which is different or more general than those previously employed to this end. These results are treated in the following section, including its subsections and Appendices A and B where all the equations solved in this paper are listed (hence, the reader wishing to get an immediate glance at them may immediately jump to these appendices). The last section, entitled "Outlook", tersely outlines further developments whose treatment is postponed to subsequent papers.

Solvable N -body problems
In this section we describe two classes of solvable N -body problems. The models of the first class are not new, hence their treatment is not elaborated beyond their identification; several models of the second class are new, hence they are fully dealt with.
In Subsection 2.1 we introduce a system of two matrix first-order ordinary differential equations (ODEs) characterizing the time-evolution of the two N × N matrices U ≡ U (t) and V ≡ V (t), and we indicate how the corresponding initial-value problem can be explicitly solved.
In Subsection 2.2 we show how -via the introduction of two appropriate ansätze -the N eigenvalues z n (t) of the matrix U (t) can be identified with the N coordinates of N unit-mass point particles whose time-evolution, generally taking place in the complex plane, is characterized by Newtonian equations of motion ("acceleration equal force", with nonlinear one-body and twobody forces). These models are thereby shown to be solvable. The first ansatz yields various models whose solvability was already known; hence their treatment is limited to the derivation of their equations of motion. The second ansatz yields several new models -as well as several previously known models -characterized by equations of motion with specific velocity-dependent one-body and two-body forces "of goldfish type" featuring several arbitrary coupling constants. They are all listed in Appendix A. The alternative class of many-body models obtained by changing the dependent variables from the N zeros of a monic polynomial of degree N to its N coefficients is discussed in Subsection 2.3; the corresponding equations of motion are listed in Appendix B.
As it is well known (see for instance [5], in particular Chapter 4, entitled "Solvable and/or integrable many-body problems in the plane, obtained by complexification"), models such as those described below, which describe the motions of N points in the complex z-plane, can be reformulated as N -body models describing the motion in a plane of N point-particles, the positions of which are identified by real 2-vectors in that plane. While we leave the elaboration of this connection to the interested reader, we feel justified by it to refer to our findings as describing "physical" (if not necessarily "realistic") many-body problems in the plane.

A solvable system of two matrix ODEs
In this subsection we discuss the following system of two matrix ODEs, satisfied by the two N × N matrices U and V : Lower case letters generally denote scalars. The 5 scalar parameters α, β, γ, η, ρ are timeindependent but a priori arbitrary, until specific restrictions on their values are explicitly mentioned. Superimposed dots always indicate time differentiations.
Remark 2.1. The factor β multiplying the term U 2 in the right-hand side of the first of the two matrix ODEs (2.1) could be eliminated (i.e., replaced by unity) by rescaling U (and by correspondingly replacing γ with an adjusted value, sayγ); then the factor η multiplying the term (U V + V U ), or the (adjusted) factorγ multiplying the term V , could also be eliminated (i.e., replaced by unity) by rescaling V . It is preferable not to do so in order to keep open the possibility to set one or the other of these parameters, η or γ, to zero (but we will never set both of them to zero, to avoid decoupling the time evolution of U (t) from that of V (t)). Note moreover that the introduction of a constant scalar term (implicitly multiplied by the N × N unit matrix I) in the right-hand side of the first equation would amount -up to a redefinition of other parameters -to adding to the matrix U the N × N unit matrix I multiplied by a new (constant) parameter, implying just a constant shift of all the eigenvalues of the matrix U , a trivial change not worth pursuing; while the introduction of two different parameters in front of the two terms U V and V U could be eliminated by adding, in the right-hand side of the first of the two matrix ODEs (2.1), the commutator [U, V ] of U and V times a convenient parameter, without any effect on the eigenvalues of U . If some of the parameters in the matrix ODEs (2.1) vanish this model may reduce to one of those that have been previously treated, see [6] (in particular its Section 4.2.2 entitled "Goldfishing") and the more recent papers [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]; we do not discard these models in the following, since we consider in any case worthwhile to exploit the unified treatment provided by the present approach based on (2.1). As for the possibility of replacing the right-hand side of the second matrix ODE (2.1) with a more general function of the matrix V , it is a possibility whose exploration is postponed to a future investigation (see Section 3). Solution of the system of two matrix ODEs (2.1). Clearly the solution of the second of the two matrix ODEs (2.1) reads To solve the first of the two matrix ODEs (2.1) it is convenient to set obtaining thereby for the N × N matrix F ≡ F (t) the following second-order linear matrix It is then plain via (2.2) that the general solution of this matrix ODE reads where F ± are two arbitrary constant N × N matrices and the two scalar functions f ± (t; v) (which of course become N × N matrices when the scalar v is replaced by the N × N matrix V 0 , see (2.5)) are two independent solutions of the scalar second-order linear ODË And it is easily seen that two independent solutions of this ODE are given by the following formulas: Here the (scalar!) function Φ(a; c; z) is the confluent hypergeometric function (see, for instance, [24]). Note that the formula (2.3) with (2.5) implies the following explicit solution of the initial-value problem for the matrix U (t): with the time-independent matrix C defined in terms of the initial data U (0) and V 0 ≡ V (0) (see (2.2)) as follows:  .7)) or of β (see (2.8a)), or values of α and ρ such that α/ρ is an integer (in which case the two functions f ± (t; v) do not provide two independent solutions of the ODE (2.6): see Section 6.7 of [24]) -which are clearly problematic (although the formulas may in these cases be reinterpreted via appropriate limiting procedures). We ignore hereafter these issues, even when listing below solvable N -body models a few of which belong to these problematic cases. Indeed in this paper we mainly limit our consideration to the identification of solvable N -body problems; these findings open the way to obtaining a rather detailed understanding of their actual behaviors (see for instance the following Remark 2.3), but such analyses exceed the scope of this paper: they should be done on a case-by-case basis, and shall perhaps be postponed to the moment when some of these N -body models evoke a specific, theoretical or applicative, interest.
Remark 2.3. If the parameter ρ is purely imaginary and the parameter α is real and negative, (of course with T real and nonvanishing, and i the imaginary unit, i 2 = −1), then clearly the matrix U (t) (see (2.8) with (2.7)) is asymptotically isochronous (i.e., asymptotically periodic with the period T independent of the initial data): indeed in this case, as t → +∞, with U + (t) given by the formula (2.8) with C = 0, hence (see (2.9a) and (2.7a)) it is periodic with period T , This observation entails of course the property of asymptotic isochrony for the solvable Nbody models identified below featuring parameters α and ρ which satisfy the conditions (2.9a).

Identif ication of solvable many-body models
The starting point is to introduce the N eigenvalues z n (t) of the matrix U (t), and the corresponding diagonalizing matrix R(t), by setting Here and hereafter indices such as n, m, , k run over the integers from 1 to N (unless otherwise indicated). Likewise we set and we introduce the N × N matrix M (t) by setting Remark 2.4. The diagonalizing matrix R(t) is defined up to right-multiplication by an arbitrary diagonal N × N matrix D(t), since the replacement of R(t) byR(t) = R(t)D(t) does not affect (2.10a). But it changes M (t) (see (2.10c)) intoM (t) = [R(t)] −1Ṙ (t) implying the following change of its diagonal elements: µ n (t) =⇒μ n (t) = µ n (t) +ḋ n (t)/d n (t), where the N quantities d n (t) are the, a priori arbitrary, elements of the diagonal matrix D(t). Hence we retain the privilege to assign at our convenience (see below) the diagonal elements µ n (t) of the matrix M (t).
It is now easily seen that via these assignments (2.10) the two matrix ODEs (2.1) get rewritten as follows: (2.11) Here and hereafter the notation [A, B] denotes the commutator of the two matrices A and B: Let us now look, componentwise, at the diagonal and off-diagonal elements of these two matrix ODEs.
The diagonal part of the first of the two ODEs (2.11) readṡ z n = αz n + βz 2 n + γy n + 2ηz n y n , (2.12a) implying y n =ż n − αz n − βz 2 n γ + 2ηz n . (2.12b) The off-diagonal part of the first of the two ODEs (2.11) reads The diagonal part of the second of the two ODEs (2.11) readṡ Hence, by time-differentiation of (2.12a) we see via (2.12b) and (2.14b) that the coordinates z n (t) satisfy the following system of N equations of motion of Newtonian type ("acceleration equal force", with one-body and two-body forces): z n = −αρz n − βρz 2 n + (α + ρ)ż n + 2βż n z n + 2ηż n ż n − αz n − βz 2 n γ + 2ηz n But in these equations of motion the role of "two-body coupling constants" is played by the quantities Y n Y n (with = n) which are in fact time-dependent. Indeed the time evolution of the off-diagonal elements Y nm of the matrix Y is determined by the off-diagonal part of the second of the two ODEs (2.11) which componentwise reaḋ yielding, via (2.10b), (2.10c), (2.12b), (2.13b) and a bit of algebra, So, in order that (2.15) become the N Newtonian equations of motion of a genuine N -body problem one must either provide some "physical interpretation" for the quantities Y n (with = n) -possibly in terms of internal degrees of freedom, an alternative we do not pursue in this paper -or find a way to "get rid" of these quantities, i.e. find a way to express them -if at all possible -via the N coordinates z m ≡ z m (t) and possibly also the N velocitiesż m ≡ż m (t). Previous experience [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23] suggest that two types of ansätze are the appropriate starting points to try and achieve this goal.

First ansatz
The first ansatz reads where we reserve the option to make a convenient assignment for the N functions g n of the coordinate z n , g n ≡ g n (z n ) ≡ g n [z n (t)]. The insertion of this ansatz in (2.16b) yields, after a bit of trivial algebra, Here and throughout the notation +((n → m)) indicates the addition of whatever comes before it, with the index n replaced by m. We now take advantage of the freedom (see Remark 2.4) to assign the diagonal elements µ n of the matrix M by setting and we moreover make the assignment with g an arbitrary constant (i.e.,ġ = 0). Thereby the system of N (N − 1) equations (2.18) gets reduced to the following, much simpler, system of only N algebraic equations: which amounts to the following 3 equations (recall that we exclude the uninteresting possibility that γ and η both vanish): This boils down to the following 3 possibilities: The corresponding solvable N -body models -obtained by inserting (2.17) with (2.20) and with these assignments of the parameters in the Newtonian equations of motion (2.15) -read, in the first 2 of these 3 cases -after conveniently setting with c an arbitrary nonvanishing constant -as follows: while in the third of these 3 cases they read But these are well-known solvable N -body problems, see for instance [5]. Hence we conclude that from the matrix system (2.1) no new solvable many-body models are obtained via the ansatz (2.17).

Second ansatz
We proceed then to consider a second ansatz, reading with g n and f n functions of the coordinate z n that we reserve to assign later: . Its insertion in (2.16b) yields, again after a bit of trivial algebra and now the assignment µ n = 0 (again justified by Remark 2.4), the following system of N (N − 1) second-order ODEs: To reduce this system we clearly must now set (2) are 3 constant parameters that we reserve to assign later. Thereby the system of N (N − 1) ODEs (2.27) gets transformed into the following system of (only) N ODEs: Note that to write this system in more compact form we employed a mixed notation, using sometimes the functions f n or f instead of their explicit expressions, see (2.28b).
To complete the task of ascertaining for which values of the (so far arbitrary) 8 constant parameters α, β, γ, η, ρ, f (0) , f (1) , f (2) the system of N (N − 1) ODEs (2.29) is satisfied we must utilize the equations of motion (2.15) which, via the ansatz (2.26) with (2.28), now read Comparison of this system of N ODEs to the system (2.29) yields the following system of N algebraic equations (note that -as it were, "miraculously" -the velocitiesż n have disappeared from these equations): It is now clear that, in order to satisfy this system -identically, i.e. for any values of the coordinates z n -one must set either So, in both cases, we hereafter only retain the freedom to assign the 2 constants a, b rather than the 3 constants Clearly in case (i) we get the following system of N algebraic equations: αρz n + βρz 2 n + [α + ρ + 2aη + bγ + 2(β + 2bη)z n ] aγ + (2aη + bγ)z n + 2bηz 2 n = 2η αz n + βz 2 n + (a + bz n )(γ + 2ηz n ) (a + bz n ), (2.34a) and likewise in case (ii) Hence in case (i) the following set of 4 nonlinear algebraic equations must be satisfied by the 7 parameters a, b, α, β, γ, η, ρ: Likewise in case (ii) the following set of 4 nonlinear algebraic equations must be satisfied by the 7 parameters a, b, α, β, γ, η, ρ: Then a trivial if rather tedious computation yields, in case (i) respectively in case (ii), the solutions reported in Table 2.1 respectively in Table 2.2 (note that for α = β = η = 0 these two cases coincide, so the corresponding results are only included in Table 2.2).   We therefore conclude that the many-body models of goldfish type characterized by the Newtonian equations of motion (2.30) (or, equivalently, (2.29)) are solvable provided either the quantities f n are expressed by (2.32a) and the 7 parameters a, b, α, β, γ, η, ρ are consistent with Table 2.1 or the quantities f n are expressed by (2.33a) and the 7 parameters a, b, α, β, γ, η, ρ are consistent with Table 2.2. There are therefore altogether 24 solvable models. Some of these models are however not new (in particular, when some parameters vanish); moreover in some cases the solution (2.8) with (2.7) of the matrix equation (2.1) is only valid in a limiting sense (see Remark 2.2). We leave these issues to whoever will be interested -possibly in specific, theoretical or applicative, contexts -in more detailed investigations of anyone of these models; our focus in this paper is rather in the unified treatment, and the simultaneous display (see the This table indicates the 5 sets of values to be assigned to the 7 parameters a, b, α, β, γ, η, ρ in (2.37) with (2.38). Asterisks indicate that the corresponding parameters can be assigned freely. Note that in every case there are 3 free parameters. appendices below), of all of them (except for the elimination, as already mentioned, of the cases with γ = η = 0). It is convenient to write the corresponding equations of motion in two different ways, depending whether the parameter η does or does not vanish.
If the parameter η vanishes, η = 0, the equations of motion read as follows: with the following assignments corresponding respectively to case (i) and case (ii).
In case (i) f n = aγ + bγz n (2.38) and the 7 parameters a, b, α, β, γ, η, ρ are restricted according to and the 6 parameters a, b, α, β, γ, ρ are restricted according to Table 2.4 (being the relevant subcase of Table 2.2). If the parameter η does not vanish, η = 0, the equations of motion are more conveniently written in terms of the dependent variables reading then as follows: x n =ẋ 2 n x n + λẋ n x n + βẋ n x n + ρ ẋ n + λ − µx n − βx 2 or equivalentlÿ and with the following assignments corresponding respectively to case (i) and case (ii). In case (i) and the 7 parameters a, b, α, β, γ, η, ρ are restricted according to Table 2.6 (being the relevant subcase of Table 2.2): These 24 Newtonian equations of motion are exhibited in Appendix A. Their solvable character is of course implied by the fact that the N coordinates z n (t) coincide with the N eigenvalues of the matrix U (t) (see (2.10a)). To ascertain the behavior of these solutions z n (t) one must in each case take account of the restrictions on the parameters characterizing these models, as detailed above, which are of course also relevant in order to identify the corresponding evolution of the matrix U (t): as implied by inserting in the explicit formula (2.8) with (2.7) -in addition to the parameters α, β, γ, η, ρ associated with the solvable N -body model under considerationthe expressions of the initial values U (0) and V (0) ≡ V 0 of the matrices U (t) and V (t) in terms of the N initial values z n (0) of the N coordinates and the N initial valuesż n (0) of the N velocities. To obtain these expressions it is useful to note that it is possible -and convenient -to assume that the diagonalizing matrix R(t) (see (2.10)) is initially just the N × N unit matrix I, implying initially (see (2.10)) The first of these two formulas provides the explicit expression of U (0) in terms of the initial data z n (0).
In the second formula the initial values y n (0) of the diagonal elements of the matrix V 0 in terms of the initial coordinates z n (0) and velocitiesż n (0) of the N particles read (see (2.12b)), while the off-diagonal elements Y nm (0) (with n = m) of the matrix V 0 are given by the ansatz (2.26) yielding with the quantities g n (0) and f n (0) given by the formulas (see (2.28)) g n (0) = 1 γ + 2ηz n (0) , with the appropriate assignments of the parameters γ, η, f (0) , f (1) and f (2) characterizing the various solvable many-body models, see above (in particular for f (0) , f (1) and f (2) see (2.38) or (2.39) or (2.42) or (2.43), as appropriate). Note moreover that the dyadic character of the off-diagonal part of the matrix V 0 = Y (0), see (2.45d), implies a simplification when one must compute functions of this matrix V 0 such as those appearing in the explicit expression (2.8) with (2.7) of U (t); this simplification becomes particularly significant when the matrix V 0 = Y (0) is altogether dyadic, Y nm (0) = v n v m , since for any dyadic matrix, say X nm = x n x m , there holds the simple formula where ϕ(x) is any (scalar) function for which the (matrix) expression ϕ(X) makes good sense. This simplification clearly happens iff and in case (ii) whenever (see (2.33) and Table 2.2)

47c)
Special cases and their (autonomous) isochronous variants. Certain special models among those identified above as solvable can be isochronized via the following change of dependent and independent variables, Here the quantities ζ n (τ ) are assumed to satisfy the Newtonian equations written above, see (2.37) with (2.38) or (2.39), of course with the new (complex ) independent variable τ replacing the time t; ω is an arbitrary real (for definiteness, positive) constant to which we associate the period T = 2π ω ; (2.49) and the number σ is adjusted so as to produce, for the dependent variables z n ≡ z n (t) (with the real independent variable t interpreted as "time") autonomous equations of motion (the special models providing the starting points for the application of this trick being appropriately selected to allow such an outcome). Since the application of this trick is by now quite standard (see, for instance, Section 2.1 entitled "The trick" of [6]), we dispense here from any detailed discussion of this approach and limit ourselves to reporting the results. This trick is only applicable to (2.37) with (2.38) in the very special cases with α = ρ = 0 and either γ = 0 or a = b = 0 (as long as one is only interested in getting autonomous equations of motion). Then the assignment σ = 1 yields the isochronous equations of motion z n = 3iωż n + 2ω 2 z n + 2βz n (ż n − iωz n ) + 2(ż n − iωz n ) Likewise the application of this trick to (2.37) with (2.39), again in the very special cases with α = ρ = 0 and either γ = 0 or a = b = 0, and again with the assignment σ = 1, yields the isochronous equations of motion z n = 3iωż n + 2ω 2 z n + 2βz n (ż n − iωz n ) Neither one of these two models is however new: see Examples 4.2.2-6 and 4.2.2-7 in [6]. An analogous treatment applied, mutatis mutandis, to (2.41) with (2.42) in the special case with ρ = 0, 2aη = bγ and 2αη = βγ (implying λ = 0), yields (again via the assignment σ = 1) the isochronous equations of motion Likewise an analogous treatment applied to (2.41) with (2.43) in the special case with ρ = 0 and either 2αη = βγ (implying λ = 0) and α = −2aη + bγ or γ = 0 (implying λ = 0) and α = 2aη, yields (again via the assignment σ = 1) the same isochronous equations of motion (up to the, merely notational, replacement of 2bη with 2bη − β).

A related class of solvable many-body models
In this subsection we consider the Newtonian equations of motion that obtain by identifying the N dependent variables of the models discussed above as the N zeros of a monic (timedependent) polynomial of degree N , and by then focussing on the time-evolution of the N coefficients of this polynomial. It is again convenient to treat separately the two cases with η = 0 and with η = 0. In the η = 0 case the starting point are the equations of motion (2.37) with (2.38) or (2.39). We then introduce the time-dependent (monic) polynomial ψ(z, t) whose zeros are the N eigenvalues z n (t) of the N × N matrix U (t): The last of these formulas introduces the N coefficients c m ≡ c m (t) of the monic polynomial ψ(z, t); of course it implies that these coefficients are related to the zeros z n (t) as follows:

51c)
and so on. The fact that the initial-value problem associated with the time evolution of the N coordinates z n can be solved by algebraic operations implies that the same solvable character can be attributed to the time evolution of the monic polynomial ψ(z, t) and of the N coefficients c m (t). The procedure to obtain the equations of motion satisfied by the N coefficients c m (t) from the N equations of motion satisfied by the N zeros is tedious but standard; a key role in this development are the identities reported, for instance, in Appendix A of [6] (but note that there are two misprints in these formulas: in equation (A.8k) the term (N + 1) inside the square brackets should instead read (N − 3); in equation (A.8l) the term N 2 inside the square brackets should instead read N (N − 2) -these misprints have been corrected in the recent paperback version of this monograph [6]). Here we limit our presentation to reporting the final result.
The equation characterizing the time evolution of the monic polynomial ψ(z, t) implied by the Newtonian equations of motion (2.37) with (2.28b) reads as follows: (2.53c) where of course the quantities f (0) , f (1) , f (2) should be expressed in terms of the other parameters as implied by (2.28b) with (2.38) or (2.39), as the case may be. Note that, while this equation, (2.52), satisfied by the function ψ(z, t) (where of course subscripted variables denote partial differentiations) might seem a linear PDE, it is in fact a nonlinear functional equation, because some of its coefficients, see (2.53), depend on the quantities c 1 and c 2 which themselves depend on ψ, indeed clearly (see (2.51)) where we used the shorthand notation ψ (j) (z, t) to denote the j-th partial derivative with respect to the variable z of ψ(z, t), Likewise, the equation characterizing the time evolution of the monic polynomial φ (x, t) implied by the Newtonian equations of motion (2.41) via the following assignment (analogous to (2.51b)), The equations of motion of Newtonian type satisfied by the N coefficients c m (t) which obtain from (2.52) hence correspond to the Newtonian equations of motion (2.37) read as follows: where c n vanishes for n < 0 and for n > N while c 0 = 1 (see (2.51)), and of course the coefficients p (j) and q (j) k are defined by (2.53). Again, this system of ODEs might seem linear, but it is in fact nonlinear because some of its coefficients depend on the dependent variables c 1 and c 2 , see (2.53). The more explicit version of these equations of motion that obtain by expressing the various coefficients in terms of the free parameters are listed in Appendix B. They are of course just as solvable as the Newtonian equations of motion satisfied by the N coordinates x n (t), see Appendix A, to which they correspond via (2.51); and in particular whenever the parameter ρ is imaginary and the parameter α is real and negative they are asymptotically isochronous with period T (see Remark 2.3).
Likewise, the equations of motion of Newtonian type satisfied by the N coefficients c m (t) which obtain from (2.57) hence correspond to the Newtonian equations of motion (2.41) read as follows: where of course again c n vanishes for n < 0 and for n > N while c 0 = 1 (see (2.55)) and of course the coefficients p (j) and q (j) k are now defined by (2.58). Again, this system of ODEs might seem linear, but it is in fact nonlinear because some of its coefficients depend on the dependent variables c 1 , c N −1 and c N , see (2.58). The more explicit version of these equations of motion that obtains by expressing the various coefficients in terms of the free parameters are listed in Appendix B. They are of course just as solvable as the Newtonian equations of motion satisfied by the N coordinates x n (t), see Appendix A, to which they correspond via (2.55); and in particular whenever the parameter ρ is imaginary and the parameter α is real and negative they are asymptotically isochronous with period T (see Remark 2.3).
Special cases and their (autonomous) isochronous variants. Certain special models among those identified above (in this subsection) as solvable can be isochronized by an analogous trick to that employed at the end of the preceding Subsection 2.2. One route to this end takes as starting point the isochronized systems of Newtonian equations of motion (2.50) and applies to them the same procedure employed above to obtain the equations of motion (2.59) with (2.53) and (2.60) with (2.58). An equivalent procedure is to apply to certain special subcases of these systems of ODEs, (2.59) with (2.53) and (2.60) with (2.58), the following change of dependent and independent variables: Here the quantities χ n (τ ) are assumed to satisfy the systems of ODEs written above, see (2.59) with (2.53) and (2.60) with (2.58), of course with the new (complex ) independent variable τ replacing the time t; ω is an arbitrary real (for definiteness, positive) constant to which we associate the period T, see (2.49); and the number σ is adjusted so as to produce autonomous ODEs for the new dependent variables c m ≡ c m (t) (with the real independent variable t interpreted as "time": the special models providing the starting points for the application of this trick being appropriately selected in order to allow such an outcome). Since the application of this trick is quite standard, we dispense here from any detailed discussion of this approach and limit ourselves to reporting the results. The ODEs that follow from (2.59) (with α = ρ = η = 0, f (0) = f (1) = f (2) = 0 and σ = 1) read as follows: (2.62a) those that follow from (2.59) (with α = ρ = η = 0, f (0) = f (1) = 0, f (2) = −β and σ = 1) read instead as follows Neither one of these two isochronous many-body problems is new. The analogous results that follows from (2.60) are instead generally new. There are then two sets of cases. The first set of isochronous models obtain from the assignment σ = 1 and read with the following restriction on the parameters: The second set of isochronous models obtain from the assignment σ = −1 and read with the following restrictions on the parameters:

Outlook
Results analogous, but somewhat more general, than those reported in this paper can be obtained by an analogous treatment based on a somewhat more general -but still solvable -system of two N × N matrix ODEs than (2.1), such as, for instance, which contains the 2 additional scalar constants ρ 0 and ρ 2 (and clearly reduces to (2.1) for ρ 0 = ρ 2 = 0). These developments will be reported in subsequent papers. Finally, let us recall that Diophantine findings can be obtained from a nonlinear autonomous isochronous dynamical system by investigating its behavior in the infinitesimal vicinity of its equilibria. The relevant equations of motion become then generally linear, but they of course retain the properties to be autonomous and isochronous. For a system of linear autonomous ODEs, the property of isochrony implies that all the eigenvalues of the matrix of its coefficients are integer numbers (up to a common rescaling factor). When the linear system describes the behavior of a nonlinear autonomous system in the infinitesimal vicinity of its equilibria, these matrices can generally be explicitly computed in terms of the values at equilibrium of the dependent variables of the original, nonlinear model. In this manner nontrivial Diophantine findings and conjectures have been discovered and proposed: see for instance the review of such developments in Appendix C (entitled "Diophantine findings and conjectures") of [6]. Analogous results obtained by applying this approach to the isochronous systems of autonomous nonlinear ODEs introduced above -and in subsequent papers -will be reported if they turn out to be novel and interesting.

B Second appendix
In this appendix we list the second series of 24 Newtonian equations of motion whose solvable character has been demonstrated in this paper; they correspond to those reported in Appendix A via the transformation among the N zeros z n and the N coefficients c m of a monic polynomial, see (2.51) and (2.55). In each case the parameters they feature (such as a, b, α, β, γ, η, ρ, as the case may be) are arbitrary constants; the (assigned) values of the other ones of these parameters (which also characterize the time-evolution of the solutions of these equations, see (2.8)), are also reported. Let us emphasize that if the parameter ρ is an imaginary number and the parameter α is real and negative, the corresponding many-body problem is asymptotically isochronous with period T , see Remark 2.3; and that isochronous many-body models are characterized by the 4 Newtonian equations of motion (2.62), (2.63) and (2.64) displayed at the end of Subsection 2.3. Let us also mention again that the equations of motion reported below are not all new; in particular all those that are linear are of course not new. Let us recall that it is always assumed that c n = 0 for n < 0 and for n > N , and c 0 = 1. η = 0, case (i) (5 models, corresponding to Table 2.3): (1) a = β = 0, α = −bγ: (2) β = ρ = 0, α = −bγ: