Hex  1.0 Hydrogen-electron collision solver
Hex Documentation
Date
1. 2. 2014

# Hex-ecs

Hex-ECS computes partial T-matrices for elastic, excitation and ionization electron-hydrogen scattering. Together with the interface program Hex-DB it offers the possibility of generating many scatterign quantities like differential or integral cross sections.

## Language and libraries

Hex is written in C++11 to make use of comfort of the modern C++ extensions, so one may need a newer compiler. Tested compilers are:

• GCC 4.8.1
• Intel C++ Composer XE 14.0 SP1 Update 1 (tested with GCC 4.8.1 headers)

Both worked with the same Makefile, just by setting the variable CPP to "g++" or "icpc". The program also uses following external packages (tested versions are given in parentheses):

• GNU Scientific Library (1.16): for Wigner coupling coefficients and some other special functions.
• SuiteSparse/UMFPACK (4.2.1/5.6.2): for sparse matrix manipulation and for a direct sparse system solver.
• HDF5 (1.8.11): for standardized binary output of intermediate results (used to speed up further simillar calculations).
• FFTW (3.3.3): for fast evaluation of Chebyshev expansions by Fast Fourier Transform.
• OpenMP (2.1): for parallelization at single machine.
• MPI (OpenMPI 1.6): for parallelization at cluster.

The next libraries are optional:

• OpenBLAS (0.2.8): Free BLAS implementation that can be compiled for a specific CPU. OpenBLAS is able to run in parallel using pthreads or OpenMP. OpenBLAS is optional because SuiteSparse can be configured to use a different BLAS implementation.
• png++ (0.2.5, requires libpng-1.5.x or older): PNG read/write interface, for debugging purposes. The code that references png++ can be excluded from compilation by the option -DNO_PNG. For convenience, this library is directly included in the source code of Hex-ecs.

All listed libraries are open-source and easily obtainable on the internet and in the repositories of some Linux distributions (openSUSE 12.3 at least).

Equations in this documentation use MathJax, which should work in every up-to-date JavaScript-enabled web browser. Tested browsers are:

• Mozilla Firefox 24.0
• Konqueror 4.10.5

## Theory

Unknown scattering state, eigenfunction of the full system hamiltonian, is split into two parts, asymptotic “incoming particle” state and the scattered part, which is a solution of the driven Schrödinger equation,

 $\left(E - \hat{H}\right) \Psi_{\mathrm{sc}} = \hat{H}_{\mathrm{int}} \Psi_{\mathrm{inc}} \ .$ (1)

The state $$\Psi_{\mathrm{inc}}$$ is a product of initial atomic state and a projectile plane wave,

 $\Psi_{\mathrm{inc}}(\mathbf{r}_1, \mathbf{r}_2) = \frac{1}{k_i r_1 r_2} P_{n_i l_i}(r_1) Y_{l_i m_i}(\mathbf{\hat{r}}) \cdot (2\pi)^{-3/2} \sum_{lm} \sqrt{4\pi} \mathrm{i}^l \hat{j}_l(k_i r_2) Y_{lm}(\mathbf{\hat{k}}_i) Y_{lm}^\ast(\mathbf{\hat{r}}_2) \ ,$ (2)

and anti/symmetrized with respect to electron exchange, according to total spin,

 $\Psi^S_{\mathrm{inc}}(\mathbf{r}_1, \mathbf{r}_2) = \Psi_{\mathrm{inc}}(\mathbf{r}_1, \mathbf{r}_2) + (-1)^S \Psi_{\mathrm{inc}}(\mathbf{r}_2, \mathbf{r}_1) \ .$ (3)

Hamiltonian $$\hat{H}$$ is a sum of free hamiltonian and interaction hamiltonian. Free hamiltonian is the hamilton operator for electron and atom separated far away,

 $\hat{H}_{\mathrm{free}} = -\frac{\nabla_1^2}{2} - \frac{\nabla_2^2}{2} -\frac{1}{r_1} \ ,$ (4)

whereas the interaction hamiltonian contains the interaction between projectile and the atomic constituents, that is

 $\hat{H}_{\mathrm{int}} = \frac{1}{r_{12}} - \frac{1}{r_2}$ (5a)

for direct case and

 $\hat{H}_{\mathrm{int}} = \frac{1}{r_{12}} - \frac{1}{r_1}$ (5b)

