Independently Optimized Orbital Sets in GRASP—The Case of Hyperfine Structure in Li I

Abstract

In multiconfiguration Dirac–Hartree–Fock (MCDHF) calculations, there is a strong coupling between the localization of the orbital set and the configuration state function (CSF) expansion used to determine it. Furthermore, it is well known that an orbital set resulting from calculations, including CSFs describing core–core correlation and other effects, which aims to lower the weighted energies of a number of targeted states as much as possible, may be inadequate for building CSFs that account for correlation effects that are energetically unimportant but decisive for computed properties, e.g., hyperfine structures or transition rates. This inadequacy can be traced in irregular or oscillating convergence patterns of the computed properties as functions of the increasing orbital set. In order to alleviate the above problems, we propose a procedure in which the orbital set is obtained by merging several separately optimized, and mutually non-orthogonal, orbital sets. This computational strategy preserves the advantages of capturing electron correlation on the total energy through the variational MCDHF method and allows to target efficiently the correlation effects on the considered property. The orbital sets that are merged are successively orthogonalized against each other to retain orthonormality. The merged orbital set is used to build CSFs that efficiently lower the energy and also adequately account for the correlation effects that are important for the property. We apply the procedure to compute the hyperfine structure constants for the $1s^{2}2s{}^2{S}_{1/2}$ and $1s^{2}2p{}^2{P}_{1/2,3/2}^o$ states in ${}^7$Li and show that it leads to considerably improved convergence patterns with respect to the increasing orbital set compared to standard calculations based on a single orbital set, energy-optimized in the variational procedure. The perspectives of the new procedure are discussed in a broader context in the summary.

1. Introduction

The fully relativistic multiconfiguration Dirac–Hartree–Fock (MCDHF) method, as implemented in the General Relativistic Atomic Structure Package (grasp) [1], has for many years been used to provide a broad user community with crucial atomic data. In astrophysics, the computed transition energies and transition data for many elements and various ionization stages have been used for abundance analysis and plasma diagnostics. In some cases, the computed transition energies are of spectroscopic accuracy, which renders them useful for the unambiguous identification of lines in recorded spectra [2,3,4,5]. Computed broadening parameters, such as hyperfine structures and isotope shifts, have been used for detailed line modeling, revealing isotopic ratios and giving clues to nuclear synthesis [6]. Further, Landé $g_J$-factors and magnetically induced transitions have been used to infer knowledge about the magnetic fields of the Sun’s atmosphere [7]. Calculated hyperfine interaction constants, field shifts, and mass isotope shifts [8,9], related to detailed electron–nucleus interactions [10], have been combined with high-precision measurements to extract nuclear information, such as magnetic dipole, electric quadrupole [11] and magnetic octupole [12] moments, changes in nuclear radii along isotopic chains [13], and nuclear deformations [14,15].

Multiconfiguration methods are conceptually simple: expand the wave functions for atomic fine-structure states in a limited basis of configuration state functions (CSFs) constructed from antisymmetrized and $jj$-coupled products of relativistic orbitals and apply the variational principle on a weighted energy average of the states to derive equations for the radial parts of the relativistic orbitals. Solve the equations iteratively in the so-called self-consistent field (SCF) procedure to obtain a radial orbital basis. In the subsequent configuration interaction (CI) calculations, expand the wave functions in a larger basis of CSFs and determine the expansion coefficients by diagonalizing the corresponding Hamiltonian matrix. Given the approximate wave functions, measurable properties, such as hyperfine structures, isotope shifts, and transition rates, are evaluated using the first-order perturbation theory as the expectation values or amplitudes of operators describing the properties [16]. Guided by rules based on the Z-dependent perturbation theory, the CSF bases and the orbitals used for their constructions are systematically increased to account for valence–valence, core–valence, core–core electron correlation, as well as spin and orbital polarizations and radial correlation of importance for both the total energy and the computed properties [16,17].

The spatial locations of the orbitals depend on the weighted energy average of the CSFs used in the SCF procedure; the challenge is to determine an orbital basis that saturates all of the above correlation effects, some of which are energetically unimportant but more important for a computed property. Failure to determine an orbital basis that is adequate for describing important correlation effects for the property at hand may lead to slow and irregular convergence as the orbital basis is enlarged or even convergence to an incorrect limit [10,18,19].

