---
jupytext:
  formats: ipynb,md:myst,py:light
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.13.8
kernelspec:
  display_name: Python 3
  language: python
  name: python3
---

```{code-cell}
:tags: [hide-cell]

import mmf_setup; mmf_setup.nbinit()
import os
from pathlib import Path
FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/'
os.makedirs(FIG_DIR, exist_ok=True)
import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL)
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
try: from myst_nb import glue
except: glue = None
```

(sec:SturmLiouville)=
# 8. Sturm-Liouville Theory

Sturm-Liouville theory describes the general properties of equations like the
time-independent Schrödinger equation, establishing the rigorous connection with linear
algebra for wavefunctions.  Here I will couch things in terms of the language of quantum
mechanics.

:::{margin}
Inserting a complete set of states in the **position basis**
\begin{gather*}
  \op{1} = \int\d{y}\ket{y}\bra{y}
\end{gather*}
we have
\begin{gather*}
  \int\braket{x|\op{H}|y}
  \underbrace{\braket{y|\psi}}_{\psi(y)}\d{y} 
  = \underbrace{\braket{x|\psi}}_{\psi(x)}E.
\end{gather*}
The operators under consideration here are local in the sense that
\begin{gather*}
  \braket{x|\op{H}|y} = \delta(x-y)H\left(\diff{}{x}, x\right)
\end{gather*}
so the equation can be expressed as in the book:
\begin{gather*}
  \mathcal{L}(x)\psi(x) = \psi(x)\underbrace{\lambda}_{E}.
\end{gather*}
:::
The idea is to find the spectrum of equations with the form
\begin{gather*}
  \op{H}\ket{\psi} = \ket{\psi}E.
\end{gather*}
The basic idea is to make $\op{H}$ self-adjoint, so that has a complete set of
independent eigenvectors $\ket{n}$ with real eigenvalues $E_n$.  Complications include:
1. The bare operator may not be self-adjoint, requiring some modification or special
   boundary conditions.
2. The spectrum might have regions with continuous eigenvalues $E_{k}$ where $k$ is part
   of a continuum and the corresponding states $\ket{k}$ are not normalizable.  The
   appropriate mathematical context for this is called a **rigged Hilbert space** (see
   e.g. {cite}`Madrid:2005` or {cite}`Ballentine:2014` for details.)

We summarize the key results:

* To check if $\op{H}$ is self-adjoint, first see if it is hermitian.  If so, then check if it is
  self-adjoint by integrating by parts.  This may introduce boundary terms which spoil
  self-adjointness.
* Sometimes self-adjointness can be restored by adding a weight function:
  \begin{gather*}
    w(x)\mathcal{L}(x)\psi(x) = w(x)\psi(x)E.
  \end{gather*}
  This does not affect the eigenvalues or eigenfunctions, but, if needed, then this must
  appear in the inner product, and the eigenfunctions will only be orthogonal **with
  this weight**.
* Once a complete set of eigenstates is found, any function can be expressed as a linear
  combination:
  \begin{gather*}
    \psi(x) = \sum_{n}c_n\psi_n(x).
  \end{gather*}

All of this can be expressed in the following general formalism, if properly
interpreted:
\begin{gather*}
  \op{H}\ket{n} = \ket{n}E_n, \tag{eigenvalue equation}\\
  \braket{m|n} = \delta_{mn}, \tag{orthonormality}\\
  \sum_{n}\ket{n}\bra{n} = \op{1}, \tag{completness}\\
  \ket{\psi} = \sum_{n}\ket{n}\underbrace{\braket{n|\psi}}_{c_n}.
  \tag{expansion}
\end{gather*}

1. If a weight $w(x)$ is needed to make $\op{H}$ self-adjoint, then it should be
   incorporated into the inner product:
   \begin{gather*}
     \braket{f|g} = \int f^*(x)g(x) w(x)\d{x}.
   \end{gather*}
:::{margin}
In physics, the convention is to keep the factors of $2\pi$ with the momentum integrals
and momentum delta-functions.  In mathematics, things tend to be more symmetric with
factors of $\sqrt{2\pi}$ appearing everywhere.
:::
2. If the spectrum is continuous, then $\delta_{mn}$ should become a Dirac delta
   function and the sums $\sum_{n}$ should become integrals.  One must treat the
   normalization with a bit of care.  For example, the following are common in physics
   for position and momentum eigenstates:
   \begin{gather*}
     \int\d{x}\; \ket{x}\bra{x} = \int\frac{\d{x}}{2\pi} \ket{k}\bra{k} = \op{1},
     \tag{completeness}\\
     \braket{x|y} = \delta(x-y), \qquad
     \braket{k|q} = (2\pi)\delta(k-q), \tag{orthogonality}\\
     \braket{x|k} = e^{\I k x}.
   \end{gather*}

The general trick is to insert the appropriate factors of $\op{1}$.  For example, the
compute the Fourier transform $\psi_k = \braket{k|\psi}$ from the wavefunction $\psi(x)
= \braket{x|\psi}$ we simply insert $\op{1}$ in the appropriate form:
\begin{gather*}
  \psi_k = \braket{k|\psi} = \braket{k|\op{1}|\psi} = 
  \int\d{x}\;\underbrace{\braket{k|x}}_{\braket{x|k}^*}\braket{x|\psi} = 
  \int\d{x}\;e^{-\I k x}\psi(x).
