Reevaluation of Stark-induced transition polarizabilities in cesium

Atomic parity violation (APV) plays an important role in probing the electroweak sector of the standard model (SM) of elementary particles at low energy. The information derived from table-top APV experiments is both complementary to and in competition with that from large-scale particle colliders (see, e.g., the review [1] and references therein). To date, the 1997 Boulder experiment [2] searching for APV in ${}^{133}\mathrm{Cs}$ remains the most accurate. A substantial body of work has been devoted to the interpretation of and the extraction of electroweak observables from the Boulder results.

In its setup, the Boulder experiment [2] employed the ¹³³Cs $6S_{1/2}\rightarrow{7S_{1/2}}$ transition, whose $E1$ amplitude nominally vanishes due to the parity selection rule. However, parity nonconserving (PNC) weak interactions between the atomic nucleus and electrons admix small components of $P_{1/2}$ states into the nominal $S_{1/2}$ states, thus opening the $E1$ channel. Using the parity-mixed multi-electron states $|6S^{\prime}_{1/2}\rangle$ and $|7S^{\prime}_{1/2}\rangle$ and the hyperfine basis (see Eq. (4) below), the APV transition amplitude may be written as

	$\displaystyle A^{\mathrm{PNC}}_{fi}$	$\displaystyle=\langle 7S_{1/2}^{\prime},\,F_{f}M_{f}|-\mathbfcal{E}_{L}\cdot{\boldsymbol{D}}|6S_{1/2}^{\prime},\,F_{i}M_{i}\rangle$
		$\displaystyle=i\mathrm{Im}(E1_{\mathrm{PNC}})\mathbfcal{E}_{L}\cdot\langle F_{f}M_{F_{f}}|\boldsymbol{\sigma}|F_{i}M_{F_{i}}\rangle\,,$		(1)

where $\mathbfcal{E}_{L}$ is the laser electric field driving the $E1$ transition, $\boldsymbol{D}$ is the electric dipole operator, and $\boldsymbol{\sigma}$ is the Pauli matrix.

Due to the smallness of $E1_{\rm{PNC}}$, which is on the order of $\sim 10^{-11}$ in atomic units, measuring the PNC transition amplitude $A^{\mathrm{PNC}}_{fi}$ directly is a formidable challenge. To overcome this difficulty, it was suggested that one uses the Stark-interference technique [3, 4], which relies on the mixing of states of opposite parities due to an externally applied electric field. The transition rate $R$ between the parity-mixed states then includes contributions from the Stark-induced $E1$, magnetic-dipole $M1$, and the PNC-induced amplitudes [5]

$$R\propto|A^{\mathrm{Stark}}_{fi}+A^{M1}_{fi}+A^{\mathrm{PNC}}_{fi}|^{2}\,.$$

(2)

Upon expansion, the right-hand side of Eq. (2) yields the Stark-PNC interference term, $2\mathrm{Re}[A^{\rm{Stark}}_{fi}(A^{\rm{PNC}}_{fi})^{*}]$, whose sign is subject to the handedness of the experiment measuring $R$. Thus, the PNC amplitude $A^{\mathrm{PNC}}_{fi}$ can be extracted from the Stark-PNC interference term by measuring the changes in $R$ under parity reversals. Based on the Stark-interference technique, the Boulder group [5] reported the following values

$$\frac{\mathrm{Im}(E1_{\mathrm{PNC}})}{\beta}=\begin{cases}-1.6349(80)\,\mathrm{mV/cm}\\ \,\,\,\,{\rm for}\,\,\,6S_{1/2},\,{F_{i}=4}\rightarrow{7S_{1/2},\,{F_{f}=3}}\,,\\ -1.5576(77)\,\mathrm{mV/cm}\\ \,\,\,\,{\rm for}\,\,\,6S_{1/2},\,{F_{i}=3}\rightarrow{7S_{1/2},\,{F_{f}=4}}\,,\end{cases}$$

(3)

where $\beta$ is the atomic vector polarizability.

A weighted average of the two values in Eq. (3) yields the nuclear-spin-independent observable, i.e., the nuclear weak charge, while their difference determines nuclear-spin-dependent effects, e.g., the nuclear anapole moment. For the extraction of these quantities, knowledge of the vector transition polarizability $\beta$ is essential and substantial attention [6, 7, 8, 9, 10, 11, 12] has been paid over the years to determining its value. Since 2000, the most accurate value of $\beta$ has been determined based upon combining a semi-empirical calculation of the hyperfine-induced magnetic-dipole $6S_{1/2}\rightarrow 7S_{1/2}$ transition amplitude $M1$ [12] with a measurement of the ratio $M1/\beta$ [11]. Another approach to estimating $\beta$ combines a calculation of the scalar polarizability $\alpha$ with the measurement of the ratio $\alpha/\beta$ [13]. The latest most accurate determination of $\beta$ was published by the Purdue group [10] who adopted the most accurate value of $\alpha/\beta=9.905(11)$ [13] and used the sum-over-state (SoS) method to calculate $\alpha$. Their calculation of $\alpha$ were carried out using experimentally and theoretically determined matrix elements and energies. Although the uncertainties of the $\alpha/\beta$ [13, 10] and $M1/\beta$ [11, 12] approaches are comparable, both approximately at the level of $0.2\%$, their central values differ by $\sim 0.7\%$ or 2.7$\sigma$. This difference critically undermines the accuracy of extracting electroweak observables from the Boulder APV measurement.

