Application of the dual-kinetic-balance sets in the relativistic many-body problem of atomic structure

# Application of the dual-kinetic-balance sets in the relativistic many-body problem of atomic structure

## Abstract

The dual-kinetic-balance (DKB) finite basis set method for solving the Dirac equation for hydrogen-like ions [V. M. Shabaev et al., Phys. Rev. Lett. 93, 130405 (2004)] is extended to problems with a non-local spherically-symmetric Dirac-Hartree-Fock potential. We implement the DKB method using B-spline basis sets and compare its performance with the widely-employed 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 parity-non-conserving amplitude of the transition, and the second-order many-body 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 many-body correlation corrections.

###### pacs:
03.65.Pm, 31.30.Jv

## I Introduction

Applications of perturbation theory in quantum mechanics require summations over a complete set of states of the lowest-order 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, so-called B-splines 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 B-spline set. Further, the variational Galerkin method is invoked and the solution of the resulting matrix eigenvalue problem produces a quasi-spectrum for the atom. In non-relativistic calculations, the lowest-energy orbitals of the resulting basis set closely agree with those of the unperturbed atom, and calculations of various properties of the low-lying states can be carried out. In particular, one could generate single-particle orbitals in some suitable lowest-order approximation, and use the resulting basis set in applications of many-body perturbation theory (MBPT) .

Application of the outlined approach to the relativistic problems brings in a complication– the appearance in the atomic quasi-spectrum of non-physical “spurious” states. These states appear in the solution of the single-particle radial Dirac equation for angular symmetry, ( orbitals). The spurious states rapidly oscillate and, moreover, spoil the mapping of the generated quasi-spectrum onto the low-energy 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 B-splines in relativistic many-body 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 low-energy 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 “kinetically-balanced” 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 built-in symmetry, the authors refer to their method as the “dual-kinetic-balance” (DKB) approach. In both methods (by contrast to the ND prescription), the spurious states were shown to be completely eliminated from the quasi-spectrum.

