An Introduction to Relativistic Theory as Implemented in GRASP

Abstract

Computational atomic physics continues to play a crucial role in both increasing the understanding of fundamental physics (e.g., quantum electrodynamics and correlation) and producing atomic data for interpreting observations from large-scale research facilities ranging from fusion reactors to high-power laser systems, space-based telescopes and isotope separators. A number of different computational methods, each with their own strengths and weaknesses, is available to meet these tasks. Here, we review the relativistic multiconfiguration method as it applies to the General Relativistic Atomic Structure Package [grasp2018, C. Froese Fischer, G. Gaigalas, P. Jönsson, J. Bieroń, Comput. Phys. Commun. (2018). DOI: 10.1016/j.cpc.2018.10.032]. To illustrate the capacity of the package, examples of calculations of relevance for nuclear physics and astrophysics are presented.

1. Introduction

Atomic structure theory has developed for over a century in synergy with computational science and techniques. Early contributions came from Douglas Hartree, who was inspired by Niels Bohr in the early 1920s to use his numerical skills on atomic systems. This resulted in his publications on “The Wave mechanics of an Atom with a Non-Coulomb Central Field” [1,2,3,4]. In these, he also proposed and used the self-consistent field (SCF) method to solve what would later be labeled the Hartree equations, where both the potentials and the wave functions are unknown. In 1930, Slater [5] and Fock [6] independently pointed out that the Hartree method did not take into account the required anti-symmetry of the wave functions under the exchange of two electrons and suggested ways to deal with this. This led Hartree [7] to develop a more practical method than the one previously proposed by Fock, which was later labeled Hartree–Fock. In these early days, the gross structures of small ions were the testing ground, leaving out correlation and relativistic effects. Inspired by this line of work, atomic systems became, and still are, fundamental testing grounds for quantum mechanics, many-body physics and fundamental processes.

As soon as the Hartree and Hartree–Fock methods were established, the interest turned to the corrections: the most important being relativistic effects and electron correlation. The theory of special relativity was incorporated into quantum mechanics by Dirac [8]. Dirac himself doubted its importance for atomic systems, which probably, together with the increased complexity involved, delayed the development of fully relativistic and computational methods. Swirles, however, published a paper on relativistic self-consistent fields [9], but they had to wait nearly 30 years before it became possible to attempt solutions.

The many-body instantaneous interactions introduce a correction to the Hartree–Fock approach, by keeping the electrons apart (in contrast to the independent particle approach) – a correlation of their motions [10]. To deal with this, the configuration interaction method was introduced, a term that appeared in 1933 [11,12], but the method had been developed earlier by, among others, Slater [13] and Condon [14]. The fully variational multiconfiguration approach, which is in focus in this paper, was first explored by Hartree, Hartree, and Swirles for O${}^{+}$ [15].

It is important to remember that the computational methods were very basic at these early stages of atomic calculations. On the machinery side, researchers used slide rules, hand calculators, and mechanical machines [16,17], while electronic computers started to become available in the late 1950s. At this time, numerical methods were readily available for the physics community, and when FORTRAN II appeared in 1958, simple calculations were possible for the first time without having to use machine code. The history of atomic structure calculations, leading to the methods described in this work, was reviewed by Hartree [18,19]. The history of configuration interaction methods is also discussed in a more recent review by Shavitt [20].

The introduction of electronic computers, accompanied by efficient numerical methods and program languages, led to the start of a new era in the late 1950s for the variational methods we will discuss in this paper, not the least through works by Hartree’s doctoral students David Mayers [21] and Charlotte Froese [22]. It is clear that Hartree continued to play an important role in the pioneering work on atomic calculations, which was recently described in a monograph [23]. The development of the non-relativistic Hartree–Fock theory for the newly introduced computers was led by Froese [24] and later extended to the multiconfiguration numerical approach [25,26]. These methods and the associated computer programs [27,28] are described in two monographs [29,30]. At the same time, the relativistic methods had a renaissance through the pioneering work of Grant [31,32], which was described in the review by the same author [33]. An effort in the early 1970s by Desclaux, Mayers, and Grant [34] led to the first multiconfiguration Dirac–Hartree–Fock (MCDHF) programs based on Slater determinants [35,36] and later on coupled configuration state functions (CSFs) as eigenfunctions to the relevant angular momenta [37] with the appropriate angular codes [38,39].

The structures of these earlier codes reflect the hardware and software then available. Since then, tremendous progress has been made in the development of the computational tools based on the MCDHF method [37,40,41,42,43] owing to improving software utilities, vastly increasing speed, memory, and overall power of computer facilities that gave the General-purpose Relativistic Atomic Structure Package (grasp) the flexibility of today’s mark. Ten years ago, the Computational Atomic Structure (CompAS) [44] group was formed. Its focus is to develop methodologies that utilize modern computational technology to meet the demands for accurate atomic data. Examples are a series of computations for transition energies that achieve spectroscopic accuracy, i.e., the accuracy that is similar to what can be obtained in experiments and that can be used to unambiguously assign observed lines in plasmas to the correct transitions [45,46,47,48,49]. Other notable achievements are the series of computations of energies and transition rates by Gaigalas and co-workers for the complex Lanthanides for opacity modeling in relation to kilonova [50,51,52]. The CompAS group is also using this opportunity to treat new challenges and applications opened up by increased computation power paired with a deeper understanding of correlation and other effects in atomic systems. Some examples of recent contributions to the front-line of atomic physics are hyperfine-quenched lines to determine isotope abundances in planetary nebula [53], magnetically induced transitions to measure field strengths in the solar corona [54,55,56], sympathetic cooling by negative ions [57], tests of models for quantum electrodynamics (QED) in complex systems with quenched correlation [58,59,60], electronic parameters relevant in nuclear physics [61,62,63], electron affinities and negative ions properties [64,65]. At the core of all activities, the CompAS group is developing state-of-the-art computer codes and methodologies for atomic structure calculations [44].

Here, we review the development of the relativistic multiconfiguration method as it applies to grasp [42,43], a suite of Fortran95 codes, adapted to run under MPI, that compute approximate wave functions of atomic states and evaluate a large number of atomic properties, such as transition energies and rates, diagonal and off-diagonal hyperfine structures [66], Landé factors and magnetic interactions [67], isotope shifts and detailed interactions with extended and deformed nuclei [68], electron densities [69], and atomic state labels in various coupling schemes [70,71].

The present paper provides the necessary theoretical support for the grasp2018 manual [72], which aims to guide users in the practical handling of the above codes. A detailed account of multiconfiguration methods and relativistic atomic physics can be found in [30,73]. To illustrate the capacity of the grasp suite of codes, we look at a few applications in nuclear physics and astrophysics. Atomic units [74] are adopted throughout.

2. Grasp Theory—Wave Functions

2.1. One-Electron Dirac Orbital Functions

For a particle in a spherically symmetric scalar potential $q\mathsf{\Phi }=V\left(r\right)$, the wave functions of the Dirac equation

$$i\hslash \frac{\partial \mathsf{\Psi }(\mathbf{r},t)}{\partial t}=\left\{c\mathit{\alpha }·\mathit{p}+c^2(\beta -I)+V\left(r\right)\right\}\mathsf{\Psi }(\mathbf{r},t)$$

(1)

have the form

$$\mathsf{\Psi }(\mathbf{r},t)=e^{-iEt/\hslash }\psi \left(\mathbf{r}\right),$$

(2)

where E is an eigenvalue of the Dirac Hamiltonian

$$\mathcal{H}\psi \left(\mathbf{r}\right)=\left\{c\mathit{\alpha }·\mathit{p}+c^2(\beta -I)+V\left(r\right)\right\}\psi \left(\mathbf{r}\right)=E\psi \left(\mathbf{r}\right).$$

(3)

Here, $\mathit{\alpha }$ and $\beta$ are the usual 4 × 4 Dirac matrices, $\mathbf{p}\equiv -i\nabla$ is the electron momentum operator and c is the speed of light (=$1/\mathit{\alpha }=137.035999084\left(21\right)a_0{E}_{\mathrm{h}}/\hslash ).$ The eigenfunctions $\psi \left(\mathbf{r}\right)$ are called orbitals. The one-electron Hamiltonian $\mathcal{H}$ commutes with angular momentum $\mathbf{j}=\mathbf{l}+\mathbf{s}$, and thus, the orbitals can be taken as eigenfunctions of ${\mathbf{j}}^2$ and $j_z$ with $jm$ as resulting good quantum numbers. Introducing the quantum number $\kappa$

$$\kappa =\left\{\begin{array}{cc}-(l+1)& \mathrm{for}j=l+1/2\hfill \\ +l& \mathrm{for}j=l-1/2\hfill \end{array}\right.,$$

(4)

the orbitals can be written in spherical co-ordinates as [73]

$${\psi }_{n\kappa m}(r,\theta ,\phi )=\frac{1}{r}\left(\begin{array}{c}{P}_{n\kappa }\left(r\right){\Omega }_{\kappa m}(\theta ,\phi )\\ i{Q}_{n\kappa }\left(r\right){\Omega }_{-\kappa m}(\theta ,\phi )\end{array}\right),$$

(5)

where $P_{n\kappa }\left(r\right)$ and $Q_{n\kappa }\left(r\right)$ are the radial functions and ${\Omega }_{\kappa m}(\theta ,\phi )$ are two-component spherical spinors built from the coupling of the spherical harmonics $Y_{l{m}_l}(\theta ,\phi )$ and the spin functions ${\chi }_{m_s}^{(1/2)}$. For these one-electron orbitals, the grasp code adopts the usual spectroscopic notation

$$1s,2s,2p-,2p,3s,3p-,3p,3d-,3d,\dots ,4f-,4f,5g-,\dots$$

where the relations (4) between l, j and $\kappa$ are exemplified for $l\le 4$ in

Table 1. Spectroscopic notation of relativistic shells.

	${\mathit{s}}_{1/2}$	${\mathit{p}}_{1/2}$	${\mathit{p}}_{3/2}$	${\mathit{d}}_{3/2}$	${\mathit{d}}_{5/2}$	${\mathit{f}}_{5/2}$	${\mathit{f}}_{7/2}$	${\mathit{g}}_{7/2}$	${\mathit{g}}_{9/2}$
	s	p-	p	d-	d	f-	f	g-	g
l	0	1	1	2	2	3	3	4	4
j	$1/2$	$1/2$	$3/2$	$3/2$	$5/2$	$5/2$	$7/2$	$7/2$	$9/2$
$\kappa$	$-1$	$+1$	$-2$	$+2$	$-3$	$+3$	$-4$	$+4$	$-5$

.

As it can be seen from this table, the $n\kappa m$ label of the relativistic orbital (5) is equivalent to $nljm$. In the non-relativistic limit ($c\to \infty$), $P_{n\kappa }\to P_{nl}$ and $Q_{n\kappa }\to 0$, respectively.

2.2. Atomic State Function and Configuration State Functions

Relativistic atomic structure theory in grasp treats the many-electron atom within Furry’s bound state interaction picture [75] of quantum electrodynamics. The formalism resembles that of non-relativistic atomic theory but is based on a Hamiltonian that uses the Dirac operator $c\mathit{\alpha }·\mathbf{p}+c^2(\beta -I)$ to describe the relativistic dynamics of electrons in the field of other charged particles, as described in Section 2.3 below.

The program structure draws from the one of the non-relativistic program atsp2K [28]. grasp assumes that the wave function of an atomic state $\mathsf{\Gamma }J{M}_J\pi$, with $\mathsf{\Gamma }$ being its identifying label, J the total angular momentum quantum number, $M_J$ the total magnetic quantum number and $\pi$ its parity, is approximated by an atomic state function (ASF), $\mathsf{\Psi }\left(\mathsf{\Gamma }J{M}_J\pi \right)$, which is a linear combination of configuration state functions (CSFs) in $jj$-coupling

$$\mathsf{\Psi }\left(\mathsf{\Gamma }J{M}_J\pi \right)\equiv |\mathsf{\Gamma }J{M}_J\pi \rangle =\sum _{\mathit{\alpha }=1}^{N_{CSF}}{c}_{\mathit{\alpha }}^{\mathsf{\Gamma }J}\mathsf{\Phi }\left({\gamma }_{\mathit{\alpha }}J{M}_J\pi \right).$$

(6)

For each CSF $\mathsf{\Phi }\left({\gamma }_{\mathit{\alpha }}J{M}_J\pi \right)$, the multi-index label ${\gamma }_{\mathit{\alpha }}$ contains all the needed information on its structure, i.e., the constituent subshells with their symmetry labels and the way their angular momenta are coupled to each other, as described in Section 2.4. The main difference with the non-relativistic program is the need to use four-component spinor orbital functions (5) to build the configurational state functions (CSF) with accompanying technical machinery to express the variational equations in computable form.

grasp can perform both variational MCDHF and configuration–interaction (CI) calculations. The latter assume that all data on the CSFs, radial orbital functions and angular coupling coefficients have previously been computed. The CI approach is therefore limited to the eigenvalue problem

$$\mathbf{H}\mathbf{c}=E\mathbf{c},$$

(7)

where $\mathbf{H}$ is the matrix of one of the chosen Hamiltonians of Section 2.3 in the CSF space spanned by a given set $\left\{\mathsf{\Phi }\left({\gamma }_{\mathit{\alpha }}J{M}_J\pi \right)\right\}$ and $\mathbf{c}={(c_1,\dots ,c_{N_{CSF}})}^t$ is the eigenvector corresponding to the eigenvalue E, the total energy. Note that the magnetic quantum number $M_J$ is irrelevant for field-free atoms, as assumed for the MCDHF and CI approaches considered in the present section. The mixing coefficients $\left\{c_{\mathit{\alpha }}^{\mathsf{\Gamma }J}\right\}$ appearing in (6) are therefore a priori $M_J$-independent. However, the electronic magnetic quantum number $M_J$ might become relevant and affect the electronic structure when external perturbations and/or hyperfine effects due to the coupling of the electronic $\mathbf{J}$ and nuclear $\mathbf{I}$ angular momenta are considered, as described in Section 3.2.

The MCDHF approach determines the orbital components $P_{n\kappa }\left(r\right),Q_{n\kappa }\left(r\right)$ of the one-electron orbitals (5) spanning the CSF space, permitting the numerical construction of the CSFs. The radial functions $P_{n\kappa }\left(r\right)$ and $Q_{n\kappa }\left(r\right)$ are defined on a grid

$$r_i=\frac{A}{Z}(exp\left[B(i-1)\right]-1),i=1,\dots ,NP,$$

(8)

where $A=2\times {10}^{-6}$, B = 0.05, and $NP=590$ are parameters that define the default grid. The default grid is usually sufficient for light to medium heavy systems. For the heaviest systems, it may be necessary to use a few thousand points, as is discussed in the accompanying grasp manual §13.5. The radial functions are obtained by solving a set of differential equations (see Section 2.7) with finite difference methods [26,73]. Just as the wave functions of a one-electron equation, the orbitals in grasp make up an orthonormal set.

2.3. Dirac-Coulomb, Dirac-Coulomb-Breit Hamiltonians and QED Corrections

In grasp, the Hamiltonian used for self-consistent calculations, see Section 2.7, is the Dirac–Coulomb Hamiltonian

$${\mathcal{H}}_{DC}=\sum _{i=1}^N\left(c{\mathit{\alpha }}_i·{\mathbf{p}}_i+c^2({\beta }_i-I)+V_{nuc}\left(r_i\right)\right)+\sum _{j>i=1}^N\frac{1}{r_{ij}}.$$

(9)

The nuclear potential $V_{nuc}\left(r\right)$ results from a nuclear charge density given by a two-parameter Fermi distribution function [76]. The nuclear charge, the two parameters of the Fermi distribution, the mass of the nucleus, the nuclear spin and the magnetic dipole and quadrupole moments are, in the grasp code suite, defined by the program rnucleus and saved in a file isodata (see accompanying manual §2.2 and §8.1).

Corrections to the above Hamiltonian can be included in configuration interaction calculations; see Section 2.8. To represent magnetic interactions and retardation effects [77,78,79], one can add the so-called transverse photon interaction, which is correct to order of ${\alpha }^2$.

$$\begin{array}{ccc}\hfill {\mathcal{H}}_{TP}=& & -\sum _{j>i=1}^N\left[\frac{{\mathit{\alpha }}_i·{\mathit{\alpha }}_{j}cos({\omega }_{ij}{r}_{ij}/c)}{r_{ij}}\right.\hfill \\ & & +\left.({\mathit{\alpha }}_i·{\nabla }_i)({\mathit{\alpha }}_j·{\nabla }_j)\frac{cos({\omega }_{ij}{r}_{ij}/c)-1}{{\omega }_{ij}^2{r}_{ij}/c^2}\right].\hfill \end{array}$$

(10)

In this expression, the ∇-operators act only on $r_{ij}$, while ${\omega }_{ij}$ represents the energy of the virtual photon exchanged between two electrons, as described in QED [78]. In the low-photon energy limit, when ${\omega }_{ij}\to 0$, the expression (10) reduces to the Breit interaction

$${\mathcal{H}}_{Breit}=-\sum _{j>i=1}^N\frac{1}{2r_{ij}}\left[{\mathit{\alpha }}_i·{\mathit{\alpha }}_j+\frac{\left({\mathit{\alpha }}_i·{\mathbf{r}}_{ij}\right)\left({\mathit{\alpha }}_j·{\mathbf{r}}_{ij}\right)}{r_{ij}^2}\right].$$

(11)

Adding the transverse photon interaction to the Dirac–Coulomb Hamiltonian gives the Dirac–Coulomb–Breit Hamiltonian

$${\mathcal{H}}_{DCB}={\mathcal{H}}_{DC}+{\mathcal{H}}_{TP}\simeq {\mathcal{H}}_{DC}+{\mathcal{H}}_{Breit}.$$

(12)

Additional important QED contributions are the self-energy (SE) correction and the vacuum polarization (VP). The SE correction, which is the result of emission and absorption of a virtual photon by the same electron, is given as a sum of one-electron corrections weighted by the fractional occupation numbers of the one-electron orbitals in the wave function. The VP correction, which is related to the creation and annihilation of virtual electron–positron pairs in the field of the nucleus, can be described by a correction to the Coulomb potential. The above QED terms [80,81,82,83] are included in the configuration interaction rci code, as implemented in the original grasp program [78], to yield the final Hamiltonian

$${\mathcal{H}}_{DCB+QED}={\mathcal{H}}_{DCB}+{\mathcal{H}}_{SE}+{\mathcal{H}}_{VP}.$$

(13)

Nuclear recoil effects [68], hyperfine interactions [66], and symmetry breaking interactions with external magnetic fields [84] are included in first-order perturbation theory; see Section 3.

2.4. Building Configuration State Functions

Adopting the following notation for a subshell containing w equivalent electrons

$${\left(nlj\right)}^w={\left(n\kappa \right)}^w,$$

a general relativistic configuration consists of m groups of equivalent electrons

$${\left(n_1{\kappa }_1\right)}^{w_1}{\left(n_2{\kappa }_2\right)}^{w_2}...{\left(n_m{\kappa }_m\right)}^{w_m},N=\sum _{i=1}^m{w}_i,$$

(14)

where $w_i$ is the occupation number of the relativistic subshell i [30].

The relativistic Hamiltonian for an N-electron system commutes with the total angular momentum operator $\mathbf{J}={\mathbf{j}}_1+\dots +{\mathbf{j}}_N$, and the solutions to the wave equation can be taken as eigenfunctions of ${\mathbf{J}}^2$ and $J_z$ with $J{M}_J$ as good quantum numbers. Atoms are fermionic systems, and wave functions are required to be antisymmetric with respect to permutations of the electron co-ordinates. Not all states arising from the angular momentum coupling $\mathbf{J}={\sum }_{i=1}^N{\mathbf{j}}_i$ are permitted, because the Pauli exclusion principle selects only fermionic states within each subshell of equivalent electrons. A configuration state function (CSF) is the simplest approximation of a many-electron wave function being both antisymmetric and an eigenfunction of ${\mathbf{J}}^2$ and $J_z$. Such a CSF can be built

• By using the well-known vector coupling techniques of angular momentum theory [85] to couple sequentially, from left to right, the subshell angular momenta $\mathbf{J}={\sum }_{i=1}^m{\mathbf{J}}_i$,

  $$(\left((J_1,J_2)J_{12}\right),j_3)J_{123},\dots ,J_{12\dots m-1},J_m)J{M}_J,$$

  associated with the m antisymmetric subshell wave functions

  $$|{\left(n_i{\kappa }_i\right)}^{w_i}{\mathit{\alpha }}_i{\nu }_i{J}_i{M}_i\rangle ,$$

  (15)
• Antisymmetrize the resulting coupled products through the permutations restricted to the exchange of electron coordinates involving different subshells [86].

In each wave function (15), $\nu$ is the seniority number [87,88], making it possible to discriminate states arising from the same relativistic subshell ${\left[j\right]}^w$ configuration having the same J-value. As shown in

Table 2. Subshell states ${\left[j\right]}^w$ are listed for $j=1/2,\dots ,9/2$, in both seniority $\left(\nu \right)$ and quasi-spin $\left(2Q\right)$ representations. An extra number $Nr$ is introduced for $j=9/2$.

