Basis set calculations of heavy atoms

Most modern calculations of many-electron atoms use basis sets of atomic orbitals. An accurate account for the electronic correlations in heavy atoms is very difficult computational problem and optimization of the basis sets can reduce computational costs and increase final accuracy. Here we suggest a simple differential ansatz to form virtual orbitals from the Dirac-Fock orbitals of the core and valence electrons. We use basis sets with such orbitals to calculate different properties in Cs including hyperfine structure constants and QED corrections to the valence energies and to the E1 transition amplitudes.


I. INTRODUCTION
Many efficient methods of atomic calculations, such as configuration interaction (CI), many-body perturbation theory (MBPT), coupled cluster (CC), or their combinations, employ basis sets of one-electron orbitals.For example, the methods CI+MBPT and CI+CC [1][2][3][4][5] use CI for valence electrons, where the electronic-correlations effects are strong, and account for weaker core-valence correlations by means of either the second order MBPT, or (linearized) CC method.In these methods, the relativistic effects are treated within the no-pair approximation for the Dirac-Coulomb-Breit Hamiltonian.On the other hand, QED effects may be approximately accounted for using the model-QED-potential approach [6][7][8].For heavy polyvalent atoms such calculations become very computationally expensive [9].That is why it is very important to develop efficient basis sets, which provide high accuracy at reasonable length.
The Dirac-Fock (DF) method serves as initial approximation for most of the atomic calculations.However in practice, it is useful to keep DF orbitals in the basis set without re-expanding them in some nonphysical basic orbitals.On the other hand, the highly excited DF orbitals are usually ineffective in accounting for the correlations between valence electrons, as their radius grows too rapidly.Because of that the B-splines [10], Sturm orbitals [11], or other simple orbitals [12,13] usually turn out to be more useful.Therefore, an efficient basis set has to include different types of orbitals.An effective method to merge two subsets of orbitals in one basis set was suggested in Ref. [14].The first subset consisted of DF orbitals for core and valence shells and the second subset consisted of B-splines, whose parameters were optimized to complement the first subset.Such mixed basis sets were tested for calculations of the valence energies of Au and Fr and sufficiently fast convergence was observed: saturation was reached for 25 splines per partial wave [14].
It is known, that accurate calculations of the hyperfine structure and parity non-conservation effects for heavy atoms are much more challenging then calculations of the energies, or transition amplitudes.These properties depend on the wave function at small distances.Correla-tions change the behaviour of the valence-electrons wave functions in this region and lead to the large corrections to the matrix elements of these operators.Because of that the convergence with respect to the number of orbitals in the basis set is typically much slower.Test calculations of the magnetic hyperfine constant A for Au did not converge even for the mixed basis set with 45 B-splines per partial wave [14].It was found that small components of the valence orbitals near the origin were not smooth.Increasing the size of the basis set decreased the amplitude of these non-physical dips or bumps, but shifted them closer to the origin, where the hyperfine interaction was particularly strong.The solution suggested in that paper was to add to the basis set DF orbitals for the ion Au M + , or, in other words, the orbitals for the V N −M potential, where N is the number of electrons in the neutral atom.Such orbitals are similar to those of the neutral atom, but are more contracted.That allowed to smoothly change the orbitals near the nucleus.The calculation of the hyperfine constants A for Au already converged for two such ionic orbitals and 24 B-splines per partial wave.To make this recipe work one need to choose M such that ionic orbitals are neither too different, nor too similar to those of the neutral atom.In the former case such orbitals become less useful, while in the latter one can run into the linear dependency problem.In Ref. [14] the orbitals for the V N −6 potential were used.Here, we suggest another method to form additional orbitals to supplement DF orbitals and B-splines.The new procedure is more formal and does not require arbitrary adjustments.We test these basis sets for calculation of the QED corrections and hyperfine constants in neutral Cs.

II. METHOD
The basis set for atomic calculations consists of subsets for different partial waves, which are defined by the relativistic quantum number κ = (l − j)(2j + 1), where l and j are the orbital and total angular momenta.The relativistic two-component radial wave functions have the form arXiv:2312.07782v1 [physics.atom-ph]12 Dec 2023 Most atomic calculations are done in the no-pair approximation.QED corrections may be included approximately by means of the model operators [7,8,15].This means that the Hamiltonian is projected on the positive energy electronic states and the basis set does not include negative energy continuum.To this end, the small components Q of the virtual basis orbitals can not be chosen arbitrary.Most often they are found from the kinetic balance condition [16]: where α ≈ 1 137 is the fine structure constant (we use atomic units throughout the paper).Below we assume that the small components are always formed with the help of Eq. ( 2) and focus on the large components P .
When some small perturbation is added to the Hamiltonian, the large component P of the electronic orbital slightly changes.These changes do not affect the asymptotic behaviour at large and small distances and the number of nodes of the function P remains the same.To account for such changes we need to add basic functions with similar properties.Let us consider a simple stretching transformation: where k = 1 + ε, |ε| ≪ 1. Expanding P in powers of ε we get: Thus, if we want to guarantee that stretched DF orbital Pnκ can be accurately expanded in the basis set, we can add a virtual orbital: Note that this function has the same asymptotic behavior at r → 0 and at r → ∞ as P nκ , however it has an extra node.

