## Abstract

The general solutions of many three-dimensional Lotka–Volterra systems, previously known to be at least partially integrable, are constructed with the aid of special functions. Examples include certain ABC and May–Leonard systems. The special functions used are elliptic and incomplete beta functions. In some cases, the solution is parametric, with the independent and dependent variables expressed as functions of a ‘new time’ variable. This auxiliary variable satisfies a nonlinear third-order differential equation of a generalized Schwarzian type, and results of Carton-LeBrun on the equations of this type that have the Painlevé property are exploited, so as to produce solutions in closed form. For several especially difficult Lotka–Volterra systems, the solutions are expressed in terms of Painlevé transcendents. An appendix on incomplete beta functions and closed-form expressions for their inverses is included.

## 1. Introduction

### (a) Background

Models of Lotka–Volterra type occur frequently in the physical and engineering sciences, as well as the biological. In any such model, there is a *d*-dimensional state vector *x*=(*x*_{1},…,*x*_{d}), a function of time *τ*, which satisfies an autonomous system of ordinary differential equations (ODEs) of the form
1.1
The overdot signifies d/d*τ*. By letting *x*_{0}≡1 and *a*_{0i}≡0, one can optionally rewrite this as , *i*=0,…,*d*. There is a paucity of closed-form solutions: when *d*≥3 and even when *d*=2, system (1.1) is usually integrated numerically rather than symbolically.

In most applications, *x*_{i}≥0 for all *i*, though it may be useful to allow the *x*_{i} and even *τ* to be complex. Such models were first introduced in population ecology [1,2], where *x* always lies in the non-negative orthant, and the terminology reflects this: *a*_{i0} is the intrinsic growth rate of the *i*th ‘species’, and depending on the signs of the off-diagonal elements of the interaction matrix , one speaks of predation, competition or mutualism. If *a*_{i0}>0 and *a*_{ii}<0, in the absence of the other *d*−1 species the *i*th species will grow logistically, but the behaviour of *x*=*x*(*τ*) in the interior of the orthant may be much more complicated.

Modelling by small-*d* Lotka–Volterra systems has occurred many times in physics. In laser physics, this ranges from the initial study of multimode coupling [3] to a treatment of Raman amplification in an optical fibre [4], in which *τ* is the distance along the fibre. Other physical applications include the modelling of the (integrable) interaction of Langmuir waves in plasmas [5] and more chaotic phenomena: the Kueppers–Lortz instability in a rotating fluid [6,7], for which *d*=3, and the Boltzmann dynamics of a small number of rarefied, spatially homogeneous gases in a background medium [8,9]. Several integrable Lotka–Volterra systems with fairly small *d*, of mathematical interest, have been obtained as spatial discretizations of the Korteweg–de Vries equation or as generalizations of the same [10,11].

In population biology, small-*d* Lotka–Volterra systems include the original and now classical *d*=2 predator–prey model of Volterra and Lotka, well known though structurally unstable; the famous May–Leonard model of *d*=3 cyclically competing species with equal intrinsic growth rates, which for certain parameter values has a stable limit cycle [12]; and generalizations that model the asymmetric competition of *d* species in a chemostat [13,14]. There is also a literature on evolutionarily stable strategies in game dynamics, modelled by differential systems of replicator type that are equivalent to Lotka–Volterra ones [15,16]. Lotka–Volterra systems with *d*≥3 may display chaotic behaviour, a fact first noted in ecological modelling [17,18].

Many non-quadratic nonlinear differential systems can be transformed to (1.1) by changes of variable [19–22]. Complex chemical reactions with mass-action kinetics provide examples. Spatially homogeneous species concentrations evolve in a Lotka–Volterra way only for certain stoichiometries [23,24]. But mass-action systems with non-quadratic polynomial nonlinearities, which can be obtained from (formal) chemical reactions in more than one way, can often be transformed to the Lotka–Volterra form.

Several elementary changes of variable deserve mention, for any *d*-dimensional system of the form (1.1). First, by scaling *τ* one can scale the column vector of intrinsic growth rates and the interaction matrix A. Additionally, by scaling *x*_{1},…,*x*_{d} independently one can scale each of the *d* columns of A. Thus, the parameter space has dimensionality not *d*^{2}+*d* but *d*^{2}−1.

If the intrinsic growth rates *a*_{10},…,*a*_{d0} are *equal* (to *a*_{*0}, say), then in terms of a transformed time *τ*′:=*e*^{a*0τ} the transformed variables *x*_{i}′:=*e*^{−a*0τ}*x*_{i} will satisfy a system of the form (1.1) but with zero growth rates. Furthermore [25], if *x*_{1},…,*x*_{d} satisfy (1.1) with zero growth rates, then the ratios will satisfy a (*d*−1)-dimensional system like (1.1), but with growth rates and interactions . Thus, Lotka–Volterra systems with equal intrinsic growth rates, for which the parameter space dimensionality is effectively *d*^{2}−*d*, should be peculiarly amenable to analysis; though even they have usually been studied not symbolically but numerically.

The object of this paper is the construction of general solutions *x*=*x*(*τ*) of systems of the form (1.1), with the aid of special functions. In §2, some ancillary results on integrable *d*=2 systems are obtained, but in §§3 and 4 the focus is on *d*=3 systems. The special functions employed are the incomplete beta function and its inverse, which sometimes reduces to an elementary or elliptic function. For many *d*=3 systems known to be at least partially integrable, general solutions are constructed with unexpected ease. By default, each of *x*,*τ* is expressed as a function of an auxiliary variable *t*, a ‘new time’. A minor but interesting example is the fully (rather than merely cyclically) symmetric May–Leonard model, which is treated in §3*c*.

In most of the integrated systems, the three species have equal intrinsic growth rates, but not in all. Several systems with unequal ones which have the Painlevé property (PP), meaning that any solution *x*=*x*(*τ*) can be extended analytically to a one-valued function on the complex *τ*-plane, are integrated in §4 with the aid of Painlevé transcendents.

Many of the interesting examples of closed-form solutions *x*=*x*(*τ*) that are presented here are made possible by the results of Carton-LeBrun [26]. She classified all nonlinear third-order ODEs of a certain type that have the PP; and the ODE satisfied by the new time *t*=*t*(*τ*) used here was fortuitously included. She also integrated many such ODEs with the aid of elliptic functions. The examples in §§3 and 4 may serve to reawaken interest in her work.

This paper includes an appendix on closed-form expressions for inverse incomplete beta functions, which should be of independent interest.

### (b) Previous results, reformulated

The following is a summary of symbolic (i.e. non-numerical) results on small-*d* Lotka–Volterra systems, which provides the context for the new results.

Lotka–Volterra solutions cannot always be expressed in closed form, i.e. in terms of elementary or ‘known’ functions, even if the behaviour of the solutions is non-chaotic; or so it is believed. An example is the classical system of Volterra and Lotka, with *d*=2 and no self-interactions (*a*_{11}=*a*_{22}=0). If the intrinsic growth rates are equal, it has a general solution in terms of elementary functions [27]. If the rates are unequal, it is nonetheless integrable (see below) and non-chaotic, but this does not imply the existence of a closed-form general solution. The case of unequal growth rates when *d*=2 has in fact been used as a test-bed for geometric numerical integration schemes [28].

Owing to the difficulty of constructing explicit solutions of systems of the form (1.1), there has been much work on the simpler problem of finding constants of the motion, i.e. first integrals. Any such is a function *I*=*I*(*x*_{1},…,*x*_{d}), preferably one-valued and well-behaved, obeying on any system trajectory. (By the overdot, the total time derivative is meant.) If such an *I* exists, any trajectory is confined to a (*d*−1)-dimensional surface *I*≡*const*., which, for example, facilitates the study of long-time behaviour. If there are *d*−*m* functionally independent first integrals, then motion is confined to an *m*-dimensional surface, obtained as an intersection. When *m*=1, the system is completely integrable (its general solution can be ‘reduced to quadratures’, though not necessarily usefully). When *m*=2, it is partially integrable and is at least non-chaotic.

For *d*=2,3, restrictions on the *d*^{2}+*d* Lotka–Volterra parameters that imply the existence of a single closed-form first integral have been investigated by many authors. In most cases, the resulting first integral is of the form
1.2
with each *f*_{k}=*f*_{k}(*x*_{1},…,*x*_{d}) a polynomial (usually though not always of degree 1 with constant term allowed). Results along this line include those of Cairó and co-workers [29–32], who in many cases used the classical technique of Darboux polynomials (DPs) [33]; and recent results based on an integrating factor technique [34]. In the case when *d*=3 and a single first integral of a form similar to (1.2) is known to exist, there has been some work on the construction of a second, functionally independent first integral of a less algebraic form, involving a quadrature [35–38].

The DP concept is fundamental. A DP *f*=*f*(*x*_{1},…,*x*_{d}) is a polynomial for which along any system trajectory, where the ‘cofactor’ *K*=*K*(*x*_{1},…,*x*_{d}) is necessarily also a polynomial, of degree 1 with constant term allowed. The algebraic surface *f*=0 is thus invariant under the Lotka–Volterra flow. If the cofactors of DPs *f*_{1},…*f*_{r} are linearly dependent, some product of the form (1.2) will be a first integral: generically, a non-trivial one. Finding a first integral of this form, with associated parametric restrictions, is facilitated by the coordinate functions *x*_{1},…,*x*_{d} always being DPs, associated with the invariant hyperplanes *x*_{i}=0, *i*=1,…,*d*, and having respective cofactors . To construct a first integral, a single additional DP *f* will suffice, if one sufficiently restricts parameters to ensure linear dependence of the cofactors of *x*_{1},…,*x*_{d},*f*.

Tables 1 and 2 are a distillation of the results of the preceding authors on degree-1 DPs *f* that yield first integrals of the form . In the two tables, *i*,*j* and *i*,*j*,*k* are arbitrary permutations of 1,2 and 1,2,3, respectively; and *f*=0 is an invariant line and plane, respectively. For both *d*=2 and 3, the cases when ≤*d* restrictions on the parameters yield a first integral of this form split into three types: I, II and III. (The numbering here differs from that of Cairó & Feix [29] but agrees with that of Hua *et al*. [39].) The primed cases can be viewed as projectively transformed versions of the unprimed ones. For instance, when *d*=2 the case *II*_{i}′ is related to the case II by .

The type-I cases are degenerate in the sense that the first integral *I* does not involve any additional DP *f* and is merely of the form . In fact, in each type-I case the restrictions on the parameters, expressed determinantally in terms of the *d*+1 column vectors , *j*=0,…,*d*, require that the matrix be singular, i.e. . Lotka–Volterra models with singular A (e.g. ones in which A is antisymmetric and *d* is odd) may have physical applications, though their relevance to population ecology has been strongly questioned [40,2]. But they will not be considered further here.

For *d*=2, in case II (which is defined by *a*_{10}=*a*_{20}=0, i.e. by the intrinsic growth rates of the two species being zero), the resulting first integral |*x*_{1}|^{l1}|*x*_{2}|^{l2}|*f*|^{l3} is easily seen to be
1.3
where . For *d*=3, the similar case *II*_{i} is defined by *a*_{j0}=*a*_{k0}=:*a*_{*0}, i.e. by the two species other than species *i* having equal growth rates, and by *a*_{ji}=*a*_{ki}=:*a*_{*i}, i.e. by their being affected equally by species *i*. A third, more obscure, condition *a*_{i0}*a*_{*i}=*a*_{*0}*a*_{ii} must also be satisfied. (It will be satisfied if all intrinsic growth rates are zero.) The resulting first integral is
1.4
The *d*=3 expression (1.4) reduces to a power of the *d*=2 first integral (1.3) when (*i*,*j*,*k*)=(3,1,2) and the third species decouples: *a*_{*i}=0.

The tables also list restrictions of type III, which are expressed in terms of the quantities
1.5
But this paper will focus on *d*=2 and especially *d*=3 systems that satisfy the less complicated restrictions of case II and case *II*_{i}, respectively.

It must be mentioned that not all Lotka–Volterra systems with *d*=2,3 and no self-interactions (*a*_{ii}=0 for *i*=1,…,*d*) fit into the DP framework of tables 1 and 2. The salient example is the classical *d*=2 model, which is defined by the two restrictions *a*_{11}=*a*_{22}=0 and is quite anomalous. It is integrable, but its first integral has the unusual form
1.6
A framework into which this model fits is only now being developed [41]. Each factor *f*_{k} in (1.6) satisfies for some degree-1 polynomial cofactor *K*_{k} and is therefore a DP or generalized DP, sometimes called a ‘second integral’ [33]. To understand (1.6), one needs the concept of the *multiplicity* of an invariant line (or algebraic curve) of a polynomial vector field [42]. The DPs *f*_{1}=*x*_{1}, *f*_{2}=*x*_{2} are associated with the invariant lines *x*_{1}=0, *x*_{2}=0, which are generically not multiple. The exponential factor *f*_{3} in (1.6) comes from the line ‘at infinity’, which like *x*_{1}=0 and *x*_{2}=0 is invariant in any Lotka–Volterra model but which when *a*_{11}=*a*_{22}=0 has a non-trivial (double) multiplicity. It therefore gives rise to an extra (generalized) DP, namely *f*_{3}.

The extent to which imposing parametric restrictions on a *d*=2 Lotka–Volterra system can produce single or multiple invariant lines, thereby engendering DPs or generalized DPs, is now fully understood, thanks to the recent work of Schlomiuk & Vulpe [43]. The number of invariant lines, counted with multiplicity, can if finite be as large as six (including the one at infinity). The line configuration 4.18 of Schlomiuk & Vulpe [43] corresponds to the classical model (*a*_{11}=*a*_{22}=0), and the configurations 4.5 and 4.1 to the cases II and III of table 1. But an integrability study of their many other configurations, most defined by severe parametric restrictions and many including invariant lines with non-trivial multiplicity, remains to be carried out. An extension to *d*=3 of their *d*=2 analysis, which is already quite lengthy, is also not yet available. The many possible configurations of invariant planes have not been fully classified, even if partial results have been obtained [31,34].

However, there is already a substantial literature on the integrability properties of *d*=3 Lotka–Volterra systems without self-interactions [35,32,44,45]. These are called ABC or *A*_{1}*A*_{2}*A*_{3} systems, since if *a*_{11}=*a*_{22}=*a*_{33}=0 and *a*_{ij}*a*_{jk}*a*_{ki}≠0 for some permutation *ijk* of 123, permuting species and scaling columns will yield an interaction matrix
1.7
If an ABC system satisfies *A*_{1}*A*_{2}*A*_{3}=−1 or *A*_{i}=1 for some *i*, and its intrinsic growth rates are suitably constrained, a first integral for it can be constructed from degree-1 DPs. This is because such a system will be of type I or type II, as defined in table 2. But there are exotic ABC systems with first integrals based on DPs of degree 2 or greater, or on generalized DPs. Being tightly restricted parametrically, they are not covered by table 2. Some of these systems were first obtained by a Painlevé analysis [46]. To some of these exotic ABC systems exotic *d*=2 systems are associated, as any *d*=3 system with equal intrinsic growth rates can be reduced to a *d*=2 system, typically with unequal rates [47].

A *d*=3 system with even tighter parametric restrictions is the May–Leonard model of cyclic competition [12]. This is a system with equal intrinsic growth rates (normally taken to be +1, though reducible to zero by a change of variables), and with the interaction matrix
1.8
It has 1-species, 2-species and 3-species fixed points, the reference being to the number of components of *x* with non-zero values. The large-time () behaviour of this system on the positive octant is well known. If *α*+*β*∈(−1,2) the 3-species fixed point (at which *x*_{1}=*x*_{2}=*x*_{3}=(1+*α*+*β*)^{−1}) is stable, and if alternatively *α*,*β*>1, the 1-species fixed points (1,0,0),(0,1,0),(0,0,1) are stable. The integrability properties of this system follow from table 2. If, for example, *α*+*β*=−1 or *α*=*β*=1, then and the system is of the degenerate type I. If the common intrinsic growth rate *a*_{10}=*a*_{20}=*a*_{30}=:*a*_{*0} is reduced to zero, and moreover *α*=*β* so that the system is fully rather than merely cyclically symmetric, then it is of type II, and specifically it is case *II*_{i} for each *i*. Hence the *α*=*β* case of the May–Leonard model is completely integrable: for each permutation *ijk* of 123, the expression (1.4) will be a first integral (cf. [48]). Normalized, this is
1.9
Each *I*_{i} is time-independent if *a*_{*0}=0, and depends exponentially on *τ* if *a*_{*0} is non-zero. The triple *I*_{1},*I*_{2},*I*_{3} are not functionally independent. Other first integrals are known [49,50].

The relationship between the integrability (partial or complete) of a differential system and its having the PP, i.e. the property that each of its solutions *x*=*x*(*τ*) has a one-valued continuation to the complex *τ*-plane (without branch points of any order), is a bit murky. It is widely believed that the integration of any system with the PP should reduce to quadratures; or, if not, that the solutions of the system should be expressible in terms of ‘known’ functions, such as the Painlevé transcendents that arise as solutions of second-order ODEs with the PP. Like the *d*=3 ABC system, the general *d*=2 Lotka–Volterra system has been subjected to Painlevé analyses [39,51]. It appears that the only such systems with the PP also have DPs from which a first integral can be constructed; thus they are at least partially integrable. But a Painlevé analysis of general *d*=3 Lotka–Volterra systems has yet to be performed. And, until now, no explicit solution involving Painlevé transcendents of any *d*=3 Lotka–Volterra system seems to have been published.

### (c) Overview of following results

In §2, we illustrate our methods by integrating the above case II of the *d*=2 Lotka–Volterra system with the aid of a ‘new time’ *t*. The resulting expressions for *x*_{1},*x*_{2} and *τ* as functions of *t*, involving an incomplete beta function, are new. This approach is equivalent to that of Goriely [36], but unlike him we do not posit a fixed form for the new time transformation *t*=*t*(*τ*), or make explicit use of the first integral (1.3). Non-zero but equal intrinsic growth rates are also treated.

In §3, we obtain our central result: the general solution of case *II*_{i} of the *d*=3 Lotka–Volterra system, initially requiring that the growth rates *a*_{10},*a*_{20},*a*_{30} be equal. The solution appears as (3.10), which expresses *x*_{1},*x*_{2},*x*_{3} and *τ* as functions of a new time *t*. The new time as a function of the old turns out to satisfy a generalized Schwarzian equation (gSE), which is integrated with the aid of the incomplete beta function. No explicit use is made of the first integral (1.4). Many systems are solved in closed form, including certain ABC and May–Leonard ones.

Section 3 relates our solution of case-*II*_{i} *d*=3 systems to the results of Carton-LeBrun [26], who classified all gSEs with the PP. Her results yield a classification of equal growth rate case-*II*_{i} systems with this property, and each has an elementary or an elliptic general solution *x*=*x*(*τ*). In §4, we express the general solutions of several *d*=3 systems with unequal growth rates in terms of Painlevé transcendents. These may be the first such solutions ever obtained.

In §5, we summarize the results of §§2–4, and make some final remarks.

## 2. Two-dimensional integration

