Elliptic Hypergeometric Laurent Biorthogonal Polynomials with a Dense Point Spectrum on the Unit Circle

Using the technique of the elliptic Frobenius determinant, we construct new elliptic solutions of the $QD$-algorithm. These solutions can be interpreted as elliptic solutions of the discrete-time Toda chain as well. As a by-product, we obtain new explicit orthogonal and biorthogonal polynomials in terms of the elliptic hypergeometric function ${_3}E_2(z)$. Their recurrence coefficients are expressed in terms of the elliptic functions. In the degenerate case we obtain the Krall-Jacobi polynomials and their biorthogonal analogs.

In the degenerated case (when both periods of elliptic functions become infinity) we obtain biorthogonal analogs of the Krall-Jacobi orthogonal polynomials. We show that these biorthogonal polynomials satisfy a 4th order differential equation which can be presented in the form of quadratic operator pencil.

Laurent biorthogonal polynomials and their basic properties
The Laurent biorthogonal polynomials LBP P n (z) appeared in problems connected with the two-points Padé approximations (see, e.g., [15]).
Let L be some linear functional defined on all possible monomials z n by the moments c n = L{z n }, n = 0, ±1, ±2 . . . .
In general the moments c n are arbitrary complex numbers. The functional L is defined on the space of Laurent polynomials P(z) = N 2 n=−N 1 a n z n where a n are arbitrary complex numbers and N 1,2 arbitrary integers: n=−N 1 a n c n .

7)
b n = d n h n h n−1 = T n+1 ∆ n−1 T n ∆ n = 0, n = 1, 2, . . . from which one can obtain expression for the normalization constant h n in terms of the recurrence parameters: There is a one-to-one correspondence between the moments c n and the recurrence coefficients b n , d n (provided restrictions b n d n = 0 are fulfilled).
We say that the LBP are regular if b n d n = 0. This condition is equivalent to the condition ∆ n ∆ (1) n = 0, n = 0, 1, . . . .
In the regular case there is a simple formula relating the biorthogonal partners Q n (z) with polynomials P n (z) [12]: Q n (z) = zP n+1 (1/z) − z n−1 P n (1/z) P n (0) .
The rescaled LBPP n (z) differ from initial LBP P n (z) only by a trivial rescaling of recurrence parameters. The momentsc n of the rescaled LBP are connected with initial moments c n by the relationc n = q n c n . Note that the rescaled biorthogonal partners Q n (z) are transformed as Q n (z) = q −n Q n (zq). (2.11) There is a connection between the LBP and the restricted relativistic Toda chain [16]. Assume that LBP P n (z; t) depend on an additional (so-called "time") parameter t. This mean that the recurrence coefficients b n (t), d n (t) become functions of the parameter t. We assume that the relatioṅ holds for all n = 0, 1, . . . . This ansatz leads to the following equations for the recurrence coefficients [16] For the corresponding moments c n (t) we have the relatioṅ c n = c n−1 , n = 0, ±1, ±2, . . . .
Another possible ansatz [16] P n (z) = −b n (P n (z) − zP n−1 (z)) leads to the equationṡ In this case we have for the moments the relatioṅ In spite of the apparent difference between equations (2.12) and (2.13), it can be shown (see, e.g., [16]) that these two systems are both equivalent to the restricted relativistic Toda chain equations. The term "restricted" in this context means that it is assumed an additional condition b 0 = 0. This means that in formulas (2.12) or (2.13) we should assume n = 0, 1, 2, . . . . For the nonrestricted relativistic Toda chain equations (2.12) or (2.13) are valid for all integer values of n = 0, ±1, ±2, . . . .

Laurent biorthogonal polynomials and QD-algorithm
The (restricted) "discrete-time" relativistic Toda chain corresponds to the following ansatz for the moments where h is an arbitrary parameter. We have the transformation formula for the corresponding Laurent biorthogonal polynomials Formulas (3.1) and (3.2) can be interpreted as Christoffel and Geronimus transformations for LBP [32]. The corresponding recurrence coefficients are transformed as [32] 3) it is assumed the argument t for the coefficients b n , d n ). These relations can be presented in a slightly different equivalent form as where we have denotedb n = b n (t + h) etc for brevity. Relations (3.4) describe so-called QDalgorithm for the two-point Padé approximation (see, e.g., [5] for details). In other words, the (restricted) discrete-time relativistic Toda chain is equivalent to the QD-algorithm for the two-point Padé approximation. Usually, this algorithm works as follows. We start from the given moments c n (t), n = 0, ±1, ±2 . . . where the dependence on "time" is trivial: and define the coefficient d 0 (t) for all t = t 0 + jh, j = 0, ±1, ±2 as The initial value t 0 is not essential, usually it is assumed that t 0 = 0, in this case we can write Assume that b 0 (t) = 0 for all t. Then at the first step we find b 1 (t) = b Then we find d n = 0 for all n and j. Then we obtain all sequences d n , n = 0, 1, 2, . . . for j = 0, ±1, ±2, . . . . There is a remarkable connection with the QD-algorithm for the ordinary orthogonal polynomials [5]. Indeed, let us introduce the monic polynomials (3.5) where the polynomials P (j) n (z) are defined as P (j) n (z) = P n (z; hj). Then relations (3.1) and (3.2) become and These relations can be interpreted as Geronimus and Christoffel transforms for the orthogonal polynomials W (j) n (z). The compatibility condition between (3.6) and (3.7) leads to the recurrence relation which describes the three-term recurrence relation for the ordinary orthogonal polynomials W (j) n (z) where the recurrence coefficients are [5] g (j) Moreover we have compatibility conditions for the coefficients e (j) Relations (3.8) coincide with those introduced by Rutishauser and describing the ordinary QDalgorithm [7]. It is easy to verify that relations (3.8) are equivalent to relations (3.4) for the two-point QD-algorithm. Thus starting from known solution P (j) n of the discrete-time relativistic Toda chain (or, equivalently, two-point QD-algorithm) we can obtain a set of the ordinary orthogonal polynomials W (j) n (z) depending on additional "time" parameter j. Note that sometimes the introduced orthogonal polynomials W Thus the orthogonal polynomials W where the linear functional τ (j) is defined by the moments

The normalization constant q (j)
n has the expression It would be instructive to interpret (3.1) and (3.2) in terms of so-called bilinear technique by using the determinantal identities. This technique is standard in the theory of integrable systems.
As a first step, we give a compressed expression to d n − b n as which can be derived from the determinantal identity, or Jacobi identity, for the Toeplitz determinant: Then the relations (3.1) and (3.2) can be transformed to the following bilinear equations, respectively, where the functions σ (j) n are defined by n is proportional to the Laurent biorthogonal polynomial P

Laurent and Baxter biorthogonal polynomials
There is an alternative (but essentially equivalent) approach to biorthogonal polynomials proposed by G. Baxter [4]. The pair P n (z), Q n (z) of the biorthogonal polynomials is defined in this approach by means of initial conditions P 0 = Q 0 = 1 and the following recurrence system where e (1,2) n are some complex coefficients. It is clear that e (1) n = −Q n+1 (0). Notation P * n (z) is standard for so-called reciprocal polynomials, i.e. P * n (z) = z n P n (1/z), Q * n (z) = z n Q n (1/z). Assume that e (1) n e (2) n (1 − e (1) n e (2) n ) = 0 (this is the nondegenerate case). Then, excluding Q * n (z) from the system (4.1) we arrive at a 3-term recurrence relation for the polynomials P n (z): n−1 ).
Conversely, assume that we have the nondegenerate Laurent biorthogonal polynomials P n (z) satisfying (2.6). We can construct their biorthogonal partners Q n (z) by (2.10). Then it is elementary to verify that polynomials P n (z), Q n (z) satisfy system (4.1) with e (1) n = −Q n+1 (0). Sometimes system (4.1) is more convenient for analysis due to apparent symmetry between polynomials P n (z), Q n (z) and corresponding coefficients e (1) n , e (2) n . Note also that the Laurent and Baxter biorthogonal polynomials in turn are equivalent to the socalled Laurent orthogonal polynomials proposed by Jones and Thron [15]. The Jones and Thron polynomials contains terms z k with both positive and negative degree k. For details of this equivalence see, e.g., [12] and [22].
There is an important special case when all the Toeplitz determinants are positive ∆ n > 0 and moreover the moments satisfy the condition c n = c −n (as usual,c n means complex conjugation of c n ). In this case the biorthogonal partners Q n (z) coincide with complex conjugated polynomials Q n (z) =P n (z) and there exists nondecreasing function σ(θ) of bounded variation on the unit circle such that the orthogonality relation holds. I.e. in this case we have polynomials P n (z) which are orthogonal on the unit circle (abbreviated as OPUC [24]). Historically, these polynomials were introduced first by Szegő [29] and are called the Szegő polynomials orthogonal on the unit circle. They satisfy the recurrence relation P n+1 (z) = zP n (z) − a n z nP n (1/z), where the coefficients a n = −P n+1 (0) are called the reflection (or Schur, or Verblunsky, . . . ) parameters. The relation (4.3) was first derived by Szegő himself [29]. The reflection parameters are complex numbers satisfying the important inequality In fact, condition (4.4) is equivalent to the condition of positive definite Toeplitz forms ∆ n > 0 or to existence of a positive measure on the unit circle providing orthogonality property (4.2). If, additionally, all the moments are real, then they satisfy condition c −n = c n . In this case the reflection parameters are real parameters satisfying the restriction −1 < a n < 1, n = 0, 1, 2, . . . . The biorthogonal partners then coincide with initial polynomials Q n (z) = P n (z). It is easy to show that the measure dσ is symmetric with respect to real axis in this case, namely the function σ(θ) satisfies the condition σ(2π − θ) + σ(θ) = const.
Recall that the Weierstrass sigma function is defined by the infinite product [1] where the product is taken over all points of the lattice s = 2mω 1 +2m ′ ω 3 , m, m ′ = 0, ±1, ±2, . . . excluding the point with m = m ′ = 0. 2ω 1 and 2ω 3 are the so-called primitive elliptic periods. It is convenient to introduce the third period 2ω 2 = −2ω 1 − 2ω 3 [1]. The Weierstrass sigma function possess quasi-periodic properties [1] where the constants η α are defined as We have v i for simplicity, moreover it is assumed that the upper limit for i, j in the products is n − 1). Formula (5.1) was obtained by Frobenius in [9]. A simple elementary method to derive formula (5.1) can be found in [3]. Frobenius and Stickelberger derived also in [8] several other explicit formulas for "elliptic determinants" in connection with the theory of rational interpolation.
Let φ k (x), ψ k (x), k = 0, 1, . . . (we assume that φ 0 = ψ 0 = 1) be two sets of functions in some argument x. Assume that there exists a linear functional L such that The linear functional L is defined on the space of functions constructed from bilinear combinations of the type with arbitrary coefficients c ik .
Introduce the following functions By construction, these functions are biorthogonal with respect to the functional L, where the normalization coefficients h n are Expanding the determinant in (5.2) over the last row we have explicit expression for the polynomial P n (x): The auxiliary determinants H n (k) are defined by canceling the kth column, i.e. Here are "generalized binomial coefficients". Similar expression can be obtained for the biorthogonal partners Q n (x) if one replaces the parameters v i with u i . In case when the sequence v j is linear with respect to j: v j = wj+ξ we obtain the conventional "elliptic binomial coefficients" [10]: is elliptic Pochhammer symbol. Note that usually the elliptic number is defined in terms of the theta function [x] = θ 1 (wx)/θ 1 (w) [10], but for our purposes these definitions are in fact equivalent.
We thus constructed an explicit system of biorthogonal functions P n (x), Q n (x) starting from the elliptic Frobenius determinant. This system can be further specified by a concrete choice of the basic functions φ n (x), ψ n (x) and the linear functional σ. Note that the idea to construct explicit families of biorthogonal functions directly from corresponding Gram determinants is due to Wilson [31]. For general biorthogonal rational functions the determinant representation can be found e.g. in [25] and [6]. Put where w is an arbitrary real parameter which is incommensurable with the real period 2ω 1 over the integers, i.e. we will assume that for any integers N 1 , N 2 . Then for the entries of the Frobenius matrix we have This matrix has the Toeplitz form. We can therefore define corresponding monic Laurent biorthogonal polynomials by the formula where the moments are defined as and the Toeplitz determinant ∆ n is defined by (5.3). As in the previous section, define the elliptic numbers [x] as and the elliptic Pochhammer symbol The elliptic hypergeometric function is defined by the formula We have Proposition 1. The Laurent biorthogonal polynomials defined by formulas (6.2) and (6.3) are expressed in terms of the elliptic hypergeometric function: is the coefficient to provide monicity P n (z) = z n + O(z n−1 ) of the polynomials P n (z).
Remark. The parameters of the elliptic hypergeometric function in our case satisfy condition and hence M = 0 in the definition of the hypergeometric function (6.4). Our definition of the elliptic hypergeometric function is in accordance with the conventional one [10,27]. The main difference is replacing the theta functions with the Weierstrass sigma functions. This replacement leads to appearance of the additional factor e M s(s−1) . Indeed, there is relation between these functions [1] (the constant factor is not essential because it is canceled in all expressions for elliptic hypergeometric series). Using this relation we can replace all sigma functions with the theta functions θ 1 (z) which leads to formula (6.4).
Now we calculate the normalization coefficients h n directly from Frobenius formula (5.1): In what follows we will assume the following restriction α = wm for any integers m. Indeed, otherwise the normalization coefficient h n becomes singular and we have a degeneration. We observe also that the determinant ∆ Formulas (6.7) and (6.8) allow us to find explicit expressions for the recurrence coefficients b n , d n .
Indeed, from (2.7) and (2.8) we have We thus obtained a new explicit example of the Laurent biorthogonal polynomials which have both explicit expression in terms of the elliptic hypergeometric function (6.5) and explicit recurrence coefficients (6.9), (6.10). As a by-product, we have also obtained a new explicit solution of the discrete-time relativistic Toda chain or, equivalently, a new explicit solution of the two-point QD-algorithm. Indeed, the recurrence coefficients b n , d n given by (6.9) and (6.10) provide an explicit elliptic solution of the two-point QD-algorithm (3.4) with t = α, h = w. In turn, using correspondence (3.5) we can obtain elliptic solution of the ordinary QD-algorithm (3.8), or equivalently, the discrete-time Toda chain solutions. As far as we know these solutions are new.
In order to find explicit (bi)orthogonality relation for these polynomials we need first the explicit Fourier expansion of the elliptic functions of the second kind. We will do this in the next section.
7 Fourier series of the elliptic functions of the second kind Assume that f (z) is the simplest elliptic function of the second kind [1] f (z) = κ σ(z + α + β) σ(z + α) e γz (7.1) with some complex parameters κ, β, α, γ. The function f (z) is quasi-periodic with respect to periods 2ω 1 , 2ω 3 : where µ 1 = e 2η 1 β+2ω 1 γ , µ 3 = e 2η 3 β+2ω 3 γ . We demand that function f (z) be purely periodic with respect to the (real) period 2ω 1 j: where j = 1, 2, . . . is an arbitrary positive integer. This leads to the condition µ j 1 = 1 or where m = 0, ±1, ±2, . . . . Note that for j > 1 we should avoid the values m = 0, ±j, ±2j, . . . because they correspond to pure 2ω 1 -periodicity. Of course, it is assumed that m and j are coprime, i.e. µ 1 is a primitive root of the unity of order j: Moreover, we assume that α = −α 0 − iα 1 , where both parameters α 0,1 are real and are restricted by conditions Conditions (7.4) mean that the parameter −α lies within the fundamental parallelogram (i.e. rectangle in our case). If α takes values beyond this parallelogram, it is possible to reduce it to canonical choice (7.4) using shifts by periods 2ω 1 , 2ω 3 . Due to quasiperiodicity property of the function f (z) this will lead only to redefining of the parameter γ. Moreover we assume that the imaginary part −α 1 of α is nonzero. This assumption is very natural if we would like to avoid singularities of the function f (z) on whole real axis. Equivalently, one can present α in the form where 0 < ν < 1 is a fixed parameter which describes the relative value of the imaginary part α 1 = −2iνω 3 with respect to the imaginary period 2ω 3 . Thus we have the function f (z) which is periodic and bounded on the whole real axis. It is possible therefore to present f (z) in terms of the Fourier series Our problem now is to calculate the Fourier coefficients A n . By definition, (the integral is well defined because by our assumptions the function f (z) has no singularities on the real axis). In order to calculate the integral in (7.7) we exploit standard method of contour integration (see, e.g., [1] for calculation of the Fourier expansion for Jacobi elliptic functions). Choose the contour Γ as the rectangle with vertices (0, 2jω 1 , 2jω 1 + 2ω 3 , 2ω 3 ) (i.e. the horizontal length is 2jω 1 and vertical length 2|ω 3 |).
We have (the contour is traversed counterclockwise) where 1 , 3 correspond to horizontal sides of the rectangle, and integrals 2 , 4 correspond to vertical sides.
Due to periodicity property f (z + 2jω 1 ) = f (z) we have 2 + 4 = 0. For the two remaining horizontal integrals we have Making the shift z → z + 2ω 3 and using quasi-periodic property (7.2) we have and thus Hence, in order to calculate the Fourier coefficient A n we need to calculate the contour integral in l.h.s. of (7.8). This can be done by standard methods of residue theory. Indeed, inside the contour Γ the function f (z) exp −2πinz T has only j simple poles located at points At z s the function f (z) has the residue r s = µ s 1 r.
Hence we have that the residue R n of the function f (z)/T exp az T inside the rectangle Γ will be If n = m mod j then R n = 0. Nonzero value of the residue will be only for n = m + jt, t = 0, ±1, ±2, . . . . In this case Comparing with (7.8) we get , n = m, m ± j, m ± 2j, . . .

and
A n = 0, if n = m mod j.
Recall that j is a fixed positive integer -the order of the root of unity µ 1 , while m is a fixed nonnegative integer (lesser than j) coprime with j. Thus for large j the nonzero coefficients A n are more rare then for small j.
There are two important simplest cases: (i) if j = 1 and m = 0. This case corresponds to the period 2ω 1 . Then the Fourier coefficients A n are nonzero for all n = 0, ±1, ±2, . . . and we have (ii) if j = 2 and m = 1. This case corresponds to the period 4ω 1 . In this case all even Fourier coefficients are zero A 2n = 0 and for the odd Fourier coefficients we have Note that in all cases the Fourier series (7.6) converges inside the strip − This results follows from standard theorems concerning asymptotic behavior of the Fourier coefficients A n and A −n for n → ∞ [1]. The parameters v 1 , v 2 are positive as follows from the inequality 0 < ν < 1. These conditions are very natural because the boundary lines Im(z) = 2ν|ω 3 | and Im(z) = 2(ν − 1)|ω 3 | of the strip pass through the poles of the function φ(z). Note that for ν = 1/2 (i.e. when the pole of the function φ(z) lies on the horizontal line Im(z) = |ω 3 |) we have the strip symmetric with respect to the real line: |Im(z)| < |ω 3 |. The latter case correspond, e.g., to the Jacobi elliptic functions sn(z; k), cn(z; k), dn(z; k) [1].

