---
jupytext:
  formats: ipynb,md:myst
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.14.0
kernelspec:
  display_name: Python 3 (phys-555)
  language: python
  name: phys-555
---

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

import mmf_setup;mmf_setup.nbinit()
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
```

# Bound States in the 1D Schrödinger Equation

Here we consider the problem of finding the energies $E$ of bound states in a 1D potential:

\begin{gather*}
  -\frac{\hbar^2}{2m}\psi''(x) + V(x)\psi(x) = E\psi(x).
\end{gather*}

:::{margin}
Such a potential is said to have **[compact support]**.
:::
To simplify the problem, let assume that the potential is zero $V(\abs{x}>R) = 0$ for
$x$ outside of some region $-R < x < R$.  This excludes some interesting examples, but
simplifies our analysis considerably.

We will also assume that $V(x)$ is sufficiently smooth so that numerical integration
techniques will converge if we choose a sufficiently small step size.  If you want to
deal with $\delta$-functions, then you must replace these with something smooth like
gaussians and then take the limit of zero width after the computation.

## Question

Given some attractive potential $V(x)$ with compact support, how can we determine the
bound state energies?

1. How can you quickly estimate the ground state energy $E_0$?
2. How can you quickly use numerical techniques to refine this?

::::{admonition} Do It!
Stop here and this about this.  Can you come up a simple problem of this type, that you
can quickly solve?  Start with something easy, but think about how you might generalize
your methods if someone gives you a similar problem in the future.

**Hints:**

:::{dropdown} Need hints for choosing a potential $V(x)$?

A finite finite square well potential is fairly simple:

\begin{gather*}
  V(x) = \begin{cases}
    V_0 & -R < x < R,\\
    0 & \text{otherwise}.
  \end{cases}
\end{gather*}

You have probably solved this at some point by matching plane waves in the various
regions.  You may recall that determining the energies requires solving a transcendental
equation, but this can be used to check results.

However, for this problem, the existence of such an analytic solution is not important.
Instead, this potential is preferred because it will probably make integrals easy to
do.  Variations that are almost as easy include piecewise linear potentials and
piecewise polynomial potentials.  These can be used to mock up almost any potential, and
so being able to quickly work with them can be very useful.

Another option might be a smooth potential like:

\begin{gather*}
  V(x) = V_0 e^{-x^2/2\sigma^2}.
\end{gather*}

This has the advantage of having a nice simple mathematical form, but might lead to
[error function]s when integrating, which could be problematic.  Strictly speaking, such
a potential does not have compact support, but we can take $R \gtrapprox 9\sigma$ large
enough that $e^{-R^2/2\sigma^2} < 10^{-16}$ or so ([machine precision] in 64-bit double
precision) which effectively behaves as a compact potential.

For testing purposes, there are also a large class of exactly solvable potential, like
the [Pöschl-Teller potential] which have exact solutions.  These might be a good start
too, though the algebra is probably too complicated for work by hand.  Numerically, they
are a good choice.
:::

:::{dropdown} Need hints for a method?

There are many options for estimating the energies:

1. **Variational Theorem:** For the ground state, how about using the variational
   theorem? 

   \begin{gather*}
     E_0 \leq \frac{\braket{\psi|\op{H}|\psi}}{\braket{\psi|\psi}}.
   \end{gather*}
   
   Getting estimates for some exited state energies is more difficult, but can be done
   in some cases (i.e. for the first excited state $E_1$ if one has parity).

2. **Perturbation Theory:** Another possibility is to write your potential as $V(x) =
   V_0(x) + \lambda V_1(x)$ where you know the solution to $V_0(x)$ and $\lambda$ is
   small so that you can use perturbation theory.  This will not work in all cases, but
   if you are free to choose the potential, then you can definitely choose an example
   where this will work.
   
3. **WKB/Semiclassical/Bohr-Sommerfeld:**  A third possibility is to use the [WKB
   approximation] or the [Bohr-Sommerfeld quantization] condition (with a Maslov
   correction of 1/2) for the classical action for a closed trajectory $\bigl(x(t),
   p(t)\bigr)$ in classical phase space with approximate energy levels $E = E_n$:

   \begin{gather*}
     \oint p\d{q} = (n + \tfrac{1}{2})\pi \hbar
   \end{gather*}
   
   If you can solve for the classical equations of motion, this is a good option.
:::

:::{dropdown} Need hints for a trial wave-function?

If you are using the variational method, then you will need to be able to compute
$\braket{\psi|\op{H}|\psi}$.  Integrals of polynomials are easy to do, so you might like
to choose wavefunctions that are piecewise linear, quadratic, etc.  If you want to do a
better job, you might want to consider the asymptotic form of the wavefunction outside
the region of support, but remember -- the goal here is to **quickly** get an estimate,
so keep it simple.

You might try something else like

\begin{gather*}
  \psi(x) \propto \abs{x-\tilde{R}}^\alpha.
\end{gather*}

Here you could vary both $\tilde{R} \approx R$ and $\alpha$.  *(After working though
with this, however, the expressions were not so simple, so I don't think I would do this
in the future.)*

If you can use a computer to do the integrals, then you can use more complex forms.

If your potential is a gaussian, then perhaps a gaussian anzatz (with a polynomial
prefactor) might be a good choice since the product of gaussians is a gaussian, and
gaussian integrals are easy to do.
:::

:::{dropdown} Need hints for integrals?

If you choose a potential and ansatz that requires gaussian integrals, then you likely
need to compute integrals of the form

\begin{gather*}
  \int \d{x} x^n e^{-x^2/2\sigma^2}.
\end{gather*}

Here are some hints about how to do this.  It is worth being able to evaluate them quickly
without references.  A fundamental result is:

\begin{gather*}
  \int_{-\infty}^{\infty}\d{x} e^{-x^2} = \sqrt{\pi}
\end{gather*}

which can be done by using the trick of computing the square of this integral, noting
that the exponent becomes $-x^2-y^2 = r^2$ and converting the integral to a 2D integral
in cylindrical coordinates.  This result might be worth memorizing.

Need more hints? Consider dimensional analysis?  This should allow you to easily compute:
  
\begin{gather*}
  I(\sigma) = \int_{-\infty}^{\infty}\d{x} e^{-x^2/2\sigma^2}
\end{gather*}

What about $x^{n}$ where $n$ is odd?  Hint: Consider parity.

If you need a more general result, consider the derivatives

\begin{gather*}
   \diff[n]{I(\sigma)}{\sigma}
\end{gather*}

and [Feynman's integration trick] of taking the derivative under the integral.
:::
::::

## Shooting

There are many numerical techniques capable of solving the Schrödinger equation, but
shooting is one of the easiest.  The idea is to start with a guess for the wavefunction
$\psi(x_i)$ and its derivative $\psi'(x_i)$ at some point $x_i$, and then integrate
forward to another point $x_f$ trying match an appropriate boundary condition.  Most
numerical techniques require the differential equation to be specified as a first-order
equation $\vect{y}'(x) = f(x, \vect{y})$, so we use a vector $\vect{y} = (\psi, \psi')$

\begin{gather*}
    \diff{}{x}\overbrace{
    \begin{pmatrix}
      \psi(x)\\
      \psi'(x)
    \end{pmatrix}
    }^{\vect{y}(x)}
    =
    \overbrace{
    \begin{pmatrix}
      \psi'(x)\\
      \frac{2m}{\hbar^2}\Bigl(V(x) - E\Bigr)\psi(x)
    \end{pmatrix}      
    }^{f(x,\vect{y}(x))}.
\end{gather*}

:::{margin}
SciPy also provides [solve_bvp], but it is significantly less straightforward to use
effectively than shooting for eigenvalue problems.  [Numerov's method] can also be
extremely efficient.
:::
This is pretty easy to implement, either by hand using [Euler's method] or (better)
[Leapfrog integration], or using an adaptive routine like SciPy's [solve_ivp].  The
remaining challenges are:

1. Choose an appropriate integration range $x\in [x_i, x_f]$.
2. Choose an appropriate initial state $\vect{y}(x_i) = [\psi(x_i), \psi'(x_i)]$.
3. Choose an appropriate boundary condition at $x_f$ to "shoot for".

::::{admonition} Do It!
Try to answer find appropriate boundary conditions.

**Hints:**

:::{dropdown} Need hints for the range $x_i < x < x_f$?

Since our potential is compact on the domain $[-R, R]$, this might be an appropriate
range.  Depending on how you answer the other questions, you might need to expand this.
How does the wavefunction behave outside the range of support?  (I.e. what is it's
asymptotic form?)  Can you use this asymptotic form to determine how large to make $x_f
= nR$ so that you get an accurate answer to some precision?

*Don't be discouraged if it is not immediately obvious what to do.  I typically start
with something reasonable, then play with it a bit to get an intuitive idea of what is
happening.  Once I have a good idea, I then revisit and refine my choice.  Thus, I might
first try $x_f =2R$ or maybe $x_f = 10R$.  Later I will come back and make this more
precise.*
:::

:::{dropdown} Need hints for the initial state $[\psi(x_i), \psi'(x_0)]$?

1. It seems like we need two conditions here, but the Schrödinger equation is linear.
   How does the solution change if we start with $[\alpha \psi(x_i), \alpha \psi'(x_i)]$
   instead?
2. Is there any point at which you could set one of $\psi(x_i)=0$ or $\psi'(x_i)=0$?
   *(For some potentials this will work very well, but for others might not be simple.)*
:::

:::{dropdown} Need hints for the shooting condition?

Why are the bound state energies $E$ quantized?  How must the wavefunction behave for
large $x$ if it represents a bound state?  What happens if $E$ is not one of the
eigen-energies?
:::
::::

## A Simple Numerical Solution

:::{important}
Do not work through this section until you have thought carefully about the problem and
how you might try to solve.  Ideally, you should try to solve it yourself, but if you
are not very comfortable coding, then come here after you have thought about the math.
:::

:::{margin}
The Wikipedia entry [Pöschl-Teller potential] uses the notation $\lambda = n$ and units
where $V(x) = -\lambda(\lambda+1)\sech^2(x)/2$ and $E_n=-n^2/2$.  Work through the
dimensional analysis to check that I have not made errors in going back to physical
units here.
:::
Let's consider a the problem of a particle in the [Pöschl-Teller potential]:

\begin{gather*}
  V(x) = -\frac{n(n+1)\hbar^2}{2ma^2}\sech^2\frac{x}{a},
\end{gather*}

which has $n$ bound states with energies

\begin{gather*}
  E_n = -\frac{n^2\hbar^2}{2ma^2}.
\end{gather*}

We will use the Leapfrog integration scheme for a fixed step size.

```{code-cell}
# :tags: [hide-input, full-width]

