Modulated Bi-Orthogonal Polynomials on the Unit Circle: The 2j-k\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2j-k$$\end{document} and j-2k\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$j-2k$$\end{document} Systems

We construct the systems of bi-orthogonal polynomials on the unit circle where the Toeplitz structure of the moment determinants is replaced by det(w2j-k)0≤j,k≤N-1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\det (w_{2j-k})_{0\le j,k \le N-1} $$\end{document} and the corresponding Vandermonde modulus squared is replaced by ∏1≤j<k≤N(ζk-ζj)(ζk-2-ζj-2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\prod _{1 \le j < k \le N}(\zeta _k - \zeta _j)(\zeta ^{-2}_k - \zeta ^{-2}_j) $$\end{document}. This is the simplest case of a general system of pj-qk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pj-qk$$\end{document} with p, q co-prime integers. We derive analogues of the structures well known in the Toeplitz case: third order recurrence relations, determinantal and multiple-integral representations, their reproducing kernel and Christoffel–Darboux sum, and associated (Carathéodory) functions. We close by giving full explicit details for the system defined by the simple weight w(ζ)=eζ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ w(\zeta )=e^{\zeta }$$\end{document}, which is a specialisation of a weight arising from averages of moments of derivatives of characteristic polynomials over USp(2N)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{USp}(2N)$$\end{document}, SO(2N)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{SO}(2N)$$\end{document} and O-(2N)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{O}^-(2N)$$\end{document}.


List of symbols w k
The k-th Fourier coefficient of w(z), k ∈ Z D r (z, Z) The (n +3)×(n +3) master matrix of 2 j −k structure with offset r ∈ Z D r (z, Z) The determinant of D r (z, Z) E s (z, Z) The (n +3)×(n +3) master matrix of j −2k structure with offset s ∈ Z E s (z, Z) The determinant of E s (z, Z) D (

r ) n
The n × n matrix of 2 j − k structure with offset r ∈ Z D The j − 2k multiple integral with weight f P n (z; r ) The monic 2 j − k polynomial of the first kind with offset r ∈ Z of degree n Q n (z; r ) The monic 2 j −k polynomial of the second kind with offset r ∈ Z of degree n R n (z; s) The monic j − 2k polynomial of the first kind with offset s ∈ Z of degree n S n (z; s) The monic j −2k polynomial of the second kind with offset s ∈ Z of degree n h (r ) n The norm of 2 j − k bi-orthogonal polynomials g (s) n The norm of j − 2k bi-orthogonal polynomials Z n (z) The (n + 1)-vector of monomials of degrees zero to n P n (z; r ) The (n + 1)-vector of P-polynomials of degrees zero to n Q n (z; r ) The (n + 1)-vector of Q-polynomials of degrees zero to n R n (z; s) The (n + 1)-vector of R-polynomials of degrees zero to n S n (z; s) The (n + 1)-vector of S-polynomials of degrees zero to n h (r ) n The (n +1)×(n +1) diagonal matrix of 2 j −k norms g (s) n The (n +1)×(n +1) diagonal matrix of j −2k norms P (r ) n and Q (r ) n The (n + 1) × (n + 1) lower triangular matrices in the LDU decomposition of D The coefficient of z in R n (z; s) The coefficient of z in S n (z; s) K n (z, Z; r ) The reproducing Kernel of the 2 j − k bi-orthogonal polynomials L n (z, Z; s) The reproducing Kernel of the j − 2k bi-orthogonal polynomials δ (r ) n and η (r ) n The recurrence coefficients associated to P n (z; r ) β (r ) n and α (r ) n The recurrence coefficients associated to Q n (z; r ) κ (s) n and ρ (s) n The recurrence coefficients associated to R n (z; s) γ (s) n and θ (s) n The recurrence coefficients associated to S n (z; s) The basis for polynomials of degree at most n consisting of P m (z; r ), m = 0, 1, . . . , n Q (r ) n The basis for polynomials of degree at most n consisting of Q m (z; r ), m = 0, 1, . . . , n R (s) n The basis for polynomials of degree at most n consisting of R m (z; s), m = 0, 1, . . . , n S (s) n The basis for polynomials of degree at most n consisting of S m (z; s), m = 0, 1, . . . , ň P n (z; r ) and P n (z; r ) Associated functions corresponding to P n (z; r ) Q n (z; r ) and Q n (z; r ) Associated functions corresponding to Q n (z; r ) R n (z; s) and R n (z; s) Associated functions corresponding to R n (z; s) S n (z; s) and S n (z; s) Associated functions corresponding to S n (z; s) F 1 and F 2 Carathéodory functions L 1 Linear difference operator annihilating P n (z; r ) L 2 Linear difference operator annihilating Q n (z −2 ; r ) L 3 Linear difference operator annihilating R n (z 2 ; s) L 4 Linear difference operator annihilating S n (z −1 ; r ) P n (z; r ) Casorati matrix annihilated by L 1 Q n (z; r ) Casorati matrix annihilated by L 2 R n (z; s) Casorati matrix annihilated by L 3 S n (z; s) Casorati matrix annihilated by L 4

Motivation
The unitary group U (N ) with Haar (uniform) measure possesses the explicit character formula of Weyl [72] , ζ l := e iθ l ∈ T, θ l ∈ (−π, π], (1.1) where {ζ 1 , . . . , ζ N } ∈ Spec(U ) and T = {ζ ∈ C : |ζ | = 1}. This also has an interpretation as the joint probability density function (PDF) (see e.g. [30,Chapter 2]) for the eigenphases of the Dyson circular ensemble (CUE). One fundamental application of this formula is to characterise averages over U ∈ U (N ) of class functions c(U ), i.e. symmetric functions of the eigenvalues of U only. An example of such functions are products N j=1 w(ζ j ) where w(ζ ) may be interpreted as a weight function or density. Introducing the Fourier components {w l } l∈Z of this weight (1.2) due to the well known Heine identity [70] E U (N ) we recognise that this is equivalent to studying Toeplitz determinants. Intimately connected with the above problem are systems of bi-orthogonal polynomials on the unit circle and its relationship to general, non-hermitian (i.e. complex weight) Toeplitz determinants. The system of bi-orthogonal polynomials {ϕ n (z),φ n (z)} ∞ n=0 1 with respect to the weight w(ζ ) on the unit circle may be defined by the orthogonality relation Such averages over the unitary group are ubiquitous in many areas of mathematical physics, in particular the gap probabilities and characteristic polynomial averages in the circular ensembles of random matrix theory [2,3,30,34], the spin-spin correlations of the planar Ising model [50,63], the density matrix of a system of impenetrable bosons on the ring [31] and probability distributions for various classes of non-intersecting lattice path problems [29].
Our study is motivated by yet another application of random matrix techniques, in particular to analytic number theory. Random matrix models have been very successful in constructing conjectures for estimating the integral moments of central values in the U (N ) family of L-functions, through the works of [22,23,51]. This program was extended by Ali Altug et al [4] to the three other families of L-functions: those characterised by the symmetry types USp(2N ), SO(2N ) and O − (2N ), and where the statistic of concern was the n-th moment of the m-th derivative of the characteristic polynomial A with A ∈ USp(2N ), SO where G denotes USp, SO, or O − , and d A is the Haar measure on G, in the regime N → ∞ and fixed n, m. Employing similar techniques to [24] they found that the leading coefficient in the large-N expansion is proportional to the n-th derivative of e u T n, (u) where T n, (u) := det g 2 j−k+ (u) 0≤ j,k≤n−1 , (1. 6) for n ≥ 0, ∈ Z is fixed by the symmetry type and u ∈ C, and where g l (u) = 1 2πi T e ζ +uζ −2 ζ l+1 dζ = 1 (l + 1) 0 F 2 ; 1 2 l + 1, 1 2 (l + 1); 1 4 u . (1.7) Clearly for the USp(2N ), SO(2N ) and O − (2N ) types the relevant moment determinant has the structure det(w 2 j−k+r ) 0≤ j,k≤n−1 , (1.8) and the corresponding joint density function has the form 1≤ j<k≤n , ζ l := e iθ l ∈ T, θ l ∈ (−π, π]. (1.9) Note that this joint density function is not real and positive, unlike the Toeplitz case (1.1), and it will become clear that meaning can only be given to more restrictive classes of weights than applies in the Toeplitz case, see Theorems 3.1 and 3.2 for precise conditions on the existence of our systems. For example the system with a constant weight or even a terminating Laurent expansion of (1.2), i.e. a finite banded moment matrix, will not exist 2 . Of course one could take as the definition of a putative system the modulus of (1.9) however this has some disadvantages and we have chosen to pursue the analytic form here. If one takes the modulus of (1.9) it factorises into three factors 1≤ j<k≤n ζ k − ζ j 2 ζ k + ζ j , (1.10) and one can see that this also applies to the joint density function 1≤ j<k≤n (1.11) and thus any distinction is lost. More significantly is that the property of complex analyticity would be lost in some aspects of the theory, and in particular the differential structures whereby one requires spectral derivatives d/dζ of the bi-orthogonal polynomials and their associated functions -such derivatives form one member of key Lax pairs of the integrable system. We will finish our discussion of these issues with a final observation on the log-gas interpretation of Eq. (1.10): it has the repulsion of coincident eigenvalues θ k = θ j , k = j with strength β = 2 but also a repulsion when θ k = θ j ± π with half this strength. It is our goal to extend the theory of Toeplitz determinants and associated bi-orthogonal polynomial systems on the unit circle to the case where they have a structure pj − qk for two positive co-prime integers, based upon (1.8) and (1.9). In the first example of such an extension we consider the p = 2, q = 1 and p = 1, q = 2 cases exclusively, however we will only employ techniques which are capable of generalisation to the other cases. Analogous systems of bi-orthogonal polynomials on the line with the Hankel structure det(h pj+qk ) 0≤ j,k≤N −1 have been studied for some time -the area was initiated by Preiser [67]; general properties of these systems were investigated by Konhauser [52], Ilyasov [46], Iserles and Nørsett [47,48]; explicit examples of cases related to Laguerre polynomials by Konhauser [53], Carlitz [16,17], Genin and Calvez [37,38], Prabhakar [66], Srivastava [69], Raizada [68]; and to the Jacobi polynomials by Madhekar and Thakare [57][58][59][60][61]. As they stand the results reported in the above works do not translate directly into the ones we seek. Systems of bi-orthogonal Laurent polynomials [71] would be expected to provide an equivalent framework to the system we treat here, however we prefer our approach because of its direct linkage to the random matrix application described earlier.
In the random matrix literature systems of bi-orthogonal polynomials on the line are known as Muttalib-Borodin ensembles following the pioneering work of [64] and [10]; Claeys and Romano have derived recurrence relations and a scalar Riemann-Hilbert problem for general classes of weights [21]. See [32,33] for some selected recent developments. Other related lines of investigation are the multiple orthogonal polynomial ensembles studied by Kuijlaars and McLaughlin [55], and by Kuijlaars [56]. In the latter work higher rank (i.e. greater than two) matrix systems of bi-orthogonal multiple polynomials -the Nikishin systems and Angelesco systems -were discussed and only the case of the Angelesco system (see Eq. (4.8)) with p = 2 species of n 1 = n 2 = n particles {x k , k = 1, . . . , n and w 2 (x) = w 1 (−x) has any correspondence with our example -it is actually proportional to the square of (1.9).
It is also worth mentioning that on the Operator Theory side, the associated operators on L 2 (T) and their restrictions on the Hardy space H 2 (T) have been a subject of research in the last 25 years or so. In a series of works [41][42][43][44][45], Mark C. Ho introduced and made some fundamental studies on the 2 j−k operators which he referred to as slant Toeplitz operators. Later, Subhash C. Arora and Ruchika Batra studied the properties of pj −k operators ( p ∈ N ≥2 ) which they referred to as generalized slant Toeplitz or pth order slant Toeplitz operators in a collection of papers [5][6][7][8]. Particularly regarding the large-size asymptotic analysis of the pj − qk determinants, it is our hope that the interplay between Operator Theory and the Riemann-Hilbert method (which is rooted in the orthogonality structures studied this work and is the subject of a future paper) could be made somewhat tractable. Examples of such interplay for Toeplitz, bordered Toeplitz and Toeplitz+Hankel determinants can be respectively found in [9,27], and the introduction of [40].

Outline
Here we summarise the outline of our study. In Sect. 2 we present the definitions of the 2 j − k and j − 2k (master) determinants and the corresponding systems of biorthogonal polynomials. In this section we also introduce the two main tools for our analysis: the Dodgson Condensation identity and the multiple integral formulae for the 2 j −k and j −2k determinants. The LDU decompositions for the 2 j −k and j −2k determinants are also derived which will be useful for proving the Christoffel-Darboux identity. In Sect. 3 we provide the bordered determinant representations and prove the existence and uniqueness of the systems of bi-orthogonal polynomials provided that the associated determinants are non-zero. We also discuss the connection of 2 j − k and j − 2k determinants and bi-orthogonal polynomials. Finally in this section we express 2 j − k and j − 2k bi-orthogonal polynomials and reproducing kernels in terms of the corresponding master determinants. In Sect. 4 we derive pure-degree and pure-offset recurrence relations for the 2 j − k and j − 2k bi-orthogonal polynomials. We also discuss equivalent Dodgson condensation identities which result in the same recurrence relations and present several mixed recurrence relations involving both 2 j − k and j − 2k bi-orthogonal polynomials. In Sect. 4.2 we derive several relationships between polynomial tails, recurrence coefficients and determinants of the 2 j − k and j − 2k systems. In Sect. 5 we prove multiple integral formulae for the 2 j − k and j − 2k determinants, bi-orthogonal polynomials, and reproducing kernels. We use some of these formulae to derive representations of Q and R-polynomials in terms of P and S-polynomials, respectively. The multiple integral formulae are also useful in deriving the Christoffel-Darboux identity in Sect. 5.3. In Sect. 6 we introduce the associated functions and derive their multiple integral formulae representations. We also find the corresponding Casorati matrices and the first order recurrence relations they satisfy. In Sect. 7 we go back to the weight relevant to the L-functions of the symmetry types USp(2N ), SO(2N ) and O − (2N ). To study a first concrete example, in this subsection we specifically study the undeformed weight when u = 0 where we find explicit formulae for the determinants, Carathéodory functions, and the 2 j − k polynomials. Eventually, in Sect. 8, we will lay out a list of important open questions and the prospects of future work. Throughout the paper, to highlight the distinguishing features of these modulated bi-orthogonal systems and for the convenience of the reader we try to make a comparison with the Toeplitz ( j − k) theory whenever a result about 2 j − k and j − 2k systems is presented, mainly by referring to [35].

Definitions and Notations
In this section we will define the objects studied in the paper and introduce the necessary notations and conventions. Throughout the paper we respectively use j and k as indices referring to the rows and the columns, and we frequently use a boldfaced letter to distinguish the determinant of a matrix with the matrix itself: det M ≡ M. Let M be an n × n matrix. By M j 1 j 2 · · · j k 1 k 2 · · · k , and M j 1 j 2 · · · j k 1 k 2 · · · k , (2.1) we respectively mean the (n − ) × (n − ) matrix obtained from M by removing the rows j i and the columns k i , 1 ≤ i ≤ , and its corresponding determinant. Although the order of writing the row and column indices is immaterial for this definition, in this work we prefer to respect the order of indices j 1 < j 2 and k 1 < k 2 if 1 < 2 . We now recall the Dodgson Condensation identity 3 (see [1,14,36] and references therein) which is an important tool for deriving many important relationships between the objects studied in this work: Definition 2.1 For fixed offset values r , s ∈ Z, we define the following (n+3)×(n+3) master matrices of 2 j − k and j − 2k structure: is the k-th Fourier coefficient of the symbol w(z) = ∈Z w l z l .
In (2.3) and (2.4) for simplicity of notation we have suppressed the dependence of D r (z, Z) and E s (z, Z) on n. Also, throughout the paper, when z and Z are not distinguished as distinct independent variables we simply drop the arguments in the notation of the master matrices: We use the determinants of master matrices D r and E s to construct the 2 j − k and j − 2k systems of bi-orthogonal polynomials in Sect. 3, but since in all of those constructions either the last row or the last column is removed, the entry in (2.3) and (2.4) never plays a role for construction of the bi-orthogonal polynomials and thus is arbitrary for the purposes of this work. Let D (r ) n and E (s) n respectively denote the n × n matrices of 2 j − k and j − 2k structure and by D (r ) n and E (s) n denote their determinants: n and E (s) n can obviously be written in terms of the determinant of master matrices D r and E s , as D (r ) n = D r n n + 1 n + 2 n n + 1 n + 2 , and E (s) n = E s n n + 1 n + 2 n n + 1 n + 2 . (2.9) Also notice that (2.11)

Definition 2.2
For an integrable function f on the unit circle, we respectively define the 2 j − k and j − 2k multiple integrals as and In Sect. 5, in particular, we show the following multiple integral representations for D (r ) n and E (s) n : (2.15) Definition 2.3 For each offset value r ∈ Z, define the 2 j − k sequences of monic polynomials {P n (z; r )} ∞ n=0 and {Q n (z; r )} ∞ n=0 , deg P n (z; r ) = deg Q n (z; r ) = n, satisfying the bi-orthogonality condition: where h (r ) n and g (s) n are the norms of the polynomials squared and dμ(ζ ) ≡ w(ζ )dζ for some weight function w(z).

LDU Decomposition of 2j − k and j − 2k Moment Matrices
The linear space of polynomials {P n (z)} ∞ n=0 , {Q n (z)} ∞ n=0 , etc. have expansions in an appropriate basis, in this case the monomial basis {z n } ∞ n=−∞ as befits a Páde approximation problem with two fixed singularities 0, ∞. We are therefore led to consider the linear transformations from this preferred basis to our orthogonal system. Let us denote the coefficients of 2 j − k and j − 2k polynomials more precisely as n,n = 1. Let us also denote We thus have where P (r ) n , Q (r ) n , R (s) n and S (s) n are the following (n + 1) × (n + 1) lower triangular matrices whose inverses are also lower diagonal with 1's on the main diagonal. Also let us denote the diagonal matrices of norms of polynomials by h (r ) n and g (r ) n : In the following result we give the LDU decompositions of the moment matrices for the 2 j − k and j − 2k systems which parallels the Toeplitz case closely, see the unpublished work [62].

Theorem 2.1 The LDU decompositions of D (r )
n+1 and E (s) n+1 are given by Proof In this proof we drop the dependence of all functions and quantities on the offsets r and s for simplicity of notation. We have which is equivalent to (2.28). For the decomposition of E (s) n+1 we consider which yields (2.29).

Bordered Determinant Representations
In this section we focus on determinantal representations of fundamental elements in the theory where the moment matrix is bordered by rows or columns containing basis vectors for the polynomial spaces. In the following result we find that the bi-orthogonal polynomials of both systems can be represented as bordered moment determinants in almost exactly the same way as for the Toeplitz case, see e.g. Eq.(2.15,16) of [35] for comparison. As well as revealing the conditions on the existence and uniqueness of these systems in a simple way, this form will be very useful in our subsequent treatment.

Theorem 3.1 If D
(r ) n = 0, the polynomials P n (z; r ) and Q n (z; r ) exist and are uniquely given by from which one can observe that h (r ) n exists and can be written as Therefore if all h (r ) exist and are non-zero for = 0, . . . , n − 1, then The existence simply follows from the fact that, if D (r ) n = 0, we can explicitly construct the system of monic bi-orthogonal polynomials P n (z; r ) and Q n (z; r ) as in (3.1) and (3.2), respectively satisfying (2.18) and (2.19). Now, let us discuss the uniqueness if D So, if D (r ) n = 0, the above linear system has a unique solution, and thus P n (z; r ) is uniquely given by (3.1). Now, assume that Q n (z; r ) = z n + n−1 =0 q (r ) n, z satisfies the orthogonality conditions (2.19). We can write the orthogonality conditions (2.19) for m = 0, 1, . . . , n − 1 as the following linear system for solving the constants q (r ) n, , 0 ≤ ≤ n − 1: So, again, if D (r ) n = 0, the above linear system can be inverted, and thus has a unique solution. Therefore, Q n (z; r ) is uniquely given by (3.2). Finally, one can directly find (3.3) using (3.1) and (2.18), or alternatively, using (3.2) and (2.19).
In an identical manner, we can prove the following Theorem about existence and uniqueness of j − 2k system of bi-orthogonal polynomials.

Theorem 3.2 If E (s)
n = 0, the polynomials R n (z; s) and S n (z; s) exist and are uniquely given by

7)
and S n (z; s) = 1 from which one can observe that g (s) n exists and can be written as Therefore if all g (s) exist and are non-zero for = 0, . . . , n − 1, then (3.10)

Connection of 2j − k and j − 2k Polynomials
The 2 j − k and j − 2k systems are intimately related by a duality through exploiting the freedom to choose suitable offsets. This equivalence is only exhibited in a formal algebraic rather than an analytical sense and reflects a mapping of the interior of the unit circle to the exterior and vice-versa. In the Toeplitz case one has a system of self-duality. Consequently the first of such duality relations involves the determinants, which are related to one another if one selects the appropriate offset values. More precisely, we have D (r ) n = E (r +n−1) n , (3.11) as these are reflections of each other across the main anti-diagonal. Indeed, for any n×n matrix M we have det M ⊥ = det M, which is due to the identity M ⊥ = A n M T A n , where A n is the n × n matrix with ones on the anti-diagonal and zeros everywhere else, M T is the transpose of the n × n matrix M, and M ⊥ is the reflection of the n × n matrix M across its main anti-diagonal. Continuing the development of the duality theme our next result furnishes further details in regard to the polynomials themselves. But first let us recall a standard notation: for any polynomial p of degree n, we denote the reciprocal polynomial by p * , that is The following identities hold between 2 j − k and j − 2k polynomials Proof Computing S * n (z; s) from (3.8) in view of (3.12) and reflection across the antidiagonal gives Now, we do n consecutive adjacent row-swaps to move the first row to the last which results in (3.13) in view of (3.1). The relationship (3.14) can be found in an identical manner.

Biorthogonal Polynomials in Terms of the Master Determinants
For many choices of triples {a, b, c} ∈ Z 3 , one can express z a f n+b (z; r + c), f ∈ {P, Q * , R, S * }, in terms of the determinants of master matrices (2.3) and (2.4). For each polynomial P, Q * , R and S * , among the choices mentioned above, we select five which turn out to be useful in connection to the specific Dodgson condensation identities used in Sect. 4. (3.35)

The Reproducing Kernel in Terms of the Master Determinants
In the same spirit as the j − k bi-orthogonal polynomials on the unit circle, we define the reproducing kernels for the 2 j − k and j − 2k systems.

Definition 3.1
The reproducing kernel for the 2 j − k and j − 2k systems respectively are respectively defined as It is easy to see that K n (z, Z; r ) and L n (z, Z; s) has the reproducing properties as described in the following equations: and Key properties flowing from these definitions are the normalisation relation (3.42) and the projection relation Identical relations apply for the L n kernel with the obvious modifications. The above definitions of the reproducing kernels follows the pattern of the Toeplitz case and the following result linking them to bordered moment determinants reinforces this similarity. Although the Toeplitz case is not quoted often it is almost identical after making the obvious modifications.

Theorem 3.4
The reproducing kernels K n (z, Z; r ) and L n (z, Z; s) can be represented in terms of the master determinants D r (z; Z) and E s (z; Z) as Proof We only give the proof for (3.44), as (3.45) can be obtained in a similar way. Define One can easily establish that K has the same reproducing properties as K . To this end, its action on monomials is given by and Therefore, we immediately get and for all z, Z ∈ C and r ∈ Z. Let us consider the following (n + 2) × (n + 2) matrices Our first result yields extensions of the second order scalar difference equations given in Eqs. (2.23,24) of [35] for the 2 j − k system.

Theorem 4.1
The third order pure-degree recurrence relations for the 2 j − k polynomials are given by Proof Let us consider the following Dodgson condensation identities (see (2.2)): The identity (4.7) can be written as So from (3.16), (3.17), (3.19), and (2.9) we have Let us also write (4.8) as Therefore from (3.16), (3.18), (3.24), and (2.9) we arrive at and thus From (4.11) we have z P n (z; r − 1) = P n+1 (z; r ) − δ (r ) n P n (z; r ). (4.15) From this, we immediately get Using (4.15) again we get and therefore On the other hand, from (4.14) we have Combining the last two equations with (4.15) and some simplification gives (4.1).
In order to prove (4.2) we need to consider the following Dodgson condensation identities: Let us write (4.20) as So from (3.26), (3.28), (3.29), and (2.9) we have (4.23) and thus from (3.3) we get Now, let us write (4.21) as  (4.26) and thus Now, let us find the pure-n recurrence relation for Q * s. From (4.24) we have From (4.27) and (4.28) we get and immediately and Using (4.30) we have Plugging this into (4.24) yields Replacing n by n + 2 and rearranging terms gives the desired recurrence relation (4.2).
Our second result gives extensions of the second order scalar difference equations given in Eqs. (2.23,24) of [35] for the j − 2k system.

Theorem 4.2
The third order pure-degree recurrence relations for the j − 2k polynomials are given by Proof In order to prove (4.34) we use the following Dodgson condensation identities (4.40) Now, lets us write (4.41) as Therefore from (3.20), (3.22), (3.23), and (2.9) we arrive at and thus and by writing (4.47) for s − 1 we obtain Combining this with (4.44) yields Plugging this into (4.48) gives Now, write (4.44) for s + 1 to get After combining the last two equations we arrive at Note that from (4.47) we have Combine the last two equations to get Rearranging terms and using ρ By writing (4.50) for n − 1 we obtain Now, we combine the last two equations to arrive at regroup terms and replace n by n + 1 to arrive at the desired recurrence relation (4.34). Finally, proving (4.35) demands using the following Dodgson condensation identities (4.60) Write (4.59) as So we have and thus from (3.9) we have Now, rewrite (4.60) as and thus Now let us find the pure degree-recurrence relations for S * 's. Writing (4.63) for n + 1 we get Now we write (4.66) for s + 1: Combining the last two equations yields Write this for n → n + 1: and write (4.63) for s → s − 1 to obtain Combining the last two equations gives Combine the last two equations to get (4.74) and rearranging terms we obtain Write (4.63) for n → n + 1 and solve for S * n+1 (z; s + 1): Now, we combine the last two equations and regroup the terms to obtain (4.35).

Remark 4.3 (Mixed Recurrence Relations)
Recall that we can replace a 2 j −k polynomial in terms of a j − 2k polynomial using Theorem 3.3. Using such replacements in several identities found in this section, we can arrive at recurrence relations involving both 2 j − k and j − 2k polynomials. To that end, (4.11) can be written as Also (4.14) can be written as one of the following three identities (4.80) In an identical way we can obtain the following mixed recurrence relations listed below and   Proof Combining (4.11) and (4.14) we will get η (r ) n P n (z; r + 2) − z P n (z; r + 1) + z P n (z; r − 1) + δ (r ) n P n (z; r ) = 0. (4.97) Also from (4.24) and (4.27) we obtain The replacement s → s + 2 in (4.99) and (4.100) respectively gives (4.95) and (4.96).

Equivalent Dodgson Condensation Identities
If one starts with after similar simplifications one arrives at But this is equivalent to (4.11) as one gets the above equation from (4.11) under the replacement r → r + 2. Similarly, If we start with D r n + 2 0 · D r 0 n + 1 n + 2 0 n + 1 n + 2 = D r 0 n + 2 0 n + 1 · D r n + 1 n + 2 0 n + 2 − D r n + 1 n + 2 0 n + 1 · D r 0 n + 2 0 n + 2 . (4.103) After simplifications, this translates to However, this is equivalent to (4.24) if we replace r by r − 1 in (4.24). Also note that reduces after simplifications to which is equivalent to (4.63) if one replaces s with s − 2 in (4.63). In a similar fashion one can find four other Dodgson Condensation identities which give rise to recurrence relations equivalent to (4.14), (4.27), (4.47), and (4.66), these Dodgson Condensation identities are respectively given by (4.112)

Relationships of Polynomial Tails, Recurrence Coefficients and Determinants
There are many redundancies amongst the various coefficients that have been defined thus far. In this subsection we provide the necessary relationships between these coefficients to remove the redundancy, but also to link these to certain polynomial coefficients, in particular the tail coefficients, and to the determinants and thus the norms of the polynomials. Such relations form an algorithmic path to compute these norms.
Proof These relationships simply follow from the definitions of the recurrence relations in terms of the 2 j − k and j − 2k determinants.

Lemma 4.6
The 2 j − k and the j − 2k systems each have only two independent recurrence coefficients due to the following interrelationships: Proof From (4.24) we have by the definition of δ (r ) n in (4.11). From (4.27) we have (4.120) Also from (4.14) we have (4.121) Combining the last two equations gives again by the definition of δ (r ) n in (4.11). The relations (4.117) and (4.118) can be proven similarly.
The following result relates the tail coefficients of the polynomials to products of certain recurrence relation coefficients.  Proof From (4.115) and (4.117) we have The following result relates the tail coefficients of the polynomials to certain ratios of determinants. This is the extension of the Toeplitz relations for the reflection or Verblunsky coefficients given by Eqs. (2.19) of [35].

Multiple Integral Representations
The joint eigenvalue PDF of (1.9) is fundamental from our viewpoint and it is from this that we construct various marginal distributions, otherwise known as n-particle correlation functions, through integration over a set of "internal" variables, the complement of n "external" variables. All aspects of the theory will have representations of this form: the determinants, the bi-orthogonal polynomials and the reproducing kernels. Thus our first result of this nature is the 0-point correlation or normalisation integral. This is the extension of the Toeplitz result, Eq.(2.2), of [35]. So, in particular we obtain (2.14) and (2.15) in view of (2.12), (2.13), (5.1) and (5.2). The relationship (3.11) can be also found using the multiple integral formulae above. To that end notice that

Theorem 5.1 The moment determinants, and therefore the normalisation of the biorthogonal polynomials, possess the multiple integral representations
which is equivalent to (3.11). More generally we have the following result. (5.11) or equivalently

Theorem 5.2 For any integrable function f we have:
(5.12)

Multiple Integral Formulae for the Bi-orthogonal Polynomials
It is well known that the bi-orthogonal polynomials in the Toeplitz system can be written as averages of characteristic polynomials over the unitary group, see Eq. (2.16,18) of [35]. In the following Theorem we give the analogous results for the 2 j − k and j − 2k polynomials. Another interpretation of such integrals is as rational modifications of the weight, known as Christoffel-Darboux-Uvarov transformations of the underlying bi-orthogonal system in the Toeplitz case, which was exhaustively treated in [73].

Theorem 5.3
We have the following multiple integral representations for 2 j − k and j − 2k polynomials: Proof We only prove the formula for P n as the other ones can be found similarly.
where we have denoted z by ζ 0 . Therefore (5.18)

Multiple Integral Formulae for the Reproducing Kernels
We will require the following important lemma, which is known as the "integration over successive variables", in a subsequent result. This is the extension of the result in the Toeplitz case, which is easy to derive following the line of Proposition 5.1.2 of [30].

(5.20)
Proof Let us use j and k respectively as indices of rows and columns. Expanding the determinant in the integrand along the last row gives det 1≤ j≤m 1≤k≤m Recalling the normalisation (3.42) and the projection relations (3.43) for the reproducing kernel we deduce By m − k − 1 successive swaps of adjacent columns we have Equation (5.19) now follows by combining the last two equations. The proof of (5.20) is similar and we do not provide the details here.

Theorem 5.5
The reproducing kernels for 2 j −k and j −2k systems have the following multiple integral representations: (5.24) and L n (z 2 , z 1 ; s) = 1

(5.25)
Proof For simplicity of notation, in this proof we denote P n (z; r ), Q n (z; r ) and K n (z 2 , z 1 ; r ) respectively by P n (z), Q n (z) and K n (z 2 , z 1 ). We have where we have temporarily denoted z 1 and z 2 respectively by ζ 0 and ζ −2 0 . Therefore . Now we apply Theorem 5.4 : . (5.26) Repeating this n − 1 times yields (5.27) and n+1 K n (z 2 , z 1 ). The equation (5.25) can be established similarly and we do not provide the details here.

Remark 5.6
One can also see the consistency of multiple integral formulae for the reproducing kernels with the definitions (3.36) and (3.37) by considering the leading asymptotic behaviour as z or Z tend to infinity. For example as z → ∞ we have and, expectedly, from the right hand side of (5.24) we find the same asymptotic behaviour: by (5.13). Similar consistency checks can be done for K n (z, Z; r ) as Z → ∞, and also for L n (z, Z; r ) as z → ∞ and separately as Z → ∞.

Theorem 5.7 We have
S (z; s + 2). (5.31) Proof These identities can be proven if one considers special arguments for the reproducing kernel. Indeed, to prove (5.28) we consider n+1 P n (z; r + 2), (5.33) where the last equality follows from (5.13). From (4.135) we have Now (5.28) follows from (5.32), (5.33), and (5.34). The identity (5.29) can be proven by considering: and , (5.37) which is obtained by (4.134). (5.30) and (5.31) can be proven similarly, respectively by considering L n (0, z; s) and L n (z, 0; s) and employing , and R n (0; s) which are consequences of (4.136) and (4.137) respectively. We obviously assume that the corresponding determinants do not vanish so that the polynomials exist (and are unique, see Theorems 3.1 and 3.2). The following theorem characterizes the number of polynomials needed to represent P n (z, r ), Q n (z, r ), R n (z, s) and S n (z, s) respectively in terms of polynomials in the bases P (r ) n , Q (r ) n , R (ŝ) n , and S (ŝ) n for various values ofr andŝ. . While if it is desired to express P n (z; r ) in terms of polynomials in the basis P (r +2m) n , then one needs to only use the following m + 1 polynomials in the subset {P n (z; r + 2m), P n−1 (z; r + 2m), . . . , P n−m (z; r + 2m)} ⊂ P (r +2m) n .

2.
If it is desired to express z P n (z; r ) as a linear combination of polynomials in P , then one needs to use all n + 2 polynomials in P . While if it is desired to express z P n (z; r ) in terms of polynomials in P (r +2m−1) n+1 , then one needs to only use the m + 1 polynomials in the subset {P n+1 (z; r + 2m − 1), P n (z; r + 2m − 1), . . . , P n+1−m (z; r + 2m − 1)} 3. If it is desired to express Q n (z; r ) as a linear combination of polynomials in Q (r +m) n , then one needs to use all n + 1 polynomials in Q

If it is desired to express S n (z; s) as a linear combination of polynomials in S
Proof Theorem 5.7 establishes the statement of this theorem for expressing P n (z, r ), Q n (z, r ), R n (z, s) and S n (z, s) respectively in terms of polynomials in P all exist and are non-zero. Let us now prove the first item for expressing P n (z; r ) in the basis P (r −2m) n via induction. The base of the induction for m = 1 clearly holds because of (5.28). Passing from the induction hypothesis to the desired result is obvious because by (5.28) we have P (z; r − 2(m − 1)) = ν=0Â ν, P ν (z; r − 2m), where the constantsÂ ν, are also described in (5.28). Thus the induction hypothesis P n (z; r ) = n =0 A m,n P (z; r −2(m −1)), immediately implies the desired result. To prove the first item for expressing P n (z; r ) in the basis P (r +2m) n we use induction again. The base of the induction clearly holds due to (5.39). Now the induction hypothesis B m−1,n P n− (z; r + 2(m − 1)), (5.44) in conjunction with In the Toeplitz case there exist relations which allow one of the polynomials to be expressed in terms of a linear combination of two of the other one, see Eq. (2.21,22) of [35]. These two relations are also the components of the first order recurrence, Eq.(2.70), in the matrix variable, Eq.(2.69). Furthermore these two relations are somewhat symmetrical. However for the present system we find the following single, unsymmetrical relation, of Q n in terms of a bilinear combination of the P n , and the analogous relation for the R n , S n system.

Remark 5.10
Since we observe that (5.48) is compatible with (4.135). A similar observation shows the compatibility of (5.49) with (4.136).

The Christoffel-Darboux Identity
Having defined the reproducing kernel via a sum over products of bi-orthogonal polynomials and deduced that these polynomials satisfy third order difference equations we seek the explicit evaluation of the sum involved in the kernel. The well known Christoffel-Darboux sum in the Toeplitz result can be found in Prop. 2.5 of [35].

Now (5.55) follows immediately by noticing
Remark 5. 12 We have found that the sum is most simply evaluated in terms of just the polynomials P n and is therefore unsymmetrical with respect to Q n . It is possible to replace certain terms in this determinant by Q n using (5.48) but not all in a symmetrical way.

Associated Functions
Our goal in this section is to define linearly independent associated functions corresponding to the polynomials P n (z; r ), Q n (z; r ), R n (z; s), and S n (z; s) each of which respectively satisfies the recurrence relations (4.1), (4.2), (4.34), (4.35) or their equivalents. To this end, let us define the associated functionš It is also convenient to write these associated functions aš P n (z; r ) = P n (z; r ) − F 1 (z; r )P n (z; r ), (6.5) Q n (z; r ) = Q n (z; r ) − F 2 (z; r )Q n (z −2 ; r ), (6.6) R n (z; s) = R n (z; s) − F 2 (z; s)R n (z 2 ; s), (6.7) S n (z; s) = S n (z; s) − F 1 (z; s)S n (z −1 ; s). (6.8) where and (6.14) Notice that due to For the same reason we have the following expressions for P n and S n P n (z; r ) = 1 2 P † n (z; r ) + P † n (−z; r ) , (6.17) where In the following theorem we show that the associated functions above satisfy the same recurrence relations as the orthogonal polynomials. To this end, let us first define the following linear difference operators n )z 2 f n (z; r ), (6.21) n z 2 f n (z; s), (6.23) and z −2 f n (z; s). Proof We have n )z 2 P n (z; r ) . (6.30) Writing z 2 = ζ 2 + (z 2 − ζ 2 ) and using (4.1) for polynomials in z and ζ we obtain where we have used which can be seen from (3.1). Indeed, using the determinantal representation of P n (z; r ) we can write T dζ 2π iζ w(ζ )ζ 2−r P n (ζ ; r ) (6.34) Now consider the following Dodgson condensation identity: which can be written as Combining this with (6.34) yields L 1 [P n (z; r )] = 0. It is obvious that the same argument, in particular, shows that L 1 [ P n (z; r )] = 0. The proof of (6.27), (6.28), and (6.29) follow from similar considerations and we do not provide the details here. We just point out that the proof of (6.28) is achieved much more easily if one first proves that is annihilated by (6.23), and then deduce (6.28) simply in view of (6.15). Following this path one needs to show that vanishes (dμ(ζ ; s) = w(ζ )ζ −s dζ ), while in the "direct" proof to show that (6.11) is annihilated by (6.23), one is required to show that the expression above is zero when we replace the (ζ + z) by (ζ + z) 2 . This is not as easy, since the presence of the ζ 2 in the integrand requires us to deal with bordered j − 2k determinants.
One is expected to find three linearly independent solutions to the third order linear difference equations L i [ f n (z; r )] = 0, i = 1, 2, 3, 4. The following corollary gives a natural fundamental set of solutions for these difference equations, which can be immediately concluded from Theorem 6.1 due to presence of z 2 , as opposed to odd functions of z, in the linear difference operators (6.21) through (6.24). Corollary 6.1.1 A fundamental set of solutions for the recurrence relations (6.21), (6.22), (6.23) and (6.24) is respectively given by P n (z; r ), P n (−z; r ), P n (z; r ) , Q n (z −2 ; r ), Q n (z; r ), Q n (−z; r ) , R n (z 2 ; s), R n (z; s), R n (−z; s) and S n (z −1 ; s), S n (−z −1 ; s), S n (z; s) , where P n (z; r ), Q n (z; r ), R n (z; s) and S n (z; r ) are given by (6.9), (6.10), (6.11), and (6.12).
In Theorem 5.3, Eqs.(5.13) through (5.16), we gave representations of the biorthogonal polynomials as averages of characteristic polynomials and in the following result we give the analogous result for the associated functions as averages of reciprocals of characteristic polynomials. This also corresponds to formulae Eq.(2.34,35) of [35] or Eq.(3.3,4) of [73] for the Toeplitz case. Theorem 6. 2 We have the following multiple integral representations for associated functions P n (z; r ), Q n (z; r ), R n (z; s), and S n (z; s) defined respectively by (6.9), (6.10), (6.11), and (6.12): By symmetry, for each choice of m we have the same object, so we set m = n + 1 and write Now, recalling (2.12) and (5.14) we observe that Now (6.37) follows by recalling (6.10) and noticing A key step in understanding the recurrence structures in the 2 j − k and j − 2k systems is their formulation as a first order 3 × 3 matrix difference equation. The matrix variable so defined would then be come a central object in the integrable system and possibly related to the solution of a rank 3 Riemann-Hilbert problem. One choice for the 2 × 2 matrix variable in the Toeplitz case was given by Eq.(2.69) along with its matrix recurrence relation Eq.(2.70) (other choices were given in Eqs.(2.73,76)) of [35].