for exchange (anti/symmetrized) case. Having the solution of equation (1), we can extract the scattering amplitude, which is a matrix element of $$\hat{H}_{\mathrm{int}}$$ between the solution $$\Psi = \Psi_{\mathrm{sc}} + \Psi_{\mathrm{inc}}$$ and some asymptotic state $$\Psi_{\mathrm{out}}$$ that is experimentally inquired. Conventionally, $$\Psi_{\mathrm{out}}$$ is also a product of (now final) atomic state, possibly different from the original, with an outgoing plane wave of the projectile. Any exchange effect are already included in the original anti/symmetrization, so they need not be considered now. The projection

 $f = -4\pi^2 \left<\Psi_{\mathrm{out}}\right|\hat{H}_{\mathrm{int}}\left|\Psi\right>$ (6)

can be also written as

 $f = -4\pi^2 \left<\Psi_{\mathrm{out}}\right|E - \hat{H}_{\mathrm{free}}\left|\Psi_{\mathrm{sc}}\right> \ ,$ (7)

where the following was used: $$\hat{H}_{\mathrm{int}} = \hat{H} - \hat{H}_{\mathrm{free}}$$, $$\Psi = \Psi_{\mathrm{sc}} + \Psi_{\mathrm{inc}}$$ and $$\hat{H}\Psi = E\Psi$$.

## Method

All computations are done in time-independent way, and the exterior complex scaling is used instead of boundary condition fitting. Radial part of sought wave-functions is expanded in a B-spline basis $$\left\{B_i\right\}_{i = 0}^{\mathrm{Nspline}-1}$$ of a given order,

 $\Psi_{\mathrm{sc}}^{LM}(\mathbf{r}_1, \mathbf{r}_2) = \sum_{l_1 l_2} \psi_{l_1 l_2}^{LM}(\mathrm{r}_1, \mathrm{r}_2) \mathcal{Y}_{l_1 l_2}^{LM}(\mathbf{\hat{r}}_1, \mathbf{\hat{r}}_2) \ ,$ (8)
 $\psi^{LM}_{\mathrm{sc},l_1 l_2}(\mathrm{r}_1, \mathrm{r}_2) = \frac{1}{r_1 r_2} \sum_{ij} \psi_{l_1 l_2,ij}^{LM} B_i(r_1) B_j(r_2) \ ,$ (9)

and when projecting the equation (1) on a bipolar spherical function $$\left<\mathcal{Y}_{l_1 l_2}^{LM}\right|$$ to get rid of angular dependence, and on a pair of B-splines to get rid of radial dependence and keep only matrix elements, one arrives at a matrix equation for components of $$\psi_{l_1 l_2}^{LM}(\mathrm{r}_1,\mathrm{r}_2)$$,

 $\left[\left(ES_{ik}S_{jl} - \frac{1}{2} D_{ik}S_{jl} - \frac{1}{2} S_{ik} D_{jl} - \frac{1}{2}l_1(l_1+1) M_{ik}^{(-2)}S_{jl} - \frac{1}{2}l_2(l_2+1) S_{ik}M_{jl}^{(-2)} + M_{ik}^{(-1)}S_{jl} + S_{ik}M_{jl}^{(-1)}\right) \delta_{l_1}^{l_1'}\delta_{l_2}^{l_2'} - \sum_{\lambda} f_{l_1l_2l_1'l_2';L}^\lambda R_{ijkl}^\lambda \right] \psi_{l_1' l_2',kl}^{LMS} = \chi_{l_1l_2,kl}^{LMS} \ ,$ (10)

or symbolically

 $\left[ \mathsf{Id}_1 \otimes \mathsf{Id}_2 \otimes \left(E\mathsf{S}\otimes\mathsf{S} - \frac{1}{2} \mathsf{D}\otimes\mathsf{S} - \frac{1}{2} \mathsf{S}\otimes\mathsf{D} - \frac{1}{2} l_1 (l_1+1) \mathsf{M}^{(-2)}\otimes\mathsf{S} - \frac{1}{2} l_2 (l_2+1) \mathsf{S}\otimes\mathsf{M}^{(-2)} + \mathsf{M}^{(-1)} \otimes \mathsf{S} + \mathsf{S} \otimes \mathsf{M}^{(-1)}\right) + \sum_{\lambda} \mathsf{f}_L^\lambda \otimes \mathsf{R}^\lambda \right] \mathsf{\psi}^{LMS} = \mathsf{\chi}^{LMS} \ .$ (10*)

Symbol $$\otimes$$ stands for Kronecker product (“flattened tensor product”) and matrices have following meanings:

• Matrix $$\mathsf{Id}_1$$ is identity of rank equal to maximal allowed angular momentum $$l_1$$. Analogically for $$\mathsf{Id}_2$$ .
• Matrix $$\mathsf{D}$$ is matrix of derivative overlaps of B-splines. It can be shown for B-splines basis which is zero at boundaries that in such case