In this work, we analyze the situation and propose an efficient and computationally cheap procedure to deal with the problems above. The new computational methodology can be used to generate an orbital basis that is suitable for describing the correlation effects that are important for a property and, at the same time, suitable to account for the energetically important correlation effects. We use the procedure to compute the hyperfine constants for the $1s^{2}2s{}^2{S}_{1/2}$ and $1s^{2}2p{}^2{P}_{1/2,3/2}^o$ states in ${}^7$Li and to demonstrate better convergence properties with respect to the increasing orbital basis, compared with standard calculations.

2. Variational Calculations

2.1. The MCDHF Method

In the multiconfiguration Dirac–Hartree–Fock (MCDHF) method, the atomic state function (ASF) for an atomic state $\mathsf{\Gamma }J{M}_J\pi$ is a linear combination of configuration state functions (CSFs)

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

(1)

where ${\gamma }_{\alpha }$ specifies the orbital occupancies and spin-angular coupling tree quantum numbers of each CSF. The CSFs are antisymmetrized and symmetry-adapted many-electron functions, built from relativistic single electron Dirac central field orbitals,

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

(2)

where $P_{n\kappa }\left(r\right)$ and $Q_{n\kappa }\left(r\right)$ are the large and small components of the radial wave function, and $\mathsf{\Omega }(\theta ,\phi )$ is the angular wave function for each component. $P_{n\kappa }\left(r\right)$ and $Q_{n\kappa }\left(r\right)$ are represented on a logarithmic grid and are required to be orthonormal within each symmetry.

The MCDHF method is energy-driven and may simultaneously target several states ${\mathsf{\Gamma }}^i{J}^i,i=1,\dots ,N_{ASF}$. The radial orthonormality condition is introduced as

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

(3)

for the orbitals with the same $\kappa$-value (${\kappa }_a={\kappa }_b=\kappa$), and applying the variational principle to the statistically weighted energy functional of the targeted states

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

(4)

where the Lagrange multipliers ${\lambda }_{ab}$ are introduced to ensure the orthonormality of the orbitals. We obtain a set of equations (see [20,21]) of the form

$$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].$$

(5)

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

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

(6)

where $V_{nuc}\left(r\right)$ is the nuclear potential. $Y(a;r)$, the direct potential, and $\overline{X}(a;r)$, the exchange potential, result from the variations of integrals associated with the Dirac–Coulomb Hamiltonian matrix elements $\langle {\mathsf{\Gamma }}^i{J}^i\left|\right|{\mathcal{H}}_{DC}\left|\right|{\mathsf{\Gamma }}^i{J}^i\rangle$ between the ASFs. Here, $w_a$ is the generalized occupation number of orbital a, and ${ϵ}_{ab}$ are the energy parameters related to the Lagrange multipliers. In the MCDHF program, as implemented in grasp [1], the orbital variational equations, together with the accompanying boundary conditions for the radial functions, are solved numerically using finite difference methods through an iterative process, until the self-consistent field (SCF) is obtained. This SCF procedure is coupled with the eigenvalue problem that provides the CSF mixing coefficients $\left\{c_{\alpha }^{\mathsf{\Gamma }J}\right\}$ of the $N_{ASF}$ atomic states considered—see [16,20,22] for details.

2.2. Localization of the Radial Orbitals and Their Dependence on the Energy Functional

The radial functions resulting from solving the variational equations strongly depend on the energy functional and thus the CSF expansion. This fact can be used to our advantage to capture some specific correlation effects that are important for a property but can also bring undesirable distortions of the wave function for the property of interest (see [23] for an in-depth discussion). As an example, we consider the ground state $1s^{2}2s^2{}^1{S}_0$ in Be I. To account for valence–valence (VV) correlation, we start with a set of radial orbitals (to be determined) $\{1s,2s,2p\text{-},2p,3s,3p\text{-},3p,3d\text{-},3d,\dots ,\}$ and allow double (D) substitutions from the outer ($2s$) orbital of the $1s^{2}2s^2$ reference configuration to the set of correlation orbitals, which eventually form an expansion

$${\mathsf{\Psi }}^{VV}\left(1s^{2}2s^2{}^1{S}_0\right)=c_0\mathsf{\Phi }\left(1s^{2}2s^2{}^1{S}_0\right)+\sum _{n\kappa ,n^{\prime }{\kappa }^{\prime }}{c}_{n\kappa ,n^{\prime }{\kappa }^{\prime }}\mathsf{\Phi }\left(1s^{2}n\kappa n^{\prime }{\kappa }^{\prime }{}^1{S}_0\right).$$