III. TEST CALCULATIONS FOR CS
To illustrate usefulness of the ansatz (5) let us consider QED corrections in neutral Cs.Firstly, we used a sufficiently large basis set to diagonalize the Dirac-Fock operator together with the model QED potential of Ref. [15,17].After that we calculated the difference between large components with and without QED corrections and compared them with the respective scaled orbitals (5).
Figure 1 shows that the latter very much resemble the former.Thus, we can expect that adding orbitals (5) to the basis set can significantly improve its quality and speeds up the convergence.We formed the basis set for Cs which included DF orbitals 1s − 7s, 2p − 7p, and 3d − 6d.For each DF orbital we added respective virtual orbital (5).To have the same number of virtual orbitals per partial wave, we formed one additional p orbital per partial wave and three additional d orbitals per partial wave using the method of Ref. [13].Thus, this basis set included 7 virtual orbitals per partial wave.This calculation was compared with calculation on the long basis set, which included these 7 plus 23 additional virtual orbitals per partial wave.Below, we refer to these basis sets as B7 and B30, respectively.We see that the QED correction for the valence 6s orbital is practically the same for both basis sets.For the orbital 1s, both the basis sets produce some non-physical oscillations, but their amplitude for the large basis set is few times smaller.
In the calculations of neutral atoms we usually explicitly need only wave functions of the valence electrons, where our short basis set seems to work nicely.To test it further we calculated QED corrections to the valence energies and to the E1 transition amplitudes.In both cases we used model QED potential [15].This model potential accounts for the QED corrections to the wave functions.For E1 amplitudes there are also QED corrections to the vertex, which are not included here.They were calculated in Refs.[18,19] and found to be much smaller.Results obtained with the basis sets B7 and B30 are compared in Table I.We see that two calculations give practically identical QED correction to the binding energy of the 6s electron.For the other partial waves the corrections to the energies are much smaller and the difference between results is about 1 -2%.For the s − p transition amplitudes the agreement is again very good, but for p − d transitions the difference between two basis sets is about 10%.On one hand, these corrections are smaller.On the other hand, in d wave we have only four orbitals formed with the ansatz (5), while other three orbitals of the basis set B7 are apparently less useful.We mentioned earlier that convergence in calculations of the hyperfine constants A is often very slow.Because of that we compared random-phase approximation (RPA) corrections for these constants obtained with basis sets B7 and B30.For the orbitals 6s, 6p 1/2 , and 6p 3/2 the difference was around 1%. II show that the basis sets, which include virtual orbitals (5) provide fast convergence for calculations of the QED and RPA corrections, which are particularly sensitive to the quality of the basis set.Below we compare QED corrections calculated with these basis sets with calculations of other groups and study potential sources of the differences.We focus primarily on the QED corrections to the E1 transition amplitudes.Next two columns give QED or RPA corrections calculated on the basis set B7 with 7 virtual orbitals per partial wave and basis B30 with 30 virtual orbitals per partial wave.Last column gives the difference between two calculations in percent.

Property
Units DF value QED correction Basis B7 Basis B30 Diff.QED corrections to the E1 amplitudes in Cs were calculated several times.Sapirstein and Cheng [19] calculated 6s 1/2 −6p 1/2 amplitude within one-determinant approximation, while calculations [20][21][22] included correlations, but neglected the vertex correction.Traditionally QED corrections to the E1 amplitude are expressed in terms of the parameter R:
There are two contributions to the QED correction R: where R PO accounts for the correction to the orbitals (perturbed orbitals) and R V is the vertex correction.
The former is supposed to be larger than the latter [7] and it can be approximated by model QED potentials.At present, the vertex correction can be calculated only within the consistent QED approach.
In Table II we compare our results for R PO with results of other groups.Sapirstein and Cheng [19] made ab initio QED calculation using local Kohn-Sham potential.They calculated both corrections and found, that the vertex correction R V was significantly smaller, R V = −0.065.Other calculations used non-local DF operator and model potential [7,23,24].Here we made calculations with model potentials [15] and [7], which we denote as STY and FG potentials respectively.
Table II shows that the results obtained with two QED potentials agree within 2 -3%.The difference between different groups is somewhat bigger, about 4 -6%.There are small differences in the model potentials described in Refs.[7] and [23,24].Incompleteness of the basis FIG. 2. Comparison of the QED calculations with the short (B7) and the long (B30) basis sets.The short basis set includes 7 orbitals per partial wave.In the s wave all of them are formed using ansatz (5), while in other partial waves additional orbitals were added using method [13].
TABLE II.Comparison of the PO QED corrections RPO (see Eq. ( 7)) to the E1 transition amplitudes in Cs with other calculations.Our calculations were done using model potentials STY [15] and FG [7].

