Orthogonal Structure on a Wedge and on the Boundary of a Square

Orthogonal polynomials with respect to a weight function defined on a wedge in the plane are studied. A basis of orthogonal polynomials is explicitly constructed for two large class of weight functions and the convergence of Fourier orthogonal expansions is studied. These are used to establish analogous results for orthogonal polynomials on the boundary of the square. As an application, we study the statistics of the associated determinantal point process and use the basis to calculate Stieltjes transforms.

We also study orthogonal polynomials on the boundary of a parallelogram. Without loss of generality, we can assume that our wedge is of the form f (x, 1)g(x, 1)w 1 (x)dx + 1 0 f (1, y)g(1, y)w 2 (y)dy. (1.2) Since is a subset of the zero set of a quadratic polynomial l 1 (x, y)l 2 (x, y), where l 1 and l 2 are linear polynomials, the structure of orthogonal polynomials on is very different from that of ordinary orthogonal polynomials in two variables [4] but closer to that of spherical harmonics. The latter are defined as homogeneous polynomials that satisfy the Laplace equation Y = 0 and are orthogonal on the unit circle, which is the zero set of the quadratic polynomial x 2 + y 2 − 1. The space of spherical polynomials of degree n has dimension 2 for each n ≥ 1 and, furthermore, one basis of spherical harmonics when restricted on the unit circle are cos nθ and sin nθ , in polar coordinates (r , θ), and the Fourier orthogonal expansions in spherical harmonics coincide with the classical Fourier series.
In Sect. 2, we consider orthogonal polynomials on a wedge. The space of orthogonal polynomials of degree n has dimension 2 for each n ≥ 1, just like that of spherical harmonics, and they satisfy the equation ∂ 1 ∂ 2 Y = 0. The main results are • An explicit expression in terms of univariate orthogonal polynomials when w 1 (x) = w 2 (x) = w(x) where w is any weight function on [0, 1] (Theorem 2.2), • Sufficient conditions for pointwise and uniform convergence (Theorem 2.4), as well as normwise convergence (Corollary 2.5), • Explicit expression in terms of Jacobi polynomials when w 1 (x) = w α,γ (x) and w 2 (x) = w β,γ (x) (Theorem 2.7), • Sufficient conditions for normwise convergence (Theorem 2.9).
In Sect. 3 we study orthogonal polynomials on the boundary of a parallelogram, which we can assume as the square [−1, 1] 2 without loss of generality. For a family of generalized Jacobi weight functions that are symmetric in both x and y, we are able to deduce an orthogonal basis in terms of four families of orthogonal bases on the wedge in Theorem 3.2. Furthermore, the convergence of the Fourier orthogonal expansions can also be deduced in this fashion, as shown in Theorem 3.3.
In Sect. 4 we use orthogonal polynomials on the boundary of the square to construct an orthogonal basis for the weight function w(max{|x|, |y|}) on the square [−1, 1] 2 . This mirrors the way in which spherical harmonics can be used to construct a basis of orthogonal polynomials for the weight function w( x 2 + y 2 ) on the unit disc. However, unlike the unit disc, the orthogonal basis we constructed are no longer polynomials in x, y but are polynomials of x, y and s = max{|x|, |y|}.
The study is motivated by applications. In particular, we wish to investigate how the applications of univariate orthogonal polynomials can be translated to multivariate orthogonal polynomials on curves. As a motivating example, univariate orthogonal polynomials give rise to a determinantal point process that is linked to the eigenvalues of unitary ensembles of random matrix theory. In Sect. 5, we investigate the statistics of the determinantal point process generated from orthogonal polynomials on the wedge, and find experimentally that they have the same local behaviour as a Coulomb gas away from the corners/edges.
In "Appendix A", we give the Jacobi operators associated with a special case of weights on the wedge, whose entries are rational. Finally, in "Appendix B" we show that the Stieltjes transform of our family of orthogonal polynomials satisfies a recurrence that can be built out of the Jacobi operators of the orthogonal polynomials, which can in turn be used to compute Stieltjes transforms numerically. This is a preliminary step towards using these polynomials for solving singular integral equations.

