format: text

text:
# Density-functional calculation of ionization energies of current-carrying atomic states

## I Introduction

Density-functional theory (DFT) hk  is based on identification of the single-body charge density $`n(𝐫)`$ as key variable in terms of which all ground-state observables of an interacting many-electron system can be described dftbook ; parryang . Although $`n(𝐫)`$ is in principle sufficient for this purpose, in practice it often turns out to be advantageous to employ additional variables. The most commonly employed such additional variable is the spin magnetization $`𝐦(𝐫)`$ (or the spin-resolved charge densities $`n_{}(𝐫)`$ and $`n_{}(𝐫)`$), leading to spin-density functional theory (SDFT) dftbook ; parryang ; gunlund . Other choices have occasionally been found useful, e.g. in solid-state physics scdft ; sdwdft .

In the present paper we are interested in orbital magnetism produced by currents forming in open-shell atoms. The current density seems a useful extra variable in this situation. In the absence of external magnetic fields ground-state currents can in principle be calculated also by means of (S)DFT. However, in practice the calculation of orbital currents and their effects on observables is impossible in standard formulations of DFT and SDFT, because no explicit prescription for calculating the true (many-body) current density is known in these theories. Of course, one can always calculate the current arising from the single-particle Kohn-Sham (KS) orbitals of SDFT or DFT, but these orbitals are constrained only to reproduce the correct charge (and spin) densities of the interacting many-body system, and there is no guarantee that the current calculated from them bears any relation to the true current.

Orbital magnetism is thus basically out of reach of conventional DFT. In view of the importance of currents, appearing either spontaneously or induced by external magnetic fields in a wide variety of many-body systems, such as atoms and molecules with open shells, nuclear-magnetic resonance, cyclotron resonance, superconductivity, and magnetism of transition metal and rare-earth compounds, it is thus desirable to develop a DFT-based approach that allows to directly address the effects of orbital currents. In the present paper we explore one such formulation of DFT, namely nonrelativistic current-density functional theory (CDFT) vr1 ; vr2 ; vr3 . Relativistic DFT strange ; rajagopal ; macdonald  in principle also provides explicit information on the current, but standard implementations of it are formulated in a spin-only version, which prohibits extraction of information on the currents. Furthermore, the formalism of relativistic DFT is considerably more complicated than that of SDFT.

CDFT is formulated in terms of the charge density $`n(𝐫)`$ and the nonrelativistic paramagnetic current density vector

$`𝐣_p(𝐫)={\displaystyle \underset{k}{}}𝐣_{p,k}(𝐫)`$
$`={\displaystyle \frac{\mathrm{}}{2mi}}{\displaystyle \underset{k}{}}\left[\psi _k^{}(𝐫)\psi _k(𝐫)(\psi _k^{}(𝐫))\psi _k(𝐫)\right].`$ (1)
This current is to be calculated from the CDFT Kohn-Sham equations vr1 ; vr2 ; vr3

$$\left[\frac{1}{2m}\left(\frac{\mathrm{}}{i}\frac{q}{c}𝐀_s(𝐫)\right)^2+V_s^c(𝐫)\right]\psi _k(𝐫)=ϵ_k^c\psi _k(𝐫),$$
(2)
where an upper index ‘$`c`$’ denotes CDFT,

$$V_s^c(𝐫)=v_s^c(𝐫)+\frac{q^2}{2mc^2}\left(𝐀_{ext}(𝐫)^2𝐀_s(𝐫)^2\right),$$
(3)
$$v_s^c(𝐫)=v_{ext}(𝐫)+v_H(𝐫)+v_{xc}^c(𝐫)$$
(4)
and