Explicit biorthogonality relation
In this section we obtain explicit biorthogonality property of the obtained Laurent biorthogonal polynomials.
To do this we need to find explicit realization of the moments c n given by formula (6.3). We note that where f (z) is the elliptic function of the second kind (7.1) (in our case κ = 1/σ(β) but the constant κ does not play any role in formulas for the polynomials P n (z) and their recurrence coefficients b n , d n ).
is an infinite set of points belonging to the unit circle |z s | = 1. Due to condition (6.1) we have that all these points are distinct z s = z t if t = s and hence they are dense on the unit circle. Thus we found explicit realization of the moments c n and hence we immediately obtain biorthogonality relation for our Laurent biorthogonal polynomials where Q n (z) are biorthogonal partners (2.2) with respect to polynomials P n (z). The Fourier coefficients A s play the role of discrete weights in this biorthogonality relation. Hence we have obtained In the periodic case f (z + 2ω 1 j) = f (z) the elliptic polynomials (6.5) P n (z) are biorthogonal (8.3) on the unit circle |z| = 1 with respect to a dense point measure with weights A s given by expression (7.9).
Note that the biorthogonal partners Q n (z) in our case can be found explicitly in terms of the elliptic hypergeometric function. Indeed, from (2.2) we see that the polynomials Q n (z) are Laurent biorthogonal polynomials corresponding to the "reflected" momentsc n = c −n . From explicit expression (6.3) it follows that the moments c −n are obtained from the moments c n by reflection of the parameters α → −α, β → −β, γ → −γ, whereas the parameter w remains unchanged (under such procedure we obtain the moments −c −n but any constant common factor in front of moments leads to the same polynomials Q n (z)). Hence we can obtain expression for the polynomials Q n (z) from the expression (6.5) for polynomials P n (z) by reflection of parameters α, β, γ: where the coefficientB n is obtained from corresponding coefficient B n (6.6) by the same reflection of the parameters α, β, γ. Thus both polynomials P n (z) and their biorthogonal partners Q n (z) have similar expressions in terms of elliptic hypergeometric function.
So far, we assumed that the function f (z) is periodic with the period 2ω 1 j. This assumption means that the parameter γ should satisfy condition (7.3). Parameters α and β are assumed to be arbitrary (with the only condition (7.4)). What happens if the function f (z) is not periodic, i.e. if the parameter γ is arbitrary? It appears that this general case can be easily reduced to the already considered. Indeed, assume that we change the parameter γ, i.e. assume that the parameters α and β remain the same butγ = γ + χ, where χ is an arbitrary complex parameter. Then it is easily seen from explicit expression (6.5) that the new Laurent biorthogonal polynomialsP n (z) are obtained by simple rescaling of the argument: P n (z) = q n P n (z/q), where q = e wχ . This corresponds to transformation of the momentsc n = ǫq n c n as seen directly from (6.3) (the common constant ǫ = e αw is inessential and can be put equal to 1).
Assume that we choose the parameter χ such that the new functioñ f (z) = κ σ(z + α + β) σ(z + α) eγ z will be periodic with the period 2ω 1 j. This means that the parameter χ should be chosen from condition (see (7.3)) where m is co-prime with j.
Then the new polynomialsP n (z) will be biorthogonal on the unit circle according to above obtained proposition: where the spectral points z s on the unit circle are given by (8.2) and the weights A s by (7.9). Note that the normalization coefficients h n remain unchanged under the rescaling transform as seen from (2.9), i.e.h n = h n . Taking into account thatQ n (z) = q −n Q n (z) (see (2.11)) we obtain from (8.5) the biorthogonal relation Relation (8.6) means that for generic values of γ polynomials P n (z) and Q n (z) are biorthogonal on the non-unit circle |z| = 1/|q| with respect to the same dense point measure.
It is interesting to note that for every integer j = 1, 2, . . . (i.e. for every period T = 2ω 1 j) we can construct corresponding circle providing biorthogonality relation (8.6). Thus there exist infinitely many orthogonality circles for different values of the integer parameter j.
For the radius r of the circle of biorthogonality we have from (8.4) (recall that we assume parameter w to be real) r = 1/|q| = e η 1 βw ω 1 |e wγ | .