\end{gather*}
Likewise, the inverse Fourier transform is
\begin{gather*}
  \psi(x) = \braket{x|\psi} = \braket{x|\op{1}|\psi} = 
  \int\frac{\d{k}}{2\pi}\braket{x|k}\braket{k|\psi} = 
  \int\frac{\d{k}}{2\pi} e^{\I k x}\psi_k.
\end{gather*}


## Self-Adjoint vs Hermitian

A hermitian operator $\op{H}$ is simply one whose matrix elements are hermitian
\begin{gather*}
  H_{mn} = \braket{m|\op{H}|n}, \qquad
  H_{mn} = H_{nm}^*, \qquad
  \mat{H} = \mat{H}^\dagger.
\end{gather*}
A self-adjoint operator is one that is its own adjoint, meaning 
\begin{gather*}
  \braket{\op{H}m|n} = \braket{m|\op{H}n}
\end{gather*}
for all states $\ket{m}$ and $\ket{n}$.  In a finite-dimensional vector space, all
hermitian matrices are self-adjoint, but the same need not be true in
infinite-dimensional Hilbert spaces.

For example, consider the momentum operator $\op{H} = \op{p}$ for a particle in a box
from $x\in [a, b]$.  The position space representation is
\begin{gather*}
  \braket{x|\op{p}\psi} = -\I\hbar \psi'(x).
\end{gather*}
To be self-adjoint, we must have
\begin{gather*}
  \braket{\op{p}f|g} = \braket{f|\op{p}g},\\
  \int_{a}^{b} [-\I\hbar f'(x)]^*g(x)\d{x} = \int_{a}^{b} f^*(x)[-\I\hbar g'(x)]\d{x},\\
  \int_{a}^{b} f'(x)^*g(x)\d{x} = -\int_{a}^{b} f^*(x)g'(x)\d{x},\\
  \int_{a}^{b} [f^*g]' \d{x} = f^*(a)g(a) - f^*(b)g(b) = 0,
\end{gather*}
for all functions $f(x)$ and $g(x)$.  This is integration by parts, and clearly this
will not be satisfied if $f(x)$ and $g(x)$ do not appropriately vanish at the endpoints
of the interval.

In this case, self-adjointness can be ensured by imposing Dirichlet boundary conditions
$f(a) = f(b)$ on our functions, as we do when we consider a particle in an infinite
square well.

:::::{admonition} Another Example of a Self-Adjoint Extension
Another interesting example of a problem that is not self-adjoint concerns the quantum
mechanics of a particle on the semi-infinite interval $r \in [0, \infty)$ with the
unbounded potential $V(r) = -m\alpha r^4 = -\hbar^2 a r^4/2m$.  This is not
self-adjoint for different reasons in that the solutions are unbounded, but there is a way to render this
well-behaved.

The key comes from considering the classical trajectories:
\begin{gather*}
  \ddot{q} = \diff{\dot{q}}{q}\dot{q} = 4 \alpha q^3,\\
  \frac{\dot{q}^2-\dot{q}_0^2}{2} = \alpha (q^4 - q_0^4).
\end{gather*}
As the particle falls away from the origin, we expect the velocity (and thus momentum)
to behave as
\begin{gather*}
  \dot{q} = \sqrt{\dot{q}_0^2 + 2\alpha(q^4 - q_0^4)} \sim \sqrt{2\alpha}q^2.
\end{gather*}
The probability of finding a particle in some region $\d{q}$ is proportional to the
amount of time it spends there, and hence, inversely proportional to the velocity.  We
thus expect
\begin{gather*}
  \abs{\psi(r)}^2 \propto \frac{1}{\dot{q}} \propto \frac{1}{r^2}.
\end{gather*}
Now notice that the probability in the tails is finite:
\begin{gather*}
  \int_{r_0}^{\infty}\frac{1}{r^2}\d{r} = \frac{1}{r_0}.
\end{gather*}
This suggests that we might be able to express our problem in a space of $L_2$ functions
with an appropriate boundary condition at infinity.

Equivalently, consider a particle starting (t=0) at rest $\dot{q}_0 = 0$ at position
$q=q_0$.  The previous condition of normalizing the wavefunction is equivalent to noting
that the particle gets to $q=\infty$ in **finite time!**:
\begin{gather*}
  t = \int_{q_0}^{\infty} \frac{1}{\sqrt{2\alpha}\sqrt{q^4 - q_0^4}}\d{q}
    = \frac{1}{q_0\sqrt{2\alpha}}\int_{1}^{\infty} \frac{1}{\sqrt{x^4 - 1}}\d{x}\\
    = \frac{1}{q_0\sqrt{2\alpha}}\frac{\sqrt{\pi}\Gamma(5/4)}{\Gamma(3/4)}
    = \frac{1}{q_0\sqrt{\alpha}}0.927\dots
\end{gather*}
*(Quick dimensions check: $[\sqrt{\alpha}] = \sqrt{E/MD^4} = 1/DT$, so
$[q_0\sqrt{\alpha}] = 1/T$.)*

One way to do this is to consider a limiting case of a hard wall
\begin{gather*}
  V(r) = \begin{cases}
    -m\alpha r^4 & r < R\\
    \infty & R \geq R.
  \end{cases}
\end{gather*}
For each $R$, the Sturm-Liouville problem will be self-adjoint because of the physical
boundary condition $\psi(R) = 0$.  The limit of $R\rightarrow \infty$ is interesting
because, as $R$ gets larger, the potential admits more and more bound states, but the
nature of these bound states for $r \ll R$ becomes universal.  This provides an example
of a renormalization-group limit cycle closely related to the [Efimov effect][].
:::::

[Efimov effect]: <https://en.wikipedia.org/wiki/Efimov_state>