$$𝐀_s(𝐫)=𝐀_{ext}(𝐫)+𝐀_{xc}(𝐫).$$
(5)
Here $`v_{ext}`$ and $`𝐀_{ext}`$ are external static electric and magnetic potentials, $`v_H`$ is, as usual, the Hartree potential, and $`v_{xc}^c`$ and $`𝐀_{xc}`$ are the exchange-correlation ($`xc`$) scalar and vector potentals of CDFT, respectively vr1 ; vr2 ; vr3 . Gauge invariance of CDFT requires that the $`xc`$ energy $`E_{xc}[n,𝐣_p]`$ depends on the current only through the so called vorticity vr1 ; vr2 ; vr3

$$\nu [n,𝐣_p](𝐫)=\times \frac{𝐣_p(𝐫)}{n(𝐫)},$$
(6)
i.e., is of the form $`E_{xc}[n,𝐣_p]=\overline{E}_{xc}[n,\nu ]`$. This dependence provides a useful constraint on approximate CDFT functionals.

Although CDFT formally solves the problem how to obtain current-related information from DFT, many practical questions remain. One is, clearly, the construction of approximate current-density functionals. A brief summary of progress in this area is given in Sec. II. Another is the actual implementation of CDFT. In practice, a fully self-consistent solution of the CDFT equations is still quite demanding. Furthermore, in many interesting situations the effect of orbital magnetism, while important for a qualitative and quantitative understanding, is relatively small, so that a fully self-consistent treatment of the orbital degrees of freedom is not always required. The question thus arises if one can put the insights, gained from the formalism of full-fledged CDFT, at work within conventional DFT, in order to achieve an improved description of orbital magnetism via a self-consistent solution of the widely implemented traditional Kohn-Sham equation.

A simple answer to this question was given by one of us in Ref. jpert , by pointing out that the CDFT Kohn-Sham equations can be written in the form of the DFT ones plus a remainder that depends explicitly on the $`xc`$ and external vector potentials. In the absence of external magnetic fields and for not too large $`xc`$ vector potentials it then suggests itself to use low-order perturbation theory in order to describe the CDFT modifications to the DFT equations. This idea has been worked out in Ref. jpert , where explicit expressions for the CDFT corrections to a number of important DFT quantities were given. The resulting approach is labeled ‘perturbative CDFT’, or simply pCDFT. A simple example of pCDFT expressions is the correction to the DFT eigenvalue spectrum $`ϵ_k^d`$, which can be cast in the form jpert

$$ϵ_k^c=ϵ_k^d\frac{q}{c}d^3r𝐣_{p,k}^{KS}(𝐫)𝐀_{xc}[n,𝐣_p^{KS}](𝐫).$$
(7)
## II Current-dependent local-density approximation

In the present paper we apply pCDFT to a study of the effect of currents in open shells on atomic ionization energies. To this end, we require an explicit expression for the $`xc`$ energy in the presence of orbital currents. In Refs. vr1 ; vr2  Vignale and Rasolt proposed an extension of the local-density approximation (LDA) of ordinary, charge-only, DFT to the case of CDFT. Their functional takes the form

$$E_{xc}[n,𝐣_p]=d^3rn(𝐫)\overline{e}_{xc}(n,\nu ),$$
(8)
where

$$\overline{e}_{xc}(n,\nu )=e_{xc}(n,0)+\frac{mk_F^0}{24\pi ^2}\left(\frac{\chi _L}{\chi _L^0}1\right)\frac{|\nu (𝐫)|^2}{n(𝐫)}.$$
(9)
Here $`k_F^0`$ is the Fermi wave vector of a non-interacting electron gas, related to the density via $`k_F^0=(3\pi ^2n)^{1/3}`$. The function $`e_{xc}(n,0)`$ is the exchange-correlation energy per particle in the absence of external magnetic fields, and can be approximated, e.g., by the usual LDA or any of the available generalized-gradient approximations (GGA’s). The second term on the right-hand side of Eq. (9) is a weak-field (linear response) expression for the current-dependent part of the functional. The functional defined by the preceeding equations has been used in several CDFT calculations ebert ; handy1 ; handy2 ; handy3  and should be good enough for a first orientation about the size and nature of current-related phenomena.

