Evaluation and Comparison of the Configuration Interaction Calculations for Complex Atoms

Configuration interaction (CI) methods are the method of choice for the determination of wave functions for complex atomic systems from which a variety of atomic properties may be computed. When applied to highly ionized atoms, where few, if any, energy levels from observed wavelengths are available, the question arises as to how a calculation may be evaluated. Many different codes are available for such calculations. Agreement between the results from different codes in itself is not a check on accuracy, but may be due to a similarity in the computational procedures. This paper reviews basic theory, which, when applied in a systematic manner, can be the basis for the evaluation of accuracy. Results will be illustrated in the study of 4s 2 4p 5 (odd) and 4s 2 4p 4 4d (even) levels in W 39+ and the transitions between them.


Introduction
Currently, tungsten (W) is an element of great interest as a possible wall material in future tokamaks, but progress is hindered by the lack of accurate spectroscopic data.In connection with a study of heavy ion impurities in fusion plasmas, Fournier [1] performed extensive collisional-radiative calculations (energy levels, radiative transition probabilities for allowed and forbidden transitions, collisional cross-sections) for candidate lines in observed spectra for W 37+ -W 47+ using the RELAC code of Klapisch et al. [2,3].This data was used by Utter et al. [4] to identify some of the lines observed in their electron beam ion trap (EBIT).Radtke et al. [5], using a similar HULLAC code [6] for their predictions and high-resolution X-ray and EUV spectroscopy for their measurements, extended the ions of interest to W 25+ -W 47+ .Similarly, Ralchenko et al. [7] observed spectra in the W 39+ -W 47+ ion range with an EBIT and confirmed the identification of some of the lines through the use of Flexible Atomic Code (FAC) [8].The tentative identification of lines has spurred others to improve upon these results using more recent variational codes that are fully relativistic and can include both Breit and quantum electrodynamic (QED) effects, as well as correlation in a many-electron system.Froese Fischer and Gaigalas [9] used GRASP2K [10] for the study of the 4p 6 4d, 4p 6 4f and 4p 5 4d 2 configurations in W 37+ and the transitions among these levels, Aggarwal, Jha and Mohan (AJM) [11] and also Aggarwal and Keenan [12] used a similar version called GRASP [13,14] for 4s 2 4p 5 , 4s 2 4p 4 4d and 4s4p 6 configurations in W 39+ , and Quinet [15] used the same code to produce a theoretical survey of 4p k and 4d k ground configurations in W 29+ -W 43+ from which, M1 transition data were computed.
In highly ionized atoms where few, if any, classified levels are available, the question of accuracy remains.Codes for atomic structure calculations do not, by themselves, determine the accuracy of the final result, but rather, the computational procedure that was adopted.A good example is the calculation reported by Jonauskas et al. [16], where GRASP2K was used to get the radial functions for orbitals, but the calculation of the configuration interaction (CI) strength for hundreds of configurations was performed using the method of global characteristics.This calculation is hardly a typical GRASP2K calculation, but the results from this calculation, referred to as GRASP2K results, were found to be in large disagreement with other methods [15].Codes have many options, and there can be significant disagreement, even when the same codes are used [11,12].
Configuration interaction methods are the method of choice for complex atomic systems [17].Many options determine the computational model.This paper describes factors that define the accuracy of the final result and strategies that can be used for assessing the accuracies.For lighter elements, energy levels are available for comparison with theory, but when these are not available, other methods must be considered.Though there are many similarities between non-relativistic methods corrected for relativistic effects (see [18] for the basic theory) and fully relativistic calculations to which Breit and QED corrections are added, this paper is concerned primarily with the latter, though, at times, the problem may be discussed in non-relativistic terms.