(7)

To account for core–valence (CV) correlation, we allow D substitutions from the reference configuration in which one substitution is from the outer orbital, and one is from the inner orbital to a set of correlation orbitals to yield

$${\mathsf{\Psi }}^{CV}\left(1s^{2}2s^2{}^1{S}_0\right)=c_0\mathsf{\Phi }\left(1s^{2}2s^2{}^1{S}_0\right)+\sum _{n\kappa ,n^{\prime }{\kappa }^{\prime }}{c}_{n\kappa ,n^{\prime }{\kappa }^{\prime }}\mathsf{\Phi }\left(1s2sn\kappa n^{\prime }{\kappa }^{\prime }{}^1{S}_0\right).$$

(8)

Finally, to account for core–core (CC) correlation, we allow D substitutions from the inner orbital of the reference configuration to the set of radial orbitals which yields

$${\mathsf{\Psi }}^{CC}\left(1s^{2}2s^2{}^1{S}_0\right)=c_0\mathsf{\Phi }\left(1s^{2}2s^2{}^1{S}_0\right)+\sum _{n\kappa ,n^{\prime }{\kappa }^{\prime }}{c}_{n\kappa ,n^{\prime }{\kappa }^{\prime }}\mathsf{\Phi }\left(2s^{2}n\kappa n^{\prime }{\kappa }^{\prime }{}^1{S}_0\right).$$

(9)

The radial orbitals with principal quantum numbers up to $n=4$, resulting from independent MCDHF calculations based on the above expansions, are shown in Figure 1.

Not surprisingly, the correlation orbitals used to build CSFs that account for valence–valence electron correlation strongly overlap the outer $2s$ orbital. The orbitals generated to build CSFs that account for core–core electron correlation strongly overlap the inner $1s$ orbital. The orbitals used to build CSFs that account for core–valence electron correlation effects overlay the region in space in which the probability distributions of both $1s$ and $2s$ overlap. Although we only show a subset of correlation orbitals, it is clear from Figure 1 that orbitals from a core–core calculation are inefficient at accounting for valence–valence correlation and vice versa.

3. Computed Properties and Their Dependence on Correlation Effects

In this section, we discuss computed properties and their dependencies on correlation effects. For simplicity, we restrict ourselves to a hyperfine structure with the $1s^{2}2s{}^2{S}_{1/2}$ and $1s^{2}2p{}^2{P}_{1/2,3/2}^o$ states in ${}^7$Li as specific examples.

3.1. Hyperfine Structure

The hyperfine structure of an atomic energy level is caused by the interaction between the electrons and the electromagnetic multipole moments of the nucleus. The Hamiltonian for the interaction may be represented as a multipole expansion [24]

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

(10)

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 [25]. The k = 1 term represents the magnetic dipole interaction and the k = 2 term represents the electric quadrupole interaction. Higher-order terms are much smaller and can often be neglected. For an N-electron atom, the electronic tensor operators are, in atomic units,

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

(11)

$${\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).$$

(12)

In the formulas above, $\alpha$ is the fine-structure constant, $\mathit{\alpha }$ is the Dirac matrix, and ${\mathbf{C}}^{\left(k\right)}$ is the renormalized spherical harmonic, a spherical tensor operator of rank k. The splitting of atomic fine-structure levels due to the hyperfine interaction is often given in terms of the interaction constants A and B

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

(13)

$$B_{\mathsf{\Gamma }J}=2Q_I{\left[\frac{J(2J-1)}{(J+1)(2J+1)(2J+3)}\right]}^{1/2}\langle \mathsf{\Gamma }J\Vert {\mathbf{T}}^{\left(2\right)}\Vert \mathsf{\Gamma }J\rangle ,$$

(14)

where I is the nuclear spin, ${\mu }_I$ is the nuclear magnetic dipole moment, and $Q_I$ is the nuclear electric quadrupole moment.

3.2. Polarization Effects

