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.


Introduction
Let Ω be a wedge on the plane that consists of two line segments sharing a common end point. With respect to a positive measure dµ defined on Ω, we study orthogonal polynomials of two variables with respect to the bilinear form f, g = Ω f (x, y)g(x, y)dµ.
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.
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 also the zero set of a quadratic polynomial. 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.
For orthogonality on the wedge, we consider two cases. In the first case w 1 (x) = w 2 (x) = w(x), where w is any weight function on [0, 1]. In the second case w 1 (x) = x α (1−x) γ and w 2 (x) = σx β (1−x) γ are Jacobi weight functions and σ is a fixed positive constant. The space of orthogonal polynomials of degree n has dimension 2, just like that of spherical harmonics, and they satisfy the equation ∂ 1 ∂ 2 Y = 0. In both cases we construct an explicit basis of orthogonal polynomials and study the Fourier orthogonal expansions, which turns out to be related to the Fourier orthogonal expansions with respect to w 1 and w 2 on [0, 1], but the connection is subtle and non-trivial.
The second case is used to 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. Furthermore, the convergence of the Fourier orthogonal expansions can also be deduced in this fashion.
Spherical harmonics can be used to construct a basis of orthogonal polynomials for the weight function w( x 2 + y 2 ) on the unit disk. Likewise, we can 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 . However, unlike the unit disk, the orthogonal basis we constructed is no longer polynomials in x, y but in polynomials of x, y and s = max{|x|, |y|}.
The paper is organized as follows. In the next section we study the case when w 1 = w 2 = w. The case when w 1 and w 2 are the Jacobi weight functions are studied in Section 3, which requires substantially more effort. Orthogonal polynomials on the boundary of square are studied in Section 4, which are used to construct orthogonal functions on the square for weight functions that depend only on 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 Section 6, 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.
2.1. 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 Let I be the polynomial ideal of R[x, y] generated by (1 − x)(1 − y). If f ∈ I, then f, g w1,w2 = 0 for all g. The bilinear form defines an inner product on Π 2 n modulus I, or equivalently, on the quotient space R[x, y]/I.
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 xy ≡ 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), Let L 2 (w) denote the L 2 space with respect to w on [0, 1]. The Fourier orthogonal expansion of f ∈ L 2 (w) is defined by 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

2.2.
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. Theorem 2.2. Let w be a weight function on [0, 1] and let φw(x) := (1 − x) 2 w(x). Define P n (x, y) = p n (w; x) + p n (w; y) − p n (w; 1), n = 0, 1, 2, . . . , Then {P n , Q n } are two polynomials in H 2 n (w) and they are mutually orthogonal. Furthermore, (2.5) P n , P n w = 2h n (w) and Q n , Q n w = 2h n−1 (φw).
Proof. Since P n (x, 1) = P n (1, x) and Q n (x, 1) = −Q n (1, x), it follows that 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 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. Then Proof. By our definition, Similarly, 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 , which proves the stated result.

3.1.
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.
Let P (α,β) n denote the usual Jacobi polynomial of degree n defined on [−1, 1]. Then P by [12, is an orthogonal polynomial of degree n with respect to (1 − x) γ x α on [0, 1], I γ,α 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 (3.1), 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.
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 γ,α n +h γ,β n for m = n. Similarly, To derive the norm of Q n , Q n , we need to use c γ,α = (γ + 1) 2 /(α + γ + 2) 2 c γ+2,α . 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 A ∼ B if c 1 B ≤ A ≤ c 2 A, where c 1 and c 2 are constants that are independent of varying parameters in A and B. By (3.1) 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 Pn ∼ 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

Consequently, it follows that
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 lose of generality that the parallelogram is the square [−1, 1] 2 .
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 and define ψ : In particular, the norm of S n f − f is bounded by those of S α+δ1,β+δ2,γ m G δ1,δ2 − G δ1,δ2 as in Theorem 3.4.
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 (4.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 Since S α+δ1,β+δ2,γ n G δ1,δ2 • ψ(x 2 , y 2 ) → G δ1,δ2 (x, y) and f (x, y) = G 0,0 (x, y) + yG 0,1 (x, y) + yG 1,0 x, y + xyG 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 disk that are constructed using spherical harmonics. However, theses 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 4.2, whereas for n ≥ 3, they are constructed in Theorem 4.2. For n ≥ k, denote by P m,2n−2k the orthogonal polynomial of degree m with respect to t 2n−2k+1 w(t) on [0, 1] and with P 0,2n−2k (s) := 1. For n ≥ 0 and 0 ≤ k ≤ n, we define 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 (x, y) = 0 for all k, i, n, then f (x, y) = 0. Indeed, by the orthogonality of polynomials on the boundary, we see that modulus the ideal. Changing order of summation shows that This completes the proof.
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 gasses and random matrix theory is no longer obvious: the interaction of the points is not Coulomb (that is, it can not 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 Figure 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 Figure 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.

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.
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 − zI)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 S Ω [P k w](z) ≈ S Ω [P 0 w](z)q k 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: q j 0 = 1 C z k q j k−1 + (A z k − zI)q j k + B z k q j k+1 = 0 for k = 1, 2, . . . , n − 1 q 1 n = 1 0 and q 2 n = 0 1 .
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 Figure B. We also find that there are still stability issues near the corner. Resolving these issues is ongoing research.