% Manuscript date: 14 August 2000
\documentstyle[preprint,pra,aps]{revtex}
\tighten
\begin{document}
\title{PRACTICAL ATOMIC PHYSICS}
\author{Yong-Ki Kim}
\address{Atomic Physics Division \\ National Institute of Standards and
Technology \\ Gaithersburg, Maryland 20899-8423 U.S.A.}
\date{\today}
\maketitle
\begin{abstract}
This lecture note describes how to calculate (a) relativistic wave
functions for neutral atoms and atomic ions, (b) atomic energy levels,
(c) bound-bound transition probabilities, (c) bound-bound plane-wave
Born excitation cross sections by electron impact, (d) bound-bound
distorted wave Born excitation cross sections by electron impact, (e)
bound-free photoionization cross sections, and (f) bound-free ionization
cross sections by electron impact. These atomic data are obtained from
the multiconfiguration Dirac-Fock (MCDF) wave functions described in
(a).
\end{abstract}
\newpage
\section{Introduction}
\label{sec1}
Many decades after the Schr\"{o}dinger equation was discovered in 1925,
we still have difficulty in calculating ``exact'' wave functions for
many-electron atoms. For all practical purposes, ``exact'' or ``almost
exact'' wave functions are available only for hydrogen-like and
helium-like systems. ``Almost exact'' solutions to the electron-impact
ionization of the hydrogen atom have been available only in the last few
years \cite{bray,mccurdy,pindzola}.
However, atomic physics data, particularly collision data, are needed in
solving many problems in other branches of physics, chemistry, and
engineering. These ``practical'' applications require ``practical''
solutions, majority of which involve neutral atoms and ions with far
more than two electrons.
The objective of this series of lectures is to describe methods to
obtain such ``practical'' solutions using ``practical'' wave functions.
Relativistic wave functions are needed to handle heavy atoms (roughly Kr
and heavier), particularly to handle intermediate coupling in heavy
atoms, highly excited states of light atoms, and highly charged ions.
These lectures will be based on the use of a set of Fortran codes
developed from the early 1970s primarily by Jean-Paul Desclaux of
Grenoble, France \cite{descl1} with help from Paul Indelicato of Paris
and myself at various stages of the code development and usage.
Theoretical details will be kept to the minimum.
The codes described in this note have been frozen in their contents
around 1992. This is why these codes are labeled DF92, MJ92, PH92, and
DW92. The backbone of the code package is the Dirac-Fock wave function
code, DF92. Other codes use the wave functions generated by DF92. The
functions of these codes are:
\newline \indent DF92: calculate single- and multi-configuration
Dirac-Fock wave functions;
\newline \indent MJ92: integrate angular functions and produce
appropriate Slater coefficients;
\newline \indent PH92: calculate bound-bound and bound-free
photoexcitation and photoionization cross sections; bound-bound
Plane-wave Born excitation cross sections by electron impact;
\newline \indent DW92: calculate bound-bound and bound-free
distorted-wave Born excitation and ionization cross sections by electron
impact.
In addition, a theoretical model called the BEB (binary-encounter-Bethe)
model has been developed \cite{kim94} to use simple data from DF92 to
calculate total cross sections for the ionization of atoms and ions by
electron impact. The BEB model is simple, analytic, and requires no
special computer codes. Spreadsheet programs on a PC, such as Lotus123
and MS Excel is sufficient.
It is customary to use atomic units in atomic and molecular
calculations. In atomic units, the electron rest mass $m$, Planck's
constant divided by 2$\pi = \hbar$, the electronic charge $e$, and the
Bohr radius $a_0$ are set to unity. All the physical quantities
associated with the ground state of the hydrogen atom, except for the
total energy, become unity in atomic units. The total energy of the
hydrogen atom is $-$0.5 a.u., i.e., 1 a.u. of energy is 27.2 eV. The
atomic unit of energy is called hartree, while the total energy of the
hydrogen atom with the {\it infinite mass} nucleus is called a rydberg
(= 0.5 hartree). The speed of light $c$ in a.u. becomes 1/$\alpha
\approx 137$, where $\alpha=e^2/\hbar c$ is the fine-structure constant.
Atomic units are used throughout this note. However, symbols
representing appropriate dimensions will be retained in expressions for
collision cross sections because cross sections are expressed in many
different units, viz., cm$^2$, \AA$^2$, barns, megabarns, $\pi a_0^2$,
\AA$^2$/eV, $\pi a_0^2/$hartree-steradian, etc. Some fundamental
constants associated with atomic units are listed in Table
\ref{tbl:table1}.
The audience of this lecture is assumed to be familiar with the basic
nonrelativistic quantum theory for atomic structure and collision. A
brief summary of the required atomic structure theory is described in
Sec. \ref{sec2}, collision theory in Sec. \ref{sec3}, the BEB model in
Sec. \ref{sec4}, the DF92 code in Sec. \ref{sec5}, MJ92 in Sec.
\ref{sec6}, PH92 in Sec. \ref{sec7}, and DW92 in Sec. \ref{sec8}. Input
data formats are described Sec. \ref{sec9}.
\section{Outline of Atomic Structure Theory}
\label{sec2}
All computer codes described in this lecture note use relativistic
Hartree-Fock (often called Dirac-Fock) wave functions. Relativistic
theories for both wave functions and collision cross sections have their
origin in nonrelativistic theories. The fundamental difference is that
nonrelativistic theories use LS coupling while relativistic theories use
jj coupling. Since the users of this note are much more familiar with
nonrelativistic notations than relativistic ones, I shall outline
underlying theories for atomic structure and collision in
nonrelativistic notations first, and then point out changes necessary to
switch to relativistic notations.
\subsection{Hartree-Fock Wave Functions}
\label{sec21}
For many-electron atoms, the only practical wave functions suitable for
various applications are Hartree-Fock (HF) wave functions. Hartree-Fock
wave functions are the best wave functions based on the independent
particle model, i.e., there is a wave function for each electron, and
the one-electron function does not contain any variable referring to
other electrons. These wave functions for individual electrons are
called {\it orbitals}. Each orbital for a bound electron in an atom
carries four quantum numbers, the principal quantum number $n$, the
angular momentum $l$, the projection $m_l$ of the angular momentum $l$,
and the projection $m_s$ of the spin angular momentum $s$. The spin
angular momentum for all individual electrons is 1/2, and hence not
explicitly included in the quantum number set.
For a many-electron atom, the one-electron quantum numbers $l,\ m_l$,
$s=$ 1/2, and $m_s$ are combined to determine the total orbital angular
momentum $\bbox{L}$, the total spin $\bbox{S}$, the total angular
momentum $\bbox{J}$, and matching projections $M_L,\ M_S$, and $M_J$.
A one-electron orbital used in the Hartree-Fock wave function has the
form:
\begin{equation}
\psi_{nlm_l}(r,\theta,\varphi) = r^{-
1}P_{nl}(r)Y_{lm_l}(\theta,\varphi),
\label{eq:hforbital}
\end{equation}
\noindent where $r^{-1}P_{nl}(r)$ is the one-electron radial wave
function, and $Y_{lm_l}(\theta,\varphi)$ are the spherical harmonics.
These one-electron functions satisfy the usual orthonormality
conditions.
Since the Schr\"{o}dinger equation does not explicitly contain
spin-dependent operators, the spin eigenfunctions are often omitted in
the definition of the one-electron orbital, Eq. (\ref{eq:hforbital}).
When necessary, the two-component spinors are used, $\hat{\alpha}$ for
spin up and $\hat{\beta}$ for spin down:
\begin{equation}
\hat{\alpha} = \left( \begin{array} {c} 1 \\ 0 \end{array} \right),
\hspace{1in} \hat{\beta} = \left( \begin{array} {c} 0 \\ 1 \end{array}
\right).
\label{eq:spinor}
\end{equation}
\noindent When one of these spinors are multiplied to Eq.
(\ref{eq:hforbital}), the resulting one-electron function with spin
projection is called a spin-orbital.
The total wave function $\Psi_{LSM_LM_S}$ is built by constructing one
or more Slater determinants from the one-electron orbitals. The Slater
determinants are weighted by Clebsch-Gordan coefficients to represent
appropriate $L,\ S,\ M_L$, and $M_S$. When the resulting $L$ and $S$
are unique, a Hartree-Fock wave function consists of one Slater
determinant. Closed-shell atoms, alkali atoms, and halogen atoms are
examples of the single-determinant case.
A simple list of occupied one-electron orbitals, such as 1s$^2$ for the
helium atom, is called a configuration. For example, 1s$^2$ 2s$^2$
2p$^2$ is the configuration for the ground state of the carbon atom. In
this case however, one must also specify $L$ and $S$ because two 2p
electrons can be coupled to have different $L$ and $S$: $^3$P, $^1$D,
and $^1$S. The superscript on the left is the spin multiplicity,
$2S+1$, and the capital letters denote $L$. The letters A, B, C, D, E,
and J and their lower case letters are not used for the orbital angular
momentum, though S/s(=0) and L/l(=8) are used. When coupling more than
two electrons, e.g., 3d$^3$, $\bbox{L}$ and $\bbox{S}$ are not
sufficient to give unique total angular quantum numbers. For such
cases, a seniority number must be assigned. For the cases of our
interest, however, seniority can be ignored. Instead, we must specify
whether an atomic level is the second, third, or fourth level with the
same $L,\ S,\ J$, and parity, usually starting from the lowest level in
energy.
A total wave function consists of Slater determinants weighted by
appropriate Clebsch-Gordan coefficients. Each Slater determinant is
normalized to unity, and the total wave function is also normalized to
unity. This normalization process also alters the coefficients in front
of each Slater determinant.
The Hartree-Fock equation is obtained by applying the variational
principle to the expectation value for the total energy. The expression
for the total energy can be reduced to a combination of one-electron
integrals involving a pair of $P_{nl}$'s and two-electron integrals
involving two to four different $P_{nl}$'s.
The spherical harmonics are analytically integrated, and appear as a
numerical constant in front of each one-electron and two-electron radial
integral. This constant, which includes the Clebsch-Gordan coefficient
to represent appropriate $L$ and $S$, the normalization constants, and
the results of integrating the spherical harmonics over $\theta$ and
$\varphi$, are known as the Slater coefficients. Tables of Slater
coefficients are available for simple cases \cite{condon}. For
complicated cases computer codes are available both in LS and jj
couplings. The function of the MJ92 code to be described later is to
calculate these Slater coefficients not only for energy levels, but also
for various matrix elements to be calculated for collision cross
sections.
When the integration over the angular variables and summation over the
spin are performed, the total energy of an atom is expressed in terms of
radial integrals and the Slater coefficients:
\begin{eqnarray}
E &=& \sum_i \sum_{a,b} I(ab) \langle P_a(r_i) \mid \hat{h_i} \mid
P_b(r_i) \rangle \nonumber \\ &+& \sum_{i,j} \sum_{a,b,c,d} \sum_{k}
R^{k}(abcd) \langle P_a(r_i)P_b(r_j) \mid \frac {r_<^{k}}{r_>^{k+1}}
\mid P_c(r_i)P_d(r_j) \rangle ,
\label{eq:etot}
\end{eqnarray}
\noindent where subscripts $i$ and $j$ refer to individual electrons,
$a,\ b,\ c$ and $d$ stand for sets of one-electron quantum numbers, $r_<
= \min (r_i,r_j)$, $r_> = \max (r_i,r_j)$, $I(ab)$ and $R^{k}(abcd)$ are
the Slater coefficients, and $\hat{h_i}$ is the radial part of the
one-electron Schr\"{o}dinger Hamiltonian
\begin{equation}
\hat{h_i} = -\frac {1}{2} \frac {d^2}{dr_i^2} + \frac{l(l+1)}{2r_i^2} -
\frac {Z}{r_i} ,
\label{eq:oneh}
\end{equation}
\noindent with the nuclear charge $Z$.
The variational principle $\delta E =$ 0 is applied to Eq.
(\ref{eq:etot}) along with the orthonormality conditions for the radial
functions to obtain an integro-differential equation for the radial
functions $P_{nl}$, which is called the Hartree-Fock equation:
\begin{equation}
\left\{ \frac {d^2}{dr^2} + \frac {2}{r} \left[ Z - Y_{nl}(r) \right] -
\frac {l(l+1)}{r^2} - \epsilon_{nl,nl} \right\} P_{nl}(r) = \frac
{2}{r}X_{nl}(r) + \Sigma_{n^{\prime}} \epsilon_{nl,n^{\prime}l}
P_{n^{\prime}l}(r) ,
\label{eq:hfeqn}
\end{equation}
\noindent where the ``effective screening'' term $Y_{a}(r)$ is
\begin{equation}
Y_a(r) = \Sigma_{b,k} A^k(ab) Y^{k}_{bb}(r) ,
\label{eq:directeq}
\end{equation}
\noindent while the ``exchange'' term $X_{a}(r)$ is
\begin{equation}
X_a(r) = \Sigma_{b\neq a,k} B^k(ab) Y^{k}_{ab}(r)P_{b}(r) ,
\label{eq:excheq}
\end{equation}
\noindent with appropriate Slater coefficients $A^k(ab)$ and $B^k(ab)$.
The function $Y^{k}_{ab}(r)$ is defined by
\begin{equation}
Y^{k}_{ab}(r) = r \int^{\infty}_{0} U^{k}(r,s)P_a(s)P_b(s)ds
\label{eq:ykdef}
\end{equation}
\noindent where the function $U^{k}(r,s)$ is given by
\begin{eqnarray}
U^k(r,s) &=& \frac {s^k}{r^{k+1}} \hspace{0.5in} {\rm for}\ r \geq s
\label{eq:ukrs} \\ &=& \frac {r^k}{s^{k+1}} \hspace{0.5in} {\rm for}\ r
< s . \label{eq:uksr}
\end{eqnarray}
The exchange term includes $P_a$ itself in the integrand, which makes
the Hartree-Fock equation an inhomogeneous equation in $P_a$. The
orthonormality conditions for the radial functions are combined with the
variational principle by introducing Lagrange multipliers $\epsilon$.
There is a diagonal $\epsilon_{aa}$ for each radial function $P_a$, and
off-diagonal Lagrange multipliers $\epsilon_{ab}$ to insure the
orthogonality between two radial functions $P_a$ and $P_b$ with the same
orbital angular momentum $l$. Equation (\ref{eq:hfeqn}) must be solved
as coupled, inhomogeneous, integro-differential equations. Typically,
about 90\% of the computational time for a Hartree-Fock equation is
spent to satisfy the exchange term and the off-diagonal Lagrange
multipliers on the right-hand side (RHS) of Eq. (\ref{eq:hfeqn}).
Equation (\ref{eq:hfeqn}) without the RHS is homogenous in $P_a$ and it
is known as the Hartree equation. There are as many Hartree or Hartree-
Fock equations as the total number of distinct radial functions $P_a$.
These equations for individual $P$'s must be solved simultaneously. The
Hartree equation looks very much like the Schr\"{o}dinger equation for
the hydrogen atom. In this case, the diagonal Lagrange multiplier
$\epsilon_{aa}$ plays the role of an energy eigenvalue, and it is often
used as an approximate binding energy of the electrons represented by
the radial function. This approximation is known as Koopmans' theorem.
However, the correct way to calculate the binding energy of an atomic
electron is to take the difference between the total energy of the atom
with and without the electron in question.
Normally, a radial function $P_{nl}$ is assigned for all $m_l$ for a
given $l$. This is called a restricted Hartree-Fock wave function.
However, in some condensed matter and chemistry applications, radial
functions are assumed to be different depending on the spin projection.
Such a radial function is called an unrestricted Hartree-Fock wave
function.
The standard method to solve the Hartree-Fock equation is to start with
a set of approximate radial functions---often from the Thomas-Fermi
model or hydrogenic functions with screened nuclear charges---and solve
the inhomogeneous, coupled second-order integro-differential equations.
More detailed discussions of the Hartree-Fock method are found in a book
by Slater \cite{slater} and a monograph by Fischer, Brage and
J\"{o}nsson \cite{fischer}.
\subsection{Electron Correlation}
\label{sec22}
Since the independent particle model is not an exact theory for
many-electron systems, often the Hartree-Fock wave functions are not
accurate enough to reproduce atomic properties. The difference between
the predictions of the non-relativistic Hartree-Fock wave functions (or
equivalent independent particle model) and the ``exact results minus
relativistic corrections'' are called {\it electron correlation
effects.} The reader is warned, however, that this is a very loose
definition, which was sufficient 50 years ago when powerful computers
and relativistic wave functions were not available. We now know much
more about relativistic effects, such as quantum electrodynamic (QED)
corrections. Dirac-Fock wave functions for all atoms in the periodic
table can easily be calculated in minutes now with powerful PC's. The
modern definition of the ``electron correlation effects'' should be the
difference between the ``best'' predictions using the independent
particle model---which includes the Dirac-Fock wave functions with QED
corrections---and the ``exact'' results, usually meaning the ``best''
experimental results.
There are several ways to account for the electron correlation effects.
One obvious approach is to introduce an empirical or {\it ab initio}
effective potential. A core polarization potential is often used to
account for the electron correlation involving core electrons. An
effective potential with several empirical parameters may replace the
RHS of the Hartree-Fock equation, Eq. (\ref{eq:hfeqn}), to make the
equation homogeneous and eliminate the time-consuming exchange term.
The Slater approximation (also known as the Kohn-Sham or Herman-Skillman
approximation) replaces the exchange term with the integral $Y_{nl}(r)$
[Eq. (\ref{eq:directeq})] multiplied by a scaled coefficient [see Sec.
\ref{sec221}].
The exchange term is often referred to as a ``nonlocal'' potential
because it is not a simple function of $r$ as the ``effective
screening'' term, which is called a ``local'' potential. The Hartree
equation that contains only local potentials can be solved rapidly with
modern computers. Moreover, because of its hydrogen-like solutions, one
can generate a complete set of eigenfunctions, which can serve as a
convenient starting point for a perturbation series. This is a method
of choice if electron correlation contributions from the continuum
states are important.
A more traditional approach known as the configuration interaction (CI)
method is to add configurations that include excited states. For
instance, the ground state of He can be described by two configurations,
1s$^2$ + 2p$^2$. This method allows more flexibility in the wave
function. Physically, one can interpret the multi-configuration wave
function mentioned above as an indication that the two electrons in He
are in the 1s$^2$ configuration and the 2p$^2$ configuration
simultaneously. The relative ratio of the two configurations is
determined again by applying the variational principle. Since the CI
method is popular and its physical interpretation is intuitive, we
explain the method in more detail in Sec. \ref{sec222}.
\subsubsection{Slater Approximation}
\label{sec221}
To avoid the difficulties associated with the exchange term Eq.
(\ref{eq:excheq}) in the Hartree-Fock equation, Slater \cite{slater}
proposed an approximation [see Appendix 22 of his book] to simplify the
exchange term. Wave functions based on Slater's approximation is known
as the Hartree-Slater wave functions. There are some variations, such
as the Herman-Skillman wave functions \cite{herman}. Later, Kohn and
Sham \cite{kohn} derived a slightly different expression for the
exchange term using a variational method.
The exchange term in the Slater approximation is given by:
\begin{equation}
\frac {2}{r}X_{nl}(r)= - \gamma \frac{3}{2} \left[ \frac{24}{\pi} \rho
(r) \right] ^{1/3},
\label{eq:slaterexch}
\end{equation}
\noindent where $\rho(r)$ is the total charge density
\begin{equation}
\rho (r) = \Sigma_{nl} N_{nl}P_{nl}(r)^2,
\label{eq:rho}
\end{equation}
\noindent with the electron occupation number $N_{nl}$ of the $nl$
orbital.
The constant $\gamma=1$ in the original Slater approximation, while
$\gamma=2/3$ in the Kohn-Sham approximation. In many applications,
including the popular Cowan code for atomic structure \cite{cowan},
$\gamma$ is taken as an adjustable parameter. Note that all these
approximations are efforts to simplify the exchange term, Eq.
(\ref{eq:excheq}) in the Hartree-Fock equation. As a result, these
approximations will lead to solutions that are not Hartree-Fock wave
functions any more, and usually lead to results inferior than the
genuine Hartree-Fock solutions.
Regardless of what is used for $\gamma$, the Slater approximation and
its variants simplify the exchange term into a local potential, which
makes the Hartree-Fock equation into a Schr\"{o}dinger equation for a
hydrogen-like system with an effective local potential. This in turn
makes it possible to generate a complete set of one-electron orbitals
covering both the discrete and continuous spectrum to be used in
many-body perturbation theories.
The Kohn-Sham approximation is often used in solid-state physics and
quantum chemistry applications as the starting point to introduce more
electron correlation, commonly known as the density functional theory.
However, to introduce many-body correlation beyond the Hartree-Fock
method, the charge density function $\rho(r)$ must be derived from wave
functions superior than Hartree-Fock wave functions.
\subsubsection{Configuration Mixing}
\label{sec222}
The additional configurations---known as the correlation
configurations---must have the same total angular quantum numbers, $L$
and $S$ or $J$ and the parity as those of the main configuration, which
are configurations normally occupied by the atomic state in question.
Extra configurations with different angular quantum numbers or parity do
not have any interaction with the occupied configurations. In
principle, there are no limits to how many configurations one can mix as
long as a computer code/computer can handle them. In practice, however,
introducing too many configurations may lead to numerical instability
and redundancy problems.
Note that ``solving'' the Hartree-Fock equations means to determine
radial functions $P_a$ for all quantum number sets $a$, whether they are
for actually occupied orbitals or correlation orbitals. Two main
methods to obtain such solutions are either to solve $P_a$ as pure
numerical functions or by expanding $P_a$ in terms of known analytic
functions (e.g., Slater-type functions or Gaussian functions) or
numerical functions (e.g., splines). The former is known as the
numerical Hartree-Fock wave functions, while the latter is known as the
basis-set Hartree-Fock wave functions.
Both methods have their advantages and disadvantages, which we shall not
discuss in detail here. At present, most atomic wave functions are of
the numerical type, while most molecular wave functions are of the
basis-set type. The numerical Hartree-Fock wave functions are usually
more accurate than the basis-set Hartree-Fock wave functions, but the
former are more prone to numerical convergence problems when a large
number of configurations are used. The basis-set wave functions have no
convergence problems but the computation takes much longer, and the
results depend on the choice of basis sets used, i.e., there is no
unique solution for a given atomic state, even when the same number of
configurations are used.
We shall illustrate the fine points of the configuration mixing method
by an example. For the ground state of Be, we use:
\begin{equation}
\Psi({\rm Be}) = a \psi({\rm 1s}^2{\rm 2s}^2) + b \psi({\rm 1s}^2{\rm
2p}^2),
\label{eq:beci}
\end{equation}
where $a$ and $b$ are ``configuration mixing coefficients'' to be
determined by applying the variational principle to the total energy
expression, which will naturally involve $a$ and $b$. We also have an
additional choice in selecting the individual orbitals $\psi({\rm
1s}^2{\rm 2s}^2)$ and $\psi({\rm 1s}^2{\rm 2p}^2)$.
One choice is to calculate them separately by solving the Hartree-Fock
equations for each of them. Then, $a$ and $b$ are solved by the
variational principle while keeping the radial functions $P_{\rm 2s}$
and $P_{\rm 2p}$ frozen. This is known as the ``CI'' wave functions
(not to be confused with the CI {\it method}, which is another name
often used for the configuration mixing method).
Another choice is to let the variational procedure determine both the
mixing coefficients and the radial functions. In practice a set of
trial radial functions are kept frozen while the mixing coefficients $a$
and $b$ are determined. Then, new $Y^k$ functions are calculated with
the new mixing coefficients and the Hartree-Fock equations are now
solved for the radial functions while the mixing coefficients are kept
frozen. This process is repeated until the mixing coefficients and the
radial functions have satisfied the convergence criteria for both. This
type of wave functions are called the ``multi-configuration''
Hartree-Fock (MCHF) wave functions. Obviously, MCHF wave functions are
more compact, flexible and superior wave functions than CI wave
functions, but MCHF wave functions will not converge well if too many
correlation configurations are introduced.
Radial functions obtained in this manner are called ``relaxed'' radial
functions, in contrast to the ``frozen'' radial functions mentioned
earlier. The {\it relaxed} radial functions peak at the same place the
radial functions in the main configuration peak. In the above example,
the ``relaxed'' 2p radial function obtained through the MCHF procedure
will peak near where the 2s orbital peaks, while the peak position of
the ``frozen'' 2p radial function will not change as a result of the
configuration maxing procedure.
Another example of the ``relaxed'' and ``frozen'' radial functions is
the {\it term-dependent} solutions. The excited configuration, 2s2p, of
Be can form either $^3$P or $^1$P states because of the two spin angular
momenta involved. One can calculate radial functions 2s and 2p either
for the specific total spin state, or for an average of the triplet and
singlet spin states, and then calculate the total energy for a specific
total spin while keeping the radial functions frozen.
The total spin specific radial functions are called ``term-dependent''
radial functions, while the spin-independent radial functions are called
``configuration average'' radial functions. The term-dependent 2p
radial functions for the triplet and singlet states of Be are very
different. For better results, term-dependent radial functions should
be used, though it takes longer to calculate term-dependent radial
functions.
\subsubsection{Dirac-Fock Equation}
\label{sec223}
The Dirac-Fock method is a relativistic version of the Hartree-Fock
method described so far. The one-electron Schr\"{o}dinger Hamiltonian,
Eq. (\ref{eq:oneh}) is replaced by the one-electron Dirac Hamiltonian
(in atomic units):
\begin{equation}
h_D = \bbox{\alpha}\cdot \bbox{p} c + c^2(\beta_D -1) -Z/r,
\label{eq:dirach}
\end{equation}
\noindent with 4$\times$4 matrix operators $\bbox{\alpha}$ and $\beta_D$
defined by
\begin{eqnarray}
\bbox{\alpha} = \left( \begin{array}{cc}
0 & \bbox{\sigma}_P \\ \bbox{\sigma}_P & 0 \end{array} \right)
\hspace{1.0in} \beta_D = \left( \begin{array}{cc}
I & 0 \\ 0 & -I \end{array} \right)
\label{eq:diracab}
\end{eqnarray}
\noindent in terms of the usual 2$\times$2 Pauli spinor
$\bbox{\sigma}_P$, and the 2$\times$2 unit matrix $I$.
The one-electron Dirac wave function associated with $h_D$ is given by
\begin{eqnarray}
\psi_{n\kappa m}(r,\theta,\varphi) = \frac {1}{r} \left( \begin{array}
{cc} P_{n\kappa}(r) \chi_{\kappa m}(\theta,\varphi) \\ iQ_{n\kappa}(r)
\chi_{-\kappa m} (\theta,\varphi) \end{array} \right),
\label{eq:diracpsi}
\end{eqnarray}
\noindent where $P_{n\kappa}$ and $Q_{n\kappa}$ are the radial functions
for the large and small components, respectively, $\kappa$ is the Dirac
quantum number that combines $l$ and $j$, and the two-component spinor
$\chi_{\kappa m}$ combines the usual spherical harmonics
$Y_{lm_l}(\theta,\varphi)$ with the two-component spinors $\hat{\alpha}$
and $\hat{\beta}$ defined earlier, Eq. (\ref{eq:spinor}).
The Dirac quantum number $\kappa$ is defined by (see Table
\ref{tbl:table2}):
\begin{equation}
\kappa = \left\{ \begin{array}{cc} -l - 1 & {\rm if}\ j = l+1/2, \\
+l & {\rm if}\ j = l-1/2, \end{array} \right.
\label{eq:kappa}
\end{equation}
\begin{equation}
j = \mid \kappa \mid - 1/2.
\label{eq:jkappa}
\end{equation}
The relativistic correction to the Coulomb repulsion between two bound
electrons, $e^2/r_{jk}$, is known as the Breit operator, $H_B$:
\begin{equation}
H_B = \Sigma_{j,k>j} \left[ \frac {\bbox{\alpha}_j\cdot\bbox{\alpha}_k}
{r_{jk}} \cos(\omega_{jk}r_{jk}) + \frac {(\bbox{\alpha} \cdot \nabla)_j
(\bbox{\alpha} \cdot \nabla)_k} {r_{jk}} \frac {\cos(\omega_{jk}r_{jk})
- 1} {\omega_{jk}^2} \right],
\label{eq:fullbreit}
\end{equation}
\noindent where $\omega=(\epsilon_j - \epsilon_k)/c$ is the frequency of
the virtual photon exchanged between the interacting electrons in the
framework of QED. The variable $\omega$ is rigorously defined only when
the unperturbed system is assumed to be a collection of non-interacting
electrons, which is not the case for many-electron atoms.
The Breit operator with $\omega \neq 0$ is called the {\it
frequency-dependent} Breit operator, while its $\omega = 0$ limit is
called the {\it zero-frequency} Breit operator. The frequency-dependent
Breit operator becomes ill-defined in multiconfiguration Dirac-Fock
(MCDF) calculations, and hence is not recommended in a routine
application of the MCDF method.
There are still unresolved fundamental questions about the use of the
Breit operator as part of the unperturbed Hamiltonian for a
many-electron system, so that the radial functions are fully influenced
by the Breit operator. Hence for routine applications, it is safer to
use the Dirac-Coulomb Hamiltonian $H_{DC}$ (often called the ``no
virtual $e^+$-$e^-$ pair Hamiltonian'') which is given by
\begin{equation}
H_{DC} = \Sigma_j \left[ h_D(j) + \Sigma_{k>j} e^2/r_{jk} \right],
\label{eq:hdc}
\end{equation}
\noindent and use $H_B$ in the first-order perturbation after the radial
functions have been determined by using $H_{DC}$.
As in the case of the nonrelativistic Hartree-Fock method, Slater
determinants are built with the Dirac orbitals $\psi_{n\kappa m}$ to get
the proper total $\bbox{J}$ using jj-coupling. The expression for the
total energy is derived by sandwiching $H_{DC}$ with the eigenfunction
of $\bbox{J}$, the spin-angular parts [$\chi$ in Eq.
(\ref{eq:diracpsi})] are integrated analytically, and the variational
principle is applied to the total energy expression in terms of radial
functions.
The Hartree-Fock equation consists of a set of second-order differential
equation for each radial function $P_{nl}$. The matching Dirac-Fock
equation consists of a pair of coupled first-order differential
equations for the large component and the small component radial
functions $P_{n\kappa}$ and $Q_{n\kappa}$ for each Dirac orbital
$\psi_{n\kappa}$. There are as many pairs as the total number of
distinct radial functions. It is customary, as is done in the
Hartree-Fock method, to assume all one-electron functions with the same
$n$ and $\kappa$ but with different $m=m_j$ share the same radial
functions.
After the application of the variational principle, one gets the
following Dirac-Fock equation:
\begin{eqnarray}
\left[ \begin{array} {cc} \frac {d}{dr} + \frac {\kappa_A}{r} - \frac
{\epsilon_{AA}}{c} & \ \ \ \frac {V_A(r)}{c} - 2c \\ - \frac
{V_A(r)}{c} & \ \ \ \frac {d}{dr} - \frac {\kappa_A}{r} - \frac
{\epsilon_{AA}}{c} \end{array} \right] \left( \begin{array} {cc} P_A(r)
\\ Q_A(r) \end{array} \right) = \nonumber \\ \Sigma_{\kappa_B=\kappa_A}
\frac{\epsilon_{A,B}}{c} \left( \begin{array} {cc} Q_B(r) \\ -P_B(r)
\end{array} \right) + \left( \begin{array} {cc} X_{Q_A}(r) \\ X_{P_A}(r)
\end{array} \right),
\label{eq:diracfock}
\end{eqnarray}
\noindent where, in analogy to the Hartree-Fock equation
(\ref{eq:hfeqn}), $A$ and $B$ represent the one-electron quantum number
set $n\kappa$, $V_A(r)$ is the ``screened nuclear charge'' term
equivalent to $(2/r)[Z-Y_A(r)]$ in Eq. (\ref{eq:hfeqn}),
$\epsilon_{A,B}$ are the Lagrange multipliers, the $X_{P,Q}(r)$ are the
exchange terms equivalent to $(2/r)X_{nl}(r)$. Both $V_A(r)$ and
$X_{P,Q}(r)$ contain relativistic Slater coefficients resulting from the
integration of spin-angular functions and the Clebsch-Gordan
coefficients to make the total wave function an eigenfunction of $J^2$.
As is common in Hartree-Fock calculations, Dirac-Fock wave functions can
also be classified into the configuration average, or term-specific wave
functions, except that the relativistic configuration average is
$(2J+1)$-weighted whereas the nonrelativistic configuration average is
$(2S+1)(2L+1)$-weighted. The total $J$ for a relativistically closed
shell (e.g., 2p$_{1/2}^2$, 3d$_{5/2}^6$) is 0, and the corresponding
wave function is a single determinant. Similarly, alkali- and
halogen-like configurations can be represented by a single Slater
determinant of relativistic one-electron orbitals.
Conversion between $J$ and $LS$ is not unique when the number of
open-shell electrons is three or more. Even in a simpler case such as
2s2p, the matching relativistic configurations are 2s2p$_{1/2}$ and
2s2p$_{3/2}$, resulting in one $J=2$ level, two $J=1$ levels, and one
$J=0$ level. The $J=2,0$ levels and one of the two $J=1$ levels belong
to the $^3$P term, while the other $J=1$ level belongs to the $^1$P
term. Since $LS$ labels are not used in the jj-coupling, the two $J=1$
levels are usually identified as the ``first'' and ``second'' $J=1$
levels, starting from the lowest one.
The 2s2p nonrelativistic configuration is a good example that a matching
relativistic representation requires a larger number of relativistic
configurations. The nonrelativistic 2p$^2$ configuration for the ground
state of carbon is another example: the equivalent relativistic
configuration consists of three---2p$_{1/2}^2$, 2p$_{1/2}$2p$_{3/2}$,
and 2p$_{3/2}^2$. The 2p$_{1/2}$ and 2p$_{3/2}$ orbitals are not
equivalent any more, and hence the 2p$_{1/2}$2p$_{3/2}$ relativistic
configuration may have all possible $J$ values, viz., $J=1$ and 2, while
2p$_{1/2}^2$ can only have $J=0$ because it is a relativistically closed
shell.
\section{Outline of Collision Theory}
\label{sec3}
\subsection{General Comments}
\label{sec31}
Atomic collision theory takes diverse formats depending on the type of
incident particles. For photon-atom collisions, one can get accurate
results if accurate target wave functions are available. In other
words, the accuracy of photoexcitation or photoionization cross sections
is entirely dependent on the ability to calculate good target wave
functions.
Theories for electron-atom collisions are divided into two classes
according to the treatment of the interaction between the incident
electron and the bound electrons in the target. Most theories treat the
interaction as a perturbation to the total Hamiltonian of the target
(without the incident electron). The other type of theories treat this
interaction on the same level as the interaction among the bound
electrons and with the nucleus. A typical example of the perturbative
theories is the plane-wave Born approximation (PWBA), while the R-matrix
method is a typical strong-coupling theory.
In general, strong-coupling theories are suitable for slow incident
electrons, while perturbative theories provide reliable answers for fast
incident electrons. It is rare for a collision theory to be valid at
both low and high incident energies. The only exception is the
binary-encounter-dipole (BED)/BEB model for electron-impact ionization
of atoms and molecules. The BED/BEB model is outlined in Sec.
\ref{sec8}.
Theories for atom-atom and ion-atom collisions are the most difficult
ones, even if accurate wave functions are known for the projectile and
the target, because the nature of the interaction between the colliding
particles strongly depends on the relative velocity of the colliding
particles. For slow collisions, the two particles form a transient
diatomic molecule, and this is one of the most difficult class of
collision problems to solve. In this section, discussions will be
limited to the theory of photon-atom and electron-atom collisions, not
only because these cross sections are needed most in many applications,
but also because the atom-atom and ion-atom collisions must be treated
case by case.
\subsection{Photon-Atom Collisions}
\label{sec32}
The interaction of an atom with radiation is introduced by adding the
vector potential $\bbox{A}$ to the kinetic energy operator $\bbox{p}$ in
the Hamiltonian, i.e., $\bbox{p}^2$ is replaced by $[\bbox{p} - (e/c)
\bbox{A}]^2$. Then, $\bbox{p}\cdot\bbox{A}$ is treated as the
perturbation. The resulting transition matrix element contains
$\bbox{p}$ as the operator, which is known as the velocity form. The
momentum operator is then replaced by $\bbox{r}$ using the commutation
relation \cite{bethe57}
\begin{equation}
[\bbox{r},H] = \frac {i\hbar}{m}\bbox{p},
\label{eq:commut}
\end{equation}
\noindent where $H$ is the usual Hamiltonian
\begin{equation}
H = \Sigma_j \left[ \frac {p_j^2} {2m} + V(\bbox{r}_j) + \Sigma_{k>j}
\frac {e^2}{r_{jk}} \right],
\label{eq:nrhamilt}
\end{equation}
\noindent where the summation is over all electrons in the target.
The transition probability for going from an upper level $n$ to a lower
level $n^{\prime}$ by emitting a photon is called Einstein's $A$
coefficient:
\begin{equation}
A_{nn^{\prime}} = \frac {4e^2\omega^3} {3\hbar c^3} \mid \langle \Psi_n
\mid \Sigma_j \bbox{r}_j \mid \Psi_{n^{\prime}} \rangle \mid ^2,
\label{acoef}
\end{equation}
\noindent where the angular frequency $\omega = (E_n - E_{n^{\prime}})
/\hbar$ with the total energies $E$. This is called the dipole or $E1$
transition matrix element in the length form. Einstein's $A$
coefficient (dimension = s$^{-1}$) is the total number of photons
emitted per second. Alternatively, the dipole oscillator strength
$f_{n^{\prime}n}$ is used for photoabsorption:
\begin{eqnarray}
f_{n^{\prime}n} &=& \frac {2m\omega g_n}{3\hbar g_{n^{\prime}}} \mid
\langle \Psi_n \mid \Sigma_j \bbox{r}_j \mid \Psi_{n^{\prime}} \rangle
\mid ^2 \nonumber \\ &=& \frac {mc^3g_n}{2e^2\omega^2g_{n^{\prime}}}
A_{nn^{\prime}},
\label{eq:fvalu}
\end{eqnarray}
\noindent where $g$ is the multiplicity of the respective state. The
above definition of oscillator strength applies only to the dipole
transition. The dipole oscillator strength is often called the f-value,
and it is dimensionless.
A very useful sum rule, the Thomas-Kuhn-Reich (TKR) sum rule for the
f-value is known:
\begin{equation}
\Sigma_n f_{n^{\prime}n} = N,
\label{fsumrule}
\end{equation}
\noindent where $N$ is the total number of electrons in an atom. The
summation must include all transitions from a given state, including all
transitions to the continuum and inner-shell excitations. The TKR sum
rule is satisfied only for the entire atom, but not for each orbital.
Since the proof of the sum rule uses the nonrelativistic Hamiltonian,
Eq. (\ref{eq:nrhamilt}), it does not hold for heavy atoms.
When the photon energy is high, non-dipole transitions become
significant. Multipole transition probabilities are calculated by
expanding the photon field into multipole electromagnetic field. The
magnetic dipole transition is called the $M1$ transition, the electric
quadrupole transition the $E2$ transition, the magnetic quadrupole
transition the $M2$ transition, and so on. Each multipole transition
has a different set of selection rules. Rigorous selection rules are in
terms of $\Delta J$ and parity. The selection rules based on spin
quantum numbers are approximate and often violated in heavy atoms where
spin-orbit interaction is significant. Selection rules for the $E1,\
M1$, and $E2$ transitions are listed in Table \ref{tbl:table3}. The
PH92 code will calculate $En$ and $Mn$ transitions for all $n$, provided
that wave functions satisfy the selection rules.
When wave functions from local potentials are used, the length and
velocity forms of the f-value and hence transition probability agree
exactly. Results from Hartree-Fock wave functions will not, because
they are computed with non-local potentials. The difference between the
length and velocity forms of the f-value {\it is not} an indication of
the accuracy of the f-value. Numerical convergence of the f-value as
the wave functions are improved is the only sure way of assuring
accuracy.
\subsection{Electron-Atom Collisions}
\label{sec33}
When the incident electron energy $T$ is only a few times the
excitation/ionization threshold, the incoming electron and the target
atom form a transient compound state, or an excited negative ion, and a
successful collision theory at such low values of $T$ must take the
correlation between the incident and target electrons into account.
This is why most electron-atom collision theories are not reliable at
low $T$.
At high $T \gtrsim 10B$, where $B$ is the binding energy of the target
electron, first-order perturbation theory such as the plane-wave Born
approximation becomes reliable, provided that good target wave functions
are used. The user must keep in mind that approximations for the
collision theory and approximations for the target wave functions are
not related. The approximations used for the incident, scattered, and
ejected electrons determine the reliability of the collision theory
part. It is important to use as accurate target wave functions as
possible before attempting to use ``fancy'' collision theories.
For modeling applications, a large number of various types of cross
sections are needed. In such cases, it is important to find
``practical'' solutions of as uniform quality as possible, rather than a
mixture of a few very accurate results and poor quality results for the
remainder. In many cases, it is very difficult to calculate accurate
target wave functions for highly excited states. The ``practical''
solutions for such cases reduce to single-configuration target wave
functions used with the plane-wave Born or distorted-wave Born
approximations for the collision part.
Electron-atom collision cross sections can be divided into three basic
categories: (a) elastic scattering, (b) excitations, and (c)
ionization. For ions in a plasma environment, (d) recombination becomes
also a significant collision process. Unlike photon-atom collisions,
there is no selection rules in electron-atom collisions, though
spin-forbidden and dipole-forbidden transitions tend to be weaker than
spin-allowed and dipole-allowed transitions. This lecture note covers
only discrete excitations and ionization.
The plane-wave Born approximation (PWBA) uses plane waves for the
incident and scattered electron. The target wave functions should be as
accurate as practical. Bethe \cite{bethe30} worked out details of the
PWBA in the early days of quantum mechanics. Although he applied the
PWBA mostly to the hydrogen atom, most of his conclusions are also valid
for many-electron atoms and molecules \cite{inokuti}.
First, he showed that the essence of the PWBA is given in its form
factor, known as the generalized oscillator strength (GOS), which is a
function of the momentum transfer from the incident electron to a target
electron but independent of $T$. The GOS is uniquely determined from
the initial and final state wave functions of the target. Hence, once
the GOS for a transition is calculated, the PWBA cross section for that
transition can be deduced from the GOS for all $T$. Of course, the
quality of a GOS will depend on the quality of the target wave functions
used, but this is the simplest and easiest collision cross section to
calculate. For modeling purposes, Hartree-Fock or slightly better wave
functions will produce ``useful'' cross sections at high $T$, beyond the
cross section peak.
Bethe showed further that the asymptotic form (i.e., high $T$) of the
PWB cross sections for excitation and ionization, which is usually
called the Bethe cross section, is given by
\begin{equation}
\sigma_{\rm Bethe} = \frac {4\pi a_0^2}{T/R} \left[ M_n^2 \ln \left(
\frac {T}{R} \right) + C_n \right] ,
\label{eq:bethe}
\end{equation}
\noindent where $a_0$ is the Bohr radius, $R$ is the rydberg energy, the
dipole quantity $M_n^2$ for excitation is
\begin{equation}
M_n^2 = f_{0n}/(E/R)
\label{eq:mn2}
\end{equation}
\noindent with the dipole oscillator strength for excitation to state
$n$ from the initial state $0$, excitation energy $E$, and $C_n$, which
is also a constant characteristic of the excitation but independent of
$T$. For a dipole-forbidden transition $M_n^2$ vanishes.
{\it The logarithmic dependence of the asymptotic cross section of a
dipole-allowed transition is valid for all inelastic collisions of
charged particles (electrons, protons, ions, pions, etc.) with atoms,
molecules, and their ions. Similarly, the lack of the logarithmic term
in the asymptotic cross section for a dipole-forbidden transition is
also valid for all incident charged particles and targets.} Classical
collision theories do not lead to this logarithmic dependence, and their
asymptotic cross sections for dipole-allowed transitions tend to be too
low at high $T$. Ionization cross sections always contain both
dipole-allowed and dipole-forbidden transitions and hence the
logarithmic dependence is always observed in asymptotic ionization cross
sections.
Another important conclusion about the PWBA is that the same formula
such as Eq. (\ref{eq:bethe}) applies to heavy ions simply by using
$T=\case{1}{2}mv^2$ where $v$ is the {\it speed} of the ion. For
example, a proton with 2 MeV in kinetic energy has the speed of an
electron with about 1 keV in kinetic energy. Hence, if the
electron-impact excitation or ionization cross section of He at $T=1$
keV is known, the same value will represent the cross section for a 2
MeV proton on He. Note that this type of scaling is valid only when the
heavy ions are moving fast. A slow positive incident ion has additional
interactions that an equal-speed electron does not have, such as charge
transfer and formation of a transient molecule with the target.
The validity of the PWB cross sections at low $T$ can be extended by
altering the denominator $T/R$ on the RHS of Eq. (\ref{eq:bethe}) as is
shown in the next section on the BEB model. The general form of the
Bethe cross section is valid not only for the integrated excitation and
ionization cross sections, but also for the singly differential (=energy
distribution of secondary electrons) and the doubly differential
(=energy and angular distributions of secondary electrons) cross
sections by replacing the constants $M_n^2$ and $C_n$ with appropriate
{\it functions}.
When the target is an ion, the incident electron is subject to a
long-range Coulomb attraction, and plane waves should not be used. An
obvious alternative is to use the continuum Coulomb function with a bare
or screened nuclear charge which is a constant. This is known as the
Coulomb-wave Born approximation (CWBA). For a many-electron target,
even for a neutral target, the screened nuclear charge seen by the
incident electron is a function of the distance from the nucleus.
Hence, another alternative, though more laborious, is to use an
$r$-dependent screened charge calculated from the charge density of the
target. This is known as the distorted-wave Born approximation (DWBA).
In the DWBA, the incident and scattered electron is represented by
partial waves expanded in terms of the angular momentum $l$, i.e., there
is no closed analytic expression for the radial part of the continuum
functions. The angular part consists of spherical harmonics as in the
bound-state wave functions.
The PH92 code calculates transition probabilities and electron-impact
PWB cross sections between bound states using bound-state wave functions
generated by the DF92 code. The DW92 code calculates electron-impact
DWB excitation cross sections between bound states using the wave
functions from the DF92 code. The computer time required for the DWB
cross sections are almost two orders of magnitude longer than the PWB
cross sections for the same excitation. Moreover, the total number of
partial waves that can be used in DWB cross sections are limited to $l
\lesssim 50$ because higher $l$ requires handling numbers smaller than
$10^{-308}$ for orbitals near the origin, which is the limit of the
``real'' numbers most computers can handle. When this limit is reached,
usually PWB cross sections become valid and can replace DWB cross
sections for higher $T$.
\section{Outline of the BED/BEB Model}
\label{sec4}
As was mentioned in the beginning of this lecture series, ``almost
exact'' solutions to the electron-impact ionization of the hydrogen atom
were obtained only in the last decade, 70 years after the discovery of
quantum mechanics. The difficulty is rooted in handling correlation of
two electrons in the continuum, while accounting for the
indistinguishability of the two electrons. The conventional ionization
theory expands the incident, scattered, and ejected electron wave
functions in partial waves, requiring in principle three independent
sets of infinitely many partial waves. This type of expansion cannot be
carried out for fast incident or ejected electrons because the required
computer resources in cpu time, storage, random memory, and the limits
on the exponent range of real numbers cannot be met in any foreseeable
future. Fortunately, the PWBA---which does not expand continuum
functions---is valid at high $T$. Unfortunately, however, most
applications require ionization cross sections for slow incident
electrons.
In 1994, Kim and Rudd \cite{kim94} developed a theory that does not
depend on partial wave expansion. They combined modified Mott cross
section for the collision of two free electrons and the asymptotic Born
cross section, Eq. (\ref{eq:bethe}) in the singly differential form,
which is referred to as the binary-encounter-dipole (BED) model. The
BED model requires the knowledge of continuum dipole oscillator
strength, $df/dE$, for each orbital in the target. However, for many
atoms and most molecules, $df/dE$ is unknown. For such cases, Kim and
Rudd adopted a model $df/dE$ that has the shape of the $df/dE$ for the
hydrogen atom, and obtained a simple, analytic expression for the
integrated ionization cross section per orbital. This simplified cross
section is called the binary-encounter-Bethe (BEB) cross section. The
BEB cross section per orbital requires only three orbital constants, the
binding energy $B$, the kinetic energy $U=\langle \bbox{p}^2/2m
\rangle$, and the electron occupation number $N$:
\begin{equation}
\sigma_{\rm BEB} = \frac {S}{t+u+1} \left[ \case {1}{2}\ln t \left( 1 -
\frac {1}{t^2} \right) + 1 - \frac {1}{t} - \frac {\ln t}{t+1} \right],
\label{eq:beb}
\end{equation}
\noindent where $S=4\pi a_0^2N(R/B)^2$, $t=T/B$, and $u=U/B$. Energies
expressed in terms of the ionization threshold $B$ is said to be in the
threshold unit. The first logarithmic term on the RHS of Eq.
(\ref{eq:beb}) comes from the leading term of the Bethe cross section,
Eq. (\ref{eq:bethe}), the last logarithmic term represents the
interference between the direct and exchange (because electrons are
fermions) scattering, and the middle term, $1-1/t$, originates from the
direct (=Rutherford scattering) and exchange collisions.
The only {\it ad hoc} expression in the BED/BEB model is the
denominator, $t+u+1$ on the RHS of Eq. (\ref{eq:beb}). Note that the
Bethe cross section, Eq. (\ref{eq:bethe}), has only t ($=T/B$) in the
denominator. Originally, the $T$ in the Bethe cross section---and in
all rigorous collision theories---came from normalizing the cross
section to the incoming electron flux per unit area perpendicular to the
incident beam direction, i.e., divide the transition rate by $v^2$ where
$v$ is the speed of the incident electron.
However, Burgess proposed to modify it with the argument that the
effective kinetic energy of the incident electron seen by the target
bound electron is $T +$ potential energy of the bound electron.
Accordingly, $T$ was replaced by $T+U+B$, or $t+u+1$ in the threshold
unit. Another way to look at Eq. (\ref{eq:beb}) is to consider the
expression in the square brackets on the RHS as the basic cross section
in an independent particle model, and the Burgess denominator $t+u+1$ as
the scaling factor to represent the correlation between the two
colliding electrons. In this way, the Burgess denominator must be
adapted to different types of targets---e.g., neutral atoms versus ions,
inner versus valence shells---while the basic cross section is unchanged
as long as an independent particle model, such as the Hartree-Fock
method, is used to describe the target.
The BEB model is very successful in reproducing known ionization cross
sections for large and small molecules from threshold to several keV in
$T$ \cite{nistweb}. Until the BEB model was developed, existing methods
for estimating molecular ionization cross sections were either empirical
or additivity rules with limited applicability. The BED/BEB model has
also successfully been applied to small atoms \cite{kim00}, and
applications to large, heavy atoms and their ions are being developed
now.
\section{Description of DF92}
\label{sec5}
The DF92 code calculates bound-state Dirac-Fock wave functions. The
wave function generated may be an eigenfunction of $J^2$ or the
configuration average, i.e., $2J+1$-weighted average of total energies
of all possible $J$ that can be constructed from the occupied electrons.
Also, the wave function generated may be a single configuration wave
function or multiconfiguration wave function.
The type of wave function to be generated is determined by the input.
The basic input data includes: (a) a title line to identify the input
data, (b) electron configurations, i.e., electron occupation numbers for
each orbital, (c) Slater integrals to specify $J^2$ eigenfunction, (d)
nuclear charge, and (e) the total number of bound electrons. The last
two numbers should be the same for a neutral atom. The relativistic
configurations include the principal quantum number $n$, orbital angular
momentum $l$, and $j$ quantum number. The row marked ``Input'' in Table
\ref{tbl:table2} shows the $lj$ combination to be used in the input data
for all computer codes described in this lecture note.
The DF92 code internally generates the appropriate total energy
expressions for configuration average wave functions. To generate an
eigenfunction of a specific $J^2$, appropriate Slater integrals, which
are generated from the MJ92 code described in Sec. \ref{sec6}, must be
provided in the input data. The configuration average is the correct
$J^2$ eigenfunction for a relativistically closed-shell configuration,
alikali-like configuration, and halogen-like configuration. When the
configurations used in an MCDF calculation can lead to more than one
eigenfunctions of the same $J^2$, as in the case of 2s2p*+2s2p, which
leads to either $^3$P$_1$ or $^1$P$_1$, the input data must specify the
order of the eigenfunction, e.g., the first or second $J=1$
eigenfunction counting from the lowest.
Accounting for electron correlation, i.e., the many-body effect beyond
the Hartree-Fock method, remains to be the most difficult challenge for
both relativistic and nonrelativistic atomic structure theories. The
MCHF and MCDF methods can easily account for the leading correlation
terms, which are configurations within the Coulomb complex. The Coulomb
complex consists of all combinations of one-electron orbitals that have
the same energy in the hydrogenic limit and the same total parity.
For example, the ground state of Mg is 3s$^2$. The Coulomb complex
consists of 3s$^2$, 3p$^2$, and 3d$^2$ in the nonrelativistic notation,
or 3s$^2$, 3p*$^2$, 3p$^2$, 3d*$^2$, and 3d$^2$ in the relativistic
notation. For highly charged ions, the relativistic Coulomb complex
reduces to 3s$^2$ and 3p*$^2$ because not only $n$ but also $j$ must be
the same to be the same energy in the Dirac theory for the hydrogen atom
\cite{weiss}. For neutral and lightly charged ions, however, the
nonrelativistic Coulomb complex should be used even in relativistic
calculations.
For the calculation of excitation energies, i.e., the energy difference
between two levels, the user should try to include as many correlation
configurations as possible, though it is difficult to obtain convergence
if too many configurations are introduced. An important concept is to
``balance'' the correlation in the initial and final states.
``Unbalanced'' calculations will lead to unreliable results. When an
accurate binding energy of an orbital is desired, the user should take
the difference in the total energies of the initial state and the final
state with one electron removed from the orbital.
Many atomic energy levels are often known to 6--8 significant figures
from experiments. Theoretical calculations of energy levels can compete
with the best experiments only in one-electron and two-electron atoms.
However, isoelectronic trends are well reproduced by relativistic
theories, and hence it is advantageous to study systematics in energy
levels along isoelectronic series \cite{baik}.
In the case of transition probabilities, best theories can compete with
the best experiments, particularly for the transitions from the ground
state to low-lying excited states. Although occasionally experimental
transition probabilities with an accuracy of $\pm1$\% for strong, dipole
and spin allowed transitions are reported in the literature, an
experimental accuracy of $\pm10$\% is common, and this is the accuracy
theory can achieve.
When many transition probabilities are needed, such as in plasma
modeling, consistency among the same kind of data are as important as
accuracy. To maintain consistency, wave functions of similar quality
should be used throughout. For this purpose, ``configuration average''
wave functions are recommended because they are easy to generate, i.e.,
without convergence difficulties. Then, the configuration average
radial functions are kept frozen and used in the CI mode to generate
eigenfunctions of $J^2$ as needed. Another advantage of using
configuration average wave functions is that they are all automatically
orthogonal to each other, an important requirement when many levels of
similar energies are involved. Schmidt orthogonalization alone will
produce mostly useless wave functions. Orbitals in a many-electron wave
function must be orthogonalized through Lagrange multipliers and the SCF
procedure.
All SCF codes for atomic wave functions require trial radial functions
at the start. The DF92 code has options to generate trial radial
functions using the Thomas-Fermi method or hydrogenic functions with
screened charges. It can also read an external input for radial
functions. For instance, configuration average orbitals can be
generated first, saved as binary files, and then used as the input for
an MCDF or CI calculation for a given $J$.
For deep inner shells, quantum electrodynamic (QED) corrections are
substantial. The DF92 code estimates these corrections automatically.
The ionization energies of inner-shell electrons become reliable only
when these QED corrections are included. In addition, the DF92 code
uses an extended nucleus: a uniform charge distribution for light nuclei
and a Fermi distribution for heavy nuclei. An option for using a point
nucleus is also available.
Aside from the total energy, the standard output from the DF92 code
lists the convergence data for each cycle of iteration, the orbital and
kinetic energies of each orbital, expectation values of $\langle r^n
\rangle$, configuration weights for a multiconfiguration or CI wave
function, and a list of various terms contributing to the total energy.
The radial functions and the total charge density are tabulated and
``punched'' as an option. The input format for the DF92 code and sample
input data are listed in Sec. \ref{sec9}.
\section{Description of MJ92}
\label{sec6}
The MJ92 code is an all-purpose code for integrating the angular part of
expressions for total energies, $En$ and $Mn$ transition probabilities,
photoionization, PWB cross sections for discrete excitations by electron
impact, and DWB cross sections for excitations and ionization by
electron impact. The integrated results are generically called Slater
coefficients.
The DF92 code assumes that the wave functions are linear combinations of
Slater determinants each with definite parity, $J$, and $M_J$.
Expectation values and matrix elements between the DF wave functions
generated by DF92 must be calculated using the same phase convention
used for both the wave functions and matrix elements. Therefore, it is
crucial that the user generates the necessary Slater coefficients using
the same phase convention in handling Clebsch-Gordan coefficients.
Countless mistakes were made in the history of atomic structure
calculations because mixed phase conventions were used. One should
never use formulas from two different books on angular momentum
coupling, unless it is proven beyond any doubt that the books use the
same phase convention.
For convenience, the input configurations for the MJ92 code are in
nonrelativistic notation, e.g., 2p$^2$. Then the code generates all
relativistic counterparts, 2p*$^2$, 2p*2p, and 2p$^2$, and evaluates
appropriate matrix elements between relativistic configurations. The
output identifies the angular coefficients and matching radial
integrals. These output data are placed in a file, which can be read by
the DF92, PH92, and DW92 codes, or if the list is short, they can be
incorporated directly into the input file for the DF92, PH92, and DW92
codes.
For an isoelectronic sequence, the same output file from the MJ92 code
can be used by ions with different nuclear charges. Also, an output
file from the MJ92 code may contain Slater coefficients for many
different energy states or transitions. Different cases are
distinguished by the first four characters on the first line of an
output data set. A general input convention for all codes described in
this lecture note is that the start of an input data set is identified
by an asterisk ``*'' in the first column, and a 3-digit number. The
first four characters of all input lines are often used as identifiers.
The MJ92 code also generates generic input data for the DF92 code
listing the appropriate relativistic configurations to be used. The
ordering of the configurations listed in this generic input data must be
kept to make the phase conventions used by the DF92 and MJ92 codes
consistent. The input format and sample input data for the MJ92 code
are listed in Sec. \ref{sec9}.
\section{Description of PH92}
\label{sec7}
The PH92 code calculates $En$ and $Mn$ transition probabilities for
discrete transitions. The code also generates appropriate continuum
partial waves for photoionization and calculates photoionization cross
sections. The partial waves, however, are calculated in the field of
frozen core orbitals, though the electron exchange interaction is
treated correctly. This is known as the continuum Dirac-Fock wave
functions. Since this is not a correlated continuum wave function, the
results for photon energies near the ionization threshold are
unreliable.
The PH92 code can also calculate plane-wave Born cross sections for
discrete excitations by electron impact. As is shown in Sec.
\ref{sec9}, the keyword in the input data identifies which type of cross
section is calculated.
To use the PH92 code, appropriate Slater coefficients must be generated
using the MJ92 code and Dirac-Fock wave functions using the DF92 code.
As was mentioned in the theory section, the accuracy of the $En$ and
$Mn$ transition probabilities is determined by the accuracy of the
Dirac-Fock wave functions used. The accuracy of the Born cross
sections, however, is governed by the approximation used in the Born
approximation in addition to the accuracy of the Dirac-Fock wave
functions used. The PWB cross sections are reliable only at high
incident electron energies $T$ beyond the cross section peak.
Also, the shape of the PWB cross sections near the threshold for
discrete excitations of ion targets is wrong. The PWB cross section
starts from zero, while the correct cross section starts from a finite
value. The distorted-wave Born (DWB) cross section described in the
next section has the correct shape near the threshold.
For neutral targets, both the PWB and DWB cross sections have the
correct shape, i.e., start from zero. The PWBA overestimates the cross
section between the threshold and the peak by as much as 50\%. For a
quick estimate, the Burgess denominator mentioned in Sec. \ref{sec4} may
be used to scale the PWB cross sections at low $T$. The DWB cross
sections should not be scaled by the Burgess denominator because the DWB
cross sections become too low when scaled by the Burgess denominator.
The input format and sample input data for the PH92 code are listed in
Sec. \ref{sec9}.
\section{Description of DW92}
\label{sec8}
The DWB cross sections take almost two orders of magnitude more cpu time
than the PWB cross sections to calculate. Therefore, DWB cross sections
should be used only for the discrete excitations of ion targets by
electron-impact. For the excitations of neutral targets, PWB cross
sections scaled by the Burgess denominator are recommended. Besides,
the DWB cross sections can be calculated only for low $T$ because of the
limitations in partial wave expansions mentioned earlier.
For ionization, DWB cross sections take several orders of magnitude
longer in cpu time than the BEB model described in Sec. \ref{sec4}. For
ionization, the correct shape near the threshold is given by both the
DWBA and the BEB model, and the BEB cross sections are more reliable
than the DWB cross sections at low $T$.
The DW92 code generates continuum partial waves for the incident,
scattered, and ejected electrons. The partial waves are generated in
the field of frozen core orbitals with proper electron exchange
interaction, as is also done in the PH92 code. As in the case of
bound-state orbitals, continuum partial waves can be made $J$-specific
(=term dependent) or in the configuration average because it is certain
that either the initial or the final state contains open-shell
configurations. Since the continuum partial waves are not true SCF wave
functions---core orbitals are frozen during the SCF process---it is not
worth elaborating on term-dependent partial waves. This is the reason
that configuration average partial waves are the default. The
term-dependent partial waves are seldom needed, and the results from
them are not necessarily more reliable.
Computational difficulties limit the use of the DWBA to low $T$, about 5
times the threshold energy in $T$. For higher $T$, either the
Coulomb-wave Born approximation or the PWBA should be used. Since the
largest differences between the DWBA and PWBA results are found in
partial waves with low $l < 5$, the user can combine DWB partial wave
cross sections for low $l$ and PWB partial wave cross sections for high
$l$.
Convergence in the DWB cross sections must be confirmed by calculating
partial wave cross sections for high $l$. Since the PWB cross sections
from the PH92 code are calculated without partial wave expansions, the
DWB cross sections can be combined with the PWB cross sections from the
PH92 code after subtracting the low $l$ partial wave cross sections
included in the DWB cross sections. For this purpose, the DW92 code
provides options to calculate the CWB and PWB cross sections in partial
wave expansion. Small differences may persist between the DWB cross
sections from the DW92 code and the PWB cross sections from the PH92
code because the latter does not include the exchange contribution,
which is significant for low $l$ partial waves. The input format and
sample input data for the DW92 code are listed in Sec. \ref{sec9}.
\centerline{\bf Exercises}
To familiarize with the DF92, MJ92, PH92, and DW92, the user should try
the following set of exercises.
{\bf Exercise 1.} Calculate Hartree-Fock wave functions for Be:
1s$^2$2s$^2$ $^1$S and 1s$^2$2p$^2$ $^1$S separately. Then, calculate
the total energy of the ground state of Be by using (a) 1s$^2$2s$^2$
only, and (b) 1s$^2$[a2s$^2$+b2p$^2$] keeping the radial functions
frozen.
{\bf Exercise 2.} Repeat step (b) above by relaxing radial functions,
i.e., use the MCDF option. Compare the total energies and $\langle r
\rangle$ of the 2s and 2p radial functions.
{\bf Exercise 3.} Calculate the total energies for Be 2s2p $^3$P and
2s2p $^1$P states using term-dependent and configuration average radial
functions. Compare total energies and $\langle r \rangle$.
{\bf Exercise 4.} Calculate $E1$ transition probabilities for the 2s2p
$^3$P $\rightarrow$ 2s$^2$ $^1S$ and 2s2p $^1$P $\rightarrow$ 2s$^2$
$^1S$ transitions for Be. Try the same transitions for Xe$^{50+}$ and
compare the ratios of the triplet and singlet transitions.
{\bf Exercise 5.} Calculate the PWB excitation cross sections for the
2s$^2$ $^1S$ $\rightarrow$ 2s2p $^3$P and $^1$P excitations for Be and
Xe$^{50+}$.
{\bf Exercise 6.} Calculate the same transitions in Ex. 5 using the
DWBA.
{\bf Exercise 7.} Calculate the BEB cross section for the ionization of
He and compare to the experimental data by Shah et al., J. Phys. B {\bf
21}, 2751 (1988). The kinetic energy of the 1s orbital of He can be
obtained by calculating the Dirac-Fock wave function of He.
\newpage
\section{Input Format for MJ92, DF92, PH92 and DW92}
\label{sec9}
%This version of input format for MJ92 and DF92 prepared by Y.-K. Kim
\subsection{Introduction}
\label{sec91}
The MJ92 code calculates angular coefficients for the total energy
expression (including the Breit operator), and the DF92 code solves
multiconfiguration Dirac-Fock (MCDF) equations. The MJ92 code also
calculates the Slater coefficients for arbitrary $En$ and $Mn$
transition probabilities, PWB cross sections, and DWB cross sections.
Although the MJ92 code can handle angular momentum coupling and
integration for several other processes, such as the hyperfine
interaction, Auger transitions, and parity nonconservation, we shall not
cover such topics in this lecture note. The features in MJ92 and PH92
for these extra processes have been omitted in the following description
of the input format.
For each set of input data, variables to be read are listed either in
upper or lower case with the following conventions:
\begin{itemize}
\item Variables listed in upper case must be given a value.
\item For variables listed in lower case default values are
automatically assigned if no value is given as input.
\item Default values are given in parentheses after the definition of
the variables.
\item Lines beginning with a \# in the first column are treated as
comments.
\end{itemize}
All input data are read in free format using the {\bf STR...} set of
subroutines. The conventions used are the following :
\begin{itemize}
\item Input stream may include an arbitrary number of data sets, one for
each case to be calculated. A data set consists of a certain number of
files each file being divided into an arbitrary number of records.
\item A file is defined as an arbitrary number of input lines followed
by a line with {\bf END} as its first 3 characters.
\item A record is defined as an arbitrary number of input data. A
record may extend over an arbitrary number of lines and more than one
record may be on the same line. The colon ``:'' after at least one blank
may be used to interrupt data list of a record and let the program
assign default values to the rest of the data in the record.
\item An input line is a set of 72 characters followed by a carriage
return.
\item The first 4 characters of each input line are considered as an
identifier that can be used to identify a set of input data.
\item An input value is defined as a string of characters {\em without
blanks}.
\item Strings are separated by at least one blank.
\item An input value and its optional identifier (see below) must be on
the same line.
\end{itemize}
For convenience each input value can be labeled with an identifier
followed by the character ``=''. Both the identifier and the character
``='' are considered as comments and are thus ignored ({\em These are
not list-directed input!)}. Valid input data are: strings of characters
(case insensitive), integers or real numbers.
In all codes, atomic states are defined as a superposition of electronic
configurations each of which is defined as a list of orbitals with their
electron occupation numbers. These configurations are given, only in
the case of the MJ92 code, in LS coupling and the program will generate
all jj configurations arising from a given LS configuration. The rest
of the codes use jj coupling for input data.
Note that:
\begin{itemize}
\item The configurations given as input data in LS coupling use 1s, 2s,
2p and so on as labels.
\item The jj configurations generated by the MJ92 use 1s, 2s, 2p*, 2p
and so on as labels. The * denotes $j = l - \case{1}{2}$, while no *
stands for $j = l + \case{1}{2}$.
\item When a certain number of orbitals are common to all
configurations, they can be given only once at the beginning as a record
with CORE in the first 4 columns. Also if some or all of these core
orbitals make up a rare gas configuration, it is sufficient to give the
symbol of the rare gas to generate the full set of orbitals. Valid
symbols are: NE, AR, KR, XE and RN.
\end{itemize}
\noindent As an example input for silicon will read:
\begin{verbatim}
In LS coupling In jj coupling
CORE NE 3S2 : CORE NE 3S2 :
C 1 3P2 : C 1 3P*2 :
END C 2 3P2 :
END
\end{verbatim}
The structure of input data are the same for all programs and include:
\begin{itemize}
\item The first input line associated with Fortran unit number 5. This
line contains information to be given only once and includes the
following parameters:
\begin{description}
\item[idebug:] Minimum output if = 0, and long output otherwise. (0)
\item[iread:] Fortran unit number to read the next sets of input data.
(5)
\item[iwri:] Fortran unit number for standard (printed) output. (6)
\item[ipun:] Fortran unit number to save formatted (=ascii) output. (7)
\item[iba:] Fortran unit number for Slater coefficients. (8)
\item[ibb:] Fortran unit number to read wave functions in binary form.
(9)
\item[ibc:] Fortran unit number to write wave functions in binary form.
(9)
\item[nblipa:] Number of output lines per page. (57)
\item[npow] Power of 10 by which the speed of light is to be multiplied.
(0)
\end{description}
\item The second input line contains the following item:
\begin{description}
\item[ttot:] If positive, overall cpu time for the job in seconds. (0.)
\end{description}
\item Specific data for each case to be calculated.
\item A line with two asterisks as the first two characters indicating
the end of the input data for the case(s) to be calculated.
\end{itemize}
\newpage
\subsection{Input Format for MJ92}
\label{sec92}
\centerline{\bf Files Used}
\begin{itemize}
\item Initial input file (initially Fortran unit number 5)
\item Optionally a second input file (defined by IREAD as the unit
number)
\item Standard output file (defined by IWRI)
\item One or two files to save results (defined by IPUN and IBA which
may be the same)
\item Two scratch files (defined by IBB and IBC)
\end{itemize}
\centerline{\bf Input Data Provided Only Once at the Beginning}
\subsubsection {\bf 1st set}
\begin{description}
\item[idebug:] Minimum output if = 0, and long output otherwise. (0)
\item[iread:] Fortran unit number to read the next sets of input data.
(5)
\item[iwri:] Fortran unit number for standard (printed) output. (6)
\item[ipun:] Fortran unit number to save formatted (=ascii) output. (7)
\item[iba:] Fortran unit number for Slater coefficients. (8)
\item[ibb:] Fortran unit number to read wave functions in binary form.
(9)
\item[ibc:] Fortran unit number to write wave functions in binary form.
(9)
\item[nblipa:] Number of output lines per page. (57)
\end{description}
\subsubsection {\bf 2nd set}
\begin{description}
\item[ttot:] If positive, overall cpu time for the job in seconds. (0.)
\end{description}
The following data sets may be in different files depending on the
values of {\bf IREAD, IBA} read in the 1st and 2nd sets above.
\vspace{0.2in}
\centerline {\bf Input Data for Each Case to Be Calculated}
\subsubsection {\bf 1st set}
\begin{description}
\item [ITITLE:] Title line with one asterisk as its first character and
a maximum of 79 characters. The rest of the input on this line is
ignored.
\end{description}
\subsubsection {\bf 2nd set}
This option is not used in MJ92. A colon ``:'' is required.
\subsubsection {\bf 3rd set}
\begin{description}
\item[KEYWORD :] Keyword for the type of expectation values to
calculate. Options are:
\subitem {ENERGY} for the total energy
\subitem {EXCITATION} for the DWB cross section for electron-impact
excitation
\subitem {IONIZATION} for the DWB cross section for electron-impact
ionization
\subitem {ELECTRIC} for electric multipole transition
\subitem {MAGNETIC} for magnetic multipole transition
\subitem {BORN} for the PWB cross section for electron-impact excitation
\subitem Only the first 4 characters are meaningful
\item [ipun1:] File unit number to save Slater coefficients. (8)
\end{description}
For photoionization, the required Slater coefficients are identical to
those for $En$ and $Mn$ discrete transitions by photon, because the only
difference is in radial functions. The required Slater coefficients for
photoionization can be generated by using the ELECTRIC or MAGNETIC
option, and rename the final-state orbitals. Also, note that the role
of the initial and final states switches. The sample input data for the
PH92 code show examples of properly labeled Slater coefficients for
photoionization.
\subsubsection {\bf 4th set}
\noindent {\em Only for ELECTRIC or MAGNETIC options:}
\begin{description}
\item[MULTPOL :] Multipolarity. For instance, MULTPOL=1 means dipole,
MULTPOL=2 means quadrupole.
\end{description}
\subsubsection {\bf 5th set}
\begin{description}
\item[JJTOT :] Twice the value of the total angular momentum $J$.
\item[mjjtot :] Twice the projection $M=J_z$ of the total angular
momentum $J$. (JJTOT)
\end{description}
\subsubsection{\bf 6th set}
\begin{description}
\item [LIST:] List of configurations in LS coupling.
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
Note that:
\begin{itemize}
\item The {\bf configurations must be in LS coupling.} The MJ92 code
will automatically generate all matching jj configurations.
\end{itemize}
\subsubsection{\bf 7th set}
\noindent {\em Only for EXCITATION or IONIZATION option}:
\begin{description}
\item [LABEL:] Label for partial waves, only the first 3 characters will
be used.
\item [LADONE:] Angular coefficients that have already been generated
for orbital angular momentum $l$ up to and including LADONE. If nothing
was done earlier, a negative number. This option is used when the
generation of Slater coefficients takes long on a slow computer.
\item [LQMAX:] Maximum value of the orbital angular momentum $l$ for the
expansion in partial waves.
\end{description}
\noindent The MJ92 code will generate matching $j$ and $\kappa$ values
corresponding to the $l$ in the above input data.
\subsubsection{\bf 8th set}
\noindent {\em For EXCITATION and IONIZATION options:} The 5th, 6th and
the 7th sets above are for the initial state and must be repeated for
the final state. For the DWB cross sections, Set 7 is associated with
the incident electron while Set 8 is associated with the scattered
electron.
\noindent {\em For ELECTRIC and MAGNETIC option:} The 4th, 5th, and the
6th sets above are for the initial state and must be repeated for the
final state. {\em For BORN,} only the 5th and 6th sets above must be
repeated for the final state.
\subsubsection{\bf 9th set}
\noindent {\em Only for IONIZATION option:} LABEL, LQDONE and LQMAX with
the same definitions as Set 7 for the ejected electron.
\subsubsection{\bf 10th set}
\noindent {\em For ENERGY, EXCITATION or IONIZATION option}: Indicate
whether to include or not the Breit operator (for EXCITATION and
IONIZATION options, only the magnetic part of the Breit operator will be
included). Must be Y or N.
\subsubsection{\bf 11th set}
\noindent {\em Only for EXCITATION or IONIZATION option}: If Y, the
Dirac-Fock equation for partial waves is reduced to the configuration
average only. If N, then extra angular coefficients are generated.
Must be Y or N.
\vspace{0.2in}
\centerline{\bf End of Input Data}
\noindent Two asterisks as the first 2 characters. This indicates the
end of all input data.
\newpage
\subsection{Input Format for DF92}
\label{sec93}
\centerline{\bf Files Used}
\label{sec931}
\begin{itemize}
\item Initial input file (Fortran unit 5)
\item Optionally a second input file (defined by IREAD as the unit
number)
\item Optionally a third input file (defined by IBA)
\item Standard output file (defined by IWRI)
\item One file to save formatted (=ascii) results (defined by IPUN)
\item One file to read radial functions in binary form (defined by IBB)
\item One file to write radial functions in binary form (defined by IBC,
may be the same as IBB)
\item One scratch file when the Breit operator is included in the SCF
process (defined internally as IBE)
\end{itemize}
\centerline{\bf Input Data Provided Only Once at the Beginning}
\label{sec932}
\subsubsection {\bf 1st set}
\begin{description}
\item[idebug:] Minimum output if = 0, and long output otherwise. (0)
\item[iread:] Fortran unit number to read the next sets of input data.
(5)
\item[iwri:] Fortran unit number for standard (printed) output. (6)
\item[ipun:] Fortran unit number to save formatted (=ascii) output. (7)
\item[iba:] Fortran unit number for Slater coefficients. (8)
\item[ibb:] Fortran unit number to read wave functions in binary form.
(9)
\item[ibc:] Fortran unit number to write wave functions in binary form.
(9)
\item[nblipa:] Number of output lines per page. (57)
\item[npow:] Power of 10 by which the speed of light is multiplied to
emulate the nonrelativistic limit $c \rightarrow \infty$.
\end{description}
\subsubsection {\bf 2nd set}
\begin{description}
\item[ttot:] If positive overall cpu time for the job in seconds. (0.)
\end{description}
The following data sets may be in different files depending on the
values of {\bf IREAD, IBA} given in the 1st set above.
\vspace{0.2in}
\centerline{\bf Input Data for Each Case to Be Calculated}
\label{sec933}
\subsubsection {\bf 1st set}
\begin{description}
\item [ITITLE:] Title line with one asterisk as its first character and
a maximum of 79 characters. {\bf Unlike the MJ92 code}, the 3-digit
number in columns 2--4 are used to identify matching Slater coefficients
in Fortran unit IBA, if IBA is defined in the 1st set above.
\end{description}
\subsubsection{\bf 2nd set}
\begin{description}
\item [LIST:] List of configurations in jj coupling.
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
Note that:
\begin{itemize}
\item {\bf Unlike in MJ92, configurations must be in jj coupling.}
\end{itemize}
\subsubsection {\bf 3rd set}
\begin{description}
\item [NEIGV :] The order of the eigenvalue for which convergence is to
be achieved. If positive it will be the NEIGV-th eigenvalue (the first
being the lowest). If negative it will be the eigenvalue having the
$|$NEIGV$|$-th configuration as its largest component.
\item [ICMUL :] Initial configuration mixing coefficients are read
(described as the 18th set) as input data if ICMUL $\neq 0$. Furthermore
if ICMUL$<0$, the weights are frozen to their initial values during all
the SCF process. ICMUL=$-2$ generates configuration average wave
functions.
\item [iprfgr :] If IPRFGR $\neq 0$, energy expression is written in the
printed output. (0)
\item [lamrq:] If LAMRQ $> 0$, initial estimates of the off-diagonal
Lagrange multipliers (ODLM) are calculated using the Rayleigh quotient
method. (0)
\item [cofrg:] Multiply initial ODLM's (if any) by the value of COFRG.
(To be read only if LAMRQ is positive.) (1.)
\end{description}
\subsubsection {\bf 4th set}
\begin{description}
\item [NZ:] Atomic number.
\item [NEL:] Number of bound electrons. NZ$>$NEL for a positive ion.
\item [NORBSC:] Orbitals to be made self-consistent. If NORBSC=0, all
orbitals are made self consistent. If $<0$, all orbitals are frozen.
If $>0$, some orbitals (list to be given as the 9th set) are frozen.
\item [iorth:] If IORTH $\neq 0$, orbitals are Schmidt orthogonalized at
the beginning and before each diagonalization of the total energy
matrix. (0)
\item [ndep:] Default option for the trial radial functions to start the
SCF process. (1)
\subitem NDEP = 1, generate trial functions from the Thomas-Fermi
potential.
\subitem = 2, generate from a local potential (to be provided as input).
\subitem = 3, hydrogenic.
\subitem = 4, read from an input file in ascii format.
\subitem = 5, read from a binary file.
\subitem $>$ 6, restart from a job terminated by time control.
\item [nlec:] Record number to read initial radial functions from a
binary file. (0)
\item [nec:] Record number to write the calculated radial functions in a
binary file. (0)
\item [ifcwf:] Convergence acceleration factors are kept fixed if IFCWF
$\neq 0$. (0)
\item[nes :] Number of trials to adjust one-electron energies while
solving the Dirac equation. (40)
\item [infini:] Not used in this version. (0)
\item [inexwf:] Defines the next radial function to be made self
consistent. If = 0, the code will choose the orbital that shows the
largest variation between successive iterations. Otherwise all
functions are solved sequentially. (0)
\item [idr1:] If $\neq 0$, radial mesh points are defined using the data
from an input binary file. Otherwise, the input data set in the 7th set
are used. (0)
\item [iprfgr:] If IPRFGR $\neq 0$, extra Slater integrals for the
total energy expression are printed only. (0)
\item [idfbeg:] If IDFBEG $= 0$, deferred correction is used from the
beginning of the SCF process. Otherwise, it is used only after the
accuracy of the radial functions is better than 0.01. (0)
\end{description}
\subsubsection{\bf 5th set}
\begin{description}
\item[NSTEP :] Number of diagonalization cycles in the SCF process. If =
0, input is replaced by 1 for a single (5 for a multi) configuration
calculation.
\end{description}
\subsubsection{\bf 6th set}
\noindent {\em Only if NSTEP $> 0$}: ISCPAR(J,N) J=1 to 3, AZPRY(1,N),
AZPRY(2,N), ACCF(1,N), ACCF(2,N), N=1 to NSTEP.
\begin{description}
\item[ISCPAR(1,.) :] Exchange potential is included only if $\neq 0$.
Furthermore, Schmidt orthogonalization is used if $< 0$. (1)
\item[ISCPAR(2,.) :] If $= 0$, omit off-diagonal Lagrange multipliers.
(1)
\item[ISCPAR(3,.) :] Number of iterations divided by number of orbitals.
(30)
\item[AZPRY(1,.) :] Nuclear charge. (=NZ)
\item[AZPRY(2,.) :] Accuracy to be satisfied by the radial functions.
(10**(NSTEP$-6-$N))
\item[ACCF(1,.) :] Mixing coefficient of the eigenvectors between
successive iterations. (1.)
\item[ACCF(2,.) :] Fraction of electron exchange to be included. (1.)
\end{description}
\subsubsection {\bf 7th set}
\noindent Option to modify radial mesh and accuracy criteria. Must be
{\bf Y} or {\bf N}. {\em Only if option is Y:}
\begin{description}
\item[H:] Mesh step size in the variable {\it t} defined as $t=\ln
(r/r_1) + ar$. (0.05)
\item[DR(1):] First tabulation point multiplied by the nuclear charge
NZ. (0.01)
\item[AMESH:] Mesh step size in the linear region, i.e. value of {\it a}
in the above definition of {\it t}. (0.01)
\item[TESTE:] Relative accuracy to be obtained for one-electron
energies. (0.000005)
\item[RAP:] Defines the accuracy to which the small component must be
matched when solving the Dirac equation. (100.)
\end{description}
\subsubsection {\bf 8th set}
\noindent Option to modify nuclear parameters. Must be {\bf Y} or {\bf
N}. If option is N, standard atomic masses are used with a uniform
nuclear charge distribution and a radius given in NUCPOT for $Z \le 45$.
For $Z > 45$, Fermi distribution is used. {\em Only if option is Y:}
\begin{description}
\item[A:] Atomic mass (used to calculate reduced mass of all particles
including electrons).
\item[FERC:] Nuclear radius in fermi. (0.)
\item[FERT:] Thickness parameter for the Fermi charge distribution of
the nucleus. (0.)
\item[NUC:] Index of the nuclear radius (1 if FERC=0, and max(NUC,41)
otherwise).
\end{description}
For the nucleus a point charge is used when FERC $\leq 0$ and a finite
one when FERC $> 0$. For a finite charge distribution, a uniform charge
inside a sphere of radius FERC is used if FERT $\leq 0$ and otherwise a
Fermi distribution written as: \\
\centerline{$Ct / ( 1. + \exp ((r-FERC)*(4\ln 3)/FERT) )$,}
$Ct$ being a normalization constant.
\subsubsection {\bf 9th set}
\noindent {\em Only if NORBSC $>0$}: Labels of the frozen orbitals and a
line with {\bf END} as the first 3 characters.
\subsubsection {\bf 10th set}
\noindent Option to modify the default options for some of the orbitals.
Must be {\bf Y} or {\bf N}.
\noindent {\em Only if option is Y}, for each orbital to be modified:
\begin{description}
\item[LABEL:] Label of the orbital to be modified.
\item[IDOI:] Type of initial orbital (see definition of NDEP in Set 6).
\item[ZH:] Only if IDOI=3, screening of the nuclear charge for
hydrogenic initial orbital. A positive number reduces the effective
nuclear charge.
\item[INDSOL(1):] Method for solving the Dirac equation. (1=normal case;
3=when frozen orbitals are involved; 4=for correlation orbitals; 5 for
using B-splines, {\bf this method may be tricky to use.}). (1)
\item[indsol(2):] If $\geq 0$, use deferred correction. If = 0, use the
solution of previous iteration to calculate deferred correction. (0)
\item[indsol(3):] If = 0, enforce number of nodes on radial functions.
(0)
\item[indsol(4):] If $\neq 0$, use the Schmidt orthogonalization before
solving the Dirac equation. (0)
\item[indsol(5):] If $\neq 0$, use the Rayleigh quotient to obtain an
estimate of the one-electron eigenvalue. If = 0, use previous value.
(0)
\item[indsol(6):] If = 0, use full exchange. If not, reduce exchange by
the coefficient given in the 6th set, {\bf ACCF(2,N)}. (0)
\item[iormax:] Maximum number of tabulation points. (Maximum value
allowed by array dimensions. This is controlled by maxpts in
length.def.)
\item[scew:] Initial value of the mixing coefficient of the orbital
between two successive iterations. (0.3)
\item[K:] Flag to distinguish between electron (0) and muon ($\neq 0$).
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\subsubsection {\bf 11th set}
\noindent Option to modify method for computing off-diagonal Lagrange
multipliers. Must be {\bf Y} or {\bf N}. {\em Only if option is Y:}
\begin{description}
\item[TESTEP:] Criteria to switch from the difference method to the sum
method. (0.001)
\end{description}
\noindent For each Lagrange multiplier:
\begin{description}
\item[LABOI(1), LABOI(2):] Labels of the orbitals involved in the
Lagrange multiplier.
\item[N:] Selection of the method to compute Lagrange multiplier. ($-1$)
\subitem N $< -1$: The difference method.
\subitem N = $-1$: The difference or sum method depending on the
occupation numbers.
\subitem N = 0: Omit Lagrange multiplier.
\subitem N = 1: The sum method.
\subitem N $> 1$: The Rayleigh quotient.
\subitem Example:
\subsubitem 0.001
\subsubitem 6D* 5D* 2
\subsubitem 5D* 6D* 0 (or try 2)
\subitem END
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\subsubsection {\bf 12th set}
\noindent {\em Only if NDEP=2}:
\begin{description}
\item [K:] Number of tabulation points for the local potential to use.
\item [VN(I) I=1 to K:] Local potential multiplied by $r$ (i.e.,
$r*V(r)$) with the format 5(1PD14.7).
\end{description}
\subsubsection {\bf 13th set}
\noindent {\em Only if IDOI=4}:
\begin{description}
\item [K:] Number of tabulation points for the orbital to read.
\item [AP(I) I=1 to 10:] Coefficients for series expansion of the large
component radial function near the origin.
\item [AQ(I) I=1 to 10:] Coefficients for series expansion of the small
component radial function near the origin.
\item [P(I) I=1 to K:] Large component radial function.
\item [Q(I) I=1 to K:] Small component radial function.
\end{description}
\noindent with the format 5(1PD14.7) for AP, AQ, P, and Q.
\subsubsection {\bf 14th set}
\noindent {\em Only if absolute value of ICMUL=1}, for each
configuration for which the weights are to be read:
\begin{description}
\item[I:] Configuration index.
\item[CW:] Initial coefficient of the configuration.
\end{description}
Note that some of the weights may already have been read from a binary
file.
\subsubsection {\bf 15th set}
\noindent{Only if the Fortran unit to read the Slater coefficients was
defined to be the same as the rest of the input data. For each line}:
\begin{description}
\item [A:] Angular coefficient.
\item [IC, JC:] Configuration indexes.
\item [K:] Multipolarity.
\item [LABEL(I) I=1 to 4:] Labels of the orbitals associated with this
radial integral.
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\subsubsection {\bf 16th set}
\begin{description}
\item [NBREIT:] Option for the Breit operator.
\subitem If $< 0$, the Breit operator is omitted.
\subitem If $= 0$, the Breit operator is calculated in the first-order
perturbation.
\subitem If $> 0$, The magnetic part of the Breit operator is included
in the SCF process at the zero-frequency limit.
\item[NRET:] If $\neq 0$, full retardation correction in the Coulomb and
Lorentz gauges are calculated.
\subitem **** This will almost double CPU time. Furthermore must be
used only for single configuration calculations. Otherwise meaningless
results may be produced. ****
\item [ibe:] Fortran unit number to write magnetic integrals when they
are used in the SCF process, i.e., when NBREIT $> 0$. (IBC+2)
\end{description}
\subsubsection {\bf 17th set}
\noindent{Only if the Fortran unit to read the Slater coefficients was
defined to be the same as the rest of the input data. For each line}:
\begin{description}
\item [A:] Angular coefficient for the magnetic interaction.
\item [B:] Angular coefficient for the retardation interaction.
\item [IC, JC:] Configuration indexes.
\item [K:] Multipolarity.
\item [LABEL(I) I=1 to 4:] Labels of the orbitals associated with this
radial integral.
\item [NEN:] If $\geq 0$, the integration of both radial variables is
from 0 to $\infty$. If $< 0$, the integration over $r_2$ is from 0 to
$r_1$.
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\subsubsection {\bf 18th set}
\begin{description}
\item[ilams:] Distance from the origin to integrate charge density for
self-energy screening, in 1/100 Compton wavelength. (30)
\end{description}
\subsubsection {\bf 19th set}
\begin{description}
\item[input(1):] If $\neq 0$, Coulomb or Breit integrals are calculated.
\item[input(2):] If $\neq 0$, one-electron integrals are calculated.
\item[input(3):] If $\neq 0$, expectation values of $\langle r^n
\rangle$ are calculated.
\item[input(4):] If $\neq 0$, one-electron radial functions are
tabulated in the standard output file.
\item[input(5):] If $\neq 0$, the total charge density is written to
Fortran unit IPUN in the ascii format. Furthermore if $> 0$, the total
charge density is also listed in the standard output file.
\item[input(6):] Number of radial functions to be written on Fortran
unit IPUN.
\item[input(7):] If $\neq 0$, the series expansion coefficients of the
radial functions near the origin are tabulated.
\end{description}
\subsubsection {\bf 20th set}
\noindent {\em Only if INPUT(1) $\neq 0$:}
\begin{description}
\item[LABELI(I) I=1 to 4:] Label of the orbitals associated with the
integral.
\item[K:] Multipolarity.
\item[NEM:] If = 0, the Coulomb integral, and Breit integral otherwise.
Furthermore if $<0$, integration over the variable S is from 0 to R
only. The integrals are defined as: $\int F(R)*UK(R,S)*G(S) dRdS$ with
\subitem
\begin{tabular}{lll}
$F = P1*P2 + Q1*Q2$ & $G = P3*P4 + Q3*Q4$ & if NEM $= 0$ \\
$F = P1*Q2$ & $G = P3*Q4$ & if NEM $\neq 0$ \\
\end{tabular}
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\subsubsection {\bf 21st set}
\noindent {\em Only if INPUT(2) $\neq 0$}:
\begin{description}
\item[LABELI(I) I=1,2:] Label of the orbitals associated with the
requested integral.
\item[IALL:] If $\neq 0$, all one-electron integrals are calculated
starting from the one defined by the first label.
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\subsubsection {\bf 22nd set}
\noindent {\em Only if INPUT(3) $\neq 0$}:
\begin{description}
\item[LABELI(I) I=1,2:] Label of the orbitals associated with the
requested integral.
\item[K:] Power of $r$.
\item[NEM:] Integral is:
\subitem $(P1*P2 + Q1*Q2)r{\sp K}$ \ \ \ \ \ if NEM $= 0$
\subitem $(P1*Q2 + P2*Q1)r{\sp K}$ \ \ \ \ \ if NEM $\neq 0$
\item[IALL:] If $>0$, ($<0$) all integrals (only the diagonal ones) are
calculated starting from the one defined by the first label.
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\subsubsection {\bf 23rd set}
\noindent {\em Only if INPUT(6) $\neq 0$}:
\begin{description}
\item[LABELI(I) I = 1 to INPUT(6):] Label of the orbitals to be written
on Fortran unit IPUN in the ascii format.
\end{description}
\centerline{\bf End of Input Data}
\noindent Two asterisks as the first 2 characters.
\newpage
\subsection{Input Format for PH92}
\centerline{\bf Files Used}
\begin{itemize}
\item Initial input file (initially Fortran unit number 5)
\item Optionally a second input file (defined by IREAD as the unit
number)
\item Optionally a third input file (defined by IBA as the unit number)
\item Standard output file (defined by IWRI)
\item One file to read wave functions in binary form (defined by IBB)
\end{itemize}
\centerline{\bf Input Data Provided Only Once at the Beginning}
\subsubsection {\bf 1st set}
\begin{description}
\item[idebug:] Minimum output if = 0, and long output otherwise. If
IDEBUG $> 9$, and IOUTPR in the 10th set is Y, last 1000 values of the
continuum wave functions are printed. (0)
\item[iread:] Fortran unit number to read the next sets of input data.
(5)
\item[iwri:] Fortran unit number for standard (printed) output. (6)
\item[ipun:] Not used. (7)
\item[iba:] Fortran unit number for Slater coefficients. (8)
\item[ibb:] Fortran unit number to read wave functions in binary form.
(9)
\item[ibc:] Not used. (9)
\item[nblipa:] Number of output lines per page. (57)
\end{description}
\subsubsection {\bf 2nd set}
\begin{description}
\item[ttot:] If positive, overall cpu time for the job in seconds. (0.)
\end{description}
The following data sets may be in different files depending on the
values of {\bf IREAD, IBA} read in the 1st and 2nd sets above.
\vspace{0.2in}
\centerline {\bf Input Data for Each Case to Be Calculated}
\subsubsection {\bf 1st set}
\begin{description}
\item [ITITLE:] Title line with one asterisk as its first character and
a maximum of 79 characters. The 3-digit number in columns 2--4 is used
to match the Slater coefficients in the Fortran unit IBA.
\end{description}
\subsubsection {\bf 2nd set}
\begin{description}
\item[KEYWORD :] Keyword for the type of expectation values to
calculate. Options are:
\subitem {GENERAL} for off-diagonal mean values of $\langle i \mid r
\mid j \rangle$.
\subitem {ELECTRIC} for electric multipole transition
\subitem {MAGNETIC} for magnetic multipole transition
\subitem {BORN} for the PWB cross section for electron-impact excitation
\subitem {PHOTOIONIZATION} for photoionization
\subitem Only the first 4 characters are meaningful
\end{description}
\subsubsection {\bf 3rd set}
\begin{description}
\item [LIST:] List of initial-state configurations in jj coupling (e.g.,
upper state for the $En$ and $Mn$ transitions, lower state for the Born
excitation cross sections).
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\subsubsection {\bf 4th set}
\begin{description}
\item [NRECI:] Record number to read the initial-state orbitals from a
binary file.
\end{description}
\subsubsection {\bf 5th set}
\begin{description}
\item [LIST:] List of final-state configurations in jj coupling (e.g.,
lower state for the $En$ and $Mn$ transitions, upper state for the Born
excitation cross sections).
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\subsubsection {\bf 6th set}
\begin{description}
\item [NRECF:] Record number to read the final-state orbitals from a
binary file.
\end{description}
\noindent{\it Only for PHOTOIONIZATION, the 7th through the 17th set}:
\subsubsection{\bf 7th set}
\begin{description}
\item [LABEL:] Label for ejected electron partial waves, only the first
3 characters will be used.
\end{description}
\subsubsection{\bf 8th set}
\begin{description}
\item [EVFREE:] Ejected photoelectron energy in eV.
\end{description}
\subsubsection{\bf 9th set}
\begin{description}
\item [AUGER:] Should be N to compute photoionization cross section.
\end{description}
\subsubsection{\bf 10th set}
\begin{description}
\item [IOUTPR:] If Y, tabulate partial waves.
\end{description}
\noindent {\it Only if IOUTPR is Y}:
\begin{description}
\item [IPSTEP:] Tabulate every IPSTEP points. (20)
\end{description}
\subsubsection{\bf 11th set}
\begin{description}
\item [MODFREE:] If Y, modify the methods for the partial waves.
\end{description}
\noindent {\it Only if MODFREE is Y}:
\begin{description}
\item [IOPTLA:] If Y, include off-diagonal Lagrange multipliers. (Y)
\item [IOPEXD:] If Y, include exchange interaction between bound and
free electrons. (Y)
\item [IORTH:] If Y, Schmidt orthogonalize before solving the Dirac-Fock
equation. (Y)
\end{description}
\subsubsection{\bf 12th set}
\begin{description}
\item [SCFCRIT:] If Y, convergence acceleration is used.
\end{description}
\noindent {\it Only if SCFCRIT is Y}:
\begin{description}
\item [COFSTO:] Initial value of mixing between two iterations (kept
frozen if no convergence acceleration). (0.9)
\item [CTEST:] SCF accuracy to be satisfied. ($10^{-5}$)
\item [ITERMA:] Maximum number of SCF iterations. (600)
\end{description}
\subsubsection{\bf 13th set}
\begin{description}
\item [WKB:] If Y, modify parameters used in the normalization of the
asymptotic Dirac-Fock solution using the WKB approximation.
\end{description}
\noindent {\it Only if WKB is Y}:
\begin{description}
\item [CTNOR:] Accuracy for partial wave normalization. ($3\times 10^{-
6}$)
\item [MAXNOR:] Maximum number of extensions (by 240 points) to reach
the asymptotic form of partial waves. (100)
\end{description}
\subsubsection{\bf 14th set}
\begin{description}
\item [NBFREE:] Number of continuum partial waves.
\end{description}
\subsubsection{\bf 15th set}
\begin{description}
\item [KFREE(I), I=1 to NBFREE:] $\kappa$ quantum number of partial
waves.
\end{description}
\subsubsection {\bf 16th set}
\begin{description}
\item [ITITLE:] Title line with one asterisk as its first character and
a maximum of 79 characters. The 3-digit number in columns 2--4 is used
to match the Slater coefficients in the Fortran unit IBA.
\end{description}
\subsubsection{\bf 17th set}
\noindent{\it For each transition matrix element, a line with (=output
from MJ92, the first four characters must be PHOT)}:
\begin{description}
\item [A:] Angular coefficient.
\item [IC, JC:] Configuration indexes.
\item [K:] Multipolarity.
\item [LABELI, LABELF:] Labels of orbitals for the matrix element.
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\vspace{0.2in}
\noindent {\bf Only for ELECTRIC and MAGNETIC}:
\subsubsection{\bf 7th set}
\begin{description}
\item [EVI:] Energy of the initial (=upper) state in eV.
\item [EVF:] Energy of the final (=lower) state in eV.
\end{description}
\subsubsection{\bf 8th set}
\begin{description}
\item [JJI:] 2$\times J$ of the initial state.
\item [JJF:] 2$\times J$ of the final state.
\end{description}
\subsubsection{\bf 9th set}
\begin{description}
\item [IPOL:] Multipolarity of the transition.
\end{description}
\subsubsection {\bf 10th set}
\begin{description}
\item [ITITLE:] Title line with one asterisk as its first character and
a maximum of 79 characters. The 3-digit number in columns 2--4 is used
to match the Slater coefficients in the Fortran unit IBA.
\end{description}
\subsubsection{\bf 11th set}
\noindent{\it For each transition matrix element, a line with (=output
from MJ92, the first four characters must be `` En '' for ELECTRIC and
`` Mn '' for MAGNETIC multipole transitions)}:
\begin{description}
\item [A:] Angular coefficient.
\item [IC, JC:] Configuration indexes.
\item [K:] Multipolarity.
\item [LABELI, LABELF:] Labels of orbitals for the matrix element.
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\vspace{0.2in}
\noindent {\bf Only for BORN}:
\subsubsection{\bf 7th set}
\begin{description}
\item [JJI:] 2$\times J$ of the initial (=lower) state.
\item [JJF:] 2$\times J$ of the final (=upper) state.
\end{description}
\subsubsection{\bf 8th set}
\begin{description}
\item [EEXC:] Excitation energy in eV. {\bf Colon (:) required after
EEXC.}
\end{description}
\subsubsection{\bf 9th set}
\begin{description}
\item [TEV:] List of incident energy in eV. The code accepts a total of
100 incident energies separated by a blank. Be careful not to use
columns 1--4.
\end{description}
\subsubsection {\bf 10th set}
\begin{description}
\item [ITITLE:] Title line with one asterisk as its first character and
a maximum of 79 characters. The 3-digit number in columns 2--4 is used
to match the Slater coefficients in the Fortran unit IBA.
\end{description}
\subsubsection{\bf 11th set}
\noindent{\it For each transition matrix element, a line with (=output
from MJ92, the first four characters must be BORN)}:
\begin{description}
\item [A:] Angular coefficient.
\item [IC, JC:] Configuration indexes.
\item [K:] Multipolarity.
\item [LABELI, LABELF:] Labels of orbitals for the matrix element.
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\vspace{0.2in}
\centerline{\bf End of Input Data}
\noindent Two asterisks as the first 2 characters. This indicates the
end of all input data.
\newpage
\subsection{Input format for DW92}
\centerline{\bf Files Used}
\begin{itemize}
\item Initial input file (initially Fortran unit number 5)
\item Optionally a second input file (defined by IREAD as the unit
number)
\item Optionally a third input file (defined by IBA as the unit number)
\item Standard output file (defined by IWRI)
\item One file to read wave functions in binary form
\item Three scratch files (Fortran units IBC=IBB, IBC+1, IBC+2)
\end{itemize}
\centerline{\bf Input Data Provided Only Once at the Beginning}
\subsubsection {\bf 1st set}
\begin{description}
\item[idebug:] Minimum output if = 0, and long output otherwise. If
IDBUG $> 9$, and IOUTPR in the 10th set is Y, last 1000 values of the
continuum wave functions are printed. (0)
\item [iread:] Fortran unit number to read the next sets of input data.
(5)
\item [iwri:] Fortran unit number for standard (printed) output. (6)
\item [ipun:] Not used. (7)
\item [iba:] Fortran unit number for Slater coefficients. (8)
\item [ibb:] Fortran unit number to read wave functions in binary form.
(9)
\item [nblipa:] Number of output lines per page. (57)
\item [npow:] Power of ten by which to multiply the speed of light. (0)
\end{description}
\subsubsection {\bf 2nd set}
\begin{description}
\item[ttot:] If positive, overall cpu time for the job in seconds. (0.)
\end{description}
The following data sets may be in different files depending on the
values of {\bf IREAD, IBA} read in the 1st and 2nd sets above.
\vspace{0.2in}
{\bf DW92 calculates cross section for one incident energy. New input
data must be submitted for each incident energy.}
\subsubsection {\bf 3rd set}
\begin{description}
\item [ITITLE:] Title line with one asterisk as its first character and
a maximum of 79 characters. The 3-digit number in columns 2--4 is used
to match the Slater coefficients in the Fortran unit IBA.
\end{description}
\subsubsection {\bf 4th set}
\noindent Debugging output requests:
\begin{description}
\item [IDEBUG1:] If $\neq 0$, stop execution after generating free
waves. (0)
\item [IDEBUG2:] If $\neq 0$, stop execution after interpolation between
the initial and final state radial meshes. (0)
\item [IDEBUG3:] If $\neq 0$, stop execution after calculating matrix
elements. (0)
\item [IDEBUG4:] If $\neq 0$, print cross section matrix elements. (0)
\item [IDEBUG5:] If $\neq 0$, maximum output during the calculation of
partial waves. (0)
\end{description}
\subsubsection {\bf 5th set}
\begin{description}
\item[KEYWORD :] Keyword for the type of cross section to calculate.
Options are:
\subitem {EXCITATION} for discrete excitations by electron impact.
\subitem {IONIZATION} for ionization by electron impact.
\subitem Only the first 4 characters are meaningful
\end{description}
\subsubsection {\bf 6th set}
\begin{description}
\item [LIST:] List of initial-state configurations in jj coupling.
\end{description}
\noindent A line with {\bf END} as the first 3 characters.
\subsubsection {\bf 7th set}
\begin{description}
\item [NRECI:] Record number to read the initial-state orbitals from a
binary file.
\item [JTARG:] 2$\times J$ of the target.
\item [LABEL:] Label for the partial waves. Only the first 3 characters
will be used.
\item [ENFEV:] Energy of the free electron in eV.
\end{description}
\subsubsection{\bf 8th set}
\begin{description}
\item [IOUTPR:] If Y, tabulate partial waves.
\end{description}
\noindent {\it Only if IOUTPR is Y}:
\begin{description}
\item [IPSTEP:] Tabulate every IPSTEP points. (20)
\end{description}
\subsubsection {\bf 9th set}
\noindent TYPE of partial waves to be calculated:
\begin{description}
\item [DW:] Distorted waves.
\item [CW:] Coulomb waves.
\item [PW:] Plane waves.
\end{description}
\noindent {\it If TYPE=CW, then read}:
\begin{description}
\item [ZEFF:] Effective nuclear charge.
\end{description}
\noindent {\bf The following 10th through 12th sets are only for DW}:
\subsubsection{\bf 10th set}
\begin{description}
\item [MODFREE:] If Y, modify the methods for the partial waves.
\end{description}
\noindent {\it Only if MODFREE is Y}:
\begin{description}
\item [IOPTLA:] If Y, include off-diagonal Lagrange multipliers. (Y)
\item [IOPEXD:] If Y, include exchange interaction between bound and
free electrons. (Y)
\item [IORTH:] If Y, Schmidt orthogonalize before solving the Dirac-Fock
equation. (Y)
\end{description}
\subsubsection{\bf 11th set}
\begin{description}
\item [SCFCRIT:] If Y, convergence acceleration is used.
\end{description}
\noindent {\it Only if SCFCRIT is Y}:
\begin{description}
\item [COFSTO:] Initial value of mixing between two iterations (kept
frozen if no convergence acceleration). (0.9)
\item [CTEST:] SCF accuracy to be satisfied. ($10^{-5}$)
\item [ITERMA:] Maximum number of SCF iterations. (600)
\end{description}
\subsubsection{\bf 12th set}
\begin{description}
\item [WKB:] If Y, modify parameters used in the normalization of the
asymptotic Dirac-Fock solution using the WKB approximation.
\end{description}
\noindent {\it Only if WKB is Y}:
\begin{description}
\item [CTNOR:] Accuracy for partial wave normalization. ($3\times 10^{-
6}$)
\item [MAXNOR:] Maximum number of extensions (by 240 points) to reach
the asymptotic form of partial waves. (100)
\end{description}
\subsubsection {\bf 13th set}
\begin{description}
\item [LQMIN:] The highest $l$ value of the partial waves for which
calculation has been completed already.
\item [LQMAX:] Maximum $l$ value of the partial waves to be calculated.
{\bf Warning: most computers cannot handle LQMAX $> 50$ because of the
smallest real number allowed (10$^{-308}$).}
\end{description}
\subsubsection {\bf 14th set}
\noindent {\it For IONIZATION only}:
\begin{description}
\item [LABEL:] Label for the scattered electron partial waves.
\item [ENFEV:] Scattered electron energy in eV.
\end{description}
\noindent {\bf For EXCITATION (IONIZATION), the sets 3 and 6 through 14
are for the incident electron (incident and scattered electron), and the
same sets must be repeated for the scattered (ejected) electron.}
\subsubsection {\bf 15th set}
\begin{description}
\item [ITITLE:] Title line with one asterisk as its first character and
a maximum of 79 characters. The 3-digit number in columns 2--4 is used
to match the Slater coefficients in the Fortran unit IBA.
\end{description}
\subsubsection {\bf 16th set}
\begin{description}
\item [BREIT:] If Y, Breit operator is included. (Appropriate Slater
coefficients from MJ92 must be available.)
\item [CONFAV:] If Y, configuration average is used. If N,
term-dependent partial waves are generated and used. In the latter
case, appropriate Slater coefficients from MJ92 must be available.
\end{description}
\centerline{\bf End of Input Data}
\noindent Two asterisks as the first 2 characters. This indicates the
end of all input data.
\vspace{0.2in}
Since the calculation of DWB cross sections involves a large number of
partial waves and hence Slater integrals, the Slater integrals are to be
read from a file separate from the input data listed above. The Slater
coefficients generated by the MJ92 code should be identified as the
Fortran unit 8. Then, the DW92 code reads the Slater coefficients and
use them in the calculation of the DWB cross sections.
\newpage
\subsection{Samples of Input Data for MJ92}
\begin{verbatim}
:
:
*881 Be-like, 2s2+2p2, total energy
:
ENERGY :
JT=0 :
CORE (1S)2 :
C 1 (2S)2 :
C 2 (2P)2 :
END
Y :
*882 Be-like, 2s2p+2p3d, J=1 total energy
:
ENERGY :
JT=2 :
CORE (1S)2 :
C 1 (2S)1 (2P)1 :
C 2 (2P)1 (3D)1 :
END
Y :
*883 Be-like, 2s2p+2p3d -> 2s2+2p2, E1
:
ELECTRIC :
MULTPOL=1
# upper state configurations first
JT=2 :
CORE (1S)2 :
C 1 (2S)1 (2P)1 :
C 2 (2P)1 (3D)1 :
END
JT=0 :
CORE (1S)2 :
C 1 (2S)2 :
C 2 (2P)2 :
END
*884 Be-like, 2s2p+2p3d <- 2s2+2p2, Born
:
BORN :
# initial state configurations first
JT=0 :
CORE (1S)2 :
C 1 (2S)2 :
C 2 (2P)2 :
END
JT=2 :
CORE (1S)2 :
C 1 (2S)1 (2P)1 :
C 2 (2P)1 (3D)1 :
END
*885 Be-like, 2s2p+2p3d <- 2s2+2p2, DW
:
EXCITATION :
0 :
C 1 (1S)2 (2S)2 :
C 2 (1S)2 (2P)2 :
END
LABL IN -1 20 :
2 :
C 1 (1S)2 (2S)1 (2P)1 :
C 2 (1S)2 (2P)1 (3D)1 :
END
LABL OUT -1 20 :
# N=no Breit; Y=configuration average
N Y :
**
\end{verbatim}
\subsection{Samples of Input Data for DF92, PH92, and DW92}
\begin{verbatim}
0 5 6 7 8 :
:
*881 C++, 2s2+2p2, ground state, MCDF
C 1 (1S)2 (2S)2 :
C 2 (1S)2 (2P*)2 :
C 3 (1S)2 (2P)2 :
END
NEIGV=1 ICMUL=0 IPRFGR=0 :
NZ=6 NOE=4 NORBSC=0
IORTH=0 NDEP=0 NLEC=0 NEC=1 :
NSTEP=0
N :
N :
N :
N :
NBRE=0 NRET=0 :
:
:
*882 C++, 2s2p+2p3d, 2nd J=1 MCDF
C 1 (1S)2 (2S)1 (2P*)1 :
C 2 (1S)2 (2S)1 (2P)1 :
C 3 (1S)2 (2P*)1 (3D*)1 :
C 4 (1S)2 (2P)1 (3D*)1 :
C 5 (1S)2 (2P)1 (3D)1 :
END
NEIGV=2 ICMUL=0 IPRFGR=0 :
NZ=6 NOE=4 NORBSC=0
IORTH=0 NDEP=0 NLEC=0 NEC=2 :
NSTEP=0
N :
N :
N :
N :
NBRE=0 NRET=0 :
:
:
**
0 5 6 7 8 :
:
# *881 Be-like, 2s2p+2p3d -> 2s2+2p2, E1
ELECTRIC :
C 1 (1S)2 (2S)1 (2P*)1 :
C 2 (1S)2 (2S)1 (2P)1 :
C 3 (1S)2 (2P*)1 (3D*)1 :
C 4 (1S)2 (2P)1 (3D*)1 :
C 5 (1S)2 (2P)1 (3D)1 :
END
NRECI=2
C 6 (1S)2 (2S)2 :
C 7 (1S)2 (2P*)2 :
C 8 (1S)2 (2P)2 :
END
NRECF=1
EVI=0.
EVF=12.69
JJI=2
JJF=0
IPOL=1
*881 Be-like, 2s2p+2p3d <- 2s2+2p2, Born
BORN :
C 1 (1S)2 (2S)2 :
C 2 (1S)2 (2P*)2 :
C 3 (1S)2 (2P)2 :
END
NRECI=1
C 4 (1S)2 (2S)1 (2P*)1 :
C 5 (1S)2 (2S)1 (2P)1 :
C 6 (1S)2 (2P*)1 (3D*)1 :
C 7 (1S)2 (2P)1 (3D*)1 :
C 8 (1S)2 (2P)1 (3D)1 :
END
NRECF=2
JJI=0
JJF=2
EEXC=12.69 :
TEV=13. 14. 15. 17. 20. 25. 30. 40. 50. :
**
0 5 6 7 8 :
:
*881 C++, 2s2p+2p3d <- 2s2+2p2, DW, T=15 eV
EXCITATION :
C 1 (1S)2 (2S)2 :
C 2 (1S)2 (2P*)2 :
C 3 (1S)2 (2P)2 :
END
FILE=1 JINI=0 LABL=IN ENGY=15. TABUL=N TYPE=DW
LGEXOR=N SCF=N WKB=N :
0 20 :
*882 C++, 2s2p+2p3d <- 2s2+2p2, DW, T'=15-12.69=3.31 eV
C 4 (1S)2 (2S)1 (2P*)1 :
C 5 (1S)2 (2S)1 (2P)1 :
C 6 (1S)2 (2P*)1 (3D*)1 :
C 7 (1S)2 (2P)1 (3D*)1 :
C 8 (1S)2 (2P)1 (3D)1 :
END
FILE=2 JFIN=2 LABL=OUT ENGY=3.31 TABUL=N TYPE=DW
LGEXOR=N SCF=N WKB=N :
0 20 :
*883 C++, 2s2p+2p3d <- 2s2+2p2, DW, lmax=20
# N=no Breit; Y=config. ave.
N Y :
**
\end{verbatim}
\acknowledgements
The computer codes described in this lecture note were developed
primarily by J. P. Desclaux with modifications by myself and P.
Indelicato. This lecture series was arranged with the financial support
from the Korea Atomic Energy Research Institute and the Asia-Pacific
Center for Theoretical Physics.
\begin{references}
\bibitem{bray}I. Bray and A. Stelbovics, Phys. Rev. A {\bf 46}, 6995
(1992); I. Bray and D. V. Fursa, Phys. Rev. A {\bf 54}, 2991 (1996).
\bibitem{mccurdy}M. Baertschy, T. N. Resigno, W. A. Issacs, and C. W.
McCurdy, Phys. Rev. A {\bf 60}, R13 (1999).
\bibitem{pindzola}F. Robicheaux, M. S. Pindzola, and D. R. Plante, Phys.
Rev. A {\bf 55}, 3573 (1977).
\bibitem{descl1}J.-P. Desclaux, Comput. Phys. Communic. {\bf 21}, 207
(1975). This is the first public documentation of his relativistic wave
function code.
\bibitem{kim94}Y.-K. Kim and M. E. Rudd, Phys. Rev. A (1994).
\bibitem{condon}E. U. Condon and G. H. Shortley, {\it The Theory of
Atomic Spectra} (Cambridge Univ. Press, London, 1959).
\bibitem{slater}J. C. Slater, {\it Quantum Theory of Atomic Structure}
(McGraw-Hill, New York, 1960), Vol. II.
\bibitem{fischer}C. F. Fischer, T. Brage, and P. J\"{o}nsson, {\it
Computational Atomic Structure} (Inst. Phys. Publishing, London, 1997);
See also C. F. Fischer, {\it The Hartree-Fock Method for Atoms} (Wiley
and Sons, New York, 1977).
\bibitem{herman}F. Herman and S. Skillman, {\it Atomic Structure
Calculations} (Prentice Hall, Englewood Cliffs, NJ, 1963).
\bibitem{kohn}W. Kohn and L. J. Sham, Phys. Rev. {\bf 140}, A1133
(1965).
\bibitem{cowan}R. D. Cowan, {\it The Theory and Atomic Structure and
Spectra} (Univ. Calif. Press, Berkeley, 1981).
\bibitem{bethe57}H. A. Bethe and E. E. Salpeter, {\it Quantum Mechanics
of One- and Two-Electron Atoms} (Plenum, New York, 1977), p. 251.
\bibitem{bethe30}H. Bethe, Ann. Physik {\bf 5}, 325 (1930).
\bibitem{inokuti}M. Inokuti, Rev. Mod. Phys. {\bf 43}, 297 (1971); M.
Inokuti, Y. Itikawa, and J. E. Turner, Rev. Mod. Phys. {\bf 50}, 23
(1978).
\bibitem{nistweb}See examples in {\it http://physics.nist.gov/ionxsec}.
\bibitem{kim00}Y.-K. Kim, W. R. Johnson, and M. E. Rudd, Phys. Rev. A
{\bf 61}, 034702 (2000).
\bibitem{weiss}A. W. Weiss and Y.-K. Kim, Phys. Rev. A {\bf 51}, 4487
(1995).
\bibitem{baik}Y.-K. Kim, D. H. Baik, P. Indelicato, and J. P. Desclaux,
Phys. Rev. A {\bf 44}, 148 (1991).
\end{references}
\begin{table}
\caption{Fundamental constants associated with atomic units.$^*$}
\begin{tabular}{ll}
electron mass $m$ & 9.109 381 88 (72) $\times 10^{-31}$\ kg \\
electron charge $e$ & -1.602 176 462 (63) $\times 10^{-19}$\ C \\
Planck's constant/2$\pi\ \hbar$ & 1.054 571 596 (82) $\times 10^{-34}$\
Js \\
speed of light $c$ & 299 792 458\ m/s (exact) \\
inverse fine-structure constant $1/\alpha = \hbar c/e^2$ & 137.035 999
76 (50) \\
Bohr radius $a_0=\hbar/me^2$ & 0.529 177 2083 (19) \AA \\
rydberg constant $R_{\infty}=\case{1}{2}mc^2\alpha^2$ & 109 737.315 685
49 (83) cm$^{-1}$ \\ & 13.605 691 72 (53) eV \\
electron volt $e/C$ in joules & 8065.544 77 (32) cm$^{-1}$
\end{tabular}
\label{tbl:table1}
\noindent $^*$ From P. J. Mohr and B. H. Taylor, J. Phys. Chem. Ref.
Data {\bf 28}, 1713 (1999). This reference contains the latest list of
recommended fundamental constants.
\end{table}
\begin{table}
\caption{Dirac quantum number $\kappa$.}
\begin{tabular}{cccccccccc}
$\kappa$ & $-$1 & 1 & $-$2 & 2 & $-$3 & 3 & $-$4 & 4 & $-$5 \\
$l$ and $j$ & s$_{1/2}$ & p$_{1/2}$ & p$_{3/2}$ & d$_{3/2}$ & d$_{5/2}$
& f$_{5/2}$ & f$_{7/2}$ & g$_{7/2}$ & g$_{9/2}$ \\
Input & s & p* & p & d* & d & f* & f & g* & g
\end{tabular}
\label{tbl:table2}
\end{table}
\begin{table}
\caption{Selection rules for $E1,\ M1$, and $E2$ transitions.$^*$}
\begin{tabular}{cccc}
& E1 & M1 & E2 \\ \tableline
(1) & $\Delta J=0,\ \pm 1$ & $\Delta J=0,\ \pm 1$ & $\Delta J=0,\ \pm
1,\ \pm 2$ \\
& (No $0 \leftrightarrow 0$) & (No $0 \leftrightarrow 0$) & (No $0
\leftrightarrow 0,\ \case{1}{2} \leftrightarrow \case{1}{2},\ \ 0
\leftrightarrow 1$) \\
\vspace{0.1in}
(2) & $\Delta M=0,\ \pm 1$ & $\Delta M=0,\ \pm 1$ & $\Delta M=0,\ \pm
1,\ \pm 2$ \\
\vspace{0.1in}
(3) & Parity change & No parity change & No parity change \\
\vspace{0.1in}
(4) & One electron jump & No electron jump & One or no electron jump \\
& $\Delta \ell = \pm 1$ & $\Delta \ell =0,\ \Delta n=0$ & $\Delta \ell
=0,\ \pm 2$ \\
\vspace{0.1in}
(5) & $\Delta S=0$ & $\Delta S=0$ & $\Delta S=0$ \\
\vspace{0.1in}
(6) & $\Delta L=0,\ \pm 1$ & $\Delta L=0$ & $\Delta L=0,\ \pm 1,\ \pm
2$ \\
& (No $0 \leftrightarrow 0$) & & (No $0\leftrightarrow 0, \
0\leftrightarrow 1$)
\label{tbl:table3}
\end{tabular}
\noindent $^*$ From R. H. Garstang, {\it Forbidden Transitions} in {\it
Atomic and Molecular Processes}, edited by D. R. Bates (Academic Press,
New York, 1962), p. 4.
\end{table}
\end{document}