The most important correlation effect for the hyperfine structure is the polarization of the atomic core due to the Coulomb interaction with the valence electrons. Of particular importance is the spin polarization of the closed s shells. If the two s-electrons in the same subshell have different spin densities at the nucleus, a contact interaction is induced in the non-relativistic limit. Since inner s-electrons have high, but canceling, spin densities at the nucleus, a very small imbalance due to the interactions with the outer electrons is sufficient to cause a contribution to the A constant comparable to that of an open valence shell. A schematic illustration is shown in Figure 2. The Coulomb exchange interaction reduces the repulsion between core and valence electrons with the same spin orientation, pulling the core electron towards the valence subshell and causing an imbalance.

For example, consider the $1s^{2}2s{}^2{S}_{1/2}$ and $1s^{2}2p{}^2{P}_{1/2,3/2}^o$ states in ${}^7$Li. In the multiconfiguration method, adhering to the non-relativistic notation, spin polarization (SP) can be described by adding a CSFs of the form $\mathsf{\Phi }(1s1s_2{(}^{3}S)2s{}^{2}S)$ and $\mathsf{\Phi }(1s1s_2{(}^{3}S)2p{}^2{P}^o)$, respectively, to yield

$${\mathsf{\Psi }}^{SP}\left(1s^{2}2s{}^{2}S\right)=c_1\mathsf{\Phi }\left(1s^{2}2s{}^{2}S\right)+c_2\mathsf{\Phi }(1s1s_2{(}^{3}S)2s{}^{2}S)$$

(15)

and

$${\mathsf{\Psi }}^{SP}\left(1s^{2}2p{}^2{P}^o\right)=c_1\mathsf{\Phi }\left(1s^{2}2p{}^2{P}^o\right)+c_2\mathsf{\Phi }(1s1s_2{(}^{3}S)2p{}^2{P}^o),$$

(16)

where $1s_2$ is allowed to be non-orthogonal to $1s$ and $2s$. For $1s^{2}2p{}^2{P}^o$, there is also an orbital polarization (OP), which is due to the interaction of the core with the non-spherical charge distribution of the valence orbital. This interaction distorts the spherical symmetry of the core, leading to additional contributions, in the non-relativistic limit, to the orbital, spin-dipolar, and quadrupole terms. In the multiconfiguration method, orbital polarization can be described by adding CSFs of the forms $\mathsf{\Phi }(1s3d_3{(}^{3}D)2p{}^2{P}^o)$ and $\mathsf{\Phi }(1s3d_4{(}^{1}D)2p{}^2{P}^o)$ to finally yield

$$\begin{array}{ccc}\hfill {\mathsf{\Psi }}^{SP+OP}\left(1s^{2}2p{}^2{P}^o\right)& =& c_1\mathsf{\Phi }\left(1s^{2}2p{}^2{P}^o\right)+c_2\mathsf{\Phi }(1s1s_2{(}^{3}S)2p{}^2{P}^o)\hfill \\ \hfill & & +c_3\mathsf{\Phi }(1s3d_3{(}^{3}D)2p{}^2{P}^o)+c_4\mathsf{\Phi }(1s3d_4{(}^{1}D)2p{}^2{P}^o),\hfill \end{array}$$

(17)

where $3d_3$ and $3d_4$ are non-orthogonal [17].

3.3. Localization of the Polarization Orbitals

The efficiency of the above representation relies on the assumption that the correlation orbitals are optimally localized to describe the polarization effects [18]. In relativistic MCDHF calculations, which require the orbitals to be orthonormal in the current version of grasp, such an orbital basis, optimally localized for describing SP and OP, can be obtained as follows:

• Perform a weighted average Dirac–Fock calculation for $1s^{2}2s{}^2{S}_{1/2}$ and $1s^{2}2p{}^2{P}_{1/2,3/2}^o$;
• Keep $1s,2s,2p\text{-},2p$ frozen and perform weighted average MCDHF calculations for $1s^{2}2s{}^2{S}_{1/2}$ and $1s^{2}2p{}^2{P}_{1/2,3/2}^o$ based on the CSF expansions formed by allowing S’s substitution from the reference configuration to a set of $s,d\text{-},d$ orbitals;
• Compute the hyperfine interaction constants and monitor the convergence as the set of $s,d\text{-},d$ orbitals is increased;
• Stop when the hyperfine interaction constants are not changing anymore.

In

