Application of the dualkineticbalance sets in the relativistic manybody problem of atomic structure
Abstract
The dualkineticbalance (DKB) finite basis set method for solving the Dirac equation for hydrogenlike ions [V. M. Shabaev et al., Phys. Rev. Lett. 93, 130405 (2004)] is extended to problems with a nonlocal sphericallysymmetric DiracHartreeFock potential. We implement the DKB method using Bspline basis sets and compare its performance with the widelyemployed approach of Notre Dame (ND) group [W.R. Johnson and J. Sapirstein, Phys. Rev. Lett. 57, 1126 (1986)]. We compare the performance of the ND and DKB methods by computing various properties of Cs atom: energies, hyperfine integrals, the paritynonconserving amplitude of the transition, and the secondorder manybody correction to the removal energy of the valence electrons. We find that for a comparable size of the basis set the accuracy of both methods is similar for matrix elements accumulated far from the nuclear region. However, for atomic properties determined by small distances, the DKB method outperforms the ND approach. In addition, we present a strategy for optimizing the size of the basis sets by choosing progressively smaller number of basis functions for increasingly higher partial waves. This strategy exploits suppression of contributions of high partial waves to typical manybody correlation corrections.
pacs:
03.65.Pm, 31.30.JvI Introduction
Applications of perturbation theory in quantum mechanics require summations over a complete set of states of the lowestorder Hamiltonian. Usually, the relevant spectrum is innumerable. In practical applications such eigenspectra are often modeled using finite basis sets, chosen to be numerically complete. Since the sets are finite, the otherwise infinite summations become amendable to numerical evaluations.
The use of a finite basis set composed of piecewise polynomials, socalled Bsplines de Boor (2001), has proven to be particularly advantageous in atomic physics and quantum chemistry applications Bachau et al. (2001). In this approach, an atom is placed in a sufficiently large cavity and the atomic wavefunctions are expanded in terms of the underlying Bspline set. Further, the variational Galerkin method is invoked and the solution of the resulting matrix eigenvalue problem produces a quasispectrum for the atom. In nonrelativistic calculations, the lowestenergy orbitals of the resulting basis set closely agree with those of the unperturbed atom, and calculations of various properties of the lowlying states can be carried out. In particular, one could generate singleparticle orbitals in some suitable lowestorder approximation, and use the resulting basis set in applications of manybody perturbation theory (MBPT) .
Application of the outlined approach to the relativistic problems brings in a complication– the appearance in the atomic quasispectrum of nonphysical “spurious” states. These states appear in the solution of the singleparticle radial Dirac equation for angular symmetry, ( orbitals). The spurious states rapidly oscillate and, moreover, spoil the mapping of the generated quasispectrum onto the lowenergy orbitals of the atom. At the same time they are required for keeping the set complete. The problem of spurious states was discussed in the literature in details, see e.g., Ref. Goldman (1985); Quiney et al. (1987); Johnson et al. (1988); Igarashi (2006), and several solutions were proposed. In the pioneering applications of the Bsplines in relativistic manybody problem by the Notre Dame group, Johnson et al. (1988) added an artificial potential spike centered at the origin to the Hamiltonian matrix. The overall effect was to shift the spurious states to higher energies thus restoring the lowenergy mapping to the physical states. We will refer to the sets generated using this prescription as the Notre Dame (ND) sets. Another solution was to use “kineticallybalanced” sets Quiney et al. (1987), which related the small and the large components of the basis set functions via the Pauli approximation. Recently, an extension of this method was proposed in Ref. Shabaev et al. (2004). Here, due to additional relations between the small and large components, the negative (Dirac sea) and positive energy spectra are treated in a symmetric fashion. To emphasize this builtin symmetry, the authors refer to their method as the “dualkineticbalance” (DKB) approach. In both methods (by contrast to the ND prescription), the spurious states were shown to be completely eliminated from the quasispectrum.
Motivated by the success of the DKB method in computing properties of hydrogenlike ionsShabaev et al. (2004); Igarashi (2006), here we investigate the suitability of the DKB method in modeling the spectrum of the (nonlocal) DiracHartreeFock (DHF) potential. We compare the performance of the ND and DKB methods by computing various properties of Cs atom: singleparticle energies, hyperfine integrals, and the paritynonconserving amplitude of the transition. We find that for properties involving matrix elements accumulated near the nucleus, the DKB method outperforms the ND approach. Otherwise, if the electronic integrals are accumulated far from the nucleus, both methods produce results of a similar quality.
In addition, we investigate a possibility of using varying number of basis set functions for different angular symmetries. Summations over intermediate states in expressions of perturbation theory are carried out both over angular quantum numbers and for fixed over radial quantum numbers. Usually, as (and ) increases, the correlation corrections due to higher partial waves become progressively smaller. Intuitively, one expects that for higher partial waves it would be sufficient to use smaller radial basis sets of lesser quality. This would reduce storage requirements for manybody calculations (for example, in implementing coupledcluster formalism ) and would speed up numerical evaluations. While such an approach is common in quantum chemistry, see, e.g., Ref.Davidson and Feller (1986), the question of building the optimal Bspline basis set was not addressed yet in relativistic manybody calculations. We illustrate optimizing the basis sets by computing the secondorder energy correction for several states of Cs.
This paper is organized as follows. First we recapitulate the Galerkintype approach to generating a finitebasis set quasispectrum for the Dirac equation in Section II. The variational method is invoked for relativistic action and the problem is reduced to solving the generalized eigenvalue problem in Section II.1. The DHF potential is specified in Section II.2. Further we specify ND and DKB sets in Section II.3 and boundary conditions in Section II.4. A numerical analysis of Cs atom is provided in Section III. In Sections III.1 and III.2 we compare the performance of the ND and the DKB sets by computing singleparticle energies, hyperfine integrals (Section III.1), and paritynonconserving amplitudes (Section III.2) using both methods. The spurious states arising from the ND method are examined in Section III.3. In Section III.4 we consider secondorder energy corrections in the DKB method and discuss a strategy of optimizing the size of the basis set.
Ii Problem setup
We are interested in solving the eigenvalue equation for the Dirac Hamiltonian, , where is the nuclear potential and is the meanfield (DiracHartreeFock) potential. is in general a nonlocal potential. Assuming that both potentials are central one may exploit the rotational invariance to parameterize the solutions as
(1) 
with being the spherical spinors. The solutions depend on the radial quantum number and the angular quantum number . The large, , and small, , radial components satisfy the conventional set of radial Dirac equations
These radial equations may be derived by seeking an extremum of the following action Johnson et al. (1988)
(2)  
where the dependent Pauli operators are defined as
Radial integrals here and below have implicit limits from to , where is the radius of the confining cavity. Boundary conditions may be imposed by adding the term to the action. The term controls the appearance of the spurious states in the quasispectrum. We will specify these two terms below.
ii.1 Reduction to the matrix form
We employ two finite basis sets and , which, since we are interested in solving the (angularlydecoupled) radial equations, may depend on the angular quantum number . We may expand the large and small components in terms of these bases
(3)  
the expansion coefficients being the same for both the large and the small components. The above expansions are plugged into the action, Eq. (2), and, further, its extremum is sought by varying the expansion coefficients. As a result one arrives at the following generalized eigenvalue equation
(4) 
where and are matrixes and is the vector of expansion coefficients in Eq.(3). The matrix elements are given by
(5) 
The matrixes entering the definition of correspond to various pieces of the radial Dirac equations,
(6)  
(7)  
(8) 
with the matrix elements of the DHF potential, , given below. The terms and arise from the boundary and “spurious state” corrections in the action. Finally, the matrix , reflects the fact that the basis sets may be nonorthonormal.
ii.2 Potentials
The nuclear potential is generated for a nucleus of a finite size. We employ the Fermi distribution with the nuclear parameters taken from Ref. Johnson and Soff (1985). As to the DiracHartreeFock potential, we employ the frozencore approximation. In this method, the calculation is carried out in two stages. First, the core orbitals are computed selfconsistently. Second, based on the precomputed core orbitals, the DHF potential is assembled for the valence orbitals and the valence orbitals are determined. In the valence part of the problem, the core orbitals are no longer adjusted. Explicitly, for a set of the angular symmetry ,
(9)  
(10)  
(11) 
with the conventionally defined multipolar contributions
(12) 
and the angular coefficient
(13) 
ii.3 ND and DKB sets
Now we specify the Notre Dame (ND) and the dualkineticbalance (DKB) basis sets. Both operate in terms of Bspline functions and first we recapitulate the relevant properties of these splines. A set of Bsplines of order is defined on a supporting grid . Usually, the gridpoints are chosen as
In our calculations the intermediate gridpoints are distributed exponentially. Bspline number of order , , is a piecewise polynomial of degree inside . It vanishes outside this interval. This property simplifies the evaluation of matrix elements between the functions of the basis set. In addition, we will make use of the fact that as , the first splines behave as (all the remaining splines are zero)
(14) 
and at , all splines vanish except for the last spline, .
The Notre Dame set is defined as
(15)  
It corresponds to an independent expansion of the large and small radial components into the Bspline basis. The DKB set involves the Pauli operators and enforces a “kinetic balance” between contributions to the components:
(16)  
Notice that, as discussed below, to satisfy the boundary conditions, we will use only a subset of the entire DKB basis.
ii.4 Spurious states and boundary conditions
With the ND set, the spurious states are shifted away to the highenergy end of the quasispectrum by adding the following action Johnson et al. (1988) for and for . This correction may be seen as arising from an artificial functionlike potential centered at the origin. Unfortunately, as shown below in numerical examples, this additional spike perturbs the behavior of the orbitals near the nucleus. (As discussed below also sets the boundary conditions at ) The DKB set does not have the spurious states at all, so that .
We need to specify boundary conditions at and at the cavity radius, . We start by discussing the boundary conditions at . For a finitesize nucleus the radial components behave as
(17)  
In the Notre Dame approach, the boundary conditions are imposed variationally by adding the boundary terms to the action integral. Varying , Ref.Johnson et al. (1988), effectively reduces to . In practice, because of the variational nature of the ND constraint, the large component, although being small, does not vanish at the origin, and the limits, (17), are not satisfied. Alternatively, FroeseFischer and Parpia (1993) proposed to impose by discarding the first Bspline of the set (this is the only Bspline that does not vanish at ). This is a “hard” constraint, since the wavefunction, Eq.(3), would vanish identically at the origin. In our calculations, because of our motivation in building the smallest possible basis set, we extend this scheme further. We exploit the powerlaw behavior of the Bsplines, Eq.(14), and match it to the small limits (17). To satisfy the matching, we need to include the Bsplines starting with the sequential number ( must be smaller than the order of the splines ).
(18) 
For and and this is equivalent to the boundary condition of Ref. FroeseFischer and Parpia (1993). For higher partial waves, however, an increasingly larger number of splines is discarded: e.g., for . One should notice that for a basis that includes partial waves , for a faithful representation of the small behavior in all the partial waves, one needs to require the order of the splines to be at least of . In particular, for , .
When the first Bspline of the set, , is included in the basis (as in the ND approach), there is another difficulty in the calculations: since it’s value does not vanish at , the matrix elements and , which contain matrix elements of , are infinite in absolute value Igarashi (2006). In practical calculations, one uses Gaussian quadratures to evaluate this integral, so the result of the integration is finite. Yet this introduces arbitrariness in the ND calculations and may be a reason for a relatively poor representation of the orbitals near the origin.
At the cavity radius, to avoid overspecifying the boundary conditions for the Dirac equation, the ND group used the boundary condition . As with the conditions at the origin, this relation was “encouraged”variationally. In our calculations (similarly to Ref. Shabaev et al. (2004)) we use the “hard” condition by removing the last Bspline from the set. Examination of the resulting orbitals reveals that the wavefunctions acquire a nonphysical inflection towards the end of the supporting grid, while the ND orbitals behave properly. Further numerical experimentation, however, shows that the inflection does not degrade the numerical quality at least for the atomic properties of the lowlying bound states of interest. At the same time, throwing away the last Bspline reduces the number of basis functions and leads to a more compact set.
To summarize, we will use the DKB basis set that includes Bsplines with sequential numbers , Eq. (18) to . We will simply refer to this choice as the DKB basis. When we refer to DKB functions for a given partial wave, it would imply larger underlying Bspline set, e.g., for waves the total number of Bsplines would be . For the ND basis , as all Bsplines participate in the expansion.
In the remainder of this paper we present results of numerical analysis for Cs atom. It is an atom with a single valence electron outside a closedshell core. As the first step we carry out finitedifference DiracHartreeFock calculations for Cs core. The core orbitals are fed into the spline code where they are used to compute the matrix elements of the potential, Eq.(9). As in Ref. Johnson et al. (1988) the numerical accuracy is monitored by comparing the resulting quasispectrum with the DHF energies from the finitedifference code.
Iii Numerical examples for Cs atom
Here we provide numerical examples involving both ND and DKB sets for Cs atom. In the two Sections immediately following we compare the performance of the ND and the DKB sets. We generate the quasispectrum using both ND and DKB sets and carry out comparisons for singleparticle energies and hyperfine integrals in Section III.1 and paritynonconserving amplitude in Section III.2. Section III.3 contains an analysis of spurious states in the ND method. In Section III.4 we analyze secondorder energy corrections and discuss a strategy of optimizing the size of the basis.
iii.1 Energies and hyperfine integrals
We compare ND and DKB quasispectrums with energies obtained using a finitedifference DHF code for the lowlying valence states in Table 1. The ND and DKB calculations were carried out using basis functions for Bsplines of order . We used a cavity of radius bohr. For the cavity of this size, only a few lowestenergy orbitals remain relatively unperturbed by the cavity. From examining Table 1, it is clear that both ND and DKB sets have a similar accuracy for energies.
In the second part of Table 1 we compare values of the radial integrals entering matrix elements of the hyperfine interaction due to the electric (EJ) and magnetic (MJ) multipolar moments of a pointlike nucleus
The angular selection rules require . We use identical integration grid for all three cases (DHF,ND,DKB) listed in the Table. The grid is sufficiently dense near the origin, it contains about 100 points inside the nucleus. The numerical integration excludes the first interval of the grid. From examining the Table we find that the DKB set outperforms the ND basis. While the ND set still recovers twothree significant figures for the magneticdipole coupling, it produces wrong results for electricquadrupole and magneticoctupole integrals. Certainly, the accuracy in the ND case improves for a larger basis set, but larger basis sets come at an additional computational cost. We carried out similar comparisons for matrix elements of the electricdipole operator. As for the energies, we find that both the ND and DKB sets perform with a similar numerical accuracy.
As we have mentioned, variationally encourages the boundary condition . However, there is no such explicit encouragement for . Here we qualitatively discuss the observed properties of the radial components near the origin for the Cs ND set. We see that, though their general behavior is to approach zero, the large component wavefunctions often have small improper inflections or oscillations near the origin. Small component wavefunctions, on the other hand, often do not even approach zero. These improper behaviors seem to be exemplified as we look at states higher in the spectrum. As we have seen here, such improper behavior of both the large and small component radial functions near the origin prove detrimental for properties accumulated near the nucleus. We conclude that while producing the results of a similar quality for matrix elements accumulated far from , the DKB set provides a better approximation to the atomic orbitals near the nucleus.
State  Set  Energy  M1 HFI  E2 HFI  M3 HFI 

