Polynomials with Multiple Zeros and Solvable Dynamical Systems including Models in the Plane with Polynomial Interactions

TL;DR

This paper develops a comprehensive method to analyze time-dependent polynomials with multiple zeros, leading to the discovery of new solvable dynamical systems, especially focusing on degree 4 polynomials with specific zero multiplicities.

Contribution

It introduces a general approach for polynomials with arbitrary zeros and multiplicities, extending previous methods to more complex cases and identifying new solvable dynamical systems.

Findings

Extended the analysis to polynomials with multiple zeros of arbitrary multiplicity.

Identified new classes of solvable dynamical systems for degree 4 polynomials.

Provided a framework for analyzing nongeneric polynomials with multiple zeros.

Abstract

The interplay among the time-evolution of the coefficients and the zeros of a generic time-dependent (monic) polynomial provides a convenient tool to identify certain classes of solvable dynamical systems. Recently this tool has been extended to the case of nongeneric polynomials characterized by the presence, for all time, of a single double zero; and subsequently significant progress has been made to extend this finding to the case of polynomials featuring a single zero of arbitrary multiplicity. In this paper we introduce an approach suitable to deal with the most general case, i. e. that of a nongeneric time-dependent polynomial with an arbitrary number of zeros each of which features, for all time, an arbitrary (time-independent) multiplicity. We then focus on the special case of a polynomial of degree 4 featuring only 2 different zeros and, by using a recently introduced additional twist of this approach, we thereby identify many new classes of solvable dynamical systems.

Equations

$${\dot{x}}_n=P^{\left(n\right)}\left(x_1,x_2\right),n=1,2,$$

$$p_N\left(z;t\right)=z^N+\sum _{m=1}^N\left[y_m\left(t\right)z^{N-m}\right]=\prod _{n=1}^N\left[z-x_n\left(t\right)\right],$$

$$y_m\left(t\right)={\left(-1\right)}^m\sum _{1\le n_1<n_2<...<n_m\le N}^N\left\{\prod _{\ell =1}^M\left[x_{n_{\ell }}\left(t\right)\right]\right\},m=1,2,...,N;$$

$${\dot{x}}_n=-{\left[\prod _{\ell =1,\ell \ne n}^N\left(x_n-x_{\ell }\right)\right]}^{-1}\sum _{m=1}^N\left[{\dot{y}}_m{\left(x_n\right)}^{N-m}\right]$$

$${\dot{x}}_n=P^{\left(n\right)}\left(x_1,x_2\right),n=1,2,$$

$$P_M\left(z;t\right)=z^M+\sum _{m=1}^M\left[y_m\left(t\right)z^{M-m}\right]=\prod _{n=1}^N\left\{{\left[z-x_n\left(t\right)\right]}^{{\mu }_n}\right\}.$$

$$M=\sum _{n=1}^N\left({\mu }_n\right).$$

$$y_1\left(t\right)=-\sum _{n=1}^N\left[{\mu }_n{x}_n\left(t\right)\right],y_M\left(t\right)={\left(-1\right)}^M\prod _{n=1}^N\left\{{\left[x_n\left(t\right)\right]}^{{\mu }_n}\right\},$$

$$P_{M,t}\left(z;t\right)=\sum _{m=1}^M\left[{\dot{y}}_m\left(t\right)z^{M-m}\right]$$

$$\sum _{m=1}^{M-\mu }\left\{{\dot{y}}_m\left(t\right)\left[\frac{\left(M-m\right)!}{\left(M-m-\mu \right)!}\right]{\left[x_n\left(t\right)\right]}^{M-m-\mu }\right\}=0,$$

$$\mu$$

$$\sum _{m=1}^{M-{\mu }_n+1}\left\{{\dot{y}}_m\left(t\right)\left[\frac{\left(M-m\right)!}{\left(M-m-{\mu }_n+1\right)!}\right]{\left[x_n\left(t\right)\right]}^{M-m-{\mu }_n+1}\right\}$$

$$n$$

$${\dot{x}}_n\left(t\right)=-{\left\{{\mu }_n!\prod _{\ell =1,\ell \ne n}^N{\left[x_n\left(t\right)-x_{\ell }\left(t\right)\right]}^{{\mu }_{\ell }}\right\}}^{-1}\cdot$$

$$\cdot \sum _{m=1}^{M-{\mu }_n+1}\left\{{\dot{y}}_m\left(t\right)\left[\frac{\left(M-m\right)!}{\left(M-m-{\mu }_n+1\right)!}\right]{\left[x_n\left(t\right)\right]}^{M-m-{\mu }_n+1}\right\},$$

$$n=1,2,...,N.$$

$$\sum _{n=1}^N\left({\mu }_m-1\right)=M-N$$

$${\dot{x}}_1=-{\left[{\mu }_1!{\left(x_1-x_2\right)}^{{\mu }_2}\right]}^{-1}\sum _{m=1}^{1+{\mu }_2}\left\{{\dot{y}}_m\left[\frac{\left({\mu }_1+{\mu }_2-m\right)!}{\left({\mu }_2-m+1\right)!}\right]{\left(x_1\right)}^{{\mu }_2-m+1}\right\},$$

$${\dot{x}}_2=-{\left[{\mu }_2!{\left(x_2-x_1\right)}^{{\mu }_1}\right]}^{-1}\sum _{m=1}^{1+{\mu }_1}\left\{{\dot{y}}_m\left[\frac{\left({\mu }_1+{\mu }_2-m\right)!}{\left({\mu }_1-m+1\right)!}\right]{\left(x_2\right)}^{{\mu }_1-m+1}\right\}.$$

$$y_1$$

$$y_3$$

$${\dot{y}}_1{\left(x_1\right)}^3+{\dot{y}}_2{\left(x_1\right)}^2+{\dot{y}}_3{x}_1+{\dot{y}}_4=0,$$

$$3{\dot{y}}_1{\left(x_1\right)}^2+2{\dot{y}}_2{x}_1+{\dot{y}}_3=0$$

$${\dot{x}}_1=-\frac{3x_1{\dot{y}}_1+{\dot{y}}_2}{3\left(x_1-x_2\right)},{\dot{x}}_2=\frac{{\left(x_2\right)}^3{\dot{y}}_1+{\left(x_2\right)}^2{\dot{y}}_2+x_2{\dot{y}}_3+{\dot{y}}_4}{{\left(x_1-x_2\right)}^3}.$$

$$\tilde{m}=1,2;\tilde{m}=1,3;\tilde{m}=1,4;\tilde{m}=2,3;\tilde{m}=2,4;\tilde{m}=3,4,$$

$$m=3,4;m=2,4;m=2,3;m=1,4;m=1,3;m=1,2.$$

$${\dot{x}}_1=-\frac{3x_1{\dot{y}}_1+{\dot{y}}_2}{3\left(x_1-x_2\right)},{\dot{x}}_2=\frac{\left(2x_1+x_2\right){\dot{y}}_1+{\dot{y}}_2}{\left(x_1-x_2\right)},$$

$${\dot{x}}_1=-\frac{3{\left(x_1\right)}^2{\dot{y}}_1-{\dot{y}}_3}{6x_1\left(x_1-x_2\right)},{\dot{x}}_2=\frac{x_1\left(x_1+2x_2\right){\dot{y}}_1-{\dot{y}}_3}{2x_1\left(x_1-x_2\right)},$$

$${\dot{x}}_1=-\frac{{\left(x_1\right)}^3{\dot{y}}_1+{\dot{y}}_4}{3{\left(x_1\right)}^2\left(x_1-x_2\right)},{\dot{x}}_2=\frac{{\left(x_1\right)}^2{x}_2{\dot{y}}_1+{\dot{y}}_4}{{\left(x_1\right)}^2\left(x_1-x_2\right)},$$

