Radial Wavefunctions and Quantum Defects
========================================
In this tutorial we show how to access quantum defects and
wavefunctions, which are used for the computation of matrix elements,
using the Python API. Some aspects of this are discussed in Appendix A
of the pairinteraction paper `J. Phys. B: At. Mol. Opt. Phys. 50, 133001
(2017) `__.
This feature is mostly used internally and therefore the interfaces are
not so user-friendly.
Ahead of our Python code, we call an IPython magic function to make the
output of plotting commands displayed inline within the notebook.
.. code:: ipython3
%matplotlib inline
Our code starts with loading the required modules for the computation.
It is irrelevant whether we use the ``pireal`` or ``picomplex`` modules
here, because we do not calculate any matrix elements.
.. code:: ipython3
# Arrays
import numpy as np
# Plotting
import matplotlib.pyplot as plt
# pairinteraction :-)
from pairinteraction import pireal as pi
Keep in mind that pairinteraction is a **Rydberg** interaction
calculator. All techniques presented below work well for the
high-:math:`n` states of Rydberg atoms, but fail miserably for low-lying
states.
While the pairinteraction software normally uses the unit system
described in the
`introduction `__,
for the wavefunction it has proven advantageous to perform the
calculation in atomic units, i.e. all units below are atomic units (in
contrast to the presentation in Appendix A of the `pairinteraction
paper `__, where we use SI
units throughout).
Coulomb wavefunctions
---------------------
Using non-relativistic quantum defect theory it is possible to find
radial wavefunctions. Therefore the Schrödinger equation is solved for a
*given* energy eigenvalue, which results in modified radial
wavefunctions with an effective (non-integer) quantum number
:math:`n^*`. The derivation is discussed in more detail in M. J. Seaton,
`Reports on Progress in Physics 46, 167
(1983) `__ (around
equation 2.77).
.. math::
r\;\Psi^{\text{rad}}_{n^*l}(r) = \frac{1}{\sqrt{n^{*2} \Gamma(n^*+l+1) \Gamma(n^*-l)}}
W_{n^*,l+1/2}\left( \frac{2 r}{n^*} \right).
where :math:`\Gamma(z)` is the Gamma function and :math:`W_{k,m}(z)` is
the so-called Whittaker function. These wavefunctions are highly
accurate for large distances and have the correct binding energy. You
can access these wavefunctions in pairinteraction using the
``Whittaker`` class.
Numerov’s method and model potentials
-------------------------------------
Model potentials
~~~~~~~~~~~~~~~~
Another method to determine the radial wavefunctions is by solving the
Schrödinger equation numerically using an *effective* Coulomb potential.
The so-called model potential is usually decomposed into three
contributions
.. math::
V_{\text{mod}}(r) = V_{\text{C}}(r) + V_{\text{P}}(r) + V_{\text{s.o.}}(r).
The first term is the charge contribution resulting from the charged
core which is screened by the filled electron shells
.. math::
% (18a) in Dalgarno1994
V_{\mathrm{C}}(r) = - \frac{1 + (Z-1)\mathrm{e}^{-\alpha_1 r} - r(\alpha_3+\alpha_4 r)\mathrm{e}^{-\alpha_2 r}}{r},
with coefficients :math:`\alpha_{1,2,3,4}` depending on the atomic
species and the orbital angular momentum :math:`l`. These coefficients
are provided with pairinteraction in an embedded database. It is
possible to substitute these coefficients with your own at runtime.
For the core polarization, only the leading dipole term is considered,
which results in:
.. math::
V_{\mathrm{P}}(r) = -\frac{\alpha_d}{2 r^4} \left[ 1 - \mathrm{e}^{-(r/r_c)^6} \right].
Here, :math:`\alpha_d` is again the core dipole polarizability and
:math:`r_c` is the effective core size, obtained by comparing the
numerical solutions with the experimentally observed energy levels.
In addition to these two terms, we add an effective expression for the
spin-orbit interaction
.. math::
V_{\mathrm{s.o.}}(r > r_c) = \frac{2 \alpha}{r^3} \mathbf{l}\cdot\mathbf{s} = \frac{\alpha}{r^3} [j(j+1) - l(l+1) - s(s+1)].
where :math:`\alpha\approx1/137` is the fine-structure constant.
Numerov’s method
~~~~~~~~~~~~~~~~
Numerov’s method can be used to solve differential equations of the form
.. math::
\frac{d^2 y}{d x^2} + g(x) y(x) = 0.
The iterative solution can be achieved as follows:
.. math::
y_{n+1} \left( 1 - \frac{h^2}{12} g_{n+1} \right)
= \left( 2 + \frac{5 h^2}{6} g_{n} \right) y_{n}
- \left( 1 - \frac{h^2}{12} g_{n-1} \right) y_{n-1}
+ \mathcal{O}(h^6).
where :math:`y_{n}=y(x_{n})`, :math:`g_{n}=g(x_{n})`,
:math:`s_{n}=s(x_{n})`, and :math:`h=x_{n+1}-x_{n}`.
To save space it is useful to reduce the sampling towards large
distances from the core as the periods of osciallations become longer
and longer. To this end we introduce the square root scaling:
.. math::
x = \sqrt{r} \;,\quad X^{\text{rad}}_{nlj}(x) = x^{3/2} \Psi^{\text{rad}}_{nlj}(r).
This scaling keeps the number of grid points between nodes of the wave
function constant. In this scaling the :math:`g(r)` in Numerov’s method
is given by
.. math::
g(r) = \frac{(2 l + 1/2)(2 l + 3/2)}{r} + 8 r (V_{\text{mod}}(r) - E).
We integrate the equation from outside to inside. Since the wavefunction
has to decay to zero at infinity we choose :math:`y_0 = 0`. Now
depending on the number of nodes of wavefunction we decide whether to
set :math:`y_1 = \pm\epsilon`, where :math:`\epsilon` is a small number.
As the inner cutoff we choose an augmented version of the classical
turning point
.. math::
r_{\text{min}} = n^2 - n \sqrt{n^2 - (l-1)^2}.
Other schemes for terminating the integration exist which are based on
identifying nodes but for large :math:`n` the augmented classical
turning point works well and is easy to implement.
Quantum defects and model potential in pairinteraction
------------------------------------------------------
The model potentials and quantum defects are stored in a database which
is shipped together with pairinteraction. For performance reasons it is
baked into the binary so that it is available in memory at all times.
The internal database can also be swapped out with a user-provided one
at runtime.
To obtain the parameters from the database we use the ``QuantumDefect``
class.
.. code:: ipython3
qd = pi.QuantumDefect("Rb", 50, 0, 0.5)
The parameters of the model potentials can be accessed like member
variables of the ``qd`` object. They have mnemonic names to mirror their
meaning in the model potentials presented above.
.. code:: ipython3
print("Core polarizability: ac =",qd.ac)
print("Effective coulomb potential")
print(" Z =",qd.Z,"(core charge)")
print(" a1 =",qd.a1)
print(" a2 =",qd.a2)
print(" a3 =",qd.a3)
print(" a4 =",qd.a4)
print("Effective core radius: rc =",qd.rc)
.. parsed-literal::
Core polarizability: ac = 9.076
Effective coulomb potential
Z = 37 (core charge)
a1 = 3.69628474
a2 = 1.64915255
a3 = -9.86069196
a4 = 0.19579987
Effective core radius: rc = 1.66242117
The effective principal quantum number in quantum defect theory is
defined as series expansion
.. math::
n^* = n - \delta_{nlj}
\quad\text{with}\quad
\delta_{nlj} = \delta_0
+ \frac{\delta_2}{(n-\delta_0)^2}
+ \frac{\delta_4}{(n-\delta_0)^4}
+ \frac{\delta_6}{(n-\delta_0)^6} .
Similarly the effective energy eigenvalues are given as
.. math::
E_{nlj} = - \frac{1}{2 R_\infty} \frac{R^*}{n^{*2}} .
The parameters :math:`\delta_{0,2,4,6,\dots}` and :math:`R^*` are also
loaded from the database, but since they are only important for the
calculation of :math:`n^*` and :math:`E_{nlj}` there is no interface to
query their values. The quantities :math:`n^*` and :math:`E_{nlj}`,
however, can be queried.
.. code:: ipython3
print("Effective quantum number: n* =",qd.nstar)
print("State energy: E(n*) =",qd.energy)
.. parsed-literal::
Effective quantum number: n* = 46.868737950193115
State energy: E(n*) = -1497.6342916553963
Even though these parameters can be accessed like member variables, they
are read-only values and attempting to change them will result in an
error.
Radial wavefunction in pairinteraction
--------------------------------------
Two types of radial wavefunctions are available in pairinteraction:
1. The numerical wavefunctions computed using Numerov’s method and the
model potentials.
2. The analytical Coulomb wavefunctions based on the Whittaker
functions.
Those two schemes can accessed by classes with the names of their
underlying method, ``Numerov`` and ``Whittaker``. Merely instatiating an
object of either class only allocates memory but does not perform any
computations. To actually run the integration, you have to call the
``integrate()`` member function.
.. code:: ipython3
n = pi.Numerov(qd).integrate()
w = pi.Whittaker(qd).integrate()
The wavefunctions which we obtain from these methods are unscaled,
i.e. they have their original scaling as defined by the calculation. The
:math:`x`-axes are square root scaled. The result returned by the
Coulomb wavefunction method is :math:`r \;\Psi^{\text{rad}}(r)`. The
result calculated by Numerov’s method is :math:`X(x)`, as defined above,
and must be multiplied by :math:`\sqrt{x}` to get
:math:`r \;\Psi^{\text{rad}}(r)`.
.. code:: ipython3
plt.xlabel("$r$ ($a_0$)")
plt.ylabel("$r^2|\Psi(r)|^2$ [a.u.]")
plt.plot(n[:,0]**2,np.abs(np.sqrt(n[:,0])*n[:,1])**2,'-',label="numeric WF")
plt.plot(w[:,0]**2,np.abs(w[:,1])**2,'--',label="Coulomb WF")
plt.legend();
.. image:: wavefunctions_files/wavefunctions_15_0.png