FD  
DKB  
ND  

FD  
DKB  
ND  

FD  
DKB  
ND  

FD  
DKB  
ND  

FD  
DKB  
ND  

FD  
DKB  
ND 
iii.2 Paritynonconserving amplitude
So far we examined properties of the individual basis set orbitals with sufficiently low energies. The real power of the finite basis set technique lies in approximating the entire innumerable spectrum by a finite size quasispectrum. This is important, for example, in computing sums over intermediate states (Green’s functions) in perturbation theory.
From the discussion of the preceding Section it is clear that the difference in quality of the ND and DKB basis sets is expected to become most apparent for the properties accumulated near the nucleus. Here, as an illustrative example, we consider the paritynonconserving (PNC) amplitude for the transition in . This amplitude appears in the second order of perturbation theory for the otherwise forbidden dipole transition and it can be represented as a sum over intermediate states
(19)  
Here and are electricdipole and weak interaction (pseudoscalar) operators, and are atomic energy levels. We will compute this expression in the singleparticle approximation. Specifically, the index runs over the entire DHF quasispectrum for partial wave, including core orbitals. The weak Hamiltonian reads
(20) 
where is the Fermi constant, is the weak charge, is the Dirac matrix (it mixes large and small components), and is the neutron density distribution. For consistency with the previous calculations the is taken to be the proton Fermi distribution of Ref. Blundell et al. (1992). Notice that the matrix elements of the weak interaction are accumulated entirely inside the nucleus.
We evaluate the sum (19) using the DKB and the ND sets with basis functions of order . The integration grids are the same in both cases and include large number of points () inside the nucleus. The PNC amplitude is conventionally expressed in units of , where is the number of neutrons in the nucleus of Cs. In these units, the results are
The finitedifference value is taken from Ref. Blundell et al. (1992). Again we note that the DKB set offers an improved performance over the Notre Dame set. Reaching the comparable accuracy in the ND approximation requires a larger basis set. For example, ND set reproduces the DKB result for the PNC amplitude.
Further insights may be gained from examining individual contributions of the intermediate states in the PNC amplitude. We plot individual contributions of the intermediate state to the PNC amplitude in Fig. 1. Both terms in Eq.(19) are included. We computed the data using ND and DKB basis sets generated in a cavity of bohr. From the plot, we observe that the dominant contribution arises from the lowlying valence states. As increases, the contributions become quickly suppressed (there is 10 orders of magnitude suppression for ). This is due to both increased energy denominators and decreased density at the nucleus for high . Comparison between the basis sets reveals that their contributions are identical until . For higher principal quantum numbers, the DKB contributions monotonically decrease, while ND contributions become irregular. Moreover, at the ND contributions start to flip signs with increasing . Generally, this oscillating pattern would lead to a deterioration of the numerical accuracy. We believe that the described irregularity is again due to the aforementioned improper behavior of orbitals near the nucleus, as the matrix elements of the weak Hamiltonian is accumulated largely in this regime.
To summarize, the DKB set is numerically complete and is well suited for carrying out practical summations over intermediate states in perturbation theory. The comparison with the ND set shows that, at least for the PNC amplitude, the ND convergence pattern becomes affected by the inaccurate representation of the orbitals near the nucleus (this is exemplified for states higher in the spectrum), while the DKB set exhibits a monotonic convergence. Additionally, the DKB set is devoid of spurious states, and therefore the incidental inclusion of spurious states in summations over intermediate states is of no concern.
iii.3 Analysis of spurious states in the ND method
Here we analyze the effect that the additional term has on spurious states which occur in the ND method. We start by taking
(21) 
with an adjustable parameter. For the case of this is equivalent to . We mentioned previously that variationally encourages the boundary condition ; this is also true for any choice of . The specific choice of effectively alters the degree of such variational “encouragement.”
We find that setting results in a single spurious state which appears as the lowest energy eigenstate for each . By subsequently increasing the value of towards we may then deduce the effect that has on these spurious states as well as the rest of the spectrum. We find that small increases in from causes the energy of the spurious state to increase, while the other states remain essentially unaffected (except for the case of near degeneracy with one of these states; this situation is discussed below). As is increased from zero, we may watch as spurious states for each first appear as the lowest energy state, then move up to the second lowest energy state, then to the third lowest energy state, etcetera. Fig. 2 displays this effect for the case of Cs.
It is also interesting to analyze the effect of the spurious state when its energy is nearly degenerate with another state in the spectrum (we will refer to this other state as the “genuine” state). As the energies approach degeneracy by varying , we observe that the spurious state begins to mix in with the genuine state. The first evidence for this is that the energy of the genuine state starts to become affected by the presence of the spurious state. Secondly, we see that the radial wavefunctions and of the genuine state begin to acquire nonphysical “bumps” that oscillate in a way that corresponds to the rapid oscillations of the respective spurious state radial functions. As the energy becomes nearly degenerate, the two states mix to such a degree that it is not even possible to define one as the “genuine state” and the other as the “spurious state.” This effect is shown in Fig. 3, where the genuine state is taken as the 4 state of Cs. As is increased such that the spurious state becomes embedded in the quasicontinuum part of the spectrum, it mixes with multiple neighboring states, and at this point it becomes difficult to track or even define the spurious state.
Presumably increasing to (the Notre Dame case) shifts the spurious states all the way to the end of the spectrum. Evidence for this is seen, for example, in the set used for Fig. 3 (, , bohr). Here the last few levels (excluding the very last level) have an energy spacing on the order of a.u., whereas the energy spacing to the final level is on the order of a.u. Furthermore, increasing past only has a substantial effect on the energy of the final level (for the energy spacing to the final level is then on the order of a.u.).
Up to this point we have been exclusively considering cases with , as these have been the only angular symmetries for which spurious states have previously been known to occur. Now we shall consider the effect that has on the cases of . As with the cases, we see that increasing past only has a substantial effect on the energy of the final level. This is an indication that a spurious state may actually lie at the end of the spectrum for angular symmetries as well. In fact, we find that setting to a negative value significantly below results in a single spurious state which appears as the lowest energy eigenstate for each . By subsequently increasing from this point, we may watch as each spurious state is shifted towards the higher energy end of the spectrum, similar to what is observed with spurious states.
Now we return briefly to subject of the matrix elements and of the ND basis, which are suspected to contribute to poor representation of orbitals near the nucleus. From Eqs. (6) and (15) we see that these include the integral , which is infinite due to the nonvanishing property of the first Bspline at . Numerically this integral is evaluated by Gaussian quadrature, producing finite values. We observe that for the (, , bohr) Cs set, increasing the numerical value of this integral to 20 times its Gaussian quadrature value results in the reappearance of spurious states as the lowest energy eigenstates for each . Simultaneously, the energy of highest energy eigenstate for each is increased substantially. Evidently the capability of to shift the spurious state to the end of the spectrum for angular symmetries is reliant upon the inaccurate numerical evaluation of this (theoretically) infinitevalue integral.
The claim made in this section of observing spurious states for angular symmetries may seem surprising at first. Shabaev et al. (2004) have proved that an independent expansion of large and small radial components with a finite set of basis functions leads to a single spurious state for angular symmetries. This proof assumes the basis functions to vanish at the origin, and the result of this proof is consistent with experience when such basis functions are employed. Arguably, connection is lost immediately with this proof because the ND set includes the first Bspline, which does not vanish at the origin. As we have seen here, the ND method depends on numerical inaccuracies in evaluating infinitevalue integrals in order to manipulate spurious states arising in cases. Because the ND method also includes the first Bspline for angular symmetries (and hence similar numerical inaccuracies), we would have no reason to discount the possibility of spurious states from occurring in these cases as well.
iii.4 Optimizing the basis set: secondorder energy correction
In the preceding sections we established that the DKB sets are more robust than the ND bases. For a comparable size of the basis set the DKB basis exhibits better numerical accuracy for properties accumulated at small radii. Or, we may say that for a fixed numerical accuracy, the DKB basis may contain smaller number of basis functions. In this section we investigate a related question: what the smallest possible basis set is for a given numerical accuracy. Keeping the set as small as possible speeds up manybody calculations that usually require multiple summations over intermediate states. Also, smaller basis sets reduce storage requirements for expansion coefficients in allorder techniques such as configuration interaction or coupledcluster methods.
Quantifying the numerical accuracy requires choosing some metric, which characterizes deviation of the selected property for a given basis from its exact value. Apparently, one should select the “metric” so that it can be easily computed and is related to the relevant atomic properties. As an example of an optimization measure, here we choose the secondorder correction to the energy of a valence electron, .
In the frozencore DHF approximation is the leading manybody correction to the energy. It is given in terms of the Coulomb integrals and single particle energies as (see, e.g., Ref. Johnson (2007)) ,
(22) 
where is the antisymmetrized Coulomb integral. The summations are carried out over core orbitals (labels and ) and virtual (noncore) orbitals (labels and ). Each summation implies summing over principal quantum numbers, angular numbers , and magnetic quantum numbers. The summation over magnetic quantum numbers may be carried out analytically and we are left with summations over radial functions.
We define a contribution of an individual partial wave , , as the difference , where stands for truncated Eq. (22); it includes summations over intermediate states (both core and virtual) with the orbital angular momentum up to . Since the calculations necessarily involve Coulomb integrals between orbitals of different angular momenta, the numerical error in is affected by the accuracy of representation of all partial waves up to .
First, in Table 2, we present results for a large set of basis functions for each partial wave. We use a sufficiently large cavity of bohr in this calculation. These results will serve as a benchmark for comparisons with the less complete (optimized) sets. We observe that the dominant (60%) contribution comes from waves, . Qualitatively, the secondorder energy correction may be described as corepolarization effect. The outer orbitals of the core are relatively “soft” and are easily polarizable. Since the core does not contain and higher partial waves, after peaking at , the partialwave contributions become suppressed as increases. Qualitatively, this suppression arises due to increased centrifugal repulsion of higher partial waves and associated smaller overlaps with core orbitals in the Coulomb integrals of Eq. (22). For example, contributes only 0.6% of the total value.
Now we would like to minimize the size of the set by choosing a different number of radial basis functions for different angular symmetries . To preserve a numerical balance between the finestructure components (for example, this may be important while recovering nonrelativistic limit) we keep the same number of functions for a given orbital angular momentum , e.g., .
In Table 2, we present an example of an optimized basis (marked as “Small set”). We also list resulting numerical errors for each partial wave by comparing with the result from the “Large set” calculations. We see that while the basisset error in higher partial waves is as large as 40%, this hardly makes any influence on the total value of the correlation energy, because of contributions of higher are suppressed. The total value of the correlation correction differs by about 1% from it’s saturated value. Considering that the correlation contribution to the energy is about 10% for the state, the less complete set would introduce only 0.1% error for the total ionization energy.
So far we discussed the correlation energy correction for the state. The optimized set remains sufficiently robust for other lowlying states as well. We have carried out a comparison similar to Table 2 for , , , and states. In all these cases the difference between the total values computed with the“Large” and “Small” sets is about 0.5%. In practical calculations one is often required to reproduce a number of properties with the same set. The atomic properties may be quite dissimilar– like hyperfinestructure interactions accumulated near the nucleus and dipole matrix elements determined by the valence region. Apparently, one has to carry out a similar loworder MBPT analysis for the relevant quantities to verify the suitability of the optimized set.
From the computational point of view, using the optimized sets speeds up the numerical evaluations. In our illustrative example, the “Large set” contains orbitals, while the “Small set” is about four times smaller ( orbitals). The resulting speedup is sizable: our computations of Eq. (22) for the state () are about 14 times faster with the optimized set, as expected due to scaling of the number of contribution in the most computationally demanding second term of Eq.(22). Similar scaling should hold for storage of expansion coefficients in allorder methods, e.g., for storing triple excitations Derevianko and Porsev (2007) one expects scaling of the storage size. Usually higher partial waves produce larger number of angular channels and the scaling should be even steeper than . Further speedup in MBPT summations and reduction in storage size may be attained by skipping a few basis set functions at the upper end of the quasispectrum. Such additional truncation of the spectrum becomes apparent from Fig. 1, where contributions to the PNC amplitude for the upper threequarters ( out of ) of the quasispectrum affect the total value below level of accuracy.
Large set  Small set  Error  