Underlying Theory
In a configuration interaction (CI) calculation, the wave function for an atomic state function (ASF), Ψ, is assumed to be a linear combination of configuration state functions (CSFs), Φ, namely: where each CSF is antisymmetric in the coordinates of the electrons and has radial and angular factors that together yield an eigenstate for total angular momentum and total spin quantum numbers LS if the calculations are performed starting with non-relativistic ls orbitals, or J and parity if orbitals are fully relativistic and coupled in the jj-coupling scheme.Thus, one configuration may generate several CSFs.In all CI methods, the vector of expansion coefficients is defined by an eigenvector of the interaction matrix, and the energy is the eigenvalue of the associated eigenvector.What differs among calculations are: (1) the configurations that define the CSFs in the expansion and (2) how the radial function of the one-electron orbitals that define the CSFs are determined.
The latter may be obtained from a model potential, a simple energy expression that depends only on the configuration in the extended average level (EAL) method of GRASP, or a fully variational method that optimizes all the radial functions that contribute to one or more ASFs of interest, such as the extended optimum level (EOL) method available in GRASP2K.
ASFs are labeled according to the composition of their wave function, where normally, the label is that of the CSF with the largest expansion coefficient.Such labels are not always unique and, ideally, when needed, should include two or more expansion coefficients, along with relative phases for unique identification, but such a complicated label is rarely adopted.For highly ionized atoms, composition in jj-coupling is usually more dominant in that the largest component represents a larger fraction of the composition.The composition of the wave function is relative to the CSFs included in a specific calculation.The least squares fitting method used by the HFR method of the Cowan code [19] achieves accurate energy levels from observed wavelengths by adjusting the values of Slater integrals.This adjustment will not correct for the missing angular symmetries in the wave function expansion.Thus, the composition for the ground state of W 37+ is 100% 4s 2 4p 6 4d in the Atomic Spectra Database (ASD) [20], whereas it is 98% in a more accurate ab initio calculation [9].
In some papers that use the early GRASP code, the ASFs of interest are the same as the set of the CSFs.Generally, when this method is used, the number of CSFs is sufficiently small, so that the entire interaction matrix can be stored in memory and a library routine used to compute all eigenvalues and their associated eigenvectors.In fact, the desired spectrum defines the CSF set, but since there is no interaction between CSFs of different J and parity, this means that permutations exist that divide the interaction matrix into distinct, non-interacting sub-blocks that require iterative methods to have initial estimates of eigenvectors of the appropriate J and parity [21].For large-scale calculations, there are several advantages to treating the sub-blocks as distinct eigenvalue problems; the calculations are faster and require less memory.GRASP2K supports both options, but the "block" version is recommended for large calculations.For large-scale calculations, the number of CSFs will vastly exceed the number of ASFs.
The accuracy of a wave function does not depend directly on the length of the expansion.Some CSFs are more important than others.In the 1980s to the 1990s, when computer memory was limited and computations were time consuming, the starting point was a single CSF that defined a Hartree-Fock wave function, and then, a few CSFs were added that accounted for the most important correlation contributions.In fact, non-relativistic Z-dependent perturbation theory defines the zero-order approximation as an expansion over the CSFs of the complex (CSFs with the same parity and principal quantum numbers) [22,23] for which the coefficients are an eigenvector of the interaction matrix.The first-order correction to the wave function then consists of all CSFs that interact with one or more CSFs in the zero-order approximation.In accurate GRASP2K calculations, including only the first-order correction, CSFs with a sufficiently large expansion coefficient form a multi-reference (MR) set.Typically, CSFs with an expansion coefficient larger than 0.1 in magnitude (or contributing 1% to the composition of the wave function) should be in the MR set.The size of this set determines both the accuracy and the computational resources needed.In practice, some CSFs in the complex may be unimportant for a particular case, and in others, CSFs outside the complex may be important, particularly in neutral atoms.
Since the interaction matrix element between CSFs that differ by more than two orbitals is zero, the interacting CSFs are contained in the set of CSFs generated by single-and double-(SD) substitutions of occupied orbitals to unoccupied or unfilled orbitals.(For states lowest in their symmetry, substitutions are excitations and often referred to as such, but for excited states, there may also be some that are de-excitations.)Thus, there are two important concepts: (1) the role of the complex, which may contain potentially strong interactions, and (2) the SD substitutions for other CSFs, whose expansion coefficients, in general, may be small, but are large in number.
An alternative to expanding the MR set is to include also second-order effects through triple (T) and quadrupole (Q) excitations, a process that may be applied to light atoms, but soon becomes unmanageable for complex ones.The process of transferring selected CSFs to the MR set, in effect, includes selected TQ excitations relative to the smaller MR set, in a first-order calculation for the larger set.In GRASP2K, the JJGEN [24] module generates CSFs for wave function expansions according to a number of rules, including SDTQ and higher excitations.
With every CI wave function, there is a set of orbitals that defines the CSF basis.The SD (or higher) process divides this set into the orbitals occupied in the MR set and other orbitals for the correction to the wave function.The latter are correlation orbitals, and associated CSFs may be large in number.In perturbation theory, summations are applied over a predetermined set of CSFs.In FAC, the needed orbitals are defined in terms of orbitals from a model potential and are essentially infinite in number, but in practice, are terminated in some manner.What is important for perturbation theory is the completeness of the basis.In variational methods, the radial functions for a set of orbitals are optimized so as to yield a stationary energy.The energy functional may be a linear combination of total energies for a set of ASFs.This is referred to as an extended optimal level calculation (EOL).A simpler strategy is to compute the radial functions from an energy functional that is defined in terms of the average energy of all CSFs in the wave function expansion of all ASFs.This method is referred to as the (EAL) method.The difference is that the latter omits the interaction between CSFs in the orbital optimization phase of the calculation and there is no large distinction between orbitals that are part of the MR set and other orbitals.In the EOL calculation, correlation orbitals are more contracted [25,26] and total energies lower than in a similar EAL calculation.
A systematic calculation [27] is one that computes wave functions of increasing expansion size that, hopefully, converge to an accurate approximation and, in the process, indicate the extent of convergence within the correlation model.This is done by increasing the set of correlation orbitals by layers, where, typically, the next layer has the extra orbital of each symmetry.Therefore, for Be, for example [26], where the complex defines the MR set as {1s 2 2s 2 , 1s 2 2p 2 } and the initial orbital set is {1s, 2s, 2p}, this could be increased to {1s, 2s, 2p, 3s, 3p, 3d}, etc.Such calculations are referred to as n = 2, n = 3 calculations, respectively, where n designates the largest principal quantum number.When orbital sets are restricted to a maximum orbital angular quantum number, they are referred to by both maximum quantum numbers as in 5f (indicating that 5g has been omitted).Therefore, a third factor affecting the accuracy of a CI calculation is: (3) the set of orbitals and the method by which they are derived.
For ground state energies, where only a single state is under consideration, the most accurate wave function in non-relativistic theory is the one with the lowest total energy.In relativistic theory, Breit and QED often raise the energy, so the comparison is more complex."Spectrum" calculations, namely calculations for energy levels relative to the ground state, are different in that only the energy differences need to be accurate.The energy of the common core largely cancels.In variational methods, it has been found advisable to have all orbitals for the core the same and treat core polarization and relaxation as a correction.Another factor in a systematic calculation is the need for maintaining a balance between the different states in the treatment of correlation, so that the energy difference is stable, converges and can be monitored.For this reason, correlation in all states should be treated in a uniform manner.Thus, if the substitution s 2 → p 2 is applied to one state, it should be applied to all states.SD excitations from a single CSF achieve a balance, but when applying the MR-SD strategy (for generality, the MR-SD method includes the case where the multi-reference set consists of only one configuration), the MR set for odd and even states should also be balanced.One strategy is to require the MR set to represent essentially the same percentage of the final composition of the wave function.This was the strategy used in determining the 2s 2 2p 2 P 3/2 -2s2p 2 4 P 5/2 excitation energy that determines the energy separation of the quartets relative to the doublets in neutral boron [28].

