El carácter del espectro del hamiltoniano cuántico

Considera la ecuación estacionaria de Schrödinger en una dimensión,
\begin{equation}{H}\psi=E\psi\hspace{0.25in}\text{con}\hspace{0.25in}{H}\equiv-\frac{\hbar^2}{2m}\frac{d^2}{dx^2}+V\end{equation} con el potencial V tal que $\displaystyle{\lim_{|x|\to\infty}V(x)=\text{cte}}$. Al tratar con problemas de juguete seguido asumimos que los eigenvalores o energía E de H están determinados por un valor entero positivo k y nos servimos denotar al k-ésimo nivel de energía como ${E_k}$. O bien, en el caso de la partícula libre por ejemplo, sabemos que aquél valor puede ser cualquier número real, o inclusive en el átomo de hidrógeno, que puede ser de ambos tipos.

Para sentar ideas, considera el siguiente esquema

uno puede asegurar que los valores ${E_1}$ y ${E_2}$ pertenecen a la parte discreta del espectro de H, sin embargo el valor ${E_3}$ (una mera etiqueta) pertenece en general a la parte continua del espectro de H.

En general se puede asegurar que el espectro de H será discreto siempre que
\begin{align}E<\lim_{|x|\to\infty}V(x)&\label{dagger}\\[0.1in]\lim_{|x|\to\infty}\psi(x)=0&\label{ddagger}\end{align} lo que a su vez no significa que, por ejemplo, si (\ref{dagger}) no se satisface, el espectro será siempre continuo. Es decir,
$$\begin{array}{lll}(\ref{dagger})\,\&\,(\ref{ddagger})&\Rightarrow&\left\{E_k\right\}_{k\in\mathbb{N}}\\[0.1in]\left\{E_k\right\}_{k\in\mathbb{N}}&\nRightarrow&(\ref{dagger})\,\&\,(\ref{ddagger})\end{array}$$ esto parece ser evidente de forma matemática, ya que ${\mathbb{Z}\subset\mathbb{R}}$, sin embargo físicamente es relevante distinguir entre los casos discreto y continuo del espectro. Cuando el espectro del hamiltoniano es discreto, se dice que los estados ${|\psi\rangle}$ son ligados (bounded states), mientras que si es continuo, se dice que los correspondientes estados son no ligados (unbounded o scattering states). Los casos en que no se satisface (\ref{ddagger}) y el espectro es discreto, se dice que se tienen estados ligados en el continuo (bounded states in the continuum), caso que seguido podría darse en nanotecnología u óptica cuántica, por ejemplo.

Una interpretación física intuitiva para la aseveración anterior puede ser una analogía con las órbitas cerradas y abiertas que se hallan en los sistemas clásicos. Supóngase un núcleo masivo y un electrón en interacción mediante un potencial coulombiano. Si se satisfacen las dos condiciones, el sistema será un átomo hidrogenoide y como es sabido, tendrá valores de energía múltiplos de un entero positivo. En cambio si alguien disparase el electrón contra el núcleo desde muy lejos y con suficiente energía cinética, la energía total sería positiva y el núcleo desviaría al electrón sin capturarlo y sin alterar su energía, que evidentemente podrá tomar cualquier valor real. El caso más simple para un espectro energético continuo es el de la partícula libre, esto es, ${V(\mathbf{r})=0,\,\forall\,\mathbf{r}\in\mathbb{R}^3}$ donde la condición de normalización ${\int\psi_m\psi_n\,d\mathbf{r}=\delta_m\delta_n}$ deja de ser válida y puede generalizarse mediante la delta de Dirac, además de que en general las eigenfunciones se salen del espacio de Hilbert ${\mathcal{L}^2}$.

Ahora bien, he asumido ambas condiciones, sin embargo la condición (\ref{ddagger}) es una puramente cuántica y se refiere a una de más características de lo que seguido llamamos funciones de onda físicamente aceptables. La condición en realidad puede aflojarse y pedir simplemente que la función de onda sea constante en ${|x|\to\infty}$. Una forma intuitiva de entender por qué es necesaria esta condición es la siguiente. Supón un sistema que satisface (\ref{dagger}) y que eliges una solución particular $\psi$ que decrece en ${x\to\infty}$. Para este propósito, sean ${x_1,\,x_2}$ con ${x_1<x_2}$ los puntos en los que E intersecta al potencial y que $\psi$ es positivo en ${x_2}$. Entonces $\psi$ debe ser apropiadamente convexa en ${x>x_2}$ de modo que (\ref{ddagger}) se satisfaga. Para que esto ocurra, de la ecuación de la forma ${\psi^{\prime\prime}=\zeta^2\psi}$, se sabe que la curvatura depende de la energía vía ${\zeta^2}$, y por tanto, debe haber sólo un cierto valor particular de la energía para esta $\psi$ tal que la curvatura sea la apropiada dado que debe satisfacer (\ref{dagger}).

Esta es básicamente la forma de identificar inmediatamente el carácter discreto del espectro energético. Ahora bien, las identificaciones físicas o la justificación como la anterior son suficientes de manera práctica, pero gracias a unas preguntas en Physics SE me he preguntado a detalle por la forma de mostrarlo matemáticamente. Lo que se me ha ocurrido de inmediato y hace ya algún tiempo, es utilizar la teoría de Sturm-Liouville.

