In electrodynamics I have seen the following:
Let $\phi$ be a solution to the Poisson equation $-\Delta \phi= \rho$, and assume that $\rho$ is compactly supported. Then we can expand $\phi$ as the following: $$ \phi(x) = \frac {q}{\|x\|} + \frac{\langle x,p\rangle}{\|x\|^3} + \frac{1}{2}\sum_{i,j=1}^n Q_{i,j}\frac{x_ix_j}{\|x\|^5}+\dots,$$ where $q=\int_\mathbb{R^n} \rho(x) dx$, $p = \int_\mathbb{R^n} x\rho(x) dx$, $Q_{i,j}=\int_\mathbb{R^n} (3x_ix_j - \|x\|^2\delta_{i,j})\rho(x) dx$.
I am interested in making this mathematically precise.
Is this an expansion of $\phi$ in terms of some orthonormal basis of $L^2(\mathbb{R^n})$? How does this basis look like? Can any solution to the Poisson equation be expanded in that way? What about harmonic functions?
It's not very hard to derive the first few terms of the expression, and one can even go for cumbersome yet exact expression for the multipole moment tensors!
We just need to use the Green's function solution for the potential
$$\phi(x)=\int d^3x'\frac{\rho(x')}{||x-x'||}$$
and then use the multipole expansion in terms of Legendre polynomials
$$\frac{1}{||x-x'||}=\sum_{n=0}^{\infty}P_n(\hat{x}\cdot\hat{x}')\frac{||x'||^n}{||x||^{n+1}}$$
and rewrite the equation in a suggestive way:
$$\phi(x)=\sum_{n=0}^{\infty}\frac{1}{||x||^{2n+1}}\int d^3x'(||x||||x'||)^nP_n(\hat{x}\cdot \hat{x}')\rho(x')$$
Let's consider the first few terms of this expansion with their associated Legendre polynomials that can be found here.
Here, it may not be so clear how to bring this object to the final form. We can rewrite however the form in brackets equivalently as
$$\int d^3x'\rho(x')\Big[3(\sum_ix_ix'_i)^2-||x||^2||x'||^2\Big]=\sum_{ij}{x_ix_j}\int d^3x'\rho(x')\Big[3(x'_ix'_j)-\delta_{ij}||x'||^2\Big]\equiv\sum_{ij}Q_{ij}x_ix_j$$
Again we rewrite the term in brackets as follows
$$\int d^3x'\rho(x')\Big[5(\sum_ix_ix'_i)^3-3||x||^2||x'||^2(\sum_ix_ix'_i)\Big]=\sum_{ijk} x_ix_jx_k\int d^3x'\rho(x')(5x_i'x_j'x_k'-3x_i||x'||^2\delta_{jk})$$
and we define the octupole moment to be
$$O_{ijk}=\int d^3x \rho(x)( 5x_ix_jx_k-3x_i||x||^2\delta_{jk})$$
we could go on and on with explicit expressions but we can write a formula for $2^n$-pole moment. Suppose the coefficients of Legendre polynomials are given by $P_n(x)=\sum_{m=0}^na_{nm}x^m$. Then we can easily show that with the intuition obtained in the examples above
$$\int d^3x'(||x||||x'||)^nP_n(\hat{x}\cdot \hat{x}')\rho(x')=\sum_{k_1,...k_n}M^{(n)}_{k_1,...k_n}x_{k_1}...x_{k_n}$$
where
$$\begin{align}M^{(n)}_{k_1,...k_n}&=\int d^3 x' (a_{n1}x_{k_1}'...x_{k_n}'+a_{n3}x_{k_1}'...x_{k_{n-2}}'||x'||^2\delta_{k_{n-1}k_n}+...\\&=\sum_{m=0}^n\int d^3 x' \rho(x')\Big[a_{nm}x_{k_1}'...x_{k_m}'||x'||^{n-m}\delta_{k_{m+1}k_{m+2} }...\delta_{k_{n-1}k_n}\Big]\\&=\sum_{m=0}^n\int d^3 x' \rho(x')\Big[a_{nm}||x'||^{n-m}\prod_{i=1}^m x'_{k_i}\prod_{i=0}^{\frac{n-m-1}{2} }\delta_{k_{m+2i+1}k_{m+2i+2} }\Big]\end{align}$$