Recently, our theory group performed the most sophisticated to date ab initio calculations of the $E1$ transition matrix elements in Cs [14]. Here, we use these newly determined $E1$ matrix elements and a SoS approach to reevaluate the scalar and vector polarizabilities $\alpha$ and $\beta$. We show that the updated value of $\beta$ agrees well with that obtained from the $M1/\beta$ method of Refs. [11, 12], thus reconciling the two alternative approaches (see Fig. 1).

The paper is organized as follows. In Sec. II, we provide a review and derivation for the second-order transition polarizabilities $\alpha$ and $\beta$. In Sec. III, we detail the numerical methods employed for the computation of these quantities. In Sec. IV, we present our numerical values and error estimates and provide a comparison with previous results. Unless stated otherwise, atomic units are used throughout.

In the presence of a DC electric field, the initial and final $S$ states of Cs admix states of opposite parities, thus enabling the otherwise forbidden $E1$ transition between the $6S_{1/2}$ and $7S_{1/2}$ states. In this section, we rederive the conventional results for the Stark-induced $E1$ transition amplitudes. The reader may refer to the original paper by Bouchiat and Bouchiat [4] for an alternative derivation.

We start by introducing the hyperfine basis

$$|n\,({IJ})FM_{F}\rangle=\sum\limits_{M_{J}M_{I}}C^{FM_{F}}_{JM_{J}IM_{I}}|n\,JM_{J}\rangle|IM_{I}\rangle\,,$$

(4)

whose members are formed by coupling electronic states $|n\,JM_{J}\rangle$ of angular momentum $\boldsymbol{J}$ and nuclear states $|IM_{I}\rangle$ of spin $\boldsymbol{I}$ to form states of definite total angular momentum $\boldsymbol{F}=\boldsymbol{I}+\boldsymbol{J}$. Here, $M_{F}$, $M_{J}$, and $M_{I}$ are the magnetic quantum numbers, $n$ stands for the remaining quantum numbers, such as the principal quantum number of the electronic state, and $C^{FM_{F}}_{JM_{J}IM_{I}}$ is the conventional Clebsch-Gordan coefficients.

In the hyperfine basis, Eq. (4), the initial and final states involved in the $i\rightarrow f$ transition are

	$\displaystyle|i\rangle$	$\displaystyle\equiv|n_{i}(IJ_{i})F_{i}M_{i}\rangle\,,$		(5a)
	$\displaystyle|f\rangle$	$\displaystyle\equiv|n_{f}(IJ_{f})F_{f}M_{f}\rangle\,.$		(5b)

In the presence of an externally applied static electric field $\mathbfcal{E}_{S}$, the initial and final states acquire the admixtures

	$\displaystyle|\delta i\rangle$	$\displaystyle=-\sum_{a\neq i}|a\rangle\frac{{\mathbfcal{E}}_{S}\cdot\boldsymbol{D}_{ai}}{\Delta{E_{ia}}}\,,$		(6a)
	$\displaystyle|\delta f\rangle$	$\displaystyle=-\sum_{a\neq f}|a\rangle\frac{{\mathbfcal{E}}_{S}\cdot\boldsymbol{D}_{af}}{\Delta{E_{fa}}}\,,$		(6b)

where $\Delta{E_{ab}}\equiv{E_{a}-E_{b}}$ and $\boldsymbol{D}_{ab}\equiv\langle a|\boldsymbol{D}|b\rangle$ is the electric dipole matrix element.

If a laser is now applied, it can drive the transition $i\rightarrow{f}$, whose Stark-induced $E1$ transition amplitude is given by

$\displaystyle A_{fi}$

$\displaystyle=-\langle\delta f|\boldsymbol{{\mathbfcal{E}}}_{L}\cdot{\boldsymbol{D}}|i\rangle-\langle f|\boldsymbol{{\mathbfcal{E}}}_{L}\cdot{\boldsymbol{D}}|\delta i\rangle={{\mathcal{E}}}_{L}{{\mathcal{E}}}_{S}a_{fi}\,,$

(7)

where $\mathbfcal{E}_{L}$ is the laser electric field. In the last step of Eq. (7), we have factored out the amplitudes of the electric fields and defined the Stark-induced transition polarizability

	$\displaystyle a_{fi}$	$\displaystyle\equiv\sum_{a\neq f}\frac{(\boldsymbol{\hat{\varepsilon}}\cdot{\boldsymbol{D}}_{fa})(\boldsymbol{\hat{e}}\cdot{\boldsymbol{D}}_{ai})}{\Delta{E_{fa}}}$
		$\displaystyle+\sum_{a\neq i}\frac{(\boldsymbol{\hat{e}}\cdot{\boldsymbol{D}}_{fa})(\boldsymbol{\hat{\varepsilon}}\cdot{\boldsymbol{D}}_{ai})}{\Delta{E_{ia}}}\,.$		(8)

