Orthogonality of 3d diffusion equation solution base functions I have a question regarding the orthogonality of the base functions for the solution of the 3d diffusion equation. My motivation is that I want to calculate the diffusion process inside a 3d hollow sphere with rotational invariant initial conditions and no-flux conditions at the surface. Since I have not found the solution process anywhere in one place, I'm posting a sketch of the process here as introduction.
Background
The diffusion equation in 3d for rotational invariant geometries is
$$
\frac{\partial c}{\partial t} = D \left( \frac{\partial^2 c}{\partial r^2} + \frac 2 r \frac{\partial c}{\partial r} \right)
$$
with concentration $c(x,t)$, diffusion constant $D>0$, time $t \ge 0$ and radial coordinate $r \ge 0$. Using $u=r \, c$ this can be transformed to the 1d diffusion equation
$$
\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial r^2} \quad .
$$
This can be solved (amongst other) by separation of variables with the solution being
$$
u(r,t)=\rho_0 + \rho_1 r + \sum\limits_{k=0}^\infty \exp\left(-D \lambda_k^2 t \right) \big( \kappa_k \cos\left( \lambda_k r \right) + \sigma_k \sin\left( \lambda_k r \right) \big)
$$
with $\rho_0$, $\rho_1$, $\lambda_k$, $\kappa_k$ and $\sigma_k$ calculated to match the boundary and initial conditions (with $\rho_0$ being ambiguous to $\kappa_i$ for $\lambda_i=0$).
Since $u=r \, c$, $u=0$ for $r=0$ at all times, leading to $\rho_0=0$ and $\kappa_k=0$. The zero flux condition means $\partial c/\partial r = 0$ at $r=R$, and using the derivative of $c=u/r$,
$$
\frac{\partial u}{\partial r} - \frac{u}{R} = 0 \quad (r=R)
$$
is the second boundary condition. Evaluating this for the solution yields
$$
\tan\left(\lambda_k R\right) = \lambda_k R
$$
which I can solve numerically for the $\lambda_k$ ($\lambda_0 = 0$, $\lambda_1 \approx 4.5/R$, $\lambda_2 \approx 7.7/R$ and so on).
Since $\lambda_0=0$, the position-related base functions can be defined as
$$
X_i(r):=\begin{cases} r/R & i = 0 \\ \sin(\lambda_i r) & i \ge 1 \end{cases}
$$
The $X_i(r)$ are orthorgonal, i.e.
$$
\frac 2 R \int\limits_0^R X_i(r) X_j(r) = \begin{cases} \frac 2 3 & i=j=0 \\ 1 - \frac{1}{2 \lambda_i R} \sin(2 \lambda_i R) & i=j\ge1 \\ 0 & i \neq j \end{cases} \quad ,
$$
so the $\sigma_i$ can be easily computed from initial conditions $c(r,t=0)=c_0(r)$ and the problem is solved.
Question
From my (engineering) point of view, the orthogonality of the solution's base functions heavily depends on the boundary conditions, since
$$
\int\limits_0^R r \, \sin(\lambda_i r) \, \mathrm{d}r = \frac{1}{\lambda_i^2} \cos(\lambda_i R) \big(\tan(\lambda_i R) - \lambda_i R \big)
$$
and
$$
\int\limits_0^R \sin(\lambda_i r) \sin(\lambda_j r) \, \mathrm{d}r = \frac{ R \lambda_i \lambda_j } {\lambda_i^2-\lambda_j^2} \cos(\lambda_i R) \cos(\lambda_j R) \left( \frac{ \tan(\lambda_i R) }{\lambda_i R} - \frac{\tan(\lambda_j R)}{\lambda_j R} \right)
\quad .
$$
So in order to become orthorgonal, either $\lambda_i R = \pi k$ with $k \in \mathbb{Z}$, or $\tan(\lambda_i R) = \lambda_i R$ (or, when ignoring $X_0$, $\tan(\lambda_i R) / \lambda_i R = \mathrm{const}$).
So my question is: What is the pattern behind this? Or is it just mere luck that my boundary conditions evaluate in a way that the base functions are orthorgonal?
 A: There is a natural weight factor that should be incorporated into the analysis, which is $r^2$ coming from the Jacobian of the spherical coordinate system. The corresponding inner product is
$$
                 (f,g)_{r^2} = \int_{0}^{R}f(r)\overline{g(r)}r^2dr.
$$
Using this inner product, the operator
$$
      Lf = \frac{d^2f}{d r^2}+\frac{2}{r}\frac{df}{dr}
$$
is selfadjoint, once standard endpoint conditions are imposed. This is because
$$
              Lf = \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{df}{dr}\right),
