---
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} ipython3
import sys
print(sys.executable)
```

```{code-cell} ipython3
:cell_style: center
:hide_input: false

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

# Designing Potentials

Consider a single particle in a 1D potential:

\begin{gather*}
  \op{H} = \frac{\op{p}^2}{2m} + V(\op{x}), \qquad
  \frac{-\hbar^2}{2m}\psi_n''(x) + V(x)\psi_n(x) = E_n\psi_n(x).
\end{gather*}

Can we design the potential $V(x)$ to obtain a desired spectrum?

## Example Problem: Spectral Gap

:::{margin}
Stated another way, can we find a potential $V(x)$ that has an arbitrarily large gap
$\Delta$ separating the low energy states from the high energy states?
:::
For example, supposed we want to consider localized solutions, so we consider solutions
embedded in an infinite square well $V(x) = \infty$ if $x<0$ or $L < x$.
Can we design a potential $V(x)$ with two low-energy states $E_0 = 0$ and $E_1 = E$ but a
spectral gap $\Delta$ so that all higher energy states $E_{n>1} > E + \Delta$ for arbitrary large $\Delta$?

Let's start with some fundamental properties of the 1D Schrödinger equation.  We know
that we can take the eigenstates $\psi_n(x) = \braket{x|\psi} = \psi_n^*(x)$ to be real,
and that $\psi_{n}(x)$ has $n$ nodes where $\psi_n(x) = 0$.  Thus, we can
construct the potential from the ground state wavefunction and energy:

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

### Infinite Square Well

We start with $V(x) = 0$.  This has the well known solution

\begin{gather*}
  \psi_n(x) = \sqrt{\frac{2}{L}}\sin(k_nx), \qquad
  kL = \pi (1+n), \\
  E_n = \frac{\hbar^2k_n^2}{2m} = \frac{\hbar^2\pi^2 (1+n)^2}{2mL^2}.
\end{gather*}

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

import scipy.linalg 
sp = scipy

# Check orthonormality of solution.
hbar = m = L = 1.2
Nx = 200
NE = 10

# This is a nice numerical trick: we choose Nx points between 0 and L
# such that they are separated by lattice space dx, but start and end
# dx/2 from the endpoints.  Thus (Nx-1)*dx + 2*dx/2 = Nx*dx = L.
# Choosing the endpoints to be dx/2 away from 0 and L improves the accuracy
# of numerical integration.
dx = L / Nx
x = dx/2 + np.arange(Nx)*dx
n = np.arange(NE)

k_n = np.pi*(1+n)[np.newaxis, :]/L
psi_n = np.sqrt(2/L) * np.sin(k_n*x[:, np.newaxis])

# Check orthonormality
assert np.allclose(psi_n.T.conj() @ psi_n * dx, np.eye(NE))
#plt.plot(x, psi_n)

def get_H_basis(L=1.0, Nx=200, hbar2_m=hbar**2/m):
    dx = L / Nx
    x = dx/2 + np.arange(Nx)*dx
    n = np.arange(Nx)
    k_n = np.pi*(1+n)/L
    Es = (hbar2_m * k_n**2/2)
    psis = np.sqrt(2/L) * np.sin(k_n[np.newaxis, :]*x[:, np.newaxis])
    
    # Not sure why we need this, but it gives an exactly orthonormal basis.
    psis[:, -1] /= np.sqrt(2)
    return x, Es, psis * np.sqrt(dx)

# 3000 -> 2292
# 2000 -> 1515
# 1000 -> 742
# 100 -> 57
# 50 -> 23
# 20 -> 4
# 10 -> 2
# Ne = 7*Nx/9 - 30
Nx = 200
Ne = int(Nx*0.77 - 20)  # Emperical number energy states accurate to eps...
x, Es, psis = get_H_basis(Nx=Nx)
assert np.allclose(psis.T.conj() @ psis, np.eye(len(x)))

# Test with half HO
# Set length so ground state fits to machine eps
a = L / np.sqrt(-np.log(np.finfo(float).eps))