Definition 6.1
The Casorati matrices (see e.g. [28]) associated with the recurrence relations (6.21) through (6.24) are defined as 52) (6.53) A necessary step in establishing the existence of a valid 3 × 3 matrix system is to compute their Casoratians and to demonstrate their boundedness and non-vanishing character for z = 0, ∞. Lemma 6.4 (Abel's Lemma, [28]) The Casoratians for the matrix systems (6.46) through (6.48) are respectively given by .

(6.57)
Proof Since the Casorati matrix P n (z; r ) satisfies the recurrence relation (6.21), its determinant satisfies the recurrence relation which follows from rewriting the objects with index n + 3 in the last row of det P n+1 (z; r ) in terms of the objects with indices n, n + 1, and n + 2 using (6.21) and performing row operations. Using (6.58) we can write det P n (z; r ) det P 0 (z; r ) = (−1) n z 2n The formulas (6.55), (6.56), and (6.57) can be established similarly.

The Symplectic and Orthogonal Example
Now we return to the specific case of 2 j − k determinants which appeared in [4]. As we discussed in Sect. 1, it is shown in [4] that Analysing this particular determinant in full details will be carried out in our future work, but in this section as a concrete first example, we are going to present explicit formulae for the determinants, orthogonal polynomials, and recurrence coefficients corresponding to the undeformed weight w(z; 0) ≡ w(z) = e z . (7.2) Notice that for this weight 3) The entries of the determinant exist if r ∈ Z ≥0 . Let us now recall a result in [65] which is all we need to evaluate D Using the functional equation for the Barnes G-function G(z + 1) = (z)G(z), we can write n−1 j=1 (1 + j) = G(n + 1). Also using the identity (2z) = π − 1 2 2 2z−1 (z) (z + 1 2 ), we write (1 + r + 2 j) = π − n 2 2 n(n+r −1) G n + r +1 2 G n + r +2 2 G r +1 2 G r +2 2 .

