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

TL;DR

This paper introduces a comprehensive method for constructing symmetric tensor polynomials in D-dimensional space that are orthonormal under various weights, generalizing classical polynomials and enabling new applications in physics and mathematics.

Contribution

It presents a unified framework for generating orthonormal tensor polynomials in multiple dimensions, including new polynomials for non-standard weights like Fermi-Dirac and Bose-Einstein distributions.

Findings

Derived explicit formulas for polynomials up to fifth order.

Validated the orthonormalization process through matching equations and coefficients.

Extended classical polynomials to new weights relevant in physics.

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.

Equations

$$\int d^D\mathbf{\xi }\omega (\mathbf{\xi }){\mathcal{P}}_{i_1\cdots i_N}(\mathbf{\xi }){\mathcal{P}}_{j_1\cdots j_M}(\mathbf{\xi })={\delta }_{NM}{\delta }_{i_1\cdots i_N|j_1\cdots j_M}.$$

$${\mathcal{P}}_{i_1\cdots i_N}(-{\xi }_{i_1},\dots ,-{\xi }_{i_k},\dots ,-{\xi }_{i_N})=(-1{)}^N{\mathcal{P}}_{i_1\cdots i_N}({\xi }_{i_1},\dots ,{\xi }_{i_k},\dots {\xi }_{i_N})$$