$\left<B_i\right| \left(-\frac{\mathrm{d}^2}{\mathrm{d} x^2}\right) \left|B_k\right> = + \int_a^b \frac{\mathrm{d}B_i}{\mathrm{d}x} \frac{\mathrm{d}B_k}{\mathrm{d}x} \equiv + D_{ik} \ .$

• Matrix $$\mathsf{S}$$ is just standard overlap matrix of the B-spline basis,

$S_{ik} = \left<B_i\right|\left.\!B_k\right> \ .$

• Matrix $$\mathsf{M}^{(\alpha)}$$ is matrix element of power of coordinate (also called integral moment in the source code).
• Matrix $$\mathsf{R}^{\lambda}$$ is matrix of four-B-spline multipole integrals for multipole $$\lambda$$,

$R_{ijkl}^\lambda = \int_a^b\int_a^b B_i(r_1) B_j(r_2) \frac{r_<^\lambda}{r_>^{\lambda + 1}} B_k(r_1) B_l(r_2) dr_1 dr_2 \ ,$

flattened so that “i” and “j” form one multi-index [ij] and the other indices the multi-index [kl].
• Matrix $$f_L^\lambda$$ is angular part of reduced matrix element,

$\left<l_1 l_2 \right|\!| \frac{1}{r_{12}} |\!\left|l_1' l_2'\right>_L = \sum_{\lambda} f_{l_1 l_2 l_1' l_2'; L}^\lambda \frac{r_<^\lambda}{r_>^{\lambda + 1}} \ .$

flattened so that “l₁” and “l₂” form one multiindex [l₁l₂] and the other two indices the second.

Finally, the symbol $$\chi^{LMS}$$ stands for projection of the right hand side, which is

 $\chi_{l_1 l_2, ij}^{LMS} = \frac{1}{k_f} \sum_\ell \mathrm{i}^\ell \sqrt{2\pi(2\ell+1)} \left\{ \left(\sum_\lambda f_{l_1 l_2 l_i \ell; L}^\lambda R_{ijkl}^\lambda - \delta_{l_1}^{l_i} \delta_{l_2}^{\ell} S_{ik}M_{jl}^{(-1)}\right) \left[P_{n_il_i}(r_1)\right]_k \left[\hat{j}_\ell(k_i r_2)\right]_l + (-1)^{S+\Pi} \left(\sum_\lambda f_{l_1 l_2 \ell l_i; L}^\lambda R_{ijkl}^\lambda - \delta_{l_1}^{\ell} \delta_{l_2}^{l_i} M_{ik}^{(-1)} S_{jl}\right) \left[\hat{j}_\ell(k_i r_1)\right]_k \left[P_{n_il_i}(r_2)\right]_l \right\}$ (11)

or symbolically

 $\chi^{LMS} = \frac{1}{k_f} \sum_\ell \mathrm{i}^\ell \sqrt{2\pi(2\ell+1)} C_{l_i m_i \ell 0}^{L M} \left\{ \left(\sum_\lambda \mathsf{f}_L^\lambda \otimes \mathsf{R}^\lambda - \mathsf{Id}_1 \otimes \mathsf{Id}_2 \otimes \mathsf{S} \otimes \mathsf{M}^{(-1)}\right) \cdot \Delta_{l_1}^{l_i} \otimes \Delta_{l_2}^{\ell} \otimes \mathsf{P}_{n_i l_i} \otimes \mathsf{j}_{\ell,k_i} + (-1)^{S+\Pi} \left(\sum_\lambda \mathsf{f}_L^\lambda \otimes \mathsf{R}^\lambda - \mathsf{Id}_1 \otimes \mathsf{Id}_2 \otimes \mathsf{M}^{(-1)} \otimes \mathsf{S}\right) \cdot \Delta_{l_1}^{\ell} \otimes \Delta_{l_2}^{l_i} \otimes \mathsf{j}_{\ell,k_i} \otimes \mathsf{P}_{n_i l_i} \right\} \ .$ (11*)

Here, the one-dimensional (!) vectors $$\Delta$$ are zero vectors with only one element equal to one at position $$l_1 = l_i$$ (etc.). P- and j- vectors are components of respective function in chosen B-spline basis.

Factor C in the expression (11*) is the Clebsch-Gordan coefficient and the zero projection of $$\ell$$-momentum reflect the deliberate choice of scattering axis along the projectile moemntum, so that the angular momentum projection is zero.

## ECS restrictions on potential