El detalle es el siguiente, siempre que exista $\zeta$ en el conjunto resolvente de un operador simétrico T que satisface el problema ${T\varphi=\lambda\varphi}$ tal que el operador resolvente ${R(T,\zeta)}$ sea compacto y simétrico, el espectro de T será discreto ${\{\lambda_k\}_{k\in\mathbb{N}}}$. Para que esto ocurra, T debe satisfacer el problema general de Sturm-Liouville regular en forma autoadjunta ${(T-\zeta)\,v=f}$ en un intervalo ${[a,b]}$ con función de peso $\rho$, i.e. en ${\mathcal{L}^2_\rho(a,b)}$. Otras propiedades bastante importantes en este caso, son que los eigenvalores ${\lambda_k}$ son únicos y constituyen un conjunto infinito y creciente, y que el conjunto de eigenfunciones ${\{\varphi_{_k}\}_{k\in\mathbb{N}}}$ forma una base ortogonal completa del espacio ${\mathcal{L}^2_\rho(a,b)}$. Es un clavado un tanto profundo para llegar aquí y desarrollarlo podría llevar un buen espacio, de cualquier modo el lector puede acudir a cualquier libro estándar de análisis funcional o en específico de la teoría de Sturm-Liouville.

La solución explícita por método de función de Green del problema regular de Sturm-Liouville en forma autoadjunta ${(T-\zeta)\,v=f}$ es
\begin{equation}v=(T-\zeta)^{-1}f=\int_a^bG(x,s)\,f(s)\,ds\end{equation} donde la función de Green
\begin{equation}G(x,s)=\begin{cases}\frac{\varphi_b(x)\,\varphi_a(s)}{W(s)},&a\leq{s}\leq{x}\\[0.1in]\frac{\varphi_a(x)\,\varphi_b(s)}{W(s)},&x\leq{s}\leq{b}\end{cases}\end{equation} se construye con las soluciones ${\varphi_a,\,\varphi_b}$ de la ecuación homogénea ${T\varphi=0}$ que satisface las respectivas condiciones de frontera (Neumann, Dirichlet, mixtas) en ${x=a,\,b}$. De este modo el resolvente asociado a T, por definición, ${R\equiv(T-\zeta)^{-1}}$ se puede relacionar de manera explícita como
\begin{equation}R\,f=\int_a^bG(x,s)\,f(s)\,ds\end{equation} donde G se interpreta como el kernel o núcleo de R. Al mostrar que el kernel de R es de cuadrado integrable, i.e. que satisface
\begin{equation}\int_a^b\int_a^b\left|G(x,s)\right|^2ds\,dx<\infty\end{equation} uno puede mostrar que R es compacto. Para esto uno simplemente debe mostrar que G es continua en el intervalo finito ${[a,b]}$. Ahora bien, por definición del operador resolvente, uno puede hacer el álgebra y ver que
\begin{equation}R\varphi=\mu\varphi\,\Longrightarrow\,T\varphi=\left(\frac{1}{\mu}+\zeta\right)\varphi\end{equation} de donde se sigue que las eigenfunciones de T son las mismas que las de R y entonces que R también será simétrico en ${[a,b]}$, además se pueden seguir las propiedades mencionadas para los eigenvalores ${\lambda\equiv\frac{1}{\mu}+\zeta}$ tomando en cuenta que para un operador compacto y simétrico en un espacio de Hilbert $\mathcal{H}$ con producto interior ${\left\langle\cdot\,\middle|\,\cdot\right\rangle_\mathcal{H}}$ los eigenvalores forman un conjunto infinito, acotado y numerable que converge a cero, $\displaystyle{\lim_{k\to\infty}\mu_k=0}$. Así pues, si la ecuación ${H\psi=E\psi}$ puede reducirse a un problema regular de Sturm-Liouville, el espectro de H será discreto.

Al querer atacar el problema de este modo sin embargo, aparecen ciertas inconveniencias que aún no me parecen tan claras con este argumento (aunque las haya procurado defender en mis respuestas de Physics SE). La primera es que el intervalo (o dominio en general) en que la ec. de Schrödinger es válida no está restringido, como pide el problema de S-L regular. Aquí me parece crucial la restricción (\ref{dagger}), ya que en general para una energía que intersecta a V en ${x=x_1,\,x_2}$ con ${x_1<x_2}$, la función de onda tendrá la forma
\begin{equation}\psi=\begin{cases}\psi_1,&x<x_1\\\psi_2,&x_1\leq{x}\leq{x}_2\\\psi_3,&x>x_2\end{cases}\end{equation} de modo que presumo que puede emplearse ${\psi_2}$ para reducir el problema original a uno de S-L regular y obtener información del espectro de H. Para esto, H debe ser autoadjunto en ${[x_1,x_2]}$, lo que puede lograrse homogeneizando las condiciones de frontera para ${\psi_2}$, por ejemplo por interpolación de una nueva función a través de la original, que satisfaga las condiciones de frontera originales. Esto me parece suficiente así visto toscamente, sin embargo me sigue pareciendo que podría haber sutilezas o algo por el estilo; quizá pronto me decida a intentarlo con detalle. La segunda, es que originalmente H se ha establecido como un operador simétrico en general en todo el espacio en que actúa. Esta inconveniencia desaparece si las funciones de onda completas $\psi$ en general satisfacen (\ref{ddagger}). Nuevamente esto me parece dudoso, aunque a la vez, visto toscamente, cierto. De igual modo hay ciertas cuestiones quizá un tanto más complejas, como el carácter ya sea acotado o libre del operador H, pues por el teorema de Hellinger-Toeplitz un operador autoadjunto es acotado, sin embargo el hamiltoniano (la energía) en general no es acotado. En el enlace anterior a Wikipedia se puede leer que en ese caso H queda definido sólo en un subconjunto denso de ${\mathcal{L}^2}$, además de que el hamiltoniano presuntamente siempre está acodado por abajo. En fin, dejo estos detalles más sutiles al lector, por mi parte procuraré profundizar un tanto más en el tema conforme tenga oportunidad o vaya aprendiendo quizá sobre teorías un tanto más fundamentales como teoría cuántica de campos, por ejemplo.