Note that $a_{fi}$ still depends on the polarization vectors $\hat{\boldsymbol{e}}\equiv\boldsymbol{\mathcal{E}}_{S}/\mathcal{E}_{S}$ and $\hat{\boldsymbol{\varepsilon}}\equiv\boldsymbol{\mathcal{E}}_{L}/\mathcal{E}_{L}$ of the DC and laser fields.

The expression for $a_{fi}$ may be cast into a form more convenient for angular reduction. To this end, one uses the recoupling identity [15]

	$\displaystyle(R^{(k_{1})}\cdot S^{(k_{1})})(U^{(k_{2})}\cdot{V^{(k_{2})}})$	$\displaystyle=\sum\limits_{Q}(-1)^{Q-k_{1}-k_{2}}$
	$\displaystyle\times\{R^{(k_{1})}\otimes U^{(k_{2})}\}^{(Q)}$	$\displaystyle\cdot{\{S^{(k_{1})}\otimes V^{(k_{2})}}\}^{(Q)}\,,$		(9)

where the operators ${P^{(k_{1})}}$, ${Q}^{(k_{1})}$, ${R}^{(k_{2})}$, and ${S}^{(k_{2})}$ are irreducible tensor operators (ITOs) of ranks $k_{1}$ and $k_{2}$. In Eq. (II), a scalar product of two rank-$k$ ITOs is understood as the following sum over their spherical components

$${P^{(k)}}\cdot{{Q}^{(k)}}=\sum_{q=-k}^{k}(-1)^{q}P^{(k)}_{q}{Q}^{(k)}_{-q}\,,$$

(10)

and a compound ITO of rank $Q$ is defined as

$$\{P^{(k_{1})}\otimes{R^{(k_{2})}}\}_{q}^{(Q)}=\sum_{q_{1}q_{2}}C^{Qq}_{k_{1}q_{1}k_{2}q_{2}}P^{(k_{1})}_{q_{1}}R^{(k_{2})}_{q_{2}}\,,$$

(11)

where $q_{1}$ and $q_{2}$ label the spherical basis components of the ITOs. The possible values of $Q$ are limited by the triangular selection rule, i.e., ${|k_{1}-k_{2}|}\leq{Q}\leq{k_{1}+k_{2}}$.

In our case of the electric dipole couplings, the polarization and dipole operators in Eq. (II) are ITOs of rank $1$. As a result, one has

	$\displaystyle a_{fi}$	$\displaystyle=\sum\limits_{Q=0}^{2}(-1)^{Q}\{\hat{\boldsymbol{\varepsilon}}\otimes\hat{\boldsymbol{e}}\}^{(Q)}\cdot\left(\sum_{a\neq f}\frac{\{\boldsymbol{D}_{fa}\otimes{{\boldsymbol{D}_{ai}}}\}^{(Q)}}{\Delta{E_{fa}}}\right.$
		$\displaystyle\left.+(-1)^{Q}\cdot\sum_{a\neq i}\frac{\{\boldsymbol{D}_{fa}\otimes{\boldsymbol{D}_{ai}}\}^{(Q)}}{\Delta{E_{ia}}}\right)\,,$		(12)

where we have used $\{\hat{\boldsymbol{\varepsilon}}\otimes\hat{\boldsymbol{e}}\}_{q}^{(Q)}=(-1)^{Q}\{\hat{\boldsymbol{e}}\otimes\hat{\boldsymbol{\varepsilon}}\}_{q}^{(Q)}$. The term in Eq. (II) with $Q=0$ corresponds to the scalar, that with $Q=1$ to the vector, and the one with $Q=2$ to the tensor (quadrupole) contributions to the transition polarizability.

To simplify Eq. (II) further, one may introduce the effective ITOs,

$$a[k]_{q}^{(Q)}\equiv\{\boldsymbol{D}\otimes{R_{k}}\boldsymbol{D}\}^{(Q)}_{q}\,,$$

(13)

with the resolvent operator $R_{k}\equiv{(E_{k}-H_{0})^{-1}}$, where $H_{0}$ stands for the unperturbed atomic Hamiltonian. Since $R_{k}$ has the spectral resolution

$\displaystyle R_{k}=(E_{k}-H_{0})^{-1}=\sum\limits_{a\neq{k}}\Delta{E_{ka}^{-1}}|a\rangle\langle a|\,,$

(14)

and is a scalar (so that the combination $R_{k}\boldsymbol{D}$ remains a rank-1 ITO), Eq. (II) may be written as

	$\displaystyle a_{fi}$	$\displaystyle=\sum\limits_{Q=0}^{2}(-1)^{Q}\{\hat{\boldsymbol{\varepsilon}}\otimes\hat{\boldsymbol{e}}\}^{(Q)}$
		$\displaystyle\cdot\left(\langle f|a[f]^{(Q)}|i\rangle+(-1)^{Q}\langle f|a[i]^{(Q)}|i\rangle\right)\,.$		(15)

which, upon applying the Wigner-Eckart theorem, further simplifies to

	$\displaystyle a_{fi}$	$\displaystyle=\sum\limits_{Q=0}^{2}w_{Q}(\hat{\boldsymbol{\varepsilon}},\hat{\boldsymbol{e}})$
		$\displaystyle\times\left[\langle f||a[f]^{(Q)}||i\rangle+(-1)^{Q}\langle f||a[i]^{(Q)}||i\rangle\right]\,,$		(16)