Orthogonal Polynomials on a Wedge
Let P 2 n denote the space of homogeneous polynomials of degree n in two variables; that is, P 2 n = span {x n−k y k : 0 ≤ k ≤ n}. Let 2 n denote the space of polynomials of degree at most n in two variables.

Orthogonal Structure on a Wedge
Given three non-collinear points, we can define a wedge by fixing one point and joining it to other points by line segments. We are interested in orthogonal polynomials on the wedge. Since the three points are non-collinear, each wedge can be mapped to by an affine transform. Since the polynomial structure and the orthogonality are preserved under the affine transform, we can work with the wedge without loss of generality. Henceforth we work only on .
Let w 1 and w 2 be two nonnegative weight functions defined on [0, 1]. We consider the bilinear form define on by Furthermore, we can choose polynomials in H 2 n (w 1 , w 2 ) to satisfy ∂ x ∂ y p = 0.
Proof Since (1−x)(1−y)P n−2 is a subset of I , the dimension of dim H 2 n (w 1 , w 2 ) ≤ 2. Applying the Gram-Schmidt process on {1, x k , y k , k ≥ 1} shows that there are two orthogonal polynomials of degree exactly n. Both these polynomials can be written in the form of p(x) + q(y), since we can use x y ≡ x + y − 1 mod I to remove all mixed terms. Evidently such polynomials satisfy ∂ x ∂ y ( p(x) + q(y)) = 0.
In the next two subsections, we shall construct an orthogonal basis of H 2 n (w 1 , w 2 ) for certain w 1 and w 2 and study the convergence of its Fourier orthogonal expansions. We will make use of results on orthogonal polynomials of one variable, which we briefly record here.
For w defined on [0, 1], we let p n (w) denote an orthogonal polynomial of degree n with respect to w, and let h n (w) denote the norm square of p n (w), The Parseval identity implies that The n-th partial sum of the Fourier orthogonal expansion with respect to w can be written as an integral where k n (w) denotes the reproducing kernel defined by

Orthogonal Structure for w 1 = w 2 on a Wedge
In the case of w 1 = w 2 = w, we denote the inner product (2.1) by ·, · w and the space of orthogonal polynomials by H 2 n (w). In this case, an orthogonal basis for H 2 n (w) can be constructed explicitly.
Then {P n , Q n } are two polynomials in H 2 n (w) and they are mutually orthogonal. Furthermore, P n , P n w = 2h n (w) and Q n , Q n w = 2h n−1 (φw). (2.5) for n ≥ 0 and m ≥ 1. Furthermore, by the orthogonality of p n (w). Similarly, The proof is completed.
Let L 2 ( , w) be the space of Lebesgue measurable functions with finite 1 2 norms. For f ∈ L 2 ( , w), its Fourier expansion is given by where P n and Q n are defined in Theorem 2.2 and The partial sum operator S n f is defined by which can be written in terms of an integral in terms of the reproducing kernel K n (·, ·), We show that this kernel can be expressed, when restricted on , in terms of the reproducing kernel k n (w; ·, ·) defined at (2.3).
It is well-known that the kernel k n (w; ·, ·) satisfies the Christoffel-Darboux formula, which plays an important role for the study of Fourier orthogonal expansion. Our formula allows us to write down an analogue of Christoffel-Darboux formula for K n (·, ·), but we can derive convergence directly.

Theorem 2.4 Let f be a function defined on . Define
Then Proof By our definition, × w(y)dy Similarly, × w(y)dy Moreover, since f e (x) from which we see that the convergence of s n (w; f e ) and s n (φw; f o ) imply the convergence of S n f .
In particular, s n (w, f e ) and s n (φw; f o ) converge to f e and f o in L 2 (w) and in L 2 (φw), respectively.
Proof By the displayed formulas at the end of the proof of the last theorem and it is easy to see that