Transition
This Because of that, it is important to check, how much QED corrections depend on the choice of the potential.Note also that parametrization (6) allows to partly compensate for the differences in the value of the initial amplitude E1 0 in different approximations.

B. Local screening potentials
We made calculations with the model QED operator [15] and three local screening potentials: the first one is the core-Hartee potential (CH) induced by the core electrons of Cs + ground-state configuration, the other two are Kohn-Sham (KS) [25] and Slater (S) [26] potentials corresponding to the ground-state configuration of the neutral Cs.The proper asymptotic behavior of the potentials KS and S is restored by including the Latter correction [27].In Table III the results obtained using the local screening potentials are compared with the Dirac-Fock results.We see that calculations of QED corrections δε with local potentials do not correlate with the DF cal-culations for all partial waves, but the s wave.The reason for this discrepancy is obvious.In the DF approximation the energy correction comes from the direct contribution of the perturbation and the contribution from the adjustment of the self-consistent field (relaxation of the core) δV SCF .The former is large for the s wave and rapidly decreases with the orbital quantum number l.At the same time, the relaxation of the core affects all valence orbitals in a similar way.Two contributions become comparable already for the valence p orbitals.For the local potentials there is only direct contribution.The KS potential reproduces valence orbital energies ε better than two other potentials, while QED corrections δε to the s wave energies for the CH and KS potentials are of the similar accuracy.Slater potential is clearly the least accurate.
Taking into account that core relaxation is very important for the p and d waves, the local approximation is clearly inapplicable for the QED corrections to the E1 transitions between p and d orbitals.The results of the calculations for the s ↔ p transitions are presented in Table IV.Again, the KS potential gives better agreement with DF, than two other potentials.In most cases the difference is less, or about 10%.The only exception is the 6s 1/2 − 7p 1/2 transition, where the difference is 24%.For the CH and S potentials the average differences are about 14% and 22% respectively.TABLE III.QED corrections δε to the binding energies ε in Cs calculated in DF approximation and with three local potentials, core-Hartree (CH), Kohn-Sham (KS), and Slater (S).All calculations were done using model QED potential [15] (STY).
Orbital Dirac-Fock core-Hartree Kohn-Sham Slater ε δε ε δε ε δε ε δε (a.u.) (cm −1 ) (a.u.) (cm −1 ) (a.u.) (cm −1 ) (a.u.) (cm −1 ) 6s We conclude that for Cs the KS potential provides higher accuracy for the QED calculations than CH and S potentials.Note, that the result of Sapirstein and Cheng [19] for the 6s 1/2 − 6p 1/2 transition cited in Table II was obtained using KS potential.Our KS value is closer to their value, then the DF one.Thus, the approximately 20% difference between ab initio calculation [19] and calculations with model QED potentials presented in Table II can be attributed in part to the local potential error.

IV. CONCLUSIONS
Many numerical methods for atomic calculations require one-electron basis sets.Here we suggested the differential ansatz to generate virtual orbitals from the DF orbitals of the core and valence electrons.The number of such orbitals cannot exceed the number of DF orbitals.More virtual orbitals can be added using B-splines, or other traditional methods.Our test calculations for Cs show good convergence of calculations of different properties, including QED corrections and hyperfine structure constants.
We also tested three local potentials, which are often used for the QED calculations in many-electron atoms and found out that for the lower partial waves the Kohn-Sham local potential demonstrated best agreement with the DF potential.For the higher partial waves the results obtained with all tested local potentials are very different from the results obtained with the non-local DF potential.These differences are caused by the absences of the contribution from the change of the self-consistent field induced by the perturbation.This effect becomes dominant for the partial waves with l > 1 making calculation of the respective QED corrections problematic.In order to improve the accuracy of the ab initio QED calculations in heavy atoms it is necessary to include core-relaxation effects by going to the next order of the perturbation theory in the electron-electron interaction.The next step in quantum-mechanical calculations of the QED corrections requires to construct effective operators for the vertex contribution.

FIG. 1 .
FIG. 1.Comparison of the differential ansatz (5) with the difference between DF orbitals for the Dirac-Coulomb Hamiltonian with and without model QED potential.

TABLE IV .
[15]arison of the PO QED corrections RPO to the s ↔ p transition amplitudes in Cs for different potentials.For local potentials we also give the difference with DF results in percent.All calculations were done using model QED potential[15](STY).