**1. Introduction**

Generalized Laguerre polynomials (GLP) are defined explicitly

$$L\_n^a(\mathbf{x}) = \sum\_{r=0}^n \frac{(-1)^r \ (1+a)\_n \ x^r}{r! \ (n-r)! \ (1+a)\_r \ '} \tag{1}$$

where *a* is a real -valued parameter, (*a*)*r* is the Pochhammer symbol

$$(a)\_r = \begin{cases} a(a+1)\dots(a+(r-1)), & r \ge 1, \\ 1, & r=0. \end{cases}$$

In confluent hypergeometric notation, we have

$$L\_n^a(\mathbf{x}) = \frac{(1+a)\_n}{n!} \,\_1F\_1(-n; a+1; \mathbf{x}). \tag{2}$$

These polynomials satisfy the second-order linear differential equation (see, for example, [1] p. 298)

$$\text{tr}\,\mathbf{D}^2 L\_n^a(\mathbf{x}) + (1+a-\mathbf{x})\,\mathbf{D}L\_n^a(\mathbf{x}) + nL\_n^a(\mathbf{x}) = 0, \quad \mathbf{D} = \frac{d}{d\mathbf{x}}.\tag{3}$$

The so-called Shively's pseudo-Laguerre polynomials *Rn*(*<sup>a</sup>*, *x*) are defined by (see, [2])

$$R\_{ll}(a, \mathbf{x}) = \frac{(a)\_{2u}}{(a)\_n} \,\_1F\_1(-n; a+n; \mathbf{x}),\tag{4}$$

1

which are related to the proper simple Laguerre polynomial (see, [2])

$$L\_n(\mathbf{x}) = \,\_1F\_1(-n; 1; \mathbf{x}),\tag{5}$$

$$R\_n(a, \mathbf{x}) = \frac{1}{(a-1)\_n} \sum\_{r=0}^n \frac{(a-1)\_{n+r} \ L\_{n-r}(\mathbf{x})}{r!}. \tag{6}$$

Shivley deduced the generating function for pseudo Laguerre polynomials of one variable as (see, [2])

$$\left(\mathbf{c}^{2t}\right)\_0 F\_1\left(-; \frac{a}{2} + \frac{1}{2}; t^2 - \mathbf{x}t\right) = \sum\_{n=0}^{\infty} \frac{R\_n(a, \mathbf{x})}{(\frac{a}{2} + \frac{1}{2})\_n} \mathbf{x} \tag{7}$$

Now, owing to the significance of the earlier mentioned work related to Laguerre polynomials, we find record that many authors became interested to study the scalar cases of the classical sets of Laguerre polynomials into Laguerre matrix polynomials. Of those authors, we mention [3–7].

Recently, the matrix versions of the classical families orthogonal polynomials such as Jacobi, Hermite, Chebyshev, Legendre, Gegenbauer, Bessel and Humbert polynomials of one variables and some other polynomials were introduced by many authors for matrices in C *N*×*N* and various properties satisfied by them were given from the scalar case. Rather than giving an exhaustive list of references, we refer the reader to the article [8]. Theory of generalized and multivariable orthogonal matrix polynomials has provided new means of analysis to deal with the majority of problems in mathematical physics which find broad practical applications. In [9,10], Subuhi Khan and others introduced the 2-variable forms of Laguerre and modified Laguerre matrix polynomials and generalized Hermite matrix based polynomials of two variables and Lie algebraic techniques. Furthermore, several papers concerning the orthogonal matrix polynomials for two and multivariables have become more and more relevant, see for example [11–17].

The section-wise treatment is as follows. In Section 2, we deals with some basic facts, notations and results to that are needed in the work. In Section 3, we define Shivleys matrix polynomials of two variables and to study their properties. The generating matrix functions, matrix recurrence relations, summation formula and operational representations these new matrix polynomials are obtained. Some special cases of the established results are also underlined as corollaries. Finally, we give some concluding remarks in Section 4.

Throughout this paper, for C *N* denote the *N*-dimensional complex vector space and C *N*×*N* denote all square matrices with *N* rows and *N* columns with entries are complex numbers, **Re**(*z*) and **Im**(*z*) denote the real and imaginary parts of a complex number *z*, respectively. For any matrix A in C *N*×*N*, *σ*(*A*) is the spectrum of *A*, the set of all eigenvalues of *A*, which will be denoted by *A*, is defined by

$$||A|| = \sup\_{x \neq 0} \frac{||Ax||\_2}{||x||\_2},$$

where for a vector *y* in C *N*, ||*y*||2 = (*yHy*) 1 2 is Euclidean norm of *y*. I and **0** stand for the identity matrix and the null matrix in C *N*×*N*, respectively.