Orthogonal Structure
We need to construct an explicit basis of H 2 n (w α,γ , w β,γ ). The case α = β can be regarded as a special case of Theorem 2.2. The case α = β is much more complicated, for which we need several properties of the Jacobi polynomials.
is an orthogonal polynomial of degree n with respect to m,n = 0 for n > m holds trivially. For m ≥ n, we need two identities of Jacobi polynomials. The first one is, see [12, (4.5.4)] or [9, (18.9.6)], and the second one is the expansion, see [9, (18.18.14)], Putting them together shows that Substituting this expression into I γ,α m,n and using the orthogonality of the Jacobi polynomials and (2.10), we conclude that, for m − 1 ≥ n, Hence, the case m > n follows. The same argument works for the case n = m.
What is of interest for us is the fact that the dependence of I γ,α m,n on n and α is separated, which is critical to prove that Q n in the next theorem is orthogonal. Theorem 2.7 Let P 0 (x, y) = 1 and, for n = 1, 2, . . ., define Then {P n , Q n } are two polynomials in H 2 n (w α,γ , w β,γ ) and (2.14) In particular, the two polynomials are orthogonal to each other if β = α. Furthermore By the identity in the previous lemma, if n > m, then P n , Q m α,β,γ = 0 since both The case n = m follows from a simple calculation. Moreover, for m = n, by the orthogonality of the Jacobi polynomials, and it is equal to h To derive the norm of Q n , Q n , we need to use c γ, The proof is completed.
Proof Since polynomials are dense on , by the Weierstrass theorem, the orthogonal basis {P n , R n } is complete, so that the Fourier orthogonal expansion converges in L 2 ( , w α,γ , w β,γ ). By the Parseval identity, Throughout this proof we use the convention where c 1 and c 2 are constants that are independent of varying parameters in A and B. By (2.10) and the fact that (n + α + 1)/n! ∼ n α , it is easy to see that h α,γ n ∼ n −1 , so that P n , P n α,β,γ ∼ n −1 , Q n , Q n α,β,γ ∼ n 2γ +3 , P n , Q n α,β,γ ∼ n γ , and, consequently, The Fourier-Jacobi coefficients of f 1 and f 2 are denoted by f 1 α,γ n and f 2 β,γ n , respectively. It follows readily that f P n ∼ f 1 We now consider the estimate for R n part. By the definition of R n , f , R n α,β,γ ∼ f , Q n α,β,γ − n γ +1 f , P n α,β,γ .
It is easy to see that so that we only have to work with the term f , Q k α,β,γ . The definition of Q k shows that 1, Q k α,β,γ = 0, which leads to the identity The proof is completed.