$$
which gives
\begin{align}
       \int (Lf)g r^{2}dr & = \int \frac{d}{dr}\left(r^2\frac{df}{dr}\right)gdr \\
    & = r^{2}\frac{df}{dr}g-\int r^{2}\frac{df}{dr}\frac{dg}{dr}dr \\
    & = r^{2}\frac{df}{dr}g-r^{2}f\frac{dg}{dr}+\int f\frac{d}{dr}\left(r^2\frac{dg}{dr}\right)dr \\
    & = r^{2}\{ f'g-fg' \}+\int f(Lg)r^2dr
\end{align}
Therefore,
$$
     (Lf,g)_{r^2} - (f,Lg)_{r^2} = r^2\{f'\overline{g}-f\overline{g}'\}|_{0}^{R}
$$
Once proper endpoint conditions at the endpoints for the functions $f,g$ in the domain,
$$
                 (Lf,g)_{r^2}=(f,Lg)_{r^2}.
$$
When dealing with a weighted space like this, $rf \in L^2(0,R)$ with unit weight iff $f\in L^2_{r^2}(0,R)$, which motivates the substitution of the original equation that leads to an unweighted space. The resulting operator is then selfadjoint in $L^2(0,R)$ with unit weight. The transformation $f\in L^2_{r^2}(0,R)\mapsto rf \in L^2(0,R)$ is a unitary substitution between the Hilbert spaces $L^2_{r^2}(0,R)$ and $L^2(0,R)$. The resulting operator is also selfadjoint in $L^2(0,R)$ once the endpoint conditions are translated. The new ODE in $L^2(0,R)$ is sandwiched between multiplication by $1/r$ and multiplication by $r$. The new operator $M : \mathcal{D}(M)\subset L^2(0,R)\rightarrow L^2(0,R)$ is
$$
                M = rL(\frac{1}{r}f),\;\;\; f \in L^2(0,R)
$$
Explicitly,
\begin{align}
Mf & = r\left[\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{d(\frac{1}{r}f)}{dr}\right)\right] \\
   & = r\left[\frac{1}{r^2}\frac{d}{dr}\left(r\frac{df}{dr}-f\right)\right] \\
   & = r\left[\frac{1}{r^2}\left(r\frac{d^2f}{dr^2}+\frac{df}{dr}-\frac{df}{dr}\right)\right] \\
   & = \frac{d^2f}{dr^2}.
\end{align}
This also exposes the possible endpoint conditions at $r=0$. In the transformed coordinates, the left endpoint condition functionals are $\Phi(f)=f(0+)$ and $\Psi(f)=f'(0+)$. That is, $\{\Phi,\Psi\}$ is a basis of left endpoint functionals. These functions may be applied to functions on the original space to obtain the transformed basis of endpoint functionals $\{\Phi_{r^2},\Psi_{r^2}\}$ defined by
\begin{align}
    \Phi_{r^2}(g) & = \lim_{r\rightarrow 0}(rg(r)) \\
    \Psi_{r^2}(g) & = \lim_{r\rightarrow 0}(\frac{d}{dr}(rg))
     =\lim_{r\rightarrow 0}(rg'(r)+g(r)).
\end{align}
To obtain a selfadjoint problem for $L^2(0,R)$ with separated endpoint conditions, you can impose one endpoint homogeneous condition at each endpoint. This is standard. The only natural condition at $r=0$ is $f(0+)=0$ because $g,Lg \in L^2_{r^2}$ implies that $\Phi_{r^2}(g)$ exists and
$$
                    g(r) \sim \Phi_{r^2}(g)\frac{1}{r}.
$$
Because $L$ is non-singular at $r=R$, the evaluation functionals for $g$ and $g'$ at $r=R$ is a basis of endpoint functionals, and these are equivalent to the the evaluations of $f$ and $f'$ at $r=R$. Let $L_{\alpha}$ denote the restriction of $L$ to the domain with $f,Lf\in L^2_{r^2}(0,R)$ satisfying endpoint conditions
$$
           \lim_{r\rightarrow 0}(rf(r))=0,\;\;\; \cos\alpha f(R)+\sin\alpha f'(R)=0.
$$
Then $(Lf,g)_{r^2} = (f,Lg)_{r^2}$ for all $f,g \in \mathcal{D}(L_{\alpha})$. It is such symmetry that guarantees the orthogonality of eigenfunctions with different eigenvalues, just as it does for Hermitian matrices; more precisely, if $Lf=\lambda f$, $Lg=\mu g$ with $f,g\in\mathcal{D}(L_{\alpha})$ and $\lambda\ne \mu$, then $\mu,\lambda$ are real because
$$
    (\lambda-\overline{\lambda})\|f\|_{r^2}=(Lf,f)_{r^2}-(f,Lf)_{r^{2}}=0,\\
    (\mu-\overline{\mu})\|g\|_{r^2}=(Lg,g)_{r^2}-(g,Lg)_{r^2}=0,
$$
and $(f,g)_{r^2}=0$ because
$$
          (\lambda-\mu)(f,g)_{r^2}=(Lf,g)_{r^2}-(f,Lg)_{r^2}=0.
$$