$$\displaystyle\delta_{i_{1}\cdots i_{N}|j_{1}\cdots j_{N}}\equiv\delta_{i_{1}j_{1}}\cdots\delta_{i_{N}j_{N}}\,+\,\mbox{permutations of $i$'s},$$

$${\delta }_{i_1\cdots i_N{j}_1\cdots j_N}\equiv {\delta }_{i_1{j}_1}\cdots {\delta }_{i_N{j}_N}+\text{mbox}allpermutations.$$

$${\mathcal{P}}_0(\mathbf{\xi })=c_0,$$

$${\mathcal{P}}_{i_1}(\mathbf{\xi })=c_1{\xi }_{i_1},$$

$${\mathcal{P}}_{i_1{i}_2}(\mathbf{\xi })=c_2{\xi }_{i_1}{\xi }_{i_2}+f_2(\xi ){\delta }_{i_1{i}_2},\text{mbox}where{f}_2(\xi )\equiv {\bar{c}}_2{\xi }^2+{c^{\prime }}_2,$$

$$\displaystyle\mathcal{P}_{i_{1}i_{2}i_{3}}(\bm{\xi})=c_{3}\,\xi_{i_{1}}\xi_{i_{2}}\xi_{i_{3}}+f_{3}(\xi)\big{(}\xi_{i_{1}}\delta_{i_{2}i_{3}}+\xi_{i_{2}}\delta_{i_{1}i_{3}}+\xi_{i_{3}}\delta_{i_{1}i_{2}}\big{)},\;\mbox{where}\;f_{3}(\xi)\equiv{\bar{c}}_{3}\xi^{2}+{c^{\prime}}_{3},$$

$$\displaystyle\mathcal{P}_{i_{1}i_{2}i_{3}i_{4}}(\bm{\xi})=c_{4}\,\xi_{i_{1}}\xi_{i_{2}}\xi_{i_{3}}\xi_{i_{4}}+f_{4}(\xi)\,\big{(}\xi_{i_{1}}\xi_{i_{2}}\delta_{i_{3}i_{4}}+\xi_{i_{1}}\xi_{i_{3}}\delta_{i_{2}i_{4}}+\xi_{i_{1}}\xi_{i_{4}}\delta_{i_{2}i_{3}}+\xi_{i_{2}}\xi_{i_{3}}\delta_{i_{1}i_{4}}+\xi_{i_{2}}\xi_{i_{4}}\delta_{i_{1}i_{3}}+\xi_{i_{3}}\xi_{i_{4}}\delta_{i_{1}i_{2}}\big{)}$$

$$\displaystyle+g_{4}(\xi)\,\delta_{i_{1}i_{2}i_{3}i_{4}},\;\mbox{where}\;f_{4}(\xi)\equiv\big{(}{\bar{c}}_{4}\xi^{2}+{c^{\prime}}_{4}\big{)},\>\mbox{and}\;g_{4}(\xi)\equiv\big{(}{\bar{d}}_{4}\xi^{4}+{d^{\prime}}_{4}\xi^{2}+{d}_{4}\big{)}.$$

$$I_N{\delta }_{i_1\cdots i_N}\equiv \int d^D\mathbf{\xi }\omega (\mathbf{\xi }){\xi }_{i_1}\cdots {\xi }_{i_N}$$

$$\displaystyle I_{2N}=\frac{\pi^{\frac{D}{2}}}{2^{N-1}\Gamma\big{(}N+\frac{D}{2}\big{)}}\int_{0}^{\xi_{max}}d\xi\,\omega(\xi)\,\xi^{2N+D-1}.$$

$${\delta }_{i_1\cdots i_N{j}_1\cdots j_N}={\delta }_{i_1\cdots i_N|j_1\cdots j_N}+\text{mbox}othertensors$$

$${\delta }_{ii\cdots i\cdots jj\cdots j}={\delta }_{ii\cdots i|jj\cdots j}+\text{mbox}othertensors.$$

$${\delta }_{ij}={\delta }_{i|j},$$

$${\delta }_{ii|jj}={\delta }_{ij}{\delta }_{ij}$$

$${\delta }_{iijj}={\delta }_{ij}{\delta }_{ij}+{\delta }_{ii}{\delta }_{jj}.$$

$${\delta }_{iiijjj}={\delta }_{ij}{\delta }_{ij}{\delta }_{ij}+{\delta }_{ii}{\delta }_{jj}{\delta }_{ij}.$$

$${\delta }_{iiiijjjj}={\delta }_{ij}{\delta }_{ij}{\delta }_{ij}{\delta }_{ij}+{\delta }_{ii}{\delta }_{jj}{\delta }_{ij}{\delta }_{ij}+{\delta }_{ii}{\delta }_{ii}{\delta }_{jj}{\delta }_{jj}$$

$$\left({\mathcal{P}}_{(M)},{\mathcal{P}}_{(N)}\right)\equiv \int d^D\mathbf{\xi }\omega (\mathbf{\xi }){\mathcal{P}}_{(M)}(\mathbf{\xi }){\mathcal{P}}_{(N)}(\mathbf{\xi },$$

$$c_0^2\int d^D\mathbf{\xi }\omega (\xi )=1,$$

$$c_0=\pm \frac{1}{\sqrt{I_0}},$$

$$c_1^2\int d^D\mathbf{\xi }\omega (\xi ){\xi }_{i_1}{\xi }_{j_1}={\delta }_{i_1{j}_1}$$

$$c_1=\pm \frac{1}{\sqrt{I_2}}.$$

$$c_0\int d^D\mathbf{\xi }\omega (\xi )\left[c_2{\xi }_{i_1}{\xi }_{i_2}+f_2(\xi ){\delta }_{i_1{i}_2}\right]=0.$$

$$I_2{c}_2+D{I}_2{\bar{c}}_2+I_0{c}_2^{\prime }=0.$$

$$\int d^D\mathbf{\xi }\omega (\xi )\left[c_2{\xi }_{i_1}{\xi }_{i_2}+f_2(\xi ){\delta }_{i_1{i}_2}\right]\left[c_2{\xi }_{j_1}{\xi }_{j_2}+f_2(\xi ){\delta }_{j_1{j}_2}\right]={\delta }_{i_1{j}_1}{\delta }_{i_2{j}_2}+{\delta }_{i_1{j}_2}{\delta }_{i_2{j}_1}.$$

$$\displaystyle c_{2}^{2}I_{4}\big{(}\delta_{i_{1}j_{1}}\delta_{i_{2}j_{2}}+\delta_{i_{1}j_{2}}\delta_{i_{2}j_{1}}+\delta_{i_{1}i_{2}}\delta_{j_{1}j_{2}}\big{)}+\big{[}2(D+2)I_{4}c_{2}\bar{c}_{2}+2I_{2}c_{2}c_{2}^{\prime}+\bar{c}_{2}^{2}D(D+2)I_{4}+2DI_{2}\bar{c}_{2}c_{2}^{\prime}+c_{2}^{\prime 2}I_{0}\big{]}\delta_{i_{1}i_{2}}\delta_{j_{1}j_{2}}$$

$$={\delta }_{i_1{j}_1}{\delta }_{i_2{j}_2}+{\delta }_{i_1{j}_2}{\delta }_{i_2{j}_1}.$$

$$c_2^2{I}_4=1\text{mbox}and,c_2^2{I}_4+2(D+2)I_4{c}_2{\bar{c}}_2+2I_2{c}_2{c}_2^{\prime }+{\bar{c}}_2^{2}D(D+2)I_4+2D{I}_2{\bar{c}}_2{c}_2^{\prime }+c_2^{\prime 2}{I}_0=0.$$

$$c_2=\pm \frac{1}{\sqrt{I_4}}.$$

$$c_2^{\prime }=-\frac{I_2}{I_0}{c}_2-D\frac{I_2}{I_0}{\bar{c}}_2.$$

Full text

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 $D\rightarrow 1$. 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 $D\rightarrow 1$ 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 ($I_{N}$). 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 $\omega(\bm{\xi})$ where the vector $\bm{\xi}\equiv(\xi_{1},\xi_{2},\cdots,\xi_{D})$ is defined. We claim here the existence of a set of orthonormal polynomials $\mathcal{P}_{i_{1}\cdots i_{N}}(\bm{\xi})$ in this space.

[TABLE]

The polynomials $\mathcal{P}_{i_{1}\cdots i_{N}}(\bm{\xi})$ are symmetric tensors in the D-dimensional Euclidean space, expressed in terms of the vector components $\xi_{i}$ and of $\delta_{ij}$. The N${}^{\mbox{th}}$ order polynomial is symmetrical in the indices $i_{1}\cdots i_{N}$, and its parity is $(-1)^{N}$.

[TABLE]

The following tensors, defined by Harold Grad Grad (1949b), are expressed as sums of products of the Kronecker’s delta function, ($\delta_{ij}=1$ for $i=j$ and 0 for $i\neq j$).

[TABLE]

and,

[TABLE]

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,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and,

[TABLE]

Therefore the Nth order polynomial is the sum of all possible symmetric tensors built from products of $\xi_{i}$ and of $\delta_{ij}$ times coefficients which are themselves polynomials in $\xi^{2}$ to maximum allowed power. This proposal yields a unique expression for the Nth order polynomial.

We define integrals $I_{N}$ which are central to the present study. They are assumed to exist and have well defined properties.

[TABLE]

Hereafter the weight function is assumed to only depend on the modulus of the vector: $\omega(\bm{\xi})=\omega(\xi)$, $\xi\equiv|\bm{\xi}|$. By symmetry it holds that $I_{2N+1}=0$ since the integral vanishes. Using the spherical integration volume, $\int d^{D}\bm{\xi}\,\omega(\bm{\xi})=D\pi^{D/2}/\Gamma(D/2+1)\int d\xi\,\xi^{D-1}\omega(\xi)$, the $I_{2N}$ integrals become,

[TABLE]

In case that $\xi_{max}=\infty$ the weight function must have the property $\omega(\bm{\xi})\rightarrow 0$ for $\xi\rightarrow\infty$ faster than any power of $\xi$. 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, $\delta_{i_{1}\cdots i_{N}|j_{1}\cdots j_{N}}$ and $\delta_{i_{1}\cdots i_{N}j_{1}\cdots j_{N}}$, 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 $I_{N}$ 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 $\delta_{i_{1}\cdots i_{N}|j_{1}\cdots j_{N}}$ has N! terms since this tensor is a sum over all possible permutations of the $i$’s under a fixed set of $j$’s. The tensor $\delta_{i_{1}\cdots i_{N}j_{1}\cdots j_{N}}$ has (2N-1)!/2N-1(N-1)!=(2N-1)(2N-3)(2N-5)…1 terms according to the arguments below. Firstly notice that the tensor $\delta_{i_{1}\cdots i_{N}j_{1}\cdots j_{N}}$ has more terms than $\delta_{i_{1}\cdots i_{N}|j_{1}\cdots j_{N}}$ and here we seek to find these remaining tensors.

[TABLE]

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 $i$ from $j$ indices. In this notation the previous expression becomes equal to,

[TABLE]

For N=1 we have that,

[TABLE]

and we express this identity with respect to the number of terms simply as $1=1$. For N=2 notice that, $\delta_{i_{1}i_{2}j_{1}j_{2}}=\delta_{i_{1}j_{1}}\delta_{i_{2}j_{2}}+\delta_{i_{1}j_{2}}\delta_{i_{2}j_{1}}+\delta_{i_{1}i_{2}}\delta_{j_{1}j_{2}}$. Since $\delta_{i_{1}i_{2}|j_{1}j_{2}}=\delta_{i_{1}j_{1}}\delta_{i_{2}j_{2}}+\delta_{i_{1}j_{2}}\delta_{i_{2}j_{1}}$ the above relation becomes, $\delta_{i_{1}i_{2}j_{1}j_{2}}=\delta_{i_{1}i_{2}|j_{1}j_{2}}+\delta_{i_{1}i_{2}}\delta_{j_{1}j_{2}}$. Therefore it holds that $\delta_{i_{1}i_{2}j_{1}j_{2}}=\delta_{i_{1}i_{2}|j_{1}j_{2}}+\delta_{i_{1}i_{2}}\delta_{j_{1}j_{2}}$. While $\delta_{i_{1}i_{2}j_{1}j_{2}}$ contains 3 components, $\delta_{i_{1}i_{2}|j_{1}j_{2}}$ has only 2, such that it holds for this decomposition that $3=2\oplus 1$. In the short notation the above relation becomes,

[TABLE]

The tensors $\delta_{ij}\delta_{ij}$ and $\delta_{ii}\delta_{jj}$ contain 2 and 1 components, respectively. Thus the short notation gives that $\delta_{ii\cdots i|jj\cdots j}\equiv\delta_{ij}\delta_{ij}\cdots\delta_{ij}$ where the products of $\delta_{ij}$ takes into account all possible permutations. The tensor $\delta_{ii\cdots ijj\cdots j}$ cannot be expressed similarly because not all combinations of $\delta_{ij}$’s are included. Therefore the need to decompose it into $\delta_{ii\cdots i|jj\cdots j}$ plus other tensors. For N=3 according to the short tensorial notation,

[TABLE]

To determine the number of components of this tensor, notice that for $\delta_{i_{1}i_{2}i_{3}j_{1}j_{2}j_{3}}$ once a pair is chosen, say $\delta_{i_{1}j_{1}}$, 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 $\delta_{ij}\delta_{ij}\delta_{ij}$ has 3! components, thus to know the number of terms in the tensor $\delta_{ii}\delta_{ij}\delta_{jj}$ we use the following argument. There are 3 components in $\delta_{ii}$, namely, $\delta_{i_{1}i_{2}}$, $\delta_{i_{1}i_{3}}$ and $\delta_{i_{2}i_{3}}$ and similarly, 3 components in $\delta_{jj}$. Once fixed $\delta_{ii}$, and $\delta_{jj}$ the tensor $\delta_{ij}$ has only one possible component left. Therefore $\delta_{ii}\delta_{ij}\delta_{jj}$ has a total of 3 times 1 times 3, that is 9 components, and the tensorial decomposition is expressed as $15=6\oplus 9$.

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

[TABLE]

The same reasoning of the previous cases is used here, namely, once a pair is fixed, say $\delta_{i_{1}j_{1}}$, 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 $\delta_{ij}\delta_{ij}\delta_{ij}\delta_{ij}$ has 4! terms, The tensor $\delta_{ii}\delta_{jj}\delta_{ij}\delta_{ij}$ has 6 possible terms for $\delta_{ii}$, $\delta_{i_{1}i_{2}}$, $\delta_{i_{1}i_{3}}$, $\delta_{i_{1}i_{4}}$, $\delta_{i_{2}i_{3}}$, $\delta_{i_{2}i_{4}}$ and $\delta_{i_{3}i_{4}}$, and the same applies for $\delta_{jj}$. Thus only 2 choices are left for $\delta_{ij}\delta_{ij}$, once the indices of $\delta_{ii}$ and $\delta_{jj}$ are fixed. Hence the total number of terms is 6 times 6 times 2, namely, 72 terms. For the tensor $\delta_{ii}\delta_{ii}\delta_{jj}\delta_{jj}$ the pairs $\delta_{ii}\delta_{ii}$ have only three terms each and so is for $\delta_{jj}\delta_{jj}$, such that the total is 3 times 3, namely, 9 terms. The tensorial decomposition is expressed as $105=24\oplus 72\oplus 9$.

IV Orthonormalization of polynomials to order N=4

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

[TABLE]

used here for $M,N=0,\cdots 4$. The Nth order polynomial is shortly referred as $\mathcal{P}_{(N)}$ such that the polynomials of Eqs.(5), (6), (7), (8), (9) are called $\mathcal{P}_{(0)}$, $\mathcal{P}_{(1)}$, $\mathcal{P}_{(2)}$, $\mathcal{P}_{(3)}$ and $\mathcal{P}_{(4)}$, respectively. The inner product between polynomials with distinct parity vanishes, $\left(\mathcal{P}_{(even)},\mathcal{P}_{(odd)}\right)=0$. 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 ($c_{K}$ for $K=0,1,2,3,4$, ${c^{\prime}}_{K}$ for $K=2,3,4$, ${\bar{c}}_{K}$ for $K=2,3,4$, $d_{4}$, ${d^{\prime}}_{4}$ and ${\bar{d}}_{4}$). Remarkably 14 equations arise from the 9 orthonormalization conditions $\left(\mathcal{P}_{(M)},\mathcal{P}_{(N)}\right)$, as seen below.

$\bullet\quad\left(\mathcal{P}_{(0)},\mathcal{P}_{(0)}\right)$

The normalization of the N=0 polynomial is,

[TABLE]

which gives that,

[TABLE]

where Eq. (5) and the definition of $I_{0}$ in Eq.(10) have been used.

$\bullet\quad\left(\mathcal{P}_{(1)},\mathcal{P}_{(1)}\right)$

[TABLE]

The definition of $I_{2}$ in Eq.(10) is invoked to obtain that,

[TABLE]

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 $c_{0}$ and $c_{1}$. The polynomial $\mathcal{P}_{(2)}$ 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.

$\bullet\quad\left(\mathcal{P}_{(0)},\mathcal{P}_{(2)}\right)$

This conditions is equivalent to $c_{0}\left(1,\mathcal{P}_{(2)}\right)=0$ since $\mathcal{P}_{(0)}$ is a constant.

[TABLE]

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

[TABLE]

$\bullet\quad\left(\mathcal{P}_{(2)},\mathcal{P}_{(2)}\right)$

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

[TABLE]

This is because while the orthonormalization is associated to the tensor $\delta_{i_{1}j_{1}}\delta_{i_{2}j_{2}}+\delta_{i_{1}j_{2}}\delta_{i_{2}j_{1}}$ the integration over $\xi_{i_{1}}\xi_{i_{2}}\xi_{j_{1}}\xi_{j_{2}}$ leads to the tensor $\delta_{i_{1}j_{1}}\delta_{i_{2}j_{2}}+\delta_{i_{1}j_{2}}\delta_{i_{2}j_{1}}+\delta_{i_{1}i_{2}}\delta_{j_{1}j_{2}}$. This difference is responsible for the onset of more than one condition. Using the tensorial formulas of appendix A, one obtains that:

[TABLE]

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

[TABLE]

The three equations are promptly solved and the coefficients $c_{2}$, ${\bar{c}}_{2}$ and ${c^{\prime}}_{2}$ are determined below.

[TABLE]

Using Eq.(23) to write $c_{2}^{\prime}$, we have

[TABLE]

Substituting $c_{1}^{\prime}$ in Eq.(27),

[TABLE]

which is a equation for $\bar{c}_{2}$. The solutions, are:

[TABLE]

The coefficient $\bar{c}_{2}$ must be a real number, and this is the case provided that $(D+2)-DJ_{2}\geq 0$. From Eq.(23), we calculate $c_{2}^{\prime}$, to obtain that,

[TABLE]

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.

$\bullet\quad\left(\mathcal{P}_{(1)},\mathcal{P}_{(3)}\right)$

[TABLE]

The integrals are calculated with the help of the tensorial formulas of appendix A, such that the above expression becomes $\left[c_{3}I_{4}+\bar{c}_{3}I_{4}(D+2)+c^{\prime}_{3}I_{2}\right]\delta_{ijjj}=0$. The tensor $\delta_{ijjj}$ is short for $\delta_{i_{1}j_{1}j_{2}j_{3}}$, according to the notation of section III. Thus we obtain that,

[TABLE]

$\bullet\quad\left(\mathcal{P}_{(3)},\mathcal{P}_{(3)}\right)$

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

[TABLE]

We stress the difference between the tensors $\delta_{iiijjj}$ and $\delta_{iii|jjj}$ 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, $\int d^{D}\bm{\xi}\omega(\xi)\;\xi_{i}\xi_{i}\xi_{i}\xi_{j}\xi_{j}\xi_{j}=I_{6}\delta_{iiijjj}$, where all permutations of $i$’s and $j$’s are taken into account. From the other side in this short notation, $\delta_{iii|jjj}=\delta_{ij}\delta_{ij}\delta_{ij}$, as explained in section III, and, as shown there, $\delta_{iiijjj}=\delta_{ij}\delta_{ij}\delta_{ij}+\delta_{ii}\delta_{ij}\delta_{jj}$. Thus the above integrals are computed with the aid of appendix A. Using this notation, Eq.(35) becomes,

[TABLE]

which gives two equations:

[TABLE]

The solution of the first equation is,

[TABLE]

Eq.(34) is used to eliminate $c^{\prime}_{3}$ from Eq.(37),

[TABLE]

This equation can be solved for $\bar{c}_{3}$.

[TABLE]

and $c^{\prime}_{3}$ is calculated by Eq.(34):

[TABLE]

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

$\bullet\quad\left(\mathcal{P}_{(0)},\mathcal{P}_{(4)}\right)$

The orthonormalization with the N=0 polynomial means that $c_{0}\left(1,\mathcal{P}_{(4)}\right)$=0.

[TABLE]

A single equation results from this integral since it can only be proportional to the tensor $\delta_{iiii}$.

[TABLE]

$\bullet\quad\left(\mathcal{P}_{(2)},\mathcal{P}_{(4)}\right)$

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

[TABLE]

The sixth order tensorial integral, $\int d^{D}\bm{\xi}\omega(\xi)\xi_{i}\xi_{i}\xi_{j}\xi_{j}\xi_{j}\xi_{j}=I_{6}\delta_{iijjjj}$ has 15 terms. This tensor can be decomposed as $\delta_{iijjjj}=\delta_{ii}\delta_{jjjj}+\delta_{ij}\delta_{ijjj}$, namely, as a sum of two other tensors, which have 3 and 12 terms, respectively. This decomposition can be formally expressed as $15=3\oplus 12$, as discussed in section III. We notice the presence of the tensor $\delta_{i_{1}i_{2}j_{1}j_{2}}\delta_{j_{3}j_{4}}+\delta_{i_{1}i_{2}j_{1}j_{3}}\delta_{j_{2}j_{4}}+\delta_{i_{1}i_{2}j_{1}j_{4}}\delta_{j_{2}j_{3}}+\delta_{i_{1}i_{2}j_{2}j_{3}}\delta_{j_{1}j_{4}}+\delta_{i_{1}i_{2}j_{2}j_{4}}\delta_{j_{1}j_{3}}+\delta_{i_{1}i_{2}j_{3}j_{4}}\delta_{j_{1}j_{2}}$, which has 18 terms, and is equal to $2\delta_{jjjj}\delta_{ii}+\delta_{iiij}\delta_{ij}$ (2 times 3 plus 12). Therefore Eq.(IV) becomes,

[TABLE]

This lead to two equations, one proportional to $\delta_{ij}\delta_{ijjj}$,

[TABLE]

and the other proportional to $\delta_{ii}\delta_{jjjj}$,

[TABLE]

respectively.

$\bullet\quad\left(\mathcal{P}_{(4)},\mathcal{P}_{(4)}\right)$

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

[TABLE]

Here we extensively use the short notation of section III, where $\delta_{iiii|jjjj}=\delta_{ij}\delta_{ij}\delta_{ij}\delta_{ij}$. There is the integral $\int d^{D}\bm{\xi}\omega(\xi)\;\xi_{i}\xi_{i}\xi_{i}\xi_{i}\xi_{j}\xi_{j}\xi_{j}\xi_{j}=I_{8}\delta_{iiiijjjj}$. It holds that $\delta_{iiii|jjjj}=\delta_{ij}\delta_{ij}\delta_{ij}\delta_{ij}+\delta_{ii}\delta_{jj}\delta_{ij}\delta_{ij}+\delta_{ii}\delta_{ii}\delta_{jj}\delta_{jj}$. There are 105 terms in $\delta_{iiiijjjj}$, 24 in $\delta_{ij}\delta_{ij}\delta_{ij}\delta_{ij}$, 72 in $\delta_{ii}\delta_{jj}\delta_{ij}\delta_{ij}$, and 9 in $\delta_{ii}\delta_{ii}\delta_{jj}\delta_{jj}$. We computed each of the integrals individually and identify them by the following notation that uses the product of their coefficients $c_{4}$, $f_{4}$ and $g_{4}$. For instance the first one is

[TABLE]

The second term is,

[TABLE]

where we have used that $\delta_{iiiijj}\delta_{jj}=\delta_{ii}\delta_{jj}\delta_{ij}\delta_{ij}+2\delta_{ii}\delta_{ii}\delta_{jj}\delta_{jj}$.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Next these integrals are introduced into Eq.(47) to obtain that,

[TABLE]

The three remaining equations are the coefficients of the three independent tensors in the above equation.

[TABLE]

[TABLE]

[TABLE]

The six equations given by the Eqs.(IV), (45), (IV), (54), (55) and (IV) can be solved to obtain the six coefficients. Nevertheless notice that coefficients $c_{4}$, $\bar{c}_{4}$ and $c^{\prime}_{4}$ are determined from a sub set of equations, namely, Eqs.(45), (54), and (55). From them it follows that,

[TABLE]

Hence one obtains that,

[TABLE]

The remaining coefficients are given by,

[TABLE]

[TABLE]

and,

[TABLE]

where

[TABLE]

[TABLE]

Notice that $\delta_{2K}=2I_{2K+2}I_{2K-2}/\Delta_{2K}^{2}$.

Hence we have determined all the coefficients in Eqs.(20), (22), (28), (29), (31), (38), (40), (41), (57), (59), (60), (61), (62), and (63). Some of the coefficients can be summarized in simple formulas for all polynomials, using a general notation:

[TABLE]

Notice that we have chosen the positive solutions for the square roots, but the negative ones would also lead to orthogonal polynomials. We have completed the orthonormalization of the first five polynomials (N=0,1,2,3,4). The procedure can be applied for higher polynomials, although it becomes increasingly laborious.

V D-dimensional polynomials for specific weights

In this section we obtain the D-dimensional polynomials for some specific weights and from them retrieve some well known D=1 polynomials. We also obtain the $I_{N}$ functions associated to some new weights.

V.1 D-dimensional Hermite polynomials

We show that the D-dimensional Hermite polynomials are straightforwardly retrieved from the above polynomials for the gaussian weight,

[TABLE]

To obtain the integrals $I_{2N}$ of Eq.(69), we note that,

[TABLE]

Then it follows from Eq.(69) that,

[TABLE]

In this limit $c_{K}=1$, $\bar{c}_{K}=0$ $c^{\prime}_{K}=-1$, $d_{4}=1$, $\bar{d}_{4}=0$ and $d^{\prime}_{4}=0$, and the polynomials of Eqs.(5), (6), (7), (8), (9) become,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and,

[TABLE]

We notice that the tensorial basis that spans the new generalized polynomials contains the basis that spans the Hermite polynomials but not vice-versa. The D-dimensional Hermite polynomials $\mathcal{P}_{i_{1}\cdots i_{N}}$ are symmetric tensors in the indices $i_{1}\cdots i_{N}$ spanned over the basis formed by the tensors,

[TABLE]

This basis is not large enough to span the new generalized polynomials, $\mathcal{P}_{i_{1}\cdots i_{N}}$, which demand a larger basis formed by the tensors

[TABLE]

whose scalar functions F\big{(}\xi\big{)} are polynomials in powers of the vector modulus, $1$, $\xi^{2}$, $\xi^{4}$, …,$\xi^{2k}$.

V.2 D-dimensional Legendre polynomials

We define D-dimensional Legendre polynomials satisfying the general orthonormal relation of Eq.(1) with $\xi_{max}=1$ and the following weight:

[TABLE]

From this it follows that,

[TABLE]

Using the expressions corresponding coefficients up to fourth order

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

These coefficients provides orthogonal polynomials for any dimension and, as we will see in subsection V.5, they give the well known D=1 Legendre polynomials for the particular case with $D=1$.

V.3 D-dimensional Chebyshev polynomials of first kind

Following the same procedure used for the Legendre polynomials, we define the D-dimensional Chebyshev polynomials of first kind with the weight $\omega(\xi)$ below and integrated under a sphere of radius $\xi_{max}$,

[TABLE]

with the corresponding integral

[TABLE]

and coefficients up to fourth order

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The polynomials with the coefficients above are orthogonal for any dimensions and, as we will see in subsection V.5 for $D=1$ they give the well known D=1 Chebyshev polynomials of first kind.

V.4 D-dimensional Chebyshev polynomials of second kind

For the D-dimensional Chebyshev polynomials of second kind, we define the following weight function and maximum radius

[TABLE]

which give the integral and coefficients below

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The corresponding polynomials are orthogonal for any dimension and, as we will see in subsection V.5, give the well known D=1 Chebyshev polynomials of second kind.

V.5 Projection of Hermite, Legendre and Chebyshev polynomials into D=1 dimension

To obtain the projection of such polynomials in D=1 dimension, it suffices to drop the index of the vector, $\xi_{i_{1}}\rightarrow\xi$, and to take that $\delta_{i_{1}i_{2}}\rightarrow 1$ since there is only one index, and so, $i_{1}=i_{2}=1$. Thus the tensors based on the Kronecker’s delta function have each of its terms equal to one, for instance, $\delta_{i_{1}j_{1}}\cdots\delta_{i_{N}j_{N}}=1$, and so, $\delta_{i_{1}\cdots i_{N}|j_{1}\cdots j_{N}}=N\,\!!$ and $\delta_{i_{1}\cdots i_{N}\,j_{1}\cdots j_{N}}=(2N-1)\,\!!/[2^{N-1}(N-1)\,\!!]$, according to section III. Hence the orthonormality condition of Eq.(1) becomes,

[TABLE]

Notice that this is not the standard normalization employed in the definition of most of D=1 orthonormal polynomials.

V.5.1 Hermite polynomials

The Hermite polynomials obtained via dimensional reduction from Eqs.(71), (72), (73), (74) and (75) are given by,

[TABLE]

The known D=1 Hermite polynomials, defined with the normalization of probability theory Gradshteyn and Ryzhik (2014), satisfy the following orthonormality condition,

[TABLE]

The first five ones are ${{\mathit{He}}_{0}(x)=1}$, ${{\mathit{He}}_{1}(x)=x}$, ${{\mathit{He}}_{2}(x)=x^{2}-1}$, ${{\mathit{He}}_{3}(x)=x^{3}-3x}$, and ${{\mathit{He}}_{4}(x)=x^{4}-6x^{2}+3}$ and they coincide exactly with the present polynomials.

V.5.2 **Legendre polynomials **

It follows by taking $D=1$ that,

[TABLE]

[TABLE]

The Legendre polynomials Gradshteyn and Ryzhik (2014) satisfy the following orthonormality condition,

[TABLE]

The first five ones are ${{\mathit{P}}_{0}(x)=1}$, ${{\mathit{P}}_{1}(x)=x}$, ${{\mathit{P}}_{2}(x)=(3x^{2}-1)/2}$, ${{\mathit{P}}_{3}(x)=(5x^{3}-3x)/2}$, and ${{\mathit{P}}_{4}(x)=(35x^{4}-30x^{2}+3)/8}$. Thus apart from the normalization they coincide with the present polynomials.

V.5.3 Chebyshev polynomials of the first kind

[TABLE]

[TABLE]

The Chebyshev polynomials of the first kind Gradshteyn and Ryzhik (2014) satisfy the following orthonormality condition,

[TABLE]

The first five ones are ${{\mathit{T}}_{0}(x)=1}$, ${{\mathit{T}}_{1}(x)=x}$, ${{\mathit{T}}_{2}(x)=2x^{2}-1}$, ${{\mathit{T}}_{3}(x)=4x^{3}-3x}$, and ${{\mathit{T}}_{4}(x)=8x^{4}-8x^{2}+1}$. Thus apart from the normalization they coincide with the present polynomials.

V.5.4 Chebyshev polynomials of the second kind

[TABLE]

[TABLE]

The Chebyshev polynomials of the second kind Gradshteyn and Ryzhik (2014) satisfy the following orthonormality condition,

[TABLE]

The first five ones are ${{\mathit{U}}_{0}(x)=1}$, ${{\mathit{U}}_{1}(x)=2x}$, ${{\mathit{U}}_{2}(x)=4x^{2}-1}$, ${{\mathit{U}}_{3}(x)=8x^{3}-4x}$, and ${{\mathit{U}}_{4}(x)=16x^{4}-12x^{2}+1}$. Thus apart from the normalization they coincide with the present polynomials.

V.6 D-dimensional Fermi-Dirac polynomials

We seek the set of polynomials orthonormal under the weight defined by the Fermi-Dirac statistical occupation number for $\bm{u}=0$.

[TABLE]

The parameter $z\equiv e^{\mu/\theta}$ is called the fugacity, where $\mu$ is the chemical potential and $\theta$ is the temperature using the so-called reduced units $(m=k_{B}=c=\hbar=1)$. For the Fermi-Dirac weight it holds that,

[TABLE]

It is interesting to consider special limits where the integral $g_{\nu}(z)$ can be explicitly calculated and the $I_{2N}$ obtained. One of such limits is when the quantum Fermi-Dirac statistics becomes the classical Maxwell-Boltzmann statistics. This is the small fugacity limit,

[TABLE]

The other interesting limit is the so-called Sommerfeld limit Pathria and Beale (2011), used for the treatment of electrons in metals where the chemical potential (Fermi energy) is much larger than the room temperature. In this limit $\mu/\theta=\ln{z}>>1$, such that the terms of high order in $(\theta/\mu)^{2k}$ can be disregarded.

[TABLE]

V.7 D-dimensional Bose-Einstein polynomials

Similarly to the previous case, we seek the set of polynomials orthonormal under a weight which is the Bose-Einstein statistical occupation number.

[TABLE]

It holds that,

[TABLE]

The are also special limits here, the first being when the quantum Bose-Einstein statistics becomes the classical Maxwell-Botlzmann one. This is the small fugacity limit,

[TABLE]

Notice that in leading order in $z$, the Fermi-Dirac and the Bose-Einstein polynomials become identical as they reduce to the Maxwell-Boltzmann polynomials at finite temperature. In fact at this limit these polynomials are just scaled versions of the D-dimensional Hermite polynomials. The other interesting limit is that of negative vanishing fugacity near to the onset of the Bose-Einstein condensate. For this we write $\alpha\equiv-\mu/\theta$, $\alpha\rightarrow 0$. Thus for $z\equiv e^{-\alpha}$, one obtains that,

[TABLE]

where $\zeta(s)$ is the Riemann zeta function, which is defined for ${\mathrm{Re}}(s)>1$. Thus for $\nu$ integer, the above expression must be replaced by,

[TABLE]

V.8 D-dimensional Graphene polynomials

Graphene is a two-dimensional sheet of carbon atoms arranged in an hexagonal lattice where electrons move with a relativistic dispersion relation. To deal with this situation we seek a set of polynomials orthonormal under a weight which the following Fermi-Dirac statistical occupation number for $\bm{u}=0$Oettinger et al. (2013).

[TABLE]

It holds that,

[TABLE]

The limit of low doping in graphene corresponds to $\mu\rightarrow 0$ or, equivalently, $z\rightarrow 1$. The doping of graphene can be chemically adjusted for instance liu11.

[TABLE]

V.9 D-dimensional Yukawa polynomials

Consider the Yukawa potential of nuclear interactions that contains the parameter $\mu$ that renders it short ranged. Assume the Yukawa potential as a weight to obtain orthonormal polynomials.

[TABLE]

To determine the coefficients that enter the polynomials, it is enough to have the integrals,

[TABLE]

VI Expansion of a function in general D-dimensional polynomials

Suppose a function $f(\bm{\xi})$ whose argument is the D-dimensional vector $\bm{\xi}$ such that a small vector $\bm{u}$ is subtracted from it. We want to expand it in terms of this small correction using the D-dimensional polynomial basis.

[TABLE]

The coefficients $\mathcal{A}_{i_{1}\,i_{2}\cdots i_{N}}$ and the polynomials $\mathcal{P}_{i_{1}\,i_{2}\cdots i_{N}}$ only depend on $\bm{u}$ and $\bm{\xi}$, respectively.

To help the understanding of the above expansion, we consider the example of a local distribution of charges, $\rho(\bm{u})$, that interact through a two-body potential, $f(\bm{u}-\bm{u^{\prime}})$, such as the Yukawa potential. We want to obtain the potential $v(\bm{\xi})$ far way from this charge distribution, namely at a distance ${\bm{\xi}}\gg{\bm{u}}$. As the potential produced by each charge at volume $d^{D}\bm{u}\,\rho(\bm{u})$ contributes with potential $f(\bm{\xi}-\bm{u})$ to the total potential, one obtains that,

[TABLE]

Therefore we propose here a novel expansion in terms of the orthonormal multipoles $\mathcal{Q}_{i_{1}\,i_{2}\cdots i_{N}}$, called in this way because of the orthonormal basis,

[TABLE]

The key ingredient to obtain these orthonormal multipoles is that the potential $f(\bm{\xi})$ itself enters in the definition of the polynomials $\mathcal{P}_{i_{1}\,i_{2}\cdots i_{N}}(\bm{\xi})$ as the weight function that renders them orthonormal. To understand this we turn to the general properties of the expansion of Eq.(125). Notice that it is of the form $f({\bm{\xi}}-{\bm{u}})=f({\bm{\xi}})\cdot S({\bm{\xi}},{\bm{u}})$, which means that the expansion $S({\bm{\xi}},{\bm{u}}\rightarrow 0)\rightarrow 1$. Thus the corrections in ${\bm{u}}$ are necessarily small.

The key ingredient that renders the above expansion advantageous is the assumption that the weight function that defines the set of orthonormal polynomials is the expanded function itself for $\bm{u}=0$,

[TABLE]

This means that the function has to fulfill the properties previously required for a weight in order to span a polynomial basis. In this case the coefficients are readily determined using the orthonormality condition of Eq.(1).

[TABLE]

Hence the above expression once fed back to Eq.(125) gives the completeness relation for the polynomials.

[TABLE]

The completness relation is given by

[TABLE]

since $f(\bm{\xi}-\bm{u})=\int d^{D}\bm{\xi^{\prime}}\,\delta^{D}(\bm{\xi^{\prime}}-\bm{\xi})f(\bm{\xi^{\prime}}-\bm{u})$.

In order to determine other properties of the coefficients, we define $\bm{\eta}\equiv\bm{\xi}-\bm{u}$ to obtain that,

[TABLE]

At this point we introduce N=0 polynomial, which is constant, into this expression.

[TABLE]

Hence the determination of the coefficients is reduced to the expansion $\mathcal{P}_{i_{1}\,i_{2}\cdots i_{N}}(\bm{\eta}+\bm{u})$ as a sum over polynomials of equal or lower order ($M\leq N$), $\mathcal{P}_{i_{1}\,i_{2}\cdots i_{M}}(\bm{\eta})$.

[TABLE]

where $\mathcal{U}_{i_{1}i_{2}\ldots i_{M}}(\bm{u})$ is a polynomial in $\bm{u}$ that multiplies the polynomial in $\bm{\eta}$ of order $N-M$. The sought coefficient is nothing but the term proportional to $\mathcal{P}_{0}(\bm{\eta})$ in this expansion. The limit $\bm{u}\rightarrow 0$ in the above expression shows that $\mathcal{U}_{0}(0)=1$ while for all other higher order tensors it holds that $\mathcal{U}_{i_{1}i_{2}\ldots i_{M}}(0)=0$. In summary the coefficients become,

[TABLE]

and $\mathcal{A}_{i_{1}i_{2}\ldots i_{N}}(\bm{u}=0)=0$ for $N\geq 1$.

Next we directly obtain the coefficients by directly expanding $\mathcal{P}_{i_{1}i_{2}\ldots i_{N}}(\bm{\eta}+\bm{u})$ up to the order N=4. The expansion of the N=0 polynomial is trivial since $\mathcal{P}_{0}(\bm{\eta}+\bm{u})=c_{0}$

[TABLE]

Expanding the N=1 polynomial, $\mathcal{P}_{i_{1}}(\bm{\eta}+\bm{u})=c_{1}(\eta_{i_{1}}+u_{i_{1}})$, so the coefficient is,

[TABLE]

Expanding the N=2 polynomial, $\mathcal{P}_{i_{1}i_{2}}(\bm{\eta}+\bm{u})=c_{2}(\eta_{i_{1}}+u_{i_{1}})(\eta_{i_{2}}+u_{i_{2}})+[\bar{c}_{2}(\bm{\eta}+\bm{u})^{2}+c^{\prime}_{2}]\delta_{i_{1}i_{2}}=\mathcal{P}_{i_{1}i_{2}}(\bm{\eta})+(c_{2}/c_{1})[u_{i_{1}}\mathcal{P}_{i_{2}}(\bm{\eta})+u_{i_{2}}\mathcal{P}_{i_{1}}(\bm{\eta})]+(2\bar{c}_{2}/c_{1})u_{i_{3}}\mathcal{P}_{i_{3}}(\bm{\eta})\delta_{i_{1}i_{2}}+\mathcal{P}_{0}(\bm{\eta})\left[c_{2}u_{i_{1}}u_{i_{2}}+\bar{c}_{2}\bm{u}^{2}\delta_{i_{1}i_{2}}\right]/c_{0}$, its N=0 coefficient gives that,

[TABLE]

using that $c_{0}^{2}=1/I_{0}$. Similarly the expansion of the N=3 polynomial, $\mathcal{P}_{i_{1}i_{2}i_{3}}(\bm{\eta}+\bm{u})$, contains the N=0 polynomial plus higher order ones that will be omitted for simplicity.

[TABLE]

Using the equations for the coefficients of polynomials, we have

[TABLE]

Therefore the third order coefficient of the expansion is

[TABLE]

Expanding the N=4 polynomial $\mathcal{P}_{i_{1}i_{2}i_{3}i_{4}}(\bm{\eta}+\bm{u})$ in terms of $\bm{\eta}$ is a laborious task as terms in $\bm{\eta}$ must be expressed again as functions of polynomials of order N. Nevertheless we only seek the N=0 term and some considerations can be applied to simplify the task. For instance, the odd terms ($\eta_{i_{1}}$, $\eta_{i_{1}}\bm{\eta}^{2}$, $\eta_{i_{1}}\eta_{i_{2}}\eta_{i_{3}}$) do not contribute to the calculation of $\mathcal{A}_{i_{1}i_{2}i_{3}i_{4}}$ and one can take that $\bm{\eta}^{2}=DI_{2}/I_{0}+\mathcal{O}(\bm{\eta})$. After some algebra, we have:

[TABLE]

Finally, the N=4 coefficient is

[TABLE]

We summarize the coefficients below, obtained after some additional algebraic manipulation. Notice that they are functions of the integrals $I_{2N}$ previously defined.

[TABLE]

The contraction between tensors $\mathcal{A}_{i_{1}i_{2}\ldots i_{N}}$ and $\mathcal{P}_{i_{1}i_{2}\ldots i_{N}}$ up to N=4 order are obtained below.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

VII Conclusions

We propose here D-dimensional symmetric tensor polynomials orthonormal under a general weight. We show that the number of coefficients of the first five ones (N=0 to 4) matches exactly the number of equations that stem from their orthonormalization, which allows for their obtainment as functions of the integrals $I_{N}$’s. In Statistical Mechanics it is well-known that the D-dimensional Hermite polynomials are a key element to solve the Boltzmann equation for classical particles, which satisfy the Maxwell-Boltzmann statistics Kremer (2010); Philippi et al. (2006); Coelho and Neumann (2016). The present generalized polynomials are applicable to semi-classical fluids where the particles obey the Bose-Einstein and Fermi-Dirac statistics Coelho et al. (2014, 2016). The proposed generalized polynomials allows for the definition of orthonormal multipoles as shown in case of the Yukawa potential. We foresee many other applications because the proposed polynomials take into account the expanded function as the weight that render them orthonormal in the D-dimensional space.

Acknowledgements.

R. C. V. Coelho thanks to FAPERJ and to the European Research Council (ERC) Advanced Grant 319968-FlowCCS for the financial support and to Hans J. Herrmann for the kind hospitality at ETH Zürich.

Appendix A Tensorial Identities involving the $I_{N}$

Here we obtain especial formulas derived from the definition of the $I_{N}$’s given in Eq.(10) and obtained from tensorial contractions of the tensor $\xi_{i_{1}}\cdots\xi_{i_{N}}$, and, consequently, from the $\delta_{i_{1}\cdots i_{N}}$.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]