Orthogonal Polynomials on the Boundary of the Square
Using the results in the previous section, we can study orthogonal polynomials on a parallelogram. Since orthogonal structure is preserved under an affine transformation, we can assume without loss of generality that the parallelogram is the square [−1, 1] 2 . For α, γ > −1, let α,γ be the weight function We consider orthogonal polynomials of two variables on the boundary of [−1, 1] 2 with respect to the bilinear form vanishes on the boundary of the square, the bilinear form defines an inner product modulo the ideal generated by this polynomial, or in the space Let BV 2 n denote the space of orthogonal polynomials in R[x, y]/I with respect to the inner product ·, · . Proposition 3.1 For n ≥ 0, the dimension of BV 2 n is given by dim BV 2 n = n + 1, n = 0, 1, 2, and dim BV 2 n = 4, n ≥ 3.
Recall that the inner product ·, · α,β,γ studied in the previous section contains a fixed parameter σ . For fixed α, β and δ 1 , δ 2 ∈ {0, 1}, we define p Theorem 3.2 For n = 0, 1, 2, a basis for BV n is denoted by Y n,i and given by For n ≥ 3, the four polynomials in BV 2 n that are linearly independent modulo the ideal can be given by for n = 2m + 1 ≥ 3. In particular, these bases satisfy the equation ∂ 2 x ∂ 2 y u = 0.
Proof The proof relies on the parity of the integrals. For example, it is easy to see that x f (x 2 , y 2 ), g(x 2 , y 2 ) = 0 and y f (x 2 , y 2 ), g(x 2 , y 2 ) = 0 for any polynomials f and g, which implies, in particular, that Y 2m,i , Y 2n+1, j = 0 for i, j = 1, 2, 3, 4. Furthermore, it is easy to see that x y f (x 2 , y 2 ), g(x 2 , y 2 ) = 0 for any polynomials f and g. Hence, Y 2m,i , Y 2k, j = 0 for i = 1, 2 and j = 3, 4. Furthermore, using the relation it is easy to see that where in the second identity, we have adjusted the normalization constants of integrals from c α,γ and c β,γ to c α+1,γ and c β+1,γ , respectively, and used our choice of σ 1,1 . Hence, with our choice of σ 0,0 and σ 1,1 , we see that Y 2m,i is orthogonal to Y 2k, j for i, j = 1, 2 and i, j = 3, 4, respectively. Similarly, by the same consideration, we obtain that which shows, with our choice of σ 0,1 and σ 1,0 , that Y 2m+1,i is orthogonal to Y 2k+1, j for i, j = 1, 2 and i, j = 3, 4, respectively. Finally, since ∂ x ∂ y p α,β n,i (x, y) = 0, we see that Y n, j = ξ(x, y)u(x) + η(x, y)v(x), where ξ and η are linear polynomial of x, y, so that it is evident that ∂ 2 x ∂ 2 y Y n, j (x, y) = 0.
In our notation, the case α = − 1 2 β = − 1 2 and γ = 0 corresponds to the inner product in which the integrals are unweighted.
We further define o (x, y), and define ψ : In particular, the norm of S n f − f is bounded by those of S Proof Using the parity of the function, it is easy to see that where we have used the fact that F e,e is even in both variables and use the change of variables in integrals as in (3.3). The similar procedure can be used in the other three cases, as G i (x, y) is even in both variables, and the result is f (x, y) = G 0,0 (x, y) + yG 0,1 (x, y) + yG 1,0 (x, y) + x yG 1,1 (x, y), the last statement is evident.

Orthogonal System on the Square
Let w be a nonnegative weight function defined on [0, 1]. Define We construct a system of orthogonal functions with respect to the inner product by making use of the orthogonal polynomials on the boundary or the square, studied in the previous section. Our starting point is the following integral identity derived from changing variables (x, y) → (sξ, sη), where B dσ denotes the integral on the boundary of the square, Our orthogonal functions are similar in structure to orthogonal polynomials on the unit disc that are constructed using spherical harmonics. However, these function are polynomials in (s, ξ, η) for the (x, y) = (sξ, sη) ∈ [−1, 1] 2 , but not polynomials in (x, y).
Let BV 2 n be the space of orthogonal polynomials on the boundary of [−1, 1] 2 with respect to the inner product which is the inner product with α = − 1 2 , β = − 1 2 and γ = 0 studied in the previous section. Let Y n,i be an orthogonal basis for BV 2 n . For n ≤ 2, they are defined by, see Theorem 3.2, Proof Changing variables x = sξ and y = η shows The second integral is zero if i = j and n − k = m − l, whereas the second integral is zero when n − k = m − l and k = l, so that Q n k,i , Q m l, j W = 0 if i = j, k = l and n = m. By definition, s n−k Y n−k,i ξ s , η s is a polynomial of degree n − k in the variable s, so that Q n k,i is a polynomial of degree n. To show that the system is complete, we show that if f , Q n k,i = 0 for all k, i, n, then f (x, y) = 0. Indeed, by the orthogonality of polynomials on the boundary, we see that modulo the ideal. Changing order of summation shows that This completes the proof.