Prospects of Future Work
In our study of 2 j − k and j − 2k bi-orthogonal polynomial systems on the unit circle we have concentrated on laying out the theory for generic class of weights and only considered the essential orthogonality structures. Even within this circumscribed area there are other important aspects that are yet to be investigated. Our purpose in this section is to bring the 2 j − k and j − 2k systems to the attention of a wider audience of mathematicians by providing a (non-exhaustive) list of open problems which we believe are significant for the development of the theory of 2 j − k and j − 2k systems, and consequently for the theory of pj − qk systems, with relatively prime p and q.
One example of this appears to be the novel form of the joint density function for the pj − qk systems which can be rewritten as the product of differences and features a primary repulsion of β = 2 as θ k −θ j → 0 along with weaker repulsions of β = 1 at as θ k − θ j → ± 2 p π, ± 2 q π, . . .. There has been a growing interest in recent years about the asymptotic aspects of structured moment determinants other than the well-known cases of Toeplitz and Hankel. Among those studies are asymptotics of Toeplitz+Hankel determinants [26,40], and recently the bordered Toeplitz determinants [9]. A successful Riemann-Hilbert formulation for 2 j − k and j − 2k systems places them among the collection of structured moment determinants for which the large size asymptotics of the determinant and the large degree asymptotics of the corresponding systems of bi-orthogonal polynomials can be investigated. In this regard, the formulation of a (3 × 3) Riemann-Hilbert problem corresponding to these bi-orthogonal polynomial systems, both as a means of founding the whole theory upon this and deriving key results but also to pave the way for a suitable Deift-Zhou analysis will be the topic of a future publication. At a later stage, the asymptotic description of these determinants can be investigated when the symbol w is of Fisher-Hartwig type, similar to what has been done for Toeplitz [26], Hankel [18,20], and Muttalib-Borodin [19] determinants.
On a different but possibly related perspective, it would be interesting to search for possible Fredholm determinant representations for the 2 j −k and j −2k determinants in the same spirit as the Fredholm determinant representations for Toeplitz determinants which could also unveil a Riemann-Hilbert formulation (see [25] where the connection was first found for Toeplitz determinants, and also [49]). Also looking for the possible 2 j − k/ j − 2k analogue of the Fredholm sine-kernel determinant representations for gap probabilities in random matrix theory seems to be another important front to be investigated (see e.g. [54] and references therein). Recalling the last paragraph of Sect. 1, it would be interesting to ask if from the viewpoint of Operator Theory, the analogue of the (strong) Szegő limit theorem can be proven for the 2 j − k/ j − 2k determinants (See. e.g. the monographs [12,13]). If this is plausible, we could hope for arriving at yet another example of close interaction between Operator Theory and Riemann-Hilbert techniques.
What is mentioned above is only a short list and obviously the full list of 2 j − k/ j − 2k analogues of Toeplitz theory can not be enumerated in full here, given the voluminous literature on the Toeplitz side. Nevertheless we point out some other obvious candidates: Rational approximations to the Carathéodory functions with a two-point (z = 0 and z = ∞) Hermite-Padé approximation of the associated functions; Quadrature problems on the unit circle and Christoffel weights; The elucidation of an analogue of the CMV matrix [15] in the Toeplitz case, i.e. the minimum banded representation of the spectral multiplication operator; and the analogue of the discrete Fredholm determinant arising from the scattering theory approach to the matrix difference system for the Toeplitz system, commonly known as the Geronimo-Case-Borodin-Okounkov identity [11,39].
Tasks which we envisage completing progressively include the representation of the derivatives of the polynomials and associated functions in terms of the fore-mentioned bases. In the context of semi-classical weights this will provide the pair of Lax operators of the integrable system. We anticipate making progress on the question raised in [4] It would be very interesting to determine whether or not there exists a differential equation arising from our formula (5) which plays a role for symplectic and orthogonal types that Painlevé III plays for unitary symmetry.