# Now divide by a factor so that the Ne state fits... explain!
# Tested emperically.
a /= np.sqrt(Nx/8)
w = hbar/m/a**2
Vx = m * (w*x)**2/2
H = np.diag(Es) + (psis.T * Vx) @ psis
E_HO = hbar*w*(2*np.arange(len(Es)) + 1 + 0.5)
plt.semilogy(abs((sp.linalg.eigvalsh(H) - E_HO) / hbar / w))
assert np.allclose(sp.linalg.eigvalsh(H)[:Ne], E_HO[:Ne])
#display(E_HO / hbar / w)
#display(sp.linalg.eigvalsh(H) / hbar / w)

def get_spectrum(V, Nx=200):
    H = np.diag(Es) + (psis.T * V(x)) @ psis
    return sp.linalg.eigvalsh(H)
```

```
from scipy.optimize import minimize
Nx = 100
x, Es, psis = get_H_basis(Nx=Nx)

Np = 10
rng = np.random.default_rng(seed=1)

bounds = np.ones((2, Np, 2))
bounds[0, :, 0] = -20.0
bounds[0, :, 1] = 20.0
bounds[1, :, 0] = 0.0
bounds[1, :, 1] = 10.0
bounds = bounds.reshape((2*Np, 2))

def pack(x, v):
    y = np.concatenate([v, dx])
    return y

def unpack(y):
    # len(y) = 2*Np
    Np = len(y) // 2
    v, dx = y[:Np], y[Np:]
    _x = np.concatenate([[0], np.cumsum(dx)]) * L / np.sum(dx)
    _v = np.concatenate([[0], v])
    return _x, _v
    
def gap(y):
    _x, _v = unpack(y)
    Vx = np.interp(x, _x, _v)
    H = np.diag(Es) + (psis.T * Vx) @ psis
    _Es = sp.linalg.eigvalsh(H)[:3]
    dEs = np.diff(_Es)
    gap = dEs[1]/dEs[0]
    print(gap, _Es)
    return -gap

y0 = rng.random(size=(2, Np))
y0[0] -= 0.5
y0[0] *= bounds[0,1]
y0[1] *= bounds[-1,1]
y0 = y0.ravel()

res = minimize(gap, x0=y0, bounds=bounds)
Vx = np.interp(x, *unpack(res.x))
H = np.diag(Es) + (psis.T * Vx) @ psis
_Es, _Vs = sp.linalg.eigh(H)
Psis = psis @ _Vs[:, :3]
plt.close('all')
plt.plot(x, Vx, x, Psis*abs(Vx.max()) + _Es[None, :3])
```

### Transfer Matrix

Now consider the potential to be a series of $S$ steps of height $V_s$ over $S$ regular intervals:

\begin{gather*}
  V(x) = \begin{cases}
    V_s & \left.\frac{sL}{S} < x < \frac{(s+1)L}{S}\right|_{s=0}^{S-1}\\
    \infty & x < 0 \text{ or } L < x.
  \end{cases}
\end{gather*}

This problem can be solved using a transfer-matrix approach.
If the potential $V(x) = V_s$ is constant for $0<x<a$, then the solution with energy
$E_s$ must satisfy:

\begin{gather*}
  \psi(x) = A e^{\I k_s x} + B e^{-\I k_s x}, \qquad
  k_s = \frac{\sqrt{2m(E_s - V_s)}}{\hbar}.
\end{gather*}

Since the Schrödinger equation is second order, to "propagate" the wavefunction
$\psi(x)$ from $x=0$ to $x=a$, we need two initial conditions at $x=0$, which we package
as an array. For example:

\begin{gather*}
  \vec{\Psi}(x) = \begin{pmatrix}
    \psi(x)\\
    \psi'(x)
  \end{pmatrix}.
\end{gather*}

Using the exact solution, and a little symbolic manipulation, we can write:

\begin{align*}
  \vec{\Psi}(x+a) &= \overbrace{\begin{pmatrix}
    \cos(k_sa) & \sin(k_sa)/k_s\\
    -k_s\sin(k_sa) & \cos(k_sa)
  \end{pmatrix}
  }^{\mat{T_s}}\vec{\Psi}(x),\\
  &= \overbrace{\begin{pmatrix}
    \cosh(\kappa_sa) & \sinh(\kappa_s a)/\kappa_s\\
    \kappa_s\sinh(\kappa_s a) & \cosh(\kappa_s a)
  \end{pmatrix}
  }^{\mat{T_s}}\vec{\Psi}(x).
\end{align*}

The second (equivalent) form is useful if $E_s < V_s$:

\begin{gather*}
  \kappa_s = \I k_s = \frac{\sqrt{2m(V_s - E_s)}}{\hbar}.
\end{gather*}

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

# Some code to test this formula.

import sympy
sympy.var('A, B, k, a, x, x_0, psi0, psia, dpsia')
dpsi0 = sympy.var(r"\psi'_0")

def psi(x):
    return A*sympy.exp(sympy.I*k*x) + B*sympy.exp(-sympy.I*k*x)

_dpsi = psi(x).diff(x)
dpsi = lambda _x: _dpsi.subs(x, _x)

eqs = [psi(x_0)-psi0, psi(x_0+a)-psia, dpsi(x_0)-dpsi0, dpsi(x_0+a)-dpsia]
res = sympy.solve(eqs, [A, B, psia, dpsia])
psia, dpsia = res[psia], res[dpsia]
T = sympy.simplify(
  sympy.Matrix([
    [psia.coeff(psi0), psia.coeff(dpsi0)],
    [dpsia.coeff(psi0), dpsia.coeff(dpsi0)]]))
display(T)
```