To show how certain Lotka–Volterra systems of the form (1.1) can be integrated parametrically with the aid of a new time variable *t*, consider a case-II *d*=2 system of table 1: one with zero intrinsic growth rates (*a*_{10}=*a*_{20}=0). If the interaction matrix A satisfies *a*_{11}≠*a*_{21} and *a*_{22}≠*a*_{12}, by scaling (redefining) the components *x*_{1},*x*_{2}, one can scale the columns of A so that
2.1
for some *a*_{1},*a*_{2}. The resulting rather symmetric system
2.2
has four invariant lines: the usual *x*_{1}=0, *x*_{2}=0 and the one at infinity, which plays no direct role here; also *x*_{1}−*x*_{2}=0, since if *x*_{1}=*x*_{2} then . Only solutions *x*=*x*(*τ*) not lying along any of these invariant lines will be considered. By examination,
2.3
is a DP-based first integral: the reciprocal of the standard case-II first integral (1.3).

Define an auxiliary (new time) variable *t* by
2.4
so that correspond to the three just-mentioned finite invariant lines. The three *t*-intervals , (−1,1), correspond to sectors in the *x*_{1},*x*_{2}-plane lying between the lines. Thus, in the first quadrant, the sectors 0<*x*_{2}<*x*_{1} and 0<*x*_{1}<*x*_{2} correspond to the *t*-intervals and . To any solution *x*=*x*(*τ*) of system (2.2), there is associated a function *t*=*t*(*τ*), taking values in one of the three *t*-intervals. The value *t*_{0} taken by *t* at the time-origin (say, at *τ*=0) determines the *t*-interval.

By calculus applied to (2.4) and system (2.2), *t*=*t*(*τ*) satisfies
2.5
which integrates to the hyper-logistic growth law
2.6
*A*≠0 being arbitrary. Hence viewed inversely as a function of *t*, *τ* is given by
2.7
Moreover, it follows by differentiating (2.4) and exploiting (2.2) that
2.8
Thus, *x*_{1},*x*_{2} as well as *τ* can be expressed as functions of the new time *t*.

The integral in (2.7) defines an increasing function of *t*: an *incomplete beta function* *τ*=*B*_{a1,a2;t0}(*t*) with parameters *a*_{1},*a*_{2}, as is explained in appendix A. It maps in an increasing way the *t*-interval containing *t*_{0} onto some *τ*-interval, which may be infinite. System (2.2) is thus solved parametrically by
2.9a
and
2.9b
(2.9b) coming from (2.8). The parameter *t* is restricted to the relevant *t*-interval.

This is a *complete* integration of the *d*=2 system (2.2), as the expressions for *x*_{1},*x*_{2} and *τ* as functions of the new time *t* involve two undetermined constants: *A* and *t*_{0}. (Other than determining the *t*-interval, the latter merely shifts *τ*.) Incomplete beta functions are supported by many software packages, so the numerical solution of any case-II *d*=2 system is easy. Growth rates *a*_{10},*a*_{20} that are non-zero but equal can readily be incorporated; see example 2.1.

One can also eliminate *t*, and at least formally express *x*_{1},*x*_{2} as functions of the original time *τ*. Let be any convenient, standardized solution of the *A*=1 case of the nonlinear ODE (2.6). (The function will map in an increasing way some *τ*-interval (*τ*_{min},*τ*_{max}), which may be infinite, onto one of the *t*-intervals , (−1,1), ; so it will depend on the choice of *t*-interval.) The general solution *x*=*x*(*τ*) of (2.2) is then given by (2.8) or (2.9b) with
2.10
which is defined if *A*(*τ*−*τ*_{0}) lies in (*τ*_{min},*τ*_{max}). The solution *x*=*x*(*τ*) thus contains two free parameters, *A* and *τ*_{0}, and implicitly a choice of *t*-interval, i.e. a choice of sector in the *x*_{1},*x*_{2}-plane.

For this *x*=*x*(*τ*) to be more than formal, an explicit expression for the (standardized) inverse incomplete beta function is needed. As appendix A explains, expressions are available for some *a*_{1},*a*_{2}. If the set {1/*a*_{1},1/*a*_{2},1/(1−*a*_{1}−*a*_{2})} is any of {1,*m*,−*m*} (for *m* a positive integer), or the inverse is expressible in terms of elementary functions; and if it is any of {2,4,4}, {2,3,6} or {3,3,3}, in terms of elliptic ones. In each case, the inverse extends to a one-valued function on the complex *τ*-plane. Hence *x*=*x*(*τ*) given by (2.8) is also one-valued, and Lotka–Volterra system (2.2) has the PP.

Table 3 gives a closed-form expression for the (standardized) inverse function in the cases . Each is one-valued on the *τ*-plane, as stated. It should also be noted that if *a*_{1},*a*_{2} are integers which are positive, or satisfy *a*_{1}*a*_{2}<0 with *a*_{1}+*a*_{2}≤0, the function will be an *algebraic* function (and hence, finite-valued on the *τ*-plane). By (2.8), *x*=*x*(*τ*) will also be finite-valued. In any of these finite-valued cases, Lotka–Volterra system (2.2) has by definition a weak form of the PP.

### Example 2.1

Consider the symmetric system
2.11
where *x*_{1},*x*_{2}>0 and *x*_{1}≠*x*_{2} for simplicity. The state variables *x*_{1},*x*_{2} are the populations of species that grow logistically (if *a*_{*0}>0, that is), and also display mutualism: the growth of each is made more rapid by the presence of the other. Consider first the case when *a*_{*0}=0, which is the case of the system (2.2).

The auxiliary variable *t*=(*x*_{1}+*x*_{2})/(*x*_{2}−*x*_{1}) satisfies the interval condition if (initially and hence subsequently) 0<*x*_{1}<*x*_{2}; suppose this to be so. By table 3, the standardized inverse incomplete beta function for the interval is
2.12
Here cn is the Jacobian elliptic function with modular parameter *m*=*k*^{2} equal to the ‘lemniscatic’ value . The domain (*τ*_{min},*τ*_{max}) is (0,*K*_{1/4}), where . Substituting (2.12) into (2.8) yields the general 0<*x*_{1}<*x*_{2} solution, parametrized by *A*>0 and *τ*_{0}; namely,
2.13
This is defined if *A*(*τ*−*τ*_{0}) lies in (0,*K*_{1/4}). In deriving (2.13) from (2.12), such elliptic-function identities as sn^{2}+cn^{2}=1 and *m* sn^{2}+*dn*^{2}=1 have been used. For this system, the first integral (2.3) is *I*=|*x*_{1}*x*_{2}|^{1/4}|*x*_{1}−*x*_{2}|^{1/2}, and by examination .

The phase portrait in the 0<*x*_{1}<*x*_{2} sector, comprising solutions of the form (2.13), is easily worked out. The endpoints *τ*_{min}=0, *τ*_{max}=*K*_{1/4} are zeroes of sn, cn, respectively. So as *A*(*τ*−*τ*_{0})→0^{+}, *t*→1^{+} and (*x*_{1},*x*_{2}) diverges to infinity, i.e. it approaches the invariant line at infinity. It simultaneously approaches the positive half of the *x*_{2}-axis (a finite invariant line). Similarly, as , and (*x*_{1},*x*_{2}) diverges to infinity while also approaching the *x*_{1}=*x*_{2} line (another invariant line). The divergences are due to the strong mutualism. Although each solution (2.13) has a real *τ*-interval of width *A*^{−1}*K*_{1/4} as its domain, it extends in a one-valued way to the complex *τ*-plane. Hence system (2.11) has the PP.

Incorporating a non-zero common growth rate *a*_{*0} is straightforward, as was explained in §1. Exponentially modifying the independent and dependent variables yields the general 0<*x*_{1}<*x*_{2} solution, parametrized by real *τ*_{1} and *C*. It is
2.14
which is defined if *e*^{a*0(τ−τ1)}−*C* lies in (0,*K*_{1/4}). Compared with the mutualism, the role played by any growth rate *a*_{*0}>0 in causing finite-time blow-up is minor.

On the complex *τ*-plane, the poles of any instance of solution (2.13) lie on a square lattice, since sn,cn,*dn* are doubly periodic, but those of any instance of (2.14) lie on an exponentially stretched lattice. The pole locations are of interest, because solutions in the complex time domain of non-integrable systems are expected to have irregular patterns of singularities [52].

The just-concluded example was rather special, as setting leads to one-valuedness on the complex *τ*-plane. It must be stressed that for generic *a*_{1},*a*_{2}, owing to the presence of branch points, the general complex-domain solution (*x*_{1}(*τ*),*x*_{2}(*τ*)) of system (2.2) will not be one-valued or even finite-valued. The distribution of its poles over its many branches remains to be explored.

However, the example made clear the behaviour of any real-domain solution (*x*_{1}(*τ*),*x*_{2}(*τ*)) in the forward-time limit, i.e. as *τ* tends to the upper endpoint *τ*_{max} of its interval of definition. In the first quadrant, either the invariant line *x*_{1}=*x*_{2} will be an attractor (as in the example), or it will be a repeller and the lines *x*_{1}=0, *x*_{2}=0 will be attractors. In population ecology, if *a*_{1},*a*_{2}>1 so that the two species compete rather than display mutualism, these two possibilities would be called competitive coexistence and exclusion [40].

The analysis in this section of the case-II *d*=2 system assumed for simplicity that *a*_{11}≠*a*_{21} and *a*_{22}≠*a*_{12}. This can be relaxed. Suppose the intrinsic growth rates are zero (*a*_{10}=*a*_{20}=0) and, say, that *a*_{11}=*a*_{21}=:*a*_{*1}. Then by examination the three usual invariant lines, namely *x*_{1}=0, *x*_{2}=0 and the one at infinity, have respective multiplicities 1,2,1 in the sense of Christopher *et al*. [42]; and there are no others. (This configuration of invariant lines is numbered 4.20 by Schlomiuk & Vulpe [43].) The line *x*_{2}=0 thus has an extra (generalized) DP associated with it. The first integral (1.3) is accordingly replaced by
2.15
exhibiting the DPs *x*_{1},*x*_{2} and, as the third factor, the extra generalized DP. It is convenient to choose as the new time variable *s*:=*x*_{1}/*x*_{2}. ODE (2.5) is then replaced by
2.16
This can be integrated by inspection, giving an expression for *τ*=*τ*(*s*) in terms of the incomplete gamma function; or in the special case when *a*_{22}=−*a*_{12}, of the error function erf. This explicit integration sheds light on the recent solution of a particular Lotka–Volterra system with *a*_{11}=*a*_{21} and *a*_{22}=−*a*_{12}, involving *erf* and *erf*^{−1} [53, §4.2]. It can evidently be generalized.

## 3. Three-dimensional integration

### (a) General case-II_{i} systems

It will now be shown that an unexpectedly large family of *d*=3 Lotka–Volterra systems of the form (1.1) can be integrated parametrically by the technique of §2. These are case-*II*_{i} systems as defined in table 2. Specifically, they are *d*=3 systems with *a*_{ji}=*a*_{ki}=:*a*_{*i}, so that the species *j*,*k* other than *i* are equally affected by species *i*, the reverse not being assumed. Without loss of generality, *i*=3 and ( *j*,*k*)=(1,2) will be taken. Initially, it will be assumed that *a*_{10}=*a*_{20}=*a*_{30}=0, i.e. that the intrinsic growth rates are zero, since incorporating a non-zero common growth rate is an easy matter.

Systems in this family include certain *ABC* and symmetric May–Leonard models, which are treated as examples in §3*b*,*c*. Besides constructing explicit solutions, one can determine from powerful results of Carton-LeBrun [26] precisely which systems in the family have the PP that each solution *x*=*x*(*τ*) extends in a one-valued way to the complex *τ*-plane. Generalized case-*II*_{i} systems which have unequal growth rates will be examined in §4.

Suppose *a*_{13}=*a*_{23}=:*a*_{*3}. If the interaction matrix A satisfies *a*_{11}≠*a*_{21}, *a*_{22}≠*a*_{12} and *a*_{33}≠*a*_{*3}≠0, by scaling (redefining) the components *x*_{1},*x*_{2},*x*_{3}, one can scale the columns of A so that
3.1
for some *a*_{1},*a*_{2};*b*_{1},*b*_{2} and *n*≠0,−1. ( will signify 1/*n*=0.) The resulting system
3.2
which reduces to the *d*=2 system (2.2) when *x*_{3}≡0, has five invariant planes: the usual *x*_{1}=0, *x*_{2}=0, *x*_{3}=0 and the one at infinity, which plays no direct role here; also *x*_{1}−*x*_{2}=0, since if *x*_{1}=*x*_{2}, then . Solutions *x*=*x*(*τ*) not lying in any of these planes are of primary interest. By examination,
3.3
is the specialization of the case-*II*_{3} first integral (1.4). Here and below, we make the definition
3.4
for convenience; note that , with meaning , i.e. 1/*n*=0.

Define an auxiliary variable *t* as in §2 by
3.5
so that correspond to the planes *x*_{1}=0, *x*_{2}=0 and *x*_{1}−*x*_{2}=0. Thus, 0<*x*_{2}<*x*_{1} and 0<*x*_{1}<*x*_{2} correspond to the *t*-intervals and . The function *t*=*t*(*τ*) associated with any trajectory *x*=*x*(*τ*) takes values in one of the intervals , (−1,1), . On any portion of *x*=*x*(*τ*) on which does not change sign, *t* can be used as an alternative parameter: a new time.

One can obtain an expression for each of as a rational function of *x*_{1},*x*_{2},*x*_{3} by repeatedly differentiating (3.5) with respect to *τ*, the equations in system (3.2) being exploited after each differentiation. By the tedious process of eliminating the three variables *x*_{1},*x*_{2},*x*_{3} from these four formulae, one obtains a single identity satisfied by , which can be written as
3.6
This nonlinear third-order ODE for *t*=*t*(*τ*) is the *d*=3 counterpart of (2.5). It will be called a *gSE*, since it is a generalization: when (i.e. *n*=−2) its first two terms constitute a Schwarzian derivative, and if additionally *b*_{1}=*a*_{1}, *b*_{2}=*a*_{2} so that its two terms vanish, it reduces to a Schwarzian equation (SE) of a type familiar from the conformal mapping of triangles [54].

The seemingly complicated gSE (3.6), with solution *t*=*t*(*τ*), can be integrated without difficulty. Suppose for simplicity that . Let , i.e. , and view *t*,*f* (rather than *τ*,*t*) as the independent and dependent variables. By applying the chain rule where *D*_{t}:=d/d*t*, one finds that
3.7
is the ODE (surprisingly, a linear one) satisfied by *f*=*f*(*t*). This is a so-called Papperitz equation, and can be greatly simplified by a substitution of a standard type. Performing the substitution reduces (3.7) to a degenerate hypergeometric equation, namely
3.8
which can be integrated by inspection. In this way, one deduces that
3.9
where *K*_{1},*K*_{2} are undetermined constants and *t*_{0} is any convenient point in the relevant *t*-interval. For the incomplete beta function *B*_{b1,b2;t0}(*t*), see appendix A.

Formulae for the inverse function *τ*=*τ*(*t*) and the alternatively parametrized trajectory *x*=*x*(*t*) of system (3.2) follow immediately from (3.9). The restriction will now be dropped: suppose instead that, on the portion of the trajectory under study, . By integrating , one obtains
3.10a
and
3.10b
In (3.10a), *t* is restricted to the relevant *t*-interval (containing any conveniently chosen point *t*_{0}), and is further restricted by the requirement that the quantity *K*_{1}+*K*_{2} *B*_{b1,b2;t0}(*t*′), within |⋅|, not change sign. The accompanying formulae (3.10b), which extend the *d*=2 formula (2.8), are an easy exercise. They come by reverting the above-mentioned expressions for as rational functions of *x*_{1},*x*_{2},*x*_{3}, to produce expressions for *x*_{1},*x*_{2},*x*_{3} as functions of . In practice, one would compute in (3.10b) by applying and to (3.10a).

The parametric solution (*τ*,*x*)=(*τ*(*t*),*x*(*t*)) displayed in (3.10) is the central result of this paper. It is a *complete* integration of the *d*=3 Lotka–Volterra system (3.2), as the expressions for *τ* and *x*_{1},*x*_{2},*x*_{3} as functions of *t* involve three undetermined constants: the coefficients *K*_{1},*K*_{2}, and also *τ*_{0} (which merely shifts *τ*). Intrinsic growth rates *a*_{10},*a*_{20},*a*_{30} that are non-zero but equal can easily be incorporated in (3.10), as when *d*=2. Solution (3.10) subsumes the *d*=2 solution (2.9). This is because if *K*_{2}=0 then *τ*=*τ*(*t*) given by (3.10a) reduces to (2.9a), and *x*_{3}=*x*_{3}(*τ*) computed from (3.10b) will be identically zero.

### Example 3.1

Consider a *d*=3 Lotka–Volterra system of the form (3.2) with *b*_{1}=*a*_{1}, *b*_{2}=*a*_{2}, so that species 3 is unaffected by species 1 and 2, though the reverse does not hold. In this case, (*τ*,*x*)=(*τ*(*t*),*x*(*t*)) can be simplified, since, if *K*_{2}≠0, equation (3.10a) integrates to
3.11
This can be converted to a formula for the inverse function *t*=*t*(*τ*). In terms of the standardized inverse beta function , as discussed in appendix A (and used in §2), the formula is
3.12
for some *A*≠0,*B*,*C*. For *a*_{1},*a*_{2} for which is expressible in terms of standard functions, such as the cases when *a*_{1}=*a*_{2}=:*a* that are listed in table 3, this allows *x*=*x*(*τ*) when *b*_{1}=*a*_{1}, *b*_{2}=*a*_{2} to be similarly expressed, by substituting (3.12) into (3.10b).

For many *a*_{1},*a*_{2}, such as those in the table, the function has a one-valued extension to the complex *τ*-plane. For such *a*_{1},*a*_{2}, if *n* is an *integer* (*n*≠0,−1) or , system (3.2) with *b*_{1}=*a*_{1}, *b*_{2}=*a*_{2} will have the PP. This is because, if *n* is an integer or , *t*=*t*(*τ*) in (3.12) and hence *x*=*x*(*τ*) will extend to the *τ*-plane with no branch-point at *τ*=−*B*/*A*, i.e. in a one-valued way. But in each such system, the third species decouples.

One can also determine which *d*=3 Lotka–Volterra systems of the form (3.2), *with no species decoupled*, have the PP; and, in fact, integrate each such explicitly. In gSE (3.6) for *t*=*t*(*τ*), there are seven terms and thus seven coefficients, the first two being 1 and . Suppose the remaining five are not determined by (*a*_{1},*a*_{2};*b*_{1},*b*_{2};*n*), as here, but are free to vary. For some choices of the remaining five, the resulting nonlinear ODE will have the PP. By a mathematically rigorous analysis Carton-LeBrun [26] determined, classified and tabulated the many possibilities, and integrated each resulting ODE. Hence all one needs to do is compare the coefficients in gSE (3.6) against her tables to determine for which (*a*_{1},*a*_{2};*b*_{1},*b*_{2};*n*) all solutions *t*=*t*(*τ*) of the gSE will extend in a one-valued way to the *τ*-plane. For each such choice, one can compute *x*=*x*(*τ*) from (3.10b), and it will extend similarly; thus the case-*II*_{3} system (3.2) as well as the gSE will have the PP.

It is not hard to see that, for gSE (3.6) to have the PP, the parameter *n* must be an integer (*n*≠0,−1 as usual) or . Otherwise, the gSE would have complex-domain solutions *t*=*t*(*τ*) with movable branch-points. The tables of Carton-LeBrun include infinite families of gSEs in which *n* is otherwise unconstrained, such as ones with *b*_{1}=*a*_{1}, *b*_{2}=*a*_{2} that were, in fact, implicitly encountered in example 3.1. They also include many gSEs with *n*=1,2,3,5 or . As she shows, each gSE with the PP can be integrated twice to yield a polynomial identity with free parameters *K*_{1},*K*_{2}, the possible identities being listed in her table VIII. Any solution *t*=*t*(*τ*) is thus reduced to quadratures. For nearly all of these integrable gSEs, solutions *t*=*t*(*τ*) and trajectories *x*=*x*(*τ*) can be expressed in terms of elliptic functions (e.g. lemniscatic or equianharmonic ones); for a few, in terms of elementary functions.

### Example 3.2

Consider the *d*=3 Lotka–Volterra system
3.13
where for simplicity *x*_{1},*x*_{2},*x*_{3}>0 and *x*_{1}≠*x*_{2}. The state variables *x*_{1},*x*_{2},*x*_{3} are the populations of three species with a common intrinsic growth rate *a*_{*0}. If *a*_{*0}>0, then in isolation species 1 and 2 grow exponentially, and species 3 grows logistically. Species 1 and 2 display mutualism, as do 1 and 3; but species 2 preys on species 3.

Since a non-zero *a*_{*0} is easy to incorporate, set *a*_{*0}=0. System (3.13) is then the case (*a*_{1},*a*_{2};*b*_{1},*b*_{2};*n*)=(0,0;1,−2;1) of (3.2). With these values, gSE (3.6) for *t*=*t*(*τ*) is a variant of Carton-LeBrun’s *n*=1 case II. Integrating twice yields a polynomial identity , i.e.
3.14
where *K*_{1},*K*_{2} are undetermined constants. (Note that as , the right side is the general solution of the Papperitz equation (3.7).) Integrating once more yields
3.15
as the general solution of the gSE. Here *c*_{1},*c*_{2} are constrained by , and *τ* can be replaced by *Aτ*+*B*, for any *A*≠0 and *B*.

Substituting (3.15) into (3.10b) yields
3.16
as the general solution *x*=*x*(*τ*) of system (3.13) with *a*_{*0}=0, it being understood that *τ* can be replaced by *Aτ*+*B*, with *x* scaled by *A*. There are three undetermined constants: *A*, *B* and the point on the hyperbola . Any solution *x*=*x*(*τ*) of this form extends to a one-valued function on the complex *τ*-plane, so system (3.13) has the PP.

The first integral *I* of (3.3) specializes to , and one can show *I*≡|*A*|^{2}×|*c*_{1}+1|/2. By examination, is another first integral, and *J*≡−*A*^{2}. The phase portrait is clear from (3.16), the denominators in which can pass through zero as *τ* varies. Each trajectory *x*=*x*(*τ*) is defined on a real *τ*-interval (*τ*_{min},*τ*_{max}). As *τ*→*τ*^{+}_{min} or *τ*→*τ*^{−}_{max}, the state *x* diverges to infinity, i.e. it approaches the invariant plane at infinity. It simultaneously approaches one of the four *finite* invariant planes *x*_{1}=0, *x*_{2}=0, *x*_{3}=0, *x*_{1}−*x*_{2}=0.

### Example 3.3

Consider the *d*=3 Lotka–Volterra system
3.17
where *r*≠0, and for simplicity *x*_{1},*x*_{2},*x*_{3}>0 and *x*_{1}≠*x*_{2}. System (3.13) is the case (*a*_{1},*a*_{2};*b*_{1},*b*_{2};*n*)=(1/*r*,−1/*r*;4/*r*,−4/*r*;1) of (3.2). With these values, gSE (3.6) for *t*=*t*(*τ*) is a variant of Carton-LeBrun’s case e.IV. The gSE can be integrated twice to yield *t*=(*u*^{r}+1)/(*u*^{r}−1), where
3.18
with *K*_{1},*K*_{2} undetermined. (To confirm this, it is easiest not to integrate the gSE itself, but rather the Papperitz equation (3.7) or preferably the degenerate hypergeometric equation (3.8), which are equivalent.) By integrating (3.18) once more, it follows that *u*∝cn(*τ*) and thus that
3.19
Here cn denotes the lemniscatic variant of the Jacobian elliptic function cn, which has modular parameter , and *C* is arbitrary; and as above, *τ* can be replaced by *Aτ*+*B*.

Substituting expression (3.19) into (3.10b) yields
3.20
as the general solution *x*=*x*(*τ*) of system (3.17), it being understood that *τ* can be replaced by *Aτ*+*B*, with *x* scaled by *A*. There are three undetermined constants: *A*, *B* and *C*. If *r* is an integer then (3.20) like (3.16) will be one-valued on the complex *τ*-plane; hence if *r* is a non-zero integer, system (3.17) like (3.13) will have the PP.

The first integral *I* of (3.3) specializes to |*x*_{1}/*x*_{2}|^{2/r}|*x*_{3}(*x*_{1}−*x*_{2})|, and one can show *I*≡|*A*|^{2}×|*r*||*C*|^{−2/r}. By examination, *J*=(*x*_{1}−*x*_{2})^{2}*x*_{3}(2*x*_{1}−2*x*_{2}−*rx*_{3}) is another first integral, and *J*≡*A*^{4}×(−*r*^{3}). In the real *τ*-domain, each solution of the form (3.20) has an interval of definition (*τ*_{min},*τ*_{max}), on which its qualitative behaviour resembles the behaviour in example 3.2.

### (b) ABC systems

An interesting application of the results of the last subsection is to a well-known family of ABC systems, which as explained in §1*b* are (generic) *d*=3 Lotka–Volterra systems without self-interactions. Any ABC system, with interaction matrix given by (1.7), is of the form
3.21
It is known that when one of *A*_{1},*A*_{2},*A*_{3} equals unity the ABC system is completely integrable in the sense that a pair of first integrals exists, one of them being of the DP type [35]. But no general solution of ABC systems of this type has previously been published. From the present point of view, the integrability comes from an ABC system with *A*_{i}=1 being a case-*II*_{i} system as defined in table 2. This facilitates a complete integration.

Suppose without loss of generality that it is *A*_{3} which equals unity, and also suppose that *A*_{2}≠0. In terms of a scaled state vector (*x*_{1},*x*_{2},*x*_{3}):=(*y*_{1},*A*_{2}*y*_{2},*y*_{3}), the system becomes
3.22
which is the case of the case-*II*_{3} system (3.2). An integration of (3.22) therefore follows at once from the complete integration (3.10). (As above, solutions that do not lie in any of the invariant planes *x*_{1}=0, *x*_{2}=0, *x*_{1}−*x*_{2}=0 are of primary interest.)

The new time *t*:=(*x*_{1}+*x*_{2})/(*x*_{2}−*x*_{1}) is used as a parameter. Along any segment of a trajectory *x*=*x*(*τ*) on which , *τ* and *x*=(*x*_{1},*x*_{2},*x*_{3}) are expressed in terms of *t* by
3.23a
and
3.23b
Here, *τ*_{0} is the initial time, at which *t*=*t*_{0}. The parameter *t* is restricted to whichever of the intervals , (−1,1), contains *t*_{0}. It is further restricted by the condition that , within |⋅|, not change sign over the integration interval. The derivatives in (3.23c) are computed by applying and to (3.23a).

The parametric general solution (3.23) is a complete integration of the ABC system with *A*_{3}=1, since the expressions for *τ* and *x*_{1},*x*_{2},*x*_{3} involve three undetermined constants: *K*_{1},*K*_{2} and also *τ*_{0} (which merely shifts *τ*). For many choices of *A*_{1},*A*_{2}, the expressions for *τ* and *x*_{1},*x*_{2},*x*_{3} in terms of *t* can be expressed using elementary functions. This is because, for many choices of *b*_{1},*b*_{2}, the incomplete beta function *B*_{b1,b2;t0}(*t*) in (3.23a) can be so expressed. The evaluation of *τ* and *x*_{1},*x*_{2},*x*_{3} is especially easy if *B*_{b1,b2;t0}(*t*) is a polynomial or rational function of *t*.

### Example 3.4

Consider the case (*A*_{1},*A*_{2},*A*_{3})=(1,1,1) of the ABC system (3.21) and hence of the normalized system (3.22), for which . In this case, *x* and *y* are identical, and the system is simply
3.24
where *ijk* signifies each cyclic permutation of 123. This ABC system has been previously studied [55, §7], though not at length. The gSE (3.6) with does not appear on Carton-LeBrun’s list, so the system does not have the PP. But a complete integration is provided by the parametric solution (3.23).

To illustrate the parametrization, consider the trajectories *x*=*x*(*τ*) in the region 0<*x*_{1}<*x*_{2}, 0<*x*_{3}, along each of which *t*>1. Definition (A.2) implies that *B*_{1,1;t0}(*t*) equals *t* plus a constant. Substituting this into (3.23) and evaluating the integral and derivatives yields
3.25a
and
3.25b
where *C* is a constant of integration; and one can also replace *τ* by *Aτ*+*B* and scale (*x*_{1},*x*_{2},*x*_{3}) by *A*. Thus, there are three undetermined constants in all.

Any trajectory *x*=*x*(*τ*) with 0<*x*_{1}<*x*_{2} and 0<*x*_{3} can be traced out parametrically by varying *t* over (*t*_{min},*t*_{max}), which is the smaller of and . As *t*→*t*^{+}_{min}, *x* tends to a point on the *x*_{2}-axis if *C*<1, and to one on the *x*_{3}-axis if *C*>1. (If *C*=1, the trajectory lies in the *x*_{2}−*x*_{3}=0 plane.) As , *x* diverges to infinity. The first integral *I* of (3.3) specializes to , and by examination *I*≡|*A*|×4, irrespective of the constant *C*. There are other first integrals, such as the triple *J*_{1},*J*_{2},*J*_{3} defined by *J*_{i}=(*x*_{i}−*x*_{j})(*x*_{k}−*x*_{i})/*x*_{i}.

To visualize the global dynamics on , it is better to use a parametrization suggested by but more symmetric than (3.25). Note first that owing to the cyclic symmetry there are six finite invariant planes: *x*_{1}=0, *x*_{2}=0, *x*_{3}=0, *x*_{1}−*x*_{2}=0, *x*_{2}−*x*_{3}=0, *x*_{3}−*x*_{1}=0. Any trajectory *x*=*x*(*τ*) not lying in any of these planes lies entirely in one of the sectors between them. To parametrize it by an ad hoc free variable *s*, rather than by *t*, use
3.26a
and
3.26b
where ± is used if ±*x*_{1}*x*_{2}*x*_{3}>0 in the sector. The quantities *s*_{1},*s*_{2},*s*_{3} are taken to be distinct and to satisfy *s*_{1}+*s*_{2}+*s*_{3}=0; they determine the trajectory, along which *J*_{1},*J*_{2},*J*_{3} are constant. It is easily verified that, as a function of *τ*, (*x*_{1},*x*_{2},*x*_{3}) formally satisfies ABC system (3.24).

The real *s*-line is divided into four intervals by *s*_{1},*s*_{2},*s*_{3}, and to trace out a single trajectory the new time *s* in (3.26) should range over only one of the four. In (3.26b), *s*=*s*_{1},*s*_{2},*s*_{3} yield points *x*=(*x*_{1},*x*_{2},*x*_{3}) that lie on the coordinate axes, and yield infinite points. Hence any generic solution *x*=*x*(*τ*) extends in , as *τ* varies over its interval of definition, either between a pair of points on distinct axes or between a point on an axis and one at infinity.

It turns out that each trajectory in follows a plane curve, which is generically a (segment of a) parabola. This is because *x*=(*x*_{1},*x*_{2},*x*_{3}) given by (3.26b) satisfies both
3.27a
and
3.27b
which are the equations of a plane and a cone. The trajectory is therefore a segment of a conic section, which is easily determined to be a parabola. This observation is new.

The general solution (3.26) is equivalent to but is much more explicit than a solution of the ABC system with (*A*_{1},*A*_{2},*A*_{3})=(1,1,1) that was obtained in [55, §7]. By manipulating (3.26b), it follows that the new time *s* used here as an alternative to *t* has the definition
3.28
The presence of the square root, like the logarithms in (3.26a), may be an underlying reason for this ABC system not having the PP.

For any ABC system in the scaled form (3.22), the new time *t*=*t*(*τ*) satisfies gSE (3.6) with . Several of the gSEs listed by Carton-LeBrun [26] as having the PP, meaning that all solutions *t*=*t*(*τ*) extend to a one-valued way to the complex *τ*-plane, have and *a*_{1}=*a*_{2}=0. They therefore give rise to ABC systems with certain *A*_{1},*A*_{2} and *A*_{3}=1 that have the PP.

### Example 3.5

Consider the case of ABC system (3.21) and hence of (3.22), for which . This case has been previously investigated [46], though not at length. With these values for (*a*_{1},*a*_{2};*b*_{1},*b*_{2};*n*), gSE (3.6) satisfied by *t*=*t*(*τ*) is a variant of Carton-LeBrun’s case e.IV.3. The solution of this gSE is
3.29
where *C*≠0 is an arbitrary constant of integration and *τ* can be replaced by *Aτ*+*B*. (This is best shown by integrating the degenerate hypergeometric equation (3.8), though one could also integrate equation (3.23a), as in the last example.) Substituting (3.29) into (3.23c) gives
3.30
where again *τ* can be replaced by *Aτ*+*B*, in which case (*x*_{1},*x*_{2},*x*_{3}) is scaled by *A*. Thus, there are three undetermined constants in all.

Equation (3.30) provides a complete integration of the ABC system with . If −1<*C*<1, the trajectory defined by (3.30) with lies in the positive octant; if *C*=0, in the plane *x*_{1}=*x*_{2}; and if *C*=±1, in the plane *x*_{3}=0. The first integral *I* of (3.3) specializes to |*x*_{1}*x*_{2}|^{1/2}|*x*_{3}/(*x*_{1}−*x*_{2})|, and by examination *I*≡|*A*|×|(*C*^{2}−1)/2*C*|.

This system is of much interest, as it is the *A*_{3}=1 member of a certain family of ABC systems (with *A*_{1}=−1/(*A*_{3}+1) and *A*_{2}=−(*A*_{3}+1)/*A*_{3} for *A*_{3}≠0,−1), which as shown by Bountis *et al*. [46] have the PP and are completely integrable. Each member has a first integral that is quadratic in *x*_{1},*x*_{2},*x*_{3} [48]. For the *A*_{3}=1 member, this extra first integral is
3.31
and by examination *J*≡*A*^{2}×16, irrespective of the constant *C*. Each trajectory provided by (3.30) lies along the intersection of a pair of surfaces *I*≡*const*. and *J*≡*const*., which are, respectively, a quartic surface and a hyperboloid. As in example 3.4, that this ABC system has first integrals that are elementary functions of *x*_{1},*x*_{2},*x*_{3} facilitates an understanding of its global dynamics.

Bountis *et al*. remark that the ABC system with can be integrated in terms of elliptic functions; but as (3.30) reveals, and suffice. This general solution also makes manifest the one-valuedness of the continuation of *x*=*x*(*τ*) to the complex *τ*-plane.

### (c) Symmetric May–Leonard systems

A particularly striking application of the results of §3*a* is to the integration of certain May–Leonard systems. As was explained in §1*b*, the three species in any May–Leonard system compete cyclically, with *α*,*β* being the strength of clockwise and counter-clockwise interactions. The systems of interest here are those with *α*=*β*, in which the species are equivalent. A general integral *y*=*y*(*τ*) will be obtained for any *α*=*β*≠0,1: see equations (3.34)–(3.35). The system with *α*=*β*=−1, which is mutualistic rather than competitive, will be solved in terms of elliptic functions and shown to have the PP. May–Leonard systems have been analysed from the Painlevé point of view [56,57], and many first integrals have been discovered [49,50]. But the results below are much stronger.

It is worth recalling the long-time () behaviour of the system state *y*=(*y*_{1},*y*_{2},*y*_{3}) of the *α*=*β* May–Leonard model, in the positive octant. If the common intrinsic growth rate *a*_{10}=*a*_{20}=*a*_{30}=:*a*_{*0} equals +1, the 3-species fixed point has *y*_{1}=*y*_{2}=*y*_{3}=1/(1+2*α*). The state *y* converges to this point if , though if *α*>1 it is the 1-species fixed points which are stable.

It will be assumed below that *a*_{*0}=0, since restoring a non-zero common growth rate is accomplished by an easy change of variables. Lotka–Volterra system (1.1), with the May–Leonard interaction matrix (1.8) and *α*=*β*, is then of the form
3.32
Suppose *α*≠0,1, and define a scaled state vector (*x*_{1},*x*_{2},*x*_{3}):=((1−*α*)*y*_{1},(1−*α*)*y*_{2},−*αy*_{3}). The system can then be written as
3.33
which is the (*a*_{1},*a*_{2};*b*_{1},*b*_{2};*n*)=((1−*α*)^{−1},(1−*α*)^{−1};1,1;−*α*) case of the case-*II*_{3} system (3.2). An integration scheme therefore follows at once from the complete integration (3.10). As above, it will yield all solutions that do not lie in any of the invariant planes *x*_{1}=0, *x*_{2}=0, *x*_{1}−*x*_{2}=0.

The new time variable *t*:=(*x*_{1}+*x*_{2})/(*x*_{2}−*x*_{1}) plays its familiar role. Along any segment of a scaled trajectory *x*=*x*(*τ*) on which , the function *t*=*t*(*τ*) satisfies
3.34
where *K*_{1},*K*_{2} are constants. This is a specialization of (3.9), because, according to definition (A.2), the relevant incomplete beta function *B*_{1,1;t0}(*t*) equals *t* plus a constant. One would obtain *t*=*t*(*τ*) by integrating equation (3.34), the maximal *τ*-interval of definition being an interval on which the quantity within |⋅| does not change sign. One would substitute *t*(*τ*) into (3.10b), i.e. into
3.35
to obtain the components of *x*=*x*(*τ*), and hence those of *y*=*y*(*τ*). This scheme is readily carried out, either numerically or (for certain *α*) symbolically. Thus, one can, *inter alia*, confirm the results of Blé *et al*. [58] on the global dynamics of *α*=*β* May–Leonard systems. The following two examples show what is involved in the construction of a closed-form general solution.

### Example 3.6

Suppose *α*=*m*/(*m*+1) with *m* a positive integer, so that . Then the exponent *α*/(*α*−1) in (3.34) equals −*m*, and is a degree-3*m* polynomial in *t*. Hence *τ* is a degree-(3*m*+1) polynomial, and its inverse *t*=*t*(*τ*) is an algebraic function. By substituting this function into (3.35), one obtains (*x*_{1},*x*_{2},*x*_{3})(*τ*) in the case when *a*_{*0}=0. For instance, if *m*=1 and , then *τ* is a quartic polynomial in *t*. By applying the quartic formula, one can express *t*=*t*(*τ*) and hence each *x*_{i}=*x*_{i}(*τ*) in terms of radicals. (The expressions are complicated and are not given here.) A radical representation is not possible when *m*>1, i.e. when etc.

But even if *t*=*t*(*τ*) for these values of *α* cannot be expressed in terms of radicals, it is algebraic, so each component *x*_{i} (computed from (3.35)) is also an algebraic function of *τ*. This implies that each *x*_{i}, if extended analytically to the complex *τ*-plane, will be finite-valued. One concludes that when *α*=*β*=*m*/(*m*+1), *m*≥1, the May–Leonard system has a weak form of the PP. This result appears to be new.

### Example 3.7

Consider the *α*=*β*=−1 May–Leonard system (with *a*_{*0}=0), i.e.
3.36
for which . For this system, the scaled state vector (*x*_{1},*x*_{2},*x*_{3}) is (2*y*_{1},2*y*_{2},*y*_{3}).

When *α*=−1, the solution *t*=*t*(*τ*) of ODE (3.34) can be expressed using the Weierstrassian elliptic function ℘(*τ*)=℘(*τ*;*g*_{2},*g*_{3}), which satisfies . In fact with these parameter values, gSE (3.6), also satisfied by *t*=*t*(*τ*), is an ODE listed by Carton-LeBrun [26], as her *n*=1 case XIII.4. Regardless of which ODE one integrates, one finds (with *C* a constant of integration) the general solution
3.37a
and
3.37b
where ℘ satisfies
3.38
On the right-hand side of (3.37a), *τ* can be replaced by *Aτ*+*B*. By substituting (3.37a) into (3.35) and using , one obtains the general solution *x*=*x*(*τ*) of the *α*=*β*=−1 system with *a*_{*0}=0; namely,
3.39
Again *τ* on the right-hand side can be replaced by *Aτ*+*B*, with *x* scaled by *A*; so there are three undetermined constants (*A*, *B* and *C*). A more symmetrical parametrization of the general solution uses the roots of the cubic polynomial in the ODE satisfied by ℘. In terms of *y*, it is
3.40a
where
3.40b
in which *e*_{1},*e*_{2},*e*_{3} satisfy *e*_{1}+*e*_{2}+*e*_{3}=0; they determine the trajectory. As (3.40a) is meromorphic on the complex *τ*-plane, the *α*=*β*=−1 May–Leonard system has the PP.

That the *α*=*β*=−1 system can be integrated using elliptic functions was noted by Brenig [59], but the explicit solution (3.40a) is new. Despite the mild similarity to example 3.4, trajectories *y*=*y*(*τ*) are easily seen to follow *non-planar* algebraic curves in . The quadratic functions *y*_{1}(*y*_{2}−*y*_{3}), *y*_{2}(*y*_{1}−*y*_{3}), *y*_{3}(*y*_{1}−*y*_{2}) are time-independent first integrals, equalling *e*_{2}−*e*_{3}, *e*_{3}−*e*_{1}, *e*_{1}−*e*_{2}. They are essentially the ones introduced as *I*_{1},*I*_{2},*I*_{3} in (1.9) (with *x* read as *y*). Starting from (3.40a), one can confirm that the additional first integrals recently found by Llibre & Valls [49] and Tudoran & Gîrban [50] are time-independent, such as the quartic one
3.41
This has the (constant) value . But the quadratic ones , defined by for each cyclic permutation *ijk* of 123, are perhaps the most fundamental of all. The reason is that, by examination, their respective values are 3*e*_{1},3*e*_{2},3*e*_{3}.

For the general solution *y*=*y*(*τ*) (or equivalently *x*=*x*(*τ*)) of the *α*=*β*=−1 May–Leonard system with a single non-zero intrinsic growth rate, see example 4.1.

## 4. Unequal growth rates and Painlevé transcendents

Section 3 treated *d*=3 Lotka–Volterra systems classified as ‘case *II*_{i}’ according to table 2, with zero intrinsic growth rates. (The extension to non-zero but equal rates is easy, as was noted.) Several *d*=3 systems with *unequal* intrinsic growth rates, which have the PP, will now be completely integrated. They are deformations of systems treated in §3 but do not fit into the framework of table 2. This is because their first integrals, though based on DPs, are not conserved: they depend exponentially on *τ*. Their solutions *x*=*x*(*τ*) turn out to involve Painlevé transcendents, which for Lotka–Volterra systems is a novelty.

Consider the system
4.1
which is the *n*=1 specialization of the case-*II*_{3} system (3.2), modified to include an intrinsic growth rate *λ* for species 3. Up to scaling, this is the generic three-species Lotka–Volterra system in which species 1,2 have zero growth rates, species 3 has equal effects on species 1,2, and species 3 has an equal but negated effect on itself. As in §§2 and 3, define an auxiliary variable *t* by *t*=(*x*_{1}+*x*_{2})/(*x*_{2}−*x*_{1}). By calculus, one deduces that
4.2
which previously appeared as equation (3.10b), in the case when *λ*=0.

The new time *t*=*t*(*τ*) satisfies a nonlinear third-order ODE which is an extended version of gSE (3.6), and which can be deduced in the same way. By tedious elimination of *x*_{1},*x*_{2},*x*_{3} from the rational expressions for in terms of *x*_{1},*x*_{2},*x*_{3}, one obtains
4.3
where gS(*a*_{1},*a*_{2};*b*_{1},*b*_{2};*n*) stands for the left-hand side of gSE (3.6). It is straightforward to rewrite the extended gSE (4.3) in the form , where
4.4
can be viewed as a time-dependent first integral. The existence of such a simple representation as is not surprising: by referring to (4.2), one can see that *J* is the familiar DP-based first integral (3.3) of §3, written in terms of (and with absolute value signs omitted, so that, in general, branch choices must be made). *J* is time-independent only if *λ*=0.

To deal with the now-present exponential dependence on *τ*, introduce a new independent variable *z*:=*e*^{λτ/c} for some as yet undetermined *c*≠0, and denote *d*/*dz* by a prime. The statement *J*=*J*_{0}*e*^{λτ} can be rewritten as a nonlinear second-order ODE for *t*=*t*(*z*), namely
4.5
where *K*=*c*^{2}*J*_{0}/*λ*^{2} is a constant of integration. Equation (4.5) is a ‘proto-Painlevé’ ODE, in the sense that, for certain *a*_{1},*a*_{2};*b*_{1},*b*_{2};*c*, it defines a Painlevé transcendent [60]. This is obscured by its singular points being , while the defining ODEs for the transcendents (traditional, due to Painlevé and Gambier) have singular points . But an equivalent ODE is
4.6
where *w*=(*t*+1)/(*t*−1)=*x*_{2}/*x*_{1} is a new, Möbius-transformed dependent variable, with *t*=(*w*+1)/(*w*−1); and, by definition, *a*_{3}:=1−*a*_{1}−*a*_{2} and *b*_{3}:=1−*b*_{1}−*b*_{2}.

### Example 4.1

For certain (*a*_{1},*a*_{2};*b*_{1},*b*_{2};*c*), ODE (4.6) for *w*=*w*(*z*) is or is reducible to a Painlevé-III (*P*_{III}) or Painlevé-V (*P*_{V}) equation. As traditionally formulated, the *P*_{III} and *P*_{V} equations have parameters *α*,*β*,*γ*,*δ*, and for *N*=*III* or *V*, *w*_{N}(*z*)=*w*_{N}(*α*,*β*,*γ*,*δ*;*z*) will denote any solution of the respective equation. By examination, the possibilities include
4.7
where *r*≠0 is arbitrary and *R*[*w*]:=−4*w*/(*w*−1)^{2}. In each case, *K*′∝*K* is a free parameter. These can be viewed as extensions to non-zero *λ* of Carton-LeBrun’s *n*=1 cases e.IV, XI.1, XIII.1, XIII.3, XIII.4. In each, by substituting *t*=(*w*+1)/(*w*−1) into (4.2), one obtains *x*=(*x*_{1},*x*_{2},*x*_{3}) in terms of the Painlevé transcendent *w*_{N}(*z*=*e*^{λτ/c}) and its derivatives. The resulting *x*=*x*(*τ*) is the general solution, as it includes three undetermined constants: besides *K*′∝*K*, the constants implicit in the choice of *w*_{III} or *w*_{V} from the two-dimensional solution space of *P*_{III} or *P*_{V}.

It should be noted that the last complete integration listed in (4.7), in terms of *w*_{V}, is of a deformed (*λ*≠0) version of the *α*=*β*=−1 May–Leonard system, which was integrated in example 3.7 in the case when *λ*=0, meaning that all intrinsic growth rates are zero. (That the *α*=*β*=−1 system has was mentioned in that example.)

The preceding results were stimulated by later work [61] of Carton-LeBrun, who studied and integrated many gSEs of an extended type that includes equation (4.3).

These results on *d*=3 Lotka–Volterra systems with unequal intrinsic growth rates are reminiscent of those of Bountis *et al*. [46], who reduced several *d*=2 systems with unequal rates to nonlinear second-order ODEs with the PP. But each of their ODEs could be integrated with the aid of elliptic functions, no Painlevé transcendents being needed.

## 5. Summary and final remarks

We have shown how to construct the general solution *x*=*x*(*τ*) of any *d*=3 Lotka–Volterra system in which species *j*,*k* are equally affected by species *i*, and the three species have equal intrinsic growth rates. Such systems had not previously been integrated, despite extensive work. The key result is equation (3.10), which gives the solution in parametric form. There *τ*,*x* are expressed with the aid of the incomplete beta function as functions of the new time *t*:=(*x*_{j}+*x*_{k})/(*x*_{k}−*x*_{j}), or equivalently of the ratio *w*:=*x*_{k}/*x*_{j}. As examples, certain ABC systems and fully symmetric (*α*=*β*) May–Leonard systems were integrated explicitly. If the integrated system has the PP, the new time is typically not needed: as a function of *τ*, the system state *x* can be expressed in terms of elementary or elliptic functions.

From the complete integration of any Lotka–Volterra system of this type, one can derive first integrals. One is the DP-based first integral of (1.4), the constancy of which defines case *II*_{i}. An additional first integral involving a quadrature was found in the ABC case by Goriely [36], and more generally by Gao [37, theorem 5]. But they did not construct general solutions. It seems more difficult to go from a pair of first integrals to a general solution, than the reverse.

The *α*=*β*=−1 May–Leonard model was treated in examples 3.7 and 4.1. If in this case the intrinsic growth rates of the three species are equal, *x*=*x*(*τ*) is expressed in terms of the Weierstrassian elliptic function ℘(*τ*); and if a single rate (*λ*) is non-zero, the fifth Painlevé transcendent *w*_{V}(*e*^{λτ}). The *α*=*β*=−1 model with equal growth rates was recently studied geometrically by Tudoran & Gîrban [50], but the elliptic general solution found here is new. The appearance of a Painlevé function is also novel.

The results of Carton-LeBrun [26,61] on nonlinear third-order ODEs with the PP, such as certain cases of the generalized Schwarzian equation (3.6) satisfied by the new time as a function of the old, proved invaluable. They are not well known, though they have been re-worked by Cosgrove [62, appendix C], as part of a programme to extend the Painlevé–Gambier classification of second-order ODEs with the property to ODEs of higher order. In several of the examples, they facilitated the construction of a non-parametric solution *x*=*x*(*τ*), which was expressed in terms of elementary or elliptic functions of *τ*. Many similar examples could be worked out. With further effort, the classification of Carton-LeBrun could be made to yield a classification of all case-*II*_{i} *d*=3 Lotka–Volterra systems with the PP.

It may be possible to extend the approach of this paper to treat integrable *d*=3 Lotka–Volterra systems of types not considered here, such as the ones of types I and III defined in table 2. Also, our techniques may be useful in constructing solutions of other small-*d* integrable quadratic differential systems. This includes small-*d* Lotka–Volterra systems of the many sorts mentioned in §1, and small-*d* quadratic systems not of the Lotka–Volterra form.

A related inverse problem suggests itself. Many functions of time *τ* have Lotka–Volterra *representations*, i.e. can be computed from solutions *x*=*x*(*τ*) of differential systems of Lotka–Volterra type. (This is stressed by Peschel & Mende [19], who even mention hardware implementations.) As §4 made clear, Painlevé functions are included. The extent to which solutions of higher order ODEs of Painlevé type, such as Chazy equations, can usefully be quadratically represented has not yet been resolved. This topic will be explored elsewhere.

## Appendix A. Incomplete beta functions and their inverses

In a normalization used here, the incomplete beta function *B*_{a1,a2;t0}(*t*), for (real) indices *a*_{1},*a*_{2}, (real) endpoints −1,1 and a (real) basepoint *t*_{0}≠−1,1, is defined by
A 1
The basepoint *t*_{0} lies in one of the intervals , (−1,1), , and by convention *τ*=*B*_{a1,a2;t0}(*t*) is defined only for *t* in that interval. It is an increasing function of *t*. It has been known since Liouville that, if any of *a*_{1},*a*_{2},*a*_{3}:=1−*a*_{1}−*a*_{2} is an integer, *B*_{a1,a2;t0}(*t*) is an elementary function. In particular, if *a*_{1},*a*_{2} are positive integers it is a polynomial; and if *a*_{1},*a*_{2} are integers with *a*_{1}*a*_{2}<0 and *a*_{1}+*a*_{2}≤0, it is a rational function. In general, it is a ‘special’ (higher transcendental) function. If *a*_{1},*a*_{2}>0, it can be expressed in terms of the traditionally normalized incomplete beta function, which is
A 2
for . Specifically,
A 3
with *C*=*C*_{a1,a2;t0} chosen so that *B*_{a1,a2;t0}(*t*_{0})=0. The traditional function is supported by many software packages, and many representations, expansions and approximations for it are known [63]. For instance,
A 4
where _{2}*F*_{1} is the Gauss hypergeometric function. Thus, for many choices of *a*_{1},*a*_{2}, closed-form expressions exist [64, §7.3]. A few have been derived in the context of hyper-logistic population growth [65].

The *inverse* incomplete beta function is also an increasing function of its argument, and is algebraic if *a*_{1},*a*_{2} are positive integers, or integers with *a*_{1}*a*_{2}<0 and *a*_{1}+*a*_{2}≤0. This inverse function is defined on some interval containing *τ*=0, and satisfies (with ) the hyper-logistic growth equation
A 5
with the initial condition . It maps the *τ*-interval onto whichever of , (−1,1), contains *t*_{0}. It is useful to write
A 6
where is any convenient, standardized solution of (A.5) that maps *some* *τ*-interval (*τ*_{min},*τ*_{max}) onto the *t*-interval containing *t*_{0}. The time-origin *τ*_{0} is determined by the condition that *B*^{−1}_{a1,a2}(−*τ*_{0})=*t*_{0}.

The function *τ*=*B*_{a1,a2;t0}(*t*) can be continued analytically from its real interval of definition to the upper half of the complex *t*-plane. (The continuation is given by (A.1) without absolute value signs, multiplied by a phase factor.) The continuation is a Schwarzian triangle function [54] that performs a Schwarz–Christoffel transformation: provided *a*_{1},*a*_{2},1−*a*_{1}−*a*_{2} are non-negative, it conformally maps the upper half of the complex *t*-plane to the interior of some triangle *ΔABC* in the *τ*-plane with vertex angles *π*(*a*_{1},*a*_{2},1−*a*_{1}−*a*_{2}), the points on the real *t*-axis being taken to the vertices *A*,*B*,*C* of the triangle. Many conformal mapping functions of this type, which are essentially incomplete beta functions, can be found in the catalogue of von Koppenfels & Stallmann [66]. For each, the inverse takes the interior of *ΔABC* to the upper-half *t*-plane.

This inverse can sometimes be given in closed form; for example, when the unordered set {1/*a*_{1},1/*a*_{2},1/(1−*a*_{1}−*a*_{2})} is any of {2,4,4}, {2,3,6} or {3,3,3}. In these cases, the inverse can be continued to the entire complex *τ*-plane from its real interval of definition, and furthermore from *ΔABC*, by applying the Schwarz reflection principle: reflecting repeatedly through the sides of the triangle. The resulting one-valued function of *τ* is elliptic (i.e. doubly periodic) and can be given explicitly. The cases {1,*m*,−*m*} (for *m* any positive integer), and are similar, but yield inverse functions *t*=*B*^{−1}_{a1,a2;t0}(*τ*) that are elementary rather than elliptic.

Thus, there are many ‘conformally distinguished’ choices for the unordered set {1/*a*_{1},1/*a*_{2}}, each of which yields an explicit formula for the inverse function . The elementary cases are {1,*m*}, {1,−*m*}, {*m*,−*m*}, , {2,2}, , , and the elliptic ones are {2,3}, {2,4}, {2,6}, {3,3}, {3,6}, {4,4}. In each of these cases, any standardized inverse function of the type used in (A.6), which has a real interval of definition (*τ*_{min},*τ*_{max}), can optionally be continued to a one-valued function on the *τ*-plane.

Table 3, which is an integral part of this appendix, gives such a standardized inverse function for . These are the conformally distinguished cases with *a*_{1}=*a*_{2}=:*a*, and are of particular interest (the many ones with *a*_{1}≠*a*_{2} are left to the reader). In each case, there are two subcases: *t*_{0}∈(−1,1) and . The inverse function *t*=*t*(*τ*) maps *τ*∈(*τ*_{min},*τ*_{max}) onto *t*∈(−1,1), respectively, ; and by (A.5) it satisfies
A 7
In the table, each *t*=*t*(*τ*) has been chosen to satisfy *t*(*τ*=0)=0, respectively, *t*(*τ*=0)=±1 (except when *a*=0). Each is increasing on its interval of definition (*τ*_{min},*τ*_{max}).

The elementary inverses listed in table 3 (for ) are obtained by inspection, i.e. by integrating ODEs (A.7). The elliptic solutions of these ODEs (for ) are less obvious, but can be derived with the aid of function identities from elliptic functions that are known to map certain triangles conformally to the upper half-plane [67,68]. For , the Jacobian elliptic functions sn,cn,*dn* that appear in the table are ‘lemniscatic’: their common modular parameter *m*=*k*^{2} equals . For , the Weierstrassian elliptic function ℘(*τ*)=℘(*τ*;*g*_{2},*g*_{3}) that appears is ‘equianharmonic’: its parameter *g*_{2} equals zero. Also, its parameter *g*_{3} equals . It satisfies the usual Weierstrassian ODE , i.e. .

For each function appearing in the table, the corresponding interval of definition *τ*∈(*τ*_{min},*τ*_{max}) is easily determined from (A.1). It can be expressed in terms of the traditionally normalized *complete* beta function *B*(*a*_{1},*a*_{2}):=*Γ*(*a*_{1})*Γ*(*a*_{2})/*Γ*(*a*_{1}+*a*_{2}). In the subcases *t*_{0}∈(−1,1) and , the interval of definition *τ*∈(*τ*_{min},*τ*_{max}) is
A 8
where
It is understood that is infinite if *a*≤0. Thus, are , as one sees in the table. Also, in the subcases the interval ±*τ*∈(0,*K*_{a}) must be interpreted as the interval if , and replaced by if *a*≤0.

The values of *a* in the table are not the only ones for which the inverse function can be expressed in closed form. It is an algebraic function if *a* is any positive integer, and if *a*=2 it can even be expressed in terms of radicals, by employing Cardan’s formula. But for any integer *a*>1, its continuation to the *τ*-plane is multiple-valued.

No table of closed-form expressions for the inverse incomplete beta function, resembling table 3, seems to have appeared previously. However, certain numerical values of the traditionally normalized version of this function (especially, for large (half-)integral values of *a*_{1},*a*_{2}) have been calculated algorithmically and tabulated, owing to their importance in statistics [69–71].

- Received November 24, 2012.
- Accepted July 23, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.