I am trying to construct a new set of basis functions on the sphere. The reason is I am working in a special coordinate system $(r,\theta_q,\phi_q)$ where $\theta_q = \theta_q(r,\theta,\phi)$ and $\phi_q = \phi_q(r,\theta,\phi)$ are nonlinear functions of the usual spherical coordinates. The coordinates $(r,\theta_q,\phi_q)$ have a geometry more suitable to my problem, but it is important to note that they are not an orthogonal coordinate system (e.g. the basis vectors in the $(r,\theta_q,\phi_q)$ system are not orthogonal).
I want to develop basis functions similar to the spherical harmonics which exhibit the geometry of this new coordinate system. I proceeded by defining $$ \alpha_l^m(\mathbf{r}) = Y_l^m(\theta_q(\mathbf{r}),\phi_q(\mathbf{r})) $$ where $Y_l^m$ are the usual spherical harmonics. The $\alpha_l^m$ exhibit the geometrical properties of my coordinate system, however they are not eigenfunctions of the Laplacian like the spherical harmonics. So in the next step, I define: $$ \alpha_l^m(\mathbf{r}) \approx \beta_l^m(\mathbf{r}) = \sum_{k=0}^K \sum_{n=-k}^k d_{lk}^{mn}(r) Y_k^n(\theta,\phi) $$ where I choose a sufficiently high truncation limit $K$, and the coefficients $d_{lk}^{mn}$ are determined from doing a truncated spherical harmonic transform of the $\alpha_l^m$. By choosing $K$ large enough, I can make the error between $\alpha_l^m$ and $\beta_l^m$ as small as I want.
Now the last problem is the functions $\beta_l^m$ are functions of radius $r$, so in a final step, I pick some radius $a$ and define my final functions, $$ \gamma_l^m(\theta,\phi) = \sum_{k=0}^K \sum_{n=-k}^k d_{lk}^{mn}(a) Y_k^n(\theta,\phi) $$ There isn't a large variation of the $\alpha_l^m$ with radius, so I can pick $a$ from a wide range of radii and still keep the geometry I want in the basis functions.
Now, the functions $\gamma_l^m$ have the nice property that they contain the geometry of my coordinate system, and they are eigenfunctions of the Laplacian as they are just linear combinations of $Y_l^m$.
However, it would be even nicer if:
1. The $\gamma_l^m$ obeyed some orthogonality relationship like the spherical harmonics
2. The $\gamma_l^m$ formed a complete set of basis functions
I would like to expand other functions in terms of the $\gamma_l^m$, such as:
$$ f(\theta,\phi) = \sum_{l=0}^L \sum_{m=-l}^l f_l^m \gamma_l^m(\theta,\phi) $$
and pick a sufficiently large $L$ so that I can get a nice approximation to $f$. If my function $f$ exhibits the same symmetry as my coordinate system $\theta_q,\phi_q$, it seems in practice that I can choose a fairly small $L$ and get a good approximation, though I don't know if that will be true for any function $f$.
Also, in order to calculate the $f_l^m$ of such an expansion, I need an orthogonality relationship for the $\gamma_l^m$. However, what I get when doing the usual integration of both sides is: $$ \int f(\theta,\phi) \gamma_k^{n*}(\theta,\phi) d\Omega = 4\pi \sum_{l=0}^L \sum_{m=-l}^l f_l^m \sum_{r=0}^K \frac{1}{2r+1} \sum_{s=-r}^r d_{kr}^{ns*} d_{lr}^{ms} $$ which doesn't seem to lead to a nice way to solve for $f_l^m$.
Does anyone know of another way to construct a set of basis functions $\gamma_l^m$, which:
Contain the geometry of my coordinate system $(r,\theta_q,\phi_q)$
Are eigenfunctions of the Laplacian in spherical coordinates,
Form a complete basis set
Have a nice orthogonality relation?