where the multipolar polarization weights $w_{Q}(\hat{\boldsymbol{\varepsilon}},\hat{\boldsymbol{e}})$ are defined as

	$\displaystyle w_{Q}(\hat{\boldsymbol{\varepsilon}},\hat{\boldsymbol{e}})$	$\displaystyle=(-1)^{Q}\sum_{q}(-1)^{q+F_{f}-M_{f}}$
		$\displaystyle\times\begin{pmatrix}F_{f}&Q&F_{i}\\ -M_{f}&-q&M_{i}\end{pmatrix}(\hat{\boldsymbol{\varepsilon}}\otimes{\hat{\boldsymbol{e}}})_{q}^{(Q)}\,,$		(17)

where $\begin{pmatrix}F_{f}&Q&F_{i}\\ -M_{f}&-q&M_{i}\end{pmatrix}$ is the $3j$ symbol.

Finally, by summing over magnetic quantum numbers, we obtain

	$\displaystyle\langle f||a[f]^{(Q)}||i\rangle$	$\displaystyle=(-1)^{F_{i}+I-J_{i}}[F_{f},Q,F_{i}]^{1/2}$
		$\displaystyle\times\begin{Bmatrix}Q&J_{f}&J_{i}\\ I&F_{i}&F_{f}\end{Bmatrix}\sum\limits_{n_{a}J_{a}}\begin{Bmatrix}Q&J_{i}&J_{f}\\ J_{a}&1&1\end{Bmatrix}$
		$\displaystyle\times\frac{\langle n_{f}J_{f}||D||n_{a}J_{a}\rangle\langle n_{a}J_{a}||D||n_{i}J_{i}\rangle}{E_{n_{f}J_{f}}-E_{n_{a}J_{a}}}\,,$		(18)

where $[J_{1},J_{2},\ldots,J_{n}]\equiv(2J_{1}+1)(2J_{2}+1)\ldots(2J_{n}+1)$. The reduced matrix elements $\langle f||a[i]^{(Q)}||i\rangle$ are given by the same formula, but with $E_{n_{i}J_{i}}$ replacing $E_{n_{f}J_{f}}$ in the energy denominator. We point out that due to the $6j$ symbols in Eq. (II), the term with $Q=2$ vanishes for $J_{f}=J_{i}=1/2$, in particular for the transition $6S_{1/2}\rightarrow{7S}_{1/2}$ of interest, as expected.

Conventionally, the Stark-induced transition polarizability $a_{fi}$ is expressed as a linear combination of the second-order scalar and vector polarizabilities [4], $\alpha$ and $\beta$

	$\displaystyle a_{fi}$	$\displaystyle=\alpha\,(\hat{\boldsymbol{e}}\cdot\hat{\boldsymbol{\varepsilon}})\delta_{F_{f}F_{i}}\delta_{M_{f}M_{i}}$
		$\displaystyle+i\beta\,(\hat{\boldsymbol{e}}\times\hat{\boldsymbol{\varepsilon}})\cdot{\langle F_{f}M_{f}|\boldsymbol{\sigma}|F_{i}M_{i}\rangle}\,.$		(19)

These two terms map into the $Q=0$ and $Q=1$ contributions in Eq. (II), respectively. In other words,

	$\displaystyle a_{fi}$	$\displaystyle=-\sqrt{3[F_{f}]}w_{0}(\hat{\boldsymbol{\varepsilon}},\hat{\boldsymbol{e}})\alpha$
		$\displaystyle-\sqrt{2}\langle F_{f}||{\sigma}||F_{i}\rangle w_{1}(\hat{\boldsymbol{\varepsilon}},\hat{\boldsymbol{e}})\beta\,,$		(20)

where, for the $S_{1/2}$ states, the reduced matrix element $\langle F_{f}||{\sigma}||F_{i}\rangle$ in the hyperfine basis (4) is given by

	$\displaystyle\langle F_{f}||{\sigma}||F_{i}\rangle$	$\displaystyle=\sqrt{6}(-1)^{I+F_{i}-1/2}$
		$\displaystyle\times\sqrt{[F_{f},F_{i}]}\begin{Bmatrix}1/2&F_{f}&I\\ F_{i}&1/2&1\end{Bmatrix}\,,$		(21)

where we have used $\langle S=1/2||{\sigma}||S=1/2\rangle={\sqrt{6}}$.

In accordance with Eq. (II), one may then write

	$\displaystyle\alpha=$	$\displaystyle-\frac{\langle f||a[f]^{(0)}||i\rangle+\langle f||a[i]^{(0)}||i\rangle}{\sqrt{3(2F_{f}+1)}}\,,$		(22a)
	$\displaystyle\beta=$	$\displaystyle-\frac{\langle f||a[f]^{(1)}||i\rangle-\langle f||a[i]^{(1)}||i\rangle}{\sqrt{2}\langle F_{f}||{\sigma}||F_{i}\rangle}\,,$		(22b)