Motivated by the success of the DKB method in computing properties of hydrogen-like ionsShabaev et al. (2004); Igarashi (2006), here we investigate the suitability of the DKB method in modeling the spectrum of the (non-local) Dirac-Hartree-Fock (DHF) potential. We compare the performance of the ND and DKB methods by computing various properties of Cs atom: single-particle energies, hyperfine integrals, and the parity-non-conserving 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 many-body calculations (for example, in implementing coupled-cluster 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 B-spline basis set was not addressed yet in relativistic many-body calculations. We illustrate optimizing the basis sets by computing the second-order energy correction for several states of Cs.

This paper is organized as follows. First we recapitulate the Galerkin-type approach to generating a finite-basis set quasi-spectrum 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 single-particle energies, hyperfine integrals (Section III.1), and parity-nonconserving 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 second-order 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 mean-field (Dirac-Hartree-Fock) potential. is in general a non-local potential. Assuming that both potentials are central one may exploit the rotational invariance to parameterize the solutions as

 unκ(r)=1r(iPnκ(r) Ωκm(^r)Qnκ(r) Ω−κm(^r)), (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

 (Vnuc(r)+VDHF(r)+c2)Pnκ(r)+c(ddr−κr)Qnκ(r) =εnκPnκ(r), −c(ddr+κr)Pnκ(r)+(Vnuc(r)+VDHF(r)−c2)Qnκ(r) =εnκQnκ(r).

These radial equations may be derived by seeking an extremum of the following action Johnson et al. (1988)

 S = 12c∫{P(r)^O−Q(r)−Q(r)^O+P(r)}dr+12∫(P(r),Q(r))^VDHF([c]cP(r)Q(r))dr (2) + 12∫Vnuc(r)(P(r)2+Q(r)2)dr−c2∫Q(r)2dr−ε12∫(P(r)2+Q(r)2)dr+ΔSbnd+ΔSspur,

where the -dependent Pauli operators are defined as

 ^Oκ±=ddr±κr.

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 quasi-spectrum. 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 (angularly-decoupled) radial equations, may depend on the angular quantum number . We may expand the large and small components in terms of these bases

 P(r) = 2N∑i=1pi li(r), (3) Q(r) = 2N∑i=1pi si(r),

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

 A→p=εB→p, (4)

where and are matrixes and is the vector of expansion coefficients in Eq.(3). The matrix elements are given by

 Aij=Dij+Vij−2c2Sij+Abndij+Aspurij. (5)

The matrixes entering the definition of correspond to various pieces of the radial Dirac equations,

 Dij =c(∫li(r)ddrsj(r)dr−∫li(r)(κr)sj(r)dr) +c(∫lj(r)ddrsi(r)dr−∫lj(r)(κr)si(r)dr), (6) Vij =∫Vnuc(r)[li(r)lj(r)+si(r)sj(r)]dr +(VDHF)ij, (7) Sij =∫si(r)sj(r)dr, (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 non-orthonormal.

### 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 Dirac-Hartree-Fock potential, we employ the frozen-core approximation. In this method, the calculation is carried out in two stages. First, the core orbitals are computed self-consistently. 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 ,

 (VDHF)ij Unknown environment 'tabular' =(VdirDHF)ij+(VexcDHF)ij, (9) (VdirDHF)ij =∑a∈core(2ja+1) ×∫v0(a,a,r)[li(r)lj(r)+si(r)sj(r)]dr, (10) (VexcDHF)ij =−∑a∈core∑L(2ja+1)ΛκLκa ×∫vL(a,j,r)[li(r)Pa(r)+si(r)Qa(r)]dr, (11)

with the conventionally defined multipolar contributions

 vL(b,a,r)=∫rL[Pa(r′)Pb(r′)+Qa(r′)Qb(r′)]dr′, (12)

and the angular coefficient

 ΛκaLκb=(jajbL−1/21/20)2. (13)

### ii.3 ND and DKB sets

Now we specify the Notre Dame (ND) and the dual-kinetic-balance (DKB) basis sets. Both operate in terms of B-spline functions and first we recapitulate the relevant properties of these splines. A set of B-splines of order is defined on a supporting grid . Usually, the gridpoints are chosen as

 t1 =t2=⋯=tk=0, tn =tn+1=⋯=tn+k=R.

In our calculations the intermediate gridpoints are distributed exponentially. B-spline 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)

 B(k)i≤k∝ri−1, (14)

and at , all splines vanish except for the last spline, .

The Notre Dame set is defined as

 lNDi(r) =⎧⎨⎩B(k)i(r)1≤i≤n0n

It corresponds to an independent expansion of the large and small radial components into the B-spline basis. The DKB set involves the Pauli operators and enforces a “kinetic balance” between contributions to the components:

 lDKBi(r) =⎧⎪⎨⎪⎩−B(k)i(r)1≤i≤n−12c^Oκ− B(k)i−n(r)n

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 high-energy end of the quasi-spectrum by adding the following action Johnson et al. (1988) for and for . This correction may be seen as arising from an artificial function-like 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 finite-size nucleus the radial components behave as

 Pnκ ∝rl+1 and Qnκ∝rl+2 for κ<0 , (17) Pnκ ∝rl+1 and Qnκ∝rl    for κ>0 .

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, Froese-Fischer and Parpia (1993) proposed to impose by discarding the first B-spline of the set (this is the only B-spline 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 power-law behavior of the B-splines, Eq.(14), and match it to the small- limits (17). To satisfy the matching, we need to include the B-splines starting with the sequential number ( must be smaller than the order of the splines ).

 imin=|κ|+1={ℓ+2,κ<0ℓ+1,κ>0. (18)

For  and and this is equivalent to the boundary condition of Ref. Froese-Fischer 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 B-spline 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 B-spline from the set. Examination of the resulting orbitals reveals that the wavefunctions acquire a non-physical 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 low-lying bound states of interest. At the same time, throwing away the last B-spline reduces the number of basis functions and leads to a more compact set.

To summarize, we will use the DKB basis set that includes B-splines 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 B-spline set, e.g., for -waves the total number of B-splines would be . For the ND basis , as all B-splines 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 closed-shell core. As the first step we carry out finite-difference Dirac-Hartree-Fock 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 quasi-spectrum with the DHF energies from the finite-difference 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 quasi-spectrum using both ND and DKB sets and carry out comparisons for single-particle energies and hyperfine integrals in Section III.1 and parity-nonconserving amplitude in Section III.2. Section III.3 contains an analysis of spurious states in the ND method. In Section III.4 we analyze second-order energy corrections and discuss a strategy of optimizing the size of the basis.

### iii.1 Energies and hyperfine integrals

We compare ND and DKB quasi-spectrums with energies obtained using a finite-difference DHF code for the low-lying valence states in Table 1. The ND and DKB calculations were carried out using basis functions for B-splines of order . We used a cavity of radius bohr. For the cavity of this size, only a few lowest-energy 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 point-like nucleus

 IEJ(nκ) =∫drrJ+1(P2nκ(r)+Q2nκ(r)), IMJ(nκ) Missing or unrecognized delimiter for \left

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 two-three significant figures for the magnetic-dipole coupling, it produces wrong results for electric-quadrupole and magnetic-octupole 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 electric-dipole 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.

### iii.2 Parity-nonconserving 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 quasi-spectrum. 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 parity-nonconserving (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

 EPNC = ∑n=2⟨7s1/2|D|np1/2⟩⟨np1/2|HW|6s1/2⟩ε6s1/2−εnp1/2 (19) + ∑n=2⟨7s1/2|HW|np1/2⟩⟨np1/2|D|6s1/2⟩ε6s1/2−εnp1/2.

Here and are electric-dipole and weak interaction (pseudo-scalar) operators, and are atomic energy levels. We will compute this expression in the single-particle approximation. Specifically, the index runs over the entire DHF quasi-spectrum for partial wave, including core orbitals. The weak Hamiltonian reads

 HW=GF√8QWρnuc(r)γ5, (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

 EFDPNC = −0.740, EDKBPNC = −0.7395(N=40,k=7), ENDPNC = −0.8546(N=40,k=7).

The finite-difference 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 low-lying 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

 Missing or unrecognized delimiter for \right (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 non-physical “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 quasi-continuum 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 non-vanishing property of the first B-spline 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) infinite-value 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 B-spline, which does not vanish at the origin. As we have seen here, the ND method depends on numerical inaccuracies in evaluating infinite-value integrals in order to manipulate spurious states arising in cases. Because the ND method also includes the first B-spline 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: second-order 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 many-body calculations that usually require multiple summations over intermediate states. Also, smaller basis sets reduce storage requirements for expansion coefficients in all-order techniques such as configuration interaction or coupled-cluster 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 second-order correction to the energy of a valence electron, .

In the frozen-core DHF approximation is the leading many-body correction to the energy. It is given in terms of the Coulomb integrals and single particle energies as (see, e.g., Ref. Johnson (2007)) ,

 E(2)v=∑abn~gabvngvnabεv+εn−εa−εb−∑mna~gvamngmnvaεm+εn−εv−εa, (22)

where is the antisymmetrized Coulomb integral. The summations are carried out over core orbitals (labels and ) and virtual (non-core) 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 second-order energy correction may be described as core-polarization 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 partial-wave 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 fine-structure components (for example, this may be important while recovering non-relativistic 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 basis-set 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 low-lying 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 hyperfine-structure interactions accumulated near the nucleus and dipole matrix elements determined by the valence region. Apparently, one has to carry out a similar low-order 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 speed-up 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 all-order 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 speed-up 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 three-quarters ( out of ) of the quasi-spectrum affect the total value below level of accuracy.

## Iv Conclusion

Calculations of certain atomic properties, such as parity-violating amplitudes, transition polarizabilities between hyperfine levels, and many-body corrections to hyperfine interactions require accurate representation of atomic orbitals at both small and intermediate electron-nucleus distances. Solving the many-body problem in high orders of perturbation theory additionally calls for an efficient representation in terms of the basis sets. The dual-kinetic-basis 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., hydrogen-like ionsShabaev et al. (2004); Igarashi (2006). Here we extended the DKB method to many-electron systems by generating the single-particle quasi-spectrum of the Dirac-Hartree-Fock potential. Several numerical example for Cs atom were presented. We showed that the DKB method outperforms the widely-employed Notre Dame B-spline 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 many-body correlation corrections.

## V Acknowledgments

We would like to thank Walter Johnson, Charlotte Froese-Fischer, and Vladimir Shabaev for discussions. The developed dual-kinetic-balance code for the Dirac-Hartree-Fock problem was partially based on codes by the Notre Dame group. This work was supported in part by the National Science Foundation.

### References

1. C. de Boor, A practical guide to splines (Springer-Verlag, New York, 2001), revised ed.
2. H. Bachau, E. Cormier, P. Decleva, J. E. Hansen, and F. Martin, Rep. Prog. Phys. 64, 1815 (2001).
3. S. P. Goldman, Phys. Rev. A 31, 3541 (1985).
4. H. M. Quiney, I. P. Grant, and S. Wilson, Phys. Scr. 36, 460 (1987).
5. W. R. Johnson, S. A. Blundell, and J. Sapirstein, Phys. Rev. A 37, 307 (1988).
6. A. Igarashi, J. Phys. Soc. Jap. 75, 114301/1 (2006).
7. V. M. Shabaev, I. I. Tupitsyn, V. A. Yerokhin, G. Plunien, and G. Soff, Phys. Rev. Lett. 93, 130405 (pages 4) (2004).
8. E. Davidson and D. Feller, Chem. Rev. 86, 681 (1986).
9. W. R. Johnson and G. Soff, At. Data and Nucl. Data Tables 33, 405 (1985).
10. C. Froese-Fischer and F. A. Parpia, Phys. Lett. 179, 198 (1993).
11. S. A. Blundell, J. Sapirstein, and W. R. Johnson, Phys. Rev. D 45, 1602 (1992).
12. W. R. Johnson, Atomic Structure Theory: Lectures on Atomic Physics (Springer, New York, NY, 2007).
13. A. Derevianko and S. G. Porsev, Eur. J. Phys. A (2007).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters