Chebyshev, Legendre, Hermite and other orthonormal polynomials in D-dimensions

# Chebyshev, Legendre, Hermite and other orthonormal polynomials in D-dimensions

Mauro M. Doria Departamento de Física dos Sólidos, Universidade Federal do Rio de Janeiro, 21941-972 Rio de Janeiro, Brazil    Rodrigo C. V. Coelho ETH Zürich, Computational Physics for Engineering Materials, Institute for Building Materials Schafmattstrasse 6, HIF, CH-8093 Zürich, Switzerland Departamento de Física dos Sólidos, Universidade Federal do Rio de Janeiro, 21941-972 Rio de Janeiro, Brazil
July 29, 2019
###### Abstract

We propose a general method to construct symmetric tensor polynomials in the D-dimensional Euclidean space which are orthonormal under a general weight. The D-dimensional Hermite polynomials are a particular case of the present ones for the case of a gaussian weight. Hence we obtain generalizations of the Legendre and of the Chebyshev polynomials in D dimensions that reduce to the respective well-known orthonormal polynomials in D=1 dimensions. We also obtain new D-dimensional polynomials orthonormal under other weights, such as the Fermi-Dirac, Bose-Einstein, Graphene equilibrium distribution functions and the Yukawa potential. We calculate the series expansion of an arbitrary function in terms of the new polynomials up to the fourth order and define orthonormal multipoles. The explicit orthonormalization of the polynomials up to the fifth order (N from 0 to 4) reveals an increasing number of orthonormalization equations that matches exactly the number of polynomial coefficients indication the correctness of the present procedure.

###### pacs:
74.20.De, 74.25.Dw,74.25.Ha

## I Introduction

The theory of orthonormal polynomials is still an unfolding branch of Mathematics and Physics Szegö (1975); Gradshteyn and Ryzhik (2014); Siafarikas (2001); Chihara (2001). It started in the nineteen century and provided the key ingredients for the following century development of Quantum Mechanics. Interestingly it was Charles Hermite who firstly introduced tensorial properties to the D-dimensional orthonormal polynomials. For the case of a gaussian weight he obtained symmetric tensorial polynomials and showed them to reduce to the well-known previously found one-dimensional Hermite polynomials by taking the limit . Since then the D-dimensional Hermite polynomials have been studied under several aspects, such as the obtainment of recurrence formulas Berkowitz and Garner (1970). Nevertheless the tensorial orthonormal polynomials have not been systematically studied for weights other than the gaussian one so far, apart from a few attempts, such as the Laguerre weight Wünsche (2015) in D=2 with the intent to apply in quantum optics Wünsche (2001); Kok and Braunstein (2001).

In this paper we propose a general method to construct D-dimensional tensorial polynomials orthonormal under an arbitrary weight. We apply this method for the first five polynomials (N=0 to 4) and determine their coefficients. The D-dimensional Hermite polynomials are retrieved as the particular case of the gaussian weight. We obtain D-dimensional generalizations of the Legendre and Chebyshev of first and second kind polynomials. By taking the limit the one-dimensional Legendre and Chebyshev of first and second kind polynomials are retrieved. We choose other weights to construct new D-dimensional tensorial polynomials, such as the Bose-Einstein, the Fermi-Dirac and also the Graphene equilibrium distribution functions. Their interest is in the search for solutions of the Boltzmann equation describing semi-classical fluids Coelho et al. (2014, 2016). In such cases the corresponding D-dimensional Euclidean space is that of the microscopic velocity. As a last example we construct the D-dimensional polynomials for the Yukawa weight, which are useful in position space to derive the concept of an orthonormal multipole series expansion.

The D-dimensional Hermite polynomials have many interesting applications in Physics ranging from Quantum Optics Kok and Braunstein (2001) to Statistical Mechanics. In the latter case they offer fundamental aid to solve the Boltzmann equation Grad (1949a); Kremer (2010); Philippi et al. (2006); Coelho and Neumann (2016) for classical particles. Indeed it was H. Grad who first used the D-dimensional Hermite polynomials to describe the microscopic velocity space of the Boltzmann equation Grad (1949b, a); Kremer (2010). The Boltzmann equation aims a statistical description of an ensemble of particles and so describes the motion of a set of particles at a scale between the microscopic and the macroscopic levels. While the microscopic level has a deterministic description of motion, since Newton’s law is applied to the individual particles, at the macroscopic level the only laws available are those of conservation of mass and momentum for many particles. For fluids and gases the macroscopic level corresponds to the continuity and to the Navier-Stokes equation, respectively, and it can be shown that both follow from the Boltzmann equation Kremer (2010); Philippi et al. (2006); Coelho et al. (2014). A few decades ago the study of the Boltzmann equation experimented a revival because of a new method developed to solve it on a lattice version of position space. Because of its simplicity this method revolutionized the way to numerically tackle problems in fluid dynamics. It became known as the Lattice Boltzmann Method Krüger et al. (2016); Succi (2001) and uses the D-dimensional Hermite polynomials to span the distribution function, which essentially gives the number of particles in a point in phase space. The Gauss-Hermite quadrature is also used in this method to perform integration in the D-dimensional space.

Recently it was found that to render the Lattice Boltzmann applicable to semi-classical fluids the weight that render the D-dimensional polynomials must be the equilibrium distribution function itself. Therefore the gaussian weight is not appropriate because it is associated to the Maxwell-Boltzmann equilibrium distribution function whereas for semi-classical fluids the particles obey the Bose-Einstein or the Fermi-Dirac equilibrium distribution functionsCoelho et al. (2016). Therefore D-dimensional polynomials orthonormal under these weights are need in such cases and so, one must go beyond the D-dimensional Hermite polynomials.

The remarkably rich tensorial structure of the D-dimensional space is the key element that allows for the existence of the present symmetric tensor polynomials orthonormal under a general weight. This rich tensorial structure was first observed by Harold Grad Grad (1949b), but never developed to obtain D-dimensional orthonormal polynomials for a general weight. Here we develop this proposal and obtain the D-dimensional orthonormal polynomials. Notice that this rich tensorial structure is not present for D=1 since the Kronecker’s delta function is trivial and equal to one. However in higher dimensions many tensors can be built as products and sums of the Kronecker’s delta function. Harold Grad was the first to notice this wealth of tensors built from Kronecker’s delta function Grad (1949b).

This paper is organized as follows. In section II we propose the general form of the D-dimensional tensorial polynomials and explicitly write the first five ones. In section III the rich tensorial properties of D-dimensional space. The explicit construction of the first five (N=0 to 4) orthonormal polynomials is carried in section IV, which means that all their coefficients are obtained as functions of some integrals over the weight (). Next we apply this general theory to specific weights in section V. The known D-dimensional Hermite polynomials are derived from the present ones and also new D-dimensional generalizations of the Legendre and Chebyshev polynomials of the first and second kinds are proposed here. The projection of such polynomials to D=1 dimensions does give the well-known Hermite in subsection V.1, Legendre (subsection V.2) and Chebyshev (subsections V.3 and  V.4) which are projected to D=1 dimension (subsection V.5). Next we consider D-dimensional polynomials orthonormal under new weights, such as Fermi-Dirac (subsection V.6), Bose-Einstein (subsection V.7), graphene (subsection V.8) and Yukawa potential (subsection V.9). Finally we expand a general function in terms of the D-dimensional polynomials in section VI, which leads to the proposal of orthonormal multipoles. We reach conclusions in section VII. Some useful tensorial identities are discussed in appendix A.

## Ii General D-dimensional Polynomials

Consider the D-dimensional Euclidean space endowed with a weight function where the vector is defined. We claim here the existence of a set of orthonormal polynomials in this space.

 ∫dDξω(ξ)Pi1⋯iN(ξ)Pj1⋯jM(ξ)=δNMδi1⋯iN|j1⋯jM. (1)

The polynomials are symmetric tensors in the D-dimensional Euclidean space, expressed in terms of the vector components and of . The N order polynomial is symmetrical in the indices , and its parity is .

 Pi1⋯iN(−ξi1,…,−ξik,…,−ξiN)=(−1)NPi1⋯iN(ξi1,…,ξik,…ξiN) (2)

The following tensors, defined by Harold Grad Grad (1949b), are expressed as sums of products of the Kronecker’s delta function, ( for and 0 for ).

 δi1⋯iN|j1⋯jN≡δi1j1⋯δiNjN+permutations of i's, (3)

and,

 δi1⋯iNj1⋯jN≡δi1j1⋯δiNjN+ all permutations. (4)

The knowledge of the number of terms in such tensors is useful and discussed in more details in section III. The first five (N=0 to 4) polynomials are given by,

 P0(ξ)=c0, (5)
 Pi1(ξ)=c1ξi1, (6)
 Pi1i2(ξ)=c2ξi1ξi2+f2(ξ)δi1i2,wheref2(ξ)≡¯c2ξ2+c′2, (7)
 Pi1i2i3(ξ)=c3ξi1ξi2ξi3+f3(ξ)(ξi1δi2i3+ξi2δi1i3+ξi3δi1i2),wheref3(ξ)≡¯c3ξ2+c′3, (8)

and,

 Pi1i2i3i4(ξ)=c4ξi1ξi2ξi3ξi4+f4(ξ)(ξi1ξi2δi3i4+ξi1ξi3δi2i4+ξi1ξi4δi2i3+ξi2ξi3δi1i4+ξi2ξi4δi1i3+ξi3ξi4δi1i2) (9)

Therefore the N order polynomial is the sum of all possible symmetric tensors built from products of and of times coefficients which are themselves polynomials in to maximum allowed power. This proposal yields a unique expression for the N order polynomial.

We define integrals which are central to the present study. They are assumed to exist and have well defined properties.

 INδi1⋯iN≡∫dDξω(ξ)ξi1⋯ξiN (10)

Hereafter the weight function is assumed to only depend on the modulus of the vector: , . By symmetry it holds that since the integral vanishes. Using the spherical integration volume, , the integrals become,

 I2N=πD22N−1Γ(N+D2)∫ξmax0dξω(ξ)ξ2N+D−1. (11)

In case that the weight function must have the property for faster than any power of . Next we shall explicitly prove the orthonormality of the first five polynomials.

## Iii D-dimensional tensors based on the Kronecker’s delta

The orthonormality condition of Eq.(1) shows a rich tensorial structure in D dimensions revealed by the following two important tensors, and , defined in Eqs.(3) and (4), respectively. The former is associated with the orthonormality condition while the latter is the totally symmetric tensor introduced in the definition of the functions given by Eq.(10). Both tensors are expressed as sums over several terms each one expressed as a product of Kronecker’s delta functions. We determine the number of terms in these two tensors. The tensor has N! terms since this tensor is a sum over all possible permutations of the ’s under a fixed set of ’s. The tensor has (2N-1)!/2(N-1)!=(2N-1)(2N-3)(2N-5)…1 terms according to the arguments below. Firstly notice that the tensor has more terms than and here we seek to find these remaining tensors.

 δi1⋯iNj1⋯jN=δi1⋯iN|j1⋯jN+ other tensors (12)

Next we determine the remaining “other tensors” in case N=0 to 4 and determine the number of components of the above tensor by induction. For this we introduce a short notation that only distinguishes indices from indices. In this notation the previous expression becomes equal to,

 δii⋯i⋯jj⋯j=δii⋯i|jj⋯j+ other tensors. (13)

For N=1 we have that,

 δij=δi|j, (14)

and we express this identity with respect to the number of terms simply as . For N=2 notice that, . Since the above relation becomes, . Therefore it holds that . While contains 3 components, has only 2, such that it holds for this decomposition that . In the short notation the above relation becomes,

 δii|jj=δijδij (15) δiijj=δijδij+δiiδjj. (16)

The tensors and contain 2 and 1 components, respectively. Thus the short notation gives that where the products of takes into account all possible permutations. The tensor cannot be expressed similarly because not all combinations of ’s are included. Therefore the need to decompose it into plus other tensors. For N=3 according to the short tensorial notation,

 δiiijjj=δijδijδij+δiiδjjδij. (17)

To determine the number of components of this tensor, notice that for once a pair is chosen, say , the previous N=2 is retrieved concerning the number of components. Since there are 5 ways to construct this first pair, the total number of components is 5 times 3, that is, 15 terms. The tensor has 3! components, thus to know the number of terms in the tensor we use the following argument. There are 3 components in , namely, , and and similarly, 3 components in . Once fixed , and the tensor has only one possible component left. Therefore has a total of 3 times 1 times 3, that is 9 components, and the tensorial decomposition is expressed as .

Finally for N=4 the and short tensorial notation gives that,

 δiiiijjjj=δijδijδijδij+δiiδjjδijδij+δiiδiiδjjδjj (18)

The same reasoning of the previous cases is used here, namely, once a pair is fixed, say , the number of terms of the remaining indices is provided by the previous N=3 case. There are 7 ways to construct this first pair, thus the total number of terms is 7 times 15, that is, 105 terms. The tensor has 4! terms, The tensor has 6 possible terms for , , , , , and , and the same applies for . Thus only 2 choices are left for , once the indices of and are fixed. Hence the total number of terms is 6 times 6 times 2, namely, 72 terms. For the tensor the pairs have only three terms each and so is for , such that the total is 3 times 3, namely, 9 terms. The tensorial decomposition is expressed as .

## Iv Orthonormalization of polynomials to order N=4

The coefficients of the polynomials are determined here in terms of the integrals . This is done for the first five polynomials by explicitly computing their inner products. We refer to Eq.(1) by the short notation,

 (P(M),P(N))≡∫dDξω(ξ)P(M)(ξ)P(N)(ξ,

used here for . The N order polynomial is shortly referred as such that the polynomials of Eqs.(5), (6), (7), (8), (9) are called , , , and , respectively. The inner product between polynomials with distinct parity vanishes, . Thus the only relevant orthonormality relations are among polynomials with the same parity (odd with odd and even with even).

The number of equations given by orthonormalization conditions must be equal to the number of free coefficients. Indeed this is the case, according to Eqs.(5), (6), (7), (8), (9). The total number of coefficients is 14 ( for , for , for , , and ). Remarkably 14 equations arise from the 9 orthonormalization conditions , as seen below.

The normalization of the N=0 polynomial is,

 c20∫dDξω(ξ)=1, (19)

which gives that,

 c0=±1√I0, (20)

where Eq. (5) and the definition of in Eq.(10) have been used.

 c21∫dDξω(ξ)ξi1ξj1=δi1j1 (21)

The definition of in Eq.(10) is invoked to obtain that,

 c1=±1√I2. (22)

The N=0 and N=1 polynomials are naturally orthogonal because they have distinct parity and to make them orthonormal is enough to normalize them which has been done above by determining the coefficients and . The polynomial has three coefficients, according to Eq.(7), and so, three equations are needed to determine them. These equations must arise from the tensorial structure of the D dimensional space.

This conditions is equivalent to since is a constant.

 c0∫dDξω(ξ)[c2ξi1ξi2+f2(ξ)δi1i2]=0. (23)

Using the tensorial formulas of appendix A and Eq.(10), it follows that,

 I2c2+DI2¯c2+I0c′2=0. (24)

Remarkably the normalization of the N=2 polynomial leads to multiple equations, in this case the following two equations.

 ∫dDξω(ξ)[c2ξi1ξi2+f2(ξ)δi1i2][c2ξj1ξj2+f2(ξ)δj1j2]=δi1j1δi2j2+δi1j2δi2j1. (25)

This is because while the orthonormalization is associated to the tensor the integration over leads to the tensor . This difference is responsible for the onset of more than one condition. Using the tensorial formulas of appendix A, one obtains that:

 =δi1j1δi2j2+δi1j2δi2j1. (26)

The independence of the two tensors leads to two independent equations, given by,

 c22I4=1and,c22I4+2(D+2)I4c2¯c2+2I2c2c′2+¯c22D(D+2)I4+2DI2¯c2c′2+c′22I0=0. (27)

The three equations are promptly solved and the coefficients , and are determined below.

 c2=±1√I4. (28)

Using Eq.(23) to write , we have

 c′2=−I2I0c2−DI2I0¯c2. (29)

Substituting in Eq.(27),

 ¯c22[D(D+2)I4−D2I22I0]+¯c2c2[2(D+2)I4−2DI22I0]+c22(I4−I22I0)=0, (30)

which is a equation for . The solutions, are:

 ¯c2=c2D(−1+Δ2),whereΔ2≡±√2(D+2)−DJ2,J2≡I22I0I4 (31)

The coefficient must be a real number, and this is the case provided that . From Eq.(23), we calculate , to obtain that,

 c′2=−c2I2I0Δ2. (32)

The orthogonalization of the first three polynomials has been concluded here. We proceed to the next order (N=3), and because of the increasing difficulty introduce a short notation for tensors, which is discussed in section III. The N=3 polynomial has three coefficients, similarly to the N=2 case, and so, three equations are needed to determine them.

 (33)

The integrals are calculated with the help of the tensorial formulas of appendix A, such that the above expression becomes . The tensor is short for , according to the notation of section III. Thus we obtain that,

 c3I4+¯c3I4(D+2)+c′3I2=0. (34)

This normalization condition gives the two other equations necessary to calculate the coefficients.

 =δi1i2i3|j1j2j3 (35)

We stress the difference between the tensors and as the reason for multiple equations from a single normalization condition. Using the short notation of section III the integral over the six vector components becomes, , where all permutations of ’s and ’s are taken into account. From the other side in this short notation, , as explained in section III, and, as shown there, . Thus the above integrals are computed with the aid of appendix A. Using this notation, Eq.(35) becomes,

 c23I6(δijδijδij+δiiδijδjj)+[2I6(D+4)c3¯c3+2I4c3c′3+I2c′32+2I4(D+2)¯c3c′3+I6(D+2)(D+4)¯c23]δiiδijδjj =δijδijδij, (36)

which gives two equations:

 c23I6=1and,c23I6+2I6(D+4)c3¯c3+2I4c3c′3+I2c′32+2I4(D+2)¯c3c′3+I6(D+2)(D+4)¯c23=0. (37)

The solution of the first equation is,

 c3=±1√I6. (38)

Eq.(34) is used to eliminate from Eq.(37),

 ¯c23(D+2)[(I6−I24I2)(D+2)+2I6]+2[(I6−I24I2)(D+2)+2I6]c3¯c3+(I6−I24I2)c23=0. (39)

This equation can be solved for .

 ¯c3=c3D+2(−1+Δ4), (40)

and is calculated by Eq.(34):

 c′3=−I4I2Δ4c3% whereΔ4=±√2(D+4)−J4(D+2),J4=I24I6I2. (41)

The N=4 polynomial of Eq.(9) sets a new level of difficulty as six equations must be obtained to calculate the six coefficients.

The orthonormalization with the N=0 polynomial means that =0.

 c0∫dDξω(ξ)[c4ξj1ξj2ξj3ξj4+f4(ξ)(ξj1ξj2δj3j4+ξj1ξj3δj2j4+ξj1ξj4δj2j3+ξj2ξj3δj1j4+ξj2ξj4δj1j3 +ξj3ξj4δj1j2)+g4(ξ)δj1j2j3ij4]=0 (42)

A single equation results from this integral since it can only be proportional to the tensor .

The integration of the N=2 with the N=4 polynomial gives that,

 ∫dDξω(ξ)[c2ξi1ξi2+f2(ξ)δi1i2][c4ξj1ξj2ξj3ξj4+f4(ξ)(ξj1ξj2δj3j4+ξj1ξj3δj2j4+ξj1ξj4δj2j3+ξj2ξj3δj1j4 +ξj2ξj4δj1j3+ξj3ξj4δj1j2)+g4(ξ)δj1j2j3j4]=0. (43)

The sixth order tensorial integral, has 15 terms. This tensor can be decomposed as , namely, as a sum of two other tensors, which have 3 and 12 terms, respectively. This decomposition can be formally expressed as , as discussed in section III. We notice the presence of the tensor , which has 18 terms, and is equal to (2 times 3 plus 12). Therefore Eq.(IV) becomes,

 (44)

This lead to two equations, one proportional to ,

 c4I6+[c′4I4+¯c4I6(D+4)]=0, (45)

and the other proportional to ,

 c4c2I6+c4[c′2I4+¯c2I6(D+4)]+2c2[c′4I4+¯c4I6(D+4)]+2[c′4c′2I2+(¯c4c′2+c′4¯c2)I4(D+2) +¯c4¯c2I6(D+4)(D+2)]+c2[d4I2+d′4I4(D+2)+¯d4I6(D+4)(D+2)]+[d4c′2I0+(d4¯c2+d′4c′2)I2D (46)

respectively.

The normalization of the N=4 polynomial is given by,

 ∫dDξω(ξ)[c4ξi1ξi2ξi3ξi4+f4(ξ)(ξi1ξi2δi3i4+ξi1ξi3δi2i4+ξi1ξi4δi2i3+ξi2ξi3δi1i4+ξi2ξi4δi1i3+ξi3ξi4δi1i2) +g4(ξ)δi1i2i3i4][c4ξj1ξj2ξj3ξj4+f4(ξ)(ξj1ξj2δj3j4+ξj1ξj3δj2j4+ξj1ξj4δj2j3+ξj2ξj3δj1j4+ξj2ξj4δj1j3 +ξj3ξj4δj1j2)+g4(ξ)δj1j2j3j4]=δi1i2i3i4|j1j2j3j4. (47)

Here we extensively use the short notation of section III, where . There is the integral . It holds that . There are 105 terms in , 24 in , 72 in , and 9 in . We computed each of the integrals individually and identify them by the following notation that uses the product of their coefficients , and . For instance the first one is

 ∫dDξω(ξ)c24[⋯]=I8c24(δijδijδijδij+δiiδjjδijδij+δiiδiiδjjδjj)