Funciones cilíndricas de Bessel

Aprovecho el tema de la entrada anterior para compartir la solución a la ecuación de Bessel por método de Frobenius; uno usualmente ve una ecuación Bessel y no tiene que resolverla de este modo, las soluciones son funciones Bessel y listo, sólo falta tomar en cuenta detalles, sin embargo siempre es bueno degustar una u otra vez la solución completita. Se tiene la ecuación diferencial
$$u^2\frac{d^2R}{du^2}+u\frac{dR}{du}+\left(u^2-\mu^2\right)R=0$$ (de Bessel) con ${R=R(u)}$ y ${\mu=0,\pm1,\pm2,\ldots}$, que evidentemente es singular en ${u=0}$; sea ${R\equiv{u^\alpha}\varrho(u)}$, entonces ya que
$$\frac{dR}{du}=u^\alpha\varrho^\prime+\alpha{u^{\alpha-1}}\varrho$$ y también
$$\frac{d^2R}{du^2}=u^\alpha\varrho^{\prime\prime}+2\alpha{u^{\alpha-1}}\varrho^{\prime}+\alpha(\alpha-1)u^{\alpha-2}\varrho$$ se satisface
$$u^{\alpha+2}\varrho^{\prime\prime}+(2\alpha+1){u^{\alpha+1}}\rho^{\prime}+\left[\left(\alpha^2-\mu^2\right)u^\alpha+u^{\alpha+2}\right]\varrho=0$$ y cuando ${u\to0}$, siendo que ${u^\alpha}$ está definido y posiblemente es no nulo, se tiene que la ecuación se satisface con
$$\alpha=\mu$$ por tanto ${R(u)=u^\mu\varrho(u)}$. De aquí entonces, se propone la solución en serie de potencias
$$R(u)=u^\mu\varrho(u)=\sum_{k=0}^\infty{a_k}u^{k+\mu}$$ de modo que debe satisfacerse
$$\sum_{k=0}^\infty{a_k}(k+\mu)(k+\mu-1)u^{k+\mu}+\sum_{k=0}^\infty{a_k}(k+\mu)u^{k+\mu}+\left(u^2-\mu^2\right)\sum_{k=0}^\infty{a_k}u^{k+\mu}=0$$ es decir,
$$\sum_{k=0}^\infty\left[(k+\mu)(k+\mu-1)+(k+\mu)-\mu^2\right]a_ku^{k+\mu}+\sum_{k=2}^\infty{a_{k-2}}u^{k+\mu}=0$$ expandiendo para ${k=0,1}$
$$\left(2\mu+1\right)a_1u^{\mu+1}+\sum_{k=2}^\infty\left[(k+\mu)(k+\mu-1)a_k+(k+\mu)a_k-\mu^2a_k+a_{k-2}\right]u^{k+\mu}=0$$ tomando la expresión para ${k\leq1}$, se tiene que ${\left(2\mu+1\right)a_1u^{\mu+1}=0}$, entonces ${a_1=0}$. Para ${k\geq2}$ entonces,
$$a_{k}=-\frac{1}{k\left(k+2\mu\right)}a_{k-2}$$ y todos los coeficientes impares son nulos, i.e. ${a_1=a_3=\ldots=a_{2n+1}=0}$. Para coeficientes pares entonces, nótese que, usando la propia recurrencia de coeficientes,
\begin{align*}a_{2n}&=-\frac{1}{4n\left(n+\mu\right)}a_{2(n-1)}\\[0.1in]&=\frac{(-1)^2}{4^2n(n-1)(n+\mu)(n-1+\mu)}a_{2(n-2)}\\[0.1in]&=\frac{(-1)^3}{4^3n(n-1)(n-2)(n+\mu)(n-1+\mu)(n-2+\mu)}a_{2(n-3)}\\&\vdots\\&=\frac{(-1)^n}{4^nn(n-1)\cdots2\cdot1\cdot(n+\mu)(n-1+\mu)\cdots(1+\mu)}a_0\end{align*} que por definición de la función Gamma,
$$\Gamma(\eta)\equiv(\eta-1)!$$ puede escribirse como
$$a_{2n}=\frac{(-1)^n\,\Gamma(\mu+1)}{4^nn!\,\Gamma(n+\mu+1)}a_0$$ entonces
\begin{align*}R(u)&=\sum_{n=0}^\infty{a_{2n}}\,u^{2n+\mu}\\[0.1in]&=2^\mu\,\Gamma(\mu+1)\,{a_0}\sum_{n=0}^\infty\frac{(-1)^n}{n!\,\Gamma(n+\mu+1)}\left(\frac{u}{2}\right)^{2n+\mu}\end{align*} por definición, se tienen las funciones Bessel,
$$J_\mu(u)=\sum_{n=0}^\infty\frac{(-1)^n}{n!\,\Gamma(n+\mu+1)}\left(\frac{u}{2}\right)^{2n+\mu}$$ por lo que, definiendo el coeficiente de la función Bessel en ${n=0}$, precisamente como
$$\frac{1}{2^\mu\,\Gamma(\mu+1)}\equiv{a_0}$$ se tiene la solución
$$R(u)=J_\mu(u)$$ que son precisamente las llamadas funciones de Bessel de primera especie y orden $\mu$. Ahora bien, la ecuación de Bessel es de segundo orden, por lo que el método de Frobenius sólo nos ha dado una de las soluciones. Si $\mu$ tomara valores no enteros, la solución más general podría simplemente expresarse como ${R(u)=\alpha\,J_\mu(u)+\beta\,J_{-\mu}(u)}$ con
$$J_{-\mu}(u)=\sum_{n=0}^\infty\frac{(-1)^n}{n!\,\Gamma(n-\mu+1)}\left(\frac{u}{2}\right)^{2n-\mu}$$ ya que ${J_\mu,\,J_{-\mu}}$ son linealmente independientes,
$$\begin{vmatrix}J_\mu&J_{-\mu}\\J^\prime_\mu&J^\prime_{-\mu}\end{vmatrix}\neq0$$ para $\mu$ no entero. De otro modo, para $\mu$ entero positivo, podemos escribir la función Gamma en términos de un factorial, de modo que
$$J_{-\mu}(u)=\sum_{n=0}^\infty\frac{(-1)^n}{n!\,(n-\mu)!}\left(\frac{u}{2}\right)^{2n-\mu}$$ y ${(n-\mu)!}$ diverge por lo que recorriendo índices, de modo que la serie inicie en ${n=\mu}$ (es decir, reemplazando n por ${n+\mu}$),
$$J_{-\mu}(u)=\sum_{n=0}^\infty\frac{(-1)^{n+\mu}}{n!\,(n+\mu)!}\left(\frac{u}{2}\right)^{2n+\mu}$$ se obtiene la conocida relación ${J_\mu=(-1)^{\mu}J_{-\mu}}$, por ello es que sólo hemos obtenido una solución de la ecuación Bessel. En este caso se introduce la función de Bessel de segunda especie y orden $\mu$ (también llamada función de Neumann o función de Weber),
$$y_\mu(u)\equiv\lim_{p\to\mu}\left[\cot(p\pi)J_p(u)-\csc(p\pi)J_{-p}(u)\right]$$ que por regla de L'Hôpital,
$$y_\mu(u)=\frac{1}{\pi}\left[J^\prime_\mu(u)-(-1)^{\mu}J^\prime_{-\mu}(u)\right]$$ (verifica la independencia lineal con ${J_\mu}$ y nota la divergencia en el origen) con lo que la solución más general para $\mu$ entero, es
$$R(u)=\alpha\,J_\mu(u)+\beta\,y_\mu(u)$$ De cualquier modo muchas situaciones físicas, como el caso de la entrada anterior, implican condiciones de frontera que anulan la constante $\beta$. Finalmente, son llamadas funciones cilíndricas (consulta además el caso de las funciones modificadas de Bessel) porque son precisamente la solución radial de la ecuación de Laplace en coordenadas cilíndricas, o en general puede ser, como para la entrada anterior, de una ecuación tipo Poisson en coordenadas cilíndricas.

