Calculation of Polarizabilities for Atoms with Open Shells

: A version of the conﬁguration interaction method for atoms with open shells (the Conﬁguration Interaction with Perturbation Theory—CIPT method, PRA 95, 012503 (2017)) is extended for calculation of static and dynamic polarizabilities. Its use is demonstrated by calculation of the polarizabilities for the ground and excited states of Er, Tm and Yb. It is proved to be an useful tool in designing a new generation of optical atomic clocks sensitive to new physics.


Introduction
Polarizabilities are important characteristics of atoms which, together with ionisation potential, determine their chemical properties. The values of static and dynamic atomic polarizabilities are also needed for designing optical clocks. Static scalar polarizabilities determine the value of the black body radiation (BBR) shift of the clock frequency, one of the major source of the clock uncertainty. For these reasons, there is a large number of works which study atomic polarizabilities experimentally and theoretically (see, e.g., the reviews [1][2][3] and references therein). In spite of this wealth of information there are still gaps in the data, mostly for excited states of many-electron atoms, including states, which are considered to be used as clock states. The aim of this paper is to develop a method which might help to close these gaps. The method is an extension of the Configuration Interaction with Perturbation Theory (CIPT) approach [4] which was developed for atoms with a large number of electrons in open shells. We consider Er, Tm and Yb atoms. All these atoms have transitions which are used or proposed to be used as clock transitions.
The 1 S 0 to 3 P o 0 transition in Yb is used in a number of labs around the world as one of the most accurate clock transition approaching relative uncertainty of 10 −18 [5][6][7][8]. There is a large number of measurements and calculations for both states of this transition. We use the data as a testing ground for our method. It was suggested in Ref. [9,10] that the electric quadrupole transition between ground and first excited state of Er can be used as a clock transition with suppressed BBR shift. It was demonstrated in Ref. [11] both experimentally and theoretically that the fine structure (FS) transition between ground and first excited state in Tm has similar features. The BBR shift is proportional to the differential polarizability of the clock states. The later was measured in Ref. [11] and turned out to be exceptionally small, ∆α 0 = −0.063(30) a.u. The authors of Ref. [11] claim that their finding "allow the development of lanthanide-based optical clocks with a relative uncertainty at the 10 −17 level". The lanthanide-based optical clocks have another important feature, they are sensitive to new physics beyond the standard model. For example, the FS clock transition in Fm is sensitive to the hypothetical time-variation of the FS constant α. The FS ∼ (Zα) 2 , which means that if α changes in time, the FS interval changes quadratically faster. The value of this so called enhancement factor , K ≈ 2, is larger than in any current optical clocks based on neutral atoms. It is also competitive with recently suggested clock transitions in Yb [12] and Au [13]. Moreover, the enhancement factors in Yb and Au are large but negative, while the enhancement factor in Tm is large and positive. Comparing two clock transitions with large values and opposite sign of K brings further enhancement of the sensitivity to the variation of α.
The lanthanide-based optical clocks also sensitive to the local Lorentz invariance (LLI) violation. It was demonstrated in Ref. [14] (see also [15]) that for strong sensitivity to the LLI violation one needs long-living state with open f -shell and large value of the total angular momentum J, J ≥ 1. It was suggested in Ref. [14] to use an excited metastable state of Yb + . However, lanthanides satisfy these conditions in the ground state.
The arguments presented above show that the lanthanide-based optical clocks deserve further study. The electronic structure of Er was studied in our previous work [16] using the above mentioned CIPT method. This included energy levels, transition amplitudes, hyperfine structure and polarizability in the ground state. The polarizability was calculated by direct summation over physical states using calculated transition amplitudes. Polarizability of the clock state was not considered. We also study polarizabilities of Er in our earlier works [9,10]. In the approximations used in these works the polarizabilities of ground and clock states were the same. The electronic structure of Tm was studied in [17] the way similar to what we did for Er. However, polarizabilities were not considered. The polarizabilities of the ground state of Er and Tm were studied theoretically and experimentally in a number of works, see review [1] and references therein. However, we are not aware about such works for excited states.
In the present paper we develop an ab initio method of calculation of polarizabilities for complex many-electron atoms. We use it to calculate static dipole polarizabilities of the ground and clock states of Yb, Er and Tm.

Polarizabilities
The second-order Stark shift of the energy of an atom in state a is given by where E is the external electric field. The polarizability α(a) is the sum of scalar and tensor terms Here J a is the total angular momentum of the atom in state a and M a is its projection on the direction of the electric field. Both polarizabilities α 0 (a) and α 2 (a) can be expressed via sums over complete sets of intermediate states involving matrix elements of the electric dipole operator F (in length form Here |a and |n are many-electron atomic states and E a and E n are corresponding energies. Note that the sums over intermediate states at fixed values of J n are the same in (3) and (4). Therefore, calculation of both polarizabilities can be reduced to the calculation of these sums We use the Dalgarno and Lewis method to perform the calculations [18]. We introduce a correction to the atomic wave function |a caused by external field F The correction can be found by solving the equation When the correction is found, the sum (5) is reduced to Then polarizabilities are calculated by substituting (8) into (3) and (4).
Note that the values of the calculated polarizabilities are sensitive to the energy intervals between the state of interest and excited states of opposite parity, see Equations (3)-(5). The calculated energy intervals are nor always sufficiently accurate in atoms with complicated electronic structure. One way around this is to use direct summation and experimental energies in (5). However, summation over physical states is possible only for few low-lying states while the Dalgarno-Lewis approach allows summation over complete set of the many-electron states, including continuum. Therefore, we use a different approach in which we look for a resonant term, i.e., the term which strongly dominates in the sum (5) due to the small energy denominator. If such a term does exist (we use the n 0 notation for the corresponding state), we rescale its value using experimental energy interval. To remove resonant term from (5) we impose the orthogonality condition |δa = |δa − δa|n 0 n 0 |. (9) This leads to a correction to the partial sum S J n This correction is then added back with the rescaling coefficient The rescaling procedure should only be used for resonant contributions. Rescaling of large number of non-dominating contributions leads to inconsistencies between experimental energies and calculated amplitudes, which may worsen the accuracy of the results.

The RPA Method
The equations in the previous section are valid for any number of electrons in the many electron states |a , |n . In practical calculations for many-electron atoms all electrons are divided into two categories, core electrons and valence electrons. This brings the need to include an important effect of core polarisation by external field, i.e., the effect of changing of atomic potential due to the effect of external field on the core electrons.This is done with the use of the so called random-phase approximation (RPA, see e.g., [19]). The RPA equation for the core has the form where H HF is the relativistic Hartree-Fock (HF) Hamiltonian, Here − → α and β are Dirac matrixes, − → p is an operator of electron momentum, V nuc (r c ) is nuclear potential calculated with Fermi distribution of nuclear charge, V core (r c ) is the self-consistent electronic potential.
Index c in (12) and (13) numerates states in the core, δψ c is the correction to the core function due to the effect of external field F, and δV F core is the correction to the core potential due to the change in all core wave functions. It is important that the HF and RPA calculations are done with the same Hamiltonian H HF . It is also important that the valence states calculated with the use of this Hamiltonian represent good initial approximation for the valence states (see next section). Therefore, we use the V N−1 approximations for the HF and RPA calculations for Er, Tm and Yb, with one 6s electron removed. Note that this leads to different definition of the core in the HF and RPA calculations on one side and configuration interaction (CI) calculations on the other side. In the former case the 4 f and 6s electrons are attributed to the core, while in the CI case they are attributed to valence space. In the end, the many electron states |a , |n in previous section are the fourteen to sixteen electron states, while the operator of external field includes core polarisation correction F → F + δV F core .

The CIPT Method
We use the CIPT method [4] to perform the calculations. It is based on the following assumptions: (a) all many-electron basis states are ordered according to their energies; (b) a relatively small number of low-energy states is sufficient to build a good approximation to the states of interest, while contributions of all other states of higher energy can be considered as a small correction; (c) off-diagonal CI matrix elements between high-energy states can be neglected. Keeping these assumptions in mind we can write the total matrix CI equation in a block form Here block A x corresponds to low-energy states, block D x corresponds to high-energy states, and blocks B x and C x correspond to cross terms. Note that since the total CI matrix is symmetric we have C x = B x , i.e., c ij = b ji . Finding X 2 from the second equation of (14) and substituting it to the first one leads to where I is unit matrix. Neglecting the off-diagonal matrix elements in D x leads to very simple structure of the (E a I − D [4] for more details).
The CIPT method described above was used in a number of calculations for many-electron atoms (see, e.g., [16,17]) proving its usefulness. In present paper we go further applying similar approach to Equation (7). Rewriting Equation (7) in a block form leads to Note that blocks A y , B y , C y , D y in (16) are different from blocks A x B x , C x , D x in (14). The later correspond to the wave function, the former -to the correction to this wave function. Since electric dipole operator F change parity, these sets of matrix blocks correspond to states of different parity. Finding Y 2 for second equation of (16) and substituting it to the first one leads to (we neglect high-order term containing FX 2 ) (A y − E a I) + B y (E a I − D y ) −1 C y Y 1 = −FX 1 .
Here again neglecting the off-diagonal matrix elements in D y leads to simple structure of the (E a I − D y ) −1 matrix, (E a I − D y ) −1 ik = δ ik /(E a − E k ). When Y 1 is found from (17) the partial sum (8) is calculated by Correcting this value by rescaling the resonant contribution using (9)-(11) is possible if corresponding state n 0 is included into the A y matrix. Other words, the resonant state n 0 should be considered as a low-energy state in the basis for the correction to the wave function |a .

Ytterbium
We used ytterbium as a testing ground for the method. Polarizabilities of ytterbium are very well studied both theoretically and experimentally (see, e.g., reviews [1,2]) motivated by the use of Yb as one of the most precise optical clock. We used the V N−1 approximations for the HF and RPA calculations with one 6s electron removed, i.e., the initial calculations were done for the 4 f 14 6s configuration of external electrons. The single-electron basis states were then calculated in the frozen V N−1 potential using the B-spline technique [20]. The 4 f and 6s electrons were attributed to valence space in the CIPT calculations. This allowed us to include configurations with excitations from 4 f subshell. However, correlations with core electrons below 4 f were not included. As was explained in our previous work [21], for the sake of calculation of polarizabilities, the mixing with states containing excitations from 4 f was not important, and Yb could be treated as a two-valence-electron atom. All most accurate calculations for Yb are done this way (see [1,2,21] and references therein). The aim of the current calculations was not to improve accuracy but to test the technique which then could be used for atoms with more complicated electronic structure, like Er and Tm (see next section).
The results of the calculations for Yb are presented in Table 1 and compared to other data. Note that we present only those other data which are sufficient to illustrate the accuracy of present calculations. For example, the recommended values from a review [1] are based on the analysis of a large amount of data. We also cite our previous calculations where detailed analysis of the uncertainty was performed. These calculations were more accurate than in the present work but the method can be used for Yb but not for Er or Tm. The Ref. [22] is cited as the only experimental paper which reports the polarizability of the clock state ( 3 P o 0 ). The ground state polarizability of Yb was strongly dominated by the single term corresponding to the 6s 2 1 S 0 to 6s6p 1 P o 0 transition. The calculated energy interval was about 15% smaller than in the experiment. As a result, the calculated polarizability was overestimated. The ab initio value was 152 a.u. When we apply the rescaling procedure (9)-(11) the polarizability went down to 143 a.u. This value is presented in the Table. It is in very good agreement with other data.
Similar to the ground state, the polarizability of the clock state of Yb was strongly dominated by the single term corresponding to the 6s6p 3 P o 0 to 5d6s 3 D 1 transition. The ab initio value of the polarizability was 378 a.u. After rescaling procedure this value goes down to 340 a.u. This was still about 10% larger than other most accurate data (see Table 1). This is a feature of the CIPT method, accuracy is lower for higher states. This is because configuration mixing was stronger for higher states which means that it was harder to build a good initial approximation for the wave function using only low-energy basis states. One way around this is to increase the size of the effective CI matrixes A x and A y (see Equations (14) and (16)). In the end, given the complexity of the targeted systems, the 10% accuracy was acceptable. We stress that the accuracy for the low-lying states is likely to be higher.  (6) 302 (14) Ref. [21] 139 (6) Ref. [1], Recommended value From 134.4(1.0) to 144.2(1.0) From 280.1(1.0) to 289.9(1.0) Ref. [22], Experimental values

Erbium and Thulium
The calculations for Er and Tm were very similar to the case of Yb. Initial HF and RPA calculations were done for the 4 f 12 6s configurations of Er and the 4 f 13 6s configuration of Tm. There were no obvious resonant contributions, therefore, the rescaling procedure was not used. Energy levels of Er and Tm, relevant to the calculation of polarizabilities, are presented in Table 2. It was demonstrated in Ref. [9,10] that calculated polarizabilities of lanthanoids were strongly dominated by the 6s to 6p transitions. Therefore, only states of the 4 f n 6s6p (n = 12 for Er and n = 13 for Tm) were included into the low-energy matrix A y (see Equation (16)). States with excitations from the 4 f subshell were also included but to the high-energy blocks B y , C y , D y , which means that off-diagonal matrix elements between these states were neglected. The results are presented in Table 3. Table 2. Energies of the lowest states of Er and Tm (cm −1 ) which contribute the most to the polarizabilities. Er states belong to the 4 f 12 6s6p configuration while Tm states belong to the 4 f 13 6s6p configuration. J n = J a , J a ± 1, where J a the the total angular momentum of the clock states; J a = 6, 4 for Er and J a = 7/2, 5/2 for Tm.   (10) 144 (15) The results for the ground state are in good agreement with the recommended values of Ref. [1]. These values were based on the analysis of the large amount of theoretical and experimental data. There was also about 10% agreement with our previous calculations [9]. In the approximation used in Ref. [9] the polarizabilities of all states of the same configuration were the same. Present calculations showed that the polarizabilities were close but not exactly the same. The experimental differential polarizability of the clock transition in Tm is ∆α 0 = −0.063(30) a.u. [11]. The accuracy of present calculations was lower than the accuracy of this result. We can probably only say that in our calculations ∆α 0 < 1 a.u. This is still in agreement with the experimental value. Both results, for Er and Tm, showed that the differential polarizability of the clock transitions was small. This means that the BBR shift of the frequencies of these transitions was strongly suppressed.

Further Developments and Applications
The method presented in this paper can be used for calculation of dynamic polarizabilities. All one needs to do is to replace E a in (16) and (17) by E a ±hω, where ω is the frequency of the laser field. One has to keep in mind that when the energy E a +hω comes close to one of the energies E n , (see Equations (5) and (6)) the rescaling of the resonant contribution might be needed for accurate results (see Equations (9)-(11)).
The method can also be used for atoms and ions with open d or p shells. In principle, it can be used for atoms with open f -shell which have more vacancies in it than the Er and Tm atoms. However, it is numerically chalanging task because of the large number of many-electron basis states generated by exciting electrons from an open f -shell. Using supercomputers with parallel calculations is adviceable in this case.
In the end we state that a method of calculating atomic polarizabilities for complicated systems with open shells has beed developed. The method is validated by calculating static polarizabilites of ground and excited clock states of Yb, Er, and Tm. Comparison with other data shows reasonably good agreement. Calculation for clock states of Er and Tm confirm an important claim that these polarizabilities are almost the same as in ground state. This makes the clock transitions to be insensitive to the BBR shift.