Many-body effects enter this functional via the ratio of the orbital susceptibilities of the interacting and the noninteracting electron gas, $`s:=\chi _L/\chi _L^0`$. This ratio has been calculated numerically by Vignale, Rasolt and Geldart vrg  for 10 values of the dimensionless density parameter $`r_s`$, which is related to the density by $`n=3/4\pi (r_sa_0)^3`$, where $`a_0`$ is the Bohr radius. To utilize these results in CDFT, Lee, Colwell and Handy (LCH) proposed the expression handy2 ; handy3

$$s_{LCH}(r_s)=(1.0+0.028r_s)\mathrm{exp}(0.042r_s)$$
(10)
as a convenient and accurate analytical interpolation through the numerical data of Ref. vrg . In the original reference handy2  this fit is claimed to have an rms error of $`1.5\times 10^3`$. This appears to be a misprint in handy2 , because on redoing the calculation we find that over the 10 data points provided by Vignale, Rasolt and Geldart vrg  the rms error of expression (10) is a little bigger, namely $`1.9\times 10^3`$. We have also constructed two alternative interpolations through the data of Ref. vrg . With the same number of fitting parameters (three) as in the LCH fit we find that the expression

$$s_3(r_s)=0.99560.01254r_s0.0002955r_s^2$$
(11)
has an rms error of $`1.2\times 10^3`$, while the 5-term expression

$`s_5(r_s)=1.10380.4990r_s^{1/3}+0.4423\sqrt{r_s}`$
$`0.06696r_s+0.0008432r_s^2`$ (12)
has an rms error of only $`2.1\times 10^4`$ over the same 10 data points. Over the range of values of $`r_s`$ spanned by these data ($`1\mathrm{}10`$) the latter expression should thus be prefered, compared to (10) or (11). In the calculations presented below the current flows mainly in a region in which $`0.1<r_s<5`$. We have therefore performed all calculations once with the LCH fit (10) and once with the above 5-term fit (12). As will be seen below, the differences between both sets of results can be considerable. Since our calculations are self-consistent only with respect to the charge density, but perturbative with respect to the current density, we expect, however, that it is only the order of magnitude of the results that is quantitatively reliable, and on this order of magnitude both parametrizations employed agree consistently.

The LCH expression and both alternative fits differ markedly for $`r_s>12`$, i.e., in the extreme low-density limit. For the present calculation this range is less important, but in view of potential future applications of CDFT in the low-density regime it should be noted that the numerical data of Ref. vrg  do not constrain the various fits in that regime. In the opposite limit, $`r_s0`$, the five-term fit (12), although on average significantly more accurate than the LCH expression, does not correctly recover the value at $`r_s=0`$, whereas the LCH fit correctly yields $`s_{LCH}(0)=1`$. However, the limit $`r_s0`$ corresponds to infinite density, and is thus rather unimportant for typical atomic physics applications. Moreover, in this limit the asymptotically exact expression vrg

$$s_{r_s0}(r_s)=1+\frac{\alpha }{6\pi }r_s\mathrm{ln}r_s+0.08483\frac{\alpha r_s}{\pi }+O(r_s^2),$$
(13)
where $`\alpha =(4/9\pi )^{1/3}`$, is available, so that there is no need for data fitting at all in this limit.

In addition to Eq. (8) with (9) several other CDFT functionals have been proposed, but none seems suitable for our purposes. The approach of Ref. handyfctnl  has not yet led to an explicit expression for the current dependence of $`𝐀_{xc}[n,𝐣_p]`$. The functional of Ref. skudlarski  displays quantum oscillations arising from Landau-level filling in the electron gas, and is thus suitable only for extended systems. The functionals of Refs. 2dim ; pade , on the other hand, were designed specifically for two-dimensional systems in the quantum Hall regime. The expression of Ref. becke , finally, is not a vorticity functional. We will return to this last functional in our discussion of the results in Table 2, below. Given the scarcity of suitable functionals, the proposal of Ref. cdftlett  to generate a CDFT functional by means of a set of integral transformations from an input SDFT functional may prove useful in the future. In the present work, however, we restrict attention to the simple linear-response LDA defined by Eqs. (8) and (9).