Table 1. Hyperfine interaction constants (in MHz) for the $1s^{2}2s{}^2{S}_{1/2}$ and $1s^{2}2p{}^2{P}_{1/2,3/2}^o$ states in ${}^7$Li from weighted average DF calculations and weighted average MCDHF calculations including spin and orbital polarization effects with four layers of $s,d\text{-},d$ orbitals, i.e., orbitals set up to $\{6s,6d\text{-},6d\}$. The nuclear parameters for ${}^7$Li are $I=3/2$, ${\mu }_I=3.256424{\mu }_N$, and $Q=-0.0406$ barn.

	$\mathit{A}({}^2{\mathit{S}}_{1/2})$	$\mathit{A}({}^2{\mathit{P}}_{1/2}^{\mathit{o}})$	$\mathit{A}({}^2{\mathit{P}}_{3/2}^{\mathit{o}})$	$\mathit{B}({}^2{\mathit{P}}_{3/2}^{\mathit{o}})$
DF	290.249	32.356	−6.469	$-0.2235$
$\{3s,3d\text{-},3d\}$	374.047	44.741	$-4.910$	$-0.2238$
$\{4s,4d\text{-},4d\}$	380.692	42.461	$-2.619$	$-0.1968$
$\{5s,5d\text{-},5d\}$	380.341	42.610	$-2.766$	$-0.1968$
$\{6s,6d\text{-},6d\}$	380.342	42.611	$-2.766$	$-0.1968$
Exp.	401.752043 ${}^a$	45.914(25) ${}^b$	−3.055(14) ${}^b$

^a Ref. [27], ^b Ref. [28].

, we report the computed hyperfine interaction constants in MHz as functions of the increasing set of $s,d\text{-},d$ orbitals. As can be seen from the table, these small calculations (the calculations based on orbitals with a value of n up to 6 include only 96 CSFs distributed over the three states) already capture the bulk of the correlation effects. Most notable is the hyperfine constant for $1s^{2}2p{}^2{P}_{3/2}^o$, which changes sign when the polarization effects are added.

The radial orbitals resulting from the polarization calculations are shown on the left panel in Figure 3, together with the corresponding spectroscopic $1s,2s,2p\text{-},$ and $2p$ orbitals. The polarization orbitals are localized in a region between core electrons and valence electrons to account for the exchange force. For comparison, radial orbitals resulting from an energy-driven calculation, further described in Section 4.1, that include CSFs obtained from SDT substitutions are shown on the right panel in Figure 3. The correlation orbitals from the energy-driven calculation are more contracted, compared with the polarization orbitals, reflecting the fact that the bulk of the correlation energy comes from the interactions within the atomic core. We infer that the orbitals from the energy-driven calculation are unsuited to capture the important polarization effects.

4. Hyperfine Interaction Constants in Different Orbital Bases

4.1. Orbital Basis from Energy-Driven Calculations

We start with weighted average MCDHF calculations for the $1s^{2}2s{}^2{S}_{1/2}$ and $1s^{2}2p{}^2{P}_{1/2,3/2}^o$ states in ${}^7$Li, based on the expansions obtained by allowing single (S), double (D), and triple (T) substitutions from the {$1s^{2}2s,1s^{2}2p$} reference configurations to an active set of orbitals. The set is increased layer-by-layer (see [16] §3.5) up to {$13s,13p,13d,13f,13g$}. The number of CSFs in the final even and odd state expansions are, respectively, 92 015 and 255 892. These calculations are energy-driven. They include all CSFs that can be generated from the orbital basis, and the orbitals are localized in space by the variational principle, which promotes those correlation effects that maximally lower the weighted energy. In Figure 4, we display, in red squares, the resulting hyperfine interaction constants A and B, in MHz, as functions of the increasing orbital set. The constants were computed using $I=3/2$, ${\mu }_I=3.256424{\mu }_N$, and $Q=-0.0406$ barn. The behavior for the $1s^{2}2s{}^2{S}_{1/2}$ is irregular, whereas for other states, it is oscillating, making it difficult to establish an accurate final value (cf. Figure 2 in [10]).

The expansions include the CSFs that describe the important polarization effects (single substitutions from the atomic core). However, the orbitals building the CSFs are not spatially localized in the correct region, as seen in Figure 3. The correlation orbitals overlap mainly the $1s$ orbital with a tendency for each layer to be either a little more contracted or a little more diffuse each time to lower the weighted energy as much as possible with the resulting oscillations of the hyperfine constants. At the same time, we observe an overall increase in the absolute values of the magnetic dipole hyperfine constants A of the J = 1/2 states, compared with the respective absolute DF values (cf. Table 1). This is consistent with the fact that correlation effects break up the closed atomic core, inducing contributions similar to those from open shells.