$\mathit{subshell}$	$\mathit{\nu }$	J	$2\mathit{Q}$	$\mathit{Nr}$	$\mathit{subshell}$	$\mathit{\nu }$	J	$2\mathit{Q}$	$\mathit{Nr}$
${\left[1/2\right]}^0$ or ${\left[1/2\right]}^2$	0	0	1			3	$5/2$	2
${\left[1/2\right]}^1$	1	$1/2$	0			3	$7/2$	2
	3	$9/2$	2
${\left[3/2\right]}^0$ or ${\left[3/2\right]}^4$	0	0	2			3	$11/2$	2
${\left[3/2\right]}^1$ or ${\left[3/2\right]}^3$	1	$3/2$	1			3	$13/2$	2
${\left[3/2\right]}^2$	0	0	2			3	$15/2$	2
	2	2	0			3	$17/2$	2
	3	$21/2$	2
${\left[5/2\right]}^0$ or ${\left[5/2\right]}^6$	0	0	3		${\left[9/2\right]}^4$ or ${\left[9/2\right]}^6$	0	0	5
${\left[5/2\right]}^1$ or ${\left[5/2\right]}^5$	1	$5/2$	2			2	2	3
${\left[5/2\right]}^2$ or ${\left[5/2\right]}^4$	0	0	3			2	4	3
	2	2	1			2	6	3
	2	4	1			2	8	3
${\left[5/2\right]}^3$	1	$5/2$	2			4	0	1
	3	$3/2$	0			4	2	1
	3	$9/2$	0			4	3	1
	4	4	1	1
${\left[7/2\right]}^0$ or ${\left[7/2\right]}^8$	0	0	4			4	4	1	2
${\left[7/2\right]}^1$ or ${\left[7/2\right]}^7$	1	$7/2$	3			4	5	1
${\left[7/2\right]}^2$ or ${\left[7/2\right]}^6$	0	0	4			4	6	1	1
	2	2	2			4	6	1	2
	2	4	2			4	7	1
	2	6	2			4	8	1
${\left[7/2\right]}^3$ or ${\left[7/2\right]}^5$	1	$7/2$	3			4	9	1
	3	$3/2$	1			4	10	1
	3	$5/2$	1			4	12	1
	3	$9/2$	1		${\left[9/2\right]}^5$	1	$9/2$	4
	3	$11/2$	1			3	$3/2$	2
	3	$15/2$	1			3	$5/2$	2
${\left[7/2\right]}^4$	0	0	4			3	$7/2$	2
	2	2	2			3	$9/2$	2
	2	4	2			3	$11/2$	2
	2	6	2			3	$13/2$	2
	4	2	0			3	$15/2$	2
	4	4	0			3	$17/2$	2
	4	5	0			3	$21/2$	2
	4	8	0			5	$1/2$	0
	5	$5/2$	0
${\left[9/2\right]}^0$ or ${\left[9/2\right]}^{10}$	0	0	5			5	$7/2$	0
${\left[9/2\right]}^1$ or ${\left[9/2\right]}^9$	1	$9/2$	4			5	$9/2$	0
${\left[9/2\right]}^2$ or ${\left[9/2\right]}^8$	0	0	5			5	$11/2$	0
	2	2	3			5	$13/2$	0
	2	4	3			5	$15/2$	0
	2	6	3			5	$17/2$	0
	2	8	3			5	$19/2$	0
${\left[9/2\right]}^3$ or ${\left[9/2\right]}^7$	1	$9/2$	4			5	$25/2$	0
	3	$3/2$	2

extracted from [89], the seniority $\nu =2$ and $\nu =6$ are needed to discriminate the two states $J=2$ (or $J=4$) arising from ${[7/2]}^4$, i.e., from any subshell ${\left(nf\right)}^4$ or ${\left(ng-\right)}^4$. If the seniority classification fails, i.e., if two levels appear with the same $\left(\nu J\right)$ values, additional quantum numbers $\alpha$ are usually introduced to unambiguously designate the state considered [90]. A systematic approach to classify the states uses the theory of Lie groups [90] in which $\alpha$ contains the irreducible representations labels of the invoked chain of groups. Unfortunately, a complete classification scheme remains an open problem [91]. Therefore, from a practical point of view, we will restrict $\alpha$ to a simple number, $Nr$, as suggested in [89,92]. As revealed by Table 2—see also Table A.5 of [73]—this case arises for configurations ${[9/2]}^4$ or ${[9/2]}^6$ for which the same seniority, $\nu =4$, is assigned to the two allowed $J=4$ (or $J=6$) values.

Due to the current restrictions of the grasp2018 package “occupied subshells with $j\ge 9/2$ are restricted to a maximum of two electrons”, $Nr$ will never be seen in the output file rcsf.out produced by the rcsfgenerate program; see accompanying manual §3.2. However, we deliberately keep the $Nr$ notation, since the spin-angular library [93] allows us to deal with the $j=9/2$ subshell without any occupation restriction.

The above two-step procedure leads to the most general form of a CSF [30]

$$\begin{array}{ccc}\hfill & & \mathsf{\Phi }\left(\gamma J{M}_J\pi \right)\equiv |\gamma J{M}_J\pi \rangle \hfill \\ & & =|{\left(n_1{\kappa }_1\right)}^{w_1}{\mathit{\alpha }}_1{\nu }_1{J}_1{\left(n_2{\kappa }_2\right)}^{w_2}{\mathit{\alpha }}_2{\nu }_2{J}_2{J}_{12}{\left(n_3{\kappa }_3\right)}^{w_3}{\mathit{\alpha }}_3{\nu }_3{J}_3{J}_{123}\dots {\left(n_m{\kappa }_m\right)}^{w_m}{\mathit{\alpha }}_m{\nu }_m{J}_{m}J{M}_J\pi \rangle .\hfill \end{array}$$

(16)

In this notation, $\gamma$ is the compact representation of a CSF, collecting all the needed information for each subshell together with all the intermediate quantum numbers that unambiguously define the CSF. The parity $\pi$ is often omitted in the CSF notation since it can be easily deduced from the collection of angular momenta hidden in $\gamma$, using $\pi ={(-1)}^{{\sum }_{i=1}^N{l}_i}$.

Each subshell wave function (15) can be built using a recursive coupling method in terms of fractional parentage coefficients (CFPs) [30,73,86]. In addition to the seniority classification, there exists an alternative representation of the same subshell wave function provided by the quasi-spin formalism [91,94,95] that offers many advantages explored in the spin-angular algebra [93,96]. In the quasi-spin classification scheme, also used in Table 2, the subshell wave function (15) is rewritten as

$$|{\left(n_i{\kappa }_i\right)}^{w_i}{\mathit{\alpha }}_i{\nu }_i{J}_i{M}_i\rangle =|{\left(n_i{\kappa }_i\right)}^{w_i}{\mathit{\alpha }}_i{Q}_i{J}_i{M}_i\rangle =|\left(n_i{\kappa }_i\right){\mathit{\alpha }}_i{Q}_i{M}_{Q_i}{J}_i{M}_i\rangle ,$$

(17)

where $(Q,M_Q)$ carries the same information as $(w,\nu )$ through the following relations [73,91],

$$\begin{array}{ccc}\hfill Q& =& (2j+1-2\nu )/4=\left(\right|\kappa |-\nu )/2,\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill M_Q& =& -(2j+1-2w)/4=-\left(\right|\kappa |-w)/2.\hfill \end{array}$$

(19)

A nice property of the quasi-spin operator Q$(Q_{+},Q_{-},Q_z)$ is that the ladder operators $Q_{\pm }$ connect subshell wave functions with occupation numbers w and $w\pm 2$ having the same seniority $\nu$. The Wigner–Eckart (WE) theorem can be applied in the space of quasi-spin for all individual subshell states much in the same way as for J-space, allowing an efficient reduction and factorization of matrix elements and CFP matrices [96].

2.5. Second Quantization and Composite Tensor Operators

With complicated CSF structures, second quantization is the most convenient way to handle the mathematics of CSFs and their matrix elements. An extensive $jj$-coupling treatment, adequate for grasp users who have no need for the technical details, can be found in ([73], §6.8, pp. 368–390). Other treatments, mainly focusing on $ls$-coupling, are the original lectures of Brian Judd [94] and Rudzikas’s monograph [91] giving the Vilnius way of doing things. The rangular code from [96] is based on the quasi-spin constructions, replacing the older formulation following Fano [86] used in earlier versions of the atsp2K [28] and grasp [43] packages.

We define the electron creation operator $a_{n\kappa m}^{†}$ as the operator that generates the Dirac orbital when acting on the vacuum state $|0\rangle$ of the electron field

$$\begin{array}{c}\hfill |n\kappa m\rangle =a_{n\kappa m}^{†}|0\rangle .\end{array}$$

(20)

This is destroyed by the annihilation operator $a_{n\kappa m}$ so that

$$\begin{array}{c}\hfill a_{n\kappa m}|n\kappa m\rangle =|0\rangle .\end{array}$$

(21)

Creation and annihilation operators anti-commute so that

$$\begin{array}{c}\hfill \{a_r^{†},a_s\}={\delta }_{rs},\{a_r^{†},a_s^{†}\}=\{a_r,a_s\}=0.\end{array}$$

(22)

Thus, interchanging neighboring pairs of operators introduces a sign change, for example $a_r^{†}{a}_s^{†}=-a_s^{†}{a}_r^{†}$ (implying $a_r^{†}{a}_r^{†}=0$) so that the algebra of creation and annihilation operators is what we need for manipulating antisymmetric states obeying the Pauli exclusion principle: one electron only in each orbital state. Products $a_1^{†}{a}_2^{†}\dots a_n^{†}|0\rangle$ in which the operators are all different therefore represent n-electron determinants from which we can generate more complex n-electron states.

Subshells of Dirac electrons are useful subdivisions of CSFs, as the set of $2j+1$ states $\left\{\right|n\kappa m\rangle ,m=-j,\dots ,+j\}$ is equivalent in the sense that each of them can be expressed in terms of the others under a rotation of axes. Focusing on a single subshell, we can drop n and consider the $\left(\genfrac{}{}{0pt}{}{2j+1}{w}\right)$ ($0\le m\le 2j+1$) equivalent antisymmetric states

$$\begin{array}{c}\hfill |m_1\dots m_w\rangle =a_{m_1}^{†}\cdots a_{m_w}^{†}|0\rangle ,M=\sum _{i=1}^w{m}_i(m_1>\dots >m_w),.\end{array}$$

(23)

In the language of group theory ([73], §6.8.2), this belongs to a reducible representation ${\mathcal{D}}^j\times \cdots \times {\mathcal{D}}^j$ (w factors). This can be decomposed as a direct sum, a Clebsch–Gordan series, of irreducible representations ${\mathcal{D}}^J$, which can be used to characterize the subshell states—see the list in Table 2. These can be further distinguished by the seniority number $\nu$, which is defined for each subshell configuration as the lowest value of w for which a particular J appears. As explained above, the few cases in which $J,\nu$ fails to identify the subshell state uniquely are labeled $Nr$ in this table.

The operators $a_m^{†},{\tilde{a}}_m={(-1)}^{j-m}{a}_{-m}$ can be regarded ([73], §6.8.3) as defining irreducible tensor operators, ${\mathit{a}}^{\kappa †},{\mathit{a}}^{\kappa }$ of rank j. From these, we can build composite tensor operators using angular momentum theory such as

$$\begin{array}{c}\hfill {\left[{\mathit{a}}^{\kappa †}\times {\mathit{a}}^{{\kappa }^{\prime }}\right]}_q^{\left(k\right)}=\sum _{m,m^{\prime }}{a}_m^{\kappa †}{a}_{m^{\prime }}^{{\kappa }^{\prime }}\langle jm,j^{\prime }{m}^{\prime }|j{j}^{\prime }kq\rangle .\end{array}$$

(24)

and more complex operators. Omitting the $\kappa$ rank in the notation of creation and annihilation operators acting on the same subshell (${\kappa }^{\prime }=\kappa$), the three operators

$$\begin{array}{cc}\hfill Q_{+}={\displaystyle \frac{\sqrt{2j+1}}{2}}{\left[{\mathit{a}}^{†}\times {\mathit{a}}^{†}\right]}_0^{\left(0\right)},& Q_{-}=-{\displaystyle \frac{\sqrt{2j+1}}{2}}{\left[\mathit{a}\times \mathit{a}\right]}_0^{\left(0\right)}\hfill \\ \hfill Q_z=-{\displaystyle \frac{\sqrt{2j+1}}{4}}& \left\{{\left[{\mathit{a}}^{†}\times \mathit{a}\right]}_0^{\left(0\right)}+{\left[\mathit{a}\times {\mathit{a}}^{†}\right]}_0^{\left(0\right)}\right\}\hfill \end{array}$$

satisfy the commutation relations of the quasi-spin vector operator $\mathit{Q}$, enabling us to classify subshell states in terms of the quantum numbers $Q,M_Q$ of (18) and (19), which are equivalent to the classification $J,\nu$ in terms of the seniority scheme. The main use of the quasi-spin classification is that $M_Q$ is related to the occupation number w so that relations between states of the subshell configurations ${\kappa }^w$ and ${\kappa }^{w\pm 1}$ can be expressed entirely in terms of a $3j$-symbol in quasi-spin.

2.6. Calculation of Matrix Elements

The above construction of the CSFs makes it possible to derive analytical expressions for matrix elements of one- and two-electron interaction operators. A one-electron irreducible tensor operator has the form

$$F_q^{\left(k\right)}=\sum _{i=1}^N{f}_q^{\left(k\right)}\left(i\right),$$

(25)

where k and q are tensor indices. Examples of one-electron tensor operators are the first part of the Dirac–Coulomb Hamiltonian and the operators describing the interactions with the nuclear dipole and quadrupole moments.

A two-electron scalar operator has the generic form

$$G_0^{\left(0\right)}=\frac{1}{2}\sum _{i\ne j}^N{g}_0^{\left(0\right)}(i,j),$$

(26)

where $g_0^{\left(0\right)}(i,j)$ is a scalar two-body interaction between two electrons

$$g_0^{\left(0\right)}(i,j)=\sum _k{g}_k(r_i,r_j){\mathbf{T}}^{\left(k\right)}\left(i\right)·{\mathbf{T}}^{\left(k\right)}\left(j\right),$$

(27)

where ${\mathbf{T}}^{\left(k\right)}\left(i\right)·{\mathbf{T}}^{\left(k\right)}\left(j\right)$ is the scalar product of two tensors ${\mathbf{T}}^{\left(k\right)}\left(i\right)$ and ${\mathbf{T}}^{\left(k\right)}\left(j\right)$ acting on electron coordinates i and j. Examples are the Coulomb interaction, the Breit interaction and the specific mass shift interaction, arising from the nuclear recoil effect.

The matrix element of an irreducible tensorial operator between two CSFs, $|\gamma J{M}_J\rangle$ and $|{\gamma }^{\prime }{J}^{\prime }{M}_{J^{\prime }}\rangle$, can be factorized thanks to the WE theorem [85],

$$\begin{array}{ccc}\hfill & & \langle \gamma J{M}_J\left|T_q^{\left(k\right)}\right|{\gamma }^{\prime }{J}^{\prime }{M}_{J^{\prime }}\rangle =\hfill \\ & & {(-1)}^{J-M_J}\left(\begin{array}{ccc}J& k& J^{\prime }\\ -M_J& q& M_{J^{\prime }}\end{array}\right)\langle \gamma J\parallel {\mathbf{T}}^{\left(k\right)}{\parallel {\gamma }^{\prime }{J}^{\prime }\rangle }^{Edmonds},\hfill \end{array}$$

(28)

where $\langle \gamma J\parallel {\mathbf{T}}^{\left(k\right)}\parallel {\gamma }^{\prime }{J}^{\prime }\rangle$ is a reduced matrix element independent of the M-quantum numbers and of the q tensorial component. Another version of the WE theorem can be found in Brink and Satchler [97]

$$\begin{array}{ccc}\hfill & & \langle \gamma J{M}_J\left|T_q^{\left(k\right)}\right|{\gamma }^{\prime }{J}^{\prime }{M}_{J^{\prime }}\rangle =\hfill \\ & & {(-1)}^{J-M_J}\sqrt{2J+1}\left(\begin{array}{ccc}J& k& J^{\prime }\\ -M_J& q& M_{J^{\prime }}\end{array}\right)\langle \gamma J\parallel {\mathbf{T}}^{\left(k\right)}{\parallel {\gamma }^{\prime }{J}^{\prime }\rangle }^{Brink-Satchler}.\hfill \end{array}$$

(29)

Equating the l.h.s of (28) and (29) provides a relation between the two definitions of the reduced matrix elements (RMEs)

$$\langle \gamma J\parallel {\mathbf{T}}^{\left(k\right)}\parallel {\gamma }^{\prime }{J}^{\prime }{\rangle }^{Edmonds}=\sqrt{2J+1}\langle \gamma J\parallel {\mathbf{T}}^{\left(k\right)}{\parallel {\gamma }^{\prime }{J}^{\prime }\rangle }^{Brink--Satchler}.$$

(30)

Many formulas in various papers may appear different due to this ambiguity. In the present work, we adopt Edmonds’ formulation (28) of the WE theorem [85] that also fits with Racah [98,99], Judd [100], Cowan [101], or Rudzikas [91], while Brink–Satchler WE theorem (29) is coherent with Rose [102], except for an extra ${(-1)}^{2k}$ phase factor (As pointed out by Judd [103], such a phase factor leaves the relation between RMEs phase-free when k is integral, as is almost always the case. k is indeed integer (even or odd) for all irreducible tensorial operators representing physical quantities, invariant under a $2\pi$ rotation. However, as observed by Judd [94] and Rudzikas [91], phase systems can be crucial for second quantization operators.).

The reduced matrix element of a one-electron operator between the CSFs, in turn, can be expressed as a weighted sum over the active orbitals of one-electron reduced matrix elements

$$\langle \gamma J\parallel {\mathbf{F}}^{\left(k\right)}\parallel {\gamma }^{\prime }{J}^{\prime }\rangle =\sum _{a,b}{\xi }_{ab;k}\langle n_a{\kappa }_a\left(1\right)\parallel {\mathbf{f}}^{\left(k\right)}\left(1\right)\parallel n_b{\kappa }_b\left(1\right)\rangle .$$

(31)

The reduced matrix element $\langle n_a{\kappa }_a\parallel {\mathbf{f}}^{\left(k\right)}\parallel n_b{\kappa }_b\rangle$ depends only on the orbitals and on the nature of the operator. It can be further reduced to a radial integral times a matrix element involving the spherical spinors of the two orbitals. The spin-angular coefficients ${\xi }_{ab;k}$ contain all information about the configuration and the angular couplings. They can be expressed in terms of CFPs and recoupling coefficients necessary to match subshell states in the two CSFs [104,105]. These coefficients are computed by routines of the grasp librang90 library, which are then called by the rangular program for MCDHF approach, the rci code for CI calculations, the rbiotransform and rtransition programs for transition properties, the rhfs and hfszeeman95 codes for hyperfine parameters and magnetic interactions, the ris4 programs for isotope shifts and the rdensity program for radial electron densities and natural orbitals, see accompanying manual §1.3 and §2.2.

For scalar two-electron operators, the application of the WE theorem (28) is trivial

$$\langle \gamma J{M}_J|G_0^{\left(0\right)}|{\gamma }^{\prime }{J}^{\prime }{M}_{J^{\prime }}\rangle =\frac{{\delta }_{J{M}_J,J^{\prime }{M}_J^{\prime }}}{\sqrt{2J+1}}\langle \gamma J\parallel {\mathbf{G}}^{\left(0\right)}\parallel {\gamma }^{\prime }{J}^{\prime }\rangle .$$

(32)

Similarly to (31), the two-body reduced matrix elements can be expressed as [73]

$$\langle \gamma J\parallel {\mathbf{G}}^{\left(0\right)}\parallel {\gamma }^{\prime }{J}^{\prime }\rangle =\sum _k\sum _{abcd}{\xi }_{abcd;k}{X}^k\left(abcd\right),$$

(33)

with

$$X^k\left(abcd\right)={(-1)}^k\langle a\parallel {\mathbf{T}}^{\left(k\right)}\parallel c\rangle \langle b\parallel {\mathbf{T}}^{\left(k\right)}\parallel d\rangle {\mathcal{R}}^k\left(abcd\right).$$

(34)

The summations over $(a,b,c,d)$ in (33) are running over the relativistic subshells $n_a{\kappa }_a$ while the summation over $(i,j)$ in (26) runs over the electron coordinates r${}_i$. The effective interaction strength, $X^k\left(abcd\right)$, is specific to the nature of the interaction and involves only the active orbitals [73,106]. It can be written in terms of a radial double integral ${\mathcal{R}}^k\left(abcd\right)$ and factors involving matrix elements of the spherical spinors of the active orbitals. The spin-angular ${\xi }_{abcd;k}$ coefficients are computed by the routines of the librang90 library, which are then called by the rangular program for MCDHF approach, the rci program for CI calculations, and the ris4 program for isotope shifts; compare to the accompanying manual §2.2.