Sampling the Associated Determinantal Point Process
Associated with an orthonormal basis q 0 (x), . . . , q N (x) is a determinantal point process, which describes N points λ 1 , . . . , λ N distributed according to is the reproducing kernel, see [1] for an overview of determinantal point processes.
In the particular case of univariate orthogonal polynomials with respect to a weight w(x), the associated determinantal process is equivalent to a Coulomb gas-that is, the points are distributed according to where Z N is the normalization constant-as well as the eigenvalues of unitary ensembles, see for example [3] for the case of an analytic weight on the real line or [8] for the case of a weight supported on [−1, 1] with Jacobi-like singularities.
In the case of our orthogonal polynomials on the wedge, the connection with Coulomb gases and random matrix theory is no longer obvious: the interaction of the points is not Coulomb (that is, it cannot be reduced to a Vandermonde determinant squared times a product of weights), nor is there an obvious distribution of random matrices whose eigenvalues are associated with the points. 1 We note that there are recent universality results due to Kroó and Lubinsky on the asymptotics of Christoffel functions associated with multivariate orthogonal polynomials [6,7], but they do not apply in our setting.
Using the algorithm for sampling determinantal point processes associated with univariate orthogonal polynomials [10], which is trivially adapted to the orthogonal polynomials on the wedge, we can sample from this determinantal point process. We use this algorithm to calculate statistics of the points. In Fig. 1, we use the sampling algorithm in a Monte Carlo simulation to approximate the probability that no eigenvalue is present in a neighbourhood of three points for α = β = γ = 0. That is, we take 10,000 samples of a determinantal point process, and calculate the distance of the nearest point to z 0 , for z 0 equal to (1, 1), (0, 1), (0.5, 1) and (0.7, 1). The plots are of a complementary empirical cumulative distribution function of these samples. This gives an estimation of the probability that no eigenvalue is in a neighbourhood of z 0 . We have scaled the distributions so that the empirical variance is one: this ensures that the distributions tend to a limit as N becomes large, which is the regime where universality is present. In Fig. 2 we plot the same statistics but for samples from the unweighted Coulomb gas on the wedge, which has the distribution 1 Z N k< j λ k − λ j 2 for λ k supported on the wedge. As this is a Vandermonde determinant squared, it is also a determinantal point process with the basis arising from orthogonalized complexvalued polynomials 1, (x +iy), (x +iy) 2 , . . . [2]. We approximate this orthogonal basis using the modified Gram-Schmidt algorithm with the wedge inner product calculated via Clenshaw-Curtis quadrature. Again, this fits naturally into the sampling algorithm of [10], hence we can produce samples of this point process. What we observe is that, while our determinantal point process is not a Coulomb gas, it appears to be in the same universality class as the Coulomb gas away from the edge and corner, as the statistics follow the same distribution. This universality class matches that of the Gaussian Unitary Ensemble, as seen in Fig. 3 where we compare the three for N = 50.

Conclusion
We have introduced multivariate orthogonal polynomials on the wedge and boundary of a square for some natural choices of weights. We have also generated a complete orthogonal basis with respect to a suitable weight inside the square. We have looked at determinantal point process statistics and observed a relationship between the resulting statistics and Coulomb gases, suggesting that, away from the corner and edge, they are in the same universality class. One of the motivations for this work is to solve singular integral equations and evaluate their solutions on contours that have corners, in other words, to generalized the approach of [11]. Preliminary work in this direction is included in "Appendix B", which shows how the recurrence relationship that our polynomials satisfy can be used to evaluate Stieltjes transforms.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
For z off the contour, we can successfully and stably calculate the Stieltjes transform using a block-wise version of Olver's algorithm, which is equivalent to solving the 2n + 1 × 2n + 1 block-tridiagonal system q 0 = 1 C z k q k−1 + (A z k − z I )q k + B z k q k+1 = 0 fork = 1, 2, . . . , n − 1 q n = 0 0 where q 0 ∈ C 1 and q k ∈ C 2 for k = 1, 2, . . . , n. Then Olver's algorithm consists of performing Gaussian elimination adaptively until a convergence criteria is satisfied. For z on or near the contour, we no longer see quick decay in the Stieltjes transform (it is no longer a minimal solution to the recurrence), hence n must be prohibitively large. Instead, we adapt Olver's algorithm in a vein similar to Miller's algorithm to allow for a non-decaying tail. We do so by calculating two additional solutions q 1 , and q 2 (with the same block-sizes as before) satisfying: . Note that we need to choose n larger than necessary to avoid the errors in the tail, and there are unresolved numerical errors if z is close to the corner While this holds true for all n, we note that in practice we need to choose n bigger than the number of coefficients in order to observe numerical stability, see Fig. 4. We also find that there are still stability issues near the corner. Resolving these issues is ongoing research.