or, as explicit sums over intermediate states,

	$\displaystyle\alpha$	$\displaystyle=\delta_{J_{f}J_{i}}\sqrt{\frac{1}{6}}\sum\limits_{n_{a}J_{a}}\frac{(-1)^{J_{a}-J_{i}}}{\sqrt{3(2J_{i}+1)}}$
		$\displaystyle\times\langle n_{f}J_{f}||D||n_{a}{J_{a}}\rangle\langle n_{a}{J_{a}}||D||n_{i}J_{i}\rangle$
		$\displaystyle\times\left(\frac{1}{E_{n_{f}J_{f}}-E_{{n_{a}J_{a}}}}+\frac{1}{E_{n_{i}J_{i}}-E_{{n_{a}J_{a}}}}\right)\,,$		(23a)
	$\displaystyle\beta$	$\displaystyle=-\frac{1}{2}\sum\limits_{n_{a}J_{a}}\begin{Bmatrix}1&J_{i}&J_{f}\\ J_{a}&1&1\end{Bmatrix}$
		$\displaystyle\times\langle n_{f}J_{f}||D||n_{a}J_{a}\rangle\langle n_{a}J_{a}||D||n_{i}J_{i}\rangle$
		$\displaystyle\times\left(\frac{1}{E_{n_{f}J_{f}}-E_{n_{a}J_{a}}}-\frac{1}{E_{n_{i}J_{i}}-E_{n_{a}J_{a}}}\right)\,.$		(23b)

These equations recover the conventional expressions in the literature, see, e.g., formulae in Ref. [16] specialized for the initial and final states of the $S_{1/2}$ character. In this case, the $E1$ selection rules fix the intermediate states to the $P_{1/2}$ and $P_{3/2}$ angular characters.

In the last section, we derived the second order transition polarizabilities $\alpha$ and $\beta$. In this section, we present the numerical methods with which these polarizabilities are calculated. Our approach is a blend of relativistic many-body methods of atomic structure and high-precision experimental values for atomic level energies.

Since the Cs atom has $55$ electrons, its electronic structure is relatively simple: it has a single valence electron outside the [Xe]-like closed-shell core. This simplicity greatly facilitates accounting for many-body effects due to the residual electron-electron interaction (correlation). In what follows, we describe several approximations of increasing complexity through which the correlation contributions to $\alpha$ and $\beta$ are computed.

The lowest-order approximation in the electron-electron interaction is the mean-field Dirac-Hartree-Fock (DHF) method, wherein each electron experiences an “averaged” influence from all other electrons (and of course the Coulomb interaction with the nucleus). Within the DHF approach, we use the “frozen-core” approximation, where atomic orbitals in the [Xe]-like closed-shell core are computed self-consistently, and the valence orbitals are determined afterward in the resulting $V^{N-1}$ DHF potential of the core.

We point out that even at this lowest-order DHF level, the intermediate states involved in calculating $\alpha$ and $\beta$, Eqs. (23), span a countable yet infinite set of bound states and an innumerable set of states in the continuum. Since this Hilbert space is infinitely large, direct numerical summations, while possible, require different numerical implementations for various many-body methods. An elegant way to handle this issue is the $B$-spline approach popularized by the Notre Dame group [17, 18, 19]. This approach generates a finite and numerically complete basis set that has been proven useful in evaluating otherwise infinite sums. In this approach, the set of eigenfunctions is a linear combination of $B$-spline functions covering a radial grid extending from the origin to $R_{\mathrm{max}}$, the radius of an artificially imposed spherical cavity. The Notre Dame approach is further refined by employing the dual-kinetic balance $B$-spline basis set [20] which helps mitigate the issue of spurious states and improve the numerical quality of orbitals both near and far away from the nucleus.

The low-$n$ orbitals from a $B$-spline finite basis set closely resemble those obtained with the conventional finite-difference techniques with a sufficiently large radial grid extent. We refer to these low-$n$ orbitals as “physical” states. As $n$ increases, this mapping deteriorates, so higher-$n$ basis orbitals often differ substantially from their finite-difference counterparts; we refer to such states as “nonphysical” states.

The value of $n$ separating the physical and nonphysical parts of the pseudospectrum primarily depends on $R_{\rm max}$ and to some extent on the number of basis functions. The dependence on the cavity’s radius is easily understood by recalling that low-$n$ orbitals decays exponentially with increasing distance from the nucleus (origin) so they cannot “know about” the existence of a cavity of sufficiently large radius. In contrast, high-$n$ orbitals have their maxima at larger distances and therefore are much more susceptible to the cavity’s presence. Our $B$-spline basis set contains $N=60$ basis functions of order $k=9$ per partial wave generated in a cavity of radius $R_{\mathrm{max}}=250\,\mathrm{a.u.}$. These parameters are chosen so that the fractional differences in the DHF eigenenergies between the basis set and the finite-difference approach for physical states ($n^{\prime}\leq 12$) are within $0.015\%$. Similarly, the basis-set values of the $E1$ matrix elements involving physical states differ from their finite-difference counterparts by less than $0.1\%$. A detailed discussion of the proper mapping of the finite basis set orbitals to the physical states may be found in Ref. [14].