Example of an Analysis
As an example, let us compare some configuration interaction calculations for the 4s 2 4p 5 (odd) and 4s 2 4p 4 4d (even) ASFs of W 39+ .The notation implies a 1s 2 2s 2 2p 6 3s 2 3p 6 3d 10 core, and most calculations have treated this core as inactive.Because this is a heavy atom, a fully relativistic approach is needed, and Breit, vacuum polarization, and self-energy corrections should be included, although, as documented by Aggarwal, Jha, and Mohan [11], the Breit correction is the most important.It will be assumed that each MR set consists of a single CSF, namely 4s 2 4p 5 (odd) and 4s 2 4p 4 4d (even).
Table 1 shows the excitations that were included by various authors in defining the interaction matrix.Some single excitations from the core or from valence orbitals to 5s, 5p, 5d were included by some authors, but these do not have a large effect, so this table is restricted to the potentially significant excitations in the complex.It also shows whether the same excitation has been applied to both the odd and even states.A blank indicates the excitation was not applied to any state.Probably because of the large number of cases being considered, Fournier [1] restricted excitations to single excitations that are balanced.Radtke applied the 4p → 4f excitation only to the odd states.Aggarwal, Jha and Mohan (AJM) [11] included some double excitations not quite in a balanced fashion, whereas Aggarwal and Keenan [12] performed a series of GRASP calculations.The first, GRASP1 (G1), has only a few excitations and is almost balanced, whereas the second, GRASP2 (G2), includes more excitations, often only in the odd state.Also included were several odd configurations, such as 4s 2 4p4d 4 , 4s 2 4p 2 4d 2 4f and 4s 2 4p 3 4d 3 , which are higher-order corrections and do not interact directly with the 4s 2 4p 5 configuration.Thus, they do not affect the wave function significantly, but greatly increase the computational effort.Finally, the present GK 1 calculation using the GRASP2K code, included all SD excitations within the n = 4 complex.
Table 1.Excitations to the n = 4 orbital set from 4s 2 4p 5 (odd) and 4s 2 4p 4 4d (even) configurations in different calculations.For each calculation, the Yes/No (Y/N) notation indicates whether the excitation has been included in the odd and even states, respectively.A blank field indicates that the excitation was not considered for any state.
In W 39+ , the ground state is 4s 2 4p 5 and the first excited configuration is 4s 2 4p 4 4d.The 4s4d → 4p 2 excitation produces the 4s4p 6 configuration, which happens to be a configuration within the spread of 4s 2 4p 4 4d.In order to include all the J = 1/2 levels, orbitals were optimized on the six lowest levels.
Table 2 shows the results from four different calculations using GRASP2K.In GK 0 , the orbitals for the two levels of 4s 2 4p 5 were optimized for the equally weighted average energy for the two odd ASFs, 2 P 1/2 and 2 P 3/2 .With the core orbitals then fixed, the 4s, 4p, 4d orbitals for 4s 2 4p 4 4d and 4s4p 6 were determined for the equally weighted energy of six J = 1/2, eight J = 3/2, eight J = 5/2, five J = 7/2 and two J = 9/2 ASFs.The GK 1 calculation proceeded in the same fashion, except that the wave function expansion included SD excitations from 4s 2 4p 5 for the odd states and from 4s 2 4p 4 4d for the even to the n = 4 set of orbitals, which includes all the important CSFs within the complex.GK 2 differs in that the excitations include SD excitations to the larger n = 5 orbital set.The table shows that the energy spectra for both GK 1 and GK 2 are very similar.These calculations both include the effect of "valence correlation", where all the occupied orbitals in the n = 4 subshells are treated as valence electrons.They differ in that GK 2 also includes excitations to the n = 5 orbitals.Such calculations have neglected the effect of core polarization.In a CI calculation, the effect of core polarization arising from the 3d 10 subshell is represented by excitations, where one of the orbitals is from the 3d subshell and the other is a valence orbital.The GK CV results include the SD excitations from both valence and 3d-core-valence excitations, with Level 21 slightly lower than Levels 22 and 23, which have the same energy in the observed spectra.
Table 2. Comparison of energy levels (in units of 1,000 cm −1 ) for 4s 2 4p 5 , 4s4p 6 , 4s 2 4p 4 4d configurations in W 39+ from four systematic calculations and the final composition in LS and jj-coupling (see the text for details).In the latter, the J-value of 4s 2 4p 4 is given in square brackets to designate the coupling.Levels not in the final order are shown in blue.In Table 2, the levels are presented in the order of the most accurate calculation (in theory), namely GK CV .Levels not in this order are depicted in blue.This shows that the order of the levels without correlation (GK 0 ) is not always correct, although some levels are very close together and the order of such levels can be expected to be uncertain.The computed spectra from GK 1 and GK 2 are almost in the final order in that Levels 22 (J = 1/2) and 21 (J = 5/2) are very close.Only with the addition of core polarization are the levels in the final order.Levels 21-23 have been observed as parts of a blended line [4].
Also included in Table 2 is the largest component of the wave function composition for each level.For neutral atoms, LS coupling describes the composition well and has the advantage that it is easy to classify E1 transitions readily into spin-allowed and spin-forbidden transitions.However, for fully relativistic calculations in jj-coupling for highly ionized atoms, the largest component in LS coupling is often less than 50%.Since states are labeled by the largest component, this may lead to situations where two different ASFs have the same label, which occurs frequently when three CSFs interact strongly.Table 2 includes the largest component in both schemes and shows that, in most cases, the largest component in the jj-coupling scheme is more dominant and, thus, is a better representation of the ASF.A more detailed study of composition for the J = 1/2 states shows that, for Level 9, the 4s4p 6 CSF accounts for 14.4% of the wave function composition and should have been included in the MR set for even states, although it would affect only the even J = 1/2 levels.This illustrates that it is not always possible to determine the best MR set before a calculation is started when observed energy levels are not available.
Table 3 compares the J-values of levels 20-25 between the different theoretical calculations.In the GRASP2K calculations, the order of these levels changed significantly when the SD correlation in the complex was added in, going from GK 0 to GK 1 .The J = 1/2 and J = 9/2 levels were affected the most, and some of their positions changed.Note that other calculations with smaller expansions that exclude excitations to CSFs with more 4d and 4f electrons (see Table 1) are closer to that of the uncorrelated order of GK 0 .The final column shows the results from a least squares fitting of observed wavelengths using the HFR code [19] for which Levels 24 and 25 are interchanged.Because of their uncertainty, Levels 24 and 25 are not included in the ASD website [20].
Table 3. Table showing the order of J-values for Levels 20 to 25 for the different calculations.Some upper levels are connected with the 2 P o 3/2 ground configuration by a direct E1 transition, so that their energy levels can be derived from observed wavelengths.Table 4 compares experimental values with theoretical values from various calculations.The calculated value that agrees most closely with the experiment is shown in blue.In all cases, the present results are the most accurate.The fact that often this is GK 2 rather than GK CV , indicates that more effects need to be included.Generally, there is reasonable agreement between F [1], AJM [11] and G1 [12], since the expansions are similar in that the excitations with the largest effect, namely 4p 2 → 4d 2 , were not considered.Both AJM and G1 were obtained through the use of the same code [13,14] and the same EAL method for obtaining radial functions for the orbitals.It should be mentioned that Aggarwal and Keenan [12] labeled Level 13 as 2 F 5/2 , rather than 2 D 5/2 .Because of the lack of balance in their excitations, particularly in G2, most transition energies became larger in going from G1 to G2, increasing the discrepancy with the experiment.[11]; e Aggarwal and Keenan [12].Exp.: experiment.
Table 5 compares transition data from methods GK 2 , GK CV and AJM along with an indicator of accuracy, δT [29].The latter is the deviation from unity of the ratio of the length and velocity form of the line strength, which can be used as an accuracy indicator in addition to the transition energy.The latter may not always be available, in which case, the length/velocity test may still be applied in a statistical sense for a set of lines.Unlike the cancellation factor as defined by Zhang et al. [30] that neglects the cancellation in the radial factor of the transition matrix element, the δT factor is sensitive to all forms of cancellation.In Table 5, for most transitions, δT is considerably smaller for GK results than for AJM.An exception is the first transition to 4 D 5/2 .This is an intercombination transition, where the line strength for a relativistic calculation is exceptionally small, due to cancellation [31], a case for which all methods may have more uncertainty.In spite of the much larger δT in this one case, the average value for these 10 transitions is 0.02, 0.025 and 0.041 for the three methods, GK 2 , GK CV and AJM.What is a concern, however, is that the average δT is larger for GK CV than for GK 2 .At the same time, the average δT for the present values are smaller than the average for AJM values, in agreement with the energy level evaluation.A better understanding of the source of the remaining discrepancy is needed.The omission of 4s4p 6 from the even MR set has not resulted in a noticeable discrepancy, because, in this particular case, the n = 4 expansion does not change, but its inclusion would result in a more reliable n = 5 calculation.
Lines from M1 transitions between the fine-structure levels of the ground state have been observed [7].Quinet [15] recently has computed many energy levels using the GRASP [13,14] code in the EAL mode with SD excitations in the n = 4 shell for the 4p k (k = 1 − 5) and 4d k (k = 1 − 9) of tungsten ions W 29+ through W 43+ .Two experimental measurements have been reported for W 39+ , one at 131.8 Å [5] and one at 134.74 Å [7] for which the 2 P 1/2 energy levels are 758,725 cm −1 and 742,170 cm −1 , respectively.Table 6 shows fine-structure splitting of a number of calculations.All are closest to the Ralchenko et al. [7] value, but the inclusion of some core-valence correlation did not bring results into better agreement with observation.At the moment, the source of the remaining discrepancy is not clear.As shown by Aggarwal, Jha and Mohan [11] and confirmed in the present work, the largest correction to the energy levels is the Breit correction, with vacuum polarization and self-energy being significantly smaller.
Table 5.Comparison of E1 transition data for transitions from the 4s 2 4p 5 2 P o 3/2 ground state to selected upper levels for GK 2 , GK CV and AJM.Experimental values of λ are provided when available [4], and the method in best agreement with that which is observed is displayed in blue.Highly ionized tungsten and similar elements are not super-heavy elements, where Breit and QED corrections are more important than the many-body effects.At the same time, they are sufficiently ionized, so that correlation effects converge fairly rapidly, requiring only a few layers of correlation orbitals.As a result, the method for obtaining the radial functions of orbitals may not be as critical as for neutral atoms.Certainly, in the present case, the determination of the orbitals for an n = 4 expansion (within the complex) from an EAL calculation in GRASP or model potential in FAC, in comparison with EOL for GRASP2K, does not appear to make a significant difference.Greater differences are expected to occur when the orbital sets are increased to n = 5 or higher.
The most accurate results were the MR-SD results obtained using GRASP2K.The size of the matrices were as follows: The inclusion of core-valence increases the size of the interaction matrices significantly, but with sparse-matrix methods and iterative methods for finding the desired eigenvalues, the calculations can readily be performed on current computers.