hbar = m = a = 1.0
n = 2

V0 = hbar**2/2/m/a**2
ns = np.arange(1, n+1)
Es = -ns**2 * V0

def V(x, V0=V0, n=n):
    return -n*(n+1)*V0/np.cosh(np.asarray(x)/a)**2

dx = a/1000

xi = -10*a
xf = 10*a
Nsteps = int(np.ceil((xf-xi)/dx))

def shoot(E, xi=xi, xf=xf, Nsteps=Nsteps):
    """Return (xs, psis).
    
    Arguments
    ---------
    E : float
        Energy
    xi : float
        Initial point.
    xf : float
        Final point
    NSteps : int
        Number of steps.
    """
    dx = (xf-xi)/Nsteps

    # Initial ansatz... can be anything since the problem is linear
    psi = 0.1
    kappa = np.sqrt(-2*m*E)/hbar

    # Correct leapfrog initial guess using asymptotic form of wavefunction
    dpsi = kappa*psi*np.exp(kappa*dx/2)

    # Why does this work as well???
    #dpsi = 0

    x = xi
    xs = [xi]
    psis = [psi]
    for _n in range(Nsteps):
        ddpsi = 2*m*(V(x) - E)*psi/hbar**2
        dpsi += ddpsi * dx
        psi += dpsi * dx
        x += dx
        xs.append(x)
        psis.append(psi)

    xs, psis = map(np.asarray, (xs, psis))
    return xs, psis