With the finite basis set, one may further facilitate the numerical evaluations of $\alpha$ and $\beta$ by splitting the summations in Eqs. (23a) and (23b) into the “core-valence” (“cv”), “main”, and “tail” contributions

	$\displaystyle\alpha$	$\displaystyle=\alpha_{\rm{cv}}+\alpha_{\rm{main}}+\alpha_{\rm{tail}}\,,$		(24a)
	$\displaystyle\beta$	$\displaystyle=\beta_{\rm{cv}}+\beta_{\rm{main}}+\beta_{\rm{tail}}\,,$		(24b)

where the cv terms correspond to summations over $2\leq{n_{a}}\leq{5}$, the main terms to summations over $6\leq{n_{a}}\leq{12}$, and the tail terms to summations over $13\leq{n_{a}}\leq{\infty}$, respectively. The cv term comes from the core particle-hole intermediate states with excitations to the valence orbital blocked by the Pauli exclusion principle [21]. The infinity in $13\leq{n_{a}}\leq{\infty}$ corresponds to the maximum number of basis set orbitals of a given angular character. We disregard the summation over Dirac negative-energy states, as their contribution in the length-gauge for dipole operators is suppressed by $\alpha_{\rm fs}^{4}$, where $\alpha_{\rm fs}\approx 1/137$ is the fine-structure constant. We have chosen the boundary $n_{a}=12$ between the main and tail terms with the convention of the earlier work, Ref. [10], in mind. Since we have carefully chosen our finite basis set so that that $B$-spline single-electron orbitals with $n_{a}\leq 12$ coincide with their finite-difference counterparts, the intermediate many-body states $|n_{a}J_{a}\rangle$ in $\alpha_{\rm{main}}$ and $\beta_{\rm{main}}$ map into physical states.

The next-level approximation is the Brueckner orbitals (BO) method which incorporates certain many-body effects beyond the DHF treatment. BOs qualitatively describe the phenomenon where the valence electron charge causes the atomic core to become polarized, thus inducing a dipole and higher-rank multipolar moments within the core. Consequently, the redistributed charges within the core attract the valence electron. Compared to the DHF approximation, the BO method improves the theory-experiment agreement for valence electron removal energies. In our work, the BO basis set is obtained by rotating the DHF set using the second-order self-energy operator, see Ref. [14] for further details.

A further improvement upon the DHF and BO methods is the random phase approximation (RPA), which is a linear response theory implemented within the mean-field framework [22, 23]. The primary function of RPA is to account for the screening of externally applied fields by the core electrons. The main advantages of the RPA formalism are that RPA is an all-order method and the RPA transition amplitudes are gauge-independent. For more details about our finite-basis-set implementation of RPA, the reader is referred to Ref. [14]. The RPA(BO) approach incorporates both the core polarization and the core screening effects. The quality of the RPA(DHF) and RPA(BO) dipole matrix elements is substantially improved over the DHF or BO methods, see, again Ref. [14].

To proceed beyond RPA(DHF) and RPA(BO), we employ the all-order relativistic many-body coupled-cluster (CC) approach, which systematically accounts for correlation contributions at each level of approximation. In our recent work [14], $E1$ matrix elements between the $6,7S_{1/2}$ and $nP_{1/2,3/2}$ states for $6\leq n\leq 12$ were computed using the CCSDpTvT method. This method incorporates single (S), double (D), and triple (T) excitation from the reference DHF state [14] in the CC formalism. The “pTvT” qualifier in CCSDpTvT refers to a perturbative treatment of core triples and a full treatment of valence triples. In addition to an accurate treatment of the many-body effects, the CCSDpTvT $E1$ matrix elements values include scaling, dressing, Breit, and QED corrections [14]. These CCSDpTvT values are the most complete theoretical determinations of the $E1$ matrix elements in Cs to date. The results are complete through the fifth order of many-body perturbation theory and include some chains of topologically-similar diagrams to all orders. As such, the CCSDpTvT method is the most theoretically complete applied to correlation effects in Cs so far. Since the finite basis set used in Ref. [14] is identical to that employed in this work, we identify the CCSDpTvT many-body states with $6\leq n\leq{12}$ with the physical states and use the CCSDpTvT matrix elements to compute the main contribution to transition polarizabilities.

We have provided an overview of the numerical approaches employed in our calculations of the transition polarizabilities $\alpha$ and $\beta$. In what follows, we present our numerical results and estimates of uncertainties.

In Table 1, we compile our numerical results for $\alpha$ and $\beta$. In addition to the DHF, BO, RPA(DHF), RPA(BO), and CCSDpTvT results, we list our values obtained from CC calculations of varying complexity. In particular, in the SD approximation, only linear singles and doubles are included, the CCSD approximation additionally incorporates nonlinear effects, the CCSDvT approximation includes full valence triples on top of CCSD, and finally CCSDpTvT(scaled) indicates a CCSDpTvT value rescaled using experimental values for the removal energies. See Ref. [14] for further details. The final values for $\alpha$ and $\beta$ are obtained by adding to the scaled CCSDpTvT values the Breit, QED, and basis extrapolation contributions to the $E1$ matrix elements, as mentioned in Sec. III. Note that the different CC approximations only apply to the $E1$ matrix elements in Eqs. (23). For the energy denominators, we have used the DHF, BO, RPA(DHF), RPA(BO) values in the corresponding approximations and experimental values for all CC approximations. Note also that the CC approximations were only used to compute the main terms, as mentioned in Sec. III. The cv and tail terms are only calculated up to RPA(BO), since (i) their contributions are much smaller than those of the main terms, (ii) full CCSDpTvT calculations are expensive, and (iii) the disparity between the high-$n$ states in the tail terms and their physical counterparts is significant. The semi-empirical result for $\beta$ is obtained by dividing our theoretically determined result for $\alpha$ by the experimentally measured ratio $\alpha/\beta=9.905(11)$ [13] (see below for further details on this point).

