Discrete-Time Goldfishing

The original continuous-time"goldfish"dynamical system is characterized by two neat formulas, the first of which provides the $N$ Newtonian equations of motion of this dynamical system, while the second provides the solution of the corresponding initial-value problem. Several other, more general, solvable dynamical systems"of goldfish type"have been identified over time, featuring, in the right-hand ("forces") side of their Newtonian equations of motion, in addition to other contributions, a velocity-dependent term such as that appearing in the right-hand side of the first formula mentioned above. The solvable character of these models allows detailed analyses of their behavior, which in some cases is quite remarkable (for instance isochronous or asymptotically isochronous). In this paper we introduce and discuss various discrete-time dynamical systems, which are as well solvable, which also display interesting behaviors (including isochrony and asymptotic isochrony) and which reduce to dynamical systems of goldfish type in the limit when the discrete-time independent variable $\ell=0,1,2,...$ becomes the standard continuous-time independent variable $t$, $0\leq t<\infty $.


Introduction
The original "goldfish" dynamical system [1,2] is characterized by the system of N Newtonian equations of motion ⋆ This paper is a contribution to the Proceedings of the Conference "Symmetries and Integrability of Difference Equations (SIDE-9)" (June 14-18, 2010, Varna, Bulgaria). The full collection is available at http://www.emis.de/journals/SIGMA/SIDE-9.html Notation 1.1. Here and hereafter N is an arbitrary positive integer (generally N ≥ 2), superimposed dots denote differentiations with respect to the independent variable t ("continuous time"), and the N dependent variables z n ≡ z n (t) may be interpreted as the coordinates of N point-like unit-mass particles -henceż n denote their velocities andz n their accelerations, consistently with the interpretation of (1.1a) as a set of Newtonian equations of motion with velocitydependent forces. The indices -such as n, m, j, k -generally run from 1 to N (below, as a convenient reminder, we often indicate this explicitly; as well as the exceptions to this rule).
Hereafter we denote as "dynamical system of goldfish type" any dynamical system characterized by Newtonian equations of motion featuring in their right-hand sides -which have, in the Newtonian context, the significance of "forces" -a velocity-dependent term such as that appearing in the right-hand side of (1.1a) (of course in addition to other terms). Let us also emphasize that the dynamical system (1.1) is the simplest model belonging to the Ruijsenaars-Schneider integrable class [3,4].
For instance a simple extension of the above model (reducing to it for ω = 0) is characterized by the Newtonian equations of motion: z n − z m , n = 1, . . . , N. (1.2a) The corresponding solution of the initial-value problem is again given by a simple rule: the N values of the dependent variables z n ≡ z n (t) at time t are related by the formula z n (t) = ζ n (t) exp(−iαωt), n = 1, . . . , N, (1.2b) to the N solutions ζ n (t) of the algebraic equation (for the unknown ζ) [ζ − z j (0)]); or, equivalently but more directly, the N values of the dependent variables z n ≡ z n (t) at time t are the N solutions z n (t) of the algebraic equation (for the unknown z) N k=1ż (1.2d) Notation 1.2. Here and hereafter i is the imaginary unit, i 2 = −1, ω is an arbitrary constant (dimensionally, an inverse time), and α is an arbitrary (dimensionless) constant. Clearly -unless both ω and αω are both imaginary, Re(ω) = Re(αω) = 0 -the time-evolution of this system takes place in the complex z-plane, i.e. the dependent variables z n ≡ z n (t) are complex ; but it may as well be viewed as describing the evolution of N point-like particles moving in the real xy-plane -whose positions at time t are characterized by the (real) Cartesian coordinates x n ≡ x n (t), y n ≡ y n (t) -by setting z n (t) = x n (t) + iy n (t); and one of the remarkable features of the resulting real model is the possibility to write its Newtonian equations of motion in covariant -i.e., rotation-invariant -form, see Chapter 4 of [4]. Hereafter we generally refer to this model, and its generalizations, see below, in their complex versions.
Remark 1.1. Clearly for ω = 0 this model, (1.2), reduces to the previous model (1.1). For ω = 0 the time evolution of this model, (1.2), depends mainly on the values of the two constants ω and αω, as displayed by its solution, see (1.2d). If both these constants are real, Im(ω) = Im(αω) = 0 (hence as well Im(α) = 0), the time evolution of this model is confined, indeed completely periodic if the real number α is rational, while if α is irrational it is multiply periodic, being a nonlinear superposition of two periodic evolutions with the two noncongruent periods T = 2π/ |ω| and T /α. Note that these outcomes obtain for generic initial data: hence, for α rational, α = q/p with q and p coprime integers (and p > 0), this system is isochronous, its generic solutions being completely periodic with period pT -or possibly with a period which is a, generally small, integer multiple of pT : indeed, when the equation (1.2d) is itself periodic with period pT , the unordered set of its N roots is clearly periodic with the same period pT , but the periodicity of the time-evolution of each individual coordinate z n (t) may then be a, generally small, integer multiple of pT due to the possibility that different roots get exchanged through the time evolution (for a discussion of this phenomenology -including a justification of the assertion that the relevant integer multiple of pT is generally small -see [5]).
On the other hand, if ω is real but α is imaginary, say αω = iγ with γ real and nonvanishing, then clearly in the remote future -i.e., as t → ∞, and up to relative corrections of order exp(−|γ|t) -all the N coordinates z n (t) tend to the origin, z n (∞) = 0, if γ < 0, while if γ > 0 they all diverge (see (1.2b) and (1.2c)).
If instead ω is not real, Im(ω) = 0, then in the remote future (i.e., as t → ∞, and up to relative corrections of order exp(−| Im(ω)|t)) the N solutions of (1.2c) become asymptotically, if Im(ω) > 0, the N solutions ζ n = ζ n (∞) of the time-independent polynomial equation of order N in ζ Note that this implies (see (1.2b)) that, if Im(ω) > 0 but αω is real, αω = ρ with ρ real and nonvanishing, then the model (1.2) is asymptotically isochronous, its generic solutions becoming, in the remote future, completely periodic with period 2π/|ρ|, up to corrections vanishing exponentially as t → ∞ (for a more detailed discussion of the notion of asymptotic isochrony see Chapter 6, entitled "Asymptotically isochronous systems", of [6]).
For ω = 0 (i.e., when the model (1.2) reduces to (1.1)) it is possible to restrict consideration to real dependent variables z n , but even then it is more interesting not to do so, so that the time evolution takes place in the plane rather than on the real line: see the remarkable behavior of this dynamical system in this case ("the game of musical chairs"), as detailed in Section 4.2.4 of [4]. Hence let us reiterate that we always consider the dependent variables z n to be complex numbers, both in the continuous-time case, z n ≡ z n (t), 0 ≤ t < ∞, and (see below) in the discrete-time case, z n ≡ z n (ℓ), ℓ = 0, 1, 2, . . . .
Another large class of solvable dynamical systems "of goldfish type" is characterized by the equations of motion z n = a 1żn + a 2 + a 3 z n − 2(N − 1)a 4 z 2 n + N m=1, m =n (z n − z m ) −1 2ż nżm + (a 5 + a 6 z n )(ż n +ż m ) + a 7 z n (ż n z m +ż m z n ) + 2 a 8 + a 9 z n + a 10 z 2 n + a 4 z 3 n , n = 1, . . . , N, featuring 10 arbitrary constants (see [4, equation (2.3.3-2)]). In this case the solvability is achieved by identifying the N dependent variables z n (t) with the N roots of a time-dependent polynomial ψ(z, t) of degree N in z satisfying a linear second-order PDE in the two independent variables z and t.
For an explanation of the origin of the name "goldfish" attributed to these models see Section 1.N of [6] and the literature cited there. In this book [6] (see in particular its Section 4.2.2, entitled "Goldfishing", and the papers referred to there) several other solvable models "of goldfish type" are reported, including isochronous ones (i.e., models featuring solutions which are completely periodic with a period independent of the initial data). A few additional models of goldfish type have been identified more recently [7,8,9].
The most remarkable aspect of these dynamical systems is their solvability, namely the possibility to solve their initial-value problems by algebraic operations, amounting generally to finding the N eigenvalues of an N × N explicitly known time-dependent matrix (see below), or equivalently to finding the N roots of an explicitly known time-dependent polynomial of degree N (see for instance (1.1b) and (1.2d)). Quite interesting is also the identification of multiply periodic, completely periodic, or even isochronous or asymptotically isochronous cases.
In the present paper we present various discrete-time dynamical systems "of goldfish type", so denoted because all these models reduce, in the limit when the discrete-time independent variable ℓ = 0, 1, 2, . . . becomes continuous, to continuous-time dynamical systems of goldfish type. All these models are moreover solvable, i.e. the solution of their initial value problems can be achieved by finding the N eigenvalues z n (ℓ) of N × N matrices explicitly known in terms of the initial data and of the discrete-time independent variable ℓ; or equivalently by finding the N roots z n (ℓ) of a polynomial, of degree N in the complex variable z, as well explicitly known in terms of the initial data and of the discrete-time independent variable ℓ. Some of these models feature interesting behaviors, even isochrony or asymptotic isochrony. Two of these models (see Subsection 2.1 and 2.2) were treated in the paper [10], which has not been published because -after it was submitted for publication but before getting any feedbacknew solvable models were identified and it was therefore considered preferable to report all these models in a single paper, this one. The main properties of each of these discrete-time models are reported in Section 2, and proven in Section 3. These properties include the display of the equations of motion of these discrete-time models, the solution of their initial-value problems, a terse discussion (for the first three models) of their behavior including the possibility that for special values of some of their parameters they possesses periodic or multiperiodic solutions or even display isochrony or asymptotic isochrony, and some mention of their continuous-time limits. Section 4 entitled "Outlook" concludes the paper: in it a general framework is outlined which might allow the identification of additional solvable discrete-time models. And some mathematical developments are confined to two appendices.
These findings are congruent with the recent surge of interest for discrete-time evolutionsin particular, such evolutions which are in some sense integrable or even solvable. Given the large body of research devoted to these topics over the last two decades, our reference to the relevant literature shall be limited to citing the following surveys: [11,12,13,14,15]. But a special mention must be made of the seminal papers by Nijhoff, Ragnisco, Kuznetsov and Pang, see [16] as well as the earlier paper [17], where discrete-time versions were introduced of the well-known integrable "Ruijsenaars-Schneider" and "Calogero-Moser" dynamical systems; as well as of the paper by Suris [18], which treats specifically a discrete-time version of the original goldfish model. Some results of these papers refer to models whose equations of motion feature trigonometric/hyperbolic or even elliptic functions, and are therefore more general than those treated in the present paper, whose equations of motion only feature rational functions (see below); on the other hand the findings reported below include more general models than those previously treated, demonstrate the solvability of these models by an approach somewhat different from those previously employed, and, most significantly, display the possible emergence of remarkable phenomenologies -including periodicity and even isochrony or asymptotic isochrony, see below -not previously identified for this kind of discrete-time dynamical systems.
Let us also mention that the approach developed below also allows to identify and investigate discrete-time variants of another class of solvable continuous-time dynamical systems, the prototype of which is characterized by the Newtonian equations of motion (with c an arbitrary constant), instead of (1.1a). But in this paper we merely indicate, at the appropriate point, how to proceed in this direction, postponing a complete treatment of this development to a separate paper. And let us finally pay tribute to Olshanetsky and Perelomov who were the first to show, more than 35 years ago, that the time-evolution of a nontrivial many-body system could be usefully identified with the evolution of the eigenvalues of a matrix itself evolving in a much simpler, explicitly solvable, manner: see [19] and their other papers referred to in Section 2.1.3.2 of [4], entitled "The technique of solution of Olshanetsky and Perelomov". The present paper extends their approach to the discrete-time context.