Conclusions
In summary, several important factors that affect the accuracy of results from different computational procedures based on CI methods have been identified.In the W 39+ example, the most important factor was the CSFs included in the wave function expansion.It needs to be remembered that when two CSFs, interact, say i and j, the ratio, c j /c i , of their expansion coefficients is, where E is an eigenvalue of the 2 × 2 interaction matrix, (H ij ).Assuming H ij ≡ 0, then, for the order, such that c i > c j , c j is large (there is a lot of mixing) either if H jj is close to E (near degeneracy) or if H ij is large.Including CSFs close in energy is not sufficient.Also important are CSFs with a strong interaction (a large H ij ).In quantum chemistry, correlation is classified as static (nondynamic) correlation or dynamic correlation [32].In atomic physics, static correlation arises from near degeneracies, whereas dynamic correlation includes the correction needed for the electron-electron cusp condition and is a short-range effect [33].The MR-SD method of generating expansions with sufficient layers of orbitals will include both types.In our example, the expansions from a single configuration (4s 2 4p 4 4d) showed that 4s4p 6 was an important component in the J = 1/2 levels, and the calculation would be improved if 4s4p 6 were added to the MR set.
A comparison with experimental energy levels (or wavelengths) for even a few levels helps the evaluation of computed wave functions, but if these are not available, the internal check on the agreement between length and velocity forms of line strengths can often be useful.

Table 4 .
Comparison of energy levels (in units of 1,000 cm −1 ) from two different experimental measurements with a number of different calculations.The best values are denoted in blue.

Table 6 .
The fine-structure splitting (in cm −1 ) of the 4s 2 4p 5 2 P o configuration as predicted by theory compared with experiment.
a SDTQ substitutions for n = 4 and SD for n = 5.