As shown in Table 1, the BO correction has a larger impact on $\beta$ than on $\alpha$, differing by 9.6% from the DHF value of $\beta$, while being only 2.6% away from the DHF value for $\alpha$. In contrast, the RPA contribution appears to be much more important for $\alpha$ than $\beta$, with the RPA(DHF) and RPA(BO) values for $\alpha$ differing from the DHF and BO values by around 20%, whereas the RPA(DHF) and RPA(BO) values for $\beta$ are only 0.1–0.5% away from the corresponding DHF and BO values. SD shifts $\alpha$ by 0.35% away from its RPA(BO) value while the SD change for $\beta$ is at 2.2%. CCSD moves the SD value for $\alpha$ by 2.5% and that for $\beta$ by 1%. Adding valence triples amounts to a 4.9% shifts for $\alpha$ and a 0.6% shift for $\beta$ while perturbative core triples give rise to a 0.07% shift for $\alpha$ and a 0.2% shift for $\beta$. Semi-empirical scaling of removal energies changes $\alpha$ by 0.2% and $\beta$ by 0.4%, while Breit, QED, and basis extrapolation corrections are at the level of 0.03% for $\alpha$ and 0.8% for $\beta$.

In Table 2, the cv, main, and tail contributions to $\alpha$ and $\beta$ in different approximations are presented explicitly. This allows us to determine the central values for our computations and estimate our uncertainties. The uncertainties in $\alpha_{\rm main}$ and $\beta_{\rm main}$ may be estimated by considering the convergence patterns of these terms across various approximations. Indeed, Fig. 2 shows the diminishing of contributions from terms of higher and higher order in many-body perturbation theory: the RPA contributions are large, the additional effects of nonlinear core singles and doubles and valence triples, although significantly smaller, are still substantial, whereas additional core triples and scaling effects are generally small. As a result of this observation, we estimate our uncertainty $\sigma_{CC}$ as half the difference between the CCSDpTvT and scaled CCSDpTvT values. This uncertainty represents missing contributions from higher-order CC diagrams. The uncertainty $\sigma_{\rm Breit+QED+basis}$ from the Breit, QED, and basis extrapolation contributions are assumed, conservatively, to be half the difference between the final and scaled CCSDpTvT values. The total uncertainties $\sigma_{\rm main}$ in $\alpha_{\rm main}$ and $\beta_{\rm main}$ are obtained by adding $\sigma_{CC}$ and $\sigma_{\rm Breit+QED+basis}$ in quadrature.

The contributions and uncertainties of the tail terms may be estimated by considering how much the DHF, BO, RPA(DHF), and RPA(BO) values for $\alpha_{\rm main}$ and $\beta_{\rm main}$ differ from the final CCSDpTvT results of these main terms. We observe that the RPA(DHF) and RPA(BO) approximations generally give better agreement with the final values, as to be expected since RPA is known to be responsible for a large portion of the electron correlation effects. As a result, we assume that the contributions from the tail terms are the average of the corresponding RPA(DHF) and RPA(BO) values, and that the uncertainties $\sigma_{\rm tail}$ are half of the corresponding RPA(DHF) and RPA(BO) differences.

Finally, since the cv terms are the same in the RPA(DHF) and RPA(BO) approaches, we take these values as our estimates for $\alpha_{\rm cv}$ and $\beta_{\rm cv}$. The uncertainties $\sigma_{\rm cv}$ in these contributions are assumed to be half the corresponding BO and RPA(BO) differences. The total uncertainties in our evaluations of $\alpha$ and $\beta$ are obtained by adding $\sigma_{\rm cv}$, $\sigma_{\rm main}$, and $\sigma_{\rm tail}$ in quadrature.

In Table 3, the main term of $\alpha$ is further broken down into contributions from intermediate states $n_{a}P_{J}$ with different $n_{a}$. This facilitates a detailed comparison between the result of this work and that of Ref. [10]. We first remind the reader that Ref. [10] estimated the contributions from $6,7P_{J}$ by using experimental values for the $E1$ matrix elements between $6,7S_{1/2}$ and these $P$ states. In contrast, we estimate the contributions from $6,7P_{J}$ by using theoretical CCSDpTvT values for the matrix elements from Ref. [14]. For $6,7P_{1/2}$, the two approaches agree quite well, reflecting the fact that the theoretical CCSDpTvT matrix elements $\langle 6,7S_{1/2}||D||6,7P_{1/2}\rangle$ from Ref. [14] are in good agreement with experiments. On the other hand, our estimates for the contributions from $6,7P_{3/2}$ disagree quite substantially with those from Ref. [10], due to tensions between theoretical values for $\langle 6,7S_{1/2}||D||6,7P_{3/2}\rangle$ from Ref. [14] and experimental results.