0  (100,11)  0.0000130  (35,7)  0.0000122  6% 
1  (100,11)  0.0020027  (35,7)  0.0019936  0.5% 
2  (100,11)  0.0105623  (30,5)  0.0105373  0.2% 
3  (100,11)  0.0039347  (25,4)  0.0039095  0.6% 
4  (100,11)  0.0007563  (15,4)  0.0007269  4% 
5  (100,11)  0.0002737  (10,4)  0.0002272  20% 
6  (100,11)  0.0001182  (10,4)  0.0000844  40% 
Total  0.0176609  0.0174912  1% 
Iv Conclusion
Calculations of certain atomic properties, such as parityviolating amplitudes, transition polarizabilities between hyperfine levels, and manybody corrections to hyperfine interactions require accurate representation of atomic orbitals at both small and intermediate electronnucleus distances. Solving the manybody problem in high orders of perturbation theory additionally calls for an efficient representation in terms of the basis sets. The dualkineticbasis set is shown here to adequately meet both these demands.
Previously the DKB basis was applied to systems with a single electron in a Coulomb field, i.e., hydrogenlike ionsShabaev et al. (2004); Igarashi (2006). Here we extended the DKB method to manyelectron systems by generating the singleparticle quasispectrum of the DiracHartreeFock potential. Several numerical example for Cs atom were presented. We showed that the DKB method outperforms the widelyemployed Notre Dame Bspline method Johnson et al. (1988) for problems involving matrix elements accumulated at small distances.
In addition, we presented a strategy for optimizing the size of the basis sets by choosing progressively smaller number of basis functions for increasingly higher partial waves. This strategy exploits suppression of contributions of high partial waves to typical manybody correlation corrections.
V Acknowledgments
We would like to thank Walter Johnson, Charlotte FroeseFischer, and Vladimir Shabaev for discussions. The developed dualkineticbalance code for the DiracHartreeFock problem was partially based on codes by the Notre Dame group. This work was supported in part by the National Science Foundation.
References
 C. de Boor, A practical guide to splines (SpringerVerlag, New York, 2001), revised ed.
 H. Bachau, E. Cormier, P. Decleva, J. E. Hansen, and F. Martin, Rep. Prog. Phys. 64, 1815 (2001).
 S. P. Goldman, Phys. Rev. A 31, 3541 (1985).
 H. M. Quiney, I. P. Grant, and S. Wilson, Phys. Scr. 36, 460 (1987).
 W. R. Johnson, S. A. Blundell, and J. Sapirstein, Phys. Rev. A 37, 307 (1988).
 A. Igarashi, J. Phys. Soc. Jap. 75, 114301/1 (2006).
 V. M. Shabaev, I. I. Tupitsyn, V. A. Yerokhin, G. Plunien, and G. Soff, Phys. Rev. Lett. 93, 130405 (pages 4) (2004).
 E. Davidson and D. Feller, Chem. Rev. 86, 681 (1986).
 W. R. Johnson and G. Soff, At. Data and Nucl. Data Tables 33, 405 (1985).
 C. FroeseFischer and F. A. Parpia, Phys. Lett. 179, 198 (1993).
 S. A. Blundell, J. Sapirstein, and W. R. Johnson, Phys. Rev. D 45, 1602 (1992).
 W. R. Johnson, Atomic Structure Theory: Lectures on Atomic Physics (Springer, New York, NY, 2007).
 A. Derevianko and S. G. Porsev, Eur. J. Phys. A (2007).