The Coulomb interaction $1/r_{ij}$

$$\frac{1}{r_{ij}}=\frac{1}{|{\mathbf{r}}_i-{\mathbf{r}}_j|}=\sum _{k=0}^{\infty }\frac{r_{<}^k}{r_{>}^{k+1}}{\mathbf{C}}^{\left(k\right)}({\theta }_i,{\varphi }_i)·{\mathbf{C}}^{\left(k\right)}({\theta }_j,{\varphi }_j),$$

(35)

where ${\mathbf{C}}^{\left(k\right)}(\theta ,\varphi )={[4\pi /(2k+1)]}^{1/2}{Y}_{kq}(\theta ,\varphi )$ is the renormalized spherical harmonic, which fits with the generic form (27), whereas the Gaunt and Breit interactions involve composite tensor operators, $X^{\left(1k\right)K}$, with more convoluted radial parts [73]. The magnetic interactions are therefore more complex but can similarly be reduced to angular factors multiplied by two-electron radial integrals (see Sections 6.4 and 6.5 of [73] for the decomposition and for individual interaction strengths, respectively).

In second quantization, operators such as the angular momentum $\mathbf{J}$ are expressed in the form

$$\begin{array}{c}\hfill \tilde{\mathbf{J}}=\sum _{\kappa ,{\kappa }^{\prime }}\sum _{m,m^{\prime }}{a}_m^{\kappa †}\langle \kappa m\left|\mathbf{J}\right|{\kappa }^{\prime }{m}^{\prime }\rangle a_{m^{\prime }}^{{\kappa }^{\prime }}\end{array}$$

(36)

operating on the ket (23) to the right and on the corresponding bra to the left. Matrix elements using CSFs as in (16) require detaching an active electron orbital, say $|n_a{\kappa }_a{m}_a\rangle$ from a subshell, matching one of the kets, $\langle \kappa m\left|\mathbf{J}\right|{\kappa }^{\prime }{m}^{\prime }\rangle$, in (36). Actions are similar with the bra. This is accomplished ([73], eq. (6.8.30)) by separating $|\kappa m\rangle$ using

$$\begin{array}{c}\hfill \langle {\kappa }^w\nu JM|a_m^{†}|{\kappa }^{w-1}\overline{\nu }\overline{J}\overline{M}\rangle =w^{1/2}\langle {\kappa }^w\nu JM|{\left\{|\kappa m\rangle ·|{\kappa }^{w-1}\overline{\nu }\overline{J}\overline{M}\rangle \right\}}_M^J\end{array}$$

(37)

where

$$\begin{array}{cc}\hfill \langle {\kappa }^w\nu JM|& {\left\{|\kappa m\rangle ·|{\kappa }^{w-1}\overline{\nu }\overline{J}\overline{M}\rangle \right\}}_M^J\hfill \\ \hfill & ={(-1)}^{J-M}\left(\begin{array}{ccc}J& j& \overline{J}\\ -M& m& \overline{M}\end{array}\right){(-1)}^{\overline{J}-j+J}{(2J+1)}^{1/2}(j^w\nu J\left\{\right|j^{w-1}\overline{\nu }\overline{J})\hfill \end{array}$$

where $(j^w\nu J\left\{\right|j^{w-1}\overline{\nu }\overline{J})$ is a CFP. The latter can also be expressed in terms of reduced matrix elements of creation and annihilation operators.

The construction (16) just couples the subshell states in the order chosen by the CSF generator. Decoupling one electron from a ket subshell introduces a $3j$-symbol, a product of the active electron creation operator, and a CFP in the subshell position. Moving the active creation operator to the end of the subshell sequence nearest to the interaction operator where it selects the right matrix element involves interchanging with other creation operators, introducing a factor $-1$ at each interchange. The interchanges introduce an overall phase factor, while the replacement of the original subshell state by the parent left behind changes the CSF coupling scheme, requiring a recoupling coefficient. Further discussion may be found in ([73], §6.9).

A powerful spin-angular algebra based on angular momentum theory, on the concept of the irreducible tensorial sets, second-quantization in a coupled tensorial form, quasi-spin formalism and Wick’s theorem, reduced (in quasi-spin space) coefficients of fractional parentage and on a generalized graphical method, is used in grasp to evaluate the ${\xi }_{ab;k}$ and ${\xi }_{abcd;k}$ coefficients appearing in (31) and (33), as described in detail in [89,96,106,107]. The corresponding spin-angular library, implemented in grasp2018, is fully documented by Gaigalas [93] in the present Special Issue. Interesting illustrations on how this spin-angular algebra is applied for some physical quantities can be found in [108,109] in $jj$- and $LSJ$-coupling, respectively.

2.7. Multiconfiguration Dirac–Hartree–Fock

In the MCDHF method, the wave function of an atomic state $\mathsf{\Gamma }J{M}_J$ is approximated by the atomic state function (6). For spectrum calculations, the orbital optimization often targets simultaneously several atomic states that may belong to different J symmetries, ${\mathsf{\Gamma }}^i{J}^i,i=1,\dots ,N_{ASF}$ (suppressing the $M_J$ and parity $\pi$ quantum numbers for brevity). In this scheme, the different ASFs belonging to the same J ($J=J^i=J^j$) are chosen to be orthonormal, so that

$${\left({\mathbf{c}}^{{\mathsf{\Gamma }}^i{J}^i}\right)}^{†}{\mathbf{c}}^{{\mathsf{\Gamma }}^j{J}^j}={\delta }_{i,j},$$

(38)

where ${\mathbf{c}}^{\mathsf{\Gamma }J}$ is the column vector collecting the mixing coefficients $\{c_{\mathit{\alpha }}^{\mathsf{\Gamma }J},\mathit{\alpha }=1,\dots ,N_{CSF}\}$ for a given ASF. The energy of the atomic state $\mathsf{\Gamma }J$ is

$$\begin{array}{ccc}\hfill E_{\mathsf{\Gamma }J}=\langle \mathsf{\Gamma }JM|{\mathcal{H}}_{DC}|\mathsf{\Gamma }JM\rangle & =& \frac{1}{\sqrt{2J+1}}\langle \mathsf{\Gamma }J\left|\right|{\mathcal{H}}_{DC}\left|\right|\mathsf{\Gamma }J\rangle \hfill \\ & =& \frac{1}{\sqrt{2J+1}}{\left({\mathbf{c}}^{\mathsf{\Gamma }J}\right)}^{†}\mathbf{H}{\mathbf{c}}^{\mathsf{\Gamma }J},\hfill \end{array}$$

(39)

where the Hamiltonian reduced matrix $\mathbf{H}$ has the elements

$$H_{\mathit{\alpha }\beta }=\langle {\gamma }_{\mathit{\alpha }}J\left|\right|{\mathcal{H}}_{DC}\left|\right|{\gamma }_{\beta }J\rangle .$$

(40)

The reduced matrix elements (40) can, as discussed above, be expressed in terms of angular coefficients and radial integrals [30,73]

$$H_{\mathit{\alpha }\beta }=\langle {\gamma }_{\mathit{\alpha }}J\left|\right|{\mathcal{H}}_{DC}\left|\right|{\gamma }_{\beta }J\rangle =\sum _{ab}{t}_{ab}^{\mathit{\alpha }\beta }I(a,b)+\sum _{abcd;k}{v}_{abcd;k}^{\mathit{\alpha }\beta }{R}^k(ab;cd),$$

(41)

where now the factors involving matrix elements of the spherical spinors are absorbed in the angular coefficients. The one-electron radial integrals are given by

$$\begin{array}{ccc}\hfill I(a,b)={\delta }_{{\kappa }_a{\kappa }_b}& & {\int }_0^{\infty }\left\{P_{n_a{\kappa }_a}\left(r\right)V_{nuc}\left(r\right)P_{n_b{\kappa }_b}\left(r\right)-c{P}_{n_a{\kappa }_a}\left(r\right)\left(\frac{d}{dr}-\frac{\kappa }{r}\right)Q_{n_b{\kappa }_b}\left(r\right)\right.\hfill \\ & & \left.+c{Q}_{n_a{\kappa }_a}\left(r\right)\left(\frac{d}{dr}+\frac{\kappa }{r}\right)P_{n_b{\kappa }_b}\left(r\right)+Q_{n_a{\kappa }_a}\left(r\right)\left(V_{nuc}\left(r\right)-2c^2\right)Q_{n_b{\kappa }_b}\left(r\right)\right\}dr.\hfill \end{array}$$

(42)

The two-electron radial double integrals, the so-called Slater integrals, are given by

$$\begin{array}{ccc}\hfill R^k(ab;cd)=& & {\int }_0^{\infty }\left[P_{n_a{\kappa }_a}\left(r\right)P_{n_c{\kappa }_c}\left(r\right)\right.\hfill \\ & & \left.+Q_{n_a{\kappa }_a}\left(r\right)Q_{n_c{\kappa }_c}\left(r\right)\right]\frac{1}{r}{Y}^k(bd;r)dr.\hfill \end{array}$$

(43)

with $Y^k$ defined by

$$\begin{array}{ccc}\hfill Y^k(ab;r)=& & r{\int }_0^{\infty }\frac{r_{<}^k}{r_{>}^{k+1}}\left[P_{n_a{\kappa }_a}\left(s\right)P_{n_b{\kappa }_b}\left(s\right)\right.\hfill \\ & & \left.+Q_{n_a{\kappa }_a}\left(s\right)Q_{n_b{\kappa }_b}\left(s\right)\right]ds.\hfill \end{array}$$

(44)

In the $Y^k$ integral, $r_{<}$ and $r_{>}$ are the smaller and larger of r and s, respectively.

Multiconfiguration methods are energy driven. Introducing the radial orthonormality condition

$${\mathcal{C}}_{ab}\equiv \int [P_{n_a{\kappa }_a}\left(r\right)P_{n_b{\kappa }_b}\left(r\right)+Q_{n_a{\kappa }_a}\left(r\right)Q_{n_b{\kappa }_b}\left(r\right)]dr-{\delta }_{n_a,n_b}=0,$$

(45)

and applying the variational principle on the statistically weighted energy functional of the targeted states ${\gamma }^i{J}^i,i=1,\dots ,N_{ASF}$

$$\mathcal{F}(\left\{c\right\},\left\{P\right\},\left\{Q\right\})=\frac{{\sum }_{i=1}^{N_{ASF}}(2J^i+1)E_{{\mathsf{\Gamma }}^i{J}^i}}{{\sum }_{i=1}^{N_{ASF}}(2J^i+1)}+\sum _{ab}{\delta }_{{\kappa }_a{\kappa }_b}{\lambda }_{ab}{\mathcal{C}}_{ab}.$$

(46)

where Lagrange multipliers ${\lambda }_{ab}$ are introduced to ensure the orthonormality of the orbitals, yields the equations for the radial functions $P_{n_a{\kappa }_a}\left(r\right),Q_{n_a{\kappa }_a}\left(r\right)$, see [30,73]

$$w_a\left[\begin{array}{ccc}V(a;r)\hfill & \hfill & -c\left[\frac{d}{dr}-\frac{{\kappa }_a}{r}\right]\hfill \\ c\left[\frac{d}{dr}+\frac{{\kappa }_a}{r}\right]\hfill & \hfill & V(a;r)-2c^2\hfill \end{array}\right]\left[\begin{array}{c}{P}_{n_a{\kappa }_a}\left(r\right)\hfill \\ Q_{n_a{\kappa }_a}\left(r\right)\hfill \end{array}\right]=\sum _b{ϵ}_{ab}{\delta }_{{\kappa }_a,{\kappa }_b}\left[\begin{array}{c}{P}_{n_a{\kappa }_a}\left(r\right)\hfill \\ Q_{n_a{\kappa }_a}\left(r\right)\hfill \end{array}\right].$$

(47)

The potential $V(a;r)$ consists of three terms

$$V(a;r)=V_{nuc}\left(r\right)+Y(a;r)+\overline{X}(a;r),$$

(48)

where the variations of the $R^k(ab;ab)$ integrals weighted with the angular coefficients and the state averaged expansion coefficients contribute to $Y(a;r)$. Variations of the other $R^k(ab;cd)$ integrals and the off-diagonal $I(a,b)$ integrals, again weighted with the angular coefficients and the state-averaged expansion coefficients, contribute to the exchange operator $\overline{X}(a;r)$ [110]. $w_a$ is the generalized occupation number of orbital a and ${ϵ}_{ab}$ are energy parameters related to the Lagrange multipliers.

In B-spline methods [111], once off-diagonal Lagrange multipliers have been dealt with, the radial functions $P_a\left(r\right)=P_{n_a{\kappa }_a}\left(r\right)$, $Q_a\left(r\right)=Q_{n_a{\kappa }_a}\left(r\right)$, are obtained as the solution of an eigenvalue problem. However, the differential equation method used by grasp, where the radial functions are defined as a vector of values on a grid, requires that the Dirac equation be rewritten as a pair of first-order differential equations [33,73,112]. In the first equation, we introduce $X_a^P\left(r\right)=\overline{X}(a;r)P_a\left(r\right)={\sum }_b{c}_{abk}(1/r)Y^k(ab;r)P_b\left(r\right)$ as a non-homogeneous term of the differential equation. Similarly, for the second equation, we let $X_a^Q\left(r\right)=\overline{X}(a;r)Q_a\left(r\right)={\sum }_b{c}_{abk}(1/r)Y^k(ab;r)Q_b\left(r\right)$. Furthermore, it is customary to have the coefficient of the highest derivative be unity, so that

$$\begin{array}{c}\hfill (V_{nuc}+Y(a;r))P_a\left(r\right)-c(\frac{d}{dr}-\frac{{\kappa }_a}{r})Q_a\left(r\right)=\frac{1}{w_a}\sum _b{ϵ}_{ab}{\delta }_{{\kappa }_a{\kappa }_b}{P}_b\left(r\right)-X^P(a;r)\\ \hfill c(\frac{d}{dr}+\frac{{\kappa }_a}{r})P_a\left(r\right)+(V_{nuc}\left(r\right)+Y(a;r)-2c^2)Q_a\left(r\right)=\frac{1}{w_a}\sum _b{ϵ}_{ab}{\delta }_{{\kappa }_a{\kappa }_b}{Q}_b\left(r\right)-X^Q(a;r)\end{array}$$

(49)

As a last step, the first equation should also be multiplied by $-1$, which considerably changes the symmetry. In the self-consistent iterative method, $Y(a;r),{ϵ}_{ab}(b\ne a)$, $X^P(a;r)$, $X^Q(a;r)$ are computed from current estimates and the orbital energy ${ϵ}_{aa}$, and updated radial functions are solutions of the differential equation.

For bound states, it is required that the solutions are square integrable. A necessary, but not sufficient, condition is that the radial amplitudes vanish as $r\to 0$ and $r\to \infty$. For a point charge nuclear model, a rigorous analysis of the radial behavior at $r=0$ leads to power series of the form

$$P_{n_a{\kappa }_a}\left(r\right)=r^{\nu }(p_0+p_{1}r+\dots ),Q_{n_a{\kappa }_a}\left(r\right)=r^{\nu }(q_0+q_{1}r+\dots ),$$

(50)

where $2\nu >-1$ and $p_i,q_i$ are determined by coupled linear algebraic equations. For a more realistic model of the nuclear charge distribution, other conditions apply [73]. The Equation (47), together with the accompanying boundary conditions, are solved iteratively on the radial grid (8) by a finite difference self-consistent field (SCF) procedure. Briefly, given initial estimates of the radial functions, the Hamiltonian matrix with elements (40) is constructed and diagonalized to give expansion coefficients of the CSFs for each of the ASFs. Improved estimates of the radial functions are then obtained by solving Equation (47). The two last steps are repeated until the energy of the states and/or the radial functions do not change anymore.

It is desirable to optimize all orbitals simultaneously, which is sometimes referred to as the “full variational” (FV) approach. However, due to numerical convergence issues, the MCDHF method often employs a layer-by-layer (LBL) strategy [113,114] (see Section 4.4, and also in §3.5 of the manual [72]), in which only the newly introduced orbitals for the layer considered are optimized while the remaining ones are kept frozen. In this context, a layer is a set of new orbitals to be optimized consisting of one orbital per angular momentum symmetry. The LBL approach is attractive as the computation time for each new layer is much shorter than the corresponding computation time of the FV approach. The price to pay for the LBL strategy is a larger active set of correlation orbitals to compensate for the lost degrees of freedom.

Although grasp focuses on atomic bound states, the MCDHF equations can generate continuum orbitals provided different boundary conditions at $r_{max}$. grasp ASFs can then be used as target states in collision calculations in codes such as R-matrix [73,115]. In the grasp suite of programs, the MCDHF equations are solved by the rmcdhf program, which reads the nuclear parameters, the CSFs list, as generated by the rcsfgenerate program and the necessary angular data, produced by rangular, from disk files. The initial estimates of the radial orbitals for the SCF procedure are generated by the rwfnestimate program and can be taken as screened hydrogenic functions, functions from a Thomas–Fermi calculation or converted non-relativistic orbitals; see the accompanying manual §3.3. For an overview of the program and file flows for the MCDHF calculations, see Figures 1 and 2 of the manual.

2.8. Configuration Interaction

The MCDHF calculations are used to generate an orbital basis. Given this basis, the final wave functions for the targeted states are obtained in relativistic configuration interaction (CI) calculations based on the Dirac–Coulomb–Breit and QED Hamiltonian (13). For that step, the CSF expansions (6) may be extended in comparison with those used in the previous MCDHF calculations to capture higher-order correlation contributions or to test the adequacy of correlation models. In CI calculations, the matrix elements of the CSFs are computed using the methods described in Section 2.6. The expansion coefficients of the CSFs in the targeted states are then obtained by diagonalizing the Hamiltonian matrix. The number of CSFs depends on the shell structure of the atomic system in question as well as on the model for electron correlation. Even for moderate calculations, the number of zeroes in the interaction exceeds the non-zero values, and sparse matrix methods are used. For accurate calculations, a large number of CSFs are required, leading to very large matrices. To handle these large matrices, the CSFs can a priori be divided into two groups. Let M be the size of the total CSF space. The first group, P, with m elements ($m\ll M$) contains CSFs that account for the major parts of the wave functions. The second group, Q, with $M-m$ elements contains CSFs that represent minor corrections. Allowing interaction between CSFs in group P, interaction between CSFs in group P and Q and diagonal interactions between CSFs in Q gives a matrix

$$\left(\begin{array}{cc}{H}^{\left(PP\right)}& H^{\left(PQ\right)}\\ H^{\left(QP\right)}& H^{\left(QQ\right)}\end{array}\right),$$

(51)

where $H_{ij}^{\left(QQ\right)}={\delta }_{ij}{E}_i^Q$. The restriction of $H^{\left(QQ\right)}$ to diagonal elements results in a huge reduction in the total number of matrix elements and corresponding computational time. The assumptions of the approximation and the connections to the method of deflation in numerical analysis are discussed in [30]. The structure of (51) is reminiscent of the second-order Brillouin–Wigner perturbation theory [116,117,118]. This form of the CI matrix has been available in the non-relativistic and relativistic multiconfiguration codes for a long time [41,119], where the P and Q groups were named “zeroth-order” and “first-order” sets, respectively. Our MCDHF/CI methods are therefore referred to as ‘Zero-First’ methods [120], i.e., ZF-MCDHF or ZF-CI.

grasp uses sparse matrix methods for storing the interaction matrix in which only non-zero off-diagonal elements are stored by column for the upper matrix, taking advantage of matrix symmetry. The order of the CSFs is not important except when initial estimates of the eigenvectors are required, in which case up to 4000 of the first CSFs are considered, depending on a program parameter. Special iterative eigensolvers are used that determine only selected eigenvalues in the lower portion of the spectrum [121] and matrix–vector multiplication with sparse matrices. Much of the efficiency of variational methods relies on the sparsity of the interaction matrix. In the grasp suite of programs, the CI calculations are performed by the rci program. Nuclear parameters, the CSFs list, as generated by rcsfgenerate and radial orbitals, as produced by rmcdhf, are read from disk files. Angular data needed to compute the Hamiltonian matrix elements are computed on the fly by calls to the routines of the librang90 library. The program and file flows associated with a CI calculation are displayed in Figures 1 and 2 in the accompanying manual. In grasp, the program rcsfzerofirst is used to partition the CSF expansion in zero- and first-order sets. The ZF method is not the default mode but is only used to handle very large expansions and matrices, as discussed in detail in section 14 of the manual [72] that is entirely devoted to strategies for ZF-MCDHF and ZF-CI.

2.9. Transformation to Different Coupling Schemes

In fully relativistic calculations, quantum labels for the targeted states are obtained in $jj$-coupling. Most often, this wave function representation is far from being pure, i.e., there is no dominant CSF whose quantum numbers can be used to label a state in a proper way. Using the methods developed by Gaigalas and co-workers [70,122], the wave function representation in $jj$-coupling is transformed to an approximate representation in $LSJ$-coupling. This representation is normally purer and better suited for labeling. One should be aware of the fact that even in $LSJ$-coupling, the labeling is not straightforward, and several components in the $LSJ$-coupling representation must be used recursively to find unique labels [123]. In grasp2018, the transformation from $jj$- to $LSJ$-coupling is completed by the jj2lsj program that also determines unique labels. Transformations to other coupling schemes, e.g., $JK$ or $LK$, are completed by the coupling program; see the accompanying manual §6.2 and §6.3. Needless to say, the programs to transform the wave functions and assigning unique labels are important parts of the grasp package.

3. Grasp Theory—Atomic Properties

Once the ASFs have been determined from MCDHF or CI calculations, measurable properties such as hyperfine structure splittings, isotope shifts and rates for transitions between two states can be computed as expectation values of one- and two-electron operators. For hyperfine structure and isotope shift, the combination of calculated electronic quantities and high precision measurements allows nuclear parameters to be extracted [124]. Atomic transition rates, on the other hand, are crucial in astrophysics and plasma physics for diagnostic purposes, e.g., for determining element abundances, temperatures, and electron densities [46].

3.1. Hyperfine Structures

The hyperfine structure of a fine structure level is caused by the interaction between the electrons and the electromagnetic multipole moments of the nucleus. The Hamiltonian for the interaction may be written as a multipole expansion

$${\mathcal{H}}_{hfs}=\sum _{k\ge 1}{\mathbf{T}}^{\left(k\right)}·{\mathbf{M}}^{\left(k\right)},$$

(52)

where ${\mathbf{T}}^{\left(k\right)}$ and ${\mathbf{M}}^{\left(k\right)}$ are spherical tensor operators of rank k in the electronic and nuclear spaces, respectively [125]. The $k=1$ term represents the magnetic dipole interaction and the $k=2$ term the electric quadrupole interaction. For an N-electron atom, the electronic tensor operators are [66,125]

$${\mathbf{T}}^{\left(1\right)}=\sum _{j=1}^N{\mathbf{t}}^{\left(1\right)}\left(j\right)=\sum _{j=1}^N-i\sqrt{2}\mathit{\alpha }{r}_j^{-2}{\left({\mathit{\alpha }}_j{\mathbf{C}}^{\left(1\right)}({\theta }_j,{\phi }_j)\right)}^{\left(1\right)}$$

(53)

$${\mathbf{T}}^{\left(2\right)}=\sum _{j=1}^N{\mathbf{t}}^{\left(2\right)}\left(j\right)=\sum _{j=1}^N-r_j^{-3}{\mathbf{C}}^{\left(2\right)}({\theta }_j,{\phi }_j).$$

(54)

The matrix elements of the nuclear tensor operators are related to the conventional nuclear magnetic dipole moment ${\mu }_I$ and electric quadrupole moment Q by

$${\mu }_I=\langle II|M_0^{\left(1\right)}|II\rangle$$

(55)

$$Q=2\langle II\left|M_0^{\left(2\right)}\right|II\rangle$$

(56)

where $|I{M}_I\rangle$ is the nuclear wave function. Values of nuclear magnetic dipole and electric quadrupole moments can be found in the recent compilation by Stone [126,127].

The hyperfine interaction couples the nuclear $\mathbf{I}$ and electronic $\mathbf{J}$ angular momenta to a total momentum $\mathbf{F}=\mathbf{I}+\mathbf{J}$ and zero-order wave functions of the coupled states can be written

$$|\mathsf{\Gamma }IJF{M}_F\rangle =\sum _{M_I,M_J}\langle IJ{M}_I{M}_J|IJF{M}_F\rangle |I{M}_I\rangle |\mathsf{\Gamma }J{M}_J\rangle .$$

(57)

If the hyperfine interaction is weak so that the interaction energy is small compared to the fine structure separation, ${\mathcal{H}}_{hfs}$ can be treated in first-order perturbation theory. The F-dependent hyperfine energy to be added to the energy of the $\mathsf{\Gamma }J$ level is then given by

$$E_{\mathsf{\Gamma }IJF}^{hfs}=\langle \mathsf{\Gamma }IJF{M}_F|\sum _{k\ge 1}{\mathbf{T}}^{\left(k\right)}·{\mathbf{M}}^{\left(k\right)}|\mathsf{\Gamma }IJF{M}_F\rangle ,$$

(58)

where

$${\mathbf{T}}^{\left(k\right)}·{\mathbf{M}}^{\left(k\right)}={(-1)}^k\sqrt{2k+1}{\left[{\mathbf{T}}^{\left(k\right)}\times {\mathbf{M}}^{\left(k\right)}\right]}_0^{\left(0\right)}.$$

(59)

Factoring out the dependence on the F quantum number, the hyperfine energies can be expressed in terms of the hyperfine interaction constants, A and B,

$$E_{\mathsf{\Gamma }IJF}^{hfs}=\frac{1}{2}{A}_{\mathsf{\Gamma }J}C+B_{\mathsf{\Gamma }J}\frac{\frac{3}{4}C(C+1)-I(I+1)J(J+1)}{2I(2I-1)J(2J-1)},$$

(60)

where

$$A_{\mathsf{\Gamma }J}=\frac{{\mu }_I}{I}\frac{1}{\sqrt{J(J+1)(2J+1)}}\langle \mathsf{\Gamma }J\parallel {\mathbf{T}}^{\left(1\right)}\parallel \mathsf{\Gamma }J\rangle ,$$

(61)

$$B_{\mathsf{\Gamma }J}=2Q\sqrt{\frac{J(2J-1)}{(J+1)(2J+1)(2J+3)}}\langle \mathsf{\Gamma }J\parallel {\mathbf{T}}^{\left(2\right)}\parallel \mathsf{\Gamma }J\rangle .$$

(62)

and $C=F(F+1)-J(J+1)-I(I+1)$. Inserting the expression for the ASFs in terms of CSFs (6) and using the results of Section 2.6, the reduced matrix elements in the hyperfine interaction constants can be computed as sums over radial hyperfine structure integrals weighted with the product of the expansion coefficients of the CSFs and the angular coefficients, as expressed by Equation (31) for one-body operators. In grasp, the hyperfine structure constants are computed using the rhfs program. The nuclear parameters, the CSF expansion, the radial orbitals, as obtained by rmcdhf, and the expansion coefficients of the CSFs, as obtained by the rmcdhf or rci programs, are read from files. Specific examples of computations of hyperfine constants are given in §6.1 of the manual. For more recent developments, using several independently optimized radial orbital sets in the computation of the hyperfine constants, see the article by Yan Ting et al. [128] in the present Special Issue.

3.2. External Magnetic Fields

Neglecting diamagnetic contributions and choosing the direction of the external magnetic field $\mathbf{B}$ in the z-direction, the interaction between the atomic electrons and the field can be written

$${\mathcal{H}}_m=({\mathbf{N}}^{\left(1\right)}+\mathsf{\Delta }{\mathbf{N}}^{\left(1\right)})·{\mathbf{B}}^{\left(1\right)}\equiv (N_0^{\left(1\right)}+\mathsf{\Delta }{N}_0^{\left(1\right)})B,$$

(63)

where $B=B_0^{\left(1\right)}=B_z$ is the magnetic field strength. The last term is the so-called Schwinger QED correction. For the Zeeman effect on hyperfine levels, the Hamiltonian for the Zeeman interaction should in principle also include the interaction with the nucleus, although it is weak and can be neglected.

For an N-electron atom, the electronic tensor operators are [84,129]

$${\mathbf{N}}^{\left(1\right)}=\sum _{j=1}^N{\mathbf{n}}^{\left(1\right)}\left(j\right)=\sum _{j=1}^N-i\frac{\sqrt{2}}{2\mathit{\alpha }}{r}_j{\left({\mathit{\alpha }}_j{\mathbf{C}}^{\left(1\right)}({\theta }_j,{\phi }_j)\right)}^{\left(1\right)},$$

(64)

$$\mathsf{\Delta }{\mathbf{N}}^{\left(1\right)}=\sum _{j=1}^N\mathsf{\Delta }{\mathbf{n}}^{\left(1\right)}\left(j\right)=\sum _{j=1}^N\frac{g_s-2}{2}{\beta }_j{\Sigma }_j,$$

(65)

where ${\Sigma }_j$ is the relativistic spin-matrix and $g_s=2.00232$ is the g factor of the electron spin corrected for QED effects [129]. If the magnetic field is weak, so that the interaction energy is small compared to the fine-structure separation, the interaction can be treated in first-order perturbation theory. A fine-structure level $\mathsf{\Gamma }J$ is then split according to   

$$\begin{array}{c}{\displaystyle \langle \mathsf{\Gamma }J{M}_J|N_0^{\left(1\right)}+\mathsf{\Delta }{N}_0^{\left(1\right)}|\mathsf{\Gamma }J{M}_J\rangle B}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}={(-1)}^{J-M_J}\left(\begin{array}{ccc}J& 1& J\\ -M_J& 0& M_J\end{array}\right)\langle \mathsf{\Gamma }J\parallel {\mathbf{N}}^{\left(1\right)}+\mathsf{\Delta }{\mathbf{N}}^{\left(1\right)}\parallel \mathsf{\Gamma }J\rangle B\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}=\frac{M_J}{\sqrt{J(J+1)(2J+1)}}\langle \mathsf{\Gamma }J\parallel {\mathbf{N}}^{\left(1\right)}+\mathsf{\Delta }{\mathbf{N}}^{\left(1\right)}\parallel \mathsf{\Gamma }J\rangle B.\hfill \end{array}$$