Sobredeterminación en una EDP parabólica

Considera la ecuación
$$\alpha^2f_{xx}+2\,\alpha\beta\,f_{xy}+\beta^2f_{yy}=g\hspace{1in}(1)$$ con ${g=g(x,y)}$ y coeficientes constantes. Ya que se tiene el discriminante ${\Delta\equiv(\alpha\beta)^2-\alpha^2\beta^2=0}$, la ecuación es parabólica y tiene una única curva característica dada por
$$\frac{dy}{dx}=\frac{\beta}{\alpha}\,\Longleftrightarrow\,d(\alpha{y}-\beta{x})=0$$ entonces sean ${\zeta\equiv\alpha{y}-\beta{x}}$, ${\eta=\beta{x}}$ las nuevas variables independientes, de modo que por regla de la cadena se tiene
$$\begin{array}{ll}\partial_x&=\zeta_x\partial_\zeta+\eta_x\partial_\eta=-\beta\,\partial_\zeta+\beta\partial_\eta=\beta(\partial_\eta-\partial_\zeta)\\[0.1in]\partial^2_{x}&=\beta^2\left(\partial^2_\eta-2\partial_{\eta\zeta}^2+\partial_\zeta^2\right)\\[0.1in]\partial_y&=\zeta_y\partial_\zeta+\eta_y\partial_\eta=\alpha\,\partial_\zeta\\[0.1in]\partial^2_y&=\alpha^2\partial^2_\zeta\\[0.1in]\partial^2_{xy}&=\alpha\beta\left(\partial^2_{\eta\zeta}-\partial^2_\zeta\right)\end{array}$$ con lo que entonces también se tiene simplemente que
$$\alpha^2\partial^2_x+2\alpha\beta\,\partial^2_{xy}+\beta^2\partial^2_y=\alpha^2\beta^2\partial^2_\eta$$ es decir, podemos escribir la ecuación original en la forma
$$\alpha^2\bar{f}_{xx}=\bar{g}\hspace{1in}(2)$$ con ${\bar{g}=\bar{g}(x,\zeta)\equiv{g\left[x,\alpha^{-1}(\zeta+\beta{x})\right]}}$ y de manera análoga para $\bar{f}$. La solución de esta ecuación en un intervalo ${[x_0,x]}$ es de la forma
$$\bar{f}(x,\zeta)=\alpha^{-2}\int_{x_0}^{x}\left[\int_{x_0}^s\bar{g}(\tau,\zeta)\,d\tau\right]\,ds+x\,u(\zeta)+v(\zeta)$$ con u, v, determinadas por las condiciones de frontera del problema, que ya sean tipo Dirichlet, Neumann o mixtas, sólo se admitirán dos (por supuesto en x). Lo interesante es que se trata de un problema dos dimensional que sólo precisa condiciones en una dimensión para tener solución única, si se considera la otra, el problema queda sobredeterminado. Esto no ocurre así, por ejemplo, si g incluye una primer derivada parcial con signo positivo en y, tornando la ecuación (2) en una ec. de calor, que requeriría el equivalente a una condición inicial en y.