## III Numerical results for ionization energies of current-carrying atomic states

After this preparatory discussion of the CDFT functional, we now return to the question of CDFT shifts with respect to the DFT eigenvalues. In general the eigenvalues obtained from solution of the KS equations have no rigorous physical meaning, although they can bear a semiquantitative relationship to the true energy spectrum kohnbeckeparr . One exception to this rule is the highest occupied eigenvalue of the KS spectrum, which is known to be the negative of the system’s ionization energy, and as such can be compared directly to experiment. Although once the system is prepared in a current-carrying state there are, in general, CDFT shifts to all KS eigenvalues, we therefore focus in the present paper on the highest occupied one, since it is most tightly connected with experiment.

In comparisons of experimental ionization energies with KS eigenvalues obtained within the local-density approximation (LDA) it typically turns out that the LDA eigenvalues are significantly off. The origin of this problem is known to be the asymptotic behaviour of the LDA $`xc`$ potential, which decays too rapidly, and leads to too weak binding of the outermost electron. Self-interaction corrections (SIC) pz , which correct the wrong LDA asymptotics, significantly improve on the LDA values for ionization energies, and typically are quite close to experimental data krieger . Once a self-interaction correction has been applied, and the main error of the LDA ionization energies removed, one can consider the effect of additional small corrections, such as the effects of orbital currents. Such current-dependent corrections to the ionization energy can become important, e.g., when ionization (or transfer of electrons during formation of chemical bonds) takes place in the presence of external magnetic fields, since such fields polarize the atom and can give rise to orbital currents.

It is in this situation, ionization or electron transfer in the presence of static external magnetic fields, where pCDFT calculations of current-induced shifts of the ionization energies are directly applicable. Within the context of density-functional theory, the possibility of small current-induced shifts in the ionization energies is also relevant for the calculation of excited states from time-dependent DFT (TDDFT). A recent systematic investigation tddft  of sources of error in excitation energies calculated from TDDFT concludes that the most important ingredient in such calculations is the ground-state $`xc`$ potential used to generate the Kohn-Sham response function, and that ‘the most important requirement for such a potential would be that its highest occupied eigenvalue reproduces the experimental ionization potential as closely as possible’ (cf. Sec. 6 of Ref. tddft ). TDDFT calculations of excitations from open shells are thus expected to sensitively depend on current-induced shifts of the ionization energies, if the excitation takes place in the presence of external magnetic fields.

A separate issue is whether a CDFT calculation of the type presented here can also be useful in the absence of any external magnetic field. An argument for such utility could run along the following lines: The negative of the highest occupied KS eigenvalue of DFT gives the true ionization energy. The negative of the highest occupied eigenvalue of CDFT also gives this energy. In the absence of external magnetic fields both values must thus be identical, if one works with the exact functionals. In practise, of course, we do not have the exact functionals available, and it becomes a meaningful question to ask which of the two approximate eigenvalues is closer to the experimental energy. Our original expectation (not consistently confirmed by the numerical data shown below) was that for atoms whose many-body ground state has a nonzero value of the total angular momentum quantum number $`L`$ the CDFT value would be better than the DFT one, if the CDFT-KS calculation is performed for a system prepared in a state with an orbital current.

Motivated by these considerations we have numerically calculated the pCDFT correction to the highest occupied KS eigenvalue for a series of atoms with unfilled $`p`$ and $`d`$ shells, prepared in a current-carrying state. Energies and radial wave functions of the unperturbed system were obtained numerically from a standard spherically averaged DFT calculation, using the basis-set-free program opmks opmks , both within the LDA and the BLYP GGA blyp . The full single-particle orbital for a current-carrying state was then obtained by multiplying the radial wave function with a spherical harmonic corresponding to a definite value of the magnetic quantum number $`m`$. Such states carry a current proportional to $`m`$. This procedure defines how we prepare a current-carrying state in the KS system. Experimentally, the preparation of an open-shell atom in a current-carrying state is achieved by applying suitable external magnetic or electric fields, as in the Zeeman and Stark effects footnote1 .