It should be noted that the observed oscillations in the hyperfine constants do not result from the layer-by-layer methodology [16] used in this work. A fully variational calculation performed with the MCDFGME code [29,30] for Li-like Be${}^{+}$, Ne${}^{7+}$, and Cl${}^{14+}$ manifested similar oscillatory behavior, although, in the fully variational calculation, the oscillations are slightly more damped than in the layer-by-layer case. The MCDFGME code has also been used to provide fully relaxed, larger-scale calculations for ${}^7$Li, as compared to the calculations from Ref. [31]. The resulting hyperfine constants are displayed as green circles in Figure 4. The wavefunctions include all SDT substitutions to orbitals with n up to 5, i.e., $\{5s,5p,5d,5f,5g\}$, and are fully relaxed. The single substitutions include also those relevant to the Brillouin theorem [32]. This leads to 1463, 1454, and 2478 configurations for the ${}^2{S}_{1/2}$, ${}^2{P}_{1/2}^o$, and ${}^2{P}_{3/2}^o$ states, respectively. It should be pointed out that it is increasingly difficult to obtain convergence in fully variational calculations when increasing the number of orbital layers, and we were unable to go beyond $n=5$. We have included orbitals up to g in both cases and thus should have comparable convergence in the angular correlation contribution. However, full relaxation of all orbitals corresponds to more correlation than in the case whereby the lower orbital layers are kept frozen. Strong oscillations are indeed observed, where most of the differences with the layer-by-layer calculations can be attributed to the fact that the fully relaxed calculations are performed separately for each state in the Optimal Level (OL) approach, whereas the layer-by-layer calculations are performed in the Extended Optimal Level (EOL) approach [33,34].

Oscillatory behavior has also been observed in earlier calculations of hyperfine structures for Sn, Zn, At, Cd, Ra, Au, Hg, Sc, Br, I, Bi, Be${}^{+}$, Fe${}^{6+}$, Ti, Y, Sc${}^{+}$, and Y${}^{+}$. Thus, we can conclude that the observed oscillations are inherent to the variational method that prioritizes the energy functional over the targeted atomic property. The latter observation can be understood from the time-independent perturbation theory by looking at the order-by-order expansion of the targeted observable (other than the energy). The terms of this expansion involve matrix elements coupling the zeroeth-order unperturbed states through not only the Hamiltonian but also the relevant (here, hyperfine) operators (cf. Equation (30) in [35]). Many of them do not appear in the energy functional on which the variational principle is applied to derive the MCDHF equations.

4.2. Polarization Orbitals Augmented to the Orbital Basis from Energy-Driven Calculations

Polarization orbitals need to be included in energy-driven layer-by-layer calculations, in order to better describe the correlation effects that are important for the hyperfine constants. However, due to the orthogonality restrictions, separately optimized and optimally localized polarization orbitals can not be directly augmented to the orbital basis from energy-driven layer-by-layer calculations. To overcome the orthogonality constraints, we use the rwfnrelabel program in grasp to relabel (change the principal quantum number that serves as the orbital identifier) the polarization orbitals so that they appear after the orbitals from the energy-driven layer-by-layer calculation. The relabeled polarization orbitals are then, using the rwfnestimate program in grasp, orthogonalized against the orbitals from the energy-driven layer-by-layer calculation and finally augmented. Then, we perform CI calculations with the CSFs obtained by allowing SDT substitutions from the {$1s^{2}2s,1s^{2}2p$} reference configurations to orbitals optimized in the energy-driven calculation. To these CSFs, which mainly describe correlation effects that are important for lowering the energy and, to a lesser extent, the effects that are important for the hyperfine structure, we add the CSFs obtained by S substitutions from the atomic core to the augmented, orthogonalized polarization orbitals. For energy-optimized orbital bases with $n<10$, four layers of orthogonalized polarization orbitals were augmented. For energy-optimized orbital bases with $n\ge 10$, due to linear dependencies in the orbital basis, only two polarization layers were augmented. The number of CSFs in the final even and odd state expansions are, respectively, 92,541 and 258,432, i.e., only negligibly larger than those from the energy-driven calculation. In Figure 4, we display, in blue crosses, the resulting hyperfine interaction constants A and B in MHz, as functions of the increasing energy-driven orbital set (to which two or four layers of polarization orbitals have been added). Now, the oscillations are almost completely damped out, and the final values of the interaction constants can already be accurately established from a limited, $n=9$, energy-driven orbital set merged with the polarization orbitals. The CPU times of the $n=13$ CI calculations are 20m42s in energy-driven calculations and 25m5s for energy-driven orbitals merged with the polarization orbitals’ calculation for the $1s^{2}2s{}^2{S}_{1/2}$ state, using 64 nodes on a cluster with the AMD EPYC 7542 32-Core Processor. The bulk of the CPU time is spent on the calculations for energy-driven orbital sets. The interaction constants, on the other hand, can be accurately established at $n=10$ during the calculation of the energy-driven orbital set merged with the polarization orbitals, for which the CPU time of the CI calculation is 6m30s on the $1s^{2}2s{}^2{S}_{1/2}$ state. That means that the efficiency can be improved by using a smaller orbital set compared with energy-driven calculations. The use of separately optimized polarization orbitals merged with energy-driven orbital sets represents a great improvement in both efficiency (the time-consuming parts are the calculations for large energy-driven orbital sets) and accuracy. To provide the final values of the interaction constants, the effects from the neglected orbitals with a higher l, as well as from the Breit interaction, QED, and nuclear recoil, have to be added [31,36,37]. This is, however, outside the scope of the present paper.

5. Summary and Conclusions

Orbitals resulting from MCDHF calculations strongly depend on the employed CSF expansion. Expansions describing, e.g., core–core correlation and other effects aiming at lowering the weighted energy of the targeted states as much as possible may yield an orbital basis that is inadequate for building CSFs that account for correlation effects that are important for other properties. An example is provided by hyperfine interaction constants. These constants are affected by spin and orbital polarizations, described by the CSFs obtained by allowing single substitutions from the atomic core to an orbital basis localized in the region of space between the core electrons and the valence electron(s). This region is energetically unimportant and not well described by the orbitals from energy-driven calculations. The orbitals from energy-driven calculations tend to overlap the orbitals of closed shells; at times they are a little more contracted, and, at other times, they are a little more expanded—a behavior that may lead to oscillations for properties other than the energy. To capture both the energetically important effects as well as the important effects for the hyperfine structure, we propose a procedure in which the polarization orbitals that are optimally localized to describe the polarization effects are relabeled and orthogonalized against orbitals from energy-driven calculations. The procedure is computationally efficient and dramatically improves the convergence properties of the computed hyperfine interaction constants as functions of the increasing orbital set.

The proposed procedure is general and can be applied to compute accurate values of other properties weakly coupled to the energy, as well as energy differences. Applied to the transition parameters from high Rydberg states, it can be used to accurately describe the outer part of the wave function that is important for the length form of the transition operator but energetically unimportant [19]. Applied to energy differences, separate calculations may be performed, targeting core–core, core–valence, and valence–valence electron correlations. The two latter orbital sets are relabeled and orthogonalized against the first one and against each other and then augmented to the first one. In the final CI calculation, the CSFs that account for core–core correlation are built from the core–core orbital set, and the CSFs that account for core–valence correlation are built from the core–core orbital set to which the core–valence orbitals were augmented. Finally, the CSFs that account for valence–valence correlation are built from the full orbital set, including the core–core, core–valence, and valence–valence orbitals. This is similar to the method described in [38], but it uses a merged and orthogonalized orbital basis. It should, in this context, be emphasized that calculations of a specific atomic property require accurate descriptions of effects that are of direct importance for that property, as well as the effects that are important for the total energy, as can be seen from the crossed second-order terms in the perturbative expression for a property (see [35] §7).

Although simple and readily available with the current tools in grasp, rwfnrelabel, and rwfnestimate [22], the present procedure has the drawback that it requires CSFs that are important for the property, or, in the case of energy separations, core–valence, and valence–valence correlation, to be built on merged orbital sets, which implies a certain overhead. An alternate and more general solution circumventing the problems above would be to employ different and mutually non-orthogonal orbital sets and deal with the non-orthogonalities by bi-orthogonal orbital transformations, as in [35,38]. This work is in progress based on the new concept of configuration state function generators (CSFGs) [39].