$${\dot{x}}_1=\frac{x_1{\dot{y}}_2+{\dot{y}}_3}{3x_1\left(x_1-x_2\right)},{\dot{x}}_2=-\frac{x_1\left(x_1+2x_2\right){\dot{y}}_2+\left(2x_1+x_2\right){\dot{y}}_3}{3{\left(x_1\right)}^2\left(x_1-x_2\right)},$$

$${\dot{x}}_1=\frac{{\left(x_1\right)}^2{\dot{y}}_2-3{\dot{y}}_4}{6{\left(x_1\right)}^2\left(x_1-x_2\right)},{\dot{x}}_2=-\frac{{\left(x_1\right)}^2{x}_2{\dot{y}}_2-\left(2x_1+x_2\right){\dot{y}}_4}{2{\left(x_1\right)}^3\left(x_1-x_2\right)},$$

$${\dot{x}}_1=-\frac{x_1{\dot{y}}_3+3{\dot{y}}_4}{3{\left(x_1\right)}^2\left(x_1-x_2\right)},{\dot{x}}_2=\frac{x_1{x}_2{\dot{y}}_3+\left(x_1+2x_2\right){\dot{y}}_4}{{\left(x_1\right)}^3\left(x_1-x_2\right)}.$$

$${\dot{y}}_1=f_1\left(y_1,y_2\right),{\dot{y}}_2=f_2\left(y_1,y_2\right).$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$$3x{f}_1\left(-4x,6x^2\right)+f_2\left(-4x,6x^2\right)=0;$$

$$f_1\left(y_1,y_2\right)={\alpha }_0+{\alpha }_1{y}_2,f_2\left(y_1,y_2\right)={\beta }_0{y}_1+{\beta }_1{\left(y_1\right)}^3,$$

$${\alpha }_0=\frac{4{\beta }_0}{3},{\alpha }_1=\frac{32{\beta }_1}{9};$$

Full text

462PolsWithMultiplZerosAndSolvDynSyst181216Rev

Polynomials with Multiple Zeros and Solvable Dynamical Systems including Models in the Plane with Polynomial Interactions

Francesco Calogeroa,b,1 and Farrin Payandeha,c,2

a Physics Department, University of Rome ”La Sapienza”, Rome, Italy

b INFN, Sezione di Roma 1

c Department of Physics, Payame Noor University (PNU), PO BOX 19395-3697 Tehran, Iran

1 francesco.calogero@roma1.infn.it, francesco.calogero@uniroma1.it

2 f_payandeh@pnu.ac.ir, farrinpayandeh@yahoo.com

Abstract

The interplay among the time-evolution of the coefficients $y_{m}\left(t\right)$ and the zeros $x_{n}\left(t\right)$ of a generic time-dependent (monic) polynomial provides a convenient tool to identify certain classes of solvable dynamical systems. Recently this tool has been extended to the case of nongeneric polynomials characterized by the presence, for all time, of a single double zero; and subsequently significant progress has been made to extend this finding to the case of polynomials featuring a single zero of arbitrary multiplicity. In this paper we introduce an approach suitable to deal with the most general case, i. e. that of a nongeneric time-dependent polynomial with an arbitrary number of zeros each of which features, for all time, an arbitrary (time-independent) multiplicity. We then focus on the special case of a polynomial of degree $4$ featuring only $2$ different zeros and, by using a recently introduced additional twist of this approach, we thereby identify many new classes of solvable dynamical systems of the following type:

[TABLE]

with $P^{\left(n\right)}\left(x_{1},x_{2}\right)$ two polynomials in the two variables $x_{1}\left(t\right)$ and $x_{2}\left(t\right)$.

1 Introduction

Notation 1.1. Hereafter $t$ generally denotes time (the real independent variable); (partial) derivatives with respect to time are denoted by a superimposed dot, or in some case by appending as a subscript the independent variable $t$ preceded by a comma; all dependent variables such as $x$, $y$, $z$ (often equipped with subscripts) are generally assumed to be complex numbers, unless otherwise indicated (it shall generally be clear from the context which of these and other quantities depend on the time $t$, as occasionally—but not always—explicitly indicated); parameters such as $a$, $b,$ $c,$ $\alpha,$ $\beta,$ $\gamma,$ $A$, etc. (often equipped with subscripts) are generally time-independent complex numbers; and indices such as $n$, $m$, $j$, $\ell$ are generally positive integers, with ranges explicitly indicated or clear from the context. $\blacksquare$

Some time ago the idea has been exploited to identify dynamical systems which can be solved by using as a tool the relations between the time evolutions of the coefficients and the zeros of a generic time-dependent polynomial [1]. The basic idea of this approach is to relate the time-evolution of the $N$ zeros $x_{n}\left(t\right)$ of a generic time-dependent polynomial $p_{N}\left(z;t\right)$ of degree $N$ in its argument $z,$

[TABLE]

to the time-evolution of its $N$ coefficients $y_{m}\left(t\right)$. Indeed, if the time evolution of the $N$ coefficients $y_{m}\left(t\right)$ is determined by a system of Ordinary Differential Equations (ODEs) which is itself solvable, then the corresponding time-evolution of the $N$ zeros $x_{n}\left(t\right)$ is also solvable, via the following $3$ steps: (i) given the initial values $x_{n}\left(0\right),$ the corresponding initial values $y_{m}\left(0\right)$ can be obtained from the explicit formulas expressing the coefficients $y_{m}$ of a polynomial in terms of its zeros $x_{n}$ reading (for all time, hence in particular at $t=0$)

[TABLE]

(ii) from the $N$ values $y_{m}\left(0\right)$ thereby obtained, the $N$ values $y_{m}\left(t\right)$ are then evaluated via the—assumedly solvable—system of ODEs satisfied by the $N$ coefficients $y_{m}\left(t\right)$; (iii) the $N$ values $x_{n}\left(t\right)$—i. e., the $N$ solutions of the dynamical system satisfied by the $N$ variables $x_{n}\left(t\right)$—are then determined as the $N$ zeros of the polynomial, see (1), itself known at time $t$ in terms of its $N$ coefficients $y_{m}\left(t\right)$ (the computation of the zeros of a known polynomial being an algebraic operation; of course generally explicitly performable only for polynomials of degree $N\leq 4$).

Remark 1-1. In this paper the term ”solvable” generally characterizes systems of ODEs the initial-values of which are ”solvable by algebraic operations”—possibly including quadratures yielding implicit solutions, generally also requiring the evaluation of parameters via algebraic operations. And let us emphasize that, because of the algebraic but nonlinear character of the relations between the zeros and the coefficients of a polynomial, it is clear that, even to relatively trivial evolutions of the $N$ coefficients $y_{m}\left(t\right)$ of a time-dependent polynomial, there correspond much less trivial evolutions of its $N$ zeros $x_{n}\left(t\right)$. On the other hand the fact that a time evolution is algebraically solvable has important implications, generally excluding that it can be ”chaotic”, indeed in some cases allowing to infer important qualitative features of the time evolution, such as the property to be isochronous or asymptotically isochronous (see for instance [2] [3] [4]). $\blacksquare$

The viability of this technique to identify solvable dynamical systems depends of course on the availability of an explicit method to relate the time-evolution of the $N$ zeros of a polynomial to the corresponding time-evolution of its $N$ coefficients. Such a method was indeed provided in [1], opening the way to the identification of a vast class of algebraically solvable dynamical systems (see also [5] and references therein); but that approach was essentially restricted to the consideration of linear time evolutions of the coefficients $y_{m}\left(t\right)$.

A development allowing to lift this quite strong restriction emerged relatively recently [6], by noticing the validity of the identity

[TABLE]