In Table 1 the results of LDA-SIC and pCDFT calculations of atomic ionization energies are listed for a number of current-carrying states. It is interesting to see that the differences between values for $`\mathrm{\Delta }ϵ^{pCDFT}`$ obtained with both parametrizations of the susceptibility data used in this work, Eq. (10) and Eq. (12), can be substantial, although the order of magnitude predicted by both is the same. On the other hand, it makes little difference for the size of $`\mathrm{\Delta }ϵ^{pCDFT}`$ whether the unperturbed orbitals used for calculating $`𝐣_{p,k}(𝐫)`$ and $`𝐀_{xc}(𝐫)`$ are obtained from LDA or GGA. The values listed in Table 1 were obtained with LDA. For comparison we also performed calculations using the BLYP GGA blyp , but the resulting changes in the pCDFT corrections are consistently smaller than the ones arising from changing the parametrization of the current-dependent part of the functional from (10) to (12). On the other hand, depending on the particular occupation of $`m`$-substates in the noninteracting system the size of both the current and the resulting energy shift can vary considerably. To explore this variation Table 1 lists values for more than one occupation for most atoms.

Comparison of the experimental ionization energies listed in Table 1 with the sum of LDA-SIC-KLI values and the pCDFT correction shows that pCDFT slightly improves agreement with experiment for several second-row elements, but worsens it for the first-row ones. The data are thus inconclusive as to whether pCDFT represents a true improvement on calculations of ionization energies that neglect currents. This is similar to the situation encountered by LCH in the CDFT calculation of nuclear shielding tensors, in Ref. handy3 , where the same functional (9) as here was employed and CDFT was found to yield rather small shifts that did not systematically improve agreement with experiment.

In Ref. handy3  this was attributed to deficiencies in the approximation for the density functional. Although the same could be true here, too, it is in our opinion more likely that the inconclusive comparison between LDA-SIC-KLI+pCDFT ionization energies with experiment is due to the fact that while the pCDFT values were obtained by explicitly assuming a current-carrying configuration in the open shell, the experimental data refer to configurations with on average equal population of all $`(2L+1)(2S+1)`$ degenerate substates belonging to the $`{}_{}{}^{2S+1}L`$ ground-state term of the atom. Such an average eliminates the orbital currents. A direct comparison of the sum of the LDA-SIC-KLI values and the pCDFT ones with available experimental data is thus not necessarily meaningful. Rather, experiments in the presence of suitable external fields, selectively populating states with a nonzero orbital current, are called for.

Even in the absence of such experiments, it is, however, still possible to compare the size of the pCDFT shifts in Table 1 with other calculations of the effects of orbital currents in open-shell atoms and with the accuracy of present-day DFT calculations. Such a comparison is presented in Table 2, in which we compare the order of magnitude of the pCDFT shifts in Table 1 with the deviation of LDA-SIC calculations from experiment, and with results from another recent calculation of effects of orbital currents in open-shell atoms. From the numbers in Table 2 it is obvious that once the effect of the self-interaction corrections has been taken into account, the order of magnitude of the effect of currents in the open shell is comparable to the remaining difference between the LDA-SIC data and experiment. As explained above, this does not mean that agreement with experiment is necessarily improved by simply adding the two contributions, but it implies that further refinement of density functionals for the calculation of the electronic structure of atoms must take the possibility of curent-dependent energy shifts into account tddft ; jcpbaerends ; becke .

