r"""Numerical Solutions to the Schrodinger equation.
This module contains some exact solutions to the Schrodinger equation
for use in testing.
"""
from __future__ import division
__all__ = ['ScarfII', 'S', 'R', 'HO_radial']
import warnings
from warnings import warn
import numpy as np
import scipy as sp
import scipy.special
from mmf.objects import StateVars, process_vars
from mmf.objects import ClassVar, Excluded, Computed
try:
from pygsl.sf import laguerre_n
except ImportError, err:
warnings.warn(
"Could not import PyGSL: %s" % (str(err),))
_eps = np.finfo(float).eps
[docs]class S(StateVars):
r"""Scarf II (hyperbolic) potential:
.. no_numpydoc
::
ScarfII(A=1,
B=0,
a=1,
m=1,
hbar=1)
.. math::
V(x) = \frac{\hbar^2}{2m a^2}\Biggl[A^2 + \left(B^2 - A^2 -
A\right)
\sech^{2}\frac{x}{a} + B\left(2A + 1\right)
\sech\frac{x}{a}\; \tanh \frac{x}{a}\Biggr].
.. rubric:: Notes
The depth of the potential is controlled by :attr:`A` while the
asymmetry is controlled by `B` (though the spectrum is unchanged).
.. plot::
:include-source:
from mmf.physics.seq import ScarfII
x = np.linspace(-3,3,100)
for B in [-1,0,1,2]:
V = ScarfII(A=1, B=B)
plt.plot(x, V(x), label='B=%i'%(B,))
plt.legend()
.. rubric:: Examples
.. plot::
:include-source:
from mmf.physics.seq import ScarfII
V = ScarfII(A=3, B=2)
x = np.linspace(-10, 10, 200)
v = V(x)
En = V.E()
plt.plot(x, v, 'k')
for n, E in enumerate(En):
x_ = x[np.where(v <= E)][[0, -1]]
line = plt.plot(x_, x_*0 + E, lw=2)
colour = line[0].get_color()
psi = V.psi(n)(x)
psi = psi/np.max(abs(psi))
plt.plot(x, 2*psi, colour + '-')
.. rubric:: Attributes
=========================== ==========
**State Variables:**
.. attribute:: ScarfII.A
Dimensionless parameter (depth).
.. attribute:: ScarfII.B
Dimensionless parameter (controls asymmetry).
.. attribute:: ScarfII.a
Length scale
.. attribute:: ScarfII.m
<no description>
.. attribute:: ScarfII.hbar
<no description>
=========================== ==========
.. rubric:: Methods
.. autosummary::
:toctree:
E
all_items
archive
archive_1
delta_0
initialize
items
psi
resume
suspend
"""
_state_vars = [
('A', 1, "Dimensionless parameter (depth)."),
('B', 0, "Dimensionless parameter (controls asymmetry)."),
('a', 1, "Length scale"),
('m', 1),
('hbar', 1),
]
process_vars()
def __call__(self, x):
A, B = self.A, self.B
z = x/self.a
return (self.hbar/self.a)**2/2/self.m*(
A**2 + ((B**2 - A**2 - A)
+ B*(2*A + 1)*np.sinh(z))/np.cosh(z)**2)
def E(self):
r"""Return the energy eigenvalues.
.. math::
E_n = \frac{\hbar^2}{2m a^2}\left[A^2 - (A - n)^2\right]
"""
E = self.A**2 - (self.A - np.arange(self.A))**2
return (self.hbar/self.a)**2/2/self.m*E
def psi(self, n):
r"""Return the n'th wavefunction.
.. math::
\psi_{n} &= a^{-1/2}(1+y^2)^{-A/2}e^{-B\tan^{-1}y}
R_{n}^{(A +1/2, -2B)}(y)\\
y &= \sinh\frac{x}{a}.
"""
def f(x):
y = np.sinh(x/self.a)
return (1/np.sqrt(self.a)/np.power(1 + y*y, self.A/2)
*np.exp(-self.B*np.arctan(y))
*R(y, n, a=self.A +1/2, b=-2*self.B))
return f
def delta_0(self, k):
r"""Return the s-wave phase shift as a function of wavevector
k.
This comes from the Flugge book 'Practical Quantum
Mechanics'. The form there is
.. math::
V(r) &= -\frac{\hbar^2\alpha^2}{2m}\lambda(\lambda - 1)
\sech^2(\alpha r)\\
\delta_0 &= \frac{\pi}{2} + \arg Q, \\
Q &= \frac{
\Gamma(ik/\alpha)2^{-ik/\alpha}}
{\Gamma\left(
\frac{\lambda + 1}{2} + i\frac{k}{2\alpha}\right)
\Gamma\left(
1 - \frac{\lambda}{2} + i\frac{k}{2\alpha}\right)
}
Note that the mass here is not the reduced mass. This is the
scattering length obtained from the Hamiltonian:
.. math::
\op{H} = \frac{-\hbar^2\nabla^2}{2m} + V(x)
whereas for the two-body problem, one solves
.. math::
\op{H} = \frac{-\hbar^2\nabla^2}{2M} + V(x)
where `M = m/2`. One must relate the coefficients with the
potential appropriately in order to equate coefficients.
To get the effective range expansion we note that
.. math::
k \cot\delta_0 = -k\tan\arg{Q} = -k\frac{\im Q}{\re Q}
= \frac{-1}{a_s} + \frac{r_{\text{eff}} k^2}{2} + \cdots
The scattering length, for example, follows from the `k=0`
limit of this expression and is computed by noting that
.. math::
\Gamma(2i0^{+}) &=
\gamma + \order(0^{+2})
- i\left(\frac{1}{2\;0^{+}} + \order(0^{+})\right)\\
\Gamma(x + i0^{+}) &= \Gamma(x) + i0^{+}\Gamma'(x) + \order(0^{+2})
Thus, in the limit `k=0` we have
.. math::
\frac{1}{a_s} = \lim_{k\rightarrow 0} - k\cot \delta_0
= \frac{\alpha}
{\gamma + \ln{2}
+ \tfrac{1}{2}\psi\left(\frac{\lambda + 1}{2}\right)
+ \tfrac{1}{2}\psi\left(1-\frac{\lambda}{2}\right)}
where $\gamma = -\psi(1)$ is the Euler-Mascheroni
constant and $\psi(x) = \Gamma'(x)/\Gamma(x)$ is the
digamma function.
To make contact with our normalizations we must equate
.. math::
\begin{aligned}
\alpha &= 1/a, &
\lambda &= A + 1, &
B &= 0
\end{aligned}
Examples
--------
The potential has resonances :math:`a_s = \infty` and
effective ranges :math:`r_{\text{eff}}` with `n`
bound states (the last one at threshold) for
.. math::
\begin{aligned}
A_n &= 2n -1, &
r_{\text{eff}} &= 2a H_{2n-1}, &
H_{2n-1} &= \sum_{j=1}^{2n-1}\frac{1}{j}
\end{aligned}
>>> V = ScarfII(A=1.0, B=0, a=2)
>>> def kcotdelta(k): return k/np.tan(V.delta_0(k))
>>> abs(kcotdelta(1e-15)) < 1e-20
True
"""
if self.B != 0:
warn("Scattering properties only known for B=0" +
"(they have not been checked for nonzero B)")
lam = self.A + 1
alpha = 1/self.a
ik = 1j*k
g = sp.special.gamma
psi = sp.special.digamma
gamma = -psi(1)
ainv = alpha/(
gamma + np.log(2)
+ 0.5*psi((lam + 1)/2)
+ 0.5*psi(1 - lam/2))
Q = (g(ik/alpha)*np.exp(-np.log(2)*ik/alpha)
/g((lam + 1)/2 + ik/2/alpha)
/g(1 - lam/2 + ik/2/alpha))
delta_0 = np.pi/2 + np.angle(Q)
return delta_0
return -k/np.tan(delta_0), ainv
def R(x, n, a, b):
r"""Return the Romanovski polynomial
.. math::
R_{m}^{(a,b)}(x) &= \frac{1}{w^{(a,b)}(x)}
\diff[m]{}{x}(1+x^2)^{m}w^{(a,b)}(x)\\
w^{(a,b)}(x) &= (1+x^2)^{-a}e^{b\tan^{-1}(x)}
Notes
-----
Computed using the recurrence
.. math::
R_{m}^{(a,b)}(x) &= [b + 2(m - a)x]R_{m-1}^{(a,b)}(x)
+ 2(m - 1)(m - a)(1+x^2)R_{m-2}^{(a-1,b)}(x),\\
R_{0}^{(a,b)} &= 1,\\
R_{1}^{(a,b)} &= 2(1-a)x + b.
"""
if n < 0:
raise ValueError('Input `n` must be a non-negative integer')
elif 0 == n:
return x*0 + 1
elif 1 == n:
return b - 2*(a - 1)*x
else:
R1 = R(x, n-1, a, b)
R2 = R(x, n-2, a-1, b)
return (2*(n - a)*x + b)*R1 + 2*(n - 1)*(n - a)*(1 + x*x)*R2
[docs]class ScarfII(StateVars):
r"""Scarf II (hyperbolic) potential:
.. math::
V(x) = \frac{\hbar^2}{2m a^2}\Biggl[A^2 + \left(B^2 - A^2 -
A\right)
\sech^{2}\frac{x}{a} + B\left(2A + 1\right)
\sech\frac{x}{a}\; \tanh \frac{x}{a}\Biggr].
Notes
-----
The depth of the potential is controlled by :attr:`A` while the
asymmetry is controlled by `B` (though the spectrum is unchanged).
.. plot::
:include-source:
from mmf.physics.seq import ScarfII
x = np.linspace(-3,3,100)
for B in [-1,0,1,2]:
V = ScarfII(A=1, B=B)
plt.plot(x, V(x), label='B=%i'%(B,))
plt.legend()
Examples
--------
.. plot::
:include-source:
from mmf.physics.seq import ScarfII
V = ScarfII(A=3, B=2)
x = np.linspace(-10, 10, 200)
v = V(x)
En = V.E()
plt.plot(x, v, 'k')
for n, E in enumerate(En):
x_ = x[np.where(v <= E)][[0, -1]]
line = plt.plot(x_, x_*0 + E, lw=2)
colour = line[0].get_color()
psi = V.psi(n)(x)
psi = psi/np.max(abs(psi))
plt.plot(x, 2*psi, colour + '-')
"""
_state_vars = [
('A', 1, "Dimensionless parameter (depth)."),
('B', 0, "Dimensionless parameter (controls asymmetry)."),
('a', 1, "Length scale"),
('m', 1),
('hbar', 1),
]
process_vars()
def __call__(self, x):
A, B = self.A, self.B
z = x/self.a
return (self.hbar/self.a)**2/2/self.m*(
A**2 + ((B**2 - A**2 - A)
+ B*(2*A + 1)*np.sinh(z))/np.cosh(z)**2)
def E(self):
r"""Return the energy eigenvalues.
.. math::
E_n = \frac{\hbar^2}{2m a^2}\left[A^2 - (A - n)^2\right]
"""
E = self.A**2 - (self.A - np.arange(self.A))**2
return (self.hbar/self.a)**2/2/self.m*E
def psi(self, n):
r"""Return the n'th wavefunction.
.. math::
\psi_{n} &= a^{-1/2}(1+y^2)^{-A/2}e^{-B\tan^{-1}y}
R_{n}^{(A +1/2, -2B)}(y)\\
y &= \sinh\frac{x}{a}.
"""
def f(x):
y = np.sinh(x/self.a)
return (1/np.sqrt(self.a)/np.power(1 + y*y, self.A/2)
*np.exp(-self.B*np.arctan(y))
*R(y, n, a=self.A +1/2, b=-2*self.B))
return f
def delta_0(self, k):
r"""Return the s-wave phase shift as a function of wavevector
k.
This comes from the Flugge book 'Practical Quantum
Mechanics'. The form there is
.. math::
V(r) &= -\frac{\hbar^2\alpha^2}{2m}\lambda(\lambda - 1)
\sech^2(\alpha r)\\
\delta_0 &= \frac{\pi}{2} + \arg Q, \\
Q &= \frac{
\Gamma(ik/\alpha)2^{-ik/\alpha}}
{\Gamma\left(
\frac{\lambda + 1}{2} + i\frac{k}{2\alpha}\right)
\Gamma\left(
1 - \frac{\lambda}{2} + i\frac{k}{2\alpha}\right)
}
Note that the mass here is not the reduced mass. This is the
scattering length obtained from the Hamiltonian:
.. math::
\op{H} = \frac{-\hbar^2\nabla^2}{2m} + V(x)
whereas for the two-body problem, one solves
.. math::
\op{H} = \frac{-\hbar^2\nabla^2}{2M} + V(x)
where `M = m/2`. One must relate the coefficients with the
potential appropriately in order to equate coefficients.
To get the effective range expansion we note that
.. math::
k \cot\delta_0 = -k\tan\arg{Q} = -k\frac{\im Q}{\re Q}
= \frac{-1}{a_s} + \frac{r_{\text{eff}} k^2}{2} + \cdots
The scattering length, for example, follows from the `k=0`
limit of this expression and is computed by noting that
.. math::
\Gamma(2i0^{+}) &=
\gamma + \order(0^{+2})
- i\left(\frac{1}{2\;0^{+}} + \order(0^{+})\right)\\
\Gamma(x + i0^{+}) &= \Gamma(x) + i0^{+}\Gamma'(x) + \order(0^{+2})
Thus, in the limit `k=0` we have
.. math::
\frac{1}{a_s} = \lim_{k\rightarrow 0} - k\cot \delta_0
= \frac{\alpha}
{\gamma + \ln{2}
+ \tfrac{1}{2}\psi\left(\frac{\lambda + 1}{2}\right)
+ \tfrac{1}{2}\psi\left(1-\frac{\lambda}{2}\right)}
where $\gamma = -\psi(1)$ is the Euler-Mascheroni
constant and $\psi(x) = \Gamma'(x)/\Gamma(x)$ is the
digamma function.
To make contact with our normalizations we must equate
.. math::
\begin{aligned}
\alpha &= 1/a, &
\lambda &= A + 1, &
B &= 0
\end{aligned}
Examples
--------
The potential has resonances :math:`a_s = \infty` and
effective ranges :math:`r_{\text{eff}}` with `n`
bound states (the last one at threshold) for
.. math::
\begin{aligned}
A_n &= 2n -1, &
r_{\text{eff}} &= 2a H_{2n-1}, &
H_{2n-1} &= \sum_{j=1}^{2n-1}\frac{1}{j}
\end{aligned}
>>> V = ScarfII(A=1.0, B=0, a=2)
>>> def kcotdelta(k): return k/np.tan(V.delta_0(k))
>>> abs(kcotdelta(1e-15)) < 1e-20
True
"""
if self.B != 0:
warn("Scattering properties only known for B=0" +
"(they have not been checked for nonzero B)")
lam = self.A + 1
alpha = 1/self.a
ik = 1j*k
g = sp.special.gamma
psi = sp.special.digamma
gamma = -psi(1)
ainv = alpha/(
gamma + np.log(2)
+ 0.5*psi((lam + 1)/2)
+ 0.5*psi(1 - lam/2))
Q = (g(ik/alpha)*np.exp(-np.log(2)*ik/alpha)
/g((lam + 1)/2 + ik/2/alpha)
/g(1 - lam/2 + ik/2/alpha))
delta_0 = np.pi/2 + np.angle(Q)
return delta_0
return -k/np.tan(delta_0), ainv
[docs]def R(x, n, a, b):
r"""Return the Romanovski polynomial
.. math::
R_{m}^{(a,b)}(x) &= \frac{1}{w^{(a,b)}(x)}
\diff[m]{}{x}(1+x^2)^{m}w^{(a,b)}(x)\\
w^{(a,b)}(x) &= (1+x^2)^{-a}e^{b\tan^{-1}(x)}
Notes
-----
Computed using the recurrence
.. math::
R_{m}^{(a,b)}(x) &= [b + 2(m - a)x]R_{m-1}^{(a,b)}(x)
+ 2(m - 1)(m - a)(1+x^2)R_{m-2}^{(a-1,b)}(x),\\
R_{0}^{(a,b)} &= 1,\\
R_{1}^{(a,b)} &= 2(1-a)x + b.
"""
if n < 0:
raise ValueError('Input `n` must be a non-negative integer')
elif 0 == n:
return x*0 + 1
elif 1 == n:
return b - 2*(a - 1)*x
else:
R1 = R(x, n-1, a, b)
R2 = R(x, n-2, a-1, b)
return (2*(n - a)*x + b)*R1 + 2*(n - 1)*(n - a)*(1 + x*x)*R2
[docs]class HO_radial(StateVars):
r"""Radial Harmonic Oscillator.
Notes
-----
In :attr:`d`-dimensions, the Hamiltonian is simply a sum over :attr:`d`
independent oscillators, $\smash{\op{H}}^{(d)} = \sum_{i=1}^{d}
\op{H}_{i}$, where each $\op{H}_{i}$ is an independent oscillator
in coordinate $x_{i}$. The energies are additive and thus:
.. math::
E_{n_{1},n_{2},\cdots,n_{d}} = \left(\sum_{i} n_{i} +
\frac{d}{2}\right)\hbar\omega.
It is useful to consider this as a central potential problem in
the coordinate $r$. In this case we recover a 1-d problem with
an additional set of quantum numbers corresponding to angular
momentum $l$. This follows from the expression of the Laplacian
in d-dimensions:
.. math::
\nabla^{\field{R}^d} = \frac{1}{r^{d-1}}\diff{}{r}r^{d-1}\diff{}{r} +
\frac{1}{r^2}\nabla^{S^{d-1}}
where $\nabla^{S^{d-1}}$ is the Laplacian defined on the unit sphere
$S^{d-1}$. This can be brought into a simple form by introducing
the radial momentum operator:
.. math::
\op{p}_{r} = \frac{-i}{r^{(d-1)/2}}\pdiff{}{r}r^{(d-1)/2}.
In terms of this we have:
.. math::
\frac{1}{r^{d-1}}\diff{}{r}r^{d-1}\diff{}{r} = -\op{p}_{r}^2 -
\frac{(d-1)(d-3)}{4r^2}.
Finally, introducing the appropriate :attr:`d`-dimensional
spherical harmonics $\mathcal{Y}_{l}(\vect{r})$ we have
.. math::
\frac{1}{r^2}\nabla^{S^{d-1}}\mathcal{Y}_{l}(\vect{r}) =
-l(l+d-2)\mathcal{Y}_{l}(\vect{r}).
Combining these gives the effective radial
Hamiltonian\footnote{One can define raising a lowering operators
to analyse this~\cite{Arik:2008ix}.}
.. math::
\begin{aligned}
\smash{\op{H}}^{(d)} &= \frac{1}{2m}\left(
\op{p}_{r}^2 + \frac{\nu_{l}^2 - \frac{1}{4}}{\op{r}^2}
+ m^2\omega^2 \op{r}^2 \right), &
\nu_{l} &= l + \frac{d-2}{2}.
\end{aligned}
The energy levels $E_{l,j}$ and degeneracies $g_{l,j}$ are thus
.. math::
\begin{aligned}
E_{l,j} &= \left(l + 2j + \frac{d}{2}\right)\hbar\omega, &
g_{l,j} &= \binom{d+l+2j-1}{l+2j}
= \frac{\Gamma(d+l+2j)}{\Gamma(l+2j+1)\Gamma(d)}.
\end{aligned}
with normalized radial eigenfunctions are $\Psi^{(l)}_{j}(r) =
R_{j}(r/L)$ where
.. math::
R^{l}_{j}(r) = (-1)^{j}\sqrt{\frac{2\Gamma(j+1)}{\Gamma(l+j+d/2)}}
e^{-r^2/2} r^{\nu+1/2}L_{j}^{\nu_l}(r^2)
and $L_{j}^{\nu_l}$ are the generalized Laguerre polynomials.
"""
_state_vars = [
('d', 3, 'Dimension of enclosing space.'),
('hbar', 1, r'$\hbar$'),
('m', 1, "Mass"),
('omega', 1, "Oscillator frequency $\omega$."),
('r0', Computed, "Oscillator length"),
]
process_vars()
def __init__(self, *v, **kw):
self.r0 = np.sqrt(self.hbar/self.m/self.omega)
def E(self, n, l):
"""Exact energies for the radial harmonic oscillator in
`d` dimensions with angular quantum number `l`."""
return self.hbar*self.omega*(2*n + l + self.d/2)
def __call__(self, r):
"""Potential."""
return r*r*self.m*self.omega**2/2
def psi(self, n, l):
r"""Return the exact eigenstates $\psi_{n,l}(r) =
R^{l}_{n}(r/r_0)$ as a function of `r`.
The HO eigenstates are related as:
.. math::
\Psi_{n,l}(r) = r^{-d/2}\psi_{n,l}(r)
Notes
-----
These are expressed in terms of the generalized Laguerre
polynomials $L(x) = L_{j}^{\nu_l}(x)$ which satisfy
.. math::
\psi(r) = r_0^{-1/2}R^{l}_{n}(r/r_0),\\
xL''(x) + (\nu_l + 1 -x)L'(x) + nL(x) = 0,\\
R^{l}_{n}(\tilde{r}) = (-1)^{n}\sqrt{\frac{2\Gamma(n+1)}{\Gamma(n + \nu_l + 1)}}
e^{-\tilde{r}^2/2} \tilde{r}^{\nu_l+1/2}L_{j}^{\nu_l}(\tilde{r}^2)
where $r_0 = \sqrt{\hbar/(m\omega)}$ is the oscillator
radius and $\nu_l = l + d/2 - 1$. The normalization is such
that
.. math::
\int\d^{3}{\vect{r}}\; \abs{\Psi(r)}^{2} =
\int_{0}^{\infty}\d{r}\;S_{d}\abs{\psi(r)}^{2} =
\int_{0}^{\infty}\d{r}\;\abs{\psi(r)}^{2} = 1.
Examples
--------
>>> import scipy.integrate
>>> V = HO_radial(m=1.2,omega=2.3,hbar=3.4,d=4.2)
>>> psi = V.psi(n=1,l=1)
>>> sp.integrate.quad(lambda r:psi(r)**2,0,np.inf)[0] # doctest: +ELLIPSIS
1.0000000...
"""
nu = l + (self.d - 2)/2
#coeff = (-1)**n*sp.sqrt(2.0*sp.misc.factorial(n)/
# sp.misc.factorial(n + nu))/np.sqrt(self.r0)
coeff = (-1)**n*sp.sqrt(2.0/self.r0*
np.exp(sp.special.gammaln(n+1)
- sp.special.gammaln(n + nu + 1)))
@np.vectorize
def R(r, r0=self.r0):
r_ = r/r0
r2 = r_*r_
return (coeff*np.exp(-r2/2 + (nu + 0.5)*np.log(r_ + _eps))*\
laguerre_n(n, nu, r2)[0])
return R