(66)

Usually, the dependence on the $M_J$ quantum number is factored out, and the energy splittings are expressed as

$$E_{\mathsf{\Gamma }J{M}_J}=g_J{M}_J\frac{B}{2},$$

(67)

where $g_J$ is the Landé factor

$$g_J=2\frac{\langle \mathsf{\Gamma }J\parallel {\mathbf{N}}^{\left(1\right)}+\mathsf{\Delta }{\mathbf{N}}^{\left(1\right)}\parallel \mathsf{\Gamma }J\rangle }{\sqrt{J(J+1)(2J+1)}}.$$

(68)

Analogously, a hyperfine level $\mathsf{\Gamma }IJF$ for which off-diagonal effects are small is then split according to

$$\begin{array}{c}{\displaystyle \langle \mathsf{\Gamma }IJF{M}_F|N_0^{\left(1\right)}+\mathsf{\Delta }{N}_0^{\left(1\right)}|\mathsf{\Gamma }IJF{M}_F\rangle B}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}=M_F\frac{F(F+1)+J(J+1)-I(I+1)}{2F(F+1)}\frac{\langle \mathsf{\Gamma }J\parallel {\mathbf{N}}^{\left(1\right)}+\mathsf{\Delta }{\mathbf{N}}^{\left(1\right)}\parallel \mathsf{\Gamma }J\rangle }{\sqrt{J(J+1)(2J+1)}}B\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}=M_F\frac{F(F+1)+J(J+1)-I(I+1)}{2F(F+1)}{g}_J\frac{B}{2}.\hfill \end{array}$$

(69)

The energy splittings of hyperfine levels are given by

$$E_{\mathsf{\Gamma }IJF{M}_F}=g_F{M}_F\frac{B}{2},$$

(70)

where the Landé factor is defined by

$$g_F=\frac{F(F+1)+J(J+1)-I(I+1)}{2F(F+1)}{g}_J.$$

(71)

In grasp, the Landé factors for fine- and hyperfine structure levels are computed using the hfszeeman95 program. The nuclear parameters, the CSF expansion, the radial orbitals, as obtained by rmcdhf, and the expansion coefficients of the CSFs, as obtained by the rmcdhf or rci programs, are read from files. Specific examples are given in §6.9 of the accompanying manual.

3.3. Isotope Shift

3.3.1. Mass Shift

The Dirac–Coulomb and Dirac–Coulomb–Breit Hamiltonians are valid under the assumptions that the nucleus is infinitely heavy. The effect of the recoil motion of the nucleus is given by the normal and specific mass shift operators [130,131]

$$\begin{array}{ccc}\hfill {\mathcal{H}}_{\mathrm{NMS}}^A& =& \frac{1}{2M}\sum _{j=1}^N\left({\mathbf{p}}_j^2-\frac{\mathit{\alpha }Z}{r_j}{\mathit{\alpha }}_j·{\mathbf{p}}_j-\frac{\mathit{\alpha }Z}{r_j}\left({\mathit{\alpha }}_j·{\mathbf{C}}^{\left(1\right)}({\theta }_j,{\phi }_j)\right){\mathbf{C}}^{\left(1\right)}({\theta }_j,{\phi }_j)·{\mathbf{p}}_j\right),\hfill \end{array}$$

(72)

$$\begin{array}{ccc}\hfill {\mathcal{H}}_{\mathrm{SMS}}^A& =& \frac{1}{2M}\sum _{j\ne k}^N\left({\mathbf{p}}_j·{\mathbf{p}}_k-\frac{\mathit{\alpha }Z}{r_j}{\mathit{\alpha }}_j·{\mathbf{p}}_k-\frac{\mathit{\alpha }Z}{r_j}\left({\mathit{\alpha }}_j·{\mathbf{C}}^{\left(1\right)}({\theta }_j,{\phi }_j)\right){\mathbf{C}}^{\left(1\right)}({\theta }_j,{\phi }_j)·{\mathbf{p}}_k\right),\hfill \end{array}$$

(73)

where M is the nuclear mass of the isotope with mass number A. Treating the operators in first-order perturbation theory, with ASFs of the Hamiltonian with infinite nuclear mass as zero-order functions, the isotope mass shift of a level $\mathsf{\Gamma }J$ is given by

$$\begin{array}{ccc}\hfill E_{\mathsf{\Gamma }J}^A-E_{\mathsf{\Gamma }J}^{A^{\prime }}=& & K_{\mathsf{\Gamma }J,\mathrm{MS}}\left(\frac{M^{\prime }-M}{M{M}^{\prime }}\right)=\hfill \\ & & (K_{\mathsf{\Gamma }J,\mathrm{NMS}}+K_{\mathsf{\Gamma }J,\mathrm{SMS}})\left(\frac{M^{\prime }-M}{M{M}^{\prime }}\right),\hfill \end{array}$$

(74)

where the level mass shift parameters $K_{\mathsf{\Gamma }J,\mathrm{NMS}}$ and $K_{\mathsf{\Gamma }J,\mathrm{SMS}}$ are defined by

$$\frac{K_{\mathsf{\Gamma }J,\mathrm{NMS}}}{M}=\langle \mathsf{\Gamma }J{M}_J|{\mathcal{H}}_{\mathrm{NMS}}^A|\mathsf{\Gamma }J{M}_J\rangle =\frac{1}{\sqrt{2J+1}}\langle \mathsf{\Gamma }J\parallel {\mathcal{H}}_{\mathrm{NMS}}^A\parallel \mathsf{\Gamma }J\rangle ,$$

(75)

$$\frac{K_{\mathsf{\Gamma }J,\mathrm{SMS}}}{M}=\langle \mathsf{\Gamma }J{M}_J|{\mathcal{H}}_{\mathrm{SMS}}^A|\mathsf{\Gamma }J{M}_J\rangle =\frac{1}{\sqrt{2J+1}}\langle \mathsf{\Gamma }J\parallel {\mathcal{H}}_{\mathrm{SMS}}^A\parallel \mathsf{\Gamma }J\rangle .$$

(76)

Rewriting the operators in tensorial form, the expectation values can be evaluated according to Equations (31) and (33) in terms of radial integrals (details are given in [132]). In the grasp suite of codes, the level mass shift electronic factors $K_{\mathsf{\Gamma }J,\mathrm{NMS}}$ and $K_{\mathsf{\Gamma }J,\mathrm{SMS}}$ are computed by the ris4 program (see §6.1 in the manual for a specific example). The three contributions associated with the three terms of (72) and (73) are reported separately.

3.3.2. Field Shift

In grasp, the nuclear potential results from a two-parameter Fermi nuclear charge distribution and the energy shift for a level from one isotope to another can in principle be obtained by performing two separate calculations and subtracting the energies. This is inconvenient and computationally expensive. A better way is to consider the level field shift (FS) as a first-order perturbation

$$\delta E_{\mathsf{\Gamma }J,FS}^{A,A^{\prime }}=-{\int }_{{\mathbf{R}}^3}\left[V^A\left(\mathbf{r}\right)-V^{A^{\prime }}\left(\mathbf{r}\right)\right]{\rho }_{\mathsf{\Gamma }J}^e\left(\mathbf{r}\right)d^3\mathbf{r},$$

(77)

where $V^A\left(\mathbf{r}\right)$ and $V^{A^{\prime }}\left(\mathbf{r}\right)$ are the potentials arising from the nuclear charge distributions of the two isotopes and

$${\rho }_{\mathsf{\Gamma }J}^e\left(\mathbf{r}\right)=\langle \mathsf{\Gamma }J{M}_J|\sum _{i=1}^N\delta (\mathbf{r}-{\mathbf{r}}_i)|\mathsf{\Gamma }J{M}_J\rangle$$

(78)

is the electron density of level $\mathsf{\Gamma }J$ of the reference isotope. Assuming that the electron density at the origin can be very well approximated with a spherically symmetric even polynomial function according to

$${\rho }_{\mathsf{\Gamma }J}^e\left(\mathbf{r}\right)=\frac{{\rho }_{\mathsf{\Gamma }J}^e\left(r\right)}{4\pi }\approx b_{\mathsf{\Gamma }J}\left(r\right)=b_{\mathsf{\Gamma }J,0}+b_{\mathsf{\Gamma }J,2}{r}^2+b_{\mathsf{\Gamma }J,4}{r}^4+b_{\mathsf{\Gamma }J,6}{r}^6,$$

(79)

and using the Laplacian operator in spherical coordinates to obtain ${\widehat{\nabla }}^2{r}^{n+2}=(n+2)(n+3)r^n$ and the Poisson’s equation ${\widehat{\nabla }}^2{V}^A\left(\mathbf{r}\right)=-4\pi {\rho }^A\left(\mathbf{r}\right)$, we obtain

$$\delta E_{\mathsf{\Gamma }J,FS}^{A,A^{\prime }}=\sum _{0\le n\le 6,n\mathrm{even}}{F}_{\mathsf{\Gamma }J,n}\delta {\langle r^{n+2}\rangle }^{A,A^{\prime }},$$

(80)

with

$$F_{\mathsf{\Gamma }J,n}=\frac{4\pi Z{b}_{\mathsf{\Gamma }J,n}}{(n+2)(n+3)},$$

(81)

where we used the expression for differences in the nuclear radial moments between isotopes A and $A^{\prime }$,

$$\delta {\langle r^{n+2}\rangle }^{A,A^{\prime }}={\langle r^{n+2}\rangle }^A-{\langle r^{n+2}\rangle }^{A^{\prime }}=\frac{1}{Z}{\int }_{{\mathbf{R}}^3}{r}^{n+2}\left[{\rho }^A\left(\mathbf{r}\right)-{\rho }^{A^{\prime }}\left(\mathbf{r}\right)\right]d^3\mathbf{r}.$$

(82)

In the expression above, Z is the atomic number and appears due to normalization. Accounting for the level isotope shifts of an upper u and a lower l state in a transition gives the experimentally measurable transition isotope shift

$$\delta {\nu }_{k,IS}^{A,A^{\prime }}=(\mathsf{\Delta }{K}_{NMS}+\mathsf{\Delta }{K}_{SMS})\left(\frac{M^{\prime }-M}{M{M}^{\prime }}\right)+\sum _{0\le n\le 6,n\mathrm{even}}\frac{\mathsf{\Delta }{F}_{k,n}}{h}\delta {\langle r^{n+2}\rangle }^{A,A^{\prime }},$$

(83)

where $\mathsf{\Delta }{K}_{NMS}=K_{u,NMS}-K_{l,NMS}$ and $\mathsf{\Delta }{K}_{SMS}=K_{u,SMS}-K_{l,SMS}$ are the differences of the mass shift parameters and the line electronic factors are given by

$$\mathsf{\Delta }{F}_{k,n}=F_{u,n}-F_{l,n}=\frac{4\pi Z\mathsf{\Delta }{b}_{k,n}}{(n+2)(n+3)}.$$

(84)

Here, $\mathsf{\Delta }{b}_{k,n}=b_{u,n}-b_{l,n}$ is the difference of the polynomial function coefficients between the upper (u) and the lower (l) level of the transition. For spherical nuclear charge distributions, grasp adopts for ${\rho }^A\left(r\right)$ the Fermi distribution model [73,76].

Using the refined treatment of the field shift, effects of the variation of the electron density (ved) over the nuclear volume as well as atomic energy shifts arising from changes in nuclear charge distributions, including nuclear deformations, can also be estimated. The level field shift (80) can be approximated by

$$\delta E_{\mathsf{\Gamma }J,FS}^{A,A^{\prime }}\approx F_{\mathsf{\Gamma }J,0}^{\mathrm{ved}}\delta {\langle r^2\rangle }^{A,A^{\prime }},$$

(85)

with