which provides a convenient explicit relationship among the time evolutions of the $N$ zeros $x_{n}\left(t\right)$ and the $N$ coefficients $y_{m}\left(t\right)$ of the generic polynomial (1). This allowed a major enlargement of the class of algebraically solvable dynamical systems identifiable via this approach: for many examples see [7] and references therein.

Remark 1-2. Analogous identities to (3) have been identified for higher time-derivatives [8] [9] [7]; but in this paper we restrict our treatment to dynamical systems characterized by first-order ODEs, postponing the treatment of dynamical systems characterized by higher-order ODEs (see Section 6). $\blacksquare$

A new twist of this approach was then provided by its extension to nongeneric polynomials featuring—for all time—multiple zeros. The first step in this direction focussed on time-dependent polynomials featuring for all time a single double zero [10]; and subsequently significant progress has been made to treat the case of polynomials featuring a single zero of arbitrary multiplicity [11]. In Section 2 of the present paper a convenient method is provided which is suitable to treat the most general case of polynomials featuring an arbitrary number of zeros each of which features an arbitrary multiplicity. While all these developments might appear to mimic scholastic exercises analogous to the discussion among medieval scholars of how many angels might dance simultaneously on the point of a needle, they do indeed provide new tools to identify new dynamical systems featuring interesting time evolutions (including systems displaying remarkable behaviors such as isochrony or asymptotic isochrony: see for instance [10] [11]); dynamical systems which—besides their intrinsic mathematical interest—are quite likely to play significant roles in applicative contexts.

Such developments shall be reported in future publications. In the present paper we focus on another twist of this approach to identify new solvable dynamical systems which was introduced quite recently [12]. It is again based on the relations among the time-evolution of the coefficients and the zeros of time-dependent polynomials [6] [7] with multiple roots (see [10], [11] and above); but (as in [12]) by restricting attention to such polynomials featuring only $2$ zeros. Again, this might seem such a strong limitation to justify the doubt that the results thereby obtained be of much interest. But the effect of this restriction is to open the possibility to identify *algebraically solvable *dynamical models characterized by the following system of $2$ ODEs,

[TABLE]

with $P^{\left(n\right)}\left(x_{1},x_{2}\right)$ two polynomials in the two dependent variables $x_{1}\left(t\right)$ and $x_{2}\left(t\right)$; hence systems of considerable interest, both from a theoretical and an applicative point of view (see [12] and references quoted there). This development is detailed in the following Section 3 by treating a specific example. In Section 4 we report—without detailing their derivation, which is rather obvious on the basis of the treatment provided in Section 3—many other such solvable models (see (4); but in some cases the right-hand side of these equations are not quite polynomial); and a simple technique allowing additional extensions of these models—making them potentially more useful in applicative contexts—is outlined in Section 5, by detailing its applicability in a particularly interesting case. Hence researchers primarily interested in applications of such systems of ODEs might wish to take first of all a quick look at these $2$ sections.

Finally, Section 6 outlines future developments of this research line; and some material useful for the treatment provided in the body of this paper is reported in $2$ Appendices.

2 Properties of nongeneric time-dependent polynomials featuring $N$

zeros, each of arbitrary multiplicity

In this Section 2 we focus on time-dependent (monic) polynomials featuring for all time $N$ different zeros $x_{n}\left(t\right)$, each of which with the arbitrarily assigned (of course time-independent) multiplicity $\mu_{n}$. They are of course defined as follows:

[TABLE]

It is plain that there exist explicit formulas—generalizing (2)—expressing the $M$ coefficients $y_{m}\left(t\right)$ in terms of the $N$ zeros $x_{n}\left(t\right)$; for instance clearly

[TABLE]

and see other examples below.

It is also plain that, while $N$ zeros $x_{n}$ and their multiplicities $\mu_{n}$ can be arbitrarily assigned in order to define the polynomial (5), this is not the case for the $M$ coefficients $y_{m}$: generally—for any given assignment of the $N$ multiplicities $\mu_{n}$—only $N$ of them can be arbitrarily assigned, thereby determining (via algebraic operations) the $N$ zeros $x_{n}$ and the remaining $M-N$ other coefficients $y_{m}.$

Remark 2-1. The generic polynomial (1), of degree $N$ and featuring $N$ different zeros $x_{n}$ and $N$ coefficients $y_{m}$, generally implies that the set of its $N$ coefficients $y_{m}$ is an $N$-vector $\vec{y}=\left(y_{1},...,y_{N}\right),$ while the set of its $N$ zeros $x_{n}$ is instead an unordered set of $N$ numbers $x_{n}$. This however is not quite true in the case of a time-dependent generic polynomial which features—as those generally considered in this paper—a continuous time-dependence of its coefficients and zeros; then the set of its $N$ zeros $x_{n}\left(0\right)$ at the initial time $t=0$ should be generally considered an unordered set, but for all subsequent time, $t>0$, the set of its $N$ zeros $x_{n}\left(t\right)$ is an ordered set, the assignment of the index $n$ to $x_{n}\left(t\right)$ being no more arbitrary but rather determined by continuity in $t$ (at least provided during the time evolution no collision of two or more different zeros occur, in which case the identities of these zeros get to some extent lost because their identities may be exchanged, becoming undetermined).

The situation is quite different in the case of a nongeneric polynomial such as (5): then zeros having different multiplicities are intrinsically different, for instance if all the multiplicities $\mu_{n}$ are different among themselves, $\mu_{n}\neq\mu_{\ell}$ if $n\neq\ell$, then clearly the set of the $N$ zeros $x_{n}$ is an ordered set (hence an $N$-vector).

We trust the reader to understand these rather obvious facts and therefore hereafter we refrain from any additional discussion of these issues. $\blacksquare$

Our task now is to identify—for the special class of nongeneric polynomials (5)— equivalent relations to the identities (3), to be then used in order to identify new solvable dynamical systems.

The first step is to time-differentiate once the formula (5), getting the relations

[TABLE]

Remark 2-2. Hereafter, in order to avoid clattering our presentation with unessential details, we occasionally make the convenient assumption that all the numbers $\mu_{n}$ be different among themselves; the diligent reader shall have no difficulty to understand how the treatment can be extended to include cases in which this simplifying assumption does not hold—indeed in the specific examples discussed below we will include in our treatment also cases in which this simplification is not valid, taking appropriate care of such cases. And we also assume—without loss of generality—that the numbers $\mu_{n}$ are ordered in decreasing order, $\mu_{n}\geq\mu_{n+1}$. $\blacksquare$

Our next step is to $z$-differentiate $\mu$ times the above formulas, firstly with $\mu=0,1,2,...,\mu_{n}-2$ and secondly with $\mu=\mu_{n}-1$; and then set $z=x_{n}$ (for each value of $n=1,2,...,N$). There clearly thereby obtain the following formulas:

[TABLE]

[TABLE]

The second set, (8b), yields the following $N$ expressions of the time-derivatives of the $N$ zeros $x_{n}\left(t\right)$ in terms of the time-derivatives of the $M$ coefficients $y_{m}\left(t\right)$:

[TABLE]

The first set, (8a), consists of linear relations among the $M$ time-derivatives $\dot{y}_{m}\left(t\right)$ of the $M$ coefficients $y_{m}\left(t\right)$: and it is easily seen, via (5b), that there are altogether

[TABLE]

such relations. So one can select $N$ quantities $\dot{y}_{m}\left(t\right)$—let us hereafter call them $\dot{y}_{\tilde{m}}\left(t\right)$—and compute, from the $M-N$ linear equations (8a), all the other $M-N$ quantities $\dot{y}_{m}\left(t\right)$ with $m\neq\tilde{m}$ as linear expressions in terms of these selected $N$ quantities $\dot{y}_{\tilde{m}}\left(t\right)$. The goal of expressing the $N$ time-derivatives $\dot{x}_{n}$ as linear equations—somehow analogous to the identities (3)—in terms of the $N$ time-derivatives $\dot{y}_{\tilde{m}}\left(t\right)$ of $N$, arbitrarily selected, coefficients $y_{\tilde{m}}\left(t\right)$ is thereby finally achieved. Indeed the task of expressing the $M-N$ quantities $\dot{y}_{m}\left(t\right)$ with $m\neq\tilde{m}$ in terms of the $N$ quantities $\dot{y}_{\tilde{m}}\left(t\right)$—and of course the $N$ zeros $x_{n}\left(t\right)$—can in principle be implemented explicitly as it amounts to solving the $M-N$ linear equations (8) for the $M-N$ unknowns $\dot{y}_{m}\left(t\right)$, with the $N$ quantities $\dot{y}_{\tilde{m}}\left(t\right)$ playing there the role of known quantities; clearly implying that the resulting expressions of the $M-N$ quantities $\dot{y}_{m}\left(t\right)$ are linear functions of the $N$ quantities $\dot{y}_{\tilde{m}}\left(t\right)$. And the insertion of these linear expressions of the $M-N$ quantities $\dot{y}_{m}\left(t\right)$ (with $m\neq\tilde{m}$) in terms of the $N$ quantities $\dot{y}_{\tilde{m}}\left(t\right)$ in the $N$ formulas (9) fulfils our goal.

The actual implementation of this development must of course be performed on a case-by-case basis, see below. In the special case with only one multiple zero—and if moreover the indices $\tilde{m}$ are assigned their first $N$ values, i. e. $\tilde{m}=1,2,...,N$—these results shall reproduce the results of the path-breaking paper [11], which was confined to the treatment of this special case.

The outcomes of these developments are detailed, in the special case with $N=2$ and $M=4$, in the following Section 3; for the motivation of this drastic restriction see below.

3 The $N=2,$ $M=4$ case

In the special case with $N=2$ the formula (9) simplifies, reading (see (5b))

[TABLE]

Let us moreover restrict attention to the case with $M=4,$ as the case with $M=3$ (implying $\mu_{1}=2,$ $\mu_{2}=1$: see Remark 2-2) has been already discussed in [10] and [12] and the case with $M=4$ is sufficiently rich (see below) to deserve a full paper.

In the case with $M=4$ there are $2$ possible assignments of the $2$ parameters $\mu_{n}$: (i) $\mu_{1}=3,$ $\mu_{2}=1;$ (ii) $\mu_{1}=\mu_{2}=2$ (see Remark 2-2, and note that we now include also the case with $\mu_{1}=\mu_{2}$).

3.1 Case (i): $\mu_{1}=3,$ $\mu_{2}=1$

In this case (5a) clearly implies the following expressions of the $4$ coefficients $y_{m}\left(t\right)$ in terms of the $2$ zeros $x_{n}\left(t\right)$:

[TABLE]

Remark 3.1-1. Note that these formulas imply that $x_{1}$ and $x_{2}$ can be computed (in fact, explicitly!) from $y_{m_{1}}$ and $y_{m_{2}}$—with $m_{1}=1,2,3$ and $m_{1}<m_{2}=2,3,4$) by solving an algebraic equation of degree $m_{2}$. $\blacksquare$

The corresponding equations (8a) read

[TABLE]

And the formulas (11) read

[TABLE]

There are now $6$ possible different assignments for the indices $\tilde{m}$:

[TABLE]

Let us list below the 6 corresponding versions of the ODEs (14):

[TABLE]

Next, let us focus to begin with—in order to explain our approach—on the first, (16a), of the $6$ formulas (16). Assume moreover that the $2$ quantities $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ evolve according to the following solvable system of ODEs:

[TABLE]

It is then clear—via the identities (12)—that we can conclude that the dynamical system

[TABLE]

While this is in itself an interesting result—to become more significant for explicit assignments of the $2$ functions $f_{1}\left(y_{1},y_{2}\right)$ and $f_{2}\left(y_{1},y_{2}\right)$ (see below)—an additional interesting development emerges if—following the approach of [12]—we now assume the $2$ functions $f_{1}\left(y_{1},y_{2}\right)$ and $f_{2}\left(y_{1},y_{2}\right)$ to be both polynomial in their $2$ arguments and moreover such that

[TABLE]

a restriction that is clearly sufficient to guarantee that the right-hand sides of the equations of motion (21) become polynomials in the $2$ dependent variables $x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ (since the numerators in the right-hand sides of the $2$ ODEs (21) are then both polynomials in the variables $x_{1}$ and $x_{2}$ which vanish when $x_{1}=x_{2}=x$ and which therefore contain the factor $x_{1}-x_{2}$).

An representative example of such functions is

[TABLE]

The conclusion is then that the dynamical system

[TABLE]

Remark 3.1-2. An equivalent—indeed more direct—way to identify the solvable dynamical system (21) as corresponding to the solvable dynamical system

[TABLE]

(see (17) and (20)), is via the relations

[TABLE]

3.2 Case (ii): $\mu_{1}=\mu_{2}=2$

In this case

[TABLE]

Remark 3.2-1. Of course a remark completely analogous to Remark 3.1-1 holds in this case as well. $\blacksquare$

The corresponding equations (8a) read

[TABLE]

and proceeding as above one easily obtains, for the $6$ assignments (15), the following systems of $2$ ODEs:

[TABLE]

[TABLE]

Hence, to the system (17), one now associates again the requirement (19); and—by making again the assignment (20a) for the system of evolution equations satisfied by $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$—one identifies again the restriction (20b), thereby concluding—via (25a)—that the polynomial system

[TABLE]

where now $a=-\beta_{0}/3$ and $b=4\beta_{1}/9$, is solvable. And the explicit solution is then quite analogous (up to simple modifications of some parameters) to that described (after eq. (21)) in the preceding Subsection 3.1.

4 Other solvable systems of 2 nonlinearly-coupled ODEs

identified via the technique described in Section 3

In this Section 4 we report a list of solvable systems of $2$ nonlinearly coupled first-order ODEs satisfied by the $2$ dependent variables $x_{1}\left(t\right)$ and $x_{2}\left(t\right)$; in each case we identify the corresponding *solvable *system of $2$ ODEs satisfied by $2$ variables $y_{\tilde{m}}\left(t\right)$ (for these, and other, notations used below see Section 3); indeed, to help the reader mainly interested in the solvable character of one of the following systems we also specify below on a case-by-case basis the information which allows to solve that specific system (we do so even at the cost of minor repetitions). Note that the majority of these models feature equations of type (4), but in a few cases the right-hand sides of these ODEs are not quite polynomial. And let us recall that in this Section 4 parameters such as $a,$ $b,$ $c$ (possibly equipped with indices) are arbitrary numbers (possibly complex).

Remark 4-1. Most of the models reported below are characterized by evolution equations of the following kind:

[TABLE]

Remark 4-2. In the following $44$ subsections we list as many solvable systems of $2$ nonlinearly-coupled first-order ODEs, most of them with polynomial right-hand sides, and we indicate how each of them can be solved. The presentation of all these models is made so as to facilitate the utilization of these findings by practitioners only interested in one of these models (or its generalization, see Section 5). Note however that not all these models are different among themselves: indeed, some feature identical equations of motion—although the method to solve them might seem different. This is demonstrated by the following self-evident identification of the following equations of motion: (32)$\equiv$(54), (33)$\equiv$(55), (40)$\equiv$(56)$\equiv$(70), (41)$\equiv$(57)$\equiv$(71), (48)$\equiv$(62)$\equiv$(68), (49)$\equiv$(63)$\equiv$(69). So in fact the list below contains only $34$ different systems of $2$ nonlinearly-coupled first-order differential equations for the $2$ time-dependent variables $x_{1}\left(t\right)$ and $x_{2}\left(t\right)$. $\blacksquare$

4.1 Model 4.(i)1.2a

[TABLE]

4.2 Model 4.(ii)1.2a

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ by (24a); and the variables $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ evolve according to (82), the explicit solution of which is given by the relevant formulas in Subsection Case A.1 of Appendix A with $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=2,$ $L=1,$ $\alpha_{0}=a,$ $\alpha_{1}=\left(3/4\right)b,$ $\beta_{0}=2a,$ $\beta_{1}=4b.$

4.3 Model 4.(i)1.2b

[TABLE]

This model (30a) (with $a_{0}=0$) is actually a special case of the more general model

[TABLE]

again with $x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ related to $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ by (12) and the variables $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ evolving according to (87) with $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=2,$ $L$ an arbitrary positive integer, $\alpha_{0}=-4a_{0},$ $\alpha_{\ell}=a_{\ell}+4b_{\ell},$ $\beta_{\ell}=2a_{\ell},$ $\gamma_{0}=-3a_{0},$ $\gamma_{\ell}=3b_{\ell}.$

4.4 Model 4.(ii)1.2b

[TABLE]

This model (31a) is actually a special case of the more general model

[TABLE]

again with $x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ related to $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ by (24a) and the variables $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ evolving according to (87) with $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=2,$ $L$ an arbitrary positive integer, $\alpha_{0}=-4a_{0},$ $\alpha_{\ell}=a_{\ell}+2b_{\ell},$ $\beta_{\ell}=2a_{\ell},$ $\gamma_{0}=-3a_{0},$ $\gamma_{\ell}=\left(3/2\right)b_{\ell}.$

4.5 Model 4.(i)1.2c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ by (12); and the variables $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=2,$ $\tilde{m}_{2}=1,$ $L=3,$ $\alpha_{0}=0,$ $\alpha_{1}=2a,$ $\alpha_{2}=2b/3,$ $\alpha_{3}=2c/9,$ $\beta_{1}=a,$ $\beta_{2}=b/3,$ $\beta_{3}=c/9,$ $\gamma_{\ell}=0.$

4.6 Model 4.(ii)1.2c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ by (24a); and the variables $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=2,$ $\tilde{m}_{2}=1,$ $L=3,$ $\alpha_{0}=0,$ $\alpha_{1}=2a,$ $\alpha_{2}=2b,$ $\alpha_{3}=2c,$ $\beta_{1}=a,$ $\beta_{2}=b,$ $\beta_{3}=c,$ $\gamma_{\ell}=0.$

4.7 Model 4.(i)1.2d

[TABLE]

4.8 Model 4.(ii)1.2d

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ by (24a); and the variables $y_{1}\left(t\right)$ and $y_{2}\left(t\right)$ evolve according to (95a), the explicit solution of which is given by the relevant formulas in Subsection Case A.3.1 of Appendix A with $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=2,$ $\alpha_{0}=-4a,$ $\alpha_{1}=8b,$ $\beta_{0}=-3a,$ $\beta_{1}=\left(9/4\right)b$. Note that this is the model treated in detail in Subsection 3.1.2, see (26).

4.9 Model 4.(i)1.3a

[TABLE]

4.10 Model 4.(ii)1.3a

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{3}\left(t\right)$ by (24a); and the variables $y_{1}\left(t\right)$ and $y_{3}\left(t\right)$ evolve according to (82), the explicit solution of which is given by the relevant formulas in Subsection Case A.1 of Appendix A with $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=3,$ $L=1,$ $\alpha_{0}=a,$ $\alpha_{1}=-b/8,$ $\beta_{0}=3a,$ $\beta_{1}=-6b.$

4.11 Model 4.(i)1.3b

[TABLE]

4.12 Model 4.(ii)1.3b

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{3}\left(t\right)$ by (24a); and the variables $y_{1}\left(t\right)$ and $y_{3}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=3,$ $L=3,$ $\alpha_{0}=-\left(16/3\right)b_{0},$ $\alpha_{\ell}=-\left(-2\right)^{-\ell}\left(a_{\ell}+2^{5-\ell}b_{\ell}\right)/3,$ $\beta_{\ell}=-\left(-2\right)^{-\ell-1}a_{\ell},$ $\gamma_{0}=-b_{0},$ $\gamma_{\ell}=-\left(-2\right)^{-\ell}b_{\ell},$ $\ell=1,2,3.$ Note that the right-hand sides of these $2$ ODEs, (38), are both polynomial iff the single parameter $b_{0}$ vanishes, $b_{0}=0$.

4.13 Model 4.(i)1.3c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{3}\left(t\right)$ by (12); and the variables $y_{1}\left(t\right)$ and $y_{3}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=3,$ $\tilde{m}_{2}=1,$ $L=3,$ $\alpha_{0}=0,$ $\alpha_{1}=3a,$ $\alpha_{2}=-3b,$ $\alpha_{3}=3c,$ $\beta_{1}=a,$ $\beta_{2}=-b,$ $\beta_{3}=c,$ $\gamma_{\ell}=0.$

4.14 Model 4.(ii)1.3c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{3}\left(t\right)$ by (24a); and the variables $y_{1}\left(t\right)$ and $y_{3}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A) with $\tilde{m}_{1}=3,$ $\tilde{m}_{2}=1,$ $L=3,$ $\alpha_{0}=0,$ $\alpha_{1}=3a,$ $\alpha_{2}=-3b/2,$ $\alpha_{3}=3c/4,$ $\beta_{1}=a,$ $\beta_{2}=-b/2,$ $\beta_{3}=c/4,$ $\gamma_{\ell}=0.$

4.15 Model 4.(i)1.3d

[TABLE]

4.16 Model 4.(ii)1.3d

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{3}\left(t\right)$ by (24a); and the variables $y_{1}\left(t\right)$ and $y_{3}\left(t\right)$ evolve according to (95a), the explicit solution of which is given by the relevant formulas in **Subsection Case A.3.2 **of Appendix A with $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=3,$ $\alpha_{0}=-8a,$ $\alpha_{1}=-16b,$ $\beta_{0}=-\left(3/2\right)a,$ $\beta_{1}=-\left(3/16\right)b$. Note that the right-hand sides of the $2$ ODEs (43) are polynomial only if $a=0.$

4.17 Model 4.(i)1.4a

[TABLE]

4.18 Model 4.(ii)1.4a

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{4}\left(t\right)$ by (24a); and the variables $y_{1}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (82), the explicit solution of which is given by the relevant formulas in **Subsection Case A.1 **of Appendix A with $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=4,$ $L=1,$ $\alpha_{0}=a,$ $\alpha_{1}=b/16,$ $\beta_{0}=4a,$ $\beta_{1}=2^{6}b=64b$.

4.19 Model 4.(i)1.4b

[TABLE]

4.20 Model 4.(ii)1.4b

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{4}\left(t\right)$ by (24a); and the variables $y_{1}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=4,$ $L=3,$ $\alpha_{0}=16b_{0};$ $\alpha_{\ell}=\left(-2\right)^{1-\ell}a_{\ell-1}+\left(-2\right)^{4-\ell}b_{\ell},$ $\beta_{\ell}=\left(-2\right)^{3-\ell}a_{\ell-1},$ $\ell=1,2,3;$ $\gamma_{\ell}=\left(-2\right)^{-2-\ell}b_{\ell},$ $\ell=0,1,2,3.$ Note that the right-hand sides of these $2$ ODEs, (47), are both polynomial only if all the $4$ parameters $b_{\ell}$ vanish, $b_{\ell}=0,$ $\ell=1,2,3$.