Results
In this section we report the main results of this paper; they are then proven in Section 3.
Notation 2.1. Hereafter the dependent variables are indicated again as z n , but they are now functions, z n ≡ z n (ℓ), of the discrete-time variable ℓ taking the integer values ℓ = 0, 1, 2, . . . ; and superimposed tildes indicate generally a unit increase of the independent variable ℓ, for instancẽ z n ≡ z n (ℓ + 1), z n ≡ z n (ℓ + 2). Hereafter δ nm is the standard Kronecker symbol, δ nm = 1 if n = m, δ nm = 0 if n = m, and underlined quantities are N -vectors, for instance z ≡ (z 1 , . . . , z N ). For the remaining notation we refer to Notation 1.1, see above, and to specific indications given case-by-case below.
As reported in this section and explained in Sections 3 and 4, the solvable models considered in this paper generally feature three equivalent versions of the second-order "equations of motion" characterizing their evolution in discrete-time. The treatment of the first model given in Subsections 2.1 and 3.1 is somewhat more detailed than that provided for the other models in the subsequent subsections where, to avoid repetitions, we often refer to the treatments provided in Subsections 2.1 and 3.1. And already in this section, as well as in Section 3, we often take advantage -to simplify the presentation of some results -of identities and lemmata collected in Appendix A.