$$F_{\mathsf{\Gamma }J,0}^{\mathrm{ved}}=F_{\mathsf{\Gamma }J,0}^{\left(0\right)\mathrm{ved}}+F_{\mathsf{\Gamma }J,0}^{\left(1\right)\mathrm{ved}}\delta {\langle r^2\rangle }^{A,A^{\prime }},$$

(86)

where $F_{\mathsf{\Gamma }J,0}^{\left(0\right)\mathrm{ved}}$ and $F_{\mathsf{\Gamma }J,0}^{\left(1\right)\mathrm{ved}}$ are estimated analytically from the parameters $(a,c)$ defining the Fermi distribution modeling the nuclear charge distribution of the reference isotope and from the four coefficients $F_{\mathsf{\Gamma }J,n}(n=0,2,4,6)$ of Equations (80) and (81). Details are given in the write-up of the relativistic isotope shift program ris4 [68]. For lighter systems, where the electronic density is essentially constant inside the nuclear volume, it is justified to only consider the first electronic factor $\mathsf{\Delta }{F}_{k,0}=\frac{2}{3}\pi Z\mathsf{\Delta }{\rho }_{\mathsf{\Gamma }J}^e\left(\mathbf{0}\right)$, and the nuclear quantity $\delta {\langle r^2\rangle }^{A,A^{\prime }}$ can be extracted from observed line shifts along isotope chains [63] and be directly compared with predictions from nuclear theory [133]. The reduction (31) for the field shift operator is relatively simple thanks to the scalar property of the radial electron density. The spin-angular coefficients are identical to the weighing factors $t_{ab}^{\mathit{\alpha }\beta }$ of the one-electron radial integrals $I(a,b)$ appearing in Equation (41) for the Hamiltonian, which is another $k=0$ tensorial operator. More details can be found in [68]—see also Section 3.4. In the grasp suite of codes, the electronic factors (81) of the level field shift (80), together with $F_{\mathsf{\Gamma }J,0}^{\left(0\right)\mathrm{ved}}$ and $F_{\mathsf{\Gamma }J,0}^{\left(1\right)\mathrm{ved}}$ of (86) are computed by the ris4 program; see §6.1 and §8.1 in the manual.

3.3.3. Total Shift

The total isotope shift in frequency of a given spectral line k is obtained from the differences of the level mass- and field shift electronic factors, which are each weighted by the adequate mass and radial nuclear factors, respectively (see Equation (83)). The auxiliary code fical of the grasp suite of codes allows this calculation. Section 12 of the manual provides a specific example of how to use ris4 together with fical to compute the effect of nuclear deformation on the frequency isotope shift for the $1s^{2}2s{}^2{S}_{1/2}-1s^{2}2p{}^2{P}_{1/2,3/2}^o$ transitions in ${}^{150,142}$Nd${}^{57+}$. Careful attention should be paid to the conventions used to define the relevant ingredients [134].

3.4. Electronic Radial Densities and Natural Orbitals

The radial electron density $D\left(r\right)$ of level $\mathsf{\Gamma }J$ described by the ASF (6) can be calculated from

$$D\left(r\right)=4\pi r^2\rho \left(r\right)=\sum _{\mathit{\alpha },\beta }{\left(c_{\mathit{\alpha }}^{\mathsf{\Gamma }J}\right)}^{*}{D}_{\mathit{\alpha }\beta }\left(r\right)c_{\beta }^{\mathsf{\Gamma }J},$$

(87)

where each density matrix element in the CSF space is given as a weighted sum of radial functions

$$D_{\mathit{\alpha }\beta }\left(r\right)=\left[\sum _{\kappa }\sum _{n^{\prime }n}{v}_{n{n}^{\prime }\kappa }^{\mathit{\alpha },\beta }{I}_{\rho }(n^{\prime }\kappa ,n\kappa ;r)\right],$$

(88)

with

$$I_{\rho }\left(n^{\prime }\kappa ,n\kappa ;r\right)\equiv \left[P_{n^{\prime }\kappa }\left(r\right)P_{n\kappa }\left(r\right)+Q_{n^{\prime }\kappa }\left(r\right)Q_{n\kappa }\left(r\right)\right].$$

(89)

The weighting coefficients $v_{n{n}^{\prime }\kappa }^{\mathit{\alpha },\beta }$ can be estimated from the ${\xi }_{ab;k=0}$ coefficients appearing in the reduction (31) of $\langle \mathit{\alpha }J\parallel {\mathbf{F}}^{\left(0\right)}\parallel \beta J\rangle$ (see [69] for more details). The radial density distribution is normalized to the number of electrons

$${\int }_0^{\infty }D\left(r\right)dr=N.$$

(90)

Although the electron density has to be evaluated within the nuclear volume for FS calculations (see Equation (77)), it is often interesting to probe the electronic radial density distribution through the entire atomic volume to investigate the role of electron correlation and relativity [69,135,136]. For the ground state of neutral beryllium, a good zeroth-order wave function is usually built on the Layzer complex (LC) $(1s^{2}2s^2+1s^{2}2p^2)$ that reflects the near-degeneracy of the two components [137] (see Section 4.2). Electron correlation can be included by building the CSF space, considering electron substitutions from this multireference. The radial electron density calculated with such an MCDHF wave function based on single and double (SD) substitutions from the LC to $n=3$ is displayed (in black) in the left panel of Figure 1. The correlation effect on the density is illustrated through the blue curve corresponding to the difference between the electron densities calculated in this correlation model, $D^{SD}\left(r\right)$ and in the two-configuration zeroth-order approximation, $D^{LC}\left(r\right)$. This difference corresponds to a correlation energy of $0.0376E_{\mathrm{h}}$, i.e., to $\simeq 0.26\%$ of the total binding energy.

The radial electronic density is often viewed as the sum over the orbital densities $I_{\rho }\left(n\kappa ,n\kappa ;r\right)$

$$D\left(r\right)=\sum _{n\kappa }{q}_{n\kappa }{I}_{\rho }\left(n\kappa ,n\kappa ;r\right)$$

(91)

weighted by the total subshell occupation numbers evaluated over the entire CSF space

$$q_{n\kappa }=\sum _{\mathit{\alpha }}{\left|c_{\mathit{\alpha }}^{\mathsf{\Gamma }J}\right|}^2{w}_{n\kappa }^{\mathit{\alpha }},$$

(92)

where $w_{n\kappa }^{\mathit{\alpha }}$ is the occupation number of subshell $n\kappa$ in the CSF $\mathsf{\Phi }\left({\gamma }_{\mathit{\alpha }}J\right)$. This is a naive picture, as demonstrated in the right panel of the same figure, that shows the importance of the off-diagonal matrix elements $D_{\mathit{\alpha },\beta }(\mathit{\alpha }\ne \beta )$ involving off-diagonal contributions $I_{\rho }\left(n^{\prime }\kappa ,n\kappa ;r\right)(n^{\prime }\ne n)$. The latter should never be omitted to obtain reliable results.

The radial density distribution (87) can be rewritten as

$$D\left(r\right)=\sum _{\kappa }\sum _{n^{\prime }n}{\rho }_{n^{\prime }n}^{\kappa }{I}_{\rho }(n^{\prime }\kappa ,n\kappa ;r),$$

(93)

with

$${\rho }_{n^{\prime }n}^{\kappa }=\sum _{\mathit{\alpha },\beta }{\left(c_{\mathit{\alpha }}^{\mathsf{\Gamma }J}\right)}^{*}{v}_{n{n}^{\prime }\kappa }^{\mathit{\alpha },\beta }{c}_{\beta }^{\mathsf{\Gamma }J}.$$

(94)

A “by-product” of the electron density are the natural orbitals (NO)