In recent work of Becke becke  the effect of orbital currents in open shells is studied from a different point of view. Instead of explicitly preparing the atom in a symmetry-broken current-carrying state and studying the resulting shifts of single-body energies, as we do here, he considers the total ground-state energies calculated from different but symmetry-equivalent configurations, which should be degenerate. Nevertheless, approximate density functionals can give rise to an artificial breaking of the degeneracy with respect to $`M`$ becke ; baerends . Once this degeneracy is broken, the corresponding ground state carries a current, just as in the present calculations. According to becke  the original degeneracy is restored if the current dependence of the $`xc`$ functional is explicitly accounted for. The current-dependent functional used in Ref. becke  is not constructed within the CDFT of Vignale and Rasolt but based on analysis of the dependence of the curvature of the $`xc`$ hole on the Kohn-Sham current and, as pointed out above, is not a vorticity functional. Implementation of the functional is done not self-consistently (as in full CDFT) or perturbatively (as in pCDFT), but in a post-LDA manner, in which orbitals obtained in a self-consistent (current-independent) LDA calculation are substituted once into the current-dependent functional footnote2 . In spite of these methodological differences between the present work and Ref. becke , Becke also finds that to within the accuracy of today’s density-functionals the current-dependent terms cannot be neglected in high-precision DFT calculations of open-shell atoms. In fact, as illustrated in the third column of Table 2, the size of the current-dependent corrections to the total energy differences obtained from subtracting the values in the last two columns of Table 2 of Ref. becke  is rather similar to the size of the current-induced corrections to $`ϵ_k`$, calculated here.

## IV Conclusions

We have performed CDFT calculations of ionization energies of open-shell atoms prepared in a current-carrying state, with the aim to illustrate the usefulness and viability of CDFT for electronic-structure calculations in the presence of orbital currents. We summarize our conclusions as follows:

(i) The perturbative approximation to CDFT, pCDFT, has been implemented for the first time and was found to be numerically feasible.

(ii) The CDFT $`xc`$ functional has been implemented in the linear-response approximation (8,9) of Vignale and Rasolt vr1 ; vr2 . Three different parametrizations for the orbital susceptibility entering this functional have been tested. The LCH 3-term parametrization was found to be slightly less accurate than claimed in the original reference (rms error $`1.9\times 10^3`$ instead of $`1.5\times 10^3`$). Two alternative parametrizations were developed. Our 3-term parametrization leads to an rms error of $`1.2\times 10^3`$, while the 5-term expression we used in our numerical calculations reduces this to $`2.1\times 10^4`$.

(iii) Orbital currents in open shells result in small but not negligible shifts of the ionization energies. Such shifts can become important, e.g., when ionization takes place in the presence of external electric or magnetic fields. The same applies to the formation and breaking of chemical bonds in the presence of such fields, and to the TDDFT calculation of excitation energies. The calculated shifts, however, do not consistently improve agreement with experiments carried out in the absence of magnetic fields. Presumably this is due to the selection of specific current-carrying states in the calculation, which are averaged over experimentally. A more conclusive test of (p)CDFT would thus require experiments carried out in the presence of magnetic fields, giving rise to a well defined current in the open shell.

(iv) The order of magnitude of the pCDFT terms shows that modern density-functional techniques have reached an accuracy at which the magnitude of small current-related effects arising in open shells is beginning to be significant compared to the remaining difference to experiment. Hence further refinement of DFT-based calculations of atomic spectra should consider the possibility of spontaneous or induced currents in open shells handy3 ; krieger ; jcpbaerends ; becke .

More research is needed to extend this analysis to other atoms and to molecules. Fully self-consistent (neither perturbative nor post-LDA) CDFT calculations would be desirable to this end, as is, obviously, construction of more reliable current-dependent $`xc`$ functionals.

Acknowledgments
This work was supported by FAPESP. One of us (KC) thanks E. K. U. Gross, G. Vignale and S. Kurth for useful discussions, and A. B. F. da Silva for hospitality at the IQSC. We would like to thank E. Engel for providing his code OPMKS on which the numerical calculations of this work were based.