First model
The first model is defined by the following second-order discrete-time equations of motion: the N values of the twice-updated variables z n ≡ z n (ℓ + 2) are given, in terms of the 2N values of the variables z m ≡ z m (ℓ),z m ≡ z m (ℓ + 1), by the N roots of the following (single) algebraic equation in the unknown z, which clearly amounts to a polynomial equation of degree N in this variable z (as it is immediately seen by multiplying this equation by the polynomial N m=1 (z − az m )). Here and below a is an arbitrary (dimensionless, nonvanishing) constant. A neater version of this formula is easily obtained by multiplying it by a and by then using the identity (A.10) with η n = az n , ζ n = a 2 z n , n = 1, . . . , N . It reads An equivalent formulation of this model is provided by the following system of N polynomial equations of degree N for the twice-updated coordinates z n ≡ z n (ℓ + 2): Again, a neater version of this formula is easily obtained by dividing it by a and by then using the identity (A.10), now with z = a 2 z n , η n = az n , ζ n = z n , n = 1, . . . , N . It reads N j=1 z j − a 2 z ñ z j − az n = (1 + a)a N −1 , n = 1, . . . , N.
And a third, equivalent version of this model is provided by the following system of N polynomial equations of degree N for the twice-updated coordinates z n ≡ z n (ℓ + 2): 3) The similarities and differences among these three sets of "equations of motion", (2.1), (2.2) and (2.3), are remarkable: let us reemphasize that they in fact yield the same evolution in discrete-time of the N coordinates z n ≡ z n (ℓ). Particularly remarkable is their similarity in the special a = 1 case, when the 3 versions (2.1b), (2.2b) and (2.3) of the equations of motion read as follows: The last of these three systems coincides with equation (1.8) of [18]. 2)is invariant under an arbitrary rescaling of the dependent variables, z n ⇒ cz n with c an arbitrary constant; including the special case c = exp(iγ) with γ an arbitrary real constant, corresponding to an overall rotation around the origin in the complex z-plane.
The solution of the initial-value problem for this model is given by the following where of course v m (0) indicates the value of v m (z,z) corresponding to the initial data z = z(0), z =z(0) ≡ z(1). A neater, equivalent formulation of this finding -obtained from (2.4) via Lemma A.4 with ζ n = z n (0)a ℓ and η m = v m (0)(a ℓ − 1)/(a − 1) -states that the N coordinates z n (ℓ) are the N solutions of the following algebraic equation in z: And another, even neater, equivalent formulation -obtained from this via the identity (A.10) with z replaced by za 1−ℓ , η k = az k (0), ζ j =z j (0) = z j (1) -states that the coordinates z n (ℓ) are the N solutions of the following algebraic equation in z: The last two equations become of course polynomial equations of degree N in z after multiplica- These formulas are also valid for a = 1 (by taking the obvious limit, i.e. replacing (a p − 1)/(a − 1) with p). If instead |a| < 1, then clearly for all (positive) values of ℓ the matrix U (ℓ) is bounded and U nm (∞) = v m (0)/(1 − a); hence for all values of ℓ the N coordinates z n (ℓ) are bounded and N − 1 of them vanish as ℓ → ∞ while one of them tends to the value (2.4b) and the identity (A.11) (with η k = a z k (0), ζ j =z j (0) = z j (1)). If a = 1 but it has unit modulus, a = exp(2πiλ) (2.6) with λ real and not integer, then clearly the matrix U (ℓ) is again, for all values of ℓ, bounded; and if moreover λ is a (strictly, i.e. non integer) rational number, with K and L two coprime integers and L > 1, then clearly the matrix U (ℓ) is periodic with period L, hence the (unordered) set of its N eigenvalues z n (ℓ) is as well periodic with period L. This shows that in this case, see (2.7), the discrete-time goldfish model, see (2.1) or (2.2) or (2.3), is isochronous. On the other hand if λ is real and irrational, then clearly the time evolution of this discrete-time dynamical system is not periodic: indeed, while the right-hand side of (2.4a) (with (2.6) and λ real and irrational) is periodic (with unit period) as a function of the real variable τ = λℓ, clearly it is not periodic as a function of the variable ℓ taking the integer values ℓ = 0, 1, 2, . . . . (We made this analysis, for convenience, referring to the coordinates z n (ℓ) as the eigenvalues of U (ℓ), see (2.4); of course an analogous discussion could be made on the basis of the alternative identification of the coordinates z n (ℓ) as the N roots of the polynomial equation (2.5) -whose similarity with (1.2d) is in any case to be noted, see below.) To explore the transition from the discrete-time independent variable ℓ to the continuous-time variable t one makes the formal replacements and (with a slight abuse of notation) with ε infinitesimal. It is then a matter of standard, if a bit cumbersome, algebra, to verify that the insertion of this ansatz, see (2.8b) and (2.8c), in (2.1b) or (2.2b) or (2.3) yields a trivial identity to order ε 0 = 1, while to order ε it reproduces (1.2a) with α = 1, reading Likewise, the discrete-time solution formula (2.5) becomes, in the continuous-time limit, which coincides with (1.2d) with α = 1. A terse outline of the derivation of these results is provided at the end of Subsection 3.1. To higher order in ε one would obtain additional relations satisfied by the solution z n (t) of this continuous-time goldfish model, which might alternatively be obtained by differentiating its equations of motion (2.9a).
Remark 2.2. Clearly, for ω real and nonvanishing, this continuous-time model, (2.9a), is isochronous: see (2.9b) and/or Remark 1.1. This is consistent with the fact that the limiting replacement (2.8b) can be considered to obtain from (2.7) -entailing isochrony of the discretetime model -by identifying εω with 2πK/L in the context of the replacement (see (2.8)) of the unit interval in the discrete-time model with the infinitesimal time interval ε to make the transition to the continuous-time case.
Remark 2.3. At every step of the discrete-time evolution the N values of the twice-updated variables z n ≡ z n (ℓ + 2) are given, in terms of the N unupdated variables z m ≡ z m (ℓ) and the N once-updated variablesz m ≡ z m (ℓ + 1), as the N roots of a polynomial, of degree N in its argument z, whose coefficients are explicitly defined in terms of the 2N unupdated and once-updated variables: see ( , provided at every step of the discrete-time evolution the appropriate identification is made of the value of each twice-updated coordinate z n ≡ z n (ℓ + 2) (among the unordered set of N values yielded by the discrete-time equations of motion) by an argument of contiguity withz n ≡ z n (ℓ + 1) and z n ≡ z n (ℓ); and likewise an appropriate identification is made by contiguity of each coordinate z n (ℓ + 1) with the corresponding coordinate z n (ℓ) (among the unordered set of N values yielded by Proposition 2.1) -these arguments of contiguity taking the place of the continuity of z n (t) as function of t applicable in the continuous-time case. But the contiguity argument breaks down if the positions at time ℓ of two different points, z n (ℓ) and z m (ℓ) with n = m, get too close to each other, corresponding to a quasi-collision, or even coincide, corresponding to an actual collision; which is however not featured by the generic solution of the discrete-time model (2.1) -nor of the standard goldfish models (1.2a) or (1.3) -clearly emerging only for a set of initial conditions z n (0),z(0) ≡ z n (1) having unit codimension in the 2N -dimensional (complex) phase space z,z.
This important remark is applicable to all the discrete-time models considered below, although it will not be repeated.
Remark 2.4. Several of the formulas written above (in this section) simplify somewhat via the following replacement of the dependent variables: In particular the 3 equivalent versions (2.1b), (2.2b) and (2.3) of the discrete time equations of motion are thereby reformulated to read and correspondingly the formula (2.5) providing the solution of the initial-value problem reads Somewhat analogous remarks are applicable to all the discrete-time models considered below; their explicit implementation is left to the interested reader.

Second model
In this subsection we treat rather tersely a discrete-time dynamical system that generalizes the discrete-time goldfish model described in the preceding Subsection 2.1. This generalization amounts to the presence of an additional free parameter, b: indeed, for b = 0 one reobtains the model treated in the preceding Subsection 2.1 (hence in this subsection we assume that b does not vanish, b = 0).
The three equivalent versions of the equations of motion of this model read as follows. The first version identifies the twice updated coordinates z n (ℓ + 2) as the N solution of the following equation in z (amounting to the identification of the N roots of a polynomial of degree N in this variable): The second and third versions consist of the following two systems: The solution of this model is provided by an analog of (the first part of) Proposition 2.1, reading as follows: The N values z n (ℓ) of the dependent variables at the discrete time ℓ are the N eigenvalues of the N × N matrix where again (see (2.4)) is now defined componentwise as follows: with the quantities v m (0) defined again as in Subsection 2.1 (see (2.4b) and the sentence following this formula). Let us recall that I is the N × N unit matrix (whose presence in (2.11a), however, might well be considered pleonastic).
As evidenced by a comparison of (2.11a) with (2.4a), the behavior of the solutions of this model (with b = 0) are less simple than those of the model discussed in the preceding Subsection 2.1. In particular a confined behavior emerges only, see (2.11a), from initial data z n (0), z n (0) ≡ z n (1) implying, via (2.11c) and (2.4b), that all the N eigenvalues of the N × N matrix aI + bV (0) have modulus not larger than unity, reading exp(−q n + 2πir n ) with the numbers q n and r n real and the N numbers q n nonnegative, q n ≥ 0. If moreover the N numbers q n all vanish and the N numbers r n are all rational, the behavior is periodic (but not isochronous, since these numbers, q n and r n , generally depend on the initial data; see (2.11c) and (2.4b)). While, if some of (but not all) the N numbers q n are positive, and none is negative, then the phenomenology we just described (corresponding to the q n 's all vanishing) emerges only asymptotically, as ℓ → ∞, up to corrections of order exp(−qℓ) with q the smallest of the nonvanishing numbers q n -provided all those r n 's are rational whose corresponding q n vanish.
Let us finally mention that, also for this second model, a transition from the discrete-time independent variable ℓ to the continuous-time variable t can be performed (as tersely outlined at the end of Subsection 3.2); but the continuous-time goldfish-type model obtained in this manner turned out to be, to the best of our knowledge, new, hence it seemed appropriate to devote a separate paper to it, see [8].

Third model
The third model is another one-parameter extension of the model treated in Subsection 2.1 (different from that treated in the preceding Subsection 2.2). Again its discrete-time equations of motion can be presented in three equivalent versions.
The first is characterized by this prescription: the twice-updated N coordinates z n ≡ z n (ℓ+2) are the N roots of the following equation in the variable z, amounting again to the determination of the N roots of a polynomial of degree N in the variable z. Here and below a + and a − are 2 arbitrary constants.
Remark 2.6. As entailed by a comparison of these discrete-time second-order equations of motion with those of the first model, see (2.1a), this third model coincides, for a − = 1, with the first model with a = a + .
A neater formulation of these equations of motion reads (after multiplication by a + , via (A.10) with ζ k = a 2 + z k , η k = a +zk ) as follows: An equivalent, second formulation of this model is provided by the following system of N polynomial equations for the twice-updated coordinates z n ≡ z n (ℓ + 2): and a neater version of these equations of motion reads (again via (A.10), but now with ζ k = z k , η k = a +zk and z replaced by a 2 And a third, also equivalent, version of these equations of motion reads as follows: (2.14) Remark 2.7. Remark 2.1 also holds for this model.
The solution of the initial-value problem for this discrete-time dynamical system is provided by the following where the two constant (i.e., ℓ-independent) N × N matrices C + and C − are defined in terms of the 2N initial data z n (0) andz(0) ≡ z n (1) by the formula with the two matrices U (0) and U (1) defined componentwise as follows: Note that we are, for simplicity, assuming that the two coupling constants a ± are different, a + = a − (see (2.15b)).
It is plain from these formulas that, if the two "coupling constants" a ± (are different and) are conveniently written as follows, with q ± and r ± real, then, if the two numbers q ± are both nonnegative, q ± ≥ 0, the time evolution of the N × N matrix U (ℓ) is bounded for all values of the discrete-time independent variable ℓ = 0, 1, 2, . . . , hence its N eigenvalues z n (ℓ) are all as well bounded (the motion is confined); if in particular the two numbers q ± both vanish, q ± = 0, and the two numbers r ± are both rational numbers, r ± = K ± /L ± with K + , L + and K − , L − coprime integers (and, for definiteness, L ± > 0), then the discrete-time evolution of the matrix U (ℓ) is periodic (with a period L independent of the initial data, being the minimum common multiple of L + and L − , hence the discrete-time dynamical system (2.12) is isochronous; while if, of the two numbers q ± , one vanishes and the other is positive, q + = 0, q − = q > 0 respectively q − = 0, q + = q > 0, and r + respectively r − are rational numbers, then the isochronous behavior (with period L + respectively L − ) only emerges asymptotically, as ℓ → ∞, up to corrections of order exp(−qℓ). While clearly if q + and q − are both positive entailing (see (2.16)) |a ± | < 0, then (see (2.15a)) the matrix U (ℓ), hence as well all it eigenvalues z n (t), vanish asymptotically (as t → ∞): z n (∞) = 0, n = 1, . . . , N . Finally let us mention the transition from this discrete-time model to its continuous-time counterpart. The treatment is completely analogous to that detailed at the end of Subsection 2.1; except that now (2.8b) must be replaced by It is then easily seen that again, at order ε 0 = 1, one gets from (2.12b) or (2.13b) or (2.14) a trivial identity, while at order ε one gets the continuous-time goldfish equations of motion (1.2a); and the solution of this model, see Proposition 2.3, reproduces in this continuous-time limit the prescription (1.2d).

Fourth model
The fourth model is also characterized by three equivalent versions of its discrete-time equations of motion. The first consists of the following prescription: the twice-updated N coordinates z n ≡ z n (ℓ + 2) are the N roots of the following equation in the variable z, where the N quantitiesĝ k (z,z) are defined as follows: Throughout this subsection, the following Subsection 3.4, and Appendix B, where α, β, γ, η and ρ are 5 arbitrary constants (but the 3 constants β, η, ρ only enter as βρ and ηρ, hence any one of these three constants could be replaced by unity without significant loss of generality); in the following we use interchangeably these constants in order to simplify some formulas. An equivalent formulation of these discrete-time equations of motion reads as follows: It is a matter of trivial algebra to rewrite these equations of motion, (2.18), as follows: And, as shown at the end of Appendix B, a neater version of this system of equations of motion then reads as follows: η z j + αβ ηz k + β = γa N −1 ηaz n + β + ηb ηz n + β , n = 1, . . . , N.
And a third, equivalent formulation of these equations of motion reads as follows: withǧ n z, z respectivelyĝ n (z,z) defined by (2.18b) respectively (2.17b).
Remark 2.8. Above and below we assume for simplicity that the parameters characterizing this model have generic values, for instance γ = 0 and γ = 1 (see (2.17a) and (2.17c)) and a = 0 (see (2.18b)). The solution of the initial-value problem for this discrete-time dynamical system is provided by the following Proposition 2.4. The N coordinates z n (ℓ) are the N eigenvalues of the N × N matrix where the N × N matrix P (ℓ 1 , ℓ 2 ) is defined as follows

21b)
and the two ℓ-independent N × N matrices A and B are defined as follows

21c)
Here and throughout we use the convention that (for arbitrary finite X j ) As for the two N × N matrices U (0) and V (0), they are defined in terms of the 2N initial data z n (0) andz n (0) ≡ z n (1) as follows: hence an explicit expression of the N × N matrix V (0) reads, componentwise, as follows: Remark 2.10. It is relevant to this expression, (2.21), of the N × N matrix U (ℓ) -whose N eigenvalues provide the N coordinates z n (ℓ) -that (2.21b) and (2.21c) entail (a n γ j + a), a n = ηv n − η β b, b n = βv n − b, n = 1, . . . , N, with v n the N (ℓ-independent) eigenvalues of the N × N matrix V (0) and Q the corresponding (ℓ-independent) diagonalizing matrix, And let us mention that, also for this fourth model, a transition from the discrete-time independent variable ℓ to the continuous-time variable t can be performed (see the end of Subsection 3.4). And, as in the case of the second model, also in this case the continuous-time goldfish model thereby obtained turned out to be, to the best of our knowledge, new. Hence it seemed appropriate to devote to this model a separate paper [9].

Proofs
In this section we prove the findings reported in the preceding Section 2.
The basic strategy to obtain all these results goes as follows. The starting point is a solvable system of two matrix first-order discrete-time ODEs, saỹ where ℓ = 0, 1, 2, . . . is the discrete-time independent variable, the two dependent variables U ≡ U (ℓ), V ≡ V (ℓ) are N × N matrices and of course superimposed tildes denote the updating of the discrete-time,Ũ ≡ U (ℓ + 1),Ṽ ≡ V (ℓ + 1). The solvable character of this matrix system entails the possibility to obtain explicitly the solution of its initial-value problem. Four cases when this is possible -corresponding to 4 simple assignments of the functions F 1 (U, V ) and F 2 (U, V ) -are treated in the following 4 subsections. Note that the two functions F 1 (U, V ), F 2 (U, V ) are assumed to depend on no other matrix besides U and V (and the unit matrix I); they may of course feature some scalar constants, and the order in which the two, generally noncommuting, matrices U and V appear in their definition is of course relevant: see below.
One assumes moreover that the N × N matrix U ≡ U (ℓ) is diagonalizable and denotes as R ≡ R(ℓ) the diagonalizing N × N matrix: where the notation z n (ℓ) for the N eigenvalues of the N ×N matrix U ≡ U (ℓ) shall be justified by the identification, see below, of these quantities with the dependent variables of the discrete-time dynamical systems introduced above.
Next we introduce the two matrices M (ℓ) and Y (ℓ) defined as follows: One then, by inserting (3.2a) and (3.3c) in (3.1), obtains the following system of two firstorder discrete-time N × N matrix evolution equations: and from these two matrix equations, by making a convenient ansatz for the two matrices M ≡ M (ℓ) and Y ≡ Y (ℓ) in terms of the 2N quantities z n ≡ z n (ℓ) andz n ≡ z n (ℓ + 1) -an ansatz which must of course be consistent with these two matrix evolution equations -one obtains a system of N second-order discrete-time evolution equations for the N coordinates z n ≡ z n (ℓ). This last step is of course only possible for special assignments, in the discrete-time matrix evolution equations (3.1), of the two matrix functions F 1 (U, V ) and F 2 (U, V ), see below.
The discrete-time dynamical system thereby obtained is then solvable, since the quantities z n ≡ z n (ℓ) are the N eigenvalues of the N × N matrix U ≡ U (ℓ) which, as solution of the, assumedly solvable, matrix evolution system (3.1), can be explicitly evaluated. How this works out is shown in detail in the following subsections: in more detail in Subsection 3.1, where the simplest case is treated.
Let us also mention, once and for all, that in the following we will conveniently assume that the matrix U is initially diagonal: implying (up to the ambiguity mentioned above, see Remark 3.1) Here and throughout I is the N × N unit matrix, i.e., componentwise, I nm = δ nm .

Solution of the first model
The point of departure to obtain the findings reported in Subsection 2.1 is the following discretetime first-order, linear, matrix system (see (3.1)): where a is an arbitrary scalar constant. Note that the second of these two ODEs entails that in this case V is a constant (i.e., ℓ-independent) N × N matrix, V (ℓ) = V (0). It is plain that the solution of the corresponding initial-value problem for the N × N matrix U reads Let us now proceed as indicated in the first part of Section 3. It is then easily seen (via (3.2) and (3.3)) that the first of the two discrete-time matrix evolution equations (3.6a) yields (see (3.4)) the matrix equation This derivation shows that this system of N 2 discrete-time equations of motion is equivalent to the solvable equation of motion (3.6a) for the N ×N matrix U ; hence it is just as solvable. Note that the dependent variables are now the N coordinates z n and the N 2 matrix elements Y nm (of which only N (N − 1) are significant, see Remark 3.2; so the number of equations and the number of dependent variables tally). To obtain a model that qualifies as discrete-time analog of the continuous-time goldfish model (2.9a) we need to distill from this system a set of only N equations of motion involving only the N coordinates z n . The standard trick to do so (see, for instance, Section 4.2.2 entitled "Goldfishing" of [6]) is to identify -if possible -an ansatz which expresses the N 2 components of the matrix Y in terms of the 2N quantities z n ,z n , yielding N equations of motion involving only the N coordinates z n ,z n and z n -to be interpreted as equations of motion of the discrete-time goldfish -and implying that the N 2 equations of motion (3.8) are all satisfied, thanks to these very equations of motion.
An educated guess for such an ansatz reads as follows: Y nm = g m , n, m = 1, . . . , N. (3.9) Note that we reserve at this stage the option to assign the N quantities g m . Via this ansatz the equations (3.8) become hence they amount to the following 2 systems, each involving only N equations: The unit in the right-hand sides could of course be replaced by an arbitrary constant c -of course the same constant in (3.10b) and (3.10c) -but this would merely entail an irrelevant rescaling of g k by c; see below.
These are now two sets of N equations, each featuring linearly the N quantities g k , that we like to eliminate in order to obtain a set of N "equations of motion" determining the twice updated coordinates z n ≡ z n (ℓ+2) in terms of the 2N coordinates z n ≡ z n (ℓ) andz n ≡ z n (ℓ+1). There are three alternative strategies to achieve this goal. One can solve the first linear system thereby obtaining g k as a function of z andz, and then insert this expression g k ≡ĝ k (z,z) in the second system; alternatively, one can solve the second linear system, thereby obtaining g k as a function ofz and z, and then insert this expression g k ≡ǧ k z, z in the first system; or one can equate the two expressions of g k obtained solving the first, respectively the second, system, i.e. writeĝ k (z,z) =ǧ k z, z . Clearly the three sets of equations of motion obtained in this manner are equivalent, i.e. they characterize the same discrete-time evolution of the N coordinates z n ≡ z n (ℓ); but they may seem quite different (indeed, see (2.1), (2.2) and (2.3)). Note that we introduced a superimposed decoration on the functionsĝ k (z,z) respectivelyǧ k z, z to emphasize that the functional dependence on their arguments is generally different, as implied by their definitions as solutions of (3.10b) respectively of (3.10c).
Let us first of all see what the first approach yields. From (3.10b) one obtains (via Lemma A.1 reported in Appendix A, with ξ k =z k , η n = az n , c = 1) the following expression of g k ≡ĝ k (z,z): The insertion of this expression of g k in (3.10c) yields the equations of motions (2.1a). Likewise, the second approach yields, from (3.10c) (again via Lemma A.1, but now with ξ k = az k , η n = z n , c = −1) the following expression of g k ≡ǧ k z, z : The second version, (2.2a), of the discrete-time equations of motion follows by inserting this expression of g k in (3.10b).
And the third approach yields, by equating (3.11) to (3.12), the third version, (2.3), of the discrete-time equations of motion.
We have seen that the solutions z n (ℓ) of these discrete-time equations of motion are provided by the eigenvalues of the N × N matrix U (ℓ), see (3.6b). To prove Proposition 2.1 we must now obtain from (3.6b) (also taking advantage of the ansatz (3.9)) the expression (2.4) of this matrix in terms of the initial data z n (0),z n (0) ≡ z n (1) of the discrete-time dynamical system. This requires that we express the two matrices U (0) and V (0) appearing in the right-hand side of (3.6b) in terms of the initial data z n (0),z n (0) ≡ z n (1).
Let us now provide a terse treatment of the transition from the discrete-time equations of motion (2.1b), which we conveniently re-write here as follows, N j=1 z n − a 2 z j z n − az j = 1 + a, n = 1, . . . , N, (3.15) to the continuous-time case, see (2.9a). It is then appropriate to treat separately the factor with j = n in the product appearing in the left-hand side of (3.15), and all the other factors with j = n. The basic equations are (2.8), entailing z n − a 2 z n = 2ε(ż n + iωz n ) + ε 2 2z n + ω 2 z n + O ǫ 3 , Hence, after a little algebra, see (2.8b). It is then clear that the insertion of these two formulas, (3.16c) and (3.16d), in (3.15) yields, at order ε 0 = 1, the trivial identity 2 = 2, and at order ε the equations of motion of the continuous-time goldfish model (2.9a).
In an analogous manner one reobtains (2.9a) from (2.2b) or from (2.3). Let us also show that (2.5), which we rewrite here conveniently as follows, yields, in the continuous-time limit, (2.9b). Indeed the relation a = 1 − iεω (see (2.8b)) entails (via the first of the three relations (2.8a)), and (via the second of the three relations (2.8c), with ℓ = 0). Via these three relations (3.17) becomes i.e. (dividing each numerator by the corresponding denominator) Clearly to order ε 0 = 1 this yields the trivial identity 1 = 1, and to order ε just the formula (2.9b). Let us end this subsection by pointing out that there is another ansatz that allows to transform the system of N 2 equations (3.8) into two separate systems of N equations, but only in the special case a = 1. This alternative ansatz reads (instead of (3.9)) But, as indicated at the end of Section 1, we postpone the treatment of the corresponding class of discrete-time dynamical systems to a separate paper.

Solution of the second model
The proof of the findings reported in Subsection 2.2 is analogous to that provided above, see Subsection 3.1, so our treatment in this subsection is quite terse, being limited to indicate the changes with respect to that reported in the preceding Subsection 3.1. Now the system of matrix evolution equations (3.6a) is generalized to read hence its solution is given by (2.11a). Clearly this evolution equation, (3.18), respectively its solution, (2.11a), reduce to (3.6a) respectively to (3.6b) when b vanishes. The rest of the treatment is analogous. (3.7a) is now generalized to read hence it yields, in place of (3.7b), In place of (3.10a) (again via the ansatz (3.9)) one now has hence in place of (3.10b) and (3.10c) one gets the two sets of N equations By solving the first set one obtains (via Lemma A.2, with ξ k =z k , η n = az n , c n = 1/(1+bz n ), and then the identity (A.7) with z = −a/b, η k = az k , ζ j =z j ) the following expression of g k ≡ĝ k (z,z): By solving instead the second set one obtains (via Lemma A.1, with ξ k = −az k , η n = − z n , c = 1 and g k replaced by g k (1 + bz k )) the following expression of g k ≡ǧ k z, z : The three versions, (2.10), of the equations of motion reported in Subsection 2.2 then follow by inserting (3.20a) in (3.19c), by inserting (3.20b) in (3.19b), and by equating (3.20a) to (3.20b).
Next, let us prove Proposition 2.2. One proceeds again in close analogy to the treatment of the preceding Subsection 3.1, hence we only mention where the treatment here differs from that provided there. It is easily seen that the expression of U (0) is the same as that given there, see (3.13), while the expression of V (0) (because one must now use Lemma A.5 with f n = 1 + bz n (0) rather than f n = 1) is now given by (2.11c). The insertion of these expressions of U (0) and V (0) in (3.18) reproduce (2.11), thereby proving Proposition 2.2.
Let us end this subsection by outlining what happens in the continuous-time limit which obtains by setting a = 1 + εbη, V (0) = εB with ε infinitesimal, and correspondingly replacing the discrete-time matrix evolution equation (3.18) with the matrix ODĖ the solution of which reads As already mentioned in Subsection 2.2, the (continuous-time) goldfish-type model obtainable by focussing appropriately on the evolution of the N eigenvalues z n (t) of this N × N matrix U (t) (evolving according to (3.21a)) was, to the best of our knowledge, new, when the solvable matrix evolution equation (3.21a) was identified as continuous-time limit of (3.18); its treatment is provided in [8].

Solution of the third model
The starting point is the following linear system of two discrete-time matrix evolution equations: where the 3 constants a ± , β are a priori arbitrary (β = 0). It is easily seen that the solution of the initial-value problem for U reads as follows (with an analogous formula for V ): and the two constant matrices C ± given by (2.15b).
Remark 3.3. This solution U (ℓ) of the initial-value problem for the discrete-time N × N matrix evolution equation (3.22a) depends only on the 2 constants a ± : see (3.22b) and (2.15b).
Indeed the system of two first-order discrete-time evolution equations (3.22a) is easily seen to correspond to the single second-order evolution equation from which the constant β has disappeared (but note that this second-order matrix ODE obtains only if β = 0; indeed if β = 0, U satisfies a first-order evolution equation, see the first of the two equations (3.22a).
We then proceed as in the first part of Section 3. It is then easily seen that the two matrix evolution equations (3.22a) become and using this formula it is easily seen that the second can be written, componentwise, as follows: At this point we use again the ansatz (3.9) for the matrix Y nm , the consistency of which is vindicated by the subsequent developments. Here the N quantities g m are again a priori arbitrary; they shall be determined as functions of the 2N un-updated and once-updated coordinates z n ≡ z n (ℓ) andz n ≡ z n (ℓ + 1), or alternatively of the once and twice updated coordinatesz n ≡ z n (ℓ + 1) and z n ≡ z n (ℓ + 2), see below. It is indeed immediately seen that via the ansatz (3.9) the system of N 2 equations (3.24) becomes hence it can be replaced by the following two separated systems of only N equations: The first, (3.25b), of these two systems defines uniquely the N quantities g k ≡ĝ k (z,z), yielding again, via (A.13), the expression (3.11) (with a replaced by a + ): Insertion of this expression in the second, (3.25c), of the two systems written just above then yields the evolution equation (2.12a). The identification of the third discrete-time dynamical system of goldfish type, see (2.12a), is thereby accomplished. The second version, (2.13a), of this model obtains by solving for g k ≡ǧ k z, z the second, (3.25c), of the two systems written above (via Lemma A.1, with ξ k = a +zk , η n = z n and c = −1/a − ), thereby obtaininǧ and by then inserting this expression of g k in (3.25b).
As for the proof of Proposition 2.3, it follows immediately from the above treatment, see in particular (2.15); there only remains to justify the identification of the two matrices U (0) and U (1), see (2.15c) and (2.15d).
The first of the two formulas, (2.15c), is just (3.5a).

Solution of the fourth model
The treatment in this subsection is rather terse, since it is analogous to that of the preceding subsections; and the notation is of course analogous. But now the starting point is the following nonlinear system of two discrete-time matrix evolution equations: featuring the 5 arbitrary constants α, β, η, ρ, γ (which, as noted in Subsection 2.4, can be reduced to 4 by taking advantage of the freedom to rescale V ). It is a standard task to see that the solution of the initial-value problem for this matrix system reads as follows: with U (ℓ) given by (2.21). We now proceed again as in the first part of Section 3. Via (3.3a) and (3.3c) we get from (3.31) the two matrix equations The first, (3.32a), of these two matrix equations entails, componentwise, hence the second, (3.32b), of these two matrix equations yields (when written componentwise) the following N 2 equations: which, as can be easily verified, can be conveniently rewritten as follows: with the two constants a and b defined by (2.17c). Next, we make the following ansatz for the matrix M nm : Here the N quantities g m are a priori arbitrary; they shall be determined as functions of the 2N un-updated and once-updated coordinates z n ≡ z n (ℓ) andz n ≡ z n (ℓ + 1), or of the once and twice updated coordinatesz n ≡ z n (ℓ + 1) and z n ≡ z n (ℓ + 2), see below. It is indeed immediately seen that via this ansatz (3.34) the system of N 2 equations (3.33b) can be rewritten as follows: hence it can be replaced by the following two separated systems of only N equations: The first of these two systems defines uniquely the N quantities g k ≡ĝ k (z,z), yielding their expression (2.17b) (see Appendix B for a proof). The second equation then yields the evolution equation (2.17a). The alternative possibility is to determine (via Lemma A.1, with ξ k = az k , η n = z n − b, c = −1/γ) the quantities g k ≡ǧ k z, z as solutions of the second system, yielding the formula (2.18b); and to then insert these expressions of g k ≡ǧ k z, z in the first system of equations. Clearly in this manner one arrives at the equations of motion (2.19a).
And a third possibility is of course to equate (2.17b) to (2.18b), see (2.20). The identification of the three variants of the equations of motion of the fourth discrete-time dynamical system of goldfish type, see Subsection 2.4, is thereby accomplished.
The proof of Proposition 2.4 follows immediately from the above treatment; and we trust that the identification in terms of the 2N initial data z n (0) and z n (1) of the two matrices U (0) and V (0), see (2.4) and (2.4), is sufficiently obvious (also in the light of the analogous treatment in the preceding subsections of this section) not to require an explicit justification here.
We end this subsection with a terse mention of the continuous-time model that obtains from that treated above (in this subsection) via the limiting transition from discrete to continuous time. The point of departure for the treatment of the continuous-time dynamical system of goldfish type obtained in this manner is the following continuous-time system of two first-order matrix evolution equationṡ U = a 1 U + a 2 V + a 3 U V,V = a 4 + a 5 V, that obtains from (3.31) via the assignments t ⇒ εℓ, U (ℓ) ⇒ U (t), V (ℓ) ⇒ V (t), α = 1 + εa 1 , β = εa 2 , η = εa 3 , ρ = εa 4 , γ = 1 + εa 5 , with ε infinitesimal. The resulting model of goldfish type was, to the best of our knowledge, new ; a detailed treatment of it is provided in [9].

Outlook
In this paper we have introduced and tersely analyzed 4 different discrete-time dynamical systems of goldfish type. The possibility to identify other discrete-time evolution equations amenable to exact treatment by variations of the methodology used in this paper is open: let us outline here an avenue to such generalizations.
Consider the system of two N × N matrix discrete-time first-order evolution equations where the two N × N matrices U ≡ U (ℓ) and V (ℓ) are the dependent variables, ℓ = 0, 1, 2, . . . is the independent discrete-time variable, S 1 , S 2 , S 3 are 3 arbitrary positive integers, F 1,s (u), arbitrary (scalar) functions of their (scalar) argument u (of course becoming N × N matrices when the scalar u is replaced by an N × N matrix). Then introduce the eigenvalues z n (ℓ) of the matrix U (ℓ), as well as the matrices Z ≡ Z(ℓ), R ≡ R(ℓ), M ≡ M (ℓ) and Y ≡ Y (ℓ), as above (see (3.2) and (3.3)). It is then plain that the matrix evolution equation (4.1a) becomes Likewise (4.1b) becomes And via the ansatz Y nm = g m (see (3.9)) this system of N 2 equations can clearly be replaced by the following two systems of N linear algebraic equations for the N quantities g k : One can then solve the first, (4.3a), respectively the second, (4.3b), of these two linear systems for the N quantities g k , getting the expressions g k =ĝ k (z,z) respectively g k =ǧ k z, z . One arrives thereby at three equivalent systems of N discrete-time second-order evolution equations for the N coordinates z n (ℓ): (i) by inserting the expressionĝ k (z,z) in (4.3b); (ii) by inserting the expressionǧ k z, z in (4.3a); (iii) by settingĝ n (z,z) =ǧ n z, z . While of course the evolution in the discrete time ℓ entailed by these equations of motions corresponds to the evolution of the eigenvalues of the matrix U (ℓ) solution of the matrix evolution equation (4.1a) (with the ansatz (3.9) properly taken into account). Hence if that matrix evolution equation, (4.1a), is solvable, the discrete-time dynamical system described by these three equivalent sets of secondorder evolution equations is as well solvable.
One has thereby identified a discrete-time solvable dynamical system. A remaining open question is the extent to which its equations of motion can be exhibited in reasonably neat form: this depends on the extent that the two quantities g k ≡ĝ k (z,z) respectively g k ≡ǧ k z, z defined as solutions of the two linear systems (4.3a) respectively (4.3b) can be expressed more explicitly than via the standard Cramer formula (ratio of two determinants).
Clearly the 4 models treated in this paper belong to this class (4.1) (the last one, however, only if ρ = 0): see (3.6a), (3.18), (3.22a) and (3.31). The interest of additional models treatable via this approach depends on the neatness of the corresponding equations of motion, which can only be investigated on a case-by-case basis.

A Appendix
In this appendix we collect various mathematical developments whose treatment in the body of the paper would interrupt the flow of the presentation.
First of all we report several mathematical identities. We consider all of them well-known, but for completeness we either prove them below, or indicate where proofs can be found. These formulas feature sets of N numbers such as ξ n or η n or ζ n ; these numbers are arbitrary but for simplicity we assume them to be distinct. The formulas of course remain valid when these numbers are not distinct, but possibly only by taking appropriate limits. Sometimes an arbitrary number z also appears.
The identity (A.1) (with z an arbitrary number) is implied by the fact that its left-hand side is a polynomial in z of degree less than N (in fact, of degree at most N − 1) which clearly has the value unity at the N points ζ n , and the right-hand side, i.e. unity, is the unique polynomial of degree less than N in z that has the value unity in N distinct points. The identity (A.2) is the special case of (A.1) with z = 0. The identity (A.3) coincides with equation (2.4.2-32) of [4] (or, as above, it is implied by the observation that its left-hand side is a polynomial in z of degree less than N the values of which at the N points ζ k coincide with the values of the polynomial z n at z = ζ k ). , with x n = ξ n , y n = η n , respectively x n = z, y n = η n , z n = ζ n ). The identities (A.8) respectively (A.9) follow from (A.7) in the limit z → ∞ respectively ζ j → ∞. The identity (A.10) follows from (A.8) and (A.7) via the trivial identity Finally, the identity (A.11) follows from (A.10) in the limit z → ∞, and the identity (A.12) is just the special case of the preceding identity (A.11) with ζ n = 0. Next we report a simple lemma (for a neat proof see for instance [18], or below, after the proof of the following Lemma A.2). To prove this formula one inserts this expression, (A.14b), of g k in (A.14a), and notes that one obtains thereby an equality provided there holds the formula Clearly this formula is implied by the identity (A.6). Lemma A.2 is thus proven. Note that, by setting c n = c in (A.14c) and using the identity (A.8) (with the dummy index k replaced by s and the index n replaced by k), one reobtains (A.13b), thereby proving Lemma A.1.
We now report, and prove, 3 other lemmata. This lemma is an immediate consequence of the preceding Lemma A.3, because the secular equation associated with the matrix (A.16a) (whose roots provide the eigenvalues) is easily seen to coincide (up to an overall, hence irrelevant, multiplicative constant) with the vanishing of the determinant in the left-hand side of (A.15a) with (A.15b) and x k = η k /(z − ζ k ). Then, for fixed m, apply Lemma A.2 with g k replaced by g k M −1 km , and c n replaced by δ nm /f n . This yields, rather immediately, the formula (A.17b), which is thereby proven.

B Appendix
In this appendix we detail the derivation of some findings for the fourth model, firstly the expression (2.17b) ofĝ n (z,z) as solution of the linear system (3.35b), and secondly the derivation of (2.19b) from (2.19a).
The insertion of this expression of σ (0) + a and of the expression (B.3b) of σ (3) in (B.3a) completes the derivation of the expression (2.17b) ofĝ n (z,z).
Finally let us tersely outline the derivation of (2.19b) from (2.19a). Firstly one uses in (2.19a) the identity (B.1); then one uses twice the identity (A.10), with ζ k = z k − b, η k = az k and with z = a(az n + b) respectively with z = −aβ/η. And the rest is trivial algebra.