Threefold symmetric Hahn-classical multiple orthogonal polynomials
Ana F. Loureiro Walter Van Assche


A. F. Loureiro and W. Van AsscheThreefold symmetric Hahn-classical multiple orthogonal polynomials 













School of Mathematics, Statistics and Actuarial Science (SMSAS), Sibson Building, 
University of Kent, Canterbury, Kent CT2 7FS (U.K.)

a.loureiro@kent.ac.uk (corresponding author)



Department of Mathematics, KU Leuven, Celestijnenlaan 200B box 2400, BE-3001 Leuven, Belgium

walter.vanassche@kuleuven.be

THREEFOLD SYMMETRIC HAHN-CLASSICAL MULTIPLE ORTHOGONAL POLYNOMIALS 
    WALTER VAN ASSCHE
    December 30, 2023
===================================================================











We characterize all the multiple orthogonal threefold symmetric polynomial sequences whose sequence of derivatives is also multiple orthogonal. Such a property is commonly called the Hahn property and it is an extension of the concept of classical polynomials to the context of multiple orthogonality. The emphasis is on the polynomials whose indices lie on the step line, also known as 2-orthogonal polynomials. We explain the relation of the asymptotic behavior of the recurrence coefficients to that of the largest zero (in absolute value) of the polynomial set. We provide a full characterization of the Hahn-classical orthogonality measures supported on a 3-star in the complex plane containing all the zeros of the polynomials. There are essentially three distinct families, one of them 2-orthogonal with respect to two confluent functions of the second kind. 
This paper complements earlier research of Douak and Maroni.  




Mathematics Subject Classification 2000: 33C45; 42C05




§ INTRODUCTION AND MOTIVATION


In this paper we investigate and characterize all the multiple orthogonal polynomials of type II that are threefold symmetric and are such that the polynomial sequence of its derivatives is also a multiple orthogonal sequence of type II. 
Within the standard orthogonality context there are only four families satisfying this property, commonly referred to as the Hahn property, and they are the Hermite, Laguerre, Jacobi and Bessel,  collectively known as the classical orthogonal polynomials. These four families of polynomials also share a number of analytic and algebraic properties. Several studies are dedicated to extensions of those properties to the context of multiple orthogonality. However, those extensions give rise to completely different sequences of multiple orthogonal polynomials. For the usual orthogonal polynomials the Hahn property is equivalent with the Bochner characterization (polynomials satisfy a second order differential equation of Sturm-Liouville type) and the existence of a Rodrigues formula. This is no longer true for multiple orthogonal polynomials since there
are families of multiple orthogonal polynomials with a Rodrigues type formula <cit.> that do not satisfy the Hahn property, and various examples of multiple orthogonal polynomials satisfy a higher-order linear recurrence
relation, which is not of Sturm-Liouville type <cit.>. Hence characterizations for multiple orthogonal polynomials based on the Hahn property, the Bochner property or the Rodrigues formula give different families of
polynomials. 

A sequence of monic polynomials { P_n}_n⩾ 0 with P_n =n is orthogonal with respect to a Borel measure μ whenever 

    ∫ P_n(x) x^k μ(x) = 0  if k=0,1,… ,n-1,  n≥ 0

and 

    ∫ P_n(x) x^n μ(x) ≠ 0  for  n⩾ 0.

Without loss of generality, often we normalize μ so that it is a probability measure. Obviously, an  forms a basis of the vector space of polynomials 𝒫. The measures described above can be represented via a linear functional , defined on 𝒫^', the dual space of 𝒫, and it is understood that the action of  over a polynomial f corresponds to ∫ f(x) μ(x). Throughout, we denote this action as 

    ⟨,f(x)⟩∫ f(x) μ(x).

The derivative of a function f is denoted by f' or / xf(x). 
Properties on 𝒫^', such as differentiation or multiplication by a polynomial, can be defined by duality. A detailed explanation can be found in <cit.> <cit.>. In particular, given g∈𝒫 and a functional ∈𝒫^', we define  

    ⟨ g(x) , f(x) ⟩⟨ , g(x)f(x) ⟩  and  ⟨ ', f(x) ⟩ - ⟨ , f' (x)⟩

for any polynomial f. 

From the definition it is straightforward that an  satisfies a second order recurrence relation 

    P_n+1(x) = (x-β_n) P_n(x) - γ_n P_n-1(x),  n⩾ 1,

with P_0(x) =1, P_1(x)=x-β_0 and γ_n≠ 0 for all integers n⩾ 1. This relation is often called the three-term recurrence relation, but we will avoid this terminology as we will be dealing with three-term recurrence relations which are of higher order. There is an important converse of this connection between orthogonal polynomials and second order recurrence relations, known as the Shohat-Favard theorem or spectral theorem for orthogonal polynomials. It states that any sequence of monic polynomials { P_n}_n⩾ 0 with P_n =n, satisfying the recurrence relation (<ref>) with γ_n≠ 0 for all n⩾1 and initial conditions P_-1(x) =0 and P_0(x)=1, is always an  with respect to some measure μ and, if, in addition, β_n∈ℝ and γ_n+1>0 for all n⩾0, then μ is a positive measure on the real line. So, basically, the orthogonality conditions and the second order recurrence relations are two equivalent ways to characterize an . 

Multiple orthogonal polynomials of type II correspond to a sequence of polynomials of a single variable that satisfy multiple orthogonality conditions with respect to r>1 measures. With the multi-index n⃗=(n_1,…,n_r) ∈ℕ^r, type II multiple orthogonal polynomials correspond to a (multi-index) sequence of monic polynomials P_n⃗(x) of degree |n⃗| = n_1+…+n_r for which there is a vector of measures 
(μ_0, …, μ_r-1) such that   

    ∫ x^k P_n⃗(x) μ_j(x) = 0 ,    0⩽ k ⩽ n_j-1

holds for every j=0,1,…, r-1.  By setting r=1, we recover the usual orthogonal polynomials. Obviously, (<ref>) amounts to the same as saying there exists a vector of r linear functionals (_0,…, _r-1) such that 

    ⟨_j, x^k P_n⃗(x) ⟩ = 0 ,    0⩽ k ⩽ n_j-1.
 

This polynomial P_n⃗ set may not exist or is not unique unless we impose some extra conditions on the r measures μ_0,…,μ_r-1, but
under appropriate conditions the polynomial set satisfies a system of r recurrence relations, relating P_n⃗ with its nearest neighbors 
P_n⃗+e⃗_k and P_n⃗-e⃗_j, where e⃗_k consists of the r-dimensional unit vector with zero entries except for the kth component which is 1 (see <cit.> for further details). The focus of the present work is  on the polynomials whose multi-indices lie on the step-line near the diagonal n⃗=(n,n,…,n), and are defined by  

    P_rn+j(x) = P_n⃗+s⃗_j(x) ,   with   s⃗_j= ∑_ℓ=0^j e⃗_ℓ,      j=0,1,…, r-1,
 
that is, 

    P_rn+j(x) = P_( n+1,…,n+1_j, n,…, n_r-j)(x)
    	,      j=0,1,…, r-1.
 

This step-line polynomial sequence {P_n}_n⩾0 satisfies a (r+1)th order recurrence relation (involving up to (r+2) consecutive terms) and it can be written as 

    [ xP_n(x) = P_n+1(x) + ∑_j=0^rζ_n,j P_n-j(x) ,      n⩾ r,;                            P_j(x)=x^j,       j=0,…, r , ]

where ζ_n,r≠ 0 for all n⩾ r. Such a step-line polynomial sequence {P_n}_n⩾0 corresponds to so-called d-orthogonal polynomials (with d=r), see <cit.> and Figure <ref> for the case of r=2.  So, basically, if there exists a vector of linear functionals (_0,…, _d-1) and a polynomial sequence {P_n}_n⩾0 satisfying 

    ⟨_j, x^k P_n(x) ⟩ = 0    and  ⟨_j, x^n P_dn+j(x) ⟩≠ 0  , 
    	   0 ⩽ k ⩽⌊n-j/d⌋,  n⩾ 0,

then the polynomials {P_n}_n⩾0 are related by (<ref>). There is a converse result, which is a natural generalization of the Shohat-Favard theorem, in the sense that if a polynomial sequence {P_n}_n⩾0 satisfies (<ref>) with ζ_n,d≠ 0 for all n⩾ d, then there is a vector of r linear functionals (_0, …, _d-1) with respect to which {P_n}_n⩾0 is d-orthogonal. Such a vector of linear functionals is not unique, but we can consider its components to be the first d elements of the dual sequence {u_n}_n⩾0 associated to {P_n}_n⩾0, which always exist and which are defined by 

    ⟨ u_n,P_m ⟩δ_n,m	,    m,n⩾ 0,

where δ_n,m represents the Kronecker symbol. The remaining elements of the dual sequence of a d- can be generated from the first ones: for each pair of integers (n,j) there exist polynomials q_n,j,ν(x) such that <cit.> 

    u_dn+j =  ∑_ν=0^d-1 q_n,j,ν(x) u_ν,   for   0⩽ j⩽ d-1,

where q_n,j,j(x) =n, q_n,j,ν(x) ⩽ n for 0⩽ν⩽ j-1 if 1⩽ j ⩽ d-1 and q_n,j,ν(x) ⩽ n -1 for j+1⩽ν⩽ d-1 if 0⩽ j⩽ d-2. Unlike the standard orthogonality (case d=1), this result only provides structural properties for the dual sequence and in practical terms it can be used in that sense. 

In this paper we mainly deal with 2-s, but our results and ideas can be extended to the d-orthogonal case. Here, we provide a characterization of all the 2-s  that are threefold symmetric and possess the Hahn property in the context of 2-orthogonality, i.e., the sequence of monic derivatives is also 2-orthogonal. The notion and properties of threefold symmetric 2-s are revised and discussed in Section <ref>, while we also bring a new result (Theorem <ref>) relating the asymptotic behavior of the recurrence coefficient to the asymptotic behavior of the absolute value of the largest zero. Section <ref> is dedicated to a detailed characterization of all the threefold symmetric 2-orthogonal polynomials satisfying the Hahn property. We bring together old and new results in a self-contained and complete analysis to fully characterize this type of sequences in terms of their explicit recurrence relations, weight functions and a differential equation of third order. In that regard, we complete the study initiated by Douak and Maroni in <cit.> by providing sufficient conditions and taking the study based only on the properties of the zeros in combination with the weights. 
We end up this coherent description with a detailed analysis of the four distinct families of threefold symmetric 2-orthogonal polynomials satisfying the Hahn property. The first case, treated in Section <ref>  corresponds to polynomials with the Appell property whose 2-orthogonality weights are supported on the three-starlike set (as in Fig.<ref>) represented via the Airy function and its derivative. These polynomials have been studied in 
<cit.>, but the explicit representation of the weights supported on a set containing all the zeros has not been considered there. In Sections <ref> and <ref> we study the second and third cases where the weights are represented via the confluent hypergeometric Kummer function of 2nd kind. These are, to the best of our knowledge, new. The fourth case is treated in Section <ref>, where the orthogonality weights are expressed via hypergeometric functions and depend on two parameters. For special choices of those parameters we recover some known polynomial families. However, the existent literature on the subject has not taken into account the explicit representation of the weights supported on a set containing all the zeros and this has been accomplished here for all the threefold symmetric polynomials satisfying the Hahn property within the context of 2-orthogonality. 

The orthogonality (<ref>) of these polynomials corresponds to an non-hermitian inner product (no conjugation is used) on three-starlike sets. The Airy function, the confluent hypergeometric function and the hypergeometric function are analytic functions, so the paths of integration can
be deformed away from the three-starlike sets without loosing orthogonality. However, the choice of the three-starlike sets is convenient because
of the positivity of the measures and the location of the zeros.

It turns out that the orthogonality measures under analysis are solutions to a second order differential equation and, as such, the recurrence coefficients for the whole set of the corresponding multiple orthogonal polynomials of type II can only be obtained algorithmically via the nearest-neighbor algorithm explained in <cit.> <cit.>. 
Some of the d-orthogonal polynomials that we found appear in the theory of random matrices, in particular in the investigation of singular values
of products of Ginibre matrices, which uses multiple orthogonal polynomials with weight functions expressed in terms of Meijer G-functions
<cit.>. For r=2 these
Meijer G-functions are hypergeometric or confluent hypergeometric functions.




§ THREEFOLD 2-


A sequence of monic polynomials {P_n}_n⩾ 0 (with P_n=n) is symmetric whenever P_n(-x) = (-1)^n P_n(x) for all n⩾0. This means that all even degree polynomials are even functions while odd degree polynomials are odd functions. Hermite and Gegenbauer polynomials are examples of symmetric polynomial sets, which also happen to be the only classical orthogonal polynomials that are symmetric. Many other examples of symmetric orthogonal polynomial sequences are around in the literature. 

The notion of symmetry of polynomial sequences has been extended and commonly referred to as d-symmetric in works by Maroni <cit.>  and followed by Ben Cheikh and his collaborators <cit.>. The case d=1 would correspond to the usual symmetric case. We believe the name is misleading, as will soon become apparent, and therefore we call it differently as it pictures better the nature of the problem. 


A polynomial sequence {P_n}_n⩾ 0 is 𝐦-fold symmetric if 

    P_n(ω x) = ω^n  P_n(x),   for any   n⩾ 0.

where 
ω= ^2 i  π/m.

By induction, property (<ref>) corresponds to 

    P_n(ω^k x) = ω^nk  P_n(x),   for any   n⩾ 0   and  k=1,2,…, m-1.
 

So, an m-fold symmetric sequence is a (m-1)-symmetric sequence in  <cit.>. 
(The symmetric sequences of Hermite and Gegenbauer polynomials are examples of twofold symmetric polynomials.) In particular, a threefold symmetric sequence {P_n}_n⩾ 0, is such that 

    P_n(ω x) =  ω^n P_n(x)   and  
    	P_n(ω^2 x) = ω^2n P_n(x)    with  ω =^2 i  π/3, n⩾0,

which corresponds to say that there exist three sequences {P_n^[j]}_n⩾ 0 with j∈{0,1,2} such that 

    P_3n+j(x) = x^j P_n^[j](x^3),    j=0,1,2.

Throughout, we will refer to {P_n^[j]}_n⩾ 0 as the diagonal components of the cubic decomposition of the threefold symmetric sequence {P_n}_n⩾ 0, which is line with the terminology adopted in a more general cubic decomposition framework in <cit.>. 

In the case of a 2- { P_n}_n⩾ 0 with respect to a vector linear functional u = (u_0,u_1), we have, as discussed in Section <ref>, 

    ⟨ u_0 , x^m P_n ⟩ =
    		{[         0       for   n⩾ 2m+1; N_0(n)≠ 0       for     n= 2m;           ]. 
          ⟨ u_1 , x^m P_n ⟩  =
    		{[         0       for   n⩾ 2m+2; N_1(n)≠ 0       for  n= 2m+1,;           ].

and there exists a set of coefficients {(β_n, α_n,γ_n)}_n⩾0 such that { P_n}_n⩾ 0 satisfies a third order recurrence relation (see <cit.>) 

    P_n+1(x) = (x-β_n) P_n(x) - α_n P_n-1(x) - γ_n-1 P_n-2(x)

with P_-2(x) = P_-1(x) =0 and P_0(x)=1. Straightforwardly from the definition, one has 

    γ_2n+1 =⟨ u_0, x^n+1 P_2n+2⟩/⟨ u_0, x^n P_2n⟩≠ 0,
    	  γ_2n+2 = ⟨ u_1, x^n+1 P_2n+3⟩/⟨ u_1, x^n P_2n+1⟩≠ 0 ,   n⩾ 0,
 
or, equivalently, 

    ⟨ u_0, x^n+1 P_2n+2⟩ = ∏_k=0^n γ_2k+1  and   ⟨ u_1, x^n+1 P_2n+3⟩ = ∏_k=0^n γ_2k+2 ,  for    n⩾ 0.

Whenever a 2- is threefold symmetric, the recurrence relation (<ref>) reduces to a three-term relation, where the β- and α-coefficients all vanish. For this type of 2-s  more can be said.  

<cit.> Let {P_n}_n⩾ 0  be a 2-orthogonal polynomial sequence with respect to the linear functional U=(u_0, u_1) satisfying (<ref>)-(<ref>). The following statements are equivalent: 


  (a) The sequence {P_n}_n⩾ 0  is threefold symmetric. 

  (b) The linear functional is threefold symmetric, that is, 

    (u_ν)_3n+j =0 ,   for  ν =0,1   and   j=1,2   with   j≠ν,

where (u)_m ⟨ u, x^m⟩ for m⩾ 0. 

  (c) The sequence {P_n}_n⩾ 0  satisfies the third order recurrence relation 

    P_n+1(x) = x P_n(x)  - γ_n-1 P_n-2(x) ,      n⩾ 2,

with P_0(x) = 1,  P_1(x)= x and P_2(x) = x^2. 



Each of the components of the cubic decomposition are also 2-s. 

<cit.> Let {P_n}_n⩾ 0 be a threefold symmetric 2-OPS. The three polynomial sequences {P_n^[j]}_n⩾ 0 (with j=0,1,2) in the cubic decomposition of {P_n}_n⩾ 0 described in  (<ref>) are 2-orthogonal polynomial sequences satisfying: 

    P_n+1^[j](x) = (x-β_n^[j]) P_n^[j](x)  -α_n^[j] P_n-1^[j](x) - γ_n-1^[j] P_n-2^[j](x) ,

where 

    [                              β_n^[j] = γ_3n-1+j+γ_3n+j+γ_3n+1+j,      n⩾ 0,; α_n^[j] = γ_3n-2+jγ_3n+j+γ_3n-1+jγ_3n-3+j+γ_3n-2+jγ_3n-1 +j
    			,      n⩾ 1,;                             γ_n^[j] = γ_3n-2+jγ_3n+jγ_3n+2+j≠ 0,      n⩾ 2.;                                                                             ]

Moreover, {P_n^[j]}_n⩾ 0 is 2-orthogonal with respect to the vector functional U^[j]=(u_0^[j],u_1^[j]) with

    u_ν^[j] = σ_3 (x^j u_3ν+j)   for each  j=0,1,2   and  ν=0,1.
 
where {u_n}_n⩾0 represents the dual sequence of {P_n}_n⩾ 0 and σ_3:𝒫' ⟶𝒫' is the operator defined by ⟨σ_3 (v), f(x) ⟩⟨ v, f(x^3) ⟩ for any v∈𝒫' and f∈𝒫.


The orthogonality measures are supported on a starlike set with three rays. 



<cit.> If γ_n>0 for n⩾ 1 in (<ref>), then there exists a vector of linear functionals  (u_0,u_1) such that the polynomials P_n defined by (<ref>) satisfy the 2-orthogonal relations (<ref>)-(<ref>). Moreover, the vector of linear functionals  (u_0,u_1) satisfies (<ref>) and there exist a vector of two measures (μ_0,μ_1) such that 

    ⟨ u_0,f(x) ⟩ = ∫_S f(x) dμ_0(x) 
          ⟨ u_1,f(x) ⟩ = ∫_S f(x) dμ_1(x)

where S represents the starlike set 

    S⋃_k=0^2 Γ_k    with     Γ_k = [0,^2π i k/3∞) ,

and the measures have a common support which is a subset of S and are invariant under rotations of 2π /3. 





Regarding the behavior of the zeros of any threefold symmetric 2-{P_n}_n⩾0, it was shown in <cit.> that between two nonzero consecutive roots of P_n+2 there is exactly one root of P_n and one root of P_n+1 and all those roots lie on the starlike set S. 


Let {P_n}_n⩾0 be a 2-OPS satisfying (<ref>). If γ_n >0, then the following statements hold: 


  (a) If x is a zero of P_3n+j, then ω^j x are also zeros of P_3n+j with ω=^2π i/3.

  (b) 0 is a zero of P_3n+j of multiplicity j when j=1,2.

  (c) P_3n+j has n distinct positive real zeros {x_n,k^(j)}_k=1^n with 0<x_n,1^(j)<… <x_n,n^(j).  

  (d) Between two real zeros of P_3n+j+3 there exist only one zero of  P_3n+j+2 and only one zero of P_3n+j+1, that is, 
x_n,k^(j+2)<x_n,k+1^(j)<x_n,k+1^(j+1)<x_n,k+1^(j+2). 


 The result is a consequence of <cit.> for the case d=2. 



So, P_3n+j with j∈{0,1,2} has n zeros (x_n,1^(j), … ,x_n,n^(j))  on the positive real line and all the other zeros are obtained by rotations of 2π/3 of (x_n,1^(j), … ,x_n,n^(j)) and 0 is a single zero for P_3n+1 and a double zero for P_3n+2. The connection between the asymptotic behavior of the γ-recurrence coefficients in (<ref>)  and the 
upper bound for the largest zero x_n,n^(j) is discussed in <cit.>, but for bounded recurrence coefficients. Here, we extend that discussion, by embracing the cases where the recurrence coefficients are unbounded with different asymptotic behavior for even and odd order indices, which will be instrumental in Section <ref>.

 If γ_n in (<ref>) are positive and, additionally, γ_2n = c_0 n^α +o(n^α) and γ_2n+1 = c_1 n^α +o(n^α) for large n, with min{c_0,c_1}⩾ 0, c=max{c_0,c_1}>0 and α⩾ 0, then largest zero in absolute value |x_n,n| behaves as 

    |x_n,n| ⩽3/2^2/3 c^1/3 n^α/3 + o(n^α/3),      n⩾ 1.


 Consider the Hessenberg matrix 

    𝐇_n = ( 
    		[     0     1     0     ⋯     0     0     0     0     0;     0     0     1     0     ⋯     0     0     0     0;   γ_1     0     0     1     0     ⋯     0     0     0;     0   γ_2     0     0     1     0     ⋯     0     0;                 ⋱                 ⋱                  ;                       ⋱                 ⋱            ;     0     0     ⋯     0 γ_n-4     0     0     1     0;     0     0     0     ⋯     0 γ_n-3     0     0     1;     0     0     0     0     ⋯     0 γ_n-2     0     0;       ])
 
so that the recurrence relation (<ref>) can be expressed as 

    𝐇_n ([   P_0(x);   P_1(x);        ⋮; P_n-1(x) ])
    	= x ([   P_0(x);   P_1(x);        ⋮; P_n-1(x) ])
    	- P_n(x) ([ 0; 0; ⋮; 1 ])

and each zero of P_n is an eigenvalue of the matrix 𝐇_n. 
The spectral radius of the matrix 𝐇_n, 

    ρ(𝐇_n) = max{ |λ| :  λ is an eigenvalue of 𝐇_n},
 
is bounded from above by  || 𝐇_n || where || · || denotes a matrix norm (see <cit.>). 
We take the matrix norm 

    || 𝐇_n ||_S =|| S^-1𝐇_n S ||_∞ = max_1⩽ i⩽ n{∑_j=1^n |(S^-1𝐇_n S)_i,j| },

where S corresponds to a non-singular matrix and (S^-1𝐇_n S)_i,j denotes the ith row and jth column entry of the product matrix S^-1𝐇_n S. In particular if S is an invertible diagonal matrix   S = diag(d_1,…, d_k, …, d_n), then 

    || 𝐇_n ||_S 
    
    	= max{d_2/d_1, d_3/d_2,    d_4 + d_1 γ_1/d_3, …, d_k + d_k-3γ_k-3/d_k-1 , …, 
    	d_n + d_n-3γ_n-3/d_n-1,   d_n-2γ_n-2/d_n}.

Setting d_k= d^k (k!)^α/3≠ 0, for some  positive constant d, gives 

    || 𝐇_n ||_S ⩽ 2^α/3(d+c/d^2) n^α/3 + o(n^α/3)   as  n→ +∞.
 
The choice of d=(2c)^1/3 gives a minimum to (d+c/d^2), so that  

    || 𝐇_n ||_S ⩽3/4^1/3( c  n^α)^1/3+ o(n^α/3)   as  n→ +∞,
 
which implies the result. 


To summarize,  Proposition <ref> combined with Lemma <ref> allow us to conclude that each P_n^[j] has exactly n real zeros {x_n,k^[j]}_k=1^n with 0<x_n,1^[j]<… <x_n,n^[j]. Moreover,  between two consecutive zeros of P_n^[j+2] there is exactly one zero of P_n^[j] and another of P_n^[j+1]: x_n,k^[j+2]<x_n,k+1^[j]<x_n,k+1^[j+1]<x_n,k+1^[j+2]. Taking into consideration the asymptotic behavior for the largest zero |x_n,n| of P_n described in Theorem <ref>,  we then conclude that 

    x_n,k^[j]⩽27/4 c  n^α + o(n^α).

Therefore, for each j∈{0,1,2}, all the zeros x_n,k^(j) of P_3n+j distinct from 0 correspond to the cubic roots of x_n,k^[j]>0, with k∈{1,…, n}, and they all lie on the starlike set S and within the disc centred at the origin with radius b= 27/4 c  n^α + o(n^α). 

Multiple orthogonal polynomials on starlike sets received attention in recent years, under different frameworks. This includes the asymptotic behavior of polynomial sequences generated by recurrence relations of the type (<ref>) when further assumptions are taken regarding specific behaviour for the γ-coefficient <cit.> or for certain type of the 2-orthogonality measures in <cit.> and <cit.>. Faber polynomials associated with hypocycloidal domains and stars have also been studied in <cit.>. 

Here, we describe all the threefold 2-s that are classical in Hahn's sense, and our study includes the  representations for the measures supported on a set containing all the zeros of the polynomial sequence. From Theorem <ref> and Proposition <ref>, the support lies on  a starlike set, that, according to  Theorem <ref> is bounded if the γ-coefficients are bounded, and unbounded otherwise. 





§ THREEFOLD SYMMETRIC HAHN-CLASSICAL 2-
 

The classical orthogonal polynomial sequences of Hermite, Laguerre, Jacobi and Bessel collectively satisfy the so-called Hahn property:   
the sequence of its derivatives is again an orthogonal polynomial sequence. In the context of 2-orthogonality, this algebraic property is portrayed as follows: 
 A monic 2- { P_n}_n⩾ 0 is "2-Hahn-classical"  when the sequence of its derivatives { Q_n}_n⩾ 0, with Q_n(x) 1/n+1 P_n+1'(x) is also a 2-. 


The study of this type of 2-orthogonal sequences was initiated in the works by  Douak and Maroni <cit.>-<cit.>. In those works (as well as in <cit.>) several properties of these polynomials were given, with the main pillars of the study being the structural properties, including the recurrence relations, satisfied by the polynomials. 
All in all, those studies encompassed the analysis of a nonlinear system of equations fulfilled by the recurrence coefficients. Douak and Maroni treated some special solutions to that system of equations, bringing to light several examples of these threefold symmetric "2-Hahn-classical" polynomials: see <cit.>. However, for those cases the support of the corresponding orthogonality measures that they found consisted of the positive real axis, which does not contain all the zeros. 
 Here, we base our analysis on the properties of the orthogonality measures and deduce the properties of the recurrence coefficients. We incorporate the works by Douak and Maroni and go beyond that by fully describing all the threefold "2-Hahn-classical" polynomials and bringing up explicitly the orthogonality measures along with the asymptotic behavior of the largest zero in absolute value as well as a Bochner type result for the polynomials (i.e., characterizing these polynomials via a third order differential equation). 

We start by observing that  the threefold symmetry of { P_n}_n⩾ 0 readily implies the threefold symmetry of { Q_n(x) 1/n+1 P_n+1'(x)}_n⩾ 0. This is a straightforward consequence of Definition <ref>, as it suffices to take single differentiation of relation (<ref>). Such property is valid for any polynomial sequence, regardless any orthogonality properties. 

Regarding threefold symmetric 2-s possessing Hahn's property, more can be said. The next result summarizes a characterization of the orthogonality measures along side with the  recurrence relations of the original sequence {P_n(x)}_n⩾ 0 and the sequence of derivatives {Q_n(x)1/n+1P'_n+1(x)}_n⩾ 0. This characterization can be found in the works <cit.>. Nonetheless, we revisit these results for a matter of completion while bringing different approaches to the original proofs, highlighting that all the structural properties for the polynomials can be derived from the corresponding 2-orthogonality measures. 


 Let {P_n}_n⩾ 0 be a threefold symmetric 2- for (u_0,u_1) satisfying the recurrence relation (<ref>) and let {Q_n(x)1/n+1P'_n+1(x)}_n⩾ 0. The following statements are equivalent: 
 


  * {Q_n(x)}_n⩾ 0 is a threefold symmetric 2-, satisfying the third-order recurrence relation 

    Q_n+1(x) = x Q_n(x) -    γ_n-1 Q_n-2(x),

with initial conditions Q_k(x)=x^k for k∈{0,1,2}.


  * The vector functional (u_0,u_1) satisfies the matrix differential equation 

 
    ( Φ[ [ u_0; u_1 ]]) '
    	+ Ψ[ [ u_0; u_1 ]]
    	= 
    	[ [ 0; 0 ]],

where 

    Φ = [ 
    [              ϑ_1        (1-ϑ_1) x; 2γ_1(1-ϑ_2)  x^2          2ϑ_2 -1;                  ]]
      and  Ψ 
    	= [ [       0       1; 2/γ_1 x       0 ]]

for some constants ϑ_1 and ϑ_2  such that

    ϑ_1,ϑ_2 ≠n-1/n ,   for all   n⩾ 1.




  * there are coefficients ϑ_1,  ϑ_2≠n-1/n, such that 𝐔=(u_0,u_1) satisfies 
 
    (  ϕ(x) u_0)” 
    		+ ( 2 /γ _1(ϑ_2 + ϑ _1-2)x^2u_0)' 
    	+ 2/γ _1(ϑ _1-2) x u_0
    	=0

 and 

	
	
	(ϑ_1-2) (2 ϑ_2-1) u_1 

	=	ϕ(x)u_0' 
			- 2 /γ_1 (ϑ_1-1) (2 ϑ_2-3)  x^2 u_0 , 
		
			      if ϑ_1≠2,
	

	
	x  u_1^' = 2 u_0^' ,  
			    if ϑ_1=2, 

where 

    ϕ(x) =  ϑ _1 (2 ϑ _2-1)
    		-2/γ _1(ϑ _1-1) 
    			(ϑ _2-1) x^3 .



  * There exists a sequence of numbers {γ_n+1}_n⩾ 0 such that 

    P_n+3(x) 	=  Q_n+3(x) 
    				+ ( (n+1) γ_n+2 - (n+3)γ_n+1) Q_n(x),

with initial conditions 
P_0(x) = Q_0(x)=1  ,   
		P_1(x) = Q_1(x)  = x   and  
	P_2(x) = Q_2(x) =x^2    .


  We prove that (a) ⇒ (d) ⇒ (b) ⇔ (c) ⇒ (a). 

In order to see (a) implies (d), we differentiate the recurrence relation satisfied by P_n to then replace the differentiated terms by its definition of Q_n. A substitution of the term x Q_n(x) by the expression provided in (<ref>) finally gives (<ref>). 

Now we prove that (d) implies (b). 
Any linear functional w in 𝒫^' can be written as 

    w = ∑_k=0^∞⟨ w, P_k ⟩ u_k.
 
Based on this, the relation Q_n(x)1/n+1P_n+1'(x) gives a differential relation between the dual sequence {v_n}_n⩾ 0 of  {Q_n}_n⩾ 0 in terms of  {u_n}_n⩾ 0 (the dual sequence of  {P_n}_n⩾ 0), which is 

    v_n' = -(n+1) u_n+1,      n⩾ 0,

whilst (<ref>) leads to 

    v_n = u_n + ( (n+1) γ_n+2-(n+3) γ_n+1) u_n+3 ,      n⩾ 0.

The choice of n=0 and n=1 in both (<ref>) and (<ref>) respectively gives

    [ [ v_0'; v_1' ]]
    	= - 
    	[ [     u_1; 2   u_2 ]]
    	  and  [ [ v_0; v_1 ]]
    	=  
    	[ [  u_0 + (  γ_2-3 γ_1) u_3; u_1 + ( 2 γ_3-4 γ_2) u_4 ]].

As explained in <cit.> (and, alternatively, in <cit.>, <cit.> and <cit.>), the elements of the dual sequence {u_n}_n⩾0  of the 2- { P_n}_n⩾ 0  can be written as 

    [ u_2n 	= E_n (x) u_0 + a_n-1(x) u_1 ,; u_2n+1  = b_n (x) u_0 + F_n(x) u_1 , ]
 
where E_n and F_n are polynomials of degree n, while a_n and b_n are polynomials of degree less than or equal to n under the assumption that a_-1(x)=0, E_0=1, F_0=1 and b_0=0. The recurrence relations (<ref>)  fulfilled by  { P_n}_n⩾ 0 yield (see <cit.>)

	   γ_2n+2 F_n+1 (x)- x F_n(x) = -a_n-1(x),            γ_2n+3 a_n+1 (x)- x a_n(x) = -F_n(x),  

	   γ_2n+2 b_n+1 (x)- x b_n(x) = -E_n(x)    and          γ_2n+3 E_n+2 (x)- x E_n+1(x) = -b_n(x),  

with initial conditions a_0 (x) =0 and  
E_1(x) = 1/γ_1 x. In particular, we obtain: 

    b_1(x)
    		= - 1/γ_2, 
    	  
    	F_1(x) = 1/γ_2x,   
    	E_2(x) = 1/γ_1γ_3  x^2 ,   and  
    	a_1(x) = -1/γ_3,

so that 

    u_2 = 1/γ_1 x   u_0  ,    
    	u_3 = - 1/γ_2 u_0 + 1/γ_2x   u_1  
    	  and   
    	u_4 = 1/γ_1γ_3  x^2   u_0 -1/γ_3  u_1.

Consequently, the first identity in (<ref>) reads as 


    [ [ v_0'; v_1' ]]
    	= - 
    	Ψ
    	[ [ u_0; u_1 ]]
    	  with   Ψ = [ [       0       1; 2/γ_1 x       0 ]].

Using (<ref>) in the second identity in (<ref>) leads to 

    [ [ v_0; v_1 ]]
    	=Φ[ [ u_0; u_1 ]]
    	  with   Φ 		
    	= [ [ ϕ_0,0 ϕ_0,1; ϕ_1,0 ϕ_1,1 ]],


where (c.f. <cit.>)

    ϕ_00(x) 
    			= 3γ_1/γ_2 , 
    	  ϕ_01(x) 
    			=  (1 -   3γ_1/γ_2) x , 
          ϕ_10(x) 
    			=   2(1/γ_1-  2γ_2/γ_3) x^2 
    	  and  ϕ_11(x) 
    			= -1+  4γ_2/γ_3.

We set ϑ_1=3γ_1/γ_2≠ 0 and ϑ_2=2γ_2/γ_3≠ 0, and now (<ref>)-(<ref>) follows after differentiating both sides of the equation in (<ref>) and then compare with (<ref>), which is a system of two functional equations in (u_0,u_1):

    ϑ_1 u_0^' + (1- ϑ_1) x  u_1^' + (2-ϑ_1) u_1 = 0 
    2/γ_1 (1-ϑ_2) x^2  u_0^' + 2/γ_1 (3-2ϑ_2) x  u_0
    			- (1-2ϑ_2)u_1^'=0.

The action of each of these equations over the monomials gives 

    n ϑ_1 (u_0)_n-1 + ( (n-1) - n ϑ_1 ) (u_1)_n = 0, 
    1γ_1( (n-1) - n ϑ_2 )  (u_0)_n+1 
    			+ (2ϑ_2 -1)n (u_1)_n-1 = 0 ,       n⩾ 0,
 
under the assumption that (u_0)_-1=(u_1)_-1=0. All the moments are well defined provided that 

    ( (n-1) - n ϑ_1 ) ( (n-1) - n ϑ_2 ) ≠ 0 ,      n⩾ 0,
 
which corresponds to condition (<ref>). Under this assumption, the threefold symmetry implies that  

    (u_0)_n = (u_1)_n+1 =0   for   n ≢0  3,
 
whilst for n≡ 0  3, we have (u_0)_n (u_1)_n+1≠ 0, and these are recursively defined by 

    2/γ_1( n+1 -( n+2) ϑ_2 ) ( n -( n+1) ϑ_1 ) (u_0)_n+3 
    	= ϑ_1(2 ϑ_2-1)(n+1)(n+2) (u_0)_n , 
       ( n+1 -( n+2) ϑ_2 ) ( n+3 -( n+4) ϑ_1 )(u_1)_n+4  
    		= ϑ_1(2ϑ_2 -1)(n+2)(n+4)   (u_1)_n+1 ,

where  (u_0)_0=(u_1)_1=1. 

Concerning the proof of  (b) ⇒ (c), we start by noting that the system of equations  (<ref>) can be written as in (<ref>) which reads as 

    C [ [ u_1';  u_1 ]]
    	= -D [ [ u_0';  u_0 ]],

with  

    C= [ [ (1-ϑ _1) x     2-ϑ _1;    2ϑ _2-1          0 ]]
    	    and      
    	D= [ [              ϑ_1                0; 2/γ_1 (1-ϑ_2)x^2  2/γ_1 (3-2ϑ_2)x ]].

Observe that C = -(2-ϑ _1)(2ϑ _2-1) is constant. The matrix C is nonsingular if ϑ_1≠ 2. For this case, a left multiplication by its adjoint matrix C^* gives 

	
	
	(θ_1-2) (2 θ_2-1) u_1' 

	=  -  2/γ_1(ϑ_1-2) (1-  ϑ_2)x^2 u_0 ' + 2/γ_1 (ϑ_1-2) (2 ϑ_2-3) x u_0 , 
		   
	   

	(θ_1-2) (2 θ_2-1)  u_1 
	= ϕ(x) u_0' - 2/γ_1 (ϑ_1-1) (2 ϑ_2-3) x^2 u_0  ,     

where ϕ(x)=Φ which is given in (<ref>). 
We differentiate (<ref>) once and compare with (<ref>) to obtain (<ref>) and (<ref>). 

When ϑ_1=2, the functional equations (<ref>) (or equivalently (<ref>)) become 

    2 u_0^' = x  u_1^' , 
    2/γ_1 (1-ϑ_2) x^2  u_0^' + 2/γ_1 (3-2ϑ_2) x  u_0
    			- (1-2ϑ_2)u_1^'=0 ,
 
and this gives (<ref>) and (<ref>).

The implication (c) ⇒ (b) can be trivially obtained by performing reciprocal operations to the ones just described. 

The proof of (b) ⇒ (a) consists of showing that { Q_n }_n⩾ 0 is 2-orthogonal for (v_0,v_1) given by (<ref>) based on  (<ref>)-(<ref>) together with the 2-orthogonality of { P_n }_n⩾ 0 with respect to (u_0,u_1) and on account of conditions (<ref>). We will omit the details which can be followed in <cit.>, with the obvious adaptations to the threefold symmetric case. 


The vector functional 𝐔=(u_0,u_1) satisfies a matrix equation of Pearson type (<ref>). This somehow mimics the properties of the  classical orthogonal polynomials. However, we refrain ourselves from calling such polynomials 2-orthogonal classical, since other characteristic properties of the classical orthogonal polynomials can be interpreted in the context of 2-orthogonality (or multiple orthogonality) and give rise to completely different sets of polynomials.




As a straightforward consequence of statement (c) in the latter result, we derive equations for the moments of the Hahn-classical vector linear functional (u_0,u_1): 
 Under the same assumptions of Theorem <ref>,   the  moments of u_0 and u_1 satisfy 

    2/γ_1[ (3n+2)( (3n+1) (ϑ _1-1) (ϑ _2-1) 
    	 -  (ϑ_2 + ϑ _1-2)   )
    	 -  (ϑ _1-2)   ](u_0)_3n+3
    
    	= (3n+2)(3n+1)ϑ _1 (2 ϑ _2-1) (u_0)_3n , 
    
    
    	(u_0)_3n+1 = (u_0)_3n+2 =0 
    		  and    
    	(u_0)_0 =1 ,

while 

    ( u_1)_3n+1 = 2 (ϑ _1-1) ( (3n+2) ϑ _2- (3n+1) ) /γ _1(ϑ _1-2)  (2 ϑ _2-1) 
    		(u_0)_3n+3
    
    		-  (3n+1)   ϑ _1 (2 ϑ _2-1) (u_0)_3n, 
       if  ϑ_1 ≠ 2
    
    	( u_1)_3n+1 = 2(3n+1)/3n+2 (u_0)_3n   if  ϑ_1 = 2
    	
    
    	(u_1)_3n = (u_1)_3n+2 =0 
    		  and    
    	(u_1)_1 =1 .


 
From (<ref>)-(<ref>) it follows 

    ⟨(  ϕ(x) u_0)” 
    		+ ( 2 /γ _1(ϑ_2 + ϑ _1-2)x^2u_0)' 
    	+ 2/γ _1(ϑ _1-2) x u_0 ,  x^n ⟩
    	=0  ,      n⩾ 0,
 
 and 

    ⟨(ϑ _1-2) (2 ϑ _2-1) u_1 ,  x^n ⟩
    	=	⟨ϕ(x)u_0^ ' 
    			- 2 /γ _1(ϑ _1-1) (2 ϑ _2-3)  x^2 u_0,  x^n ⟩ , 
    			      if ϑ_1≠2,
    	
    
    	  x  u_1^ ' = 2  ⟨  u_0^ ' ,  x^n ⟩,  
    			   if ϑ_1=2,
 
which, on account of (<ref>), leads to (<ref>). 


Similarly, by taking into account the operations defined in (<ref>), we deduce that  relations (<ref>) imply
 (<ref>).


The "2-Hahn-classical" polynomials satisfy a third order recurrence relation (<ref>) and the  γ-recurrence coefficients have a specific rational structure. 
Here, we show that such expression for the γ-coefficients actually characterizes the  threefold symmetric "2-Hahn-classical" polynomials, by proving the reciprocal condition found in the work by Douak and Maroni  <cit.>. For a matter of completion, we obtain the γ-coefficients directly from the functional equations (<ref>), rather than from algebraic manipulations on the recurrence relations. 

 Let { P_n }_n⩾ 0 be a monic 2-. Then  { P_n }_n⩾ 0 satisfies  (<ref>) with 
 
    γ_n+2
    	= n+3/n+1  n(ϑ_n -1) +1  / (n+4)(ϑ_n+1 -1) +1 γ_n+1 ,

where 

    ϑ_n = (1-(-1)^n/2) 1-n+1/2 (1-ϑ_1)/1-n-1/2 (1-ϑ_1)
    			+( 1+(-1)^n/2)1-n/2 (1-ϑ_2)/1-(n/2-1) (1-ϑ_2) ,
    			     n⩾ 1,

with ϑ_1,ϑ_2 subject to (<ref>), 
if and only if  the sequence { Q_n(x)1/n+1 P_n+1'(x) }_n⩾ 0 is 2-orthogonal satisfying the recurrence relation (<ref>) with 

    γ_n 
    	= n/n+2ϑ_n γ_n+1,   for   n⩾ 1.



 We start by proving that if { Q_n }_n⩾ 0 satisfies the recurrence relation  (<ref>) with the γ-coefficients given by (<ref>), then the γ-coefficients in the recurrence relation of  { P_n }_n⩾ 0 are given by (<ref>)-(<ref>). The assumption means that { Q_n }_n⩾ 0 is 2-orthogonal and therefore { P_n }_n⩾ 0 is Hahn-classical. According to Theorem <ref>, { P_n }_n⩾ 0 is 2-orthogonal with respect to a vector functional (u_0,u_1) satisfying the differential equation (<ref>) and  { Q_n }_n⩾ 0 is 2-orthogonal for  𝐕=(v_0,v_1) given by (<ref>)-(<ref>). The relation Q_n(x) = 1/n+1P_n+1'(x) combined with the properties (<ref>) implies 

    ⟨  v_0, x^n Q_2n⟩ = - 1/2n+1⟨(x^n v_0)' ,  P_2n+1⟩

and 

    ⟨  v_1, x^n Q_2n+1⟩ = - 1/2n+2⟨(x^n v_1)' ,  P_2n+2⟩

for any n⩾ 0. In the latter identities we replace (v_0',v_1') and (v_0,v_1) by the respective expressions given in (<ref>) and (<ref>), 
to obtain 

    ⟨  v_0, x^n Q_2n⟩ 
    			= (2n+1)^-1 (1-nϕ_0,1'(0))  ⟨ u_1,  x^nP_2n+1⟩ 
    			,    for   n⩾ 0,
    	
       ⟨  v_1, x^n Q_2n+1⟩ 
    			= (2n+2)^-1(ψ'(0) - n ϕ_1,0”(0)/2) ⟨ u_0,   x^n+1 P_2n+2⟩
    			,    for   n⩾ 0,

which are the same as 

    ⟨  v_0, x^n Q_2n⟩ 
    			= (2n+1)^-1 (1-(1-ϑ_1)n)  ⟨ u_1,  x^nP_2n+1⟩ 
    			,    for   n⩾ 0,
       ⟨  v_1, x^n Q_2n+1⟩ 
    			= (2n+2)^-12/γ_1(1 - (1-ϑ_2)n ) ⟨ u_0,   x^n+1 P_2n+2⟩
    			,    for   n⩾ 0.

Taking into account (<ref>),  from the latter we obtain (<ref>) if we set ϑ_n as in  (<ref>). Observe that 

    ϑ_2n+1 =(n+1)  ϑ_1-n/nϑ_1 - (n-1)  and  ϑ_2n+2 = (n+1)ϑ_2 - n/n ϑ_2 - (n-1),      n⩾ 0,
 
and from this, we can clearly see that ϑ_n is actually a solution to the Riccati equation

    ϑ_n+3 + 1/ϑ_n+1 = 2,   for   n⩾ 0.


Conversely, 
if { P_n }_n⩾ 0  satisfies (<ref>) with the γ-coefficients given by (<ref>) where ϑ_n is given by (<ref>), then we prove that { Q_n }_n⩾ 0  satisfies a third order recurrence relation (<ref>) with γ-coefficients as in 
(<ref>). 
We differentiate the recurrence relation (<ref>)  once and then take into account the definition of Q_n(x)1/n+1 P'_n+1(x) obtain the structural relation 

    P_n (x) = (n+1) Q_n(x) - n x Q_n-1(x) + (n-2) γ_n-1 Q_n-3.
 
In (<ref>), we replace P_n+1,   P_n and P_n-2 by the expressions provided in the latter identity to obtain  

    (n+2) Q_n+1(x)
    		= 2(n+1) x Q_n(x) - n x^2 Q_n-1(x) + 2(n-2) γ_n-1 x Q_n-3(x)
    
    		-(n+1)(γ_n-1+γ_n) Q_n-2(x) 
    		- (n-4)γ_n-3γ_n-1Q_n-5(x),

which is valid for any n⩾ 0, under the assumption that Q_-n(x)=0. As { Q_n }_n⩾ 0 is a basis for 𝒫, there are coefficients 
ξ_n+1,ν such that 

    x Q_n(x) = ∑_ν=0^n+1ξ_n+1,ν Q_ν(x),

where ξ_n+1,n+1+ν=0 for any ν⩾ 0 and ξ_n+1,n+1=1 because  Q_n is monic. Based on this, we can also write

    x^2 Q_n(x) = ∑_ν=1^n+2 ∑_σ=ν-1^n+1ξ_n+1,σξ_σ+1,ν Q_ν(x) 
    				+  ∑_σ=0^n+1ξ_n+1,σξ_σ+1,0 Q_0(x).

The threefold symmetry of { P_n }_n⩾ 0 readily implies that of { Q_n }_n⩾ 0 and this means that 
ξ_n,ν=0 whenever n+ν≠ 0  3, so that ξ_n,n-1,  ξ_n,n-2,   ξ_n,n-4,   ξ_n,n-5, … are all equal zero. 
We replace the terms xQ_n,  xQ_n-3 and x^2Q_n-2 by the respective expressions given by (<ref>) and (<ref>) in the relation (<ref>) to obtain 

    (n+2) Q_n+1(x)
    		= 2(n+1) ∑_ν=0^n+1ξ_n+1,ν Q_ν(x)  + 2(n-2) γ_n-1∑_ν=0^n-2ξ_n-2,ν Q_ν(x)
    
    		- n ( ∑_ν=1^n+1 ∑_σ=ν-1^nξ_n,σξ_σ+1,ν Q_ν(x) 
    				+  ∑_σ=0^nξ_n,σξ_σ+1,0 Q_0(x))
     
    		-(n-1)(γ_n-1+γ_n) Q_n-2(x) 
    		- (n-4)γ_n-3γ_n-1Q_n-5(x).


For n⩾ 2, we  equate  the coefficients of Q_n-2, giving

    0= 2(n+1) ξ_n+1,n-2 +  2(n-2) γ_n-1    -n ξ_n,n-3 -(n-1)(γ_n-1+γ_n) .

In particular, for n=2 the latter becomes 

    0 = 4ξ_3,0-γ_1-γ_2  ,
 
and, because of (<ref>), we have 
γ_1 = 1/3( 4ϑ_1 -3)γ_2
so that 

    ξ_3,0 = ϑ_1/3γ_2.
 
If we assume that for some n⩾ 3, 

    ξ_n,n-3 = n-2/nϑ_n-2γ_n-1,
 
then (<ref>) implies 

    (n+2) ξ_n+1,n-2 = ((n-2) ϑ_n-2 -(n-3) )γ_n-1    +(n-1)γ_n ,
 
which, after writing  γ_n-1 in terms of γ_n via (<ref>), reads as 

    ξ_n+1,n-2 = n-1/n+1ϑ_n-1γ_n  ,
 
and therefore we conclude that 

    ξ_n+3,n = n+1/n+3ϑ_n+1γ_n+2  ,     for all  n⩾ 0.

Equating the coefficients of Q_n-5 in (<ref>) leads to 

    0=  (n+2) ξ_n+1,n-5    + ( 2(n-2) γ_n-1- n ξ_n,n-3)  ξ_n-2,n-5
    
    		- n   ξ_n,n-6 
    		- (n-4)γ_n-3γ_n-1,

for n ≥ 5.
The particular choice of n=5 in the latter identity becomes 

    0= 7 ξ_6,0    + (2  
    		-  ϑ_3) ϑ_1γ_2γ_4  - γ_2γ_4,
 
after replacing ξ_5,2 and ξ_3,0 by the expressions provided by (<ref>). Using the identity (<ref>) for n=4, we then conclude that
ξ_6,0  =0. Now, suppose that ξ_n,n-6 =0 for some n⩾ 6.  Identity (<ref>) tells that 

    (n+2) ξ_n+1,n-5    = -  ( 2(n-2) γ_n-1- n ξ_n,n-3)  ξ_n-2,n-5 
    		+ (n-4)γ_n-3γ_n-1 ,
 
which, because of (<ref>), becomes 

    (n+2) ξ_n+1,n-5    = - (n-4) ( 2  -  ϑ_n-2) ϑ_n-4γ_n-1γ_n-3
    		+ (n-4)γ_n-3γ_n-1 ,
 
and hence, due to (<ref>), we conclude 

    ξ_n+6,n = 0 ,      for all   n⩾ 0.
 
Now, equating the coefficients of Q_n+1-3j for j⩾ 3 in (<ref>) gives

    0 
    		= (n+2)  ξ_n+1,n+1-3j    + (  2(n-2) γ_n-1   - n ξ_n,n-3) ξ_n-2,n+1-3j
    
    		- n ∑_σ=n-3j^n-6ξ_n,σξ_σ+1,n+1-3j ,
 
which is also 

    0 
    		= (n+2)  ξ_n+1,n+1-3j    + (  2(n-2) γ_n-1   - n ξ_n,n-3) ξ_n-2,n+1-3j
    
    		- n ∑_σ=2^jξ_n,n-3σξ_n+1-3σ,n+1-3j .
 
Suppose that ξ_n,n-3k=0 for each k=2,3,…, j-1 and n⩾ 3k, then the latter becomes 

    (n+2)  ξ_n+1,n+1-3j     
    	= n ξ_n,n-3j ,
 
which yields 


    ξ_n+1,n+1-3j
    	= 3j (3j+1)/(n+2)(n+1)ξ_3j,0.

In particular this implies that 

    ξ_3j+1,1   = 3j /(3j+2)ξ_3j,0   and   ξ_3j+2,2   = 3j +1 /(3j+3)ξ_3j,0 .


Finally we compare the coefficients of Q_0 in (<ref>) to obtain 

    0 
    		= (n+2)   ξ_n+1,0    + 2(n-2) γ_n-1ξ_n-2,0 
    		- n ( ∑_σ=0^n-1ξ_n,σξ_σ+1,0),
 
for n ≥ 6.
Recall that ξ_k,0 =0 for k≠ 0 3, so that the latter identity simplifies to 

    0 
    		= (3n+3)   ξ_3n+3,0    + 6n γ_3nξ_3n,0 
    		- (3n+2) ( ∑_σ=1^n-2ξ_3n+2,3σ+2ξ_3σ+3,0) -n ξ_3n+2,2ξ_3,0 ,

for n ≥ 6.
We have already seen that ξ_6,0=0. Proceeding by induction, we promptly observe that if ξ_3j,0=0 for j=0,1,…, n, then 
identities (<ref>)–(<ref>) allow us to conclude from (<ref>) that ξ_3(n+1),0=0 and this implies 

    ξ_n,k=0,     for  0⩽ k ⩽ n-4.
 
As a result, we conclude that 

    xQ_n = Q_n+1 + ξ_n+1,n-2 Q_n-2,   with   ξ_n+1,n-2 =  n-1/n+1ϑ_n-1γ_n≠ 0,     for   n⩾ 2,
 
with Q_j(x) = x^j for j=0,1,2. This means that {Q_n}_n⩾ 0 is 2-orthogonal and threefold symmetric. 




In  <cit.> Douak and Maroni have highlighted several properties of threefold-symmetric (therein referred to as "2-symmetric") 2-classical polynomials including a differential equation of third order. Here we show that such differential equation actually characterizes these polynomials, bringing to the theory a Bochner-type characterization. 
 
 Let {P_n}_n⩾ 0 be a  threefold symmetric 2-orthogonal polynomial sequence satisfying (<ref>). 
The sequence {P_n}_n⩾ 0 is Hahn-classical  if and only if each P_n satisfies  

    (a_n x^3- b_n) P_n”' + c_n x^2 P_n” + d_nx P_n' = e_n P_n ,

where 


    a_n    =    
    		(ϑ_n- 1) (ϑ_n+1 - 1) 
     
    	b_n    =   γ _n ((n-1) ϑ _n-1-n+2) (n ϑ _n-n+1) ((n+1) ϑ _n+1-n)/n (n+1) 
     
    	c_n    =   ϑ_nϑ_n+1-1 -(n-3) ( ϑ_n -1)( ϑ_n+1-1)
     
    	d_n    =    n  ϑ_n+1 - (n-1)  ϑ_n (2  ϑ_n+1-1)
    
    	e_n    =    n  ϑ_n+1,


for all n ≥ 1, 
and    ϑ_n given in (<ref>) with initial values ϑ_1=3(γ_1+γ_2)/4γ_2 and ϑ_2 = 3 (γ _1+γ _2)/10 γ _3+4/5 subject to ϑ_1,ϑ_2∉{n-1/n :  n⩾ 1} .

 
The necessary condition was proved in <cit.>, where Douak and Maroni have shown that  threefold-symmetric (therein referred to as "2-symmetric") 
2-orthogonal polynomials satisfying Hahn's property are solutions to (<ref>) under the definitions  (<ref>) and (<ref>)-(<ref>)  with 

    b_n =  γ _n+3((n+3) ϑ _n-n+2) ((n+4) ϑ _n+1-(n+3)) ((n+5) ϑ _n+2-(n+3))/(n+3)(n+4)  ,
 
where  ϑ_n is given in (<ref>). It turns out that b_n can be written as in (<ref>) because, according to Theorem <ref>, the γ-recurrence coefficients are recursively given by (<ref>) which implies 

    γ_n+3
    	= (n+4)(n+3)/n(n+1)
    ×( (n+1)(ϑ_n+1-1) +1  ) ( n(ϑ_n -1) +1  )( (n-1)(ϑ_n-1-1) +1  )/( (n+5)(ϑ_n+2 -1) +1 )( (n+4)(ϑ_n+1 -1) +1 )( (n+3)(ϑ_n -1) +1 ) γ_n.
 

Conversely, suppose that the 2-orthogonal polynomial sequence {P_n}_n⩾ 0 fulfils  (<ref>) under the definitions (<ref>)-(<ref>) with ϑ_n given in (<ref>). Observe that ϑ_n satisfies the Riccati equation (<ref>) with initial conditions ϑ_1=3(γ_1+γ_2)/4γ_2 and ϑ_2 = 3 (γ _1+γ _2)/10 γ _3+4/5, both assumed to be different from k-1/k for all integers k⩾ 1. 
For values of n⩾ 3, we consider the expansion P_n(x) = x^n + λ_n x^n-3 + …  and insert it in the differential equation (<ref>). We equate the coefficients of x^n and this gives 

    λ_n  =- (n-2) (n-1) n b_n/3 ((n-3) ϑ _n-n+4) ((n-2) ϑ _n+1-n+3)
    	 ,      n⩾ 0.
 
On the other hand, from the recurrence relation (<ref>) we deduce 

    γ_n = λ_n+1 - λ_n+2,    n⩾ 0.
 
After combining the latter two expressions we obtain  

    γ_n 
    	=    - (n-1) n (n+1) b_n+1/3 ((n-2) ϑ _n+1-n+3) ((n-1) ϑ _n+2-n+2)
       +  n (n+1) (n+2) b_n+2/3 ((n-1) ϑ _n+2-n+2) (n ϑ _n+3-n+1) ,

where b_n is given in (<ref>). Based on (<ref>), we consider the following substitutions in the latter expression 

    ((n-2) ϑ _n+1-n+3) = ((n-1) ϑ _n-1-(n-2)) /ϑ_n-1

and

    (n ϑ _n+3-(n-1)) = ((n+1) ϑ _n+1-n) /ϑ_n+1 ,

to find

    γ_n 	=n ϑ _n((n+2) ϑ _n+2-(n+1))
    	(	- γ_n+1(n-1)    ϑ _n-1((n+1) ϑ _n+1-n) /3  (n+2) ((n-1) ϑ _n-1-(n-2)) 
    
    	 +  γ_n+2(n+1)  ϑ _n+1((n+3) ϑ _n+3-(n+2))/ 3(n+3) (n ϑ _n-(n-1) ) ) ,
 
which is the same as 

    γ_n/n((n+3) ϑ _n-(n+2))	
     = 
      γ_n+2 (n+1) ((n+4) ϑ _n+1-(n+3))/ 3(n+3) (n ϑ _n-(n-1) ) 
    
     - γ_n+1(n-1)     ϑ _n-1((n+1) ϑ _n+1-n)  /3  (n+2) ((n-1) ϑ _n-1-(n-2))

after taking the following substitutions (derived from  (<ref>))

    ϑ _n((n+2) ϑ _n+2-(n+1)) = ((n+3) ϑ _n-(n+2)) ,
 
and

    ϑ _n+1((n+3) ϑ _n+3-(n+2)) =((n+4) ϑ _n+1-(n+3)) .
 
We subtract γ_n+1/  (n+2) ((n-1) ϑ _n-1-(n-2)) from both sides of (<ref>) and this leads to

    -3 γ_n/n((n+3) ϑ _n-(n+2))γ_n+1( 1 - 	  γ_n+1/γ_nn((n+3) ϑ _n-(n+2))/  (n+2) ((n-1) ϑ _n-1-(n-2)) )
    		
        = ( 
    			 1 -  γ_n+2/γ_n+1 (n+1) ((n+4) ϑ _n+1-(n+3))/ (n+3) (n ϑ _n-(n-1) ) ) ,
 
which implies 

    (∏_ℓ=1^n (-3) γ_ℓ/ℓ((ℓ +3) ϑ _ℓ-(ℓ +2))γ_ℓ +1) 
    	( 1 - 	  γ_2/γ_1(4 ϑ _1-3)/  3 )
    
    	=  ( 
    			 1 -  γ_n+2/γ_n+1 (n+1) ((n+4) ϑ _n+1-(n+3))/ (n+3) (n ϑ _n-(n-1) ) ).
 
The assumption on the initial value for ϑ_1 readily implies the left-hand side of the latter equality to be zero and therefore we conclude 

    0 =  ( 
    			 1 -  γ_n+2/γ_n+1 (n+1) ((n+4) ϑ _n+1-(n+3))/ (n+3) (n ϑ _n-(n-1) ) )
    	 for all   n⩾ 1.
 
Now, Theorem <ref> ensures the 2-orthogonality of the sequence  { Q_n(x)1/n+1 P_n+1'(x) }_n⩾ 0, which means that {P_n}_n⩾ 0 is Hahn-classical. 




After a single differentiation of (<ref>), we obtain a differential equation for the polynomials Q_n and these satisfy: 

    (a_n x^3- b_n) Q_n”' + (c_n+3a_n) x^2 Q_n” + (d_n+2c_n) x Q_n' = (e_n-d_n) Q_n,      n⩾ 0,
 
with a_n,b_n,c_n,d_n and e_n given by  (<ref>)–(<ref>). 





Theorem <ref> shows that a threefold symmetric 2-orthogonal polynomial sequence is Hahn-classical if and only if the γ-recurrence coefficients in 
(<ref>) can be written as (<ref>), provided that ϑ_1,ϑ_2≠n/n+1 for all positive integers n. 


If ϑ_1,ϑ_2⩾ 1, then γ_n and γ_n  are both positive for all integers n⩾ 1. 
Furthermore, from (<ref>)-(<ref>) we readily see that γ_n and γ_n are two rational functions in n, both having the same asymptotic behavior: 

    γ_n = c n^α +o(n^α)  as   n→ +∞.

Hence, as a result of Theorem <ref> together with Theorem <ref> and Theorem <ref> the two linear functionals associated with a Hahn-classical threefold symmetric 2- admit the following integral representation. 
 
 Let {P_n}_n⩾0 be a threefold symmetric and 2-  with respect to the vector linear functional (u_0,u_1) fulfilling (<ref>)-(<ref>). If ϑ_1,ϑ_2 ⩾ 1, then {P_n}_n⩾0  satisfies the recurrence relation (<ref>) with γ_n>0 and 
u_0 and u_1 admit the integral representation

    ⟨ u_k, f ⟩ = 
    1/3( ∫_0^b f(x) 𝒰_k(x)   x
    		+ ω^2k-1∫_0^bω f(x) 𝒰_k(  ω^2 x)  x
    		+  ω^1-2k∫_0^bω^2 f(x) 𝒰_k(  ω x)  x)
    		,
 
(with k=0,1)  for any polynomial f, with  ω = ^ 2π i/3   and  b=lim_n→∞( 27/4γ_n),  provided that there exist two twice differentiable functions 𝒰_0 and 𝒰_1 mapping [0,b)  to ℝ   such that 𝒰_0 is solution to 

    ( ϕ(x) 𝒰_0 (x))” 
    		+ ( 2(ϑ_2 + ϑ _1-2) /γ _1x^2 𝒰_0(x))' 
    		 	+ 2 (ϑ _1-2) /γ _1x 𝒰_0(x) = λ_0 g_0(x) ,

and  𝒰_1 is given by


	
	
	(ϑ_1-2) (2 ϑ_2-1) 𝒰_1(x)

	=	
		ϕ(x)𝒰_0'(x)
			- 2 (ϑ_1-1) (2 ϑ_2-3)x^2 /γ_1 𝒰_0(x) 
			+λ_1 g_1(x)  ,
		
			    if ϑ_1≠ 2, 
			
 
			
	x 𝒰_1^'(x) = 2   𝒰_0^'(x) +λ_1 g_1(x)  ,
	    if ϑ_1=2, 


where ϕ(x) =  ϑ _1 (2 ϑ _2-1)-2 (ϑ _1-1) 
			(ϑ _2-1) /γ _1x^3 ,  
and satisfying

    lim_x→ b f(x) ^l / x^l𝒰_0(x) = 0  
    		,   for any  l∈{0,1}  and   f∈𝒫,
          ∫_0^b𝒰_0(x)  x=1,

where λ_k is a complex constant (possibly zero) and g_k a function whose moments vanish identically on the support of 𝒰_k.


 
According to Theorem <ref>, the threefold symmetry of a 2- ensures the existence of two orthogonality measures μ_0 and μ_1 (respectively defined by the vector functional (u_0,u_1)) supported on a starlike set S on the three rays of the complex plane. On the other hand, Theorem <ref> tells that u_0 satisfies (<ref>) and u_1 fulfils (<ref>)-(<ref>). Based on these equations, we seek an integral representation for both u_0 and u_1 via a weight function 𝒰 such that (<ref>) holds for any polynomial f. Such a representation readily ensures the threefold symmetry of (u_0,u_1). Identity (<ref>)  is the same as  

    ⟨( ϕ(x)u_0)” 
    		+ ( 2(ϑ_2 + ϑ _1-2) /γ _1x^2u_0)' 
    		 	+ 2 (ϑ _1-2) /γ _1x u_0  , f ⟩=0, 
    			    ∀ f∈𝒫,
 
which, because of (<ref>), reads as 

    ⟨ u_0, ϕ(x) f”(x) 
    			-   2(ϑ_2 + ϑ _1-2) /γ _1x^2 f'(x) 
    			+ 2 (ϑ _1-2) /γ _1x f(x)
    			 ⟩=0,     ∀ f∈𝒫,

with 

    ϕ(x) = ϑ _1 (2 ϑ _2-1)-2 (ϑ _1-1) 
    			(ϑ _2-1) /γ _1x^3.

We seek a function 𝒲_0, at least twice differentiable, defined on an open set D containing the piecewise differentiable curve 𝒞 containing all the zeros of {P_n}_n⩾0 that is a subset of S=Γ_0∪Γ_1 ∪Γ_2 (represented in Fig.<ref>)  and   
such that 

    ⟨ u_0,f ⟩ = ∫_𝒞 f(x) 𝒲_0(x)  x

holds for every polynomial f. 
In the light of Theorem <ref>,  
the support of the weight function is 𝒞=Γ_0∪Γ_1∪Γ_2
where Γ_1 and Γ_2 are the two straight lines starting at some point bω and bω^2 (respectively) and ending at the origin, while Γ_0 corresponds to the straight line on the positive real axis starting at the origin and ending at b. Here b is a positive real number or can represent a point at infinity. 
Thus, from (<ref>), the weight function 𝒲_0 we seek must be such that 

    ∫_𝒞𝒲_0(x) (ϕ(x) f”(x) 
    			-   2(ϑ_2 + ϑ _1-2) /γ _1x^2 f'(x) 
    			+ 2 (ϑ _1-2) /γ _1x f(x) )
    			  x
    			= 0, ∀ f∈𝒫.

We use (complex) integration by parts to deduce 

    ∫_𝒞(  (  ϕ(x) 𝒲_0(x) )”  
    			+(  2(ϑ_2 + ϑ _1-2) /γ _1x^2  𝒲_0(x) )' 
    			+ 2 (ϑ _1-2) /γ _1x  𝒲_0(x) ) f(x)
    			  x
    	
        + ∑_j∈{1,2}.  ( f'(x) ϕ(x) 𝒲_0(x) +   f(x)( 2(ϑ_2 + ϑ _1-2) /γ _1x^2  𝒲_0(x) 
    			- (ϕ(x) 𝒲_0(x))'))|_ω^j b^b
       		= 0 ,

 which must hold for every polynomial f. This condition is fulfilled if each of the following conditions hold: 
  
 
  (a) 𝒲_0 is a solution of the second order differential equation 

    (  ϕ(x) 𝒲_0(x) )”  
    			+(  2(ϑ_2 + ϑ _1-2) /γ _1x^2  𝒲_0(x) )'  
    			+ 2 (ϑ _1-2) /γ _1x 𝒲_0(x)
    	= λ_0 g_0(x),

for some constant λ_0 (possibly zero) and a function g_0 representing the null linear functional on the vector space of polynomials supported on 𝒞, that is 

    ∫_𝒞 g_0(x) f(x)  x =0 ,     for every polynomial  f.



  (b) The solution 𝒲_0 of (<ref>) satisfies the boundary conditions

    .  ( f'(x) ϕ(x) 𝒲_0(x) +   f(x)( 2(ϑ_2 + ϑ _1-2) /γ _1x^2  𝒲_0(x) 
    			- (ϕ(x) 𝒲_0(x))'))|_ω^j b^b=0,

with j=1,2,  for every polynomial f. 


  (c) All the moments of 𝒲_0 coincide with those of u_0, that is, 

    (u_0)_n 
    	 = ∫_S𝒲_0(x) x^n   x   ,      n⩾ 0.

This corresponds essentially to (<ref>), insofar as {x^n}_n⩾ 0 forms a basis for 𝒫. In particular, and in the light of Proposition <ref>, the threefold symmetry of the linear functional u_0 means that (u_0)_n = 0 if n≠ 0 3.
  


Observe that if y is a solution of (<ref>), then so are the functions ω^j  y(ω^j x) for j=1,2.  Taking into consideration the threefold symmetry of u_0, it follows that 

    𝒲_0(x) = {[         1/3𝒰_0(x)                if      x∈Γ_0=[0,b),; -1/3ω^2𝒰_0(ω^2 x)                if            x∈Γ_1,;     -1/3ω𝒰_0(ω x)                if            x∈Γ_2,;                 0                if  x∉𝒞=Γ_0∪Γ_1∪Γ_2. ].

where 𝒰_0: [0,b) →ℝ is an at least twice differentiable function that is a  solution of  the differential equation (<ref>) satisfying the conditions (<ref>), so that (<ref>) holds. Here,  Γ_1 and Γ_2 are the two straight lines starting at some point bω and bω^2 (respectively), with ω = ^2π i /3,  and ending at the origin, while Γ_0 corresponds to the straight line on the positive real axis starting at the origin and ending at b. Hence  (<ref>) is proved for k=0. 

Similarly, an integral representation for the linear functional u_1 can be obtained from (<ref>) by seeking a differentiable function  𝒲_1 defined on 𝒞 such that 

    ⟨ u_1, f(x) ⟩ = ∫_𝒞 f(x) 𝒲_1(x)  x,  for every polynomial  f.

If ϑ_1≠ 2, then u_1 satisfies (<ref>) and this gives

    ⟨ u_1, f ⟩ 
    	=  (ϑ _1-2)^-1(2 ϑ _2-1)^-1⟨ϕ(x) u_0' - 2 (ϑ _1-1) (2 ϑ _2-3) /γ _1x^2 u_0, f ⟩.

Based on the properties (<ref>), the latter becomes

    ⟨ u_1, f ⟩ 
    	=  (ϑ _1-2)^-1(2 ϑ _2-1)^-1⟨ u_0 , - ( ϕ f)' - 2 (ϑ _1-1) (2 ϑ _2-3) /γ _1x^2 f ⟩.

Taking into consideration (<ref>), we then have 

    ⟨ u_1, f ⟩ 
    	=
    	  - (ϑ _1-2)^-1(2 ϑ _2-1)^-1
    ×∫_𝒞𝒲_0(x)  
    	(( ϕ(x) f(x))' - 2 (ϑ _1-1) (2 ϑ _2-3 ) /γ _1x^2 f(x) )  x.

We perform integration by parts on the first term of the integral, to obtain 

    (ϑ _1-2) (2 ϑ _2-1)  ⟨  u_1, f ⟩
    [                                         =∑_j∈1,2( .  ϕ(x) 𝒲_0(x)f(x) |_bω^j^b); + ∫_𝒞   f(x) 
    	(  ϕ(x)   𝒲_0'(x) - 2 (ϑ _1-1) (2 ϑ _2-3 ) /γ _1x^2𝒲_0(x) )  x. ]

The vanishing conditions (<ref>) at the end points of the contour readily imply that the integrated term vanishes identically for every polynomial f considered. As a consequence, a weight function associated with u_1 corresponds to 

    𝒲_1 (x ) 
    	= 
    	(ϑ _1-2)^-1(2 ϑ _2-1)^-1(  ϕ(x)   𝒰'(x) - 2 (ϑ _1-1) (2 ϑ _2-3 ) /γ _1x^2𝒲_0(x) ) 
    	
    
    	+ λ_1 g_1(x),

where λ_1 represents a complex constant (possibly zero) and g_1 a function representing the null linear functional. 
By construction, this choice for the weight function guarantees that (<ref>) is well defined and, in particular, all the moments satisfy the threefold symmetry property described in Proposition <ref>.
So, we conclude that (<ref>) also holds for k=1 with 𝒰_1 given by (<ref>) when  ϑ_1≠2. 

If ϑ_1=2, then using an entirely analogous approach we deduce (<ref>)  from  (<ref>). 


 An alternative to the representation (<ref>) stated in Theorem <ref> is 

    ⟨ u_k,f ⟩ = ∫_𝒞 f(x) 𝒲_k(x)  x,      (k=0,1),

where 𝒲_k(x) is given in (<ref>). 

The weight functions 𝒲_0 and 𝒲_1 are threefold symmetric in the sense that they satisfy the following rotational invariant property: 

    ω^j 𝒲_k(ω^j x) = 𝒲_k(x),      j=0,± 1,±2, … ,  for each  k=0,1.
 








We consider the following equivalence relation between two polynomial sequences {P_n}_n⩾ 0 and {B_n}_n⩾ 0: 

    {P_n}_n⩾ 0 ∼ {B_n}_n⩾ 0 iff∃  a∈ℂ\{0} ,  b∈ℂ such that   
    	B_n(x) = a^-n P_n(ax+b),

for all  n⩾ 0.
In this case, if  {P_n}_n⩾ 0 satisfies (<ref>) then {B_n}_n⩾ 0 satisfies 

    B_n+1(x) = (x-β_n-b/a) B_n(x) - α_n/a^2B_n-1(x) - γ_n-1/a^3 B_n-2(x) ,      n⩾ 1,

with initial conditions B_-2(x) = B_-1(x) =0 and B_0(x)=1.


Observe that ϑ_n is a solution of a Riccati equation (<ref>) for which ϑ_n=1 is a trivial solution. Depending on the initial conditions 
ϑ_1 and ϑ_2, there are four sets of independent solutions which in fact give rise to four equivalence classes of the threefold symmetric 2-Hahn-classical polynomials. In other words, up to a linear transformation of the variable, there are at most four distinct families of threefold symmetric 2-orthogonal Hahn-classical polynomials, which we single out: 


   Case A: ϑ_1 =1= ϑ_2. This implies that ϑ_n =1 for all n⩾ 0. 



    Case B_1: ϑ_1 ≠1 but ϑ_2=1 so that by setting ϑ_1 = μ+2/μ+1 it follows that 

    ϑ_2n-1 = n+μ+1/n+μ  and  ϑ_2n =  1 ,      n⩾ 1.





   Case B_2: ϑ_1 =1 but ϑ_2≠1 so that by setting ϑ_2 = ρ +2/ρ+1 it follows that

    ϑ_2n-1 = 1   and  ϑ_2n =  n+ρ+1/n+ρ ,      n⩾ 1.




   Case C: ϑ_1 ≠1 and ϑ_2≠ 1 and hence by setting ϑ_1 = μ+2/μ+1 and ϑ_2 = ρ +2/ρ+1 it follows that 

    ϑ_2n-1 = n+μ+1/n+μ  and  ϑ_2n =  n+ρ+1/n+ρ ,      n⩾ 1.



All these cases were highlighted in <cit.>, where expressions for the recurrence coefficients for the three components of the cubic decomposition of each of these threefold symmetric Hahn-classical polynomials were deduced. Only some of these polynomials were studied in detail: those in case A in <cit.> and a few subcases of Case C in <cit.>. The representations for the orthogonality measures that we describe in the next sections are new. The support of these orthogonality measures contain all the zeros of the polynomial sequences. Even for particular cases that already appeared in the literature, the integral representations for the orthogonality measures were either not given or given on the positive real line (with oscillating terms), which is only part of the starlike set S. 
In the next subsections we fully describe all these cases in detail.  

As observed in <cit.>, the following limiting relations take place 

    case C B_1 case A ,

and also 

    case C B_2 case A .

It turns out that cases B_1 and B_2 are related to each other by differentiation, as explained in Section <ref> and Section <ref>.



 §.§ Case A
   

In the light of Theorem <ref>, for this choice of initial conditions one has ϑ_n=1 for all integers n⩾ 1, so that the γ-recurrence coefficients are  given by 
γ_n+1
	=(n+1)(n+2)γ_1/2
 and  
γ_n+1 = γ_n+1 for n⩾ 0. As a consequence, Q_n(x) = P_n(x) for all n⩾0 and, for this reason, the 2- {P_n}_n⩾ 0 is an Appell sequence.

Since the 2-orthogonality property is invariant under any linear transformation, we can set γ_1=2, and, with this choice, recalling (<ref>)-(<ref>) it follows  that


    γ_n+1 = γ_n+1=  (n+2)(n+1),      n⩾0,
 

so that 

    [ P_n+1(x) = xP_n(x) -  n(n-1)P_n-2(x),        n⩾ 2,;          P_0(x)=1,  P_1(x) =x   and  P_2(x) = x^2. ]




Figure <ref> illustrates the behavior of the largest zero in absolute value (plot generated in Mathematica): 
bounded from above by the curve y=3/2^2/3 n^2/3 suggested by 
Theorem <ref>.

Regarding the integral representation of the corresponding orthogonality measures, we have: 

 The threefold symmetric polynomial sequence {P_n}_n⩾ 0 defined by the recurrence relation (<ref>) is 2-orthogonal with respect to (u_0,u_1) admitting the integral representation  (<ref>), where 

    𝒰_0(x) = Ai(  x)
    	  and  𝒰_1(x) = - Ai'(  x),

and b=+∞. Moreover, the sequence {Q_n(x)1/n+1P_n+1'(x)}_n⩾ 0 coincides with {P_n}_n⩾ 0 (i.e., {P_n}_n⩾ 0 is an Appell sequence). 


Here and in what follows Ai and Bi are the Airy functions of the first and second kind (see <cit.>), respectively. 

 
Under the assumptions, by virtue of Theorem <ref> and Theorem <ref>, this polynomial sequence is 2-orthogonal with respect to (u_0,u_1)  satisfying  (<ref>)-(<ref>) which reads as 

    {[ u_0” - x u_0=0 ,;   u_1  =	- u_0', ].
 
from which we conclude that the corresponding sequence of the moments {(u_0)_n}_n⩾0 and  {(u_1)_n}_n⩾0 satisfy

    (u_0)_n+3 = 	  (n+1)(n+2)  (u_0)_n   and    
    	(u_1)_n = n (u_0)_n-1 ,

with initial conditions (u_0)_0=1 and (u_0)_1=(u_0)_2=(u_1)_0=0. This implies 

    [                                      (u_0)_3n                                             =                                (3n)!/3^n (n!)          and    
    	(u_0)_3n+1 =(u_0)_3n+2 =0 ,;                                    (u_1)_3n+1                                             =                              (3n+1)!/3^n (n!) and    
    	(u_1)_3n =(u_1)_3n+2 =0 ,      n⩾ 0. ]


According to Theorem <ref>,  u_0 and u_1 admit the representation  (<ref>)  provided that there exist two twice differentiable functions 𝒰_0 and 𝒰_1 from ℝ  to ℝ   that are solutions to 

    {[ 𝒰_0” (x) 
    		 	 - x  𝒰_0(x) = λ_0 g_0(x),;   -𝒰_1(x)
    	=	
    		 𝒰_0'(x) + λ_1 g_1(x)  , ].

and satisfying

    lim_x→ω^j∞ f(x) ^l / x^l𝒰_0(x) = 0  
    		,   for    j,l∈{0,1,2}  and   f∈𝒫,
          ∫_0^∞𝒰_0(x)  x=1 ,

where λ_k ∈ℂ (possibly zero) and g_k(x)≠0 are rapidly decreasing functions, locally integrable, representing the null functional (k=0,1). For λ_0=0, the general solution of the first equation in (<ref>) can be written as

    y(x) = c_1 Ai(  x) + c_2 Bi(  x) ,

for arbitrary constants c_1,c_2. Observe that (<ref>) is realized if we take c_2=0 (see <cit.>). Since 
 <cit.>

    ∫_ 0^+∞ x^n Ai(x)   x = Γ(n+1)/3^n/3+1Γ(n/3+1) < +∞,     for all  n⩾ 0,

the result now follows if we take c_1=1 and λ_1=0. 






  §.§.§ The differential equation
 

Under these assumptions (<ref>) becomes 

    -  P_n”' (x) + x P_n'(x) =n P_n(x),      n⩾ 0,
 
whose general solution is given by 

    P_n(x) 
    	   =    c_0  12-n/31/3,2/3x^3/9
    	+ c_1  x   12-n-1/32/3,4/3x^3/9
    	+ c_2   x^2  12-n-2/34/3,5/3x^3/9,

with integration constants c_0, c_1 and c_2. Observe that for each n there is one polynomial solution and we have 

    P_3n(x) = P_n^[0] (x^3)    with    P_n^[0] (x) = (-9)^n (1/3)_n (2/3)_n 12-n1/3,2/3x/9 , 
        P_3n+1(x)= x  P_n^[1] (x^3)    with    P_n^[1] (x)  = (-9)^n (2/3)_n (4/3)_n 12-n2/3,4/3x/9 , 
        P_3n+2(x)= x^2  P_n^[2] (x^3)    with    P_n^[2] (x)  = (-9)^n (4/3)_n (5/3)_n 12-n4/3,5/3x/9.


	



  §.§.§ The cubic decomposition


In <cit.> the cubic decomposition of an Appell 2-orthogonal sequence has been highlighted. This also happens to be threefold symmetric, and the corresponding polynomials  were called Hermite-type 2-orthogonal polynomials. The three 2-s  {P_n^[j]}_n⩾ 0 (j=0,1,2) in the cubic decomposition  of {P_n}_n⩾ 0 are 2-orthogonal with respect to weights involving modified Bessel functions of the second kind, studied in <cit.>. 

Following Lemma <ref>, each of the three polynomial sequences {P_n^[j]}_n⩾ 0, with  j∈{0,1,2}, is 2-orthogonal and satisfies the third order recurrence relation (<ref>), where  


    [                                β_n^[j] =  3 (j^2+6 j n+j+9 n^2)+9 n+2,      n⩾ 0,;                         α_n^[j] = 3 (j+3 n-2) (j+3 n-1)^2 (j+3 n)
    			,      n⩾ 1,; γ_n^[j] = (j+3 n-2) (j+3 n-1) (j+3 n) (j+3 n+1) (j+3 n+2) (j+3 n+3) ,       n⩾ 2.;                                                                                   ]

It should be noted that the Appell polynomials {P_n}_n⩾ 0 and the components in the cubic decomposition were treated in <cit.> under a different normalisation, namely by considering 2/γ_1 =9.  However the integral representation provided in  <cit.> is supported on the positive real axis, and therefore different from the one given here. In <cit.> the authors have also studied these polynomials, where the focus was put on algebraic properties, including generating functions. 

 This Appell  sequence {P_n}_n⩾ 0  has already appeared in the literature in other contexts, often not recognized as 2-s. For instance, it is linked to the Vorob'ev-Yablonski polynomials associated with rational solutions of the second Painlevé equations <cit.>. 
Furthermore, in <cit.> Widder studied the so-called Airy transform defined as follows

    u(x,t)= ∫_-∞^+∞ f(y) 1/(3t)^1/3Ai(y-x/3t^1/3)  y.

Up to a scaling, it maps the sequence of monomials to the Appell 2- {P_n}_n⩾ 0 <cit.>. For this reason, in <cit.> (and also <cit.>) they have been referred to as "Airy polynomials". 


Figure <ref> corroborates the statement in Proposition <ref>: the zeros on the positive real axis of three consecutive polynomials interlace and all the other non-zero zeros are rotations of 2π/3 of them. 












 §.§ Case B_1
   

In this case, ϑ_2n=1 and ϑ_2n-1 = n+μ+1/n+μ, under the constraint that μ≠ -n, for all n⩾1. Therefore, the threefold symmetric sequence {P_n(·;μ;γ_1)}_n⩾ 0 satisfies the recurrence relation (<ref>), where 

    γ_2n +1 = (n+1)(2n+1) (μ+2)/(3n+μ+2)γ_1, 
       γ_2n +2  =  (n+1)(2n+3)(n+μ+1) (μ+2)/(3n+μ+2)(3n+μ+5)γ_1,    n⩾ 0.

Observe that γ_1 is a mere scaling factor and, from this point forth, we set γ_1 = 2/3(μ+2). The discussion is therefore about on 
{P_n(·;μ) P_n(·;μ,2/3(μ+2))}_n⩾ 0. The general case  {P_n(·;μ;γ_1)}_n⩾ 0 can be deduced after a linear transformation of the variable: 

    P_n(x;μ;γ_1) =a^-n P_n(ax;μ),   with  a= (23(μ+2)γ_1)^1/3,       n⩾ 0.

The two sequences are equivalent, under (<ref>).

Hence, {P_n(·;μ)}_n⩾ 0 satisfies  (<ref>) with γ_n+1γ_n+1(μ) given by 

    γ_2n +1 =2/3 (n+1)(2n+1)  /(3n+μ+2),   γ_2n +2 =  2/3 (n+1)(2n+3)(n+μ+1) /(3n+μ+2)(3n+μ+5)  ,      n⩾ 0,

and it is 2-orthogonal for (u_0,u_1)  for which 

    1/3  u_0” + x^2 u_0'- (μ -2) x u_0  = 0
    	    and  
    		u_1
    			=	- (μ +2)/μ( u_0' + 3x^2 u_0)
    			   if  μ≠0 , 
    
    		x u_1' = 2 u_0'   if  μ=0.

Consequently, 

    (u_0)_3n 
    		= (1/3)_n (2/3)_n  /(μ +2/3)_n  and   (u_0)_3n+1 =(u_0)_3n+2 =0 ,

and 

    (u_1)_3n+1 
    		= (2/3)_n (4/3)_n/(μ+5/3)_n  and   (u_1)_3n =(u_1)_3n+2 =0.

Observe that in order to have γ_n+1>0 for all n⩾ 0, the constraint μ>-1 needs to be imposed. 

Particular choices of the parameter μ make the γ-coefficients linear functions in n, namely: 

    for  μ=-1/2:   γ_2n (-1/2) = 4n/27,   γ_2n +1(-1/2)= 4(n+1)/9,

whilst 

    for  μ=1:   γ_2n  (1)= 2(2n+1)/27,   γ_2n +1(1)= 2(2n+1)/9.


For each μ>-1, we have 

    γ_2n+1(μ) = 4n/9 +  2(5-2μ)/27𝒪(1)
    		  whilst  γ_2n+2(μ) = 4n/27  +  2(7+2μ)/81𝒪(1)       n→+∞.
 
In the light of Theorem <ref>,   an upper bound for the  absolute value of the largest zero |x_n,n| of P_n is given by (<ref>) with c=4/9 and α=1 so that we obtain 

    |x_n,n| ⩽  3^1/3 n^1/3 + o(n^1/3),      n⩾ 1.
 
The accuracy of the result is illustrated in Figure <ref>, where we have only plotted the positive zeros of P_3n(x), but similar results hold for the zeros of P_3n+1(x) and P_3n+2(x), which are just rotations of the positive zeros.




The zeros of three consecutive polynomials interlace and lie on the three-starlike set S : Figure <ref> illustrates this. 


 In the light of Theorem <ref>, the two orthogonality weights can be expressed via the confluent hypergeometric function of the second kind 𝐔(a,b;x), which  admits the following integral representation <cit.>

    𝐔(a,b;x) = 1/Γ(a)∫_0^∞ t^a-1 (t+1)^-a+b-1^-t x  t ,
 
provided that (a)>0 and | (x) |< π/2, whilst 

    𝐔(0,b;x) = 1,
 
and one has the identity  𝐔(a,b;x) = x^1-b𝐔(a-b+1,2-b;x). 



 Let μ>-1. The threefold symmetric polynomial sequence {P_n(·:μ)}_n⩾ 0 defined by the recurrence relation (<ref>), with γ_n given by (<ref>), is 2-orthogonal with respect to (u_0,u_1) admitting the integral representation  (<ref>) where 

    𝒰_0(x) 𝒰_0(x;μ) 
    		 = 3  Γ(μ+23)/Γ(13)Γ(23) ^-x^3𝐔(μ3,23;  x^3) ,
    	 
       𝒰_1(x)𝒰_1(x;μ) 
    		 =  9    Γ(μ +5/3) /Γ(1/3) Γ(2/3)
    			   	x^2  e^-x^3 𝐔(μ3+1,53,x^3) ,

and b=+∞, 
where 𝐔 represents the Kummer (confluent hypergeometric) function of second kind.




 Under the assumptions, it follows that u_0 and u_1 satisfy (<ref>) and {P_n(·:μ)}_n⩾ 0 satisfies the recurrence relation (<ref>) with γ_n>0 for all n⩾1. 
In the light of Theorem <ref>, there exists a function 𝒲_0 and a domain 𝒞 containing all the zeros of {P_n}_n⩾0 so that 

    ⟨ u_0, f(x) ⟩  = ∫_𝒞 f(x) 𝒲_0(x)  x ,

is valid for every polynomial f. By virtue of Theorem <ref> and the asymptotic behavior of γ_n for large n, the curve that contains all the zeros of {P_n}_n⩾0 corresponds to the starlike set S in Fig. <ref>. 
According to Theorem <ref>,  u_0 and u_1 admit the representation  (<ref>),  provided that there exist two twice differentiable functions 𝒰_0 and 𝒰_1 from ℝ  to ℝ   that are solutions to 

    1/3𝒰_0”(x)+ x^2 𝒰_0'(x) -  (μ -2) x 𝒰_0(x)   = λ_0 g_0(x)
    		    if  μ>-1,
    		
    
    𝒰_1(x)
    	=	 2/3𝒰_0'(x)+2 x^2 𝒰_0(x) + λ_1 g_1(x)     if  μ≠ 0  and μ>-1,
    	
    
    	x 𝒰_1^'(x) = 2 𝒰_0^'(x) +λ_1 g_1(x) 
    	   if  μ=0,

and satisfying

    lim_x→ω^j∞ f(x) ^l / x^l𝒰_0(x) = 0  
    		,   for    j,l∈{0,1,2}  and   f∈𝒫,
          ∫_0^∞𝒰_0(x)  x=1 ,

where λ_k ∈ℂ (possibly zero) and g_k≠0 are rapidly decreasing functions, locally integrable, representing the null functional (k=0,1). For λ_0=0, the general solution of the first equation in (<ref>) can be written as

    y(x) = c_1   112-μ/32/3 -x^3  +  c_2   x  111-μ/34/3 -x^3  ,

with c_1,c_2 two arbitrary constants.  Based on  <cit.> we have 

    ^-x^3𝐔(μ3,23,x^3) = Γ(1/3)/Γ(μ+1/3)112-μ/32/3 -x^3  
    	+  μΓ(-1/3)/3Γ(μ/3+1)  x  111-μ/34/3 -x^3  ,

which is valid when b is not an integer. Observe that (see <cit.>)

    G(x;μ) = ^-x^3𝐔(μ3,23,x^3)∼^-x^3 x^-μ∑_s=0^∞(μ3)_s(μ+13)_s/s!(-z)^-s ,

and also that <cit.> 

    / x G(x;μ) = -3x^2  ^-x^3𝐔(μ3,53,x^3),

so that, for both l=0,1, we obtain lim_x→∞ f(x) ^l G(x;μ) / x^l = 0 for every polynomial f. 

According to <cit.>, the following identity 

    ∫_0^+∞  t^b-1𝐔(a,c; t) ^-t  t 
    	= Γ(b)Γ(b-c+1)/Γ(a+b-c+1),

holds, provided that  (b)>max(0, (c)-1). Therefore

    ∫_0^+∞ x^n ^-x^3𝐔(μ3,23; x^3)  x
    	=1/3 ∫_0^+∞  t^n+13-1 ^-t𝐔(μ3,23;t)   t
    	= 1/3 Γ(n+13)Γ(n+23 )/Γ(n+μ+2/3) ,

for n⩾ 0, which, upon the substitution n→ 3n, becomes 

    ∫_0^+∞ x^3n^-x^3𝐔(μ3,23;  x^3)  x
    	 = Γ(13)Γ(23) /3  Γ(μ+23)  (13)_n (23 )_n/  (μ+23)_n    ,      n⩾ 0.


Consequently, the particular choices of 

    c_1 = 3  Γ(μ+23)/Γ(23)Γ(μ+1/3)  and   
    	c_2 
    	 = -3μΓ(μ+23) /Γ(μ/3+1)Γ(13)

in the general solution (<ref>) gives the function 
𝒰_0(x;μ) in (<ref>), which 
meets the requirements (<ref>)-(<ref>). Furthermore, we have 

    (u_0)_n = (∫_Γ_0 x- ω^2  ∫_Γ_1 x-ω∫_Γ_2 x)( x^n  Γ(μ+23)/Γ(13)Γ(23) ^-x^3𝐔(μ3,23;  x^3)),
 
which matches with  (<ref>). 

Now, for the integral representation of the linear functional u_1 we consider (<ref>) with λ_1=0 
and this gives 

    [                                                        𝒰_1(x;μ)                                                               =                              - (μ +2)/μ( 𝒰_0' (x)+ 3x^2 𝒰_0(x));                                                                                                                               = 9    Γ(μ +5/3) /Γ(1/3) Γ(2/3)
    			   	x^2  e^-x^3 𝐔(μ3+1,53,x^3);                                                                 ]   if   μ≠ 0, 
    𝒰_1(x;μ)  = 9    Γ(5/3) /Γ(1/3) Γ(2/3)Γ( 2/3,x^3)
    	   if  μ=0,
 
where, for the first line of the latter identity we have used <cit.> as well as the contiguous relation <cit.>. 
The second line of the latter identity corresponds to the first line after taking μ=0 because of <cit.>, so that we obtain (<ref>). 
Here, Γ(α,z) is the incomplete Gamma function: Γ(α,z) = ∫_z^+∞ t^α-1^-t  t provided that α>0. 
Moreover, based on (<ref>), we have 

    ∫_0^+∞ x^3n+1𝒰_1(x;μ)  x = (2/3)_n (4/3)_n/(μ+5/3)_n, 
    	      μ>-1,

which coincides with the moments of u_1 given in (<ref>).




  §.§.§ Particular cases
 

For the following choices of the parameter μ the expressions  (<ref>)-(<ref>)  relate to other functions. Namely we have: 



  * μ=-1/2, then 

    𝒰_0(x;-1/2) =  3/2√(3/π) e^-x^3 U(-1/6,2/3,x^3)
      and  𝒰_1(x;-1/2) =9 √(3) e^-x^3/2 x K_1/3(x^3/2)/4 π,

where K_ν is the modified Bessel function of the second kind (see <cit.>). 


  * μ=0, then 

    𝒰_0(x;0) 
    		 = 3   /Γ(13)  ^-x^3  and  𝒰_1(x;0) 
    		 =  9    Γ(5/3) /Γ(1/3) Γ(2/3)Γ(2/3;x^3).



  * μ=1, then

    𝒰_0(x;1) =  3 √(3) e^-x^3/2√(x) K_1/6(x^3/2)/2
       π ^3/2  and  𝒰_1(x;1) = 9 √(3) e^-x^3 U(2/3,1/3,x^3)/2 π.

 

  * μ=2, then  

    𝒰_0(x;2) =  √(3)Γ(1/3) /2 πΓ(13,x^3)
      and  𝒰_1(x;2) =2 √(3)Γ(1/3) /πΓ(23,x^3).



  





  §.§.§ The differential equation



Following Lemma <ref>, the polynomial P_n(x;μ) is a solution of the differential equation 

    -2/3y”'(x)
    	+  2 x^2 y”(x)+ 2 x (μ +3/4((-1)^n+3)-n/2) y'(x)
    	=2n(μ +n/2+3 (-1)^n/4+5/4) y(x) ,
 
whose general solution can be written as 

    y(x) 
       =   	c_1  22-n/3,n/6+(-1)^n/4+μ/3+5/121/3,2/3x^3
    	+ 	c_2    x  221/3-n/3,n/6+(-1)^n/4+μ/3+3/42/3,4/3x^3
          	+ c_3   x^2  222/3-n/3,n/6+(-1)^n/4+μ/3+13/124/3,5/3x^3 .


For each positive integer n, there is only one (monic) polynomial solution, but it is not always the same: it depends on whether n equals 0,1 or 2 3. 
To be precise, we have 

    P_3n(x;μ)        P_n^[0](x^3;μ) 
       =   (-1)^n  (1/3)_n (2/3)_n/(n/2+1/4 (-1)^3 n+μ/3+5/12)_n22-n,n/2+(-1)^n/4+μ/3+5/121/3,2/3x^3 , 
    
    	P_3n+1(x;μ)       x P_n^[1](x^3;μ)
        =     x  (-1)^n   (2/3)_n (4/3)_n/(n/2+1/4 (-1)^n+1+μ/3+11/12)_n22-n,μ/3+n/2+1/4 (-1)^n+1+11/122/3,4/3x^3, 
    
    	P_3n+2(x;μ)       x^2 P_n^[2](x^3;μ) 
       =    x^2 (-1)^n  (4/3)_n (5/3)_n /(n/2+1/4 (-1)^n+2+μ/3+17/12)_n22-n, μ/3+n/2+1/4 (-1)^n+2+17/124/3,5/3x^3 .

 These polynomial sequences {P_n^[k](·;μ)}_n⩾ 0, with k∈{0,1,2}, are precisely the 2-s in the cubic decomposition of {P_n(·;μ)}_n⩾ 0. In fact, from Lemma <ref>, these three 2-s are not threefold symmetric and satisfy the recurrence relation (<ref>),  whose recurrence coefficients are given in the Appendix for completeness.  These coefficients have been computed in  <cit.>, for a different choice of the "free" parameter γ_1. We have included them here for a matter of completeness. 







  §.§.§ The sequence of derivatives


Following Theorem <ref> along with Theorem <ref>,  {Q_n(·;μ)1/n+1 P_n+1'(x;μ)}_n⩾ 0 is 2-orthogonal for the vector functional (v_0,v_1), which admits an integral representation via two weight functions (𝒲_0(·;μ),𝒲_1(·;μ)) with support on the three-starlike set S which are  given by 

    [ [ 𝒲_0(x;μ); 𝒲_1(x;μ) ]]
    	 =[[ μ+2/μ +1  -x/μ +1;        0        1;          ]]
    [ [ 𝒲_0(x;μ); 𝒲_1(x;μ) ]].

Proposition <ref> allows us to conclude 

    ⟨ v_k,f ⟩ 
    	=  (∫_Γ_0 x- ω^2  ∫_Γ_1 x-ω∫_Γ_2 x)( f(x)  𝒱_k(x;μ)),     ∀ f∈𝒫,    k=0,1,

where 

    𝒱_0(x;μ)    =   μ+2/μ +1𝒰_0(x;μ)  -x/μ +1𝒰_1(x;μ)  =  𝒰_0(x;μ+3) ,
    𝒱_1(x;μ)    =   𝒰_1(x;μ).

Furthermore,  {Q_n(·;μ)}_n⩾ 0 satisfies the second order recurrence relation (<ref>) where 

    [                              γ_2n γ_2n (μ) =2/3 n (2 n+1)/(3 n+2+μ) ,      n⩾ 1,; γ_2n+1 γ_2n +1(μ) =2/3(n+1) (2 n+1) (n +2+μ) /(3 n+2+μ) (3 n+5+ μ)  ,      n⩾ 0.;                                                                                  ]


For each integer n⩾ 0, the polynomial Q_n(x;μ) is a solution to the differential equation 

    -2/3 y^(3)(x)+2 x^2 y”(x) -1/2 x (-4 μ +3 ((-1)^n-5)+2 n) y'(x)
    	= 1/2 n (4 μ +2   n-3 (-1)^n+11) y(x),

from which we deduce 

    Q_3n(x;μ)    =   (-1)^n (1/3)_n (2/3)_n/(n/2-(-1)^n/4 +μ/3+11/12)_n 22-n,n/2-(-1)^n/4 +μ/3+11/121/3,2/3x^3, 
    
    	Q_3n+1(x;μ)    =    x  (-1)^n (2/3)_n (4/3)_n/(μ/3+n/2+(-1)^n/4 +17/12)_n 22-n,μ/3+n/2+(-1)^n/4 +17/122/3,4/3x^3 ,
     
    	Q_3n+2(x;μ)    =    x^2  (-1)^n (4/3)_n (5/3)_n/(μ/3+n/2-(-1)^n/4 +23/12)_n 22-n,μ/3+n/2-(-1)^n/4 +23/124/3,5/3x^3
    		,     n⩾ 0.




 It turns out that the forthcoming case B_2 is closely connected to this case B_1: the sequence of derivatives in case B_1 belongs to case B_2. 




 §.§ Case  B_2
  


For this family  ϑ_2n+1=1 and we can write ϑ_2n=n+ρ+1/n+ ρ, 
under the relation ϑ_2 = ρ +2/ρ +1 with ρ≠-n.  As a result, the threefold symmetric Hahn-classical sequence 
{P_n(·;ρ;γ_1)}_n⩾ 0 satisfies the recurrence (<ref>) with 

    [                     γ_2n  =n(2n+1) (ρ+3)/(3n+ρ)γ_1,      n⩾ 1,; γ_2n +1 = (n+1)(2n+1)(n+ρ) (ρ+3)/(3n+ρ)(3n+ρ+3)γ_1,      n⩾ 0. ]

Analogously to the former cases, the parameter γ_1 is redundant for the study, and therefore we can choose a representative value for γ_1. 
Here we set γ_1 = 2/3(ρ+3). Hence the 2- {P_n(·;ρ) P_n(·;ρ;2/3(ρ+3))}_n⩾ 0 satisfies (<ref>)  with 

    [                     γ_2n γ_2n (ρ)  =2n(2n+1) /3(3n+ρ) ,      n⩾ 1,; γ_2n +1γ_2n+1 (ρ)= 2(n+1)(2n+1)(n+ρ) /3(3n+ρ)(3n+ρ+3) ,      n⩾ 0. ]

Furthermore, the 2- {Q_n(x;ρ)1/n+1P_n+1'(x;ρ)}_n⩾ 0 satisfies (<ref>) with 

    [ γ_2n γ_2n (ρ) =  2 n (2 n+1)  (n+ρ +1)/3(3 n+ρ ) (3 n+ρ +3),      n⩾ 1,;                γ_2n +1γ_2n (ρ)=2 (n+1) (2 n+1) /3(3 n+ρ +3),      n⩾ 0. ]


Observe that the recurrence coefficients uniquely determine a 2-. A comparison between (<ref>) and (<ref>) 
and between (<ref>) and (<ref>) shows that 

    Q_n^case B_2(x;μ) = P_n^case B_1(x;μ+1) ,     for all  n⩾ 0,

while 

    Q_n^case B_1(x;μ) = P_n^case B_2(x;μ+2) ,     for all  n⩾ 0,

 where the notation P_n^case B_1 and P_n^case B_2 is used for the threefold 2-s defined by the recurrence coefficients (<ref>), in case B_1, and  (<ref>), in case B_2, respectively. Likewise, Q_n^case B_1 and Q_n^case B_1 denote the monic derivatives of P_n^case B_1 and P_n^case B_2, respectively. 

From this, we conclude that the  the 
polynomial sequences {Q_n^case B_1(x;μ)1/n+1/ xP_n+1^case B_1(x;μ)}_n⩾ 0
and {Q_n^case B_2(x;ρ)1/n+1/ xP_n+1^case B_2(x;ρ)}_n⩾ 0 are also Hahn-classical. Furthermore, we have 

    1/(n+2)(n+1)^2/ x^2P_n+2(x;μ)
    	= P_n(x;μ+3) ,

in both cases B_1 and B_2.
These observations also serve as an alternative to the proof of the following result:

  Let ρ>0. The threefold symmetric polynomial sequence {P_n(·: ρ)}_n⩾ 0 defined by the recurrence relation (<ref>) with γ_n given by (<ref>) is 2-orthogonal with respect to (u_0,u_1) admitting the integral representation  (<ref>) where 

    𝒰_0(x) 𝒰_0(x; ρ) 
    		= 3  Γ(ρ +3/3)/Γ(1/3) Γ(2/3)^-x^3𝐔(ρ +1/3,2/3,x^3)	,
    
    
          𝒰_1(x)𝒰_1(x; ρ) 
    	=9  Γ(ρ +3/3) /Γ(1/3) Γ(2/3)
    	x^2 ^-x^3𝐔(ρ +1/3,5/3,x^3),

and b=+∞, 
where 𝐔 represents the Kummer (confluent hypergeometric) function of second kind.

 The proof is entirely analogous to the proof of Proposition <ref>.


In fact, we have 

    𝒰_0^case B_2(x; ρ) 
    	= 𝒰_0^case B_1(x; ρ+1) 
    	  and  𝒰_1^case B_2(x; ρ) 
    	= 𝒰_1^case B_1(x; ρ-2),

provided that ρ>1. 





 §.§ Case C
   

 Here we have 

    ϑ_2n-1 = n+μ+1/n+μ  and  ϑ_2n =  n+ρ+1/n+ρ ,      n⩾ 1.
 
so that  (<ref>)  reads as  

    [       γ_2n  =  2n(2n+1)(n+μ) /(3n+μ-1)(3n+μ+2)(3n+ρ) (μ+2)(ρ+3)γ_1/2,      n⩾ 1,; γ_2n +1 = 2(n+1)(2n+1)(n+ρ)  /(3n+μ+2)(3n+ρ)(3n+ρ+3) (μ+2)(ρ+3)γ_1/2,      n⩾ 0.;                                                                                  ]
 

The Hahn-classical polynomial sequence obtained under these assumptions clearly depends  on the pair of parameters (μ,ρ)  and on γ_1, so it makes sense to incorporate this information, and therefore we shall refer to it as {P_n(·;μ,ρ;γ_1)}_n⩾ 0. 
Similar to the precedent cases, without loss of generality, it suffices to study the sequence for a particular choice of γ_1. 
A scaling of the variable would then reproduce all the other sequences {P_n(·;μ,ρ;γ_1)}_n⩾ 0 within this equivalence class. Hence, 
we set γ_1 = 2 (μ+2)(ρ+3), and the analysis is on  

    {P_n(·;μ,ρ) P_n(·;μ,ρ; 2 (μ+2)(ρ+3))}_n⩾ 0,
 
which satisfies the recurrence relation (<ref>) where the γ-coefficients are given by

    [       γ_2n γ_2n  (μ,ρ) =  2n(2n+1)(n+μ) /(3n+μ-1)(3n+μ+2)(3n+ρ) ,      n⩾ 1,; γ_2n +1γ_2n  (μ,ρ)= 2(n+1)(2n+1)(n+ρ)  /(3n+μ+2)(3n+ρ)(3n+ρ+3)  ,      n⩾ 0.;                                                                              ]

Observe that γ_n+1>0 for all n⩾ 0 provided that μ+1, ρ>0 and that

    γ_n = 4/27 +o(1) ,   as   n→ +∞.
 
Theorem <ref> ensures that the largest zero of P_n(x;μ,ρ) in absolute value is always less than or equal to 1. Figure <ref> illustrates the curve defined by these zeros. 




Regarding the orthogonality measures for the polynomial sequence {P_n(·;μ,ρ)}_n⩾ 0  we have: 

 Let μ>-1 and ρ>0. The threefold symmetric polynomial sequence {P_n(·;μ, ρ)}_n⩾ 0 defined by the recurrence relation (<ref>) with γ_n given by (<ref>)  is 2-orthogonal with respect to (u_0,u_1), admitting the integral representation  (<ref>) with b=1 and  

    𝒰_0(x) 𝒰_0(x;μ, ρ) 
    		=  3Γ(μ+2/3)Γ(ρ/3+1)/Γ(1/3)Γ(2/3)Γ(μ +ρ +2/3)  
    	(1-x^3)^μ+ρ-1/3 μ/3, ρ+1/3μ+ρ+2/31-x^3, 
    	 
          𝒰_1(x)𝒰_1(x; μ,ρ) 
    	= 3 Γ(μ+5/3)Γ(ρ/3+1)/Γ(2/3)Γ(4/3)Γ(μ +ρ +2/3)  x^2 (1-x^3)^μ+ρ-1/3μ/3+1, ρ+1/3μ+ρ+2/31-x^3.
 



 According to Theorem <ref>, {P_n(·;μ, ρ)}_n⩾ 0 is 2-orthogonal with respect to the vector of linear functionals (u_0,u_1) fulfilling the distributional equations 

    (  1- x^3)u_0 ” +  x^2 (μ +ρ -4) u_0 '-  (μ -2) (ρ -1) x u_0  =0, 
    	
       μ/(μ +2)   u_1 
    	   =  ( x^3-1)u_0 ' 
    	-  (ρ -1) x^2  u_0    if  μ≠ 0, 
    
    	x u_1' = 2 u_0'	
    		   if  μ= 0,

which yields 

    [ (u_0)_3n = (1/3)_n (2/3)_n  /(μ +2/3)_n (ρ +3/3)_n
    	,   
    	(u_0)_3n+1 =  (u_0)_3n+2 =0,      n⩾0, ]
 as well as  

    [ (u_1)_3n+1 = (2/3)_n (4/3)_n  /(μ +5/3)_n (ρ +3/3)_n ,
    	   
       	(u_1)_3n =  (u_1)_3n+2 =0,      n⩾0, ]

because (u_0)_0=(u_1)_1=1 and (u_0)_1=(u_0)_2=(u_1)_0=(u_1)_2=0.

Since lim_n→∞γ_n = 4/27, then by virtue of Theorem <ref>, u_0 and u_1  admit the integral representation 
(<ref>) with b=1, provided that there exists a pair of functions 𝒰_0 and 𝒰_1 with support on [0,1) satisfying (<ref>)-(<ref>) and subject to the boundary condition (<ref>). To be precise, under the assumptions, we seek 𝒰_0(x;μ,ρ) and 𝒰_1(x;μ,ρ) defined on defined on [0,1) such that 

    (  1- x^3)𝒰_0”(x;μ,ρ)+  x^2 (μ +ρ -4) 𝒰_0'(x;μ,ρ)-  (μ -2) (ρ -1) x 𝒰_0(x;μ,ρ)  = λ_0 g_0(x),
    		 
       μ/(μ +2)  𝒰_1(x;μ,ρ)
    	   =  ( x^3-1)𝒰_0'(x;μ,ρ) 
    	-  (ρ -1) x^2 𝒰_0(x;μ,ρ) + λ_1 g_1(x)     if  μ≠ 0, 
    
    	x 𝒰_1^'(x;μ,ρ)= 2𝒰_0^'(x;μ,ρ)+ λ_1 g_1(x)   if  μ= 0,

where λ_j, for j∈{0,1}, are constants (possibly zero) and g_j are functions representing the null linear functional,  i.e., 

    ∫_𝒞 g_j(x) x^n  x = 0,      n⩾ 0,

and, in addition, 

    lim_x→ 1 f(x) ^l / x^l𝒰_0(x;μ,ρ) = 0  
    		,   for any  l∈{0,1}  and   f∈𝒫,
          ∫_0^1 x^3n𝒰_0(x;μ,ρ)  x
    		=  (1/3)_n (2/3)_n  /(μ +2/3)_n (ρ +3/3)_n
    		,     n⩾ 0.


For λ_0=0, the differential equation (<ref>) becomes 

    (  1- x^3) 𝒰_0”(x;μ,ρ)+  x^2 (μ +ρ -4) 𝒰_0'(x;μ,ρ)-  (μ -2) (ρ -1) x 𝒰_0(x;μ,ρ)  =0 ,
 
whose general solution can be written as 

    𝒰_0(x;μ,ρ) = c_1    2-μ/3,1-ρ/32/3x^3 
    	+ c_2   x  1-μ/3, 2-ρ/34/3x^3 ,

with integration constants c_1 and c_2 which must be chosen such that (<ref>)-(<ref>) hold.  
By taking 

    c_2=  -Γ(2/3) Γ(μ +1/3)
       Γ(ρ +2/3)/Γ(4/3) Γ(μ/3) Γ(ρ +1/3) c_1,
 
we have 

    𝒰_0(x;μ,ρ)    =    c_1Γ(μ+1/3) Γ(ρ+2/3)/Γ(1/3)
          (   Γ(4/3)/Γ(μ+1/3) Γ(ρ+2/3) 2-μ/3,1-ρ/32/3x^3 
    	- Γ(2/3)  x /Γ(μ/3) Γ(ρ+1/3) 1-μ/3, 2-ρ/34/3x^3)
       =    c_1Γ(μ+1/3) Γ(ρ+2/3)/Γ(1/3)
    	(1-x^3)^μ+ρ-1/3 μ/3, ρ+1/3μ+ρ+2/31-x^3 ,

where, for the latter identity, we have used <cit.>. In fact, for this choice of c_2, it follows that  lim_x→ 1^-𝒰_0(x;μ,ρ) =0, because 

    21a,bc1= Γ(c) Γ(c-a-b)/Γ(c-a)Γ(c-b)
 
is valid whenever (c-a-b)>0. Moreover, condition (<ref>) is fulfilled, and, in addition, we successively have 

    ∫_0^1 x^3n𝒰_0(x;μ,ρ)  x 
    	   =    c_1Γ(μ+1/3) Γ(ρ+2/3)/Γ(1/3)∫_0^1 x^3n (1-x^3)^μ+ρ-1/3 μ/3, ρ+1/3μ+ρ+2/31-x^3  x
       =    c_1Γ(μ+1/3) Γ(ρ+2/3)/3 Γ(1/3)∫_0^1 (1-x)^n-2/3 x^μ+ρ-1/3 μ/3, ρ+1/3μ+ρ+2/3x ,  x
       =    c_1Γ(2/3)Γ(μ+1/3) Γ(ρ+2/3)Γ(μ +ρ +2/3)/3  Γ(μ+2/3)Γ(ρ/3+1)(1/3)_n (2/3)_n /(μ+2/3)_n (ρ/3+1)_n,

where, for the last identity, we have used <cit.>. If we take

    c_1 = 3  Γ(μ+2/3)Γ(ρ/3+1)/Γ(2/3)Γ(μ+1/3) Γ(ρ+2/3)Γ(μ +ρ +2/3)

then (<ref>) is fulfilled. As a result, the first orthogonality measure can be represented as in (<ref>) (case k=0) where 𝒰_0 is given by (<ref>). 

Now, in order to obtain the second orthogonality measure, we take λ_1=0 in (<ref>), which involves a derivative of 
𝒰_0(x;μ,ρ), that, according to <cit.>, corresponds to 

    𝒰_0'(x;μ,ρ) = -  3Γ(μ+2/3)Γ(ρ/3+1)/Γ(1/3)Γ(2/3)Γ(μ +ρ +2/3)  (μ+ρ-1)
    	x^2 (1-x^3)^μ+ρ+2/3-2
    ×μ/3, ρ+1/3μ+ρ+2/3 -11-x^3 .
 
For μ≠0, it follows from (<ref>) that 

    μ/(μ +2)  𝒰_1(x;μ,ρ)
    	=   9Γ(μ+2/3)Γ(ρ/3+1)/Γ(1/3)Γ(2/3)Γ(μ +ρ +2/3)  x^2 (1-x^3)^μ+ρ-1/3
           ( (μ+ρ-13)
    	  μ/3, ρ+1/3μ+ρ-1/31-x^3
    	 -  (ρ-13) μ/3, ρ+1/3μ+ρ+2/31-x^3)
    
    	=    3μΓ(μ+2/3)Γ(ρ/3+1)/Γ(1/3)Γ(2/3)Γ(μ +ρ +2/3)  x^2 (1-x^3)^μ+ρ-1/3μ/3+1, ρ+1/3μ+ρ+2/31-x^3 ,

where, for the last identity we have used <cit.>

    (c-1)   21a,bc-11-z-(-a+c-1)   21a,bc1-z
    	= a    21a+1,bc1-z.
 
As a result, we obtain the expression (<ref>) for 𝒰_1(x;μ,ρ) and we have 

    ∫_0^1 x^3n+1𝒰_1(x;μ,ρ)  x 
    	= (2/3)_n (4/3)_n  /(μ +5/3)_n (ρ +3/3)_n
    	= (u_1)_3n+1,      n⩾ 0.
 
When μ=0, then (<ref>) reads as 

    𝒰_1^'(x;0,ρ)
    	= -18  Γ(ρ/3+1) /Γ(1/3) Γ(ρ -1/3) x (1-x^3)^ρ -4/3 ,
 
so that 

    𝒰_1 (x;0,ρ)
    	 = K -9 x^2 Γ(ρ/3+1) /Γ(1/3) Γ(ρ -1/3)212/3,4-ρ/35/3x^3,
 
for some integration constant K. 
We set K= 6Γ(ρ/3+1)Γ(2 /3)/Γ(1 /3)Γ(ρ+1 /3), use 
identity <cit.> to obtain 

    𝒰_1 (x;0,ρ)
    	 = 6Γ(ρ/3+1)/Γ(1/3)Γ(ρ +2/3)  x^2 (1-x^3)^ρ-1/3 1, ρ+1/3ρ+2/31-x^3,
 
which is the same as (<ref>) with μ=0. 





  §.§.§ Differential equation and the cubic decomposition of {P_n(x;μ,ρ)}_n⩾ 0


Let us assume μ>-1 and ρ>0. For each positive integer n, the polynomial P_n(x;μ,ρ) satisfies the differential equation

    [                                               (x^3-1) y^(3)(x)+x^2 (μ +ρ +5) y”(x); + 1/8 x (8 μρ +14 μ -6 n^2-4 n (μ +ρ +2)-3 (-1)^n (2 μ -2 ρ +1)+18
       ρ +27) y'(x);                       = 1/16 n (4 μ +2 n+3 (-1)^n+5) (2 n-3 (-1)^n+4   ρ +3) y(x). ]
 
The general solution of the latter equation can be written as 

    [                                                                    P_n(x;μ,ρ)           = c_1  32-n/3,n/6+(-1)^n/4+μ/3+5/12,n/6-(-1)^n/4 +ρ/3+1/41/3,2/3x^3;                                                                                   +   c_2 x  321/3-n/3,n/6+(-1)^n/4+μ/3+3/4,n/6-(-1)^n/4+ρ/3+7/122/3,4/3x^3;                                                                               + c_3 x^2  322/3-n/3,n/6+(-1)^n/4+μ/3+13/12,n/6-(-1)^n/4+ρ/3+11/124/3,5/3x^3. ]
 
which has only one polynomial solution for each n⩾ 0. Similar to the precedent cases, for each positive integer n, there is only one (monic) polynomial solution and it differs depending on whether n equals 0,1 or 2 3. To be precise, we have 

    P_3n(x;μ,ρ)       P_n^[0](x^3;μ,ρ) 
       =   (-1)^n (1/3)_n (2/3)_n/(n/2+1/4 (-1)^3
       n+μ/3+5/12)_n (n/2-1/4 (-1)^3 n+ρ/3+1/4)_n
          · 32-n,n/2+1/4 (-1)^3 n+μ/3+5/12,n/2-1/4 (-1)^3 n+ρ/3+1/41/3,2/3x^3 , 
    
    	 P_3n+1(x;μ,ρ)       x P_n^[1](x^3;μ,ρ) 
       =     x  (-1)^n (2/3)_n (4/3)_n  /(n/2-1/4 (-1)^3
       n+μ/3+11/12)_n (n/2+1/4 (-1)^3 n+ρ/3+3/4)_n
          · 32-n,n/2-1/4 (-1)^3 n+μ/3+11/12,n/2+1/4 (-1)^3 n+ρ/3+3/42/3,4/3x^3 ,
       	 
    
    	P_3n+2(x;μ,ρ)      x^2 P_n^[2](x^3;μ,ρ) 
       =    x^2 (-1)^n (4/3)_n (5/3)_n  /(n/2+1/4 (-1)^3
       n+μ/3+17/12)_n (n/2-1/4 (-1)^3 n+ρ/3+5/4)_n
          · 32-n,n/2+1/4 (-1)^3 n+μ/3+17/12,n/2-1/4 (-1)^3 n+ρ/3+5/44/3,5/3x^3 .

  
  Here, the sequences {P_n^[k](·;μ,ρ)}_n⩾ 0, with k∈{0,1,2}, are precisely the 2-s in the cubic decomposition of {P_n(·;μ,ρ)}_n⩾ 0. From Lemma <ref>, these three 2-s are not threefold symmetric and satisfy the recurrence relation (<ref>). 
These coefficients have been computed in  <cit.>, for a different choice of the "free" parameter γ_1. We have included them in the Appendix  
for completeness, where we have used the software Mathematica. 



  §.§.§ The sequence of derivatives


The 2- {Q_n(x;μ,ρ)1/n+1 P_n+1'(x;μ,ρ)}_n⩾ 0 satisfies the relation (<ref>) with 

    [           γ_2n γ_2n (μ,ρ)
    				= 2 n (2 n+1) (n+ρ +1)/(μ +3 n+2) (3 n+ρ ) (3 n+ρ +3)
    				 ,      n⩾ 1,; γ_2n +1γ_2n +1 (μ,ρ)
    				= 2 (n+1) (2 n+1) (μ +n+2)/(μ +3 n+2) (μ +3 n+5) (3 n+ρ +3)
    				,      n⩾ 0,;                                                                                                       ]

whose expressions are derived from (<ref>) and (<ref>). 
A straightforward comparison between (<ref>) and (<ref>) readily shows that 

    γ_n+1 (μ,ρ) = γ_n+1 (ρ+1,μ+2),      n⩾ 0.

This, combined with the uniqueness of a polynomial sequence defined by the recurrence relation  (<ref>) implies 
Q_n(x;μ,ρ) = P_n(x;ρ+1,μ+2)
which means  

    1/n+1 P_n+1'(x;μ,ρ)=P_n(x;ρ+1,μ+2)
    	,      n⩾ 0.




  §.§.§ Particular cases
 

The so called Humbert polynomials introduced in <cit.> correspond to {P_n(x;3ν-1/2,3ν/2)}_n⩾ 0, up to a scaling of the variable. In fact, by setting μ=3ν-1/2 and ρ=3ν/2, this 2- 
satisfies

    P_n+2(x;3ν-12,3ν2) 
    
    	= x P_n+1(x;3ν-12,3ν2)  - 4/27n (n+1) (3 ν +n-1)/(ν +n-1) (ν +n) (ν +n+1)P_n-1(x;3ν-12,3ν2) ,

with initial conditions P_0=1,  P_1(x)=x and P_2(x)=x^2. Baker <cit.> gave explicit expressions for this particular case via _3F_2 hypergeometric series and these coincide with those obtained above. Pincherle polynomials, introduced in <cit.>, are a particular case obtained with ν=1/2, which received special attention. For instance, they were discussed in <cit.> and several properties were further analyzed in <cit.> with the focus of the study on the algebraic properties, including generating functions, as well as the analysis of the components arising in the cubic decomposition, but integral representations for the orthogonality measures were not included. 
Prior to this work, Douak and Maroni <cit.>  analyzed in detail the case where ν=1 (and therefore μ=1 and ν=3/2), which gives a threefold symmetric sequence with constant γ-coefficients, namely 

    P_n+2(x;1,32) = x P_n+1(x;1,32)   - 427 P_n-1(x;1,32) ,

with the same initial conditions. In other words, this means that the associated sequence coincides with the original one, and, regarding the nature of the problem, they referred to this polynomial sequence as Chebyshev type polynomials. Therein, an integral representation for the orthogonality functionals was given, with support on an interval on the real line (which obviously does not contain all the zeros of the polynomial sequence  <cit.>), which is, from our point of view, a bit artificial. Following Proposition <ref>, the polynomial sequence {P_n(x;1,32)}_n⩾0 is 2-orthogonal with respect to (u_0,u_1) which admit the integral representation  (<ref>) with b=1 and  

    𝒰_0(x;1, 32) 
    
    
    
    	=  9√(3)/4π( ( 1+√(1-x^3))^1/3- ( 1- √(1-x^3))^1/3), 
    	
       𝒰_1(x; 1, 32) 
    
    
    
    		  =  27 √(3)/8 π((√(1-x^3)+1)^2/3-(1-√(1-x^3))^2/3),

after using <cit.>. 


Another particular case that made its appearance in the literature corresponds to the case where μ=-1/2 and ρ=0, and therefore γ_n+2=γ_n=4/27 and γ_1=4/9.  In fact, {P_n(x;-1/2,0)}_n⩾ 0  is a particular case of Faber polynomials, named after the author of <cit.>, and this polynomial sequence with almost constant coefficients has been studied in <cit.> and <cit.>. The latter paper was essentially devoted to the study of the zeros, while the former was mainly dedicated to the algebraic properties as well as a representation of the measure on an interval on the real line. Here, we combine the two approaches in a more general setting. So,  {P_n(x;-1/2,0)}_n⩾ 0  
is 2-orthogonal with respect to (u_0,u_1) which admit the integral representation  (<ref>) with b=1 and  

    𝒰_0(x;-12,0) 
    	=3 √(3)( (1-√(1-x^3))^1/3+(1+√(1-x^3))^1/3)/4 π√(1-x^3) , 
    	
       𝒰_1(x; -12,0) 
    		  = 9 √(3)((1-√(1-x^3))^2/3+(1+√(1-x^3))^2/3)/8 π√(1-x^3).


The penultimate and the latter particular cases are related to each other: if we recall (<ref>) we readily see that 

    1/n+1 P_n+1'(x;-12,0) = P_n(x;1,32),     n⩾ 0.
 
It turns out that, the sequence {R_n(x;-1/2,0)}_n⩾ 0 of the associated polynomials of  {P_n(x;-1/2,0)}_n⩾ 0, defined by 

    R_n(x;-12,0) ⟨ u_0, P_n(x;-12,0)-P_n(t;-12,0)/x-t⟩
 
actually coincides with {P_n(x;1,32)}_n⩾ 0, as discussed in  <cit.>. 





§ ACKNOWLEDGMENTS


We thank the anonymous referees for their careful reading and for the suggestions and corrections, which led to an improvement of the paper.  The research of WVA was supported by FWO research grant G.0864.16N and EOS project PRIMA 30889451. 





00

ABV A. I. Aptekarev, A. Branquinho, W. Van Assche, 
Multiple orthogonal polynomials for classical weights, 
Trans. Amer. Math. Soc.  355  (2003),  no. 10, 3887–3914. 
AKS A. I. Aptekarev, V. A. Kalyagin, and E. B. Saff, 
Higher-order three-term recurrences and asymptotics of multiple orthogonal polynomials, 
Constr. Approx. 30 (2009), no. 2, 175–223.
AKVI A.I. Aptekarev, V.A. Kaliaguine,  J. Van Iseghem, 
The genetic sum's representation for the moments of a system of Stieltjes functions and its application, Constr. Approx. 16 (2000), 487–524.
BDS D. Babusci, G. Dattoli, D. Sacchetti, 
The Airy transform and the associated polynomials, 
Cent. Eur. J. Phys. 9 (2011), no. 6, 1381–1386. 
Baker B.B. Baker, 
On the relation between Pincherle's polynomials and the hypergeometric function,
Proc. Edinb. Math. Soc. 39 (1920), 58–62. 
BCBR11 Y. Ben Cheikh Y, N. Ben Romdhane, 
On d-symmetric classical d-orthogonal polynomials, 
J. Comput Appl Math. 236 (2011), 85–93.
BCD00 Y. Ben Cheikh and K. Douak, 
On 2-orthogonal polynomials related to the Bateman J_n^u,v-function, 
Meth. Appl. Anal. 7 (2000), no. 4, 641–662. 
BRMG16 N. Ben Romdhane and M. Gaied, 
A generalization of the symmetric classical polynomials: Hermite and Gegenbauer polynomials, 
Integral Transforms Spec. Funct. 27 (2016), 227–244.
BR08 N. Ben Romdhane, 
On the zeros of d-symmetric d-orthogonal polynomials, 
J. Math. Anal. Appl. 344 (2008), 888–897.
CM03 P. Clarkson, E. Mansfield, 
The second Painlevé equation, its hierarchy and associated special polynomials, 
Nonlinearity 16 (2003), no. 3, R1–R26.
CouWVA J. Coussement, W. Van Assche, 
Differential equations for multiple orthogonal polynomials with respect to classical weights: raising and lowering operators, 
J. Phys. A: Math. Gen. 39  (2006),  no. 13, 3311–3318.
DelLop S. Delvaux, A. López, 
High-order three-term recursions, Riemann-Hilbert minors and Nikishin systems on star-like sets, 
Constr. Approx. 37 (2013), 383–453.
Douak96 K. Douak, 
The relation of the d-orthogonal polynomials to the Appell polynomials, 
J. Comput. Appl. Math. 70 (1996), 279–295. 
DM92 K. Douak and P. Maroni, 
Les polynômes orthogonaux classiques de dimension deux, 
Analysis 12, (1992) 71–107.
DM95 K. Douak and P. Maroni, 
Une charactérisation des polynômes orthogonaux classiques de dimension d, 
J. Approx. Theory 82 (1995), 177–204.
DM97 I K. Douak and P. Maroni, 
On d-orthogonal Tchebyshev polynomials, I, 
Appl. Numer. Math. 24, (1997), 23–53.
DM97 II K. Douak and P. Maroni, 
On d-orthogonal Tchebyshev polynomials, II, 
Meth. Appl. Anal. 4, (1997), no. 4, 404–429.
Faber G. Faber, 
Über polynomische Entwickelungen, 
Math. Annalen 57 (1903), 389–408. 
FHVA G. Filipuk, M. Haneczok, W. Van Assche, 
Computing recurrence coefficients of multiple orthogonal polynomials, 
Numer. Algorithms  70  (2015),  no. 3, 519–543. 
GR I. S. Gradshteyn and I. M. Ryzhik, 
Table of Integrals, Series, and Products, 
Edited by A. Jeffrey and D. Zwillinger, Academic Press, New York, 7th edition, 2007.
HeSaff M.X. He, E.B. Saff, 
The zeros of Faber polynomials for an m-cusped hypocycloid, 
J. Approx. Theory 78 (1994), 410–432. 
Hum P. Humbert, 
Some extensions of Pincherle's polynomials, 
Proc. Edin. Math. Soc. 39 (1920), 21–24.
HJ R. A. Horn  and C. R. Johnson, 
Matrix Analysis, 
Cambridge University Press, 1985. 
KuijlSti A. B. J. Kuijlaars, D. Stivigny,
Singular values of products of random matrices and polynomials ensembles,
Random Matrices: Theory Appl. 3 (2014), no. 3, 1450011, 22. pp.
LOHerm I. Lamiri, A. Ouni, 
d-Orthogonality of Hermite type polynomials, 
Appl. Math. Comput. 202 (2008) 24–43. 
LO I. Lamiri, A. Ouni, 
d-orthogonality of Humbert and Jacobi type polynomials, 
J. Math. Anal. Appl. 341 (2008), 24–51.
LG11 A. López García, 
Asymptotics of multiple orthogonal polynomials for a system of two measures supported on a starlike set, J. Approx. Theory 163 (2011), 1146–1184.
Mar89 P. Maroni, 
L'orthogonalité et les récurrences de polynômes d'ordre supérieur  à deux, 
Ann. Fac. Sci. Toulouse 10, (1989), no. 1, 105–139.
MarAlg P. Maroni, 
Une théorie algébrique des polynômes orthogonaux. 
Application aux polynômes orthogonaux semi-classiques, 
in `Orthogonal Polynomials and their Applications' (Erice, 1990), Eds. C. Brezinski et al.,
IMACS Ann. Comput. Appl. Math. 9, Baltzer, Basel (1991), pp. 95–130.
Mar92 P. Maroni, 
Two-dimensional orthogonal polynomials, their associated sets and the co-recursive sets, 
Numer. Algorithms 3 (1992), 299–312.
MarVar P. Maroni, 
Variations around classical orthogonal polynomials. Connected problems, 
J. Comput. Appl. Math. 48 (1993), 133–155.
Mar99 P. Maroni, 
Semi-classical character and finite-type relations between polynomial sequences, 
Appl. Numer. Math. 31 (1999), 295–330.
MarMes P. Maroni and T.A. Mesquita, 
Cubic decomposition of 2-orthogonal polynomial sequences, 
Mediterr. J. Math. 10 (2013), 843–863.
DLMF F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, and B. V. Saunders,
NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/, Release 1.0.16 of 2017-09-18. 
Pin S. Pincherle, 
Una nuova estensione delle funzione sferiche,
Mem.  Accad. R. Bologna 5 (1890), 337–362. 
SoaVal O. Vallée, M. Soares, 
Airy Functions and Applications to Physics, 
Imperial College Press, London, 2nd edition, 2010.
VAneigh11 W. Van Assche, 
Nearest neighbor recurrence relations for multiple orthogonal polynomials, 
J. Approx. Theory 163 (2011), 1427–1448.
VAY W. Van Assche, S.B. Yakubovich, 
Multiple orthogonal polynomials associated with MacDonald functions, 
Integral Transform Special Funct. 9 (2000), 229–244.
VIseg  J. Van Iseghem, 
Vector orthogonal relations, vector QD-algorithm, 
J. Comput. Appl. Math. 19 (1987), 141–150.
Widder D. V. Widder, 
The Airy Transform, 
Amer. Math. Monthly 86 (1979),  271–277.









§ APPENDIX




As discussed in Lemma <ref>, the polynomial sequences  {P_n^[j](x)}_n⩾ 0, with j=0,1,2, arising from the cubic decomposition of threefold-symmetric 2-Hahn classical polynomials P_n(x) are 2-orthogonal polynomials. The expressions for the corresponding recurrence coefficients can be obtained from the expressions of the γ-coefficients. So, in the light of Lemma <ref>, we computed the expressions of the recurrence coefficients for each of the cubic components of the polynomial sequences discussed in Cases B_1 and C, and this is detailed below. We have included them here for a matter of completion, as these polynomial sequences {P_n^[j](x)}_n⩾ 0 were already described in <cit.> up to a linear change of variable. 





Case B_1





    β_2n^[0] (μ) = 2 (μ +3 n (3 μ +9 n (2 μ +14 n+3)-2)-1)3 (μ +9 n-1) (μ +9 n+2)
       β_2n+1^[0] (μ) =  2 (19 μ +3 n (21 μ +9 n (2 μ +10 n+17)+77)+32)/3 (μ +9 n+2) (μ +9 n+8)
       α_2n^[0] (μ) = 4 n (3 n-1) (6 n-1) (-3 μ ^2-5 μ +702 n^3+27 (8 μ -9) n^2+3 (6 (μ -3) μ -17)
       n+8)/3 (μ +9 n-4) (μ +9 n-1)^2 (μ +9 n+2)
       α_2n+1^[0] (μ) =  4 (2 n+1) (3 n+1) (6 n+1) (2 μ +3 (μ ^2+117 n^3+9 (4 μ +9) n^2+(3 μ  (μ
       +6)+5) n)-5)/3 (μ +9 n-1) (μ +9 n+2)^2 (μ +9 n+5)
       γ_2n^[0] (μ) =  8 n (2 n+1) (3 n-1) (3 n+1) (6 n-1) (6 n+1) (μ +3 n-1) (μ +3 n) (μ +3 n+1)/3 (μ +9
       n-4) (μ +9 n-1)^2 (μ +9 n+2)^2 (μ +9 n+5)
       γ_2n+1^[0] (μ) = 8 (n+1) (2 n+1) (3 n+1) (3 n+2) (6 n+1) (6 n+5)/3 (μ +9 n+2) (μ +9 n+5) (μ +9 n+8)





    β_2n^[1] (μ) = 8 (μ -1)+6 n (9 μ +9 n (2 μ +10 n+7)+5)3 (μ +9 n-1) (μ +9 n+5)
       β_2n+1^[1] (μ) =  2 (31 μ +3 n (27 μ +9 n (2 μ +14 n+31)+202)+143)/3 (μ +9 n+5) (μ +9 n+8)
       α_2n^[1] (μ) =  4 n (36 n^2-1) (-4 μ +3 n (3 (μ -2) μ +117 n^2+36 (μ -1)
       n-10)+4)/3 (μ +9 n-4) (μ +9 n-1)^2 (μ +9 n+2)
       α_2n+1^[1] (μ) = 4 (2 n+1) (3 n+1) (3 n+2)/3 (μ +9 n+2) (μ +9 n+5)^2 (μ +9 n+8)
    ((μ +2) (9 μ +37)+702 n^3+27 (8 μ +43) n^2+3 (6 μ 
       (μ +13)+187) n)
       
       γ_2n^[1] (μ) =   8 n (2 n+1) (3 n+1) (3 n+2) (6 n-1) (6 n+1)/3 (μ +9 n-1) (μ +9 n+2) (μ +9 n+5)
       γ_2n+1^[1] (μ) = 8 (n+1) (2 n+1) (3 n+1) (3 n+2) (6 n+5) (6 n+7) (μ +3 n+1) (μ +3 n+2) (μ +3 n+3)/3
       (μ +9 n+2) (μ +9 n+5)^2 (μ +9 n+8)^2 (μ +9 n+11)



    β_2n^[2] (μ) = 20 (μ +2)+6 n (15 μ +9 n (2 μ +14 n+17)+58)3 (μ +9 n+2) (μ +9 n+5)
       β_2n+1^[2] (μ) = 2 (46 μ +3 n (33 μ +9 n (2 μ +10 n+27)+209)+170)/3 (μ +9 n+5) (μ +9 n+11)
       α_2n^[2] (μ) =  4 n (3 n+1) (6 n+1) (μ +3 (μ ^2+234 n^3+9 (8 μ +17) n^2+(6 μ  (μ +5)+7)
       n)-10)/3 (μ +9 n-1) (μ +9 n+2)^2 (μ +9 n+5)
       α_2n+1^[2] (μ) =  4 (2 n+1) (3 n+2) (6 n+5) /3 (μ +9 n+2) (μ +9 n+5)^2 (μ +9 n+8)
    (2 (μ +2) (3 μ +10)+351 n^3+54 (2 μ +11) n^2+3 (3 μ 
       (μ +14)+98) n)
       
       γ_2n^[2] (μ) = 8 n (2 n+1) (3 n+1) (3 n+2) (6 n+1) (6 n+5) (μ +3 n) (μ +3 n+1) (μ +3 n+2)/3 (μ +9
       n-1) (μ +9 n+2)^2 (μ +9 n+5)^2 (μ +9 n+8)
       γ_2n+1^[2] (μ) =  8 (n+1) (2 n+1) (3 n+2) (3 n+4) (6 n+5) (6 n+7)/3 (μ +9 n+5) (μ +9 n+8) (μ +9 n+11)





Case C




    β _2n^[0](μ,ρ)= 2 (2 n+1) (3 n+1) (6 n+1)/(9 n+μ +2) (9
       n+ρ +3)-4 n (3 n-1) (6 n-1)/(9 n+μ -1) (9 n+ρ -3)
       β _2n+1^[0](μ,ρ)= 4 (n+1) (3 n+2) (6 n+5)/(9 n+μ +8) (9
       n+ρ +6)-2 (2 n+1) (3 n+1) (6 n+1)/(9 n+μ +2) (9 n+ρ +3)
       α _2n^[0](μ,ρ)=  6 n (6 n-2) (6 n-1) /(9 n+μ -4) (9 n+μ -1)^2 (9 n+ρ -3)^2 (9 n+ρ )((6 n-1) (3 n+μ -1)
       (3 n+ρ -1) 
     
       	+(6 n-3) (9 n+μ -1) (3 n+ρ -2) (3 n+ρ -1)/9 n+ρ -6
    	+(6   n+1) (3 n+μ -1) (3 n+μ ) (9 n+ρ -3)/9 n+μ +2) 
    
       α _2n+1^[0](μ,ρ)
     	= 6 (2 n+1) (6 n+1) /(9 n+μ
       +2)^2 (9 n+μ +5) (9 n+ρ ) (9 n+ρ +3)^2(2 (3 n+μ +1) (3
       n+ρ ) (3 n+1)^2 
    
       +6 n (3 n+μ ) (3 n+μ +1) (9 n+ρ +3) (3 n+1)/9 n+μ
       -1
       +(3 n+2) (6 n+2) (9 n+μ +2) (3 n+ρ ) (3 n+ρ +1)/9 n+ρ +6) 
    
       γ _2n^[0](μ,ρ)= 6 n (6 n-2) (6 n-1) (6 n+1) (6 n+2) (6 n+3)
       (μ +3 n-1) (μ +3 n) (μ +3 n+1)/(μ +9 n-4) (μ +9 n-1)^2 (μ +9 n+2)^2 (μ +9 n+5)
       (9 n+ρ -3) (9 n+ρ ) (9 n+ρ +3)
       γ _2n+1^[0](μ,ρ)= 6 n (6 n-2) (6 n-1) (6 n+1) (6 n+2) (6 n+3)
       (μ +3 n-1) (μ +3 n) (μ +3 n+1)/(μ +9 n-4) (μ +9 n-1)^2 (μ +9 n+2)^2 (μ +9 n+5)
       (9 n+ρ -3) (9 n+ρ ) (9 n+ρ +3)




    β _2n^[1](μ,ρ)= 4 (2 n+1) (3 n+1) (3 n+2)/(9 n+μ +5) (9
       n+ρ +3)+2 n-72 n^3/(9 n+μ -1) (9 n+ρ )
       β _2n+1^[1](μ,ρ)= 2 (n+1) (6 n+5) (6 n+7)/(9 n+μ +8) (9
       n+ρ +9)-4 (2 n+1) (3 n+1) (3 n+2)/(9 n+μ +5) (9 n+ρ +3)
       α _2n^[1](μ,ρ)= 6 n (6 n-1) (6 n+1)/(9 n+μ -1)^2 (9 n+μ +2) (9 n+ρ -3)
       (9 n+ρ )^2(6 n (3 n+μ ) (3
       n+ρ -1)
       
     
       +(6 n+2) (9 n+μ -1) (3 n+ρ ) (3 n+ρ -1)/9 n+ρ +3+(6 n-2) (3
       n+μ -1) (3 n+μ ) (9 n+ρ )/9 n+μ -4) 
    
       α _2n+1^[1] (μ,ρ)
     	=6 (2 n+1) /(9 n+μ +2) (9 n+μ +5)^2 (9 n+ρ +3)^2 (9 n+ρ +6)
    ×(6 (2 n+1) (3 n+1) (3 n+2) (3 n+μ +1) (3 n+ρ +1)
       
     
       +(3 n+2) (6 n+1) (6 n+2) (9 n+μ +5) (3 n+ρ ) (3 n+ρ +1)/9
       n+ρ
    
       +(3 n+1) (6 n+4) (6 n+5) (3 n+μ +1) (3 n+μ +2) (9 n+ρ +3)/9 n+μ
       +8)
    
       γ _2n^[1](μ,ρ)= 6 n (6 n-1) (6 n+1) (6 n+2) (6 n+3) (6 n+4) (3
       n+ρ -1) (3 n+ρ ) (3 n+ρ +1)/(μ +9 n-1) (μ +9 n+2) (μ +9 n+5) (9 n+ρ -3) (9
       n+ρ )^2 (9 n+ρ +3)^2 (9 n+ρ +6)
       γ _2n+1^[1](μ,ρ)= 6 n (6 n-1) (6 n+1) (6 n+2) (6 n+3) (6 n+4)
       (3 n+ρ -1) (3 n+ρ ) (3 n+ρ +1)/(μ +9 n-1) (μ +9 n+2) (μ +9 n+5) (9 n+ρ -3)
       (9 n+ρ )^2 (9 n+ρ +3)^2 (9 n+ρ +6)



    β _2n^[2](μ,ρ)= 2 (2 n+1) (3 n+2) (6 n+5)/(9 n+μ +5) (9
       n+ρ +6)-4 n (3 n+1) (6 n+1)/(9 n+μ +2) (9 n+ρ )
       β _2n+1^[2](μ,ρ)= 4 (n+1) (3 n+4) (6 n+7)/(9 n+μ +11) (9
       n+ρ +9)-2 (2 n+1) (3 n+2) (6 n+5)/(9 n+μ +5) (9 n+ρ +6)
       α _2n^[2](μ,ρ)= 6 n (6 n+1) (6 n+2)/(9 n+μ -1) (9 n+μ +2)^2 (9 n+ρ )^2
       (9 n+ρ +3)
    ((6 n+1) (3 n+μ )
       (3 n+ρ )
    
       +(6 n-1) (9 n+μ +2) (3 n+ρ -1) (3 n+ρ )/9 n+ρ -3
    
       +(6 n+3) (3
       n+μ ) (3 n+μ +1) (9 n+ρ )/9 n+μ +5)
    
       α _2n+1^[2](μ,ρ)= 6 (2 n+1) (6 n+5)/(9 n+μ
       +5)^2 (9 n+μ +8) (9 n+ρ +3) (9 n+ρ +6)^2
    (2 (3 n+μ +2) (3
       n+ρ +1) (3 n+2)^2
       
    
       +6 (n+1) (9 n+μ +5) (3 n+ρ +1) (3 n+ρ +2) (3 n+2)/9 n+ρ
       +9
    
       +(3 n+1) (6 n+4) (3 n+μ +1) (3 n+μ +2) (9 n+ρ +6)/9 n+μ +2)
    
       γ _2n^[2](μ,ρ)= 6 n (6 n+1) (6 n+2) (6 n+3) (6 n+4) (6 n+5)
       (μ +3 n) (μ +3 n+1) (μ +3 n+2)/(μ +9 n-1) (μ +9 n+2)^2 (μ +9 n+5)^2 (μ +9 n+8)
       (9 n+ρ ) (9 n+ρ +3) (9 n+ρ +6)
       γ _2n+1^[2](μ,ρ)=6 n (6 n+1) (6 n+2) (6 n+3) (6 n+4) (6 n+5)
       (μ +3 n) (μ +3 n+1) (μ +3 n+2)/(μ +9 n-1) (μ +9 n+2)^2 (μ +9 n+5)^2 (μ +9 n+8)
       (9 n+ρ ) (9 n+ρ +3) (9 n+ρ +6)