fig, ax = plt.subplots()
xs = np.linspace(xi, xf, 200)
ax.plot(xs, V(xs), label="$V$")
V_min = np.min(V(xs))
for n, En in enumerate(Es):
    dE = 1.001
    xs, psis = shoot(E=En)
    max_psi = np.where(abs(V(xs)) > 0.1, abs(psis), 0).max()
    for E, ls in [(En/dE, ':'), (En, '-'),  (En*dE, '--')]:
        xs, psis = shoot(E=E)
        ax.plot(xs, E + psis*V0/max_psi, c=f"C{n}", ls =ls, 
                label=fr"$\psi$: $E={E/V0:.3f}V_0$")
    plt.axhline(En, ls='--', c='y')
    
ax.legend()
ax.set(xlabel="$x$", ylabel=r"$E + \psi(x)$", ylim=(1.1*V_min, 0.1*abs(V_min)));
```

### To Ponder

* Notice that unless we tune the energy precisely, it diverges exponentially.  Why?

[compact support]: https://en.wikipedia.org/wiki/Support_(mathematics)#Compact_support
[error function]: https://en.wikipedia.org/wiki/Error_function
[machine precision]: https://en.wikipedia.org/wiki/Machine_epsilon
[Pöschl-Teller potential]: https://en.wikipedia.org/wiki/Pöschl–Teller_potential
[Feynman's integration trick]: https://www.cantorsparadise.com/richard-feynmans-integral-trick-e7afae85e25c
[Euler's method]: https://en.wikipedia.org/wiki/Euler_method
[Leapfrog integration]: https://en.wikipedia.org/wiki/Leapfrog_integration
[solve_ivp]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html
[solve_bvp]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html
[Numerov's method]: https://en.wikipedia.org/wiki/Numerov%27s_method
[WKB approximation]: https://en.wikipedia.org/wiki/WKB_approximation
[Bohr-Sommerfeld quantization]: https://en.wikipedia.org/wiki/Old_quantum_theory