Complex Absorbing Potentials¶
This example illustrates the use of a complex absorbing potential (CAP) with the qdyn library for the calculation of resonance energies and lifetimes.
The effective Hamiltonian of a simple model system, which obeys a shape resonance, is diagonalised using two different CAPs: a quadratic one (following the spirit of Ref. [3] and a more stable transmission free complex absorbing potential (TFCAP) introduced in Ref. [4].
Additionally you can see in the second part of this notebook how to modify config files using the QDYN Python library and how to run several instances of a Fortran program in parallel using the subprocess module.
[1]:
import matplotlib.pyplot as plt
import numpy as np
Background¶
Here we use a CAP to determine the lifetime of resonance states by solving the stationary Schrödinger equation. In general, a CAP can be used to avoid reflections of the wavefunction at the boundaries of the grid, which can also be useful in time-dependent calculations.
In this example we consider the scattering of a particle with unit mass by a repulsive potential, described by the model system introduced in Ref. [5]. It consists of a single potential curve \(V\) on a radial coordinate \(r\),
\begin{equation} V(r) = 7.5 \, r^2 \, e^{-r} \end{equation}
This potential allows for a shape resonance at \(E_{\mathrm{res}}\approx3.4\) with \(\Gamma\approx0.026\) and is shown in the following figure.
[2]:
r = np.linspace(0, 25, 200)
V = 7.5 * r**2 * np.exp(-r)
fig, ax = plt.subplots()
ax.hlines(3.426, 0, 1.2)
ax.plot(r, V, '-')
ax.set_xlim(0, 25)
ax.set_ylim(0, 5)
ax.set_xlabel(r'$r$')
ax.set_ylabel(r'$V(r)$')
ax.set_title(r'Model system')
plt.show()
The lifetime of the shape resonance can be obtained by adding a CAP to the end of the grid, causing a non-hermitian effective Hamiltonian
\begin{equation} H' = T + V(r) - i \eta W(r) \, \Theta(r - r_c) \;. \end{equation}
The parameter \(\eta\) controls the strength of the CAP and the cut-off parameter \(r_c\) in the step function \(\Theta\) allows to ensure, that the CAP does not disturb the coordinate range, where the physical potential is important.
The eigenvalues \(\mathcal{E}_i\) of \(H'\) are complex. Their real part corresponds to the physical eigenenergies \(E_i\), whereas for resonances the imaginary part is connected to their width \(\Gamma\) as
\begin{equation} \mathcal{E}_{\mathrm{res}} = E_{\mathrm{res}} - i \frac{\Gamma}{2} \;. \end{equation}
In general, the eigenvalues \(\mathcal{E}_i\) depend on \(\eta\) and \(r_c\). However, a resonance eigenvalue converges towards the physical resonance energy and width for a certain choice of the CAP parameters.
For a more detailed introduction to CAPs and their numerical properties see Ref. [3] or chapter 3 and B of Ref. [6].
Complex Eigenenergies¶
First we need to compile the program
[3]:
%%bash
make
make: 'complex_diag' is up to date.
To obtain the complex eigenvalues of the effective Hamiltonian run the program with the provided runfolders.
We start with the quadratic CAP, which has the form \(W(r) = (r-r_c)^2\).
[4]:
! ./complex_diag run_quadratic
*** Read config file run_quadratic/config ***
*** Done reading config file ***
*** Initializing system ***
*** Initializing grid ***
Initializing grid for label `` as 1D cartesian grid without mapping
*** Initializing pulses ***
*** Initializing dynamical generator ***
*** Done with initialization ***
eigensystem for non-coupled system
surf 1: eigenvalues are shifted down by E_0 = 0.000000000000E+00
The transmission free CAP is defined as \(W(r) = \left(\frac{1}{1-x^2} + \frac{1}{1+x^2} - 2\right)\) with \(x = \frac{r-r_c}{r_{\mathrm{max}}-r}\).
Using the TFCAP, in principle the correct parameters for the calculation of a resonance energy can be estimated by
\begin{equation} \eta = \frac{4 E_{\mathrm{min}}}{C^2} \;, \end{equation}
with \(C = 2.62206\) and \(E_{\mathrm{min}}\) beeing the lowest scattering energy of interest. Similar, the absorption range \(D = r_{\mathrm{max}}-r_c\) is given by the de Broglie wavelength at this energy, i.e.
\begin{equation} D = \frac{2\pi}{\sqrt{2\mu E_{\mathrm{min}}}} \;. \end{equation}
In practice you can use these values as a starting point, but you probably have to perform some manual parameter optimisation, too.
[5]:
! ./complex_diag run_tfcap
*** Read config file run_tfcap/config ***
*** Done reading config file ***
*** Initializing system ***
*** Initializing grid ***
Initializing grid for label `` as 1D cartesian grid without mapping
*** Initializing pulses ***
*** Initializing dynamical generator ***
*** Done with initialization ***
eigensystem for non-coupled system
surf 1: eigenvalues are shifted down by E_0 = 0.000000000000E+00
Next, we read in the obtained eigenvalues and plot them in the complex plane.
[6]:
E_quad = np.loadtxt('run_quadratic/eigenvals_rc005_eta000030348E-6.dat',
usecols=(1,2))
E_tfcap = np.loadtxt('run_tfcap/eigenvals_rc005_eta001890445E-6.dat',
usecols=(1,2))
[7]:
fig, [ax0, ax1] = plt.subplots(1, 2, figsize=(12,4))
ax0.plot(E_quad[:,0], E_quad[:,1], 'p')
ax0.set_xlim(0,10)
ax0.set_ylim(-6,0)
ax0.set_xlabel(r'$\mathrm{Re}(\mathcal{E})$')
ax0.set_ylabel(r'$\mathrm{Im}(\mathcal{E})$')
ax0.set_title('Quadratic CAP')
ax1.plot(E_tfcap[:,0], E_tfcap[:,1], 'p')
ax1.set_xlim(0,10)
ax1.set_ylim(-6,0)
ax1.set_xlabel(r'$\mathrm{Re}(\mathcal{E})$')
ax1.set_ylabel(r'$\mathrm{Im}(\mathcal{E})$')
ax1.set_title('TFCAP')
plt.show()
The CAP rotates the eigenvalues into the lower complex plane. The distribution of the eigenvalues depends on the functional form of the CAP and the parameters \(\eta\) and \(r_c\).
If the CAP parameters are tuned correctly, resonance states reveal themselves as isolated eigenvalues between the spectral string and the real axis.
You can easily adjust the CAP strength and the absorption range in the config files and observe their effect on the spectrum.
Converging Resonance States¶
The complex eigenvalues of the effective Hamiltonian depend on the parameters of the CAP. However, a resonance eigenvalue converges towards the physical resonance energy and width for a certain choice of parameters.
Usually one fixes the cut-off parameter \(r_c\) and varies the CAP strength \(\eta\). By this, you obtain \(\eta\)-trajectories of the eigenvalues in the complex plane. Eigenvalues belonging to resonance states then show a point of stability. This corresponds to the optimal value of \(\eta\) (which in principle is different for every resonance energy).
In the following, we calculate \(\eta\)-trajectories by running the complex diagonalisation multiple times with different CAP strengths. For this purpose, we create additional temporary runfolder and use the QDYN Python library to generate the config files.
If the QDYN Python library is not installed, you can do so with pip install qdyn
[8]:
from parallel_run import parallel_diag
We start with \(\eta=0\) and increase the CAP strength exponentially.
[9]:
factor = 1e-5
slope = 1.10
n = np.arange(0, 150, 1)
etas = factor * (slope**n - 1.) / (slope - 1.)
[10]:
runfolder = 'run_quadratic' # or 'run_tfcap'
parallel_diag(etas, runfolder, 10)
Running 150 processes using 10 cpus...
Done.
After the calculation is complete, we extract the lowest resonance from the data.
We do this by looking at the probability amplitude of the wavefunctions at small distances, since resonance wavefunctions have a larger probability amplitude in the region of the potential well.
After the diagonalisation, the program complex_diag integrates the eigenfunctions over the interval \([0, 5]\). The resulting value is also written to the file that contains the eigenenergies. In order to find the lowest resonance, we just need to look for the maximum of this inner population.
[11]:
Eres = np.zeros((len(etas), 2))
for i, eta in enumerate(etas):
filename = runfolder+'/'+'eigenvals_rc005_eta{:09}E-6.dat'.format(int(round(eta*1e+6)))
ReE, ImE, pop = np.loadtxt(filename, usecols=(1,2,3), unpack=True)
ires = pop.argmax()
Eres[i] = (ReE[ires], ImE[ires])
Now, we can plot the \(\eta\)-trajectory of the lowest shape resonance in the complex plane.
[12]:
fig, ax = plt.subplots()
ax.plot(Eres[:,0], Eres[:,1], 'o')
ax.set_xlabel(r'$\mathrm{Re}(\mathcal{E})$')
ax.set_ylabel(r'$\mathrm{Im}(\mathcal{E})$')
ax.set_title(r'$\eta$-trajectory of lowest resonance')
plt.show()
One can see, how the eigenvalue starts on the real axis for \(\eta=0\) and moves into the complex plane for increasing CAP strength until it reaches a stable point. This corresponds to the optimal CAP strength. When \(\eta\) further increases, the eigenvalue starts to move again.
The found resonance energy and width do not match the known values exactly. This is probably related to the basis set used.