Esto ilustra bien que puede sobredeterminarse un problema si no hay cuidado en tratar el carácter del mismo. La pregunta interesante es en realidad, ¿qué sucede si en general ${\alpha=\alpha(x,y),\;\beta=\beta(x,y)}$ en la ec. (1)? La ecuación podría resolverse únicamente de forma numérica, por lo que sería conveniente responder antes a la pregunta hecha; aparentemente nada impide emplear nuevamente la ecuación para la curva característica y regla de la cadena.

Ortogonalidad y normalización de las funciones Bessel de primera especie

Considérense constantes ${\lambda_v,\,\lambda_w}$ tales que se satisfacen las ecuaciones Bessel
$$\label{b1}r^2\frac{d^2}{dr^2}J_\mu\left(\lambda_v{r}\right)+r\frac{d}{dr}J_\mu\left(\lambda_vr\right)+\left(\lambda_v^2\,r^2-\mu^2\right)J_\mu\left(\lambda_vr\right)=0\hspace{0.25in}\ldots\hspace{0.25in}(1)$$ y
$$\label{b2}r^2\frac{d^2}{dr^2}J_\mu\left(\lambda_w{r}\right)+r\frac{d}{dr}J_\mu\left(\lambda_wr\right)+\left(\lambda_w^2\,r^2-\mu^2\right)J_\mu\left(\lambda_wr\right)=0\hspace{0.25in}\ldots\hspace{0.25in}(2)$$ con la condición de que ${\lambda_w\mathcal{R}}$ y ${\lambda_v\mathcal{R}}$ sean ceros de ${J_\mu}$. Entonces multiplicando la ec. (1) por ${J_\mu(\lambda_wr)}$ y la ec. (2) por ${J_\mu(\lambda_vr)}$ y restando término a término,
$$(\lambda_w^2-\lambda_v^2)r J_\mu\left(\lambda_vr\right)J_\mu\left(\lambda_wr\right)=J_\mu\left(\lambda_wr\right)\frac{d}{dr}\left[r\frac{d}{dr}J_\mu\left(\lambda_vr\right)\right]-J_\mu\left(\lambda_vr\right)\frac{d}{dr}\left[r\frac{d}{dr}J_\mu\left(\lambda_wr\right)\right]$$ entonces también, integrando en $r$,
$$(\lambda_w^2-\lambda_v^2)\int_0^\mathcal{R}r J_\mu\left(\lambda_vr\right)J_\mu\left(\lambda_wr\right)dr=\int_0^\mathcal{R}J_\mu\left(\lambda_wr\right)\frac{d}{dr}\left[r\frac{d}{dr}J_\mu\left(\lambda_vr\right)\right]dr-\int_0^\mathcal{R}J_\mu\left(\lambda_vr\right)\frac{d}{dr}\left[r\frac{d}{dr}J_\mu\left(\lambda_wr\right)\right]dr$$ donde
$$\int_0^\mathcal{R}J_\mu\left(\lambda_wr\right)\frac{d}{dr}\left[r\frac{d}{dr}J_\mu\left(\lambda_vr\right)\right]dr=J_\mu(\lambda_wr)r\frac{d}{dr}J_\mu(\lambda_vr)\bigg|_0^\mathcal{R}-\int_0^\mathcal{R}r\frac{d}{dr}J_\mu\left(\lambda_wr\right)\frac{d}{dr}J_\mu\left(\lambda_vr\right)\,dr$$ y de manera análoga para el siguiente término, por lo tanto,
$$(\lambda_w^2-\lambda_v^2)\int_0^\mathcal{R}r J_\mu\left(\lambda_vr\right)J_\mu\left(\lambda_wr\right)dr=J_\mu(\lambda_wr)r\frac{d}{dr}J_\mu(\lambda_vr)\bigg|_0^\mathcal{R}-J_\mu(\lambda_vr)r\frac{d}{dr}J_\mu(\lambda_wr)\bigg|_0^\mathcal{R}$$ es decir,
$$\int_0^\mathcal{R}r J_\mu\left(\lambda_vr\right)J_\mu\left(\lambda_wr\right)dr=\frac{1}{(\lambda_w^2-\lambda_v^2)}\left[rJ_\mu(\lambda_wr)\frac{d}{dr}J_\mu(\lambda_vr)-rJ_\mu(\lambda_vr)\frac{d}{dr}J_\mu(\lambda_w{r})\right]_0^\mathcal{R}$$ que en general, para ${\lambda_v\mathcal{R},\,\lambda_w\mathcal{R}}$ dos raíces distintas de ${J_\mu}$ y ${\mu>-1}$, se verifica que el lado derecho es nulo, i.e. se prueba la ortogonalidad de las funciones Bessel de primera especie (en el espacio ${L^2\left[0,\mathcal{R};r\right]}$).