Another noticeable feature of Table 3 is the significant difference between our values for the contributions from $n_{a}P_{J}$ with $n_{a}=8,\ldots,12$ and those of Ref. [10]. This discrepancy is due to the fact that Ref. [10] used for matrix elements $\langle 6,7S_{1/2}||D||n_{a}P_{J}\rangle$ ($n_{a}=8,\ldots,12$) theoretical values from Ref. [26], which computed them in the SDpT approximation. In Table 3, we present our SD values for the contributions with $n_{a}=8,\ldots,12$. One observes that our SD values generally agree with those used by Ref. [10], with the small deviations coming from the pT contributions and the fact that Ref. [26] used a different basis from ours, with less accurate mapping to the physical states.

We next point out that Ref. [10] estimated the tail contribution by first computing $\alpha_{\rm tail}$ in the DHF approximation, then rescaling this DHF result based on the fact that the DHF values for contributions from $n_{a}=8,\ldots,12$ are $\sim$ 30% higher than the more accurate SDpT values. In this work, we adopt a slightly different method mentioned earlier, where we estimate $\alpha_{\rm tail}$ by averaging the corresponding RPA(DHF) and RPA(BO) values. This approach stems from the observation that for individual contributions to $\alpha_{\rm main}$ from $n_{a}=8,\ldots,12$, the average of our RPA(DHF) and RPA(BO) values agree well with the final CCSDpTvT results. Reassuringly, our final value for $\alpha_{\rm tail}$ is in good agreement with that of Ref. [10]. We note also that while our DHF value for $\alpha_{\rm cv}$ agrees with that of Ref. [10], we choose to estimate this term using the RPA(DHF) and RPA(BO) methods, since these are more complete theoretical treatments.

From Table 3, the origin of the difference between our estimate for $\alpha$ and that of Ref. [10] is also clear. Out of the total disagreement of 2.50 a.u., 1.86 (74%) comes from the disagreement between experimental and theoretical values for $\langle 6,7S_{1/2}|D||6,7P_{J}\rangle$, 0.46 (18%) originates from our use of the CCSDpTvT instead of SDpT values for the main contributions with $n_{a}=8,\ldots,12$, and the remaining 0.20 (8%) comes from the cv contribution.

We close by noting that, as may be observed from the lower panel of Fig. 2, the computation of $\beta_{\rm main}$ does not converge as well with increasingly complex approximations as that for $\alpha_{\rm main}$ (upper panel of Fig. 2). Indeed, whereas the final uncertainty in $\alpha_{\rm main}$ is at 0.075%, the final uncertainty in $\beta_{\rm main}$ is about six times worse, at 0.42%. This may be understood by noting that in $\beta$, contributions from the $nP_{3/2}$ intermediate states add with an opposite sign to those from $nP_{1/2}$, whereas in $\alpha$, contributions from $nP_{1/2}$ and $nP_{3/2}$ add with the same sign, due to the prefactor $(-1)^{J_{a}-J_{i}}$ in Eq. (23a). Since the $nP_{1/2}$ and $nP_{3/2}$ states are degenerate in the nonrelativistic limit, $\beta$ is nonzero solely due to relativistic effects and is thus suppressed compared to $\alpha$. This may also be understood from the observation that the matrix element of a rank-1 tensor (the vector polarizability) between the $L=0$ states (the $S$ states in the nonrelativistic limit) vanishes due to the angular selection rules, while the same matrix element between the $S_{1/2}$ ($J=1/2$) states does not. The cancellation of terms and the resulting suppression of $\beta$ render SoS computations of the vector polarizability less reliable than those for $\alpha$. An improved evaluation of $\beta$ involves, as in previous works, combining our theoretically determined value of $\alpha=-266.30(21)$ with the experimentally measured ratio [13] $\alpha/\beta=9.905(11)$ to obtain $\beta=26.887(38)$. This semi-empirical value (in bold) for $\beta$ is also presented in Table 1. It differs from the value of $\beta=27.139(42)$ of Ref. [10] by 0.94% or 4.4$\sigma$ while is only 0.26% or 1.1$\sigma$ away from the $M1/\beta$ value of $\beta=26.957(51)$ of Ref. [12]. A comparison between our new value for $\beta$ with previous results is presented in Fig. 1. We conclude that our determination of $\beta$ brings the two alternative approaches ($\alpha/\beta$ and $M1/\beta$) into an essential agreement.

Finally, a weighted average of our value for $\beta$ and that of Ref. [12] results in

$$\beta=26.912(30)\,.$$

This is the most accurate determination of the vector transition polarizability in Cs to date. Since these two values (ours and that of Ref. [12]) were obtained using different methods, potential cross-correlation effects are anticipated to be suppressed when taking the weighted average. Note that taking weighted average over all the values in Fig. 1 would be incorrect, as all the values on the left panel are statistically correlated.

We thank D. Elliott for a discussion. This work was supported in part by the U.S. National Science Foundation grants PHY-1912465 and PHY-2207546, by the Sara Louise Hartman endowed professorship in Physics, and by the Center for Fundamental Physics at Northwestern University.