Exterior complex scaling of right hand side poses a serious problem for typical (not exponentially decreasing) potentials. One of the factors in the right hand side is the Riccati-Bessel functions $$\hat{j}$$, which exponentially diverges under ECS transformation. To avoid this, potential $$\hat{H}_{\mathrm{int}}$$ is artificially truncated at (or before) the turning point $$R_0$$, which has several consequences for the numerical construction:

• In equations (6)-(7), the radial integration is done only up to $$r_1,r_2 = R_0$$ and not further.
• Matrices $$\mathsf{S}, \mathsf{M}^{(-1)}, \mathsf{R}^\lambda$$ in (11), (11*) are to be computed, again, for $$r_1,r_2 \le R_0$$. These are referenced as “truncated overlap matrices” in the source code.

## Cross section

As was said above, the scattering amplitude is

 $f = -4\pi^2 \left<\Psi_{\mathrm{out}}\right|E - \hat{H}_{\mathrm{free}}\left|\Psi_{\mathrm{sc}}\right>_{R_0}\ ,$ (12)

where the subscript $$R_0$$ means, that the radial integration is done only for radii less than $$R_0$$ (because the original matrix, $$\hat{H}_{\mathrm{int}}$$, has to be truncated at such distance to avoid far-region divergence. The outgoing – detected – wavefunction $$\Psi_{\mathrm{out}}$$ has form simillar to (2), with replaced initial to final quantum numbers. Substituting such expansion into the equation (12) and once again using zero boundary trait of chosen B-spline basis when doing per parts integration one easily arrives at the formula for cross section

 $\sigma = \frac{4}{k_i k_f} \sum_{\ell L L'} C_{l_f m_f \ell 0}^{L m_f} C_{l_f m_f \ell 0}^{L' m_f} \left| \psi_{l_f \ell, ij}^{LMS} W[P]_i S[j]_j + \psi_{l_f \ell, ij}^{L'MS} S[P]_i W[j]_j \right|^2 \ ,$ (13)

where $$W[P]_i$$ and $$W[j]_j$$ stand for wronskian (evaluated at $$R_0 - \varepsilon$$) of the i-th B-spline and the (final) hydrogenic or Riccatti-Bessel function, and $$S[P]_i$$ and $$S[j]_i$$ stand for (truncated) overlap integrals of the i-th B-spline and the (final) hydrogenic or Riccatti-Bessel function.

## Implementation in the code

The code runs along the following outline:

• The program is initialized. Parameters from the command line are stored for later use in the class CommandLine.
• The optional multi-process parallel environment (MPI) is set up. All parallel communication is mediated by the class Parallel. The parallelization by MPI is used whenever the switch –mpi is present on the command line. The number of processes (= "communicator rank") is given as an argument to the MPI launcher. E.g.
      mpirun -np 4 hex-ecs --mpi
• Note that the parallelization by OpenMP is used always and does not use any classes, nor does it depend on MPI in any way. The number of OpenMP threads is completely independent on the communicator size and can be specified by the environment variable OMP_NUM_THREADS. The number of threads on every machine used in the MPI network can be different.
• From here on, all the steps are executed by every machine in the MPI scope.
• The input file specified on the command line (–input/-i) is found and read. If the name was not given on command line, the default "hex.inp" is used. If the file doesn't exist, the program aborts. The values from the input file are contained in the object InputFile.
• A B-spline basis is created according to the input (class Bspline).
• A list of available angular states is assembled. The number of these states is restricted by $$L$$, $$\Pi$$ and $$n_L$$.
• A preconditioner is selected in accord with the user's choice. It is stored only in the form of a pointer to the base class PreconditionerBase, whose methods are reimplemented in all derived classes (specific preconditioners). The constructor is called, which will in all cases trigger also the computation of the radial integrals that are needed by all preconditioners. The radial integrals are managed by the class RadialIntegrals.
• For all impact energies:
• The preconditioner for the set of equations is updated for this energy (method "update").
• For all initial states:
• Check that at least some allowed angular states contribute to this combination of $$L$$, $$\Pi$$, $$n_L$$ and $$l_i$$. If not, skip this inital state.
• Check that the solution for this initial state and this impact energy hasn't been already computed. If it has been, skip this initial state.
• The B-spline expansion of the right hand side is computed (method "rhs" of the chosen preconditioner class).
• Run the preconditioned conjugate gradients solver. If there already is a solution for the previous energy, use it as an initial guess.
• Save the solution to disk.
• For all initial and final states:
• Compute T-matrices for all available partial waves.
• Save the T-matrices to SQL files.
• Evaluate and print the integral cross sections.