Since the wavefunction and derivative must be continuous across different regions, we
can transfer the wavefunction over a series of intervals by simply multiplying the
matrices $\mat{T}=\mat{T_{S-1}}\cdots\mat{T_2}\mat{T_1}\mat{T_0}$.


### Delta Functions

This transfer matrix approach can also be used with Dirac delta-function potentials
$V(x) = α δ(x)$.  Across a delta-function, the wavefunction remains
continuous, but the derivative can change, resulting in a kink.  This can be found by
integrating the Schrödinger equation across the delta-function, which we take at $x=0$
for simplicity

\begin{gather*}
  \int_{-ε}^{ε}\d{x}\left\{
    \frac{-\hbar^2}{2m}Ψ''(x) + α δ(x) Ψ(x) - EΨ(x)
  \right\} = \\
  \frac{-\hbar^2}{2m}\bigl(Ψ'(ε) - Ψ'(-ε)\bigr) + α Ψ(0) = 0.\\
  Ψ'(ε) = Ψ'(-ε) + \frac{2mα Ψ(0)}{\hbar^2}.
\end{gather*}

Hence, the transfer matrix across a delta-function is:

\begin{gather*}
  \mat{T_{δ}} = \begin{pmatrix}
    1 & 0\\
    \frac{2mα}{\hbar^2} & 1 
  \end{pmatrix}
\end{gather*}





:::{note}
The more common approach for the transfer-matrix is to express

\begin{gather*}
  \psi_s(x) = \begin{pmatrix}
    e^{\I k_s x} & e^{-\I k_s x}
  \end{pmatrix}
  \overbrace{
  \begin{pmatrix}
    A_s\\
    B_s
  \end{pmatrix}}^{\vec{\phi}_s},\\
  \psi_s'(x) = \begin{pmatrix}
    \I k_s e^{\I k_s x} & -\I k_s e^{-\I k_s x}
  \end{pmatrix}
  \vec{\phi}_s.
\end{gather*}

Matching the boundary conditions at $x=a$ between intervals $s$ and $s+1$ gives

\begin{gather*}
  \psi_s(a) = \psi_{s+1}(a), \qquad
  \psi'_s(a) = \psi'_{s+1}(a),\\
  \begin{pmatrix}
    e^{\I k_s a} & e^{-\I k_s a}
  \end{pmatrix}\vec{\phi}_s = 
  \begin{pmatrix}
    e^{\I k_{s+1} a} & e^{-\I k_{s+1} a}
  \end{pmatrix}\vec{\phi}_{s+1}, \\
  \begin{pmatrix}
    \I k_s e^{\I k_s a} & -\I k_s e^{-\I k_s a}
  \end{pmatrix}\vec{\phi}_s = 
  \begin{pmatrix}
    \I k_{s+1}e^{\I k_{s+1} a} & -\I k_{s+1}e^{-\I k_{s+1} a}
  \end{pmatrix}\vec{\phi}_{s+1}.