Positivity of the measure and polynomials orthogonal on the unit circle
Return to the case when the function f (z) is periodic with the period 2ω 1 j and consider an important special case when all the Fourier coefficients of the function f (z) are nonnegative A n ≥ 0. In this case all spectral points z s belong to the unit circle |z s | = 1 and the measure on the unit circle is a positive nondecreasing function. We have 0 < h < 1. Thus for n → −∞ we have It is seen that for positivity of A n one should have 2πiR 0 = κ 0 , where κ 0 is a positive parameter, and for the real part of α we have the condition Now for for n → ∞ we have In this case we should have necessarily It is easily seen that conditions (9.1) and (9.2) are also sufficient and so we have the Proposition 3. The Fourier coefficients are positive (up to inessential common factor) if and only if the real parts of parameters α, β satisfy conditions (9.1) and (9.2). In this case the expression for the Fourier coefficients can be presented in the form is a positive parameter (as usual by Im(β) we denote the imaginary part of β).
In this case we have positive dense point measure on the unit circle. It is well known that when the measure dσ is positive on the unit circle then biorthogonal polynomials become the orthogonal polynomials on the unit circle [29,11,24]. In this case the moments c n satisfy the restriction c −n =c n and moreover all the Toeplitz determinants are positive ∆ n > 0, n = 1, 2, . . . .
The property c −n =c n can be verified directly from the definition (6.3) if the parameters α, β satisfy conditions: (here β 1 is an arbitrary real parameter). In this special case the obtained polynomials satisfy the Szegő recurrence relation (4.3). The reflection parameters a n are calculated as a n = −P n+1 (0) and using already found explicit formula (6.5) for polynomials P n (z) we have a n = −B n+1 , where B n is given by (6.6) (with α, β satisfying restrictions (9.4)). From general theory it follows that in this case the reflection parameters should satisfy the restriction |a n | < 1. This property is not obvious from explicit expression for a n in terms of elliptic Pochhammer symbols.
If, in addition to positivity of A n , we demand that the discrete measure should be symmetric with respect to the real axis we then obtain the condition A −n = A n for all n = 0, 1, 2, . . . . It is easily verified from explicit expression (9.3) that this is possible only for j = 1 and j = 2. In the first case, when j = 1 the period T = 2ω 1 and necessarily ν = 1/2 and κ 1 = 1, so that (9.5) But the Fourier coefficients with expression (9.5) correspond to the Jacobi elliptic function dn(z; k) [30]. In this case the moments are c n = dn(wn; k) and indeed satisfy the property c −n = c n ; the reflection parameters are very simple: a n = dn(w(n + 1); k) for the even n and a n = −cn(w(n + 1); k) for the odd n.
In the second case, i.e. when j = 2 we have the period T = 4ω 1 and necessarily ν = 1/2 and κ 1 = h −1 , so that These Fourier coefficients correspond to the Jacobi elliptic function cn(z; k) [30]. Again the moments c n = cn(wn; k) satisfy the desired property c −n = c n and we have the polynomials orthogonal on the unit circle with simple reflection parameters: a n = cn(w(n + 1); k) for the even n and a n = −dn(w(n + 1); k) for the odd n.
These two explicit cases of OPUC with dense point measure were first considered in [33]. Now we see that there exists much wider class of explicit elliptic OPUC with positive dense measure on the unit circle. This class of OPUC contains essentially 3 arbitrary continuous parameters: w, Im(α) = −α 1 , Im(β). We thus have two additional parameters with respect to the only parameter w in [33]. Note however, that if one demands that OPUC were real (i.e. they have real reflection parameters a n and moments c n ) then nothing more general than "cn-" and "dn-"polynomials considered in [33] appear.