4.21 Model 4.(i)1.4c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{4}\left(t\right)$ by (12); and the variables $y_{1}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=4,$ $\tilde{m}_{2}=1,$ $L=3,$ $\alpha_{0}=0,$ $\alpha_{1}=4a,$ $\alpha_{2}=4b,$ $\alpha_{3}=4c,$ $\beta_{1}=a,$ $\beta_{2}=b,$ $\beta_{3}=c,$ $\gamma_{\ell}=0.$

4.22 Model 4.(ii)1.4c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{4}\left(t\right)$ by (24a); and the variables $y_{1}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=4,$ $\tilde{m}_{2}=1,$ $L=3,$ $\alpha_{0}=0,$ $\alpha_{1}=4a,$ $\alpha_{2}=4b,$ $\alpha_{3}=4c,$ $\beta_{1}=a,$ $\beta_{2}=b,$ $\beta_{3}=c,$ $\gamma_{\ell}=0.$

4.23 Model 4.(i)1.4d

[TABLE]

4.24 Model 4.(ii)1.4d

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{1}\left(t\right)$ and $y_{4}\left(t\right)$ by (24a); and the variables $y_{1}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (95a), the explicit solution of which is given by the relevant formulas in Subsection Case A.3.3 of Appendix A with $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=4,$ $\alpha_{0}=2^{4}a=16a;$ $\alpha_{1}=2^{8}b=256b,$ $\beta_{0}=2^{-2}a=a/4,$ $\beta_{1}=2^{-6}b=b/64.$ Note that the right-hand sides of these $2$ ODEs, (50), are not polynomial.

4.25 Model 4.(i)2.3a

[TABLE]

4.26 Model 4.(ii)2.3a

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{2}\left(t\right)$ and $y_{3}\left(t\right)$ by (24a); and the variables $y_{2}\left(t\right)$ and $y_{3}\left(t\right)$ evolve according to (82), the explicit solution of which is given by the relevant formulas in Subsection Case A.1 of Appendix A with $\tilde{m}_{1}=2,$ $\tilde{m}_{2}=3,$ $L=1,$ $\alpha_{0}=2a,$ $\alpha_{1}=2b,$ $\beta_{0}=3a,$ $\beta_{1}=\left(81/2\right)b.$ Note that the right-hand sides of these ODEs is not polynomial.

4.27 Model 4.(i)2.3b

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{2}\left(t\right)$ and $y_{3}\left(t\right)$ by (12); and the variables $y_{2}\left(t\right)$ and $y_{3}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=2,$ $\tilde{m}_{2}=3,$ $L=3,$ $\alpha_{0}=0~{},$ $\alpha_{1}=2a,$ $\alpha_{2}=\left(2/3\right)b,$ $\alpha_{3}=\left(2/9\right)c,$ $\beta_{1}=3a,$ $\beta_{2}=b,~{}\beta_{3}=c/3,~{}\gamma_{\ell}=0.$

4.28 Model 4.(ii)2.3b

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{2}\left(t\right)$ and $y_{3}\left(t\right)$ by (24a); and the variables $y_{2}\left(t\right)$ and $y_{3}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=2,$ $\tilde{m}_{2}=3,$ $L=3,$ $\alpha_{0}=0~{},$ $\alpha_{1}=2a,$ $\alpha_{2}=2b,$ $\alpha_{3}=2c,$ $\beta_{1}=3a,$ $\beta_{2}=3b,~{}\beta_{3}=3c,~{}\gamma_{\ell}=0.$

4.29 Model 4.(i)2.3c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{2}\left(t\right)$ and $y_{3}\left(t\right)$ by (12); and the variables $y_{2}\left(t\right)$ and $y_{3}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=3,$ $\tilde{m}_{2}=2,$ $L=3,$ $\alpha_{0}=0~{},$ $\alpha_{1}=3a,$ $\alpha_{2}=-3b,$ $\alpha_{3}=3c,$ $\beta_{1}=2a,$ $\beta_{2}=-2b,$ $\beta_{3}=3c,$ $\gamma_{\ell}=0.$

4.30 Model 4.(ii)2.3c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{2}\left(t\right)$ and $y_{3}\left(t\right)$ by (24a); and the variables $y_{2}\left(t\right)$ and $y_{3}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=3,$ $\tilde{m}_{2}=2,$ $L=3,$ $\alpha_{0}=0~{},$ $\alpha_{1}=3a,$ $\alpha_{2}=-\left(3/2\right)b,$ $\alpha_{3}=\left(3/4\right)c,$ $\beta_{1}=2a,$ $\beta_{2}=-b,$ $\beta_{3}=c/2,$ $\gamma_{\ell}=0.$

4.31 Model 4.(i)2.4a

[TABLE]

4.32 Model 4.(ii)2.4a

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{2}\left(t\right)$ and $y_{4}\left(t\right)$ by (24a); and the variables $y_{2}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (82), the explicit solution of which is given by the relevant formulas in Subsection Case A.1 of Appendix A with $\tilde{m}_{1}=2,$ $\tilde{m}_{2}=4,$ $L=2,$ $\alpha_{0}=2a,$ $\alpha_{1}=2b,$ $\alpha_{2}=-2c,$ $\beta_{0}=4a,$ $\beta_{1}=144b,$ $\beta_{2}=-2^{6}3^{4}c=-5184c.$ Note that the right-hand sides of these ODEs are not polynomial, except for the trivial case with $b=c=0.$

4.33 Model 4.(i)2.4b

[TABLE]

4.34 Model 4.(ii)2.4b

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{2}\left(t\right)$ and $y_{4}\left(t\right)$ by (24a); and the variables $y_{2}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=2,$ $\tilde{m}_{2}=4,$ $L=3,$ $\alpha_{0}=36b_{0};$ $\alpha_{\ell}=2a_{\ell-1}+36b_{\ell},$ $\beta_{\ell}=4a_{\ell-1},$ $\ell=1,2,3;$ $\gamma_{\ell}=2b_{\ell},$ $\ell=0,1,2,3.$ Note that the right-hand sides of these $2$ ODEs, (61), are both polynomial only if all the $4$ parameters $b_{\ell}$ vanish, $b_{\ell}=0,$ $\ell=0,1,2,3$.

4.35 Model 4.(i)2.4c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{2}\left(t\right)$ and $y_{4}\left(t\right)$ by (12); and the variables $y_{2}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=4,$ $\tilde{m}_{2}=2,$ $L=3,$ $\alpha_{0}=0~{},$ $\alpha_{1}=4a,$ $\alpha_{2}=4b,$ $\alpha_{3}=4c,$ $\beta_{1}=2a,$ $\beta_{2}=2b,~{}\beta_{3}=2c,~{}\gamma_{\ell}=0.$

4.36 Model 4.(ii)2.4c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{2}\left(t\right)$ and $y_{4}\left(t\right)$ by (24a); and the variables $y_{2}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=4,$ $\tilde{m}_{2}=2,$ $L=3,$ $\alpha_{0}=0~{},$ $\alpha_{1}=4a,$ $\alpha_{2}=4b,$ $\alpha_{3}=4c,$ $\beta_{1}=2a,$ $\beta_{2}=2b,~{}\beta_{3}=2c,~{}\gamma_{\ell}=0.$

4.37 Model 4.(i)2.4d

[TABLE]

4.38 Model 4.(ii)2.4d

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{2}\left(t\right)$ and $y_{4}\left(t\right)$ by (24a); and the variables $y_{2}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (95a), the explicit solution of which is given by the relevant formulas in Subsection Case A.3.1 of Appendix A with $\tilde{m}_{1}=2,$ $\tilde{m}_{2}=4,$ $\alpha_{0}=-36a;$ $\alpha_{1}=-1296b,$ $\beta_{0}=-2a,$ $\beta_{1}=-9b.$ Note that the right-hand sides of these $2$ ODEs, (65), are not polynomial.

4.39 Model 4.(i)3.4a

[TABLE]