Para obtener la constante de normalización , sea ${\lambda\equiv\lambda_v}$ la única raíz local de ${J_\mu}$, de modo que por regla de L'Hôpital cuando ${\lambda_w\to\lambda}$,
\begin{align*}\lim_{\lambda_w\to\lambda}\int_0^\mathcal{R}r J_\mu\left(\lambda{r}\right)J_\mu\left(\lambda_wr\right)dr&=\int_0^\mathcal{R}r J_\mu^2\left(\lambda{r}\right)dr\\[0.1in]&=\lim_{\lambda_w\to\lambda}\frac{\mathcal{R}}{(\lambda_w^2-\lambda^2)}\left[J_\mu(\lambda_w\mathcal{R})\frac{d}{dr}J_\mu(\lambda{r})\bigg|_\mathcal{R}-J_\mu(\lambda\mathcal{R})\frac{d}{dr}J_\mu(\lambda_wr)\bigg|_\mathcal{R}\right]\\[0.1in]&\stackrel{\text{L'Hôpital}}{=}\lim_{\lambda_w\to\lambda}\frac{\mathcal{R}}{2\lambda_w}\,\frac{d}{d\lambda_w}J_\mu(\lambda_w\mathcal{R})\frac{d}{dr}J_\mu(\lambda{r})\bigg|_\mathcal{R}\end{align*} ahora bien, considérese la siguiente relación
$$\frac{d}{d\lambda_w}J_\mu(\lambda_w\mathcal{R})=\frac{\mu}{\lambda}J_\mu(\lambda_w\mathcal{R})-\mathcal{R}\,J_{\mu+1}(\lambda_w\mathcal{R})$$ que podemos obtener de la siguiente forma: La función generadora de ${J_\mu}$ está dada por
$$\mathcal{G}(\lambda{x},t)\equiv\mathrm{e}^{(\lambda{x}/2)(t-1/t)}=\sum_{\mu=-\infty}^\infty\,J_\mu(\lambda{x})t^\mu$$ cuya derivada parcial respecto a t es
\begin{align*}\frac{\partial\mathcal{G}(\lambda{x},t)}{\partial{t}}&=\frac{1}{2}\lambda{x}\left(1+\frac{1}{t^2}\right)\mathrm{e}^{(\lambda{x}/2)(t-1/t)}\\[0.1in]&=\frac{1}{2}\lambda{x}\left(1+\frac{1}{t^2}\right)\sum_{\mu=-\infty}^\infty\,J_\mu(\lambda{x})t^\mu\\[0.1in]&=\frac{1}{2}\lambda{x}\sum_{\mu=-\infty}^\infty\,\left[J_\mu(\lambda{x})t^\mu+J_\mu(\lambda{x})t^{\mu-2}\right]\\[0.1in]&=\frac{1}{2}\lambda{x}\sum_{\mu=-\infty}^\infty\,\left[J_{\mu-1}(\lambda{x})+J_{\mu+1}(\lambda{x})\right]\,t^{\mu-1}\\[0.1in]&=\sum_{\mu=-\infty}^\infty\,\mu\,J_\mu(\lambda{x})t^{\mu-1}\end{align*} La definición de la función generadora para ${J_\mu}$ surge de un desarrollo en serie de Laurent, que es único, por lo que igualando coeficientes de potencias de $t$,
$$J_{\mu-1}(\lambda{x})+J_{\mu+1}(\lambda{x})=\frac{2\mu}{\lambda{x}} J_\mu(\lambda{x})$$ De manera análoga, derivando parcialmente en x,
\begin{align*}\frac{\partial{\mathcal{G}(\lambda{x},t)}}{\partial x}&=\frac{1}{2}\lambda\left(t-\frac{1}{t}\right)\mathrm{e}^{(\lambda{x}/2)\left(t-1/t\right)}\\[0.1in]&=\frac{1}{2}\lambda\left(t-\frac{1}{t}\right)\sum_{\mu=-\infty}^\infty\,J_\mu(\lambda{x})t^\mu\\[0.1in]&=\frac{1}{2}\lambda\sum_{\mu=-\infty}^\infty\,\left[J_\mu(\lambda{x})t^{\mu+1}-J_\mu(\lambda{x})t^{\mu-1}\right]\\[0.1in]&=\frac{1}{2}\lambda\sum_{\mu=-\infty}^\infty\,\left[J_{\mu-1}(\lambda{x})-J_{\mu+1}(\lambda{x})\right]t^{\mu}\\[0.1in]&=\sum_{\mu=-\infty}^{\infty}\frac{d}{dx}\left[J_\mu(\lambda{x})\right]\,t^\mu\end{align*} se sigue que
$$J_{\mu-1}(\lambda{x})-J_{\mu+1}(\lambda{x})=\frac{2}{\lambda} \frac{d}{dx}\left[J_\mu(\lambda{x})\right]$$vrestando término a término la primer ecuación obtenida de esta segunda,
$$-2J_{\mu+1}(\lambda{x})=\frac{2}{\lambda} \frac{d}{dx}\left[J_\mu(\lambda{x})\right]-\frac{2\mu}{\lambda{x}} J_\mu(\lambda{x})$$ es decir
$$\frac{d}{dx}J_\mu(\lambda{x})=\frac{\mu}{x} J_\mu(\lambda{x})-\lambda\,J_{\mu+1}(\lambda{x})$$ $\blacksquare$

entonces se tiene
$$\int_0^\mathcal{R}r J_\mu^2\left(\lambda{r}\right)dr=-\frac{\mathcal{R}^2}{2\lambda}\,J_{\mu+1}(\lambda\mathcal{R})\frac{d}{dr}J_\mu(\lambda{r})\bigg|_\mathcal{R}$$
Considérese ahora la relación
$$J_{\mu+1}(\lambda\mathcal{R})=\frac{\mu}{\lambda\mathcal{R}}J_\mu(\lambda\mathcal{R})-\frac{1}{\lambda}\frac{d}{dr}J_\mu(\lambda{r})\bigg|_\mathcal{R}$$ que no es más que la relación utilizada anteriormente para la derivada en una variable distinta, con lo que se tiene finalmente que
$$\int_0^\mathcal{R}r J_\mu^2\left(\lambda{r}\right)dr=\frac{\mathcal{R}^2}{2\lambda^2}\,\left[\frac{d}{dr}J_\mu(\lambda{r})\bigg|_\mathcal{R}\right]^2$$ que es precisamente la constante de normalización de las funciones Bessel.

0^0=1

Quieres resolver la ecuación diferencial
$$f^{\prime\prime}(x)+f^{\prime}(x)+x^{-1}f(x)=0$$ con ${x\in\mathbb{R}}$, por lo que, sabiendo que la singularidad en ${x=0}$ es regular, propones una solución de la forma
$$f(x)=x^{\beta}F(x)$$ entonces obtienes que
$$x^{\beta}F^{\prime\prime}+\left[2\beta\,x^{\beta-1}+x^\beta\right]F^\prime+\left[\beta(\beta-1)x^{\beta-2}+(\beta+1)\,x^{\beta-1}\right]F=0$$ que multiplicado por ${x^2}$ ambos lados y cuando ${x\to0}$,
$$\beta(\beta-1)x^{\beta}F=0$$ asumiendo que, por ejemplo, ${0^{0}=1}$ y ya que ${F\neq{0}}$ en general, podemos considerar entonces el caso en que ${\beta(\beta-1)=0}$ y obtenemos la solución aceptable ${\beta=1}$; y en efecto, la solución puede hallarse en la forma ${x\,F(x)}$. Traigo esto a colación, pues en mecánica cuántica este razonamiento ha sido de ayuda para trabajar el problema de una partícula en un campo central.

De hecho me parece una manera fea de llegar a un resultado correcto. En clase se nos ha justificado el asumir que ${x^{\beta}}$ puede ser distinto de cero, ya que "cualquier número elevado a la cero es igual a uno"

Se nos ocurre notar que
$$0^0\to\lim_{z\to0^+}0^z=\lim_{z\to0^+}0=0$$ y también notamos que,
$$0^0\to\lim_{z\to0^+}z^0=\lim_{z\to0^+}1=1$$ esto es prácticamente equivalente a decir que "cero a la cualquier potencia es cero" y luego decir que "cualquier número a la potencia cero es uno". Bueno, en realidad el decir que ${0^0}$ está indeterminado es en general correcto, ya que una función ${f(x,y)=y^x}$ es discontinua en ${(0,0)}$, aunque quizá no lo sea tanto el decir que está indefinido, pues de hecho suele definirse ${0^0=1}$. La discusión es bastante antigua, y suele concluirse que la definición ${0^0=1}$ es tanto útil como consistente, y aunque seguido se considere también como un feo parche, lo usamos todo el tiempo en simplificaciones y demás al asumir ${a^0=1}$ para cualquier a.

Este tipo de cosas seguramente están muy presentes en los matemáticos, y al menos la gran mayoría de los que he conocido suelen ser muy cuidadosos con el lenguaje que emplean (recuerdo en alguna ocasión una larga discusión porque a alguien se le ocurrió decir "sacar factor común").

En el caso de interés supongo que hubiera sido más claro y suficiente aseverar que ${x^\beta}$ está definido para todo $\beta$ (al menos en $\mathbb{R}$) ignorando su valor para alguna de ellas y ${x=0}$, siendo posiblemente no nulo. De este modo, naturalmente se parte a obter los valores de ${\beta}$ que anulan la expresión, sin ser de interés particular el caso ${0^0}$. Aunque por supuesto naturalmente también a veces no somos capaces de conformarnos con simplemente ignorar el valor de una cantidad y considerarla posiblemente no nula. Acá dejo una entrada de askamathematician que he encontrado hace algún tiempo en el que se mencionan algunas razones por las cuales se define ${0^0=1}$.

Sobre el ciclo límite del oscilador de van der Pol

(...) Esta ecuación (van der Pol) es de la forma
\begin{equation}\ddot{x}+\mu(x^2-x_0^2)\,\dot{x}+{\omega_0}^2x=0\label{vdpeq}\end{equation} donde $\mu$ es un parámetro positivo de pequeño valor. Los sistemas descritos por una ecuación de van der Pol presentan una propiedad muy interesante. Si la amplitud ${|x|}$ excede el valor crítico ${|x_0|}$, el coeficiente de ${\dot{x}}$ será positivo y el sistema estará amortiguado. En cambio si ${|x|<|x_0|}$, habrá amortiguamiento negativo, es decir, la amplitud del movimiento irá en aumento. (...)

Jerry B. Marion, Dinámica clásica de las partículas y sistemas. Ed. Reverté.

Con este extracto del libro de Marion & Thornton se introduce el concepto de ciclo límite. Ésto es de la primera (segunda en inglés) edición en español (y no sé si única), en la que incluso se ilustra el ciclo límite del oscilador de VdP en el espacio fase pasando por ${\pm{x_0}}$ en ${\dot{x}=0}$. Entonces uno va, escribe unas líneas de código en Mathematica, y...

Lo que hace el siguiente código es resolver la ecuación (\ref{vdpeq}) numéricamente con ${\omega_0=1}$ implícito para dos condiciones iniciales distintas (sólo elegidas de modo que una solución comience dentro y otra fuera del ciclo límite), con $\mu=0.1$ y $x_0=1$ y mostrar el gráfico de ambas soluciones en el espacio fase.

Bueno, pues es claro que la amplitud del ciclo límite no es ${x_0}$, como se ilustra de manera errónea en la primera edición en español del libro de Marion, y como uno se precipitaría a pensar a partir del texto, sin embargo el enunciado del libro necesariamente debería ser correcto, ya que la ecuación de VdP  parece describirlo así, entonces, ¿qué es lo que está sucediendo?

Este comportamiento del oscilador VdP puede volverse accesible intuitivamente al considerar que dependiendo del signo que tome el amortiguamiento, el valor de éste dependerá también de $x$ (o más exactamente del factor $x^2-x_0^2$) y no sólo de $\mu$ como sucede en un oscilador amortiguado simple. Puede apreciarse más fácilmente cómo afecta esto y cómo se comporta el oscilador en ambos amortiguamientos al aumentar el parámetro $\mu$ en el espacio fase.

Considera el siguiente código con ${\mu=1,\;x_0=2}$

de donde puede pensarse que hay algo así como cierta simetría a partir de ${x=\pm{x_0}}$, pues es precisamente en donde el amortiguamiento cambia de 'carácter' (amortiguamiento negativo en la región roja y positivo fuera de ésta); si en estas regiones la ecuación de VdP 'tomara' un valor constante en el amortiguamiento respectivo, entonces sería clara la distinción de ambos tipos y el ciclo límite pasaría exactamente por ${x=\pm{x_0}}$ en ${\dot{x}=0}$. Así pues, aunque el enunciado en el libro es correcto, la ilustración no lo es.

Si uno se pone a jugar con el código, resulta curioso además que la consecuencia de ésto sea que el ciclo límite pase aproximadamente por ${x=\pm2x_0}$ en ${\dot{x}=0}$ (que tenga una amplitud de ${2x_0}$); de hecho en el código lo he tomado en cuenta en los dominios de los gráficos.

Considera, por ejemplo, el caso ${\mu=5,\;x_0=9}$, graficando la solución en el tiempo

la solución se va directamente al ciclo límite, que tiene apreciablemente una amplitud aprox. de ${2x_0}$. Acá sólo he dado argumentos heurísticos, pero uno seguramente podría pasárselo bien intentando realmente demostrarlo.

Yo he incluido a propósito un ejemplo con $\mu$ relativamente grande; de cualquier modo aquí se muestra que la amplitud es de hecho ${2x_0}$ cuando ${\mu\ll{1}}$, ¿puedes probar cuál es la amplitud del ciclo límite para cualquier $\mu$?