$$\begin{array}{c}\hfill \left\{\begin{array}{ccc}{\tilde{P}}_{\mathit{\alpha }\kappa }\left(r\right)& =& \sum _n{u}_n^{\mathit{\alpha },\kappa }{P}_{n\kappa }\left(r\right)\\ {\tilde{Q}}_{\mathit{\alpha }\kappa }\left(r\right)& =& \sum _n{u}_n^{\mathit{\alpha },\kappa }{Q}_{n\kappa }\left(r\right),\end{array}\right.\end{array}$$

(95)

that diagonalize the density matrix $\mathbf{\rho }$${}^{\kappa }$

$${\mathbf{U}}^{\kappa †}{\mathbf{\rho }}^{\kappa }{\mathbf{U}}^{\kappa }={\tilde{\mathbf{\rho }}}^{\kappa }.$$

(96)

The NOs have the attractive feature that they concentrate the expansion coefficients to relatively fewer CSFs, resulting in zero weights for others; see [29]. As an example, we take $1s^{2}2s^2{}^{1}S$ in Be I. In

Table 3. Mixing coefficients for the $n=4$ active space valence correlation expansion of the $1s^{2}2s^2{}^1{S}_0$ ground state of Be, using two different orbital bases: the MCDHF-optimized orbitals and their corresponding natural orbitals (MCDHF/NO).

CSF	MCDHF	MCDHF/NO
2s (2)	0.953738	0.953740
2s (1) 3s (1)	−0.001117	0.000000
2s (1) 4s (1)	−0.001846	0.000000
2p (2)	0.242750	0.242750
2p-(2)	0.171674	0.171674
2p (1) 3p (1)	0.000254	0.000000
2p (1) 4p (1)	0.000302	0.000000
2p-(1) 3p-(1)	0.000178	0.000000
2p-(1) 4p-(1)	0.000214	0.000000
3s (2)	−0.039770	−0.039787
3s (1) 4s (1)	−0.001052	0.000000
3p (2)	0.004905	0.004922
3p-(2)	0.003467	0.003479
3p (1) 4p (1)	−0.000333	0.000000
3p-(1) 4p-(1)	−0.000237	0.000000
3d (2)	−0.013120	−0.013134
3d-(2)	−0.010712	−0.010723
3d (1) 4d (1)	0.000530	0.000000
3d-(1) 4d-(1)	0.000432	0.000000
4s (2)	−0.004103	−0.004089
4p (2)	0.001628	0.001611
4p-(2)	0.001150	0.001138
4d (2)	−0.002808	−0.002794
4d-(2)	−0.002291	−0.002280
4f (2)	0.004766	0.004766
4f-(2)	0.004127	0.004127

, we display the expansion coefficients for an RCI calculation based on a CSF expansion built on MCDHF-optimized orbitals up to $n=4$ and accounting for valence–valence electron correlation. Redoing the RCI calculation with CSFs built on NOs gives exactly the same wave function (and energy), but it changes the expansion coefficients, as seen in the right-most column in the table. Now, many of the CSFs have zero expansion coefficients.

The use of NOs is not common in Atomic Physics, unlike in Quantum Chemistry. However, a recent study [138] suggests that natural orbitals might be a promising alternative to the LBL optimization scheme usually adopted in GRASP calculations (see Section 4) to avoid convergence problems. NOs could also prove useful in connection with the treatment of Lagrange multipliers during the orbital optimization process [139]. In grasp, the radial density function as well as the NOs are computed using the rdensity program. The nuclear parameters, the CSF expansion, the radial orbitals, as obtained by rmcdhf, and the expansion coefficients of the CSFs, as obtained by the rmcdhf or rci programs, are read from files. The radial density function is written to file in the format that makes it easy to print. The NOs are written to file in the same format as the radial wave functions from rmcdhf and can be directly used by the rci program. Examples of how to use the rdensity program to produce the radial density function and the NOs are given in §6.8 in the manual.

3.5. Radiative Transition Properties

The atomic states are coupled to the radiation field. Following Grant [73,140], the rate A for spontaneous emission from an upper state ${\mathsf{\Gamma }}^{\prime }{J}^{\prime }{M}_{J^{\prime }}$ to any of the $2J+1$ states $\mathsf{\Gamma }J{M}_J,M_J=-J,-J+1,\dots ,J$ of a lower energy level $\mathsf{\Gamma }J$ is given by

$$A=\frac{2\omega }{c}\frac{1}{(2L+1)(2J^{\prime }+1)}|\langle \mathsf{\Gamma }J\parallel {\mathbf{O}}^{\left(L\right)}\parallel {\mathsf{\Gamma }}^{\prime }{J}^{\prime }{\rangle |}^2,$$

(97)

where ${\mathbf{O}}^{\left(L\right)}={\sum }_{j=1}^N\mathbf{o}{\left(j\right)}^{\left(L\right)}$ is the electric and magnetic multipole transition operator of rank L and $\omega$ is the angular frequency of the transition. Instead of giving the rate, the strength of a transition is often expressed in terms of the dimensionless oscillator strength

$$gf=\frac{1}{\omega }\frac{1}{(2L+1)}|\langle \mathsf{\Gamma }J\parallel {\mathbf{O}}^{\left(L\right)}\parallel {\mathsf{\Gamma }}^{\prime }{J}^{\prime }{\rangle |}^2.$$

(98)

Inserting the CSF expansions for the left and right-hand ASFs, the reduced matrix element of the transition operator is

$$\langle \mathsf{\Gamma }J\parallel {\mathbf{O}}^{\left(L\right)}\parallel {\mathsf{\Gamma }}^{\prime }{J}^{\prime }\rangle =\sum _{\mathit{\alpha },\beta }{c}_{\mathit{\alpha }}^{\mathsf{\Gamma }J}{c}_{\beta }^{{\mathsf{\Gamma }}^{\prime }{J}^{\prime }}\langle {\gamma }_{\mathit{\alpha }}J\parallel {\mathbf{O}}^{\left(L\right)}\parallel {\gamma }_{\beta }^{\prime }{J}^{\prime }\rangle .$$

(99)

The reduced matrix elements between the CSFs can, according to Equation (31), be expressed in terms of spin-angular coefficients and operator strengths as

$$\langle {\gamma }_{\mathit{\alpha }}J\parallel {\mathbf{O}}^{\left(L\right)}\parallel {\gamma }_{\beta }^{\prime }{J}^{\prime }\rangle =\sum _{a,b}{\xi }_{ab;L}\langle n_a{\kappa }_a\left(1\right)\parallel {\mathbf{o}}^{\left(L\right)}\left(1\right)\parallel n_b{\kappa }_b\left(1\right)\rangle .$$

(100)

The one-particle reduced matrix elements above can be further decomposed as

$$\langle n_a{\kappa }_a\parallel {\mathbf{o}}^{\left(L\right)}\parallel n_b{\kappa }_b\rangle =\langle j_a\parallel {\mathbf{C}}^{\left(L\right)}\parallel j_b\rangle M_{a,b}^{e,m}(\omega ;G_L),$$

(101)

where $M_{a,b}^{e,m}(\omega ;G_L)$ is a quantity that can be expressed as a linear combination of three radial integrals involving combinations of $P_a\left(r\right)$, $Q_a\left(r\right)$, $P_b\left(r\right)$, $Q_b\left(r\right)$ weighted by a spherical Bessel function $j_L(\omega r/c)$ (see [73], section 8.2.1). For electric multipoles

$$M_{a,b}^e(\omega ;G_L)=M_{a,b}^e(\omega ;0)+G_L{M}_{a,b}^l\left(\omega \right),$$

(102)

where $M_{a,b}^e(\omega ;0)$ is the Coulomb gauge integral and $M_{a,b}^l\left(\omega \right)$ is the longitudinal part. $G_L$ is the gauge parameter [140,141]. Setting $G_L=0$, we obtain the matrix element in the Coulomb gauge. Choosing $G_L=\sqrt{(L+1)/L}$, we obtain the matrix element in the Babushkin gauge.

The Racah algebra used to express the reduced transition matrix elements appearing in (97) and (98) in terms of spin-angular coefficients and operator strengths assumes that the initial and final atomic states are built from a common and orthonormal orbital set. This is a severe restriction, since a high-quality wave function requires that orbitals be optimized for a specific electronic state. However, for very general initial and final state CSF expansions, the only requirement being that they are closed under orbital de-excitation, built on different and separately optimized orbital sets, it is possible to change the wave function representation in such a way that the standard Racah algebra can be used for the evaluation of the matrix elements in the new representation (see Olsen et al. [142]). The procedure for the calculation of the transition matrix element can be summarized as follows.

1.

Perform MCDHF or CI calculations for the initial and< the final states, where the radial one-electron orbital sets $\left\{(P_{n\kappa },Q_{n\kappa })\right\}$ and $\left\{(P_{n^{\prime }\kappa }^{\prime },Q_{n^{\prime }\kappa }^{\prime })\right\}$ of the lower $\mathsf{\Gamma }J{M}_J$ and upper ${\mathsf{\Gamma }}^{\prime }{J}^{\prime }{M}_{J^{\prime }}$ state wave functions are not assumed to be the same.

2.

For each $\kappa$, compute the radial orbital overlap matrix. Transform the two radial one-electron orbital sets

$$\left\{(P_{n\kappa },Q_{n\kappa })\right\}\to \left\{({\tilde{P}}_{n\kappa },{\tilde{Q}}_{n\kappa })\right\},\left\{(P_{n^{\prime }\kappa }^{\prime },Q_{n\kappa }^{\prime })\right\}\to \left\{({\widetilde{P^{\prime }}}_{n^{\prime }\kappa },{\widetilde{Q^{\prime }}}_{n^{\prime }\kappa })\right\}$$

(103)

to a biorthonormal basis, i.e., a basis such that

$$\langle ({\tilde{P}}_{n\kappa },{\tilde{Q}}_{n\kappa })|({\widetilde{P^{\prime }}}_{n^{\prime }\kappa },{\widetilde{Q^{\prime }}}_{n^{\prime }\kappa })\rangle ={\delta }_{n,n^{\prime }}.$$

(104)

The orbital transformation in effect changes the CSFs, and we have

$$\left\{\mathsf{\Phi }\left({\gamma }_{\mathit{\alpha }}J\right)\right\}\to \left\{\tilde{\mathsf{\Phi }}\left({\gamma }_{\mathit{\alpha }}J\right)\right\},\left\{\mathsf{\Phi }\left({\gamma }_{\beta }^{\prime }{J}^{\prime }\right)\right\}\to \left\{\tilde{\mathsf{\Phi }}\left({\gamma }_{\beta }^{\prime }{J}^{\prime }\right)\right\}.$$

(105)

The orbital transformation is followed by a counter transformation of the CI expansion coefficients

$$\left\{C_{\mathit{\alpha }}^{\mathsf{\Gamma }J}\right\}\to \left\{{\tilde{C}}_{\mathit{\alpha }}^{\mathsf{\Gamma }J}\right\},\left\{C_{\beta }^{{\mathsf{\Gamma }}^{\prime }{J}^{\prime }}\right\}\to \left\{{\tilde{C}}_{\beta }^{{\mathsf{\Gamma }}^{\prime }{J}^{\prime }}\right\}$$

(106)

to leave the total wave functions invariant, i.e.,

$$\sum _{\mathit{\alpha }}{C}_{\mathit{\alpha }}^{\mathsf{\Gamma }J}\mathsf{\Phi }\left({\gamma }_{\mathit{\alpha }}J\right)\equiv \sum _{\mathit{\alpha }}{\tilde{C}}_{\mathit{\alpha }}^{\mathsf{\Gamma }J}\tilde{\mathsf{\Phi }}\left({\gamma }_{\mathit{\alpha }}J\right),\sum _{\beta }{C}_{\beta }^{{\mathsf{\Gamma }}^{\prime }{J}^{\prime }}\mathsf{\Phi }\left({\gamma }_{\beta }^{\prime }{J}^{\prime }\right)\equiv \sum _{\beta }{\tilde{C}}_{\beta }^{{\mathsf{\Gamma }}^{\prime }{J}^{\prime }}\tilde{\mathsf{\Phi }}\left({\gamma }_{\beta }^{\prime }{J}^{\prime }\right)$$

(107)

3.

Use standard Racah algebra to compute the transition matrix elements in the new biorthonormal representation.

In grasp, the biorthogonal transformation is performed with the program rbiotransform. The CSF expansions, required to be closed under de-excitation (see accompanying manual §3.8), the radial orbitals, as obtained by rmcdhf, and expansion coefficients obtained by rmcdhf or rci are read from file. The transformed radial orbitals and expansion coefficients are written back to file. The computation of transition rates and weighted oscillator strengths, in both Babushkin and Coulomb gauges for electric multipoles, are completed with rtransition. The CSF expansions, the transformed radial orbitals and expansion coefficients, along with relevant energies needed to compute the transition frequencies $\omega$, are read from file. The program and file flow leading up to the transition calculations are displayed in Figures 1 and 2 in the manual.

Transition matrix elements in, respectively, the Coulomb and Babushkin gauges, approach matrix elements in the velocity and length forms in the non-relativistic limit. The resulting rates and weighted oscillator strengths for the electric multipole transitions computed in different gauges are known to be different when approximate wave functions are used. It is generally found that parameters in the Babushkin gauge, which puts more weight on the outer parts of the wave functions, are more accurate for excitations from lower lying levels, but those obtained with the Coulomb gauge, which puts more weight on the inner parts of the wave functions, are more accurate for excited states [143].

For highly accurate wave functions, that predict the observed transition energies to fractions of a percent, a statistical analysis by Ekman et al. [144] suggests that the quantity

$$dT=\frac{|A_C-A_B|}{max(A_C,A_B)},$$

(108)

where $A_C$ and $A_B$ are the rates computed in the Coulomb and Babushkin gauges, can be used to estimate the uncertainty of the computed rates. The uncertainty should be interpreted in a statistical sense and applied to groups of transitions [114]. For convenience, the $dT$ parameter is written to file by the rtransition program along with the transition rates and weighted oscillator strengths in both Coulomb and Babushkin gauges. Specific examples of transition calculations are given in Sections 6.1, 6.3, 6.4, 9, 10, and 11 of the manual. In §11.3, the programs are used to study the change of transition rates along an iso-electronic sequence.

3.6. Unexpected Transitions

The transition operator is a one-electron operator, but occasionally, transitions appear to be the result of “two-electron jumps” [14]. An example is the $2s3d{}^{1}D-2p3s{}^1{P}^o$ transition in Be i. This anomaly occurs when there is extensive mixing of the wave function in either the initial or final state. In the present case, the ${}^{1}D$ wave function has a significant $2p^2{}^{1}D$ component that contributes to the electric dipole transition labeled $2s3d{}^{1}D-2p3s{}^1{P}^o$ through the $2p^2{}^{1}D-2p3s{}^1{P}^o$ CSF transition matrix element. This two-electron jump, also classified as a “two-electron-one-photon” (TEOP) process [145], is therefore a consequence of electron correlation in the wave functions describing the initial and/or final levels. Note that TEOP transitions may also appear at the uncorrelated level of approximation [146] if the orbital relaxation resulting from the independent optimization of the initial and final states is large enough to allow significant radial non-orthogonalities. An example is the E1 radiative decay from $1s^{2}2s3p{}^3{P}_0^o$ to $1s^{2}2p^2{}^3{P}_1$ in Be-like ions for which the E1 radial integral $\langle 2s\left|\mathbf{r}\right|2p\rangle$ is multiplied by the non-zero radial overlap $\langle 3p|2p\rangle$ integral involving the spectator subshells. TEOP transitions may also appear for magnetic dipole or electric quadrupole transitions [147].

Other a priori unexpected radiative transitions can be induced by symmetry-breaking perturbations due to, for example, an external magnetic field (magnetic-field-induced transitions, MIT) and/or hyperfine interaction (hyperfine-induced transitions, HIT) and give rise to new lines in the spectrum [148]. Depending on the magnetic field strength, B, and the nuclear spin, I, the total Hamiltonian and the perturbed eigenstates can be expressed as a linear combination of unperturbed ASFs (6) as follows,

$$\begin{array}{ccc}\hfill & & B\ne 0,I=0:\mathcal{H}={\mathcal{H}}_{DCB}+{\mathcal{H}}_m,|\tilde{\mathsf{\Gamma }}{M}_J\rangle =\sum _{\mathsf{\Gamma }J}{d}_{\mathsf{\Gamma }J}|\mathsf{\Gamma }J{M}_J\rangle \hfill \\ & & B=0,I\ne 0:\mathcal{H}={\mathcal{H}}_{DCB}+{\mathcal{H}}_{hfs},|\tilde{\mathsf{\Gamma }}F{M}_F\rangle =\sum _{\mathsf{\Gamma }J}{d}_{\mathsf{\Gamma }J}|\mathsf{\Gamma }IJF{M}_F\rangle \hfill \\ & & B\ne 0,I\ne 0:\mathcal{H}={\mathcal{H}}_{DCB}+{\mathcal{H}}_m+{\mathcal{H}}_{hfs},|\tilde{\mathsf{\Gamma }}I{M}_F\rangle =\sum _{\mathsf{\Gamma }JF}{d}_{\mathsf{\Gamma }JF}|\mathsf{\Gamma }IJF{M}_F\rangle ,\hfill \end{array}$$

(109)

where ${\mathcal{H}}_{DCB}$ is the Dirac–Coulomb–Breit Hamiltonian. ${\mathcal{H}}_m$ and ${\mathcal{H}}_{hfs}$ are given by Equations (63) and (52). Additional quantum numbers are incorporated in the labels of the ASFs as needed for the different cases. The mixing coefficients, d, can be determined through, e.g., perturbation theory or by solving the full eigenvalue problem:

$$\mathbf{H}\mathbf{d}=E\mathbf{d}$$

(110)

where $\mathbf{H}$, taking the $B\ne 0$, $I\ne 0$ case as an example, is the matrix with elements:

$$\begin{array}{ccc}\hfill H_{\mathsf{\Gamma }JF,{\mathsf{\Gamma }}^{\prime }{J}^{\prime }{F}^{\prime }}=\langle \mathsf{\Gamma }IJF{M}_F|{\mathcal{H}}_{DCB}& +& {\mathbf{T}}^{\left(1\right)}·{\mathbf{M}}^{\left(1\right)}+{\mathbf{T}}^{\left(2\right)}·{\mathbf{M}}^{\left(2\right)}\hfill \\ \hfill & +& ({\mathbf{N}}^{\left(1\right)}+\mathsf{\Delta }{\mathbf{N}}^{\left(1\right)})·{\mathbf{B}}^{\left(1\right)}|{\mathsf{\Gamma }}^{\prime }I{J}^{\prime }{F}^{\prime }{M}_F\rangle .\hfill \end{array}$$

(111)

Formulas for the relevant matrix elements can be found in [84]. In grasp, the above matrix elements are computed using the hfszeeman95 program [67]. The nuclear parameters, the CSF expansion as produced by the rcsfgenerate, the radial wave functions as produced by rmcdhf, and the expansion coefficients of the CSFs as obtained by the rmcdhf or rci programs are read from file. The hfszeeman95 program diagonalizes the interaction matrix to give the energies and expansion coefficients in Equation (109). In external magnetic fields, the latter quantities are dependent on B and, given the mithit tool of hfszeeman95, it is possible to map out the energy structure as a function of B beyond the weak Zeeman– or Paschen–Back limits [149].

In the long wavelength approximation, the transition rate for an electric dipole (E1) transition between an ${\tilde{\mathsf{\Gamma }}}^{\prime }{M}_F^{\prime }$ and a lower $\tilde{\mathsf{\Gamma }}{M}_F$ magnetic hyperfine structure substate is given by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill A({\tilde{\mathsf{\Gamma }}}^{\prime }{M}_F^{\prime }\to \tilde{\mathsf{\Gamma }}{M}_F)=\frac{64{\pi }^4{e}^2{a}_0^2}{3h{\lambda }^3}\sum _q|\langle \tilde{\mathsf{\Gamma }}{M}_F|O_q^{\left(1\right)}|{\tilde{\mathsf{\Gamma }}}^{\prime }{M}_F^{\prime }{\rangle |}^2\end{array}\end{array}$$

(112)

where $\lambda =2\pi c/\omega$ is the wavelength and $O_q^{\left(1\right)}$ is the magnetic component of the electric dipole transition operator ${\mathbf{O}}^{\left(1\right)}$. Substituting Equation (109) into Equation (112) and using the WE theorem together with an uncoupling of F, the rate (in s${}^{-1}$) can be rewritten as

$$\begin{array}{c}\begin{array}{c}{\displaystyle A({\tilde{\mathsf{\Gamma }}}^{\prime }{M}_F^{\prime }\to \tilde{\mathsf{\Gamma }}{M}_F)=\frac{2.02613\times {10}^{18}}{{\lambda }^3}\sum _q|\sum _{\mathsf{\Gamma }JF}\sum _{{\mathsf{\Gamma }}^{\prime }{J}^{\prime }{F}^{\prime }}{d}_{\mathsf{\Gamma }JF}{d}_{{\mathsf{\Gamma }}^{\prime }{J}^{\prime }{F}^{\prime }}^{\prime }}\hfill \\ \times \sqrt{(2F^{\prime }+1)(2F+1)}{(-1)}^{F-M_F}\left(\begin{array}{ccc}F& 1& F^{\prime }\\ -M_F& q& M_F^{\prime }\end{array}\right)\hfill \\ \times {(-1)}^{I+J^{\prime }+F+1}\left\{\begin{array}{ccc}J& F& I\\ F^{\prime }& J^{\prime }& 1\end{array}\right\}\langle \mathsf{\Gamma }J\left|\right|{\mathbf{O}}^{\left(1\right)}\left|\right|{\mathsf{\Gamma }}^{\prime }{J}^{\prime }\rangle {|}^2.\hfill \end{array}\hfill \end{array}$$

(113)

where the wavelength $\lambda$ is expressed in ångström. The rates for unexpected symmetry breaking transitions are computed by the mithit program. The program reads matrix elements produced by hfszeeman95, and then, it constructs and diagonalizes the interaction matrix for a specified B value to determine the expansion coefficients in (109). The needed transition matrix elements are computed from the output of the rtransition_phase, which is a modification of the rtransition program, that gives also phase information. In Section §6.9 of the manual, there is an example of how to use the above programs to compute hyperfine and magnetically induced transitions $2s2p{}^3{P}_0^o-2s^2{}^1{S}_0$ in Ni XXV.

4. Selection of CSFs

Starting with the concept of electron correlation, this section aims to guide the user of grasp into the important process of selecting CSFs in order to build the CSF space. This selection ultimately determines the accuracy of computed transition energies and rates as well as other properties.

4.1. Electron Correlation

Electron correlation is one of the most important concepts in computational atomic physics. It can be defined as the effects beyond the single CSF, or Dirac–Fock, approximation. Electron correlation is often divided into static correlation and dynamic correlation [30]. Static correlation is a long-range rearrangement of the single CSF electron charge distribution due to the strong interaction with a set of CSFs that can be formed from configurations built from orbitals with the same principal quantum numbers as the ones that occupy the reference state and where we may think of orbitals with the same principal quantum numbers as being closely degenerate. The set of strongly interacting CSFs is known as the multireference (MR) and includes, in the general case, also important CSFs formed from configurations other than those above. Whereas it is natural to see the MR as a set of important CSFs defining the lowest order approximation of the wave function, the MR is sometimes also thought of as a set of important configurations with the CSFs followed by angular couplings. This is the view taken in the accompanying manual.

Dynamic correlation is a short-range effect that arises from the singularity of the $1/r_{ij}$ electron–electron interaction near points of coalescence where $r_{ij}=0$ and has a wave function cusp condition associated with it [150]. These are not isolated points but include entire regions of space. The more likely regions are those where the probability of finding a pair of electrons is the highest. Dynamic correlation does in part depend on the coupling conditions of the state. In a state where the spins of the electrons are aligned, e.g., $1s2p{}^3{P}_1^o$, the wave function is, by the Pauli exclusion principle, zero whenever ${\mathbf{r}}_1={\mathbf{r}}_2$. By continuity, the probability for the two electrons to closely coincide is small, and the dynamic electron correlation is fairly minor. For $1s2p{}^1{P}_1^o$, the spins do not align, and there is a finite probability for the electrons to coincide, and thus, the dynamic electron correlation is larger.

For many-electron systems, the largest contributions to electron correlation come from pairs of electrons which occupy the same region in space. Thus, there are large contributions from each doubly occupied orbital, with smaller additions from orbital pairs that occupy different shells. Just as for the static correlation, the dynamic correlation can be accounted for by expansions over CSFs, and the effect should be to mimic the cusp behavior of the exact wave function at points of electron coalescence.

4.2. Z-Dependent Perturbation Theory

Static and dynamic correlation can be understood in terms of Z-dependent perturbation theory, which, for simplicity, we apply in non-relativistic theory. Introduce a new variable $\rho =Zr$ in which the Hamiltonian becomes

$$\mathcal{H}=Z^2\left({\mathcal{H}}^{\left(0\right)}+Z^{-1}V\right),$$

(114)

where

$${\mathcal{H}}^{\left(0\right)}=\sum _{i=1}^N\left(-\frac{1}{2}{\nabla }_i^2-\frac{1}{{\rho }_i}\right),V=\sum _{i>j}^N\frac{1}{{\rho }_{ij}}.$$

(115)

Schrödinger equation now reads

$$\left({\mathcal{H}}^{\left(0\right)}+Z^{-1}V\right)\mathsf{\Psi }=\left(Z^{-2}E\right)\mathsf{\Psi }.$$

(116)

In this form, $1/Z$ appears as a perturbation parameter. Expanding the wave function

$$\mathsf{\Psi }={\mathsf{\Psi }}^{\left(0\right)}+Z^{-1}{\mathsf{\Psi }}^{\left(1\right)}+Z^{-2}{\mathsf{\Psi }}^{\left(2\right)}+\cdots$$

(117)

and the energy

$$E=Z^2\left(E^{\left(0\right)}+Z^{-1}{E}^{\left(1\right)}+Z^{-2}{E}^{\left(2\right)}+Z^{-3}{E}^{\left(3\right)}+\cdots \right)$$

(118)

in terms of $1/Z$ and inserting in Equation (116) gives

$$\begin{array}{ccc}\hfill ({\mathcal{H}}^{\left(0\right)}-E^{\left(0\right)}){\mathsf{\Psi }}^{\left(0\right)}& =& 0\hfill \\ \hfill ({\mathcal{H}}^{\left(0\right)}-E^{\left(0\right)}){\mathsf{\Psi }}^{\left(1\right)}& =& (E^{\left(1\right)}-V){\mathsf{\Psi }}^{\left(0\right)}\hfill \\ \hfill ({\mathcal{H}}^{\left(0\right)}-E^{\left(0\right)}){\mathsf{\Psi }}^{\left(2\right)}& =& (E^{\left(1\right)}-V){\mathsf{\Psi }}^{\left(1\right)}+E^{\left(2\right)}{\mathsf{\Psi }}^{\left(0\right)}.\hfill \end{array}$$

(119)

The solutions of the first equation are products of hydrogenic orbitals.

Let $\left|\right\{nl\}\gamma LS\rangle$ be a configuration state function constructed from products of hydrogenic orbitals. Here, $\left\{nl\right\}=\{n_1{l}_1,n_2{l}_2,\dots ,n_N{l}_N\}$ is the set of N principal and orbital quantum numbers that define the configuration (14) and $\gamma$ denotes the complete set of the coupling tree quantum numbers specifying unambiguously the considered configuration state. Then

$${\mathcal{H}}^{\left(0\right)}|\left\{nl\right\}\gamma LS\rangle =E^{\left(0\right)}|\left\{nl\right\}\gamma LS\rangle$$

(120)

with

$$E^{\left(0\right)}=-\sum _{i=1}^N\frac{1}{2n_i^2}.$$

(121)

Since $E^{\left(0\right)}$ is independent of the orbital quantum numbers, we infer that different configurations may lead to the same energy, i.e., $E^{\left(0\right)}$ is degenerate. According to first-order perturbation theory for degenerate states, ${\mathsf{\Psi }}^{\left(0\right)}$ is a linear combination of the degenerate configuration state functions $|\left\{n{l}^{\prime }\right\}{\gamma }^{\prime }LS\rangle$

$${\mathsf{\Psi }}^{\left(0\right)}=\sum _{l^{\prime }{\gamma }^{\prime }}{c}_{l^{\prime }{\gamma }^{\prime }}|\left\{n{l}^{\prime }\right\}{\gamma }^{\prime }LS\rangle .$$

(122)

The expansion coefficients are components of an eigenvector of the interaction matrix $\langle \left\{n{l}^{\prime }\right\}{\gamma }^{\prime }LS\left|V\right|\left\{nl\right\}\gamma LS\rangle$, with $E^{\left(1\right)}$ as the corresponding eigenvalue. However, only configurations with the same parity interact, and the linear combination is over all CSFs with the same set of principal quantum numbers and the same parity. This set of CSFs is referred to as the complex by Layzer [137]. A relativistic generalization of the non-relativistic Z-dependent perturbation theory has been derived for many-electron atoms [151], but discussing the zero-/first-order structure of the non-relativistic wave function is enough for our purpose.

The first-order correction ${\mathsf{\Psi }}^{\left(1\right)}$ is a solution of Equation (119) orthogonal to ${\mathsf{\Psi }}^{\left(0\right)}$. It can be expanded as a linear combination of normalized intermediate configuration state functions $|{\gamma }_{v}LS\rangle$ belonging to ${\mathcal{H}}^{\left(0\right)}$ but outside the complex. Then

$${\mathsf{\Psi }}^{\left(1\right)}=\sum _v\frac{|{\gamma }_{v}LS\rangle \langle {\gamma }_{v}LS\left|V\right|{\mathsf{\Psi }}^{\left(0\right)}\rangle }{E^{\left(0\right)}-E_{{\gamma }_{v}LS}},$$

(123)

where $E_{{\gamma }_{v}LS}=\langle {\gamma }_{v}LS|{\mathcal{H}}^{\left(0\right)}|{\gamma }_{v}LS\rangle$. Substituting Equation (122) into (123) and interchanging the orders of summation, we have

$${\mathsf{\Psi }}^{\left(1\right)}=\sum _{l^{\prime }{\gamma }^{\prime }}{c}_{l^{\prime }{\gamma }^{\prime }}\sum _v\frac{|{\gamma }_{v}LS\rangle \langle {\gamma }_{v}LS\left|V\right|\left\{n{l}^{\prime }\right\}{\gamma }^{\prime }LS\rangle }{E^{\left(0\right)}-E_{{\gamma }_{v}LS}}.$$

(124)

The mixing coefficient, $c_{l^{\prime }{\gamma }^{\prime }}$, are thus weight factors in the sum over intermediate configuration state functions $|{\gamma }_{v}LS\rangle$ interacting (i.e., having non-zero matrix elements) with configuration state functions in the complex.

4.3. Classification of Correlation Effects

The zero-order wave function ${\mathsf{\Psi }}^{\left(0\right)}$ is obtained as a linear combination of CSFs in the complex. It gives a gross description of the system, and it accounts for the major part of the long-range static electron correlation. The first-order correction ${\mathsf{\Psi }}^{\left(1\right)}$ is a linear combination of CSFs that interacts with the CSFs in the complex. It accounts for additional long-range electron correlation and the major part of the short-range dynamic correlation.

Assume for simplicity that there is only one configuration state function $\left|\right\{nl\}\gamma LS\rangle$ in the complex. Noting that V is a two-body operator, we infer that CSFs interacting with $\left|\right\{nl\}\gamma LS\rangle$ are of two types: those that differ by a single electron (single substitution S) and those that differ by two electrons (double substitution D). The former can be further subdivided into

1.

Those that differ from $\left|\right\{nl\}\gamma LS\rangle$ by one principal quantum number, but retain the same spin and orbital angular coupling. These configuration states are part of radial correlation.

2.

Those that differ by one principal quantum number and also differ in their coupling. Often the only change is the coupling of the spins, in which case the configuration states are part of spin-polarization.

3.

Those that differ in the angular momentum of exactly one electron and are accompanied by a change in orbital angular coupling of the configuration state and possibly also the spin coupling. These represent orbital polarization.

The sums over CSFs that differ in two electrons can also be classified. Let $\{a,b,c,..\}$ be occupied orbitals in $\left|\right\{nl\}\gamma LS\rangle$ and $\{v,v^{\prime },..\}$ be orbitals in a so-called active set. Then, the double substitution $ab\to v{v}^{\prime }$ (also called “replacement”), generates CSFs in the expansion for ${\mathsf{\Psi }}^{\left(1\right)}$. The function defined by CSFs from all double replacements from $ab$ is called a pair-correlation function (PCF), and it corrects for the cusp in the wave function associated with this electron pair. The PCFs from all electron pairs correct for the main part of the dynamic correlation. There is another and more general classification which distinguishes orbital replacements from valence and from core orbitals:

1.

If $ab$ are orbitals for outer electrons, the replacement represents outer or valence correlation.

2.

If a is a core orbital but b is an outer orbital, the effect represents the polarization of the core and is referred to as core–valence correlation.

3.

If both orbitals are from the core, the replacement represents core–core correlation.

4.4. The Active Set Approach

Z-dependent perturbation theory is not appropriate for practical calculations, but it provides a very useful guide for how the initial DF approximation can be improved in MCDHF or CI calculations in order to capture most of the correlation energy. The zero-order wave function ${\mathsf{\Psi }}^{\left(0\right)}$ accounting for the major part of the static correlation is an expansion over CSFs with large interactions with the CSF of interest: either those that are nearly degenerate or those with a large interaction matrix element. These CSFs define the MR set. As is discussed in the manual §4.2, the selection of the MR is far from trivial and requires a number of exploratory calculations. In addition, in order to account for dynamic correlation, the wave function $\mathsf{\Psi }$ should include CSFs generated by SD substitutions of orbitals from each CSF of the MR set, with orbitals in an active set (AS) in a process referred to as the active set approach.

As an example, we look at the $1s^{2}2s^{2}2p^{6}3s^2{}^1{S}_0$ ground state and the $1s^{2}2s^{2}2p^{6}3s3p$${}^3{P}_{0,1,2}^o$ and ${}^1{P}_1^o$ excited states in Mg I. Starting with the ground state, the $3s^2$, $3p^2$ and $3d^2$ configurations (in non-relativistic notation) are formed by orbitals with the same principal quantum numbers, and these configurations are closely degenerate. An MR for the ground state could consist of the CSFs that can be formed from these three configurations. If we now turn to the excited states, the $3s3p$ and $3p3d$ configurations are formed by orbitals with the same principal quantum numbers, and these configurations are closely degenerate. An MR in this case could consist of the CSFs that can be formed from these configurations. However, it turns out that $3s4p$ is important, and thus, a more suitable MR should consist of CSFs also from the latter. The CSFs in the MR and the corresponding expansion coefficients are shown in

Table 4. CSFs and corresponding expansion coefficients for the MR of the $1s^{2}2s^{2}2p^{6}3s^2{}^1{S}_0$ ground state and the $1s^{2}2s^{2}2p^{6}3s3p{}^3{P}_1^o$, ${}^1{P}_1^o$ excited states in Mg I.

CSFs	$3s^2{}^1{S}_0$
1s (2) 2s (2) 2p-(2) 2p (4) 3s (2)	0.964240
1s (2) 2s (2) 2p-(2) 2p (4) 3p (2)	0.214715
1s (2) 2s (2) 2p-(2) 2p (4) 3p-(2)	0.152334
1s (2) 2s (2) 2p-(2) 2p (4) 3d (2)	−0.023696
1s (2) 2s (2) 2p-(2) 2p (4) 3d-(2)	−0.019299
CSFs	$3s3p{}^3{P}_1^o$	$3s3p{}^1{P}_1^o$
1s (2) 2s (2) 2p-(2) 2p (4) 3s (1) 3p-(1)	0.811480	0.756818
1s (2) 2s (2) 2p-(2) 2p (4) 3s (1) 3p (1)	−0.571905	0.532680
1s (2) 2s (2) 2p-(2) 2p (4) 3p (1) 3d-(1)	0.082429	−0.233472
1s (2) 2s (2) 2p-(2) 2p (4) 3p (1) 3d (1)	−0.061506	0.191522
1s (2) 2s (2) 2p-(2) 2p (4) 3p-(1) 3d-(1)	0.046401	−0.166818
1s (2) 2s (2) 2p-(2) 2p (4) 3s (1) 4p-(1)	0.032508	0.142870
1s (2) 2s (2) 2p-(2) 2p (4) 3s (1) 4p (1)	−0.025260	0.063628

for the ground state and the two excited $J=1$ states, i.e., the $1s^{2}2s^{2}2p^{6}3s3p{}^3{P}_1^o,{}^1{P}_1^o$ states in $LSJ$-notation. One thing to note is that the excited states are impure in $jj$-coupling, with heavy mixing between the $|3s3p-\rangle$ and $|3s3p\rangle$ CSFs. Another thing to note is that there is more electron correlation in $3s3p{}^1{P}_1^o$, with non-aligned spins, than in $3s3p{}^3{P}_1^o$, with aligned spins (compare Section 4.1), which is reflected in the larger mixing coefficients for the CSFs in the MR for the latter.

Valence correlation is accounted for by considering CSFs obtained from double substitutions from the outer orbitals of each of the CSFs in the MR to orbitals in an active set of orbitals. Core–valence correlation with the $2s^{2}2p^6$ shell is accounted for by considering CSFs obtained from a single excitation from $2s^{2}2p^6$ and a single substitution from one of the outer electrons in each of the CSFs in the MR to orbitals in an active set of orbitals. Core–core correlation within the $2s^{2}2p^6$ shell is accounted for by considering CSFs obtained from double excitations from the $2s,2p$ orbitals of each of the CSFs in the MR to orbitals in an active set of orbitals. All the generated CSFs should interact with at least one CSF in the MR. Included in the general expansions above are also the CSFs obtained by single replacements, although there is no clear classification in valence, core–valence or core correlation effects. In the active set approach, the active set of orbitals is systematically extended in a layer-by-layer fashion until some sort of convergence of the computed properties is achieved.

For energies, higher-order corrections are captured by including CSFs that interact with the CSFs in the zero- and first-order wave function. In practice, this is the same as including some CSFs that can be generated by triple and quadruple (TQ) orbital substitutions from the orbitals of each of the CSFs in the MR to an active set. The number of CSFs generated in this way grows very rapidly with the size of the active set, and for this reason, one has to limit the size of the latter. A way to include the most important higher-order correlation effects is to increase the MR set by adding CSFs for a certain portion $p={\sum }_{\mathit{\alpha }\subset MR}{c}_{\mathit{\alpha }}^2$ of the wave function composition [136]. The overall accuracy of the wave function increases as the MR set accounts for a larger portion of the wave function.

4.5. CSF Expansions for Energy Differences

We are often interested in determining energy separations between different levels. In these cases we may, in the first approximation, define closed inner subshells as inactive and consider correlation only between the outer valence electrons. The rationale for this is that the correlation energy in the core, although large in an absolute sense, to a great extent cancels when computing energy level differences or the energy relative to the ground state. However, the presence of outer valence electrons polarizes the core. The effect of this polarization is represented by core–valence correlation that increases the binding of the valence electrons to the core. In the case of a single electron, this is reflected in a contraction of the orbital, which has a large effect on other computed properties. Generally, energy separations are much improved if core–valence correlation is included. For larger atomic systems it is not always clear which subshells should be inactive and which should be part of the active core, for which core–valence effects are to be considered. For each new system, this needs to be systematically investigated. A good starting point for analyzing the situation is to plot the radial parts of the core and valence orbitals and look at the overlap between the different orbitals. If the overlap is large, then core–valence effects are likely to be important. In grasp, the rwfnplot program is used for plotting, and examples are given in §7.3 of the manual.

4.6. The Active Set Approach as Implemented in GRASP

In grasp, lists of CSFs are generated using the active set approach as implemented in the rcsfgenerate program. The user is required to specify the configurations (non-relativistic notation) in the MR along with information whether the orbitals in the configurations are inactive (i), i.e., no substitutions are allowed, active (*), i.e., unrestricted substitutions are allowed, or have a minimal occupation (m). To generate a valence expansion for the ground state in Mg, the specification would be

   1s(2,i)2s(2,i)2p(6,i)3s(2,*)

Given an orbital set, the number of substitutions, SD in this case, and the resulting J value, a set of CSFs is generated. To generate an expansion accounting for valence and core–valence correlation with $2s^{2}2p^6$, the specification is

   1s(2,i)2s(2,i)2p(6,5)3s(2,*)

   1s(2,i)2s(2,1)2p(6,i)3s(2,*)

   1s(2,i)2s(2,i)2p(6,5)3p(2,*)

   1s(2,i)2s(2,1)2p(6,i)3p(2,*)

   1s(2,i)2s(2,i)2p(6,5)3d(2,*)

   1s(2,i)2s(2,1)2p(6,i)3d(2,*)

Again, given an orbital set, the number of substitutions (SD in this case), and the resulting J value, a set of CSFs is generated. A detailed description of the rcsfgenerate program is given in Sections 4.3 and 4.4, as well as in Sections 5.1–5.10 of the manual.

The rcsfgenerate program gives a number of configurations as determined by the substitution rules, the active set of orbitals and the specified J values. Given these configurations, a set of CSFs follows by applying the angular couplings described in Section 2.4. Not all of these CSFs interact with the CSFs formed by the configurations in the MR. To include only CSFs that interact, one has to run the rcsfinteract program. It should be remembered that the Z-dependent perturbation theory underlying the active set approach is applicable mainly for ionized systems. For neutral and near neutral systems, it may be necessary to include all generated CSFs and not only the ones interacting with the CSFs formed by the configurations in the MR; see [152] for a discussion. Examples of how to use the rcsfinteract program are given in the manual §5.5.

5. Examples of Applications

We illustrate in the present section the use of the grasp suite of programs through a few specific examples of atomic structure calculations. As explained in Section 2.8, the correlation models used for the MCDHF orbital optimization may differ from those of the subsequent CI calculations performed by the rci program. For each considered case, the details of the computational strategies can be found in the original publications.

5.1. Determination of the Nuclear Quadrupole Moment $Q{(}^{67}$Zn)

One of the most accurate methods to determine nuclear quadrupole moments Q is to combine measured nuclear quadrupole hyperfine interaction constants $B_{\mathsf{\Gamma }J}$ with calculated reduced electric matrix elements [124]. We illustrate the method by determining improved nuclear quadrupole moments of zinc, which is the second most abundant essential trace element in the human body, after iron. For zinc, the standard value cited in the 2016 review of Stone [127] was $Q{(}^{67}Zn)=0.150\left(15\right)$ b, which was derived from measurements on the $4s4p{}^3{P}_{1,2}^o$ states, where magnetic hyperfine interaction constants were used to estimate the reduced electric matrix elements [153]. The value above can be improved by combining measurements by Byron et al. and Lurio [154,155] with electronic factors from large-scale atomic calculations.

Table 5. A (MHz), $B/Q$ (MHz/b), and Q (in barn) values as functions of the increasing active set of orbitals for the $4s4p{}^3{P}_1^o$ and $4s4p{}^3{P}_2^o$ states in ${}^{67}$Zn i. $I^{\pi }=5/2^{-}$ and ${\mu }_{\mathrm{expt}}=0.875479\left(9\right){\mu }_N$. The Q-values are extracted from the relation $Q=B_{\mathrm{expt}}/(B/Q)$, where the experimental values are $B_{\mathrm{expt}}{(}^3{P}_1^o)=-18.782\left(8\right)$¹ MHz and $B_{\mathrm{expt}}{(}^3{P}_2^o)=35.806\left(5\right)$² MHz.

Active Set	$4\mathit{s}4\mathit{p}{}^3{\mathit{P}}_1^{\mathit{o}}$				$4\mathit{s}4\mathit{p}{}^3{\mathit{P}}_2^{\mathit{o}}$
	${\mathit{N}}_{\mathbf{CSFs}}$	A (MHz)	$\mathit{B}/\mathit{Q}$ (MHz/b)	Q (barn)	${\mathit{N}}_{\mathbf{CSFs}}$	A (MHz)	$\mathit{B}/\mathit{Q}$ (MHz/b)	Q (barn)
MCDHF-SrDT-SP
DHF	2	473.40	$-100.373$	0.1098	1	$419.93$	192.924	0.1159
VV+CV
$5s5p4d4f$	$1592$	$558.02$	$-131.036$	$0.1433$	$2122$	$483.71$	$254.975$	$0.1404$
$6s6p5d5f5g$	$11932$	$590.45$	$-146.084$	$0.1286$	$16961$	$507.74$	$280.708$	$0.1276$
$7s7p6d6f6g6h$	$48574$	$610.80$	$-150.997$	$0.1244$	$71610$	$529.87$	$290.233$	$0.1234$
$8s8p7d7f7g7h$	$128264$	$613.17$	$-152.617$	$0.1231$	$191495$	$532.46$	$292.535$	$0.1220$
$9s9p8d8f8g8h$	$267998$	$617.02$	$-154.391$	$0.1217$	$402586$	$536.97$	$296.441$	$0.1208$
Liu et al. [157]		$605.9$	$-150.7$	$0.1247$
Expt.		$609.086\left(2\right)$1				$531.987\left(5\right)$2

¹ Byron et al. [154]. ² Lurio [155].

shows the magnetic and electric hyperfine interaction constants from systematic MCDHF and CI calculations by Bieroń et al. [156]. The Dirac–Hartree–Fock (DHF) wave functions consist, for the $J=1$ and $J=2$ states, of two and one CSFs, respectively, and they are the simplest approximation of the atomic system. The larger calculations result from a single reference model in which configurations are generated by single, double and triple substitutions from $3p^{6}3d^{10}4s4p$ to yield $3p^{6}3d^{9}nl{n}^{\prime }{l}^{\prime }{n}^{\prime \prime }{l}^{\prime \prime }$ and $3p^{5}3d^{10}nl{n}^{\prime }{l}^{\prime }{n}^{\prime \prime }{l}^{\prime \prime }$, in non-relativistic notation, with $nl,n^{\prime }{l}^{\prime },n^{\prime \prime }{l}^{\prime \prime }$ in systematically larger orbital sets, which are denoted by the number of orbitals of the same symmetry. These configurations are augmented by configurations formed from restricted single substitutions from all core subshells with s-symmetry. We see that the calculated magnetic dipole interaction constants, A, are substantially improved and, for the largest orbital set, they are quite close to the experimental values. The quantities $B/Q$ show reasonable convergence patterns with respect to the increasing orbital set, and the extracted Q from the combined calculations and experiments for the $J=1$ and $J=2$ states agree very well. The values in Table 5, along with a number of other calculations, probing other electron correlation effects such as core–core correlation deep down in the atomic core, allowed Bieroń et al. [156] to give an updated and improved value of the quadrupole moment along with an uncertainty estimate $Q{(}^{67}$Zn) = 0.122(10) b. The calculation by Bieroń et al. illustrates the power of combined measurements and calculated electronic quantities to extract nuclear data.

5.2. Determination of Changes in Nuclear Radii

Isotope shift (IS) measurements are valuable sources for information about changes in the nuclear charge radii and distributions along isotope sequences [158,159]. Adding one or more neutrons makes the nucleus not only heavier but may also modify its charge distribution, sometimes considerably [160]. During the last decade, optical laser spectroscopy has emerged as an important alternative to earlier X-ray measurements in muonic atoms or electron-scattering experiments [161]. Using optical laser spectroscopy, measurements became possible for quite long sequences of stable and radioactive isotopes.

However, the determination of the electronic part of the IS through ab initio atomic calculations is still a challenge due to the large cancellation between the upper and lower level mass shifts, with a resulting loss of significant digits for the differences, $\mathsf{\Delta }{K}_{NMS}$ and $\mathsf{\Delta }{K}_{SMS}$, of the level mass shift parameters. For the differences to be accurate, the parameters for the upper and lower states must be very well converged with respect to the increasing active set of orbitals. Furthermore, the result is also sensitive to electron correlation deep down in the atomic core and higher-order correlation effects [162]. The differences of the level mass shift parameters can be theoretically attained with relatively high precision and reliability for few-electron ions [163], but it is desirable to have more accurate results for the larger neutral or near-neutral atoms and ions. To illustrate the usefulness of IS measurements combined with atomic structure calculations as an interface to nuclear physics, we follow Li et al. [164] and look at the $2s{}^2{S}_{1/2}-2p{}^2{P}_{1/2}^o$ and $2s{}^2{S}_{1/2}-2p{}^2{P}_{3/2}^o$ transitions in Nd${}^{57+}$. Computed mass shift differences and field shift factors from MCDHF and CI calculations are shown in

Table 6. Relativistic mass shift $\mathsf{\Delta }{K}_{MS}$ (in GHz u) and field shift F (in MHz/fm${}^2$) factors for the $2s{}^2{S}_{1/2}-2p{}^2{P}_{1/2}^o$ and $2s{}^2{S}_{1/2}-2p{}^2{P}_{3/2}^o$ transitions in Nd${}^{57+}$ from Li et al. [164]. For comparison, individual relativistic normal mass shift $\mathsf{\Delta }{K}_{NMS}$ and specific mass shift $\mathsf{\Delta }{K}_{SMS}$ coefficients (in GHz u) are also included. The results in the second row, labeled MCDHF, were obtained in the MCDHF model, with the largest size of the active set ($n=5$). The numbers in square brackets represent powers of 10.

Model	$2\mathit{s}{}^2{\mathit{S}}_{1/2}-2\mathit{p}{}^2{\mathit{P}}_{1/2}^{\mathit{o}}$				$2\mathit{s}{}^2{\mathit{S}}_{1/2}-2\mathit{p}{}^2{\mathit{P}}_{3/2}^{\mathit{o}}$
	$\mathsf{\Delta }{\mathit{K}}_{\mathit{NMS}}$	$\mathsf{\Delta }{\mathit{K}}_{\mathit{SMS}}$	$\mathsf{\Delta }{\mathit{K}}_{\mathit{MS}}$	$\mathit{F}$	$\mathsf{\Delta }{\mathit{K}}_{\mathit{NMS}}$	${\mathit{K}}_{\mathit{SMS}}$	$\mathsf{\Delta }{\mathit{K}}_{\mathit{MS}}$	$\mathit{F}$
DHF	−1083[1]	−8227[2]	−8336[2]	−7903[3]	−8721[1]	−8768[2]	−9640[2]	−8215[3]
MCDHF ($n=5$)	−1342[1]	−8196[2]	−8331[2]	−7929[3]	−8589[1]	−8761[2]	−9620[2]	−8203[3]
CI + Breit	−1449[1]	−8196[2]	−8341[2]	−7885[3]	−8577[1]	−8775[2]	−9632[2]	−8157[3]
Kozhedub et al. [165]	−1641.8[1]	−8180.90[2]	−8345.08(25)[2]		−8573.3[1]	−8769.29[2]	−9626.62(25)[2]
Brandau et al. [166]				−7520[3]				−7810[3]

. The MCDHF calculation includes CSFs that can be formed from configurations from single, double, and triple substitutions to an active set consisting of orbitals with principal quantum numbers up to $n=5$, at which level a high degree of convergence of the computed parameters has been attained. The CI calculation includes the Breit interaction, which affects mainly the normal mass shift of the $2s{}^2{S}_{1/2}-2p{}^2{P}_{1/2}^o$ transition. Included in Table 6 are results from Kozhedub et al. [165]. Although the individual contributions from the normal and specific mass shifts are somewhat different, especially for the $2s{}^2{S}_{1/2}-2p{}^2{P}_{1/2}^o$ transition, the total mass shifts are very consistent. The table also includes experimentally extracted F values by Brandau et al. [166]. The lack of agreement between theory and experiment is, as discussed in [164], due to neglected effects from the variation of the electron density over the nuclear volume as well as effects from the deformation of the nuclei. Through a reformulation of the field shift, the latter effects can now be reliably computed by the isotope shift program [68]; see §12.2 in the accompanying manual. Adding these effects, the agreement between theory and experiment is highly satisfactory, see Table 6 of Ekman et al. [68], demonstrating that differences in nuclear radii can indeed be accurately extracted by combining experimental isotope shift with computed mass shift parameters.

5.3. Spectroscopic Data for Astrophysics—Al-like Ions

Large-scale multiconfiguration calculations have been shown to be capable of supplying transition data of spectroscopic accuracy for a wide range of astrophysically important ions [46]. Spectroscopic accuracy, in this context, means that the calculations are capable of supplying transition energies with an accuracy high enough to directly aid line identification in observed astrophysical or experimental spectra or to weed out incorrect line assignments. Furthermore, the calculations provide transition rates, with uncertainties for stronger transitions of just a few percent. We demonstrate the capacity of these types of calculations in the case of Al-like iron, Fe XIV, which is an important ion for diagnosing the solar corona [167]. In active region spectra, or irradiance spectra of the active Sun, Fe XIV lines within the $n=3$ complex are amongst the most prominent ones in the extreme ultraviolet (EUV) region of the spectrum. In such cases, Fe XIV lines provide the best density diagnostics in the EUV. Transitions from the $n=4$ levels are instead prominent in the soft X-rays, where accurate theoretical energies are needed to identify the lines in the observed spectra.

The calculations by Ekman et al. [167] targeted the 360 lowest states belonging to the 30 configurations $3s^2\{3l,4l,5l\}$, $3p^2\{3d,4l\}$, $3s\{3p^2,3d^2\}$, $3s\{3p3d,3p4l,3p5s,3d4l^{\prime }\}$, $3p3d^2$, $3p^3$ and $3d^3$ with $l=0,1,\dots ,n-1$ and $l^{\prime }=0,1,2$. The starting point for the calculations was an MR set including the 16 even parity configurations $3s^2\{3sd,4sd,5sdg\}$, $3s\{3p^2,3d^2,3p4p,3p4f,3d4s\}$, $3p^2\{3d,4sd\}$ and $3d^3$ and the 14 odd parity $3s^2\{3p,4pf,5pf\}$, $3s3p\{3d,4sd,5s\}$, $3s3d4p$, $3p^3$, $3p^2\left\{4pf\right\}$ and $3p3d^2$. Configurations, from which the CSFs are formed, were then obtained by SD substitutions from occupied subshells of the configurations in the MR to orbitals in an increasing active set. No substitutions were allowed from the $1s$ subshell, and at most, one substitution was allowed from $n=2$ subshells. The first three layers in the active set contain orbitals up to $nl=5g,6h$ and $7i$, respectively. Subsequent calculations included only SD substitutions from the $n>2$ valence subshells up to $nl=8k$ and $9k$ active orbitals. The resulting expansions consisted of approximately 2 640 000 and 2 450 000 CSFs distributed over the $J=1/2,3/2,\cdots ,11/2$ angular symmetries for even and odd parity, respectively. The wave function expansions were in terms of $jj$-coupled CSFs. To adhere to the labeling used by the astrophysics community, the wave function was transformed to a representation in $LSJ$-coupling. The excitation energies for the lowest 40 states from the final CI calculation are shown in

Table 7. Energies in cm${}^{-1}$ and $LS$-compositions for the first 40 levels in Fe XIV. $E_{CI}$ Ekman et al. [167], $E_{NIST}$ NIST Atomic Spectra Database (2013) [168] and $\mathsf{\Delta }E$ difference between $E_{CI}$ and $E_{NIST}$. Indices “a” and “b” are used to differentiate between identical configurations which share the same coupling and leading $LS$-percentage composition. The first number in the $LS$-compositions is the expansion coefficient for the leading configuration and $LSJ$ term in column 3.

No.	Level	$\mathbf{LS}$-Composition	${\mathit{E}}_{\mathbf{CI}}$	${\mathit{E}}_{\mathbf{NIST}}$	$\mathbf{\Delta }\mathit{E}$
1	$3s^{2}3p{}^2{P}_{1/2}^{\circ }$	0.97	0	0	0
2	$3s^{2}3p{}^2{P}_{3/2}^{\circ }$	0.97	18 855	18 853	2
3	$3s3p^2{(}_2^{3}P){}^4{P}_{1/2}$	0.98	225 086	225 114	−28
4	$3s3p^2{(}_2^{3}P){}^4{P}_{3/2}$	0.99	232 777	232 789	−12
5	$3s3p^2{(}_2^{3}P){}^4{P}_{5/2}$	0.97	242 372	242 387	−15
6	$3s3p^2{(}_2^{1}D){}^2{D}_{3/2}$	0.86 + 0.11 $3s^{2}3d{}^{2}D$	299 402	299 242	160
7	$3s3p^2{(}_2^{1}D){}^2{D}_{5/2}$	0.85 + 0.11 $3s^{2}3d{}^{2}D$ + 0.02 $3s3p^2{(}_2^{3}P){}^{4}P$	301 627	301 469	158
8	$3s3p^2{(}_0^{1}S){}^2{S}_{1/2}$	0.75 + 0.21 $3s3p^2{(}_2^{3}P){}^{2}P$	364 945	364 693	252
9	$3s3p^2{(}_2^{3}P){}^2{P}_{1/2}$	0.75 + 0.22 $3s3p^2{(}_0^{1}S){}^{2}S$	388 711	388 510	201
10	$3s3p^2{(}_2^{3}P){}^2{P}_{3/2}$	0.95 + 0.02 $3p^2{(}_2^{1}D){}^{1}D3d{}^{2}P$	396 687	396 512	175
11	$3s^{2}3d{}^2{D}_{3/2}$	0.86 + 0.11 $3s3p^2{(}_2^{1}D){}^{2}D$ + 0.02 $3p^2{(}_0^{1}S){}^{1}S3d{}^{2}D$	473 231	473 223	8
12	$3s^{2}3d{}^2{D}_{5/2}$	0.86 + 0.11 $3s3p^2{(}_2^{1}D){}^{2}D$ + 0.02 $3p^2{(}_0^{1}S){}^{1}S3d{}^{2}D$	475 215	475 202	13
13	$3p^3{(}_3^{2}D){}^2{D}_{3/2}^{\circ }$	0.64 + 0.27 $3s3p{}^{3}P3d{}^2{D}^{\circ }$ + 0.04 $3p^3{(}_1^{2}P){}^2{P}^{\circ }$	576 599	576 383	216
14	$3p^3{(}_3^{2}D){}^2{D}_{5/2}^{\circ }$	0.69 + 0.29 $3s3p{}^{3}P3d{}^2{D}^{\circ }$	580 450	580 233	217
15	$3p^3{(}_3^{4}S){}^4{S}_{3/2}^{\circ }$	0.92 + 0.03 $3p^3{(}_3^{2}D){}^2{D}^{\circ }$ + 0.02 $3p^3{(}_1^{2}P){}^2{P}^{\circ }$	589 023	589 002	21
16	$3s3p{}^{3}P3d{}^4{F}_{3/2}^{\circ }$	0.96	641 955
17	$3p^3{(}_1^{2}P){}^2{P}_{1/2}^{\circ }$	0.80 + 0.13 $3s3p{}^{3}P3d{}^2{P}^{\circ }$ + 0.05 $3s3p{}^{1}P3d{}^2{P}^{\circ }$	642 591	642 310	281
18	$3s3p{}^{3}P3d{}^4{F}_{5/2}^{\circ }$	0.98	646 042	645 988	54
19	$3p^3{(}_1^{2}P){}^2{P}_{3/2}^{\circ }$	0.71 + 0.13 $3s3p{}^{3}P3d{}^2{P}^{\circ }$ + 0.04 $3s3p{}^{1}P3d{}^2{P}^{\circ }$	646 119	645 409	710
20	$3s3p{}^{3}P3d{}^4{F}_{7/2}^{\circ }$	0.98	651 972	651 946	26
21	$3s3p{}^{3}P3d{}^4{F}_{9/2}^{\circ }$	1.00	660 304	660 263	41
22	$3s3p{}^{3}P3d{}^4{P}_{5/2}^{\circ }$	0.65 + 0.28 $3s3p{}^{3}P3d{}^4{D}^{\circ }$ + 0.02 $3s3p{}^{3}P3d{}^2{D}^{\circ }$	690 311	690 304	7
23	$3s3p{}^{3}P3d{}^4{D}_{3/2}^{\circ }$	0.60 + 0.38 $3s3p{}^{3}P3d{}^4{P}^{\circ }$	692 653	692 662	−9
24	$3s3p{}^{3}P3d{}^4{D}_{1/2}^{\circ }$	0.87 + 0.12 $3s3p{}^{3}P3d{}^4{P}^{\circ }$	694 140	694 168	−28
25	$3s3p{}^{3}P3d{}^4{D}_{7/2}^{\circ }$	0.98	703 341	703 393	−52
26	$3s3p{}^{3}P3d{}^4{P}_{1/2}^{\circ }$	0.87 + 0.12 $3s3p{}^{3}P3d{}^4{D}^{\circ }$	703 826	703 750	76
27	$3s3p{}^{3}P3d{}^4{D}_{5/2}^{\circ }$	0.70 + 0.27 $3s3p{}^{3}P3d{}^4{P}^{\circ }$	704 114	704 114	0
28	$3s3p{}^{3}P3d{}^4{P}_{3/2}^{\circ }$	0.60 + 0.39 $3s3p{}^{3}P3d{}^4{D}^{\circ }$	704 202	704 209	−7
29	$3s3p{}^{3}P3d{}^2{D}_{3/2}^{\circ }$	0.51 + 0.27 $3s3p{}^{1}P3d{}^2{D}^{\circ }$ + 0.17 $3p^3{(}_3^{2}D){}^2{D}^{\circ }$	717 296	717 195	101
30	$3s3p{}^{3}P3d{}^2{D}_{5/2}^{\circ }$	0.48 + 0.25 $3s3p{}^{1}P3d{}^2{D}^{\circ }$ + 0.16 $3p^3{(}_3^{2}D){}^2{D}^{\circ }$	717 937	717 861	76
31	$3s3p{}^{3}P3d{}^2{F}_{5/2}^{\circ }$	0.65 + 0.32 $3s3p{}^{1}P3d{}^2{F}^{\circ }$	745 214	744 965	249
32	$3s3p{}^{3}P3d{}^2{F}_{7/2}^{\circ }$	0.65 + 0.33 $3s3p{}^{1}P3d{}^2{F}^{\circ }$	760 089	759 814	275
33	$3s3p{}^{3}P3d{}^2{P}_{3/2}^{\circ }$	0.77 + 0.14 $3p^3{(}_1^{2}P){}^2{P}^{\circ }$ + 0.06 $3s3p{}^{1}P3d{}^2{D}^{\circ }$	807 347	807 113	234
34	$3s3p{}^{3}P3d{}^2{P}_{1/2}^{\circ }$	0.86 + 0.12 $3p^3{(}_1^{2}P){}^2{P}^{\circ }$	815 394
35	$3s3p{}^{1}P3d{}^2{F}_{7/2}^{\circ }$	0.64 + 0.33 $3s3p{}^{3}P3d{}^2{F}^{\circ }$	817 790	817 593	197
36	$3s3p{}^{1}P3d{}^2{F}_{5/2}^{\circ }$	0.64 + 0.32 $3s3p{}^{3}P3d{}^2{F}^{\circ }$	820 795	820 601	194
37	$3s3p{}^{1}P3d{}^2{P}_{1/2}^{\circ }$	0.90 + 0.04 $3p^3{(}_1^{2}P){}^2{P}^{\circ }$	839 715	839 492	223
38	$3s3p{}^{1}P3d{}^2{P}_{3/2}^{\circ }$${}_a$	0.38 + 0.37 $3s3p{}^{1}P3d{}^2{D}^{\circ }$ + 0.08 $3s3p{}^{3}P3d{}^2{D}^{\circ }$	840 967	840 775	192
39	$3s3p{}^{1}P3d{}^2{P}_{3/2}^{\circ }$${}_b$	0.52 + 0.26 $3s3p{}^{1}P3d{}^2{D}^{\circ }$ + 0.08 $3s3p{}^{3}P3d{}^2{D}^{\circ }$	843 862	843 656	206
40	$3s3p{}^{1}P3d{}^2{D}_{5/2}^{\circ }$	0.67 + 0.17 $3s3p{}^{3}P3d{}^2{D}^{\circ }$ + 0.12 $3p^3{(}_3^{2}D){}^2{D}^{\circ }$	844 618	844 477	141

together with observed energies from the NIST database [168]. As can be seen from the table, the uncertainty of the calculated energies is at the order of 0.02%. Note that levels 38 and 39 have the same dominant terms in the wave function expansions. To give a single unique label, one has to revert to the recursive method of Froese Fischer and Gaigalas [123]. Whereas many of the lower levels have been identified, only a few levels with $n=4$ are known.

Table 8. Energies (in cm${}^{-1}$) for a few selected $n=4$ levels in Fe XIV. $E_{CI}$ are energies from [167], $E_{NIST}$ from NIST [168] and $E_{obs}$ from [169], along with three new tentative (T) values. Values in parentheses are percentage differences with respect to $E_{CI}$.

Pos	Configuration $\mathit{LS}$.	J	$\mathit{\pi }$	${\mathit{E}}_{\mathit{CI}}$	${\mathit{E}}_{\mathit{NIST}}$(%)	${\mathit{E}}_{\mathit{obs}}$(%)
101	$3s^{2}4s{}^2{S}_{}$	1/2	+	1 427 550	1 435 020 (0.5206)	1 426 965 (0.04) T
125	$3s^{2}4p{}^2{P}_{}$	1/2	−	1 541 937	1 568 840 (1.715)	1 541 394 (0.03)
128	$3s^{2}4p{}^2{P}_{}$	3/2	−	1 548 618	1 574 010 (1.613)	1 548 258 (0.02)
136	$3s3p{}^{3}P4s{}^2{P}_{}$	1/2	−	1 690 299	-	1 689 695 (0.0004)
150	$3s3p{}^{3}P4p{}^4{P}_{}$	3/2	+	1 795 164	-	1 795 032 (0.007) T
152	$3s3p{}^{3}P4p{}^4{P}_{}$	5/2	+	1 802 686	-	1 802 292 (0.02) T
181	$3s3p{}^{3}P4d{}^4{D}_{}$	5/2	−	1 930 871	-	1 933 758 (0.15)
184	$3s3p{}^{3}P4d{}^4{D}_{}$	7/2	−	1 935 340	-	1 938 452 (0.16)

compares calculated energies for some states with $n=4$ with available data from NIST and assignments by Del Zanna [169]. As can be seen from Table 8, certain energies in the NIST database may be incorrect. The agreement with Del Zanna is however excellent, confirming also the tentative assignments.

The calculations by Ekman et al. also provide rates for all E1, M1, and E2 transitions. In

Table 9. Multipoles (MT), transition energies $\mathsf{\Delta }E$, wavelengths $\lambda$ and transition probabilities A in s${}^{-1}$ in Fe XIV for selected transitions. $A_{CI}$ Ekman et al. [167], $A_{MCHF}$ Froese Fischer and Tachiev [170], $A_{MR-MP}$ Santana et al. [171], and $A_{NIST}$ NIST Atomic Spectra Database (2013) [168]. Accuracy estimates $dT$ have been computed based on transition probabilities in Babushkin and Coulomb gauges. The numbers in square brackets are powers of 10.

Upper Level	Lower Level	MT	$\mathbf{\Delta }\mathit{E}$ (cm${}^{-1}$)	$\mathit{\lambda }$ (Å)	${\mathit{A}}_{\mathit{CI}}$	$\mathit{dT}$	${\mathit{A}}_{\mathit{MCHF}}$	${\mathit{A}}_{\mathit{MR}-\mathit{MP}}$	${\mathit{A}}_{\mathit{NIST}}$
$3s^{2}3p{}^2{P}_{3/2}^o$	$3s^{2}3p{}^2{P}_{1/2}^o$	M1	18 854	5303.740	6.019[1]			6.016[1]
		E2	18 854	5303.740	1.474[−2]	0.004		1.466[−2]
$3s3p^2{(}_2^{3}P){}^4{P}_{1/2}$	$3s^{2}3p{}^2{P}_{1/2}^o$	E1	225 086	444.274	2.657[7]	0.062	2.620[7]	2.230[7]
	$3s^{2}3p{}^2{P}_{3/2}^o$	E1	206 231	484.891	9.777[6]	0.059	1.013[7]	8.693[6]
		M2	206 231	484.891	2.801[−1]			2.892[−1]
$3s3p^2{(}_2^{3}P){}^4{P}_{3/2}$	$3s^{2}3p{}^2{P}_{1/2}^o$	E1	232 776	429.596	5.851[5]	0.024	5.187[5]	4.833[5]
		M2	232 776	429.596	2.169[0]			2.193[0]
	$3s^{2}3p{}^2{P}_{3/2}^o$	E1	213 922	467.459	6.323[6]	0.083	5.908[6]	5.458[6]
		M2	213 922	467.459	6.644[−2]			6.838[−2]
	$3s3p^2{(}_2^{3}P){}^4{P}_{1/2}$	M1	7 690	13 002.904	1.007[1]			1.014[1]
		E2	7 690	13 002.904	2.050[−5]	0.001		1.003[−5]
$3s3p^2{(}_2^{3}P){}^4{P}_{5/2}$	$3s^{2}3p{}^2{P}_{3/2}^o$	E1	223 517	447.392	2.714[7]	0.061	2.491[7]	2.256[7]	2.5[7] ${}^c$
		M2	223 517	447.392	1.633[0]			1.685[0]
	$3s^{2}3p{}^2{P}_{1/2}^o$	M2	242 372	412.589	1.433[0]			1.423[0]
	$3s3p^2{(}_2^{3}P){}^4{P}_{1/2}$	E2	17 285	5 785.065	6.979[−3]	0.000		3.274[−3]
	$3s3p^2{(}_2^{3}P){}^4{P}_{3/2}$	M1	9 595	10 421.769	1.411[1]			1.411[1]
		E2	9 595	10 421.769	4.907[−4]	0.000		4.976[−4]
$3s3p^2{(}_2^{1}D){}^2{D}_{3/2}$	$3s^{2}3p{}^2{P}_{1/2}^o$	E1	299 401	333.999	2.426[9]	0.019	2.460[9]		2.3[9] ${}^b$
	$3s^{2}3p{}^2{P}_{3/2}^o$	E1	280 547	356.446	7.560[7]	0.003	8.669[7]		7.5[7] ${}^c$
$3s3p^2{(}_2^{1}D){}^2{D}_{5/2}$	$3s^{2}3p{}^2{P}_{3/2}^o$	E1	282 772	353.642	1.954[9]	0.027	1.998[9]		1.9[9] ${}^b$
$3s3p^2{(}_0^{1}S){}^2{S}_{1/2}$	$3s^{2}3p{}^2{P}_{1/2}^o$	E1	364 944	274.014	1.782[10]	0.011	1.716[10]		1.8[10] ${}^b$
	$3s^{2}3p{}^2{P}_{3/2}^o$	E1	346 090	288.942	1.082[9]	0.002	1.631[9]		1.2[9] ${}^c$
$3s3p^2{(}_2^{3}P){}^2{P}_{1/2}$	$3s^{2}3p{}^2{P}_{1/2}^o$	E1	388 711	257.261	1.279[10]	0.012	1.511[10]		1.4[10] ${}^b$
	$3s^{2}3p{}^2{P}_{3/2}^o$	E1	369 856	270.375	2.090[10]	0.012	2.144[10]		2.1[10] ${}^b$
$3s3p^2{(}_2^{3}P){}^2{P}_{3/2}$	$3s^{2}3p{}^2{P}_{1/2}^o$	E1	396 687	252.088	7.427[9]	0.009	7.902[9]		7.6[9] ${}^a$
	$3s^{2}3p{}^2{P}_{3/2}^o$	E1	377 832	264.667	3.254[10]	0.012	3.429[10]		3.38[10] ${}^a$

${}^a$ estimated uncertainty ≤25%, ${}^b$ estimated uncertainty ≤50%, ${}^c$ estimated uncertainty >50%.

, some of their rates are compared with results from other calculations and values from NIST. There is in general a good agreement between values from the different calculations. The quantity $dT$ defined by Equation (108) indicates that the uncertainty of the calculations by Ekman et al. is around a few percent, improving the available data from NIST.

5.4. Impact of External Magnetic Fields on Hyperfine Spectra: The 4205 Å Line in Eu II

Given a set of isotopic data (the nuclear spin, magnetic dipole and nuclear quadrupole moments) and a specified magnetic field strength, the hfszeeman95 program [67] allows the computation of complete hyperfine-Zeeman interaction matrices in a representation of the unperturbed ASFs, which were pre-calculated with grasp. From these interaction matrices, new perturbed eigenstates, labeled by $M_J$, F, or $M_F$ depending on the specified magnetic field and nuclear parameters (see Equation (109)) can be obtained through standard perturbation theory or from smaller configuration–interaction calculations, which in turns allows for determining radiative properties between these new eigenstates [148]. The perturbed eigenstates and associated radiative properties are conveniently computed with the mithit tool associated with the hfszeeman95 program.

As an example of results from such a computation, Figure 2 shows how the 4205 Å  fine-structure transition in ${}^{151,153}$Eu${}^{+}$ splits up into an complex pattern of lines when considering both isotopes together, including the symmetry breaking perturbations due to the hyperfine interaction with the nuclear moments together with the interaction with an external magnetic field set to 6 kG. In this computation, the atomic structure of the relevant states of Eu${}^{+}$ was first determined through standard MCDHF and RCI calculations with grasp, with special attention to spin-polarization in the electronic core represented via single substitutions from the inner s-shells. The hfszeeman95 calculation takes all off-diagonal effects into account, most notably, in this case, the mixing between different hyperfine states due to the external magnetic field. It is striking how complex a single fine-structure transition may become if these effects are taken into consideration.

6. Summary and Conclusions

Computers and atomic structure calculations have advanced to a level where large scale calculations in many cases produce transition energies that directly aid line identification in observed experimental or astrophysical spectra. Further methodological progress is however needed to meet the challenges from the complex Lanthanides, where the number of CSFs increases extremely rapidly with the increasing active set of orbitals [113]. The cancellation effects for the mass shift parameters in heavy and super-heavy neutral or near-neutral systems also pose problems. There is however good hope that these problems can be overcome by reduction techniques, grouping CSFs together to account for correlation in the atomic core without increasing the dimensions of the Hamiltonian matrix beyond the level that can be handled on large systems [173]. Furthermore, we can expect to see methods utilizing the fact that angular integration is independent of the principal quantum numbers, potentially decreasing the computational time for CI calculations based on large active sets [152].

While computational methodologies and codes improve, so does documentation and instructions for efficient use. The computer codes now come with detailed manuals describing different modes of use, guiding new groups of users to produce reliable atomic data that meet the set-up requirements [43]. We would like to conclude the present paper that provides the theoretical support of the grasp2018 manual [72], by mentioning some present and forthcoming theoretical, methodological, algorithmic and code developments that illustrate the dynamics of the GRASP community:

(i)

Investigating new algorithms and methods to test and improve the numerical accuracy and stability of the variational method [139,174];

(ii)

Designing efficient CSF generators that drastically reduce the computational load of MCDHF and RCI calculations [152];

(iii)

Utilizing fast biorthogonal transformation techniques to handle non-orthogonalities and allow the CSF spaces to be built from several separately optimized and mutually non-orthogonal orbital sets [152,173,175];

(iv)

Extending the grasp2018 hyperfine codes [66] to treat the magnetic octupole [176] and electric hexadecapole hyperfine interactions, as well as the Bohr–Weisskopf effect [177], which allowed to resolve the nuclear-octupole-moment puzzle in ${}^{173}$Yb [178,179];

(v)

Lifting current restrictions on maximum occupation numbers of two for orbitals $j=9/2$ in the CSF list generation to fully exploit the available spin-angular library;

(vi)

Searching for original methods that open promising perspectives for performing rigorous QED calculations within the grasp framework [180];

(vii)

Implementing the model-QED-operator approach of Shabaev et al. [181,182,183,184] in grasp and testing the results through systematic comparisons with other QED approaches and observations [58,59,60,185];

(viii)

Coding and implementing other QED approaches [186,187,188,189] in the grasp suite of programs [190].

Finally, we would like to acknowledge fruitful collaboration [191,192,193,194,195,196,197] with the authors and users of the sibling mdfgme code [36,198], which is based on the same relativistic multiconfiguration theory, but using different approaches for solving the variational equations and for evaluation of the spin-angular factors.