4.40 Model 4.(ii)3.4a

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{3}\left(t\right)$ and $y_{4}\left(t\right)$ by (24a); and the variables $y_{3}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (82), the explicit solution of which is given by the relevant formulas in Subsection Case A.1 of Appendix A with $\tilde{m}_{1}=3,$ $\tilde{m}_{2}=4,$ $L=1,$ $\alpha_{0}=3a,$ $\alpha_{1}=\left(3/16\right)b,$ $\beta_{0}=4a,$ $\beta_{1}=64b.$

4.41 Model 4.(i)3.4b

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{3}\left(t\right)$ and $y_{4}\left(t\right)$ by (12); and the variables $y_{3}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=4,$ $\tilde{m}_{2}=3,$ $L=3,$ $\alpha_{0}=0~{},$ $\alpha_{1}=4a,$ $\alpha_{2}=4b,$ $\alpha_{3}=4c,$ $\beta_{1}=3a,$ $\beta_{2}=3b,~{}\beta_{3}=3c,~{}\gamma_{\ell}=0.$

4.42 Model 4.(ii)3.4b

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{3}\left(t\right)$ and $y_{4}\left(t\right)$ by (24a); and the variables $y_{3}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=4,$ $\tilde{m}_{2}=3,$ $L=3,$ $\alpha_{0}=0~{},$ $\alpha_{1}=4a,$ $\alpha_{2}=4b,$ $\alpha_{3}=4c,$ $\beta_{1}=3a,$ $\beta_{2}=3b,$ $\beta_{3}=3c,$ $\gamma_{\ell}=0.$

4.43 Model 4.(i)3.4c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{3}\left(t\right)$ and $y_{4}\left(t\right)$ by (12); and the variables $y_{3}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=3,$ $\tilde{m}_{2}=4,$ $L=3,$ $\alpha_{0}=0~{},$ $\alpha_{1}=3a,$ $\alpha_{2}=-3b,$ $\alpha_{3}=3c,$ $\beta_{1}=4a,$ $\beta_{2}=-4b,~{}\beta_{3}=4c,~{}\gamma_{\ell}=0.$

4.44 Model 4.(ii)3.4c

[TABLE]

$x_{1}\left(t\right)$ and $x_{2}\left(t\right)$ are related to $y_{3}\left(t\right)$ and $y_{4}\left(t\right)$ by (24a); and the variables $y_{3}\left(t\right)$ and $y_{4}\left(t\right)$ evolve according to (87), the explicit solution of which is given by the relevant formulas in Subsection Case A.2 of Appendix A with $\tilde{m}_{1}=3,$ $\tilde{m}_{2}=4,$ $L=3,$ $\alpha_{0}=0~{},$ $\alpha_{1}=3a,$ $\alpha_{2}=-\left(3/2\right)b,$ $\alpha_{3}=\left(3/4\right)c,$ $\beta_{1}=4a,$ $\beta_{2}=-2b,$ $\beta_{3}=c,$ $\gamma_{\ell}=0.$

5 Extensions

In this Section 5 we tersely indicate the possibility to generalize the class of solvable models listed in the preceding Section 4, by outlining the procedure to do so in just one case, that detailed in **Subsection 4.3 **(see (30a)), in fact just the special case of it with $a_{0}=a_{1}=a_{3}=b_{1}=b_{3}=0$, so that its equations of motion read as follows:

[TABLE]

Remark 5-1. Note that the right-hand sides of the $2$ ODEs (72) are homogeneous polynomials of second degree, the coefficients of which satisfy of course the condition

[TABLE]

as implied by Remark 4-1. Moreover—as clearly implied by (72a)—the $2$ homogeneous second-degree polynomials in the right-hand sides of the 2 ODEs characterizing this model feature a common zero: they both vanish when $X=0$, namely when $x_{2}=-3x_{1}$. $\blacksquare$

Analogous extensions of other models treated in this paper shall be performed by practitioners interested in these systems of ODEs in the context of specific applications (see Section 6).

Remark 5-2. Note that, via a by now well-known trick (see, for instance, [7]) corresponding to the following time-dependent change of both independent and dependent variables,

[TABLE]

the autonomous system (72) gets replaced by the following, also autonomous, system:

[TABLE]

Here $a$ is an arbitrary (time-independent) parameter; and note that if this parameter $a$ is purely imaginary, $\mathop{\mathrm{R}e}\left[a\right]=0,$ $\mathop{\mathrm{I}m}\left[a\right]\neq 0,$ then this dynamical system (75) is generally doubly periodic**;** or even—if $a_{2}/b_{2}$ is a real rational number—isochronous, namely then all its solutions are completely periodic with a period (an integer multiple of $T=2\pi/\left|a\right|$) independent of the initial data: see Remark A.2-1 and, if need be, [7] [3]. $\blacksquare$

Because of this remarkable fact, in the remaining part of this Section 5 we limit, for simplicity, consideration to the special case (72), by investigating its extension which obtains via the following linear reshuffle of the $2$ dependent variables $x_{1}\left(t\right)$ and $x_{2}\left(t\right)$:

[TABLE]

[TABLE]

with the $6$ parameters $a_{n\ell}$, $n=1,2,$ $\ell=1,2,3$ explicitly expressed in terms of the $4$ arbitrary parameters $A_{nm},$ $n=1,2,$ $m=1,2,$ and the $2$ arbitrary parameters $a_{2}$ and $b_{2}$ (see (72c)) as follows:

[TABLE]

Remark 5.3. The fact that the $6$ parameters $a_{n\ell}$ which characterize the system (77) can be (explicitly!) expressed, see (78), in terms of $6$ a priori arbitrary parameters—the $4$ parameters $A_{nm},$ see (76), and the $2$ parameters $a_{2}$ and $b_{2}$ (see (72))—might seem to imply that this system (77) can be reduced by algebraic operations to the algebraically solvable system (72)— hence that it is itself *algebraically solvable—*for any generic assignment of its $6$ parameters $a_{n\ell}$, $n=1,2,$ $\ell=1,2,3$. That this is not the case is however implied by the observation that the property of the system (72)—to feature in the right-hand sides of its $2$ ODEs $2$ polynomials themselves featuring a common zero (see Remark 5-1)—is then clearly also featured by the generalized system (77) (we like to thank François Leyvraz for this very useful observation). Hence only (at most) $5$ of the $6$ parameters $a_{n\ell}$ ($n=1,2;$ $\ell=1,2,3$) can be arbitrarily assigned, since these $6$ parameters are constrained by the condition

[TABLE]

which is easily seen to correspond to the requirement that the right-hand sides of the $2$ ODEs (77) (with $n=1,2$) feature a common zero. $\blacksquare$

Remark 5-4. Let us finally emphasize that the trick reported in Remark 5-1 is just as applicable to the more general system (77), implying—via the ansatz

[TABLE]

The relevance of this dynamical system, (80b), in many applicative context is exemplified by too many contributions to allow reporting a full bibliography; we record here just one such paper which lists $11$ references and contains the remarkable assertion that the system (80b) ”is not solvable explicitly except in certain simple cases” [13].

6 Outlook

In this final Section 6 we tersely outline future developments of the findings reported in this paper.

There is of course the possibility to treat cases with $M>4$ (see Section 3).

There is the possibility to iterate the procedure leading to the identification of new solvable systems (as described in this paper): see for this kind of development [14] and Chapter 6 of [7].

Another natural development is to treat analogous dynamical systems evolving in *discrete *rather than continuous time. For progress in this direction see [15].

Another extension is to treat systems characterized by second-order rather than first-order differential equations, including models characterized by Newtonian equations of motion (”accelerations equal forces”); and in the cases in which these equations of motion are derivable from a Hamiltonian, an additional interesting development is the treatment of the corresponding time-evolutions in the context of quantal rather than classical mechanics.