\end{gather*}

Solving for $\vec{\phi}_{s+1}$ we have (see {cite:p}`Walker:1992`: James S. Walker was a
professor in the WSU Physics department):

\begin{gather*}
  \vec{\phi}_{s+1} = 
  \frac{1}{2}
  \begin{pmatrix}
    1 + \frac{k_s}{k_{s+1}} & 1 - \frac{k_s}{k_{s+1}}
    1 - \frac{k_s}{k_{s+1}} & 1 + \frac{k_s}{k_{s+1}}
  \end{pmatrix}
  \vec{\phi}_{s}
  =
  \frac{1}{2}
  \begin{pmatrix}
    1 + \frac{\kappa_s}{\kappa_{s+1}} & 1 - \frac{\kappa_s}{\kappa_{s+1}}
    1 - \frac{\kappa_s}{\kappa_{s+1}} & 1 + \frac{\kappa_s}{\kappa_{s+1}}
  \end{pmatrix}
  \vec{\phi}_{s}
\end{gather*}






Then

\begin{gather*}
  \psi_s(x) = A_se^{\I k_s x} + B_se^{-\I k_s x}\\
  \psi'_s(x) = \I k_s (A_se^{\I k_s x} - B_se^{-\I k_s x}),
\end{gather*}

so we can write:

\begin{gather*}
  \vec{\Psi}_s(x) = \begin{pmatrix}
    e^{\I k_s x} &  e^{-\I k_s x}\\
    \I k_s e^{\I k_s x} &  - \I k_s e^{-\I k_s x}
  \end{pmatrix}
  \vec{\phi}_s
\end{gather*}



Thus, we can express $\psi(a)$ and $\psi'(a)$ in terms of $\psi(0)$ and $\psi'(0)$:

\begin{gather*}
  \psi(0) = A + B\\
  \psi'(0) = \I k (A - B)\\
  \psi(a) = Ae^{\I k a} + Be^{-\I k a}\\
  \psi'(a) = \I k (Ae^{\I k a} - Be^{-\I k a})
\end{gather*}

:::






If $V(x) = V_s$ is
constant for $x_0 < x < x+a$ and we know $\psi_s(x_0)$ and $\psi_s'(x_0)$, then we can
integrate the Schrödinger equation to obtain $\psi_s(x)$ throughout the interval.

To see how this works, consider the interval from $x_0 = 0$ to $x=0$.  The solution must
have the form:



\begin{gather*}
    
\end{gather*}


Let $\psi(x)$


Write the wavefunction
corresponding to energy $E=E_n$ in each interval as:

\begin{gather*}
  \psi_n(x) = \begin{cases}
    \psi_s(x) & \left.\frac{sL}{S} < x < \frac{(s+1)L}{S}\right|_{s=0}^{S-1}\\
    0 & x < 0 \text{ or } L < x.
  \end{cases}, \\
  \psi_s(x) = A_s \sin(k_s x) + B_s \cos(k_s x), \qquad
  k_s = \frac{\sqrt{2m (E_n - V_s)}}{\hbar}.
\end{gather*}


If the wavefunction is $




## Delta-Functions

Here is one solution -- a symmetric infinite square well with three delta function
potentials,

\begin{gather*}
  V(x) = \begin{cases}
    \alpha \delta\Bigl(\abs{x}-a\Bigr) + \beta \delta(x) & \abs{x} < L\\
    \infty & L < \abs{x} 
  \end{cases}, \qquad
  a = \frac{L}{2}.
\end{gather*}

This is easily solved with our transfer matrix approach:

```{code-cell} ipython3
import numpy as np
from ipywidgets import interact

a = hbar = m = 1.0
Nx = 256

E4 = hbar**2 * np.pi**2 / 2 / m / a**2

@interact(E=(0, 2*E4, 0.01))
def get_psi(E=E4, alpha=10.0, beta=20.0):
    k = np.sqrt(2*m*abs(E))/hbar
    if E < 0:
        def Ts(x):
            c, s = np.cosh(k*x), np.sinh(k*x)
            return np.array(
                [[c, s/k], 
                 [k*s, c]])
    else:
        def Ts(x):
            c, s = np.cos(k*x), np.sin(k*x)
            return np.array(
                [[c, s/k], 
                 [-k*s, c]])
    Ta = np.array(
        [[1, 0],
         [2*m*alpha/hbar**2, 1]])
    Tb = np.array(
        [[1, 0],
         [2*m*beta/hbar**2, 1]])

    Psis = [[0, 1]]
    Psis.append(Ta @ Ts(a) @ Psis[-1])
    Psis.append(Tb @ Ts(a) @ Psis[-1])
    Psis.append(Ta @ Ts(a) @ Psis[-1])
    Psis.append(Ts(a) @ Psis[-1])
    x = np.linspace(0, a, Nx//4)[:-1]
    psi = np.concatenate([np.einsum('abx,b->ax', Ts(x), _Psi)[0] 
                          for _Psi in Psis[:-1]])
    x = np.concatenate([x + _n*a for _n in range(len(Psis)-1)]) - 2*a
    psi = np.concatenate([psi, [(Ts(a) @ Psis[-2])[0]]])
    x = np.concatenate([x, [2*a]])

    fig, ax = plt.subplots()
    ax.plot(x, psi)
    ax.axhline([0], ls='--', c='k')
```

```{code-cell} ipython3

```

```{code-cell} ipython3
def get_Ts(E, x=a):
    k = np.sqrt(2*m*abs(E))/hbar
    if E == 0:
        Ts = np.array(
            [[1, x], 
             [0, 1]])
    elif E < 0:
        c, s = np.cosh(k*x), np.sinh(k*x)
        Ts = np.array(
            [[c, s/k], 
             [k*s, c]])
    else:
        c, s = np.cos(k*x), np.sin(k*x)
        Ts = np.array(
            [[c, s/k], 
             [-k*s, c]])
    return Ts

@np.vectorize
def f(E=1.0, alpha=1.0, beta=2.0):
    Ts = get_Ts(E=E, x=a)
    Ta = np.array(
        [[1, 0],
         [2*m*alpha/hbar**2, 1]])
    Tb = np.array(
        [[1, 0],
         [2*m*beta/hbar**2, 1]])

    Psi = Ts @ Ta @ Ts @ Tb @ Ts @ Ta @ Ts @ [0, 1]
    return Psi[0]

@interact
def go(alpha=10.0, beta=10.0):
    Es = np.linspace(0, 10, 256)
    fig, ax = plt.subplots()
    ax.plot(Es, f(Es, alpha=alpha, beta=beta))
    ax.axhline([0], ls='--', c='k')
    ax.axvline([E4])
    ax.set(ylim=(-1, 1))
```

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

import sympy
sympy.var('a, alpha, beta, x, hbar, m, E')
#sympy.var('E', positive=True)
#a = hbar = m = 1
hbar2_m = hbar**2/m
a = hbar2_m = 1

def get_Ts(E=E, a=a, V=0):
    k = sympy.sqrt(2*(V-E)/hbar2_m)
    c, s = sympy.cos(k*a), sympy.sin(k*a)
    return sympy.Matrix([
        [c, s/k],
        [-k*s, c]])
        
def get_Td(a):
    return sympy.Matrix([
        [1, 0],
        [2*a/hbar2_m, 1]])

Ta, Tb, Ts = get_Td(alpha), get_Td(beta), get_Ts(E, a=a)
psi_R = (Ts * Ta * Ts * Tb * Ts * Ta * Ts * sympy.Matrix([[0],[1]]))[0,0]
res = sympy.simplify(psi_R)
res
```

```{code-cell} ipython3

```

```{code-cell} ipython3

```