"Classical" property of LBP
Assume that P n (z) are arbitrary Laurent polynomials satisfying 3-term RR P n+1 (z) + d n P n (z) = z(P n (z) + b n P n−1 (z)) with some recurrence coefficients b n , d n .
For any sequence µ n , n = 0, 1, 2, . . . of complex numbers such that µ 0 = 0 we define the linear operator D which acts on the space of polynomials in the argument z by the rule Dz n = µ n z n−1 .
Then it is clear that the operator D sends any polynomial of degree n to a polynomial of degree n − 1 and moreover D{1} = 0. In this sense the operator D can be called as a generalized derivative operator. If µ n = n then D = ∂ z coincides with the ordinary derivative operator with respect to the variable z.
In [32] we considered the case of the ordinary classical LBP (i.e. with respect to the operator ∂ z ) and derived necessary and sufficient conditions for existence of such polynomials. It appears that there exists many different types of such classical LBP. The simplest ones are the LBP constructed by Hendriksen and van Rossum [12]. The latter have explicit expression in terms of the Gauss hypergeometric function. Now return to our elliptic LBP P n (z) = P n (z; α, β, γ, w) (we indicate dependence on parameters α, β, γ, w for convenience) and consider the operator D with µ n defined as .
Then from explicit representation (6.5) it is elementary to verify that the operator D transforms these Laurent biorthogonal polynomials to the same family but with the sole parameter β changed: DP n (z; α, β, γ, w) = µ n P n−1 (z; α, β − α, γ, w) which means that our elliptic polynomials P n (z; α, β, γ, w) are indeed "classical" polynomials with respect to the operator D.
One can consider the "µ-exponential" function E µ (x) which is a formal solution of the operator equation Clearly we have a solution of the operator equation (10.3) in terms of the formal series In case of the elliptic µ n (10.2) we have Obviously for an arbitrary complex parameter γ we have We can also introduce "even" and "odd" µ-exponential functions which are µ-analogs of the hyperbolic "cosh" and "sinh" functions We have obvious relations These functions both have the same property with an arbitrary parameter γ, and hence for arbitrary parameters β 0 , β 1 the function f (x) = β 0 C µ (γx) + β 1 S µ (γx) is a formal solution of the operator equation Note the obvious "intertwining" property of these functions: We can use these properties in order to construct formal generating functions for the cn and dn-circle polynomials. Indeed, let P (C) n (z) and P (D) n be cn and dn-circle polynomials corresponding to the choices α = −ω 3 , β = −ω 2 and α = −ω 3 , β = ω 1 respectively. As shown in [33] these polynomials satisfy intertwining properties and hence where one can choose µ n = sn(wn)/sn(w). Construct the generating functions for the polynomials P n (z) as the formal series We have obviously where notation D z means that the operator D acts only on the variable z. As a consequence This property means that both functions F (C) (z; t) and F (D) (z; t) can be expressed in terms of "even" and "odd" µ-exponential functions with respect to the variable z: with some functions ξ i (t), η i (t), i = 1, 2. For these functions we have relations η 0 (t) = ξ 1 (t), η 1 (t) = ξ 0 (t) which follow easily from intertwining relations. Remaining functions ξ 0 (t), ξ 1 (t) can be found as follows. Put z = 0. Then from definition (10.5) we have where B We thus have explicit expressions for the functions ξ 0 (t), ξ 1 (t) in terms of formal series As shown in [33] the coefficients B are corresponding reflection parameters. We can thus present expressions for ξ 0 (t) and ξ 1 (t) in a more explicit form: dn(w(2s + 1))t 2s+1 sn 2s+1 (w) sn(w)sn(2w) · · · sn(w(2s + 1)) , cn(w(2s + 1))t 2s+1 sn 2s+1 (w) sn(w)sn(2w) · · · sn(w(2s + 1)) .
Hence we were able to find explicitly the generating functions for the OPUC P

Rational limit of the elliptic functions and corresponding Laurent biorthogonal polynomials
In this section we consider the rational limit of the elliptic functions, when both periods 2ω 1 , 2ω 2 tend to infinity. In this case for the Weierstrass functions we have simple formulas [1] σ(z) = z, ζ(z) = 1/z, ℘(z) = 1/z 2 .
The biorthogonal partners have the expression In order to find the orthogonality measure for these polynomials we first note that the moments (11.1) can be rewritten in the form where the moments c (0) n = 1/(n + α) = lim β→∞ c n correspond to a special case of the "classical" Laurent biorthogonal polynomials considered by Hendriksen and Van Rossum [12]. The moments c (0) n correspond to the recurrence coefficients .
The moments (11.5) differ from the moments c (0) n only by the constant term β −1 . This means that corresponding orthogonality measure on the unit circle has an additional concentrated mass at z = 1. Thus the biorthogonality relation for polynomials (11.2) looks as where the weight function is Note that inserting the concentrated mass at point z = 1 can be performed by the Geronimus transform of the classical Hendriksen-van Rossum polynomials (see [32] for details).
In more details DP n (z; α, β) = µ n P n−1 (z; α, β − α), where P n (z; α, β) stands for LBP (11.2) with explicit dependence on the parameters α, β. Note the operator D in this case does not coincide with the ordinary derivative operator ∂ z . Hence the polynomials (11.2) provide one of the simplest nontrivial examples of the "classical" LBP with respect to a nonclassical "derivative" operator D.
Consider also a special case when both α and β are purely imaginary parameters: Then we can put (obviously the moments c n are defined up to a unnecessary common constant factor) by the expression c n = n + i(s 1 + s 2 ) n + is 1 from which the conditionc n = c −n follows which means that the corresponding polynomials P n (z) satisfy the Szegő recurrence relation (4.3) and one can expect that these polynomials are orthogonal on the unit circle with respect to a positive measure. According to the general theory this occurs if and only if the reflection parameters a n = −P n+1 (0) satisfy the restriction |a n | < 1 for all n. from (11.2) we see that a n−1 = −B n . Then After simple calculations we get |a n−1 | 2 = s 2 1 n 2 + s 2 1 (s 1 n + s 2 ) 2 + n 2 (s 1 n + s 2 ) 2 = 1 + ξ 2 n 1 + η 2 n , where η n = n/s 1 , ξ n = n/(ns 1 + s 2 ).
These polynomials coincide with a special class of Jacobi polynomials which are orthogonal on the interval [0, 1] with the weight function w(x) = x α+j : 1 0 x α+j W (j) n (x)W (j) m (x)dx = q n δ nm .
Hence, the polynomials W (j) n with nonzero β correspond to adding a concentrated mass M = β −1 at the endpoint x = 1 of orthogonality interval for the Jacobi polynomials. I.e. the weight function for the polynomials W Such polynomials are called the Krall-Jacobi polynomials (see, e.g., [19]). These polynomials have a remarkable property: they are eigenfunctions of a fourth-order differential operator [19,20].
If we now substitute j → j − n into above formulas, we return to the Laurent biorthogonal polynomials P (j) n (z) which satisfy a differential equation of the 4th degree which can be presented in the form n 2 L 2 + nL 1 + L 0 P (j) n (z) = 0, (11.12) where L 0 is a differential operator of the 4th order, L 1 of the 3-rd order and L 2 of the second order. The operators L 0 , L 1 , L 2 do not depend on the parameter n. We thus see that the Laurent biorthogonal polynomials P (j) n (z) given by (11.2), satisfy the quadratic operator pencil equation (11.12) (with respect to the "eigenvalue" parameter n). When β = ∞ we have generalized eigenvalue problem (11.10) which is equivalent to a linear operator pencil.
Hence the Laurent biorthogonal polynomials (11.12) can be considered as biorthogonal analogs of the Krall-Jacobi orthogonal polynomials. They possess many useful properties of the Krall-Jacobi polynomials including the 4th order differential equation they satisfy. The main difference, however, is that in the biorthogonal case we have quadratic operator pencil equation (11.12) instead of usual eigenvalue problem (11.11).
We thus see that many nontrivial properties of the elliptic Laurent biorthogonal polynomials are manifested already in the rational limit. Loosely speaking, one can say that elliptic biorthogonal polynomials (6.5) are elliptic analogs of the Krall-Jacobi polynomials.