And yet another extension is to Partial Differential Equations (PDEs) rather than ODEs.

There is finally the vast universe of applications, including to cases in which the systems of evolution equations can be shown—via their solvability—to feature remarkable properties such as isochrony [2] [3] or asymptotic isochrony [4].

7 Acknowledgements

FP likes to thank the Physics Department of the University of Rome ”La Sapienza” for the hospitality from April to November 2018 (during her sabbatical), when the results reported in this paper were obtained. FC likes to thank Robert Conte and François Leyvraz for very useful discussions in the context of the Gathering of Scientists on ”Integrable systems and beyond” hosted by the Centro Internacional de Ciencias (CIC) in Cuernavaca, Mexico, from November 19th to December 14th, 2018.

8 Appendix A: Three useful classes of solvable systems of $2$

nonlinear first-order ODEs for the $2$ variables $y_{\tilde{m}}\left(t\right)$

The findings reported in this Appendix A are not new; they are displayed here to facilitate the reader of the new findings reported in the body of this paper.

Notation A-1. In this Appendix A we indicate with the notation $y_{\tilde{m}_{1}}\left(t\right)$ and $y_{\tilde{m}_{2}}\left(t\right)$—with $\tilde{m}_{1,2}=1,2,3,4$ and $\tilde{m}_{1}\neq\tilde{m}_{2}$ (and for the significance of the superimposed tilde see the last part of Section 2)—the $2$ dependent variables which satisfy the ”solvable” system of $2$ nonlinearly-coupled ODEs

[TABLE]

with the $2$ functions $f_{\tilde{m}_{1}}\left(y_{\tilde{m}_{1}},y_{\tilde{m}_{2}}\right)$ and $f_{\tilde{m}_{2}}\left(y_{\tilde{m}_{1}},y_{\tilde{m}_{2}}\right)$ assigned—conveniently for our treatment in this paper (see Section 2 above)—so that the system (81) is ”solvable”. The precise meaning of the term ”solvable” shall be clear from the following.

In this Appendix A $\alpha_{\ell},$ $\beta_{\ell},$ $\gamma_{\ell}$ are a priori arbitrary time-independent parameters, and $L$ is an a priori arbitrary nonnegative integer.

The selection of the specific systems of $2$ ODEs considered below is of course motivated by the treatment in the body of this paper, see in particular Sections 2 and 3. $\blacksquare$

8.1 Case A.1

[TABLE]

Each of these $2$ ODEs can be integrated via one quadrature, that can be performed explicitly after some purely algebraic operations. Indeed, to integrate the first of these $2$ ODEs one must first of all identify—via an algebraic operation—the $L\tilde{m}_{2}+1$ zeros $\bar{y}_{n}$ (assumed below, for simplicity, to be all different among themselves) of the polynomial in its right-hand side,

[TABLE]

next one must identify the $L\tilde{m}_{2}+1$ ”residues” $r_{n}$ defined by the ”partial fraction decomposition” formula

[TABLE]

—another algebraic operation, indeed one that can be performed explicitly; and finally one integrates the ODE getting the (generally implicit; but not always, see below) result

[TABLE]

which characterizes the solution $y_{\tilde{m}_{1}}\left(t\right)$ corresponding to the initial datum $y_{\tilde{m}_{1}}\left(0\right).$

Of course an analogous procedure characterizes—for the second ODE (82)—the solution $y_{\tilde{m}_{2}}\left(t\right)$ corresponding to the initial datum $y_{\tilde{m}_{2}}\left(0\right).$

For $L=1$ the initial-value problem for these $2$ ODEs can be solved explicitly, since the solution of the initial-value problem for the ODE

[TABLE]

8.2 Case A.2

[TABLE]

The solution of the first of these $2$ ODEs, (87a), has been already discussed above, see Subsection Case A.1; hence in this Subsection Case A.2 we need to consider only the second ODE (87b).   And since we are mainly interested in the case when the right-hand side of this ODE is polynomial, we shall limit our consideration below to the $3$ subcases with $\tilde{m}_{1}=1$ and to the single case with $\tilde{m}_{1}=2,$ $\tilde{m}_{2}=4;$ except in the special case with all parameters $\gamma_{\ell}$ vanishing, $\gamma_{\ell}=0$, which we treat separately firstly (since it is an intermediate step to solve the more general case).

In this special case the ODE (87b) reads

[TABLE]

And it is then easily seen that the solution of the initial-value problem of the (more general) ODE (87b) reads

[TABLE]

More explicit solutions can be easily obtained in the following cases:

[TABLE]

[TABLE]

Remark A.2-1. In particular, it is easily seen that the system of ODEs (72) discussed in Section 5 (see **Subsection 4.3 **with $a_{0}=a_{1}=a_{3}=b_{1}=b_{3}=0$) implies that the system of ODEs (87), which then reads

[TABLE]

[TABLE]

[TABLE]

8.3 Case A.3

[TABLE]

The most convenient technique to solve the system of $2$ ODEs (95a) is by noticing to begin with that it implies for the variable $y_{\tilde{m}_{1}}\left(t\right)$ the second-order ODE

[TABLE]

which is an autonomous Newtonian equation of motion (”acceleration equal force”) and is of course solvable by quadratures (as discussed in more detail in the following special cases).

And of course once $y_{\tilde{m}_{1}}\left(t\right)$ is known, $\dot{y}_{\tilde{m}_{2}}=\left[\dot{y}_{\tilde{m}_{1}}\left(t\right)-\alpha_{0}\right]/\alpha_{1}$ is as well known (see the first of the $2$ ODEs (95a)).

8.3.1 Case A.3.1

In this case $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=2,$ or $\tilde{m}_{1}=2,$ $\tilde{m}_{2}=4,$ so that (95b) reads

[TABLE]

here and below $w\left(t\right)$ stands for $y_{\tilde{m}_{1}}\left(t\right)$ respectively $y_{\tilde{m}_{2}}\left(t\right),$ in the $2$ cases $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=2,$ respectively $\tilde{m}_{1}=2,$ $\tilde{m}_{2}=4$.

It is easily seen that the solution of the initial-value problem of this ODE reads as follows:

[TABLE]

[TABLE]

where of course $u\left(t\right)$ stands for $y_{2}\left(t\right)$ respectively $y_{4}\left(t\right).$

Here sn$\left(z,k\right),$ cn$\left(z,k\right),$ dn$\left(z,k\right)$ are the $3$ standard Jacobian elliptic functions.

8.3.2 Case A.3.2

In this case $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=3,$ so that (95b) reads

[TABLE]

It is thus seen that in this case $y_{1}\left(t\right)$ is a hyperelliptic function.

8.3.3 Case A.3.3

In this case $\tilde{m}_{1}=1,$ $\tilde{m}_{2}=4,$ so that (95b) reads

[TABLE]

It is thus again seen that in this case $y_{1}\left(t\right)$ is a hyperelliptic function.

9 Appendix B

In this Appendix B we display—for the case $M=4$ and $N=2$—the expressions of the time-derivatives $\dot{y}_{\tilde{m}}\left(t\right)$ of the coefficients $y_{\tilde{m}}\left(t\right)$ in terms of the time-derivatives $\dot{y}_{m}\left(t\right)$ of the coefficients $y_{m}\left(t\right)$ and of the zeros $x_{n}\left(t\right)$, for the $2$ cases with $\mu_{1}=3,$ $\mu_{2}=1$ respectively $\mu_{1}=\mu_{2}=2$ (for this terminology, see Section 2).

In case (i), with $\mu_{1}=3,$ $\mu_{2}=1,$ these relations read as follows:

[TABLE]

In case (ii), with $\mu_{1}=\mu_{